Нижегородский государственный университет им. Н.И.Лобачевского
Опубликован: 30.05.2014 | Доступ: свободный | Студентов: 292 / 32 | Длительность: 11:26:00

Самостоятельная работа 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, и занесите их в таблицу.

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

Таблица 9.2. Время работы базовой версии (в секундах)
N 60 000 000 120 000 000 180 000 000 240 000 000
Версия V0 17,002 34,004 51,008 67,970

Версия 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.
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

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

Svetlana Svetlana
Svetlana Svetlana

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