Самостоятельная работа 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, время работы версии 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, время работы версии 6.1 фактически совпадает с временем работы версии 6. В данном случае эта техника не дала результата, но это не означает, что так будет всегда.

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

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

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

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

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

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

Как видно из табл. 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 ускорение растет с ростом объема данных, от 7,59 до 11.27 на 240 млн. образцов.