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

С.И. Бастраков, Р.В. Донченко, И.Б. Мееров, А.Н. Половинкин

Нижегородский государственный университет им. Н.И. Лобачевского

Статья опубликована в журнале "Научно-технический вестник Санкт-Петербургского государственного университета информационных технологий, механики и оптики", 2010, № 5(69). – C. 95–100.

Рассматриваются вопросы оптимизации вычислений в прикладных расчетных программах на языке С. Основной целью оптимизации в данном случае является эффективная утилизация имеющихся вычислительных ресурсов. Изложение проводится на примере решения задачи определения справедливой цены опциона европейского типа. Эффективность предложенного подхода подтверждена вычислительными экспериментами.

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

Введение

Оптимизация вычислений для уменьшения времени работы прикладных программ расчетного характера – актуальная задача современной науки и техники. Решением этой задачи активно занимаются ученые и инженеры по всему миру [1–3]. Существенной проблемой является то, что прикладные программисты часто не вполне владеют приемами оптимизации ПО, а системные программисты плохо разбираются в предметных областях, что затрудняет процесс оптимизации. В настоящей работе авторы ставят целью продемонстрировать методику уменьшения времени вычислений, не выходя за рамки многоядерных архитектур и языка программирования высокого уровня, оставаясь в пределах доступного исследователю программно-аппаратного окружения. Изложение ведется на примере языка С, подавляющая часть описанных приемов с тем же успехом может быть применена к программам на языке Fortran (чаще всего именно С или Fortran используются при реализации численных методов). Для иллюстрации приемов оптимизации рассматривается задача из области финансовой математики – определение справедливой цены опциона европейского типа. Процесс оптимизации вычислений иллюстрируется фрагментами программного кода, полученный результат подтверждается вычислительными экспериментами.

Постановка задачи и метод решения

Цена европейского опциона, основанного на облигации, может быть вычислена по следующим формулам [4, 5]:
где t – текущий момент времени; s – срок исполнения опциона; P(s, T) – цена в момент времени s облигации со сроком погашения T; К – цена исполнения опциона; r(t) – краткосрочная безрисковая ставка в момент времени t. Искомая величина – цена опциона C = C(0, s) . (Дополнительная информация об использованных понятиях может быть получена из работы [4]). Математическая постановка задачи, метод решения и вычислительная схема подробно описаны в [5]. Приведем краткое описание схемы проведения вычислений на концептуальном уровне.
Задача решается методом Монте-Карло.

Шаг 1. Инициализация параметров начальными значениями.
Шаг 2. Вычисления для определения цены опциона методом Монте-Карло.

В цикле от 0 до количества повторений в методе Монте-Карло выполнить:
Шаг 2.1. Пусть k – число шагов в разностной схеме решения стохастического дифференциального уравнения (СДУ), описывающего поведение ставки r(t) от момента заключения опционного контракта (t = 0) до истечения срока его исполнения (t = s). На данном шаге генерируется k нормально распределенных псевдослучайных чисел.
Шаг 2.2. Решение СДУ методом Эйлера. Используются результаты шага 2.1. Вычисляются k значений процентной ставки r(t) на рассматриваемом интервале времени.
Шаг 2.3. Определение цены облигации P(s,T) .
Шаг 3. Усреднение полученных данных: вычисление цены опциона C = C(0, s) по формуле (1).

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

В рамках данной работы для демонстрации приемов оптимизации по скорости создано несколько версий ПО. Каждая следующая версия улучшает предыдущую, сокращая время работы. Реализация выполнена на языке программирования ANSI C. Кроме того, разработана отдельная версия с использованием NVIDIA CUDA, оптимизированная под графические процессоры NVIDIA.

Параметры задачи и окружение. В тестовой задаче выполнялось определение цены 30 опционов с разным значением цены исполнения, использовалось 50000 повторений в методе Монте-Карло и вещественная арифметика двойной точности. Количество опционов было выбрано из тех соображений, чтобы время работы базовой версии составляло порядка 300 секунд. Вычислительный эксперимент производился с использованием следующего программно-аппаратного окружения: 2x Intel Xeon e5520 (4 ядра, 2.27 Ггц), 16 Гб ОЗУ, NVIDIA Tesla C1060 (240 ядер, 1.30 ггц), openSUSE Linux 11.1 (x86_64), kernel 2.6.27.29-0.1-default, glibc 2.9, gcc 4.3.2, Intel С++ Compiler for Linux 11, Intel Math Kernel Library 10.2.2, NVIDIA CUDA Toolkit 2.2.

