Опубликован: 30.05.2014 | Уровень: для всех | Доступ: платный | ВУЗ: Нижегородский государственный университет им. Н.И.Лобачевского

Самостоятельная работа 4: Оптимизация расчетов на примере задачи вычисления справедливой цены опциона Европейского типа

Версия 6. Эквивалентные преобразования. Вычисление квадратного корня

Обратим внимание на еще один момент. Замена деления умножением может привести к существенному выигрышу производительности. Еще одна возможная оптимизация – замена выражения вида 1/sqrtf() на вызов функции invsqrtf(). Проверим, проделал ли компилятор это преобразование самостоятельно.

__declspec(noinline) void GetOptionPricesV6(float *pT,
  float *pK, float *pS0, float *pC)
{
  int i;
  float d1, d2, erf1, erf2, invf;
  float sig2 = sig * sig;

#pragma simd
  for (i = 0; i < N; i++)
  {
    invf = invsqrtf(sig2 * pT[i]);
    d1 = (logf(pS0[i] / pK[i]) + (r + sig2 * 0.5f) * 
         pT[i]) * invf;
    d2 = (logf(pS0[i] / pK[i]) + (r - sig2 * 0.5f) * 
         pT[i]) * invf;
    erf1 = 0.5f + 0.5f * erff(d1 * invsqrt2);
    erf2 = 0.5f + 0.5f * erff(d2 * invsqrt2);
    pC[i] = pS0[i] * erf1 - pK[i] * expf((-1.0f) * r *
            pT[i]) * erf2;
  }
}

Добавьте в массив указателей GetOptionPrices, новую функцию.

tGetOptionPrices GetOptionPrices[9] = 
{ GetOptionPricesV0, GetOptionPricesV1, GetOptionPricesV2,
  GetOptionPricesV3, GetOptionPricesV4, GetOptionPricesV5,
  GetOptionPricesV6 };

Соберите программу и проведите эксперименты.

На описанной ранее инфраструктуре авторы получили следующие времена.

Таблица 9.8. Время работы версий с 0 до 6 (в секундах)
N 60 000 000 120 000 000 180 000 000 240 000 000
Версия V0 17,002 34,004 51,008 67,970
Версия V1 16,776 33,549 50,337 66,989
Версия V2 2,871 5,727 8,649 11,230
Версия V3 0,522 1,049 1,583 2,091
Версия V4 0,521 1,036 1,566 2,067
Версия V5 0,527 1,047 1,580 2,085
Версия V6 0,538 1,071 1,614 2,133

Как видно из 9.8, время работы версии 6 отличается от версий с 3 по 5 в пределах погрешности измерений. Тем не менее, в любом случае проделанные нами оптимизации небесполезны и в случае менее "интеллектуального" компилятора дали бы эффект. Таким образом, именно эту версию мы возьмем за базу в дальнейших экспериментах, включая запуски на сопроцессоре Xeon Phi.

Версия 6.1. Выравнивание данных

Еще одна потенциально полезная программная оптимизация связана с выравниванием обрабатываемых данных, поскольку невыровненные данные обрабатываются процессором медленнее. С большой вероятностью либо с этим справится компилятор, либо ущерб будет не так велик. Тем не менее, давайте убедимся. Для гарантированного выравнивания заменим операторы выделения/освобождения памяти (new/delete) на вызовы функций memalign() 1 Обратите внимание, что функция memalign() не доступна под ОС Windows. Вы можете использовать __mm_malloc() и __mm_free(). /free(). Остальной код не изменится, так же как и ключи сборки.

int main(int argc, char *argv[])
{
  pT = (float *)memalign(32, 4 * N * sizeof(float));
//  pT  = new float[4 * N];

  ...

  free(pT);
//  delete [] pT;
  return 0;
}

Рекомендуемая величина выравнивания (первый параметр в функции memalign()) зависит от длины используемых регистров. Для векторного расширения SSE – это 16, для AVX – 32, и для сопроцессора Xeon Phi – 64.

Кроме того, полезно добавить в код перед циклом еще одну директиву:

#pragma vector aligned

Эта директива подскажет компилятору, что данные, используемые в цикле, выровнены и могут быть использованы соответствующие команды для работы с памятью. Заметим, что если мы "обманем" компилятор, программа будет аварийно завершаться при попытке доступа к невыровненным данным.

Соберите программу и проведите эксперименты. На описанной ранее инфраструктуре авторы получили следующие времена.

Таблица 9.9. Время работы версий с 0 до 6.1 (в секундах)
N 60 000 000 120 000 000 180 000 000 240 000 000
Версия V0 17,002 34,004 51,008 67,970
Версия V1 16,776 33,549 50,337 66,989
Версия V2 2,871 5,727 8,649 11,230
Версия V3 0,522 1,049 1,583 2,091
Версия V4 0,521 1,036 1,566 2,067
Версия V5 0,527 1,047 1,580 2,085
Версия V6 0,538 1,071 1,614 2,133
Версия V6.1 0,539 1,072 1,617 2,135

