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

Нулевая (базовая) версия

Одна из типичных ошибок, которые допускают многие программисты на языке C (не только начинающие), – смешивание типов float и double при работе с числами c плавающей запятой. Оставляя за рамками обсуждения вопросы точности получаемого результата, рассмотрим, как это влияет на скорость работы, если данные представлены типом float, а используемые математические функции вызываются в форме, работающей с типом double.

В представленном выше коде функции GetOptionPrice() мы использовали, например, функцию logf(). Если программист вызовет вместо нее функцию log(), результат будет следующий. Функция log() в языке C принимает аргумент типа double и возвращает число типа double. Вследствие этого, во-первых, элемент массива pT[i] будет преобразован из типа float в тип double и, во-вторых, далее все вычисления будут вестись с типом double (с преобразованием типов при необходимости) и лишь в самом конце при присваивании результата в элемент массива p1[0] будет выполнено обратное преобразование в тип float. Учитывая, что работа с числами типа float во многих случаях происходит быстрее, чем с типом double (в особенности заметна эта разница становится при использовании векторизации), подобные преобразования сказываются на скорости работы кода отрица тельно. Д авайте посмотрим, насколько велика будет эта разница в нашем случае.

Добавьте в файл main.cpp следующую функцию.

__declspec(noinline) void GetOptionPricesV0(
  float *pT, float *pK, float *pS0, float *pC)
{
  int i;
  float d1, d2, p1, p2;

  for (i = 0; i < N; i++)
  {
    d1 = (log(pS0[i] / pK[i]) + (r + sig * sig * 0.5) *
          pT[i]) / (sig * sqrt(pT[i]));
    d2 = (log(pS0[i] / pK[i]) + (r - sig * sig * 0.5) * 
          pT[i]) / (sig * sqrt(pT[i]));
    p1 = cdfnormf(d1);
    p2 = cdfnormf(d2);
    pC[i] = pS0[i] * p1 - pK[i] * 
            exp((-1.0) * r * pT[i]) * p2;
  }
}

Теперь модифицируем функцию main(). Добавим тип данных – указатель на функцию с прототипом, соответствующим функции GetOption-PricesV0(), объявим массив указателей GetOptionPrices, и по номеру модификации будем вызывать требуемую функцию.

#include "omp.h"

...
// Начало и конец счета
double start, finish;
// Время вычислений
double t;

typedef void (*tGetOptionPrices)(float *pT, float *pK, 
  float *pS0, float *pC);
tGetOptionPrices GetOptionPrices[9] = 
{ GetOptionPricesV0 };

int main(int argc, char *argv[])
{
  int i, version;

  if (argc < 2)
  {
    printf("Usage: <executable> size version [#of_threads]\n");
    return 1;
  }
  N = atoi(argv[1]);
  version = atoi(argv[2]);
  if (argc > 3)
    numThreads = atoi(argv[3]);

  pT  = new float[4 * N];
  pK  = pT + N;
  pS0 = pT + 2 * N;
  pC  = pT + 3 * N;
  for (i = 0; i < N; i++)
  {
    pT[i] = T;
    pS0[i] = S0;
    pK[i] = K;
  }

  float res = GetOptionPrice();
  printf("%.8f;\n", res);

  omp_set_num_threads(numThreads);
  start = omp_get_wtime();
  GetOptionPrices[version](pT, pK, pS0, pC);
  finish = omp_get_wtime();
  t = finish - start;
  printf("v%d: %.8f; %lf\n", version, pC[0], t);

  delete [] pT;
  return 0;
}

Отметим, что для удобства измерение времени выполнено с использованием функции OpenMP (учитывая, что в дальнейшем именно с помощью этой технологии мы будем разрабатывать параллельную версию).

Также отметим, что память под массивы, хранящие данные, мы выделяем одной операцией, настраивая далее указатели pT, pK, pS0 и pC.

Для корректной сборки программы в командную строку нужно добавить ключ -openmp:

icc -O2 -openmp ...

Соберите программу, проведите эксперименты, увеличивая число образцов вдвое от 60 000 000 до 240 000 000, и занесите их в таблицу.

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

Версия 1. Исключение ненужных преобразований типов

Теперь приведем код предыдущей версии в порядок, вызывая правильные функций logf(),sqrtf() и expf() а также правильно указывая константы в коде (все знают, что число 1.0 без суффикса f будет рассматриваться как число типа double?).

__declspec(noinline) void GetOptionPricesV1(float *pT, 
  float *pK, float *pS0, float *pC)
{
  int i;
  float d1, d2, p1, p2;

  for (i = 0; i < N; i++)
  {
    d1 = (logf(pS0[i] / pK[i]) + (r + sig * sig * 0.5f) *
          pT[i]) / (sig * sqrtf(pT[i]));
    d2 = (logf(pS0[i] / pK[i]) + (r - sig * sig * 0.5f) * 
          pT[i]) / (sig * sqrtf(pT[i]));
    p1 = cdfnormf(d1);
    p2 = cdfnormf(d2);
    pC[i] = pS0[i] * p1 - pK[i] * 
            expf((-1.0f) * r * pT[i]) * p2;
  }
}

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

tGetOptionPrices GetOptionPrices[9] = 
{ GetOptionPricesV0, GetOptionPricesV1 };

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

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

Как видно из таблицы 9.3, время работы версии 1 в сравнении с базовой изменилось очень мало. Причина такого результата в функции cdfnormf(), точнее в ее крайне медленной работе, а также в том, что время работы cdfnormf() и cdfnorm() не отличается. На ее фоне выигрыш от исключения ненужных преобразований типов практически не заметен. В других случаях разница может быть очень значительной. Так, авторы данной работы сталкивались с ситуацией, когда разница времени достигала трех раз (!).