Базовая версия – элементы реализации. Базовая версия представляет собой неоптимизированную реализацию описанной ранее вычислительной схемы. Для генерации псевдослучайных чисел на шаге 2.1 применяется второй вариант преобразования Бокса–Мюллера [6]; равномерно распределенные на (0,1) числа получаются при помощи генератора MCG59 (реализован по описанию в [7]). Общее время работы исходной версии составляет 327 с.

Базовая версия – профилировка и оптимизирующий компилятор.
Первоначальным этапом оптимизации отлаженного работоспособного программного кода является его профилировка. В данной работе использовалась консольная версия Intel Vtune Performance Analyzer. Рассмотрим результаты профилирования базовой версии (компилятор gcc, ключ оптимизации – O2).

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

Описание функцииИмя функцииNTSelfTime
Моделирование r(t) (Шаг 2.2)pow442 34285,80%
Генерация псевдослучайных чисел
(Шаг 2.1)
log18 5413,60%
generateGaussian11 9852,32%
sincos3 9000,76%
Вспомогательные вычисленияrun17 4593,39%

 

Из результатов профилирования видно, что основную вычислительную нагрузку составляют моделирование процентной ставки (~86%) и генерация псевдослучайных чисел (~7%). Прежде чем приступить к программной оптимизации, попробуем добиться прироста производительности при помощи использования оптимизирующих компиляторов. В рамках данной работы используется один из лучших компиляторов C/C++ – Intel C++ Compiler for Linux 11.

Без каких-либо изменений кода программы общее время работы составляет 99,08 с, что в 3,3 раза меньше, чем для версии, скомпилированной gcc.

Оптимизированная версия – предварительные вычисления. Выполним профилировку базовой версии, откомпилированной при помощи Intel C++ Compiler.

Описание функцииИмя функцииNT (gcc)NT (Intel)SelfTime
Моделирование r(t) (Шаг 2.2)pow.L442 342114 99373,40%
Генерация псевдослучайных чисел
(Шаг 2.1)
log.L18 5418 4965,42%
__libm_sse2_sincos3 9004 1672,66%
Вспомогательные вычисленияrun17 459 + 11 98526 86917,15%

Как видно из профиля, время, затраченное на вычисление pow, уменьшилось почти в 4 раза, время вычисления log также уменьшилось почти в 2 раза. Можно заметить, что из списка «исчезла» функция generateGaussian, при этом увеличилось время работы функции run. Это вызвано тем, что функция generateGaussian была встроена компилятором, что привело к уменьшению времени работы (17 459 + 11 985 > 26 869).

Рассмотрим фрагмент кода, выполнение которого занимает 73,4% времени работы программы.

r[0] = f0t; 
phi1 = 0; 
phi2 = 0; 
chi1 = 0;
for (i = 0; i < numSteps; i++) 
{
    dr = (kappa * (f0t - r[i]) + phi1 + phi2 + kappa * chi1) * dt +
                       sigma1 * pow(r[i], alpha) * z1[i] + 
                       sigma2 * pow(r[i], beta) * z2[i];
    dphi1 = sigma1 * sigma1 * pow(r[i], 2 * alpha) * dt;
    dphi2 = (sigma2 * sigma2 * pow(r[i], 2 * beta) - 2 * kappa * phi2) * dt;
    dchi1 = phi1 * dt + sigma1 * pow(r[i], alpha) * z1[i];
    r[i + 1] = r[i] + dr; 
    phi1 += dphi1; 
    phi2 += dphi2; 
    chi1 += dchi1;
}

Несложный анализ приведенного фрагмента показывает, что количество вызовов «тяжеловесной» функции pow избыточно: достаточно вычислить величины pow(r[i], alpha), pow(r[i], beta) один раз и затем использовать вычисленные значения, а также их квадраты вместо pow(r[i], 2*alpha), pow(r[i], 2*beta). Кроме того, во всех вхождениях значение pow(r[i], alpha) умножается на sigma1, а pow(r[i], beta) – на sigma2, поэтому целесообразно вычислить данные произведения заранее. Ниже приведен усовершенствованный участок кода:

r[0] = f0t; 
phi1 = 0; 
phi2 = 0; 
chi1 = 0;
for (i = 0; i < numSteps; i++) 
{
    double c0, c1;
    c0 = sigma1 * pow(r[i], alpha); 
    c1 = sigma2 * pow(r[i], beta);
    r[i+1] = r[i] + (kappa * (f0t - r[i]) + phi1 + phi2 + kappa * chi1) * dt +
                c0 * z1[i] + c1 * z2[i];
    chi1 += phi1 * dt + c0 * z1[i];
    phi1 += c0 * c0 * dt; 
    phi2 += (c1 * c1 - 2 * kappa * phi2) * dt;
}

Общее время работы составляет 64,86 с, что в 1,5 раза быстрее предыдущей версии.