Как видно из таблицы 9.9, время работы версии 6.1 фактически совпадает с временем работы версии 6. В данном случае эта техника не дала результата, но это не означает, что так будет всегда.

Версия 6.2. Пониженная точность

В п. 4.1 мы указали, что использовать в рассматриваемой задаче тип double нет необходимости. Сейчас мы пойдем еще дальше. В данной задаче точность типа float также избыточна (вряд ли нам могут понадобится больше 4 знаков после запятой для вычисленной цены опциона) и ее можно еще понизить, ускорив тем самым вычисления. Добиться этого несложно, использую ключи компилятора при сборке программы.

Добавьте в командную строку опции, влияющие на точность вычисления значений математических функций.

icc ... -fimf-precision=low -fimf-domain-exclusion=31 

Соберите программу и проведите эксперименты.

На описанной ранее инфраструктуре авторы получили следующие времена.

Таблица 9.10. Время работы версий с 0 до 6.2 (в секундах)
N 60 000 000 120 000 000 180 000 000 240 000 000
Версия V0 17,002 34,004 51,008 67,970
Версия V1 16,776 33,549 50,337 66,989
Версия V2 2,871 5,727 8,649 11,230
Версия V3 0,522 1,049 1,583 2,091
Версия V4 0,521 1,036 1,566 2,067
Версия V5 0,527 1,047 1,580 2,085
Версия V6 0,538 1,071 1,614 2,133
Версия V6.1 0,539 1,072 1,617 2,135
Версия V6.2 0,438 0,871 1,314 1,724

Как видно из табл. 11.10, время работы версии 6.2 примерно на 23% меньше, чем версии 6.1.

Версия 7. Распараллеливание

Распараллеливание вычислительного цикла в функции GetOption-PricesV6() не представляет никакого труда в силу отсутствия каких-либо зависимостей между итерациями. Достаточно указать перед ним прагму omp parallel for и локализовать все переменные, в которые происходит запись.

__declspec(noinline) void GetOptionPricesV7(float *pT, 
  float *pK, float *pS0, float *pC)
{
  int i;
  float d1, d2, erf1, erf2, invf;
  float sig2 = sig * sig;

#pragma simd
#pragma omp parallel for private(invf, d1, d2, erf1, erf2)
  for (i = 0; i < N; i++)
  {
    invf = invsqrtf(sig2 * pT[i]);
    d1 = (logf(pS0[i] / pK[i]) + (r + sig2 * 0.5f) * 
         pT[i]) * invf;
    d2 = (logf(pS0[i] / pK[i]) + (r - sig2 * 0.5f) * 
         pT[i]) * invf;
    erf1 = 0.5f + 0.5f * erff(d1 * invsqrt2);
    erf2 = 0.5f + 0.5f * erff(d2 * invsqrt2);
    pC[i] = pS0[i] * erf1 - pK[i] * expf((-1.0f) * r *
            pT[i]) * erf2;
  }
}

Добавьте в массив указателей GetOptionPrices, новую функцию.

tGetOptionPrices GetOptionPrices[9] = 
{ GetOptionPricesV0, GetOptionPricesV1, GetOptionPricesV2,
  GetOptionPricesV3, GetOptionPricesV4, GetOptionPricesV5,
  GetOptionPricesV6, GetOptionPricesV7 };

Учитывая, что ключ -openmp уже был нами поставлен в командную строку при сборке ранее, изменений в ключах сборки не понадобится.

Соберите программу и проведите эксперименты.

На описанной ранее инфраструктуре авторы получили следующие времена.

Таблица 9.11. Время работы версий с 0 до 7 (в секундах), 16 ядер
N 60 000 000 120 000 000 180 000 000 240 000 000
Версия V0 17,002 34,004 51,008 67,970
Версия V1 16,776 33,549 50,337 66,989
Версия V2 2,871 5,727 8,649 11,230
Версия V3 0,522 1,049 1,583 2,091
Версия V4 0,521 1,036 1,566 2,067
Версия V5 0,527 1,047 1,580 2,085
Версия V6 0,538 1,071 1,614 2,133
Версия V6.1 0,539 1,072 1,617 2,135
Версия V6.2 0,438 0,871 1,314 1,724
Версия V7 0,058 0,084 0,126 0,153

Как видно из таблицы 9.11 ускорение растет с ростом объема данных, от 7,59 до 11.27 на 240 млн. образцов.

Svetlana Svetlana
Svetlana Svetlana

Здравствуйие! Я хочу пройти курс Введение в принципы функционирования и применения современных мультиядерных архитектур (на примере Intel Xeon Phi), в презентации самостоятельной работы №1 указаны логин и пароль для доступ на кластер и выполнения самостоятельных работ, но войти по такой паре логин-пароль не получается. Как предполагается выполнение самосоятельных работ в этом курсе?

Егор Кузьмин
Егор Кузьмин
Россия, г. Москва
Вера Борисова
Вера Борисова
Россия