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

Понятие опциона и справедливой цены

Опцион – производный финансовый инструмент – контракт между сторонами $$P_1$$ и $$P_2$$, который дает право стороне $$P_2$$ в некоторый момент времени $$t$$ в будущем купить у стороны $$P_1$$ или продать стороне $$P_1$$ акции по цене K, зафиксированной в контракте. За это право сторона $$P_2$$ выплачивает фиксированную сумму (премию) C стороне $$P_1$$. При этом $$K$$ называется ценой исполнения опциона (страйк, strike price), а C – ценой опциона.

В данной работе рассматривается простейший вариант опциона – колл-опцион европейского типа на акцию. Основная идея заключения контракта состоит в игре двух лиц – $$P_1$$ и $$P_2$$. Вторая сторона выплачивает некоторую сумму C и в некоторый момент времени $$T$$ (срок выплаты, maturity, зафиксирован в контракте) принимает решение: покупать акции по цене K у первой стороны или нет. Решение принимается в зависимости от соотношения цены $$S_T$$ и $$K$$. Если $$S_T$$ < $$K$$, покупать акции не выгодно, первая сторона получила прибыль C, а вторая – убыток C. В случае, если $$S_T$$ > $$K$$, вторая сторона покупает у первой акции по цене K, в ряде случаев получая прибыль (в зависимости от соотношения между C и $$S_T-K$$).

Основная проблема заключается в расчете справедливой цены такого опционного контракта – цены, при которой наблюдается баланс выигрыша/проигрыша каждой из сторон. Логично определить такую цену как средний выигрыш стороны $$P_2$$: $$C=E(e^{-rT}(S_T-K)^+)$$

В формуле (4) $$(S_T-K)^+$$ равно либо разности между ценой акции в момент $$T$$ и ценой исполнения (если цена акции больше), либо нулю, $$E$$ – математическое ожидание ($$S_T$$ подвержена действию случайных факторов, смотри формулу (3)). Умножение на экспоненту соответствует т.н. дисконтированию – переводит 1 у.е. в момент времени $$t=T$$ к моменту времени $$t = 0$$, что соответствует инфляции с процентной ставкой r.

Метод вычисления справедливой цены опциона

Для вычисления справедливой цены опциона по формуле (4) известно аналитическое решение, которое справедливо при сделанных нами ранее предположениях.

Аналитическое решение описывается формулой Блэка–Шоулса для вычисления цены опциона в момент времени t = 0:

$$ C=S_0F(d_1)-Ke^{-rT}F(d_2)$$ $$d_1=ln \frac {S_0} K+\frac {\left(r+\frac {\sigma^2} 2 \right) T} {\sigma \sqrt T}$$ $$d_2=ln \frac {S_0} K+\frac {\left(r-\frac {\sigma^2} 2 \right) T} {\sigma \sqrt T}$$ где $$F$$ – функция стандартного нормального распределения.

Именно эту формулу мы и будем использовать в дальнейших расчетах. Что же здесь считать? – спросите вы. Действительно, на первый взгляд все выглядит элементарным. На деле не все так просто, что мы и увидим далее.

Оценивание набора опционов

Конечно, для вычисления цены одного опциона не нужна высокопроизводительная вычислительная техника. На практике организации, работающие на финансовых рынках, вычисляют цены гигантского количества разных опционов, которые можно выпустить в конкретных рыночных условиях. Учитывая, что время финансовых расчетов существенно влияет на скорость принятия решений, каждая секунда на счету. Поэтому сокращение времени оценивания набора опционов является достаточно важной задачей.

Построим схему информационных зависимостей для одного опциона (рис 9.1).

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

Программная реализация

Точка отсчета. Расчет по аналитической формуле

Программную реализацию начнем с написания функции, выполняющей расчет по рассмотренной в разделе 2.3 формуле.

Создадим пустой файл main.cpp и напишем в нем заготовку функции main().

int numThreads = 1;
int N = 60000000;

int main(int argc, char *argv[])
{
  int 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]);

  // Здесь будут вызовы функций для разных способов расчета

  return 0;
}

Здесь version – номер модификации, эксперименты с которой мы будем проводить в конкретный момент времени, N – число образцов, numThreads – число потоков в параллельной версии.