Оптимизированная версия – использование инструкций SIMD. Выполним профилирование разработанной версии программной реализации (столбец Intel+) для поиска путей дальнейшей оптимизации.

Описание функцииИмя функцииNT (Intel)NT (Intel+)SelfTime
Моделирование r(t) (Шаг 2.2)pow.L114 99367 36065,58%
Генерация псевдослучайных чисел
(Шаг 2.1)
log.L8 4968 4428,22%
generateGaussianвстроена8 8138,58%
__libm_sse2_sincos4 1674 3184,20%
Вспомогательные вычисленияrun26 86912 27311,95%

В результате уменьшения количества вычислений суммарное время работы pow (NT Intel+) уменьшилось почти в 2 раза по сравнению с предыдущей версией (NT Intel). Время работы остальных функций практически не изменилось, существенное ускорение функции run является мнимым – компилятор решил не встраивать функцию generateGaussian. Исходя из результатов анализа текущей версии кода, на следующем этапе оптимизации было решено использовать векторизацию. Необходимо отметить, что многие циклы векторизуются компилятором Intel автоматически при включении опции -O2. При этом эвристика компилятора предсказывает, ожидается ли выигрыш от векторизации. Для устранения возможных ошибок прогноза существует директива #pragma vector always, действующая для ближайшего цикла.

Рассмотрим приведенный выше фрагмент кода. Очевидно, что цикл не может быть векторизован из-за большого количества зависимостей внутри итерации и между итерациями. Тем не менее, возможности для оптимизации все еще остаются.

Попробуем «вручную» векторизовать вызовы pow – самой трудоемкой функции в данной программе. Для этого создадим массивы из двух элементов и цикл из двух итераций. Иногда компилятор откажется от векторизации такого короткого цикла, поэтому воспользуемся директивой #pragma vector always. Ответ на вопрос о том, был ли векторизован цикл, можно получить из отчета компилятора о векторизации.

r[0] = f0t; 
phi1 = 0; 
phi2 = 0; 
chi1 = 0;
for (i = 0; i < numSteps; i++) 
{
    double c[2], sigma[2] = {sigma1, sigma2}, power[2] = {alpha, beta};
    int k;
#pragma vector always
    for (k = 0; k < 2; k++)
        c[k] = sigma[k] * pow(r[i], power[k]);
…
}

В результате изменений время работы данной версии (Intel++) составляет 54,8 с, что на 20% меньше, чем время работы предыдущей версии.

Оптимизированная версия – использование библиотек. Выполним профилировку для поиска возможностей дальнейшей оптимизации.

Описание функцииИмя функцииNT (Intel+)NT (Intel++)SelfTime
Моделирование r(t) (Шаг 2.2)__svml_pow2.A67 36053 08261,20%
Генерация псевдослучайных чисел
(Шаг 2.1)
log.L8 4428 55510,06%
generateGaussian8 8138 7269,86%
__libm_sse2_sincos4 3184 3134,22%
Вспомогательные вычисленияrun12 27311 13112,83%

Как видно из профиля, используется векторная функция вычисления степени __svml_pow2.A, ее время работы почти на 20% меньше, чем в предыдущей версии.

Практически все предыдущие оптимизации были направлены на уменьшение времени работы функции возведения в степень pow. В итоге генерация случайных чисел стала занимать существенное время от общего – около 25%, что является поводом для оптимизации и этого этапа вычислений, для чего используется оптимизированная под процессоры Intel библиотека Intel Math Kernel Library (MKL), в состав которой, в частности, входят генераторы псевдослучайных чисел. Так, для генерации нормально распределенных псевдослучайных чисел может быть использована функция vdRngGaussian. Рассмотрим профиль приложения после изменений.

Описание функцииИмя функцииNT (Intel+)NT (Intel_MKL)SelfTime
Моделирование r(t) (Шаг 2.2)__svml_pow2.A53 08253 73670,59%
Генерация псевдослучайных чисел
(Шаг 2.1)
_vmldSinCos_2621 59410 42013,69%
_vdRngGaussianBoxMuller2
_vmldLn_26
_vmldSqrt_26
__vsldBRngMCG59
Вспомогательные вычисленияrun11 13110 82914,23%

Как видно из профиля, время генерации случайных чисел уменьшилось примерно в 2 раза. Общее время работы составляет 48,35 с, что на 13% меньше, чем ранее.

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

#pragma omp parallel for
for (i = 0; i < numOptions; i++) {
    …
    getAveragePrice_MKL(&data_MKL_openMP[i]);
}

Время работы программы на 8 вычислительных ядрах составило 6,98 с, что в 46,87 раз меньше, чем в базовой версии.

Ниже приведены диаграммы, иллюстрирующие ускорение для параллельной версии (рис. 1) и результаты оптимизации (рис. 2).

