CC BY
28
38
Евразийский Союз Ученых (ЕСУ) # 7 (16), 2015 | ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ
9. Пухлий В.А. Об одном подходе к решению краевых задач математической физики. - Дифференциальные уравнения, 1979, том 15, №11, с.2039-2043.
10. Pfeifer H. Die amplitude des primaren Echos bei der Hahnschen Spin-Echo-Methode. - Annalen der Physik, Band 17. Heft 1, 1955, s.23-27.
11. Das T.P., Saha A.K. Математический анализ экспериментов по спиновому эхо методом Хана. - Phys. Rev., vol.93, 1954, p.749.
12. Waters G.S., Philips G. Новый метод измерения магнитного поля Земли. - Geophys. prospecting., vol.4, 1956, p.1.
13. Zimmerman J.R., Williams D. Использование атома водорода в фотомагнитометрической съемке. - Oil and Gas, vol. 88, 1955.
АНАЛИЗ ЭФФЕКТИВНОСТИ ПРИМЕНЕНИЯ ТЕХНОЛОГИИ CUDA ДЛЯ ЧИСЛЕННОГО РЕШЕНИЯ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ В ЧАСТНЫХ ПРОИЗВОДНЫХ ПО СХЕМЕ КРАНКА НИКОЛСОНА
Вахлаева Клавдия Павловна
Канд. физ. -мат. наук, доцент кафедры информатики и программирования Саратовского государственного
университета им. Н.Г. Чернышевского, г. Саратов Сынкова Дарья Александровна
Бакалавр по направлению «Математическое обеспечение и администрирование информационных систем»,
г. Саратов Фаюстов Денис Сергеевич
Студент 4 курса факультета компьютерных наук и информационных технологий Саратовского государственного
университета им. Н.Г. Чернышевского, г. Саратов
АННОТАЦИЯ
Приведены результаты разработки программы расчета цены конвертируемой облигации, как численного решения дифференциального уравнения в частных производных второго порядка по схеме Кранка-Николсона с использованием метода циклической редукции на гетерогенной вычислительной системе, сочетающей в себе графический и центральный процессоры. Проведен анализ эффективности применения технологии CUDA для реализации схемы Кранка-Николсона.
Ключевые слова: CUDA, GPGPU, схема Кранка-Николсона, метод циклической редукции
Введение
В последнее время активно развивается направление высокопроизводительных вычислений и, в частности, вычисления общего назначения на графических процессорах - GPGPU (General-Purpose computing on GPU) [1]. Возможность переносить интенсивные параллельные вычисления на графический процессор, а последовательные выполнять на центральном процессоре, делает актуальной задачу разработки параллельных программ, согласованных с гетерогенной вычислительной средой, сочетающей в себе графический и центральный процессоры. В рамках данного направления необходимы новые средства разработки GPGPU-приложений. Одним из таких средств является технология CUDA (Compute Unified Device Architecture), разработанная компанией NVIDIA [7].
В данной работе технология CUDA используется для реализации численного метода решения дифференциального уравнения в частных производных второго порядка, описывающего изменение цены конвертируемой
облигации, для задачи из области финансовой математики. Решение дифференциального уравнения выполняется с помощью метода конечных разностей и сводится к многократному решению систем линейных алгебраических уравнений с трехдиагональными матрицами, содержащими одинаковые элементы на каждой диагонали в отдельности.
Целью данной работы является реализация численного метода решения дифференциального уравнения в частных производных для гетерогенной вычислительной среды на базе центрального и графического процессоров.
Постановка задачи и метод решения
Математической моделью рассматриваемой задачи является модель Блэка-Шоулза [5], представляющая собой параболическое уравнение в частных производных относительно стоимости конвертируемой облигации. Жизненный цикл конвертируемой облигации делится на конечное число стадий или временных интервалов в соответствии с рисунком 1.
fo fi f* -ft+i
-------1 I I I I I I I I 1 I I I III-------------------------------►
Ч 4 r - C- -
Рисунок 1. Временная шкала жизни конвертируемой облигации.
Пусть существует m* моментов выплат процентов держатель принимает решение о конвертации (t;,i = в пределах времени исполнения T конвертируемой обли- 0,1,2,..., n; n > m* ). Множество моментов выплат явля-гации (tjj, k = 1,2,..., m* ) и n моментов времени, когда ется подмножеством моментов принятия решений.
Евразийский Союз Ученых (ЕСУ) # 7 (16), 2015 | ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ
39
Пусть Г; - процент, выплачиваемый держателю облигации в момент времени t*:
Г;
ркЛ = tk {0, ti Ф tk
(i = 0,1,.., n; k = 1,2,.., m*),
где rk - заданный параметр для каждого момента времени.
Конвертируемая облигация имеет два важных параметра - номинальная стоимость K и цена конвертации в одну акцию IC.
Задача состоит в том, чтобы определить оптимальную цену конвертируемой облигации на каждой стадии, при которой держатель получает максимальную выгоду в случае конвертации.
Цена конвертируемой облигации Ck = Ck(Vt,t) на каждом временном интервале [tk, tk+1] удовлетворяет дифференциальному уравнению в частных производных:
dCk , 1 d2Ck
dt + 2а Vt
dVt2
dCk
+ (rf-D)V^^-rfCk = 0,
где Vt - курс акций, rf - безрисковая процентная ставка, D - дивиденд по акциям и ст - волатильность.
Данное уравнение предполагает проход от времени исполнения t = T до текущего момента t = t0.
С помощью замены переменных zt = lnVt дифференциальное уравнение с переменными коэффициентами из модели Блэка-Шоулза преобразуется к уравнению с постоянными коэффициентами:
dCk , 1 ^2 d2Ck
1 + - ст'
dt 2 dz2
+ (rf- D -
1,4 dCk
2CT2)dZTrA = a
Для однозначности решения дифференциального уравнения, которым описывается изменение цены конвертируемой облигации, необходимо установить начальное распределение цены конвертируемой облигации на последней стадии:
Cn-1(Vtn,tn) = max{^Vtn, K(1 + rn)},
и на промежуточных стадиях разбиения:
Ck-i(Vtk,tn) = max{-Vtk,Krk + Ck(Vtk,tk)},k e {1, ...,n - 1},
а также задать граничные условия:
lim Ck(Vt, t) = Ke- rf(tk+i t) ( \ ' rk+m+:Le rf(tk+m+i tk+i) + e rf(tn tk+i)),
7t^0 \Z-i^m = 0 )
где t e [tk,tk+i],k = 0, ...,n - 1;
K
lim Ck(Vt, t) = — Vte-D(tk+1-t:),
Vt^+W I c
где t e [tk,tk+i],k = 0, ...,n - 1.
Полученное уравнение решается численно методом конечных разностей по схеме Кранка-Николсона [2, 3, 5]. При построении сеточной области рассматривается диапазон значений стоимости акций от нижней границы
V=0.01 до верхней границы V=10000.0 (z = lnV и "z = lnV ) и для интервала времени от срока исполнения t = T до текущего момента t = 0.
Для каждой временной стадии задается равномерная сетка с помощью двух положительных целых чисел M и J, определяющих разбиение по оси времени и по оси стоимости акций zt в соответствии с рисунком 2.
Рисунок 2. Сетка для построения схемы Кранка-Николсона
На каждой временной стадии M раз решаются системы линейных алгебраических уравнений (СЛАУ) размерности (J - 1) х (J - 1) с трехдиагональной матрицей, содержащей одинаковые элементы на каждой диагонали в отдельности [2-4].
СЛАУ с трехдиагональной матрицей, на каждой диагонали которой стоят одинаковые значения, решается методом циклической редукции, временная сложность которого составляет O(J).
Алгоритм метода циклической редукции может быть эффективно распараллелен для архитектуры CUDA за счет присущего ему внутреннего параллелизма. При
массивно-параллельном вычислении в технологии CUDA каждая отдельная нить будет выполнять часть общего объема вычислений циклической редукции, пропорциональную log2 J [6].
Схема распараллеливания
Реализация алгоритма вычисления цены конвертируемой облигации в виде одно- и многопоточных приложений для исполнения на центральном процессоре подробно рассматривается в [2]. Многопоточные версии выполнены для систем с общей памятью с использованием технологии OpenMP и высокоуровневой библиотеки
40
Евразийский Союз Ученых (ЕСУ) # 7 (16), 2015 | ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ
Intel TBB (Treading Building Blocks) [2]. Алгоритм вычисления цены конвертируемой облигации по схеме Кранка-Николсона на каждом временном интервале может быть разделен на блоки и описан следующим образом:
• вычислить начальные значения цены конвертируемой облигации на последней стадии и на промежуточных стадиях разбиения;
- для каждого момента времени от M — 1 до 0 в интервале [tk, tk+i]:
- получить СЛАУ вида Ax = B, где A - трехдиагональная матрица с одинаковыми элементами на каждой из диагоналей:
Ах =
ао ai
а_1 ао ai
а_ i ао
ai
а i ао ^1
а_ i
aoJ
(J-l)x(J-l)
Ху-1)Х1 = В',
- решить СЛАУ методом циклической редукции;
• определить оптимальное значение цены конвертируемой облигации по вектору х.
Наиболее эффективное распараллеливание приведенного выше алгоритма можно получить, распараллелив следующие функции:
• вычисление начальных значений на последней стадии и на промежуточных стадиях разбиения;
• перемножение матрицы на вектор;
• решение СЛАУ методом циклической редукции. Метод циклической редукции решает трехдиаго-
нальную систему за два этапа: прямой ход и обратный. Суть прямого хода заключается в последовательном исключении из соседних уравнений сначала переменных с нечетными номерами, затем с номерами кратными двум, но не кратными 4, и так далее. Каждый шаг редукции уменьшает число уравнений вдвое, при этом полученная система будет также трехдиагональной и можно повторять процесс до тех пор, пока дальнейшая редукция будет невозможна, после чего начинается обратная подстановка. Операции прямого хода, порождающие редуцированные системы, и операции обратной подстановки, независимы по данным и могут быть распараллелены;
• поиск оптимального значения цены конвертируемой облигации как минимального элемента массива.
Особенности реализации
Коды функций, предназначенных для распараллеливания с применением технологии CUDA, были оформлены в виде ядер CUDA-программы (kernel) и выполнялись на GPU как набор многих одновременно работающих нитей (threads).
Код для всех ядер CUDA-программы разрабатывался с использованием типа double, чтобы не снижать точность расчетов. Трехдиагональные матрицы и правые столбцы систем линейных алгебраических уравнений формируются непосредственно в памяти GPU. Вектор результатов копируется из памяти GPU в память центрального процессора, после чего ресурсы освобождаются.
Далее представлен фрагмент программы, содержащий ядро CydeReductionMethod_step_forwardO для прямого хода циклической редукции.
___global_void CycleReductionMethod_step_forward(double* cuda_f,
double* cuda_a, double* cuda_b, double* cuda_c,
int step, int start, int j, double a_0,
double a_1, double a_2, int n, int elementsNum)
{
int tid = getGlobalIdx_2D_3D();
if (tid >= elementsNum) return;
int blockSize = blockDim.x*blockDim.y*blockDim.z;
if (j == 0 && tid == 0) { cuda_a[0] = a_0; cuda_b[0] = a_1; cuda_c[0] = a_2;
}
__syncthreads();
double alpha = -cuda_a[j] / cuda_b[j]; double beta = -cuda_c[j] / cuda_b[j]; if (tid == 0) {
cuda_a[j + 1] = alpha * cuda_a[j];
cuda_b[j + 1] = cuda_b[j] + 2 * alpha * cuda_c[j];
cuda_c[j + 1] = beta * cuda_c[j];
}
for (int i = tid; i < elementsNum; i += blockSize) { int k = start * (i + 1);
cuda_f[k] = alpha * cuda_f[k - step] + cuda_f[k] + beta * cuda_f[k + step]; }
}
Для определения индекса нити вызывается функция getGlobalIdx_2D_3D(). В случае если индекс нити выходит за пределы результирующего вектора, данная нить не выполняет никаких действий, для этого в коде предназначена специальная проверка.
Результаты вычислительных экспериментов Эксперименты проводились при следующих значениях входных параметров: число временных стадий n =100, число разбиений по оси времени на каждой стадии M =100, число разбиений по оси стоимости акций J,
Евразийский Союз Ученых (ЕСУ) # 7 (16), 2015 | ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ
41
задающее размерность трехдиагональной матрицы (J — 1) х (J — 1) в системе линейных алгебраических уравнений, кратно степени двойки J = 2q, q £ И.
Вычислительные эксперименты проводились с использованием инфраструктуры, представленной в таблице 1. Для CUDA-программы измерялось время расчета с учетом всех накладных расходов на запуск нитей.
Основная вычислительная сложность схемы Кранка-Николсона определяется трудоемкостью методов решения СЛАУ. Результаты исследования зависимости времени выполнения программ, применяющих метод циклической редукции, на CPU и с использованием GPU от размерности трехдиагональной матрицы приведены в таблице 2.
Таблица 1
Тестовая инфраструктура
Центральный процессор (CPU) Intel(R) Core(TM) i5-3450 CPU @ 3.10 GHz
Память 8 ГБ
Операционная система Microsoft Windows 7 Ultimate, 64-bit (SP1)
Среда разработки Microsoft Visual Studio 2012
GPU NVIDIA GeForce GTX 670
Число CUDA-ядер 1344
Память GPU 2048 МБ GDDR5
Шина PCI-E 16x 3.0
Скорость передачи данных памяти 6008 МГц
Версия драйвера CUDA 347.62
Таблица 2
Время выполнения программ, применяющих метод циклической редукции, на CPU и с использованием GPU___
J Время работы реализации на CPU (сек) Время работы реализации с использованием GPU (сек)
128 0.017 1.127
512 0.039 1.030
1024 0.113 1.323
2048 0.130 1.422
4096 0.255 1.673
16384 1.056 2.654
32768 2.574 3.805
65536 5.647 5.878
131072 10.503 10.002
262144 28.954 18.331
524288 52.856 33.426
1048576 123.336 66.078
Как видно, начиная с J =131072, которому соответствует q = 17, реализация метода циклической редукции с применением CUDA показывает лучшее время решения. Заключение
В рамках данной работы на базе центрального и графического процессоров разработана CUDA-программа, реализующая численный метод решения дифференциального уравнения в частных производных по схеме Кранка-Николсона с применением метода циклической редукции.
Проведенное исследование времени выполнения программы показало эффективность применения технологии CUDA для ускорения решения СЛАУ с трехдиагональными матрицами методом циклической редукции при размерностях матрицы больших J = 217.
Список литературы
1. Боресков А. В., Харламов А. А., Марковский Н. Д. и др. Параллельные вычисления на GPU. Архитектура и программная модель CUDA / А. В. Боресков и др. Предисл.: В. А. Садовничий. - М.: Изд-во Московского университета, 2012. 336 с.
2. Гергель В. П., Баркалов К. А., Мееров И. Б. и др. Параллельные вычисления: технологии и численные методы: учебное пособие. В 4 т. Т. 4. - Н. Новгород: Изд-во Нижегородского госуниверситета, 2013. 369 с.
3. Мееров И. Б., Никонов А. С., Русаков А. В., Шишков А. В. Реализация одного алгоритма нахождения цены конвертируемой облигации. // Технологии
Microsoft в теории и практике программирования. Материалы конференции / Под ред. проф. Р. Г. Стронгина. - Н. Новгород: Изд-во Нижегородского госуниверситета, 2008. - С. 236-241.
4. Мееров И. Б., Никонов А. С., Русаков А. В., Шишков А. В. Параллельная реализация одного алгоритма нахождения цены конвертируемой облигации для систем с общей памятью // Технологии Microsoft в теории и практике программирования. Материалы конференции / Под ред. проф. В. П. Гер-геля. - Н. Новгород: Изд-во Нижегородского госуниверситета, 2009. - С. 287-292.
5. Косяков М. С. Применение технологии CUDA для ускорения расчета цен опционов европейского типа сеточным методом / М. С. Косяков, Д. Н. Шинка-рук, А. В Торопов, Ю. А. Шполянский // Известия высших учебных заведений. Приборостроение, 2012. - Т.55, вып.10. - С. 13-19.
6. Косяков М. С. Анализ эффективности применения технологии CUDA для решения систем линейных уравнений с трехдиагональными матрицами в задачах расчета цен опционов / М. С. Косяков, Д. Н. Шинкарук, А. В Торопов, Ю. А. Шполянский // Известия высших учебных заведений. Приборостроение, 2012. - Т.55, вып.10. - С. 20-25.
7. Параллельные вычисления CUDA [Электронный
ресурс] URL: http://www.nvidia.ru/object/cuda-
parallel-computing-ru.html. (дата обращения: 24.06.2015).