Следующим шагом напишем функцию GetOptionPrice() для расчета цены опциона по формуле (5). Отметим, что здесь и во всех остальных функциях мы используем одинарную точность, то есть тип float. В силу специфики данной задачи двойная точность в ней явно избыточна (например, в связи с тем, что входные данные изначально приходят с точностью, для которой типа float достаточно). Более того, далее мы увидим, что точность можно понизить дополнительно.

Функция GetOptionPrice() предваряется объявлением необходимых констант, используемых в расчете. Обращаем внимание на необходимость согласованного объявления данных – поскольку обычно процентная ставка указывается в годовом исчислении, то и время также должно быть в годах.

// Волатильность (% годовых 0.2 -> 20%)
const float sig = 0.2f;
// Процентная ставка (% годовых 0.05 -> 5%)
const float r   = 0.05f;
// Срок исполнения опциона (в годах)
const float T   = 3.0f;
// Цена акции в момент времени t=0 (у.е.)
const float S0  = 100.0f;
// Цена исполнения - цена покупки, зафикс. в опционе (у.е.)
const float K   = 100.0f;

float GetOptionPrice()
{
  float C;
  float d1, d2, p1, p2;

  d1 = (logf(S0 / K) + (r + sig * sig * 0.5f) * T) /
       (sig * sqrtf(T));
  d2 = (logf(S0 / K) + (r - sig * sig * 0.5f) * T) /
       (sig * sqrtf(T));
  p1 = cdfnormf(d1);
  p2 = cdfnormf(d2);
  C  = S0 * p1 - K * expf((-1.0f) * r * T) * p2;

  return C;
}

Нетривиальная часть формулы (5) – функция стандартного нормального распределения F. Для ее вычисления можно написать собственный довольно несложный код, но учитывая, что это задача не входит напрямую в данную работу, мы использовали готовое решение – функцию cdfnormf().

Наконец, последнее, что осталось сделать в этой части лабораторной работы, – вызвать функцию GetOptionPrice() из функции main() с выводом результата на печать.

int main(int argc, char *argv[])
{
  int 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]);

  float res = GetOptionPrice();
  printf("%.8f;\n", res);
  // Здесь будут вызовы функций для разных способов расчета

  return 0;
}

Для сборки разработанного кода используем следующую командную строку:

icc -O2 main.cpp -o option_prices

Выбор структуры хранения

Обсудим вопрос о хранении данных в нашей программе. Казалось бы, в чем проблема? Нужно завести четыре массива (три для хранения исходных данных и один для хранения результатов), а также несколько отдельных переменных. Дело в том, что порядок расположения данных в памяти существенно влияет на производительность программы. Так, во многих задачах, подобных той, что мы решаем, возникает дилемма: использовать набор массивов (паттерн SoA – structure of arrays) или массив структур (паттерн AoS – array of structures). В первом случае данные будут размещены в памяти так: сначала целиком разместится первый массив, потом целиком второй и т.д. Во втором случае сначала будут размещены все данные, относящиеся к первому объекту предметной области, потом – ко второму и т.д.

Вопрос о том, какой из подходов лучше, не имеет однозначного ответа. В некоторых случаях можно дать рекомендации. Так, например, интуитивно понятно, что в случае, если мы работаем с M массивами, из которых три используются в 95% случаев, а оставшиеся M–3 – в 5% случаев, имеет смысл отдельно хранить три массива и оставшиеся M–3 массива. Такой подход позволит локализовать в памяти те данные, с которыми мы чаще всего работаем.

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

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

Одна из типичных ошибок, которые допускают многие программисты на языке 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() не отличается. На ее фоне выигрыш от исключения ненужных преобразований типов практически не заметен. В других случаях разница может быть очень значительной. Так, авторы данной работы сталкивались с ситуацией, когда разница времени достигала трех раз (!).

Версия 2. Эквивалентные преобразования: cdfnorm() vs. erf()

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

Используем формулу $$cdfnorm(x)=0.5+0.5erf \left( \frac x {\sigma \sqrt T \right)}$$

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

__declspec(noinline) void GetOptionPricesV2(float *pT, 
  float *pK, float *pS0, float *pC)
{
  int i;
  float d1, d2, erf1, erf2;

  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]));
    erf1 = 0.5f + 0.5f * erff(d1 / sqrtf(2.0f));
    erf2 = 0.5f + 0.5f * erff(d2 / sqrtf(2.0f));
    pC[i] = pS0[i] * erf1 - pK[i] * expf((-1.0f) * r * 
            pT[i]) * erf2;
  }
}

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

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

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

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