Реализация на графическом процессоре. В настоящий момент в научной и научно-популярной литературе большое внимание уделяется программированию на графических процессорах (GPU). Основной акцент в таких публикациях часто делается на пиковую производительность GPU, многократно превосходящую пиковую производительность центрального процессора (CPU). При этом решаемые задачи чаще всего характеризуются большим количеством однотипных вычислений над однотипными данными. При таких условиях (идеальное распараллеливание) GPU, как правило, демонстрируют свою эффективность. В данной работе мы попробовали перенести имеющийся (достаточно хорошо распараллеливающийся) код на GPU с использованием NVIDIA CUDA, уделяя основное внимание производительности.

В выполненной реализации каждому опциону соответствует блок потоков (каждый из которых назначается на свой мультипроцессор): при этом каждый поток внутри блока выполняет свою часть итераций метода Монте-Карло. Такое распределение работы позволяет разместить необходимые для вычисления опциона общие данные в так называемой разделяемой (общей) памяти, имеющей низкую латентность. Роль центрального процессора в данном случае ограничивается подбором конфигурации запуска вычислительного ядра и копированием входных данных в память GPU, а выходных – обратно. Как и в случае реализации для CPU, использовались вычисления двойной точности. Время работы приложения на GPU составило 6,43 с, что сопоставимо с временем работы на CPU – 6,98 с (в обоих случаях рассмотрены параллельные реализации).

Конечно, одной задачи недостаточно, чтобы делать далеко идущие выводы. Однако идея переноса расчетных программ на GPU как способ быстрого и «бесплатного» увеличения их производительности представляется как минимум дискуссионной из-за следующих проблем. Оптимизация по скорости на GPU требует существенного опыта подобной работы и дополнительной подготовки: написание программ, которые эффективно используют преимущества архитектуры графического процессора без глубокого знания данной архитектуры является на сегодняшний день трудновыполнимым. Так, при наличии навыков работы с NVIDIA CUDA и работающего на CPU исходного кода перенос сравнительно несложного приложения (около 740 строк, включая исходные и заголовочные файлы с функциями вычислений, тестирования, ввода-вывода, подробными комментариями) потребовал времени, сравнимого со временем написания данного кода с нуля. Дополнительно положение осложняется отсутствием возможностей для полноценной отладки (в настоящий момент компания NVidia выпустила отладчик с ограниченной функциональностью). Архитектура GPU накладывает серьезные ограничения на схему распараллеливания; кроме того, чем больше сложность вычислительной части, тем менее эффективно удается использоватьвычислительную мощь устройства из-за относительной нехватки регистров и общей памяти. Таким образом, представляется, что эффективное использование GPU + NVIDIA CUDA возможно далеко не во всех задачах и требует наличия в коллективе опытных системных программистов.

Заключение

В работе на примере конкретной задачи исследуются вопросы оптимизации времени работы прикладного ПО для современных вычислительных систем с общей памятью. Суть предлагаемого подхода состоит в эффективном использовании алгоритмической, компиляторной, программной оптимизации, а также специализированных средств – оптимизирующих компиляторов, профилировщиков, математических библиотек. В результате применения рассмотренных методов и средств оптимизации время работы последовательной версии программы сократилось в 6,77 раза, а после распараллеливания – в 46,87 раза при выполнении на 8 вычислительных ядрах.

Работа выполнена в лаборатории «Информационные технологии» (ITLab) ВМК ННГУ (www.itlab.unn.ru)

Литература

  1. Gerber R., Bik A., Smith K., Tian X. The Software Optimization Cookbook Second Edition. High Performance Recipes for IA 32 Platforms // Intel Press, 2005.
  2. Fog A. Optimizing software in C++: An optimization guide for Windows, Linux and Mac platforms. – 2009 [Электронный ресурс]. – Режим доступа: http://www.agner.org/optimize/optimizing_cpp.pdf, своб.
  3. Касперски К. Техника оптимизации программ. Эффективное использование памяти. – СПб: BHV,2003. – 464 c.
  4. Халл Дж. Опционы, фьючерсы и другие производные финансовые инструменты. – М.: Вильямс, 2007. – 1056 с.
  5. Inui K., Kijima M. A Markovian framework in multi-factor HJM models // Journal of financial and quantitative analysis. – 1998. – V. 33. – № 3. – Р. 423–440.
  6. Box G.E.P. and M.E. Muller. A Note on the Generation of Random Normal Deviates // Ann. Math. Stat. – 1958. – V. 28. – Р. 610–611.
  7. Intel MKL. Vector Statistical Library Notes, revision -019, section 8.4.4.

Additional Resources

Intel AI Developer Program 

Лаборатория Intel в ННГУ

For more complete information about compiler optimizations, see our Optimization Notice.