Как видно из таблицы 9.4, время работы версии 2 почти в 6 раз меньше, чем версии 1. В дальнейших модификациях будем использовать именно функцию erff().

Версия 3. Векторизация: использование ключевого слова restrict

Возможности алгоритмической оптимизации на этом в основном исчерпаны. Перейдем к векторизации расчетного цикла. Прежде всего, необходимо сменить режим компиляции, указав в командной строке ключ –mavx, в противном случае компилятор будет использовать наборы команд SSE, но не AVX.

Проведя эксперименты с собранной версией, можно убедиться, что время работы функции GetOptionPricesV2() не изменится. Чтобы понять, в чем дело, укажем при сборке ключ -vec-report3, "попросив" компилятор выдать отчет о векторизации. Цифра "3" соответствует степени подробности отчета, соответствующей выдаче информации о векторизованных/невекторизованных циклах, а также причинах отсутствия векторизации. Заметим, что в последних версиях компилятора появился новый, шестой тип отчета (-vec-report6), дополнительно дающий некоторые рекомендации относительно векторизации.

Отчет о векторизации приведен ниже (рис 9.2). К сожалению, в данном случае он не слишком информативен. Заметим, что при компиляции данного кода при помощи Intel C/C++ Compiler for Windows в отчете можно увидеть фразу о том, что возможны зависимости между массивами. В целом, наш опыт говорит о том, что диагностики компиляторов Intel для C/C++ чаще всего позволяют догадаться, что именно мешает векторизации.

Вернемся к вопросу о том, как "уговорить" компилятор векторизовать цикл.

Объяснить компилятору, что опасности нет, можно разными способами: используя ключевое слово restrict в объявлении формальных параметров функции, с помощью директивы #pragma ivdep, с помощью директивы #pragma simd перед циклом. Кроме перечисленного, можно указать специальную опцию, которая подскажет компилятору, что массивы в программе никогда не пересекаются, но эта опция действует на всю программу целиком и при неаккуратном использовании может "сломать" код.

Рассмотрим сначала первый вариант. Ключевое слово restrict в объявлении указателей говорит о том, что в данную область памяти не будут обращаться при помощи других указателей. Это дает компилятору информацию о том, что массивы pT, pK, pS0 и pC не пересекаются, а, следовательно, запись pC[i] не оказывает никакого влияния на значения в других массивах, то есть можно безопасно менять порядок вычислений, организуя их векторно.

__declspec(noinline) void GetOptionPricesV3(
  float * restrict pT, float * restrict pK, 
  float * restrict pS0, float * restrict pC)
{
  int i;
  float d1, d2, erf1, erf2;

  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]));
    erf1 = 0.5f + 0.5f * erff(d1 / sqrtf(2.0f));
    erf2 = 0.5f + 0.5f * erff(d2 / sqrtf(2.0f));
    pC[i] = pS0[i] * erf1 - pK[i] * expf((-1.0f) * r * 
            pT[i]) * erf2;
  }
}

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

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

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

icc -O2 –restrict -mavx ...

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

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

Как видно из таблицы 9.5, время работы версии 3 примерно в 5.5 раз меньше, чем версии 2. Процессор, на котором проводились эксперименты, способен с использованием AVX выполнять операции одновременно над 8 элементами типа float. Почему код ускорился в меньшее число раз, чем должен, предлагаем выяснить читателям самостоятельно, собрав профиль работы функции, использую Intel VTune Amplifier XE (обратите внимание на вызовы математических функций и примените закон Амдаля).

Версия 4. Векторизация: использование директивы simd

Второй способ векторизации связан с использованием перед циклом директивы #pragma ivdep (ignore vector dependencies). Она подсказывает компилятору, что с нашей точки зрения массивы в цикле не пересекаются. Заметим, что в случае, когда мы ошиблись, компилятор иногда может установить наличие зависимости и проигнорировать нашу подсказку, но чаще всего мы получим некорректно работающий код. Данной директивой надо пользоваться с осторожностью. Заметим, что эта директива иногда используется в паре с #pragma vector always – подсказкой компилятору о том, что с нашей точки зрения векторизация, если она возможна, будет эффективна. Как мы упоминали ранее в лекции и практике по векторизации, иногда векторизация не приведет к ускорению, более того, может возникнуть замедление. Как правило, это происходит в тех ситуациях, когда данные лежат в памяти не в том порядке, в котором их нужно запаковывать в векторные регистры. В таких случаях накладные расходы на подготовку данных и запись результа тов могут превысить выигрыш от векторного выполнения операций. Однако компилятор не всегда прав. Иногда он считает векторизацию неэффективной, но мы хотим попробовать. Для этого и нужна #pragma vector always.

В настоящий момент основным становится третий способ подсказки компилятору – использование новой директивы #pragma simd. Возможности компилятора по векторизации цикла при использовании данной директивы значительно выше, кроме того, существует целый ряд настроек, которых не было ранее. Заметим, что ответственность за корректность результата в этом случае полностью ложится на программиста. Использовать данную директиву нужно с осторожностью.

__declspec(noinline) void GetOptionPricesV4(float *pT,
  float *pK, float *pS0, float *pC)
{
  int i;
  float d1, d2, erf1, erf2;

#pragma simd
  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]));
    erf1 = 0.5f + 0.5f * erff(d1 / sqrtf(2.0f));
    erf2 = 0.5f + 0.5f * erff(d2 / sqrtf(2.0f));
    pC[i] = pS0[i] * erf1 - pK[i] * expf((-1.0f) * r * 
            pT[i]) * erf2;
  }
}

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

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

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

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

Как видно из таблицы 9.6, время работы версии 4 почти не отличается от времени версии 3. Впрочем, этого и следовало ожидать.

Версия 5. Вынос инвариантов из цикла

Попытаться еще немного ускорить работу функции GetOption-PricesV4() можно, сэкономив не тех вычислениях, которые можно вынести за цикл. В данном случае это расчет 1.0f / sqrtf(2.0f) .

Вычислим это выражение отдельно и внесем в код константу вида

const float invsqrt2 = 0.707106781f;

Теперь используем ее в теле функции, заменив попутно деление на умножение.

__declspec(noinline) void GetOptionPricesV5(float *pT,
  float *pK, float *pS0, float *pC)
{
  int i;
  float d1, d2, erf1, erf2;

#pragma simd
  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]));
    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
};

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

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

Как видно из таблицы 9.7, время работы версии 5 практически совпадает с временами версий 3 и 4. В данном случае, видимо, компилятор проделал необходимые преобразования кода самостоятельно. Убедиться в этом можно, изучив ассемблерный листинг программы после сборки (укажите в командной строке ключ -Fa ).

Заметим, что хотя в данном случае мы ничего не выиграли, вынос инвариантов – полезная техника оптимизации, которая, как правило, делает код более понятным и логичным, а также нередко приводит к выигрышу производительности (компилятор далеко не всегда может осуществить подобные преобразования сам).

Версия 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. Выравнивание данных

Еще одна потенциально полезная программная оптимизация связана с выравниванием обрабатываемых данных, поскольку невыровненные данные обрабатываются процессором медленнее. С большой вероятностью либо с этим справится компилятор, либо ущерб будет не так велик. Тем не менее, давайте убедимся. Для гарантированного выравнивания заменим операторы выделения/освобождения памяти ( Обратите внимание, что функция 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 млн. образцов.

Версия 7.1. Прогрев

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

Обсудим, является ли честным такой способ оценки производительности. Что касается исключения накладных расходов на создание потоков, сомнений нет. Дело в том, что в реальных программах мы всегда можем организовать однократное создание потоков, в результате чего все последующие параллельные расчеты в рамках конкретной программы будут лишены соответствующих накладных расходов. Однако есть еще один важный момент. Дополнительное уменьшение времени может быть получено за счет так называемого "прогрева" кэша. При повторном вызове функции GetOptionPricesV7() часть необходимых ей данных может остаться в кэше, что может ускорить ее работу. Хорошо это или плохо? С одной стороны – плохо или, скорее, не совсем честно. Известно, что типичная ошибка – проведение в рамках одного процесса экспериментов с однопоточной, двухпоточной и четырехпоточной версиями программы, работающими с одними и теми же данными. В результате обычно получается сверхлинейное ускорение, которое немедленно исчезает, когда экспери менты ставятся корректно – каждый в своем процессе. С другой стороны, в реальных программах выполняется много расчетов, и данные постепенно подгружаются в кэш.

В целом, наша рекомендация выглядит следующим образом:

Проведем эксперименты с функциями GetOptionPricesV6() и GetOptionPricesV7().

Модифицируем функцию main() следующим образом.

int main(int argc, char *argv[])
{
  ...
  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);

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

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

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

Как видно из таблицы 9.12, ускорение версии 7.1 по отношению к версии 6.3 существенно выше, чем было у версии 7 по отношению к 6.2, и составляет от 12.54 (60 млн. образцов) до 13,61 (240 млн. образцов). Также видно, что версия 6.3 примерно на 7,5% быстрее, чем версия 6.2, а версия 7.1 – на треть быстрее, чем версия 7 (кроме эксперимента с 60 млн. образцов, где разница составляет более 70%, что неудивительно для такого малого объема вычислений).

Эксперименты на сопроцессоре Xeon Phi

Мощь сопроцессора Xeon Phi проявляется при использовании большого числа потоков (от 60 до 240), при этом каждое отдельное ядро существенно уступает в производительности ядрам процессора Xeon. Таким образом, сравнение имеет смысл проводить, начиная с версии 7. Однако, чтобы картина была полной, рекомендуем читателям провести эксперименты и с начальными версиями. В частности, можно увидеть, какие результаты дает векторизация кода при переходе от версии 2 к версии 3.

Здесь же мы приведем результаты экспериментов, начиная с версии 6, взятой нами за базу для процессора Xeon.

Для сборки программы под сопроцессор Xeon Phi добавьте в командную строку компилятора ключ -mmic.

Также не забудьте увеличить величину выравнивания в функции memalign() с 32 до 64.

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

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

Как видно из таблицы 9.13 на сопроцессоре Xeon Phi переход к вычислениям с пониженной точностью (версия 6.2) дает выигрыш не на 23%, как это было на процессоре, а почти в 2,3 раза. Также видно, что более значительный прирост (порядка 60%) дает прогрев, что неудивительно, принимая во внимание тот факт, что мы создаем на порядок больше потоков. Наконец, можно отметить, что времена работы версии 6.3. на сопроцессоре практически сравнялись с временами работы на процессоре. Теперь посмотрим на результаты работы параллельных версий.

Как видно из табл. 9.14-9.16, версия 7 ускоряется весьма плохо. Как и для центрального процессора, накладные расходы на работу с потоками оказываются сопоставимыми с общим временем работы программы. В то же время за вычетом расходов на работу с потоками (версия 7.1) ускорение получается вполне неплохим (от 50,5 до 60,4). Также видим, что при запуске в 120 потоков ускорение выше, чем при запуске в 60 потоков, что согласуется с техническими особенностями исполнения кода на Xeon Phi.

Однако результаты на Xeon Phi можно улучшить.

Версия 8. Оптимизация работы с кэш-памятью

В коде функций GetOptionPricesVX() мы работаем с 4-мя массивами (pT, pK, PS0, pC). При этом 3 из них используются только для чтения, а один (pC) для записи. Обратим внимание на тот факт, что значения, которые мы записываем в длинный массив pC, в цикле никак не используются. Таким образом, их кэширование вряд ли имеет смысл (массив pC относится к nontemporal data). Для записи напрямую в память, минуя кэш, результатов вычислений, идентифицированных нами как nontemporal data, предназначены так называемые streaming stores, использование которых приводит к уменьшению накладных расходов. В данном конкретном примере мы видим, что цикл имеет огромное число итераций, в нем выполняется сравнительно мало арифметических операций относительно операций с памятью, на Xeon Phi одновременно работают сотни потоков. Совокупность этих факторов приводит к тому, что мы прокачиваем через шину объем данных, сравнимый, а то и превосходящий ее пропускную способность. Использование streaming stores для массива pC позволит нам сэкономить по одной операции чтения, нужной для поддержания системы кэш-памяти в согласованном состоянии, на каждой итерации. Проверим этот факт.

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

#pragma simd
#pragma vector nontemporal
#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, GetOptionPricesV8
};

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

Как видно из табл. 9.17-9.19, времена работы версии 8 по сравнению с версией 7.1 ухудшились на 60 ядрах и улучшились на 120 и 240. Ускорение на 120 и 240 потоках также выросло.

Наконец, для сравнения приведем результаты работы версии 8 на центральном процессоре (16 ядер) и сопроцессоре (120 ядер).

Как видно из таблицы 9.20, время работы параллельной версии с прогревом на сопроцессоре примерно в 4,6 раза лучше, чем на центральном процессоре.

Также отметим, что применения прагмы vector nontemporal почти не ускоряет код на центральном процессоре (сравните времена версий 7.1 и 8), что вполне объяснимо – используется всего 16 потоков, шина не перегружается.

Дополнительные задания