Алгоритмы с «Длинными» векторами решения сеточных уравнений явных разностных схем Текст научной статьи по специальности «Математика»

АЛГОРИТМЫ С «ДЛИННЫМИ» ВЕКТОРАМИ РЕШЕНИЯ СЕТОЧНЫХ УРАВНЕНИЙ ЯВНЫХ РАЗНОСТНЫХ СХЕМ

Воротникова Д.Г., Головашкин Д.Л.

Институт систем обработки изображений РАН, Самарский государственный аэрокосмический университет имени академика С.П. Королёва (национальный исследовательский университет) (СГАУ)

Аннотация

Предложены два варианта алгоритмов с «длинными» векторами для решения сеточных уравнений явных разностных схем, позволяющих задействовать одновременно максимальное количество ядер СЦЭА даже для сеточной области небольшой размерности. На примерах разностного решения уравнений теплопроводности и Максвелла продемонстрирована эффективная реализация предложенного подхода. Произведено сравнение предложенных авторами алгоритмов, реализованных при помощи библиотеки CUBLAS, со свободно распространяемыми пакетами В-СЛЬМ и ОрепСштеп!

Ключевые слова: уравнение теплопроводности, уравнения Максвелла, векторные алго-

ритмы, разностные схемы, CUBLAS.

Введение

Ведущей современной тенденцией развития вычислительной техники следует признать наращивание производительности не за счёт увеличения тактовой частоты центрального процессора (проблема «кремниевого тупика» [1]), а посредством новых архитектурных решений (многопоточность, векторные сопроцессоры, многоядерность и др.). Наиболее популярным направлением в этой области авторы настоящей работы признают организацию вычислений общего назначения на графических процессорах (GPGPU [2]), относящихся к модели SIMD по классификации Флинна [3]. Реализация основных методов вычислительной математики в рамках упомянутой модели обсуждалась достаточно активно на протяжении последней четверти века [4, 5], хотя данная проблема и не находилась в центре внимания в силу ограниченной доступности соответствующего аппаратного обеспечения (суперкомпьютеры серии Cray [6]). В настоящее время в связи с широким распространением вычислений на графических процессорах (бытовых, таких как видеокарты серии NVIDIA GeForce [7], и профессиональных -NVIDIA Tesla [8]) создание специализированных алгоритмов и программных продуктов с учётом особенностей модели SIMD становится весьма актуальным.

Интересуясь теорией разностных схем, авторы отмечают традиционное предпочтение неявных схем явным, особенно заметное в отечественной литературе [9], [10] и объясняемое худшей устойчивостью последних в большинстве случаев. Действительно, необходимость налагать более густые сеточные области для сохранения устойчивости в случае решения параболических дифференциальных уравнений приводит к росту числа операций при вычислениях по явным схемам, несмотря на их простоту по сравнению с неявными. Однако с недавним появлением доступной аппаратной базы для организации векторных вычислений указанный недостаток нивелируется эффективностью графических процессоров в силу простоты векторизации вычислений по явным схемам и сложностью по неявным. Последнее утверждение относится к прямым методам решения сеточных

уравнений неявных схем (прогонки, циклической редукции и др.) [5] и значительно ослабляется при рассмотрении итерационных методов [4]. К сожалению, применение последних подразумевает известный уровень квалификации исследователя (необходимо удачно выбрать сам метод, начальное приближение, задать точность), что ограничивает использование неявных схем в инженерной практике.

К настоящему времени известно множество коммерческих и свободно распространяемых пакетов, реализующих вычисления по явным схемам на графических процессорах. Авторы представленного исследования пользуются разработками ОрепСшггеШ: [11] и В-СЛЬМ [12] для численного решения уравнений теплопроводности и Максвелла. Изучение открытого кода данных пакетов, основанного на «коротковекторном» представлении (по классификации Ортеги [4]) привело к идее повысить эффективность вычислений переходом к «длинным» векторам, позволяющим задействовать ресурсы видеокарты более полно.

Возможность упомянутого перехода ранее продемонстрирована на примере двух «длинновекторных» алгоритмов метода Якоби для решения сеточных уравнений простейшей неявной схемы для уравнения Пуассона, представленных в работе [4]. Развитие этой методики организации вычислений в случае явных разностных схем (для уравнений теплопроводности и Максвелла) представлено далее.

1. Векторные алгоритмы для уравнения теплопроводности

В качестве первого примера для исследования методов, основанных на использовании «длинных» векторов, рассмотрим разностное решение двумерного однородного линейного нестационарного уравнения теплопроводности:

dU (d2U d2U

-= a I —- + ——

dt { dx2 8y2

(1)

где U(x, y, t) - распределение температуры в пространстве и времени, x e[0,X ] ,y е[0,У ] ,t e[0,T ],

a - коэффициент температуропроводности. Пусть граничные условия соответствуют условиям Дирихле (температура на границах области равна нулю), а начальное условие имеет вид U(x,y,0) = ф(x,y). Простейшая явная разностная схема [9] для этой задачи на сеточной области

= {(X >yj' tk): X =ih yj = jh (К =hy=h),

tk = к x, i = 0: N +1, j = 0: N +1, к = ii^J имеет вид:

Uj = a(uj + Ui+i,j -4Ukj +U^j +Ui!j+i) + Ukj, (2)

где a = a (x / h2).

Далее нам удобнее будет представлять (2) в следующем итерационном виде, следуя рекомендациям А.А. Самарского [9]:

Uk+1 = a AUk + Uk, (3)

(T B B T2

где A =

B2

Bn

TN-1 BN-1

B„

T

± лг

т тк ^

и - вектор, полученный из двумерного массива

т тк

и1 чередованием строк

ик = (ик ик ик ик ик V

В силу равенства коэффициентов при ик++х , ик.-1 , и— , матрица А - симметричная, Б. - диаго-

U,

i_1j ' i_1j

нальная, хранящая коэффициенты при Ukj_ и Ukj+1 , а Tj - трёхдиагональная, в которой хранятся коэффициенты при Ukj (на главной диагонали), Uk_1j и Uk+1j.

1.1. Алгоритм с «короткими» векторами

Наиболее очевиден вариант алгоритма с использованием «коротких» векторов, представляющих собой строки или столбцы сеточной области. Запишем алгоритмическую реализацию выражения (2) в форме следующего псевдокода для последующего анализа: for j = 2: N +1

for i = 2: N +1

Unx, j) = a (иы, j _1) + иш (i +1 j) +

+ U hsl (i _ 1, j) + U lasl (i, j +1)) + + (1 _ 4a)U lasl (i, j);

end end

U at = Unexl.

В алгоритме U last и Unext представляют собой двумерные массивы размера (N+2)*(N+2) элементов,

содержащие значения сеточных функций. В большинстве современных векторных пакетов, например в OpenCurrent, используется именно такой способ представления данных. В такой форме записи алгоритма заметно, что он может быть векторизован по столбцам (или строкам) массива U¡ast, то есть операции сложения, умножения могут быть выполнены не поэлементно, а по отношению к столбцам указанных массивов с использованием нотации из [5] следующим образом.

for j = 1: N Unex, (1: N, j) = «(Uas, (0: N -1, j) +

+ Ulcat (2: N +1, j) + Ulcat (1: N, j +1)+ (4) + U lcat (1: N, j -1)) + (1 - 4a) Ua, (1: N, j);

end

Uiast = Umx, .

В (4) запись U iast (1: N,j) означает обращение к j-му столбцу массива Ulast без первого и последнего элементов, содержащих граничные значения, Ulast (2: N+1,j) -это тот же столбец, сдвинутый на одну позицию вниз.

Главным недостатком данного алгоритма при реализации на видеокарте является загрузка не всех микропроцессоров при пересчёте сеточной функции для небольших значений N. К примеру, чтобы загрузить полностью видеокарту последнего поколения NVIDIA GeForce GTX TITAN Z [7], нам необходима матрица, строки или столбцы которой имеют размер не менее 5760 элементов.

1.2. Первый алгоритм с «длинными» векторами

Следуя (4), рассмотрим следующий способ хранения массива, содержащего значения сеточной функции. Пусть V - это одномерный массив длины (N+2)2, содержащий значения сеточной функции, хранящиеся по строкам, показанный на рис. 1. В отличие от ранее введённого U в V входят граничные значения.

Сумма используется дважды

„Ш)________ГФП-^шц т

[ Г(АГ+3)_____ V{N+4) Lr(jV+5)l V(N+6)

V(2N+5) ~ JV(2N+6)KV(2N+7) V(2N+8f] V(3N+7) T(3N+S) ~[VJ3N+£^V(3N+10)

V(N+2)) V(2(N+2)) V(3(2N+2)) V(4(3N+2))

V((N+2)(N+2))

ГР+2)(ЛГ+/)Н) У((Ы+2ХМ+1у-2) ......

Рис. 1. Схема хранения в виде одномерного массива и попарные суммы

Разностной схеме для уравнения теплопроводности (2) при вычислении лапласиана соответствует дифференциальный шаблон «крест». На рис. 1 обозначены светло- и тёмно-серым выделением элементы вектора V, участвующие в расчёте новых значений ¥(N+4) и ¥^N+7). Очевидно, что попарная сумма элементов иш(3, 2) и иш(2, 3), соответствующих ¥(N+5) и ¥^N+6), показанная на рис. 1 двунаправленной стрелкой, подсчитывается дважды в алгоритме с короткими векторами, так как используется для нахождения двух элементов сеточных функций. Благодаря использованию «длинных» векторов вместо строк или столбцов матрицы, мы можем сократить общее число операций,

необходимое для расчёта значений сеточной функции, путём введения вспомогательного вектора Т размерности (Ж+2)(Ж+1)+1, который будет хранить в себе результаты попарного сложения.

Таким образом, алгоритм подсчёта сеточной функции на каждом 1-м временном шаге можно записать в следующем виде:

1. Заполнение вектора Т попарными суммами:

Т (2:( N +1)( Ы + 2)) ^ (2:( N +1)( Ы + 2)) +

-V(N + 3:(N + 2)2 -1).

2. Подсчёт значений для следующего временного слоя:

V (N + 4: (N +1)( N + 2)-1) =

= а(т (2: N (N + 2)-1) + Т (N + 5: N (N + 2)-1)) +

+V (N + 4: (N +1)( N + 2)-1).

3. Восстановление граничных значений.

Следующий временной слой пишется поверх предыдущего.

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

Переход к работе с «длинными» векторами позволил перейти от цикла, повторяющегося N раз (количество столбцов матрицы), в котором выполняется 3 операции векторного сложения и одна операция saxpy с векторами длины N, к двум операциям сложения векторов длины (N +1)(N + 2) - 2 и одному выполнению saxpy для векторов длины N(N + 2) -1, что (как будет показано далее) приводит к более полному использованию аппаратных возможностей видеокарты и, как следствие, ускорению вычислений.

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

1.3. Второй алгоритм с «длинными» векторами

Второй вариант - векторизация с использованием gaxpy, основной операции в (3). Главная диагональ матрицы А в (3) соответствует центральному узлу дифференциального шаблона. Коэффициент при нём состоит из двух слагаемых: коэффициента при сеточной функции и*} из конечной разности по пространству и единицы из конечной разности по времени. Первая и минус первая диагонали соответствуют коэффициентам при Ц+1 (первая) и Ц-1 (минус первая). Главная диагональ не содержит нулевых элементов в отличие от первой и минус первой, где бло-

ки неплотно примыкают друг к другу (по одному нулю в каждой диагонали):

Tj =

- 4а 1 1 1 - 4а 1

1 1 - 4а

(5)

Две другие блочные диагонали (N-1 и -N+1) формируются коэффициентами при Ukj+1 и Ukj-l , расположенными на главных диагоналях указанных блоков. Таким образом, мы получаем матрицу с пятью ненулевыми диагоналями: главная длины N х N, первая и минус первая диагонали длины N х N -1 и N-1 и -N+1 диагонали длины N х (N -1).

В работе [13] авторами настоящей статьи предложены векторные алгоритмы умножения трёхдиаго-нальной матрицы на столбец при помощи GPU. Основная операция, используемая там, - покомпонентное умножение векторов. В общем виде умножение трёхдиагональной матрицы, главная диагональ которой имеет размер N х N, на вектор z = Ax, может быть записано следующим образом в нотации из [13].

z = Diag0. * X % умножение главной диагонали на вектор х.

z (1: N2-1) = z (1: N2-1) + Diag, (1: N2-1). * x (2: N2);

% работа с первой диагональю.

z (2: N2 ) = z (2: N2) + Diag-l (2: N2). * x (1: N2 -1);

% работа с минус первой диагональю.

m = N-1; t = (m- 1)N2 -m(m-1)/2; z (1: N2 - m ) = z (1: N2 - m) +

+ DiagN-1 (t +1: t + N2 - m). * x (m +1: N2);

% работа с N-1 диагональю.

z (m +1: N2 ) = z (m +1: N2) +

+ Diag-N+l (t + 1: t + N2 - m). * x (1: N2 - m);

% работа с -N+1 диагональю.

Здесь операция « .* » обозначает покомпонентное умножение векторов.

В данном алгоритме мы также работаем с векторами наибольшей длины, а не столбцами (строками). В процессе работы алгоритма мы дважды выполняем операцию saxpy для векторов длины Nx(N-1) и длины NxN-1 и один раз производится покомпонентное умножение векторов длины NxN. В данном варианте алгоритма с «длинными» векторами не происходит затирание, а, следовательно, и восстановление граничных значений.

2. Алгоритм FDTD-метода с «длинными» векторами

Рассмотрим двумерный вариант FDTD-метода, широко применяемого [15] в дифракционной нанофо-тонике, вычислительной электродинамике и в других

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

дну

1<>

ее,

dEx dHz dEx dy

dt dz "dt dEx _ dHz dHy

(6)

dt

dy dz

где Ex, Hy, Hz - проекции векторов напряжённо-стей электромагнитного поля, ц0 , е0 - магнитная и электрическая константы, е - диэлектрическая проницаемость. Определим область D (0 < t < T, 0 <y < Ly, 0 < z < Lz), на которую наложим сеточную область Dh. В этом случае проекцию Ex определим в узлах {(tn,ym,zk): t„=nht, n = 0, N, N=T/ht, ym=mhy, m = 0,M, M=L/hy, zk=khz, k = 0, K, K=Lz/hz}, проекцию Hy -{(tn+0,5,ym, zk+0,5): tn+0,5=(n+0,5)ht, n = 0, N -1, ym=mhy, m = 1,M -1, zk+0,5=(k+0,5)hz, k = 0, K -1}, проекцию Hz - {(tn+0,5,ym+0,5, zk): tn+0,5 = (n+0,5)ht, n = 0, N -1, ym+0,5 = (m+0,5) hy, m = 0,M -1, zk=khz, k = 0,K-1}. Упрощая задачу, выделим на D подобласть Dc (0 < t < T, L <y <R, B <z < U), при работе с которой нет необходимости в заботе о значениях сеточных функций на краях D, задав на ней Dch. Будем считать,

что Ex"m,k определена в узлах{(4,ym, zk): tn=nht, n = 0, N, N=T/ht, ym=mhy, m = L,R, L=Li/hy, R=Lr/hy, zk = khz, k = B,U , U=LU/ hz, B = Lbh}, Hy^^s -

{(tn+0,5,Ут• Zk+01,5): tn+0,5 = (n+0,5)hfe n = 0, Ж -1, ym = mhy, m = L +1, R-1, zk+0,5 = (k+0,5)hz, k = B,U-1}, Hzm+0,5, k

- {(^и+0,5,Ут+0,5, ¿к): 4+0,5= (п + 0,5)^, п = 0, N -1, Ут+о,5 = (т+0,5)^,, т = ^-1, гк = к = В + 1,Ц-1}. В предложенной области и подобласти индексы т, к

обозначают узлы по пространству, п - по времени. Расстояния между узлами задаются пространственными (h и hz) и временным (^ ) шагами сетки. Сеточное значение диэлектрической проницаемости (е]к) характеризует изучаемый оптический элемент.

Тогда явная разностная схема Yee для системы уравнений (6) выглядит как:

I0

Hv"+1/2 - НУп-1/2 Fx" - Fx"

У m,k+1/2 11ym,k+1/2 ^xm,k+1 ^xm,k

Hzn+1/2 - Hzn-1/2 Fxn - Fxn

m+1/2,k m+1/2,k "m+1,k -L"i'm,k

(7)

Fr"+1 - Fr" Hzn+1/2 - Hzn+1/2

hz"+1/2 - Hz"+1/2

Hm,k+1/2 Hm,k-1/2

Рассмотрим данную явную двухслойную схему в нотации (3) и на её основе запишем первые два выражения (6) в итерационном виде, следуя рекомендациям А.А. Самарского [9]:

(8.1)

Hyk+1 = a1 AExk + Hyk,

f-T T -4 -4

где A :

-T T J2 J2

-T T

Hzk+1 = a2 BExk + Hzk

(8.2)

Í-T 0

где B

-T2 0

-T

■*■ XT

0

TN

N У

a1 = (ht/ |0hz), a2 = (ht/ |0hy) и Tt - диагональная матрица, хранящая значения коэффициентов при Exk . Отбросим краевые значения, так как при добавлении PML данные граничные условия будут считаться с помощью других уравнений, отличных от тех, что используются для пересчёта в середине сетки.

Как и в случае с уравнением теплопроводности, мы приходим к работе с диагональными матрицами, диагонали которых состоят из коэффициентов при Exm,k+1 и Ex"mk для (8.1) и Exm+U и Ex"mk для (8.2), следовательно, алгоритм сводится к описанному в пункте 1.3 алгоритму. Третье уравнение из (7) также считается при помощи «длинных» векторов.

3. Вычислительные эксперименты

В качестве аппаратного обеспечения использован ПК с операционной системой Debian 7.0 (wheezy) и установленными драйверами CUDA Toolkit 4 и компилятором gcc. Также была использована библиотека CUBLAS, являющаяся аппаратно-ориентированной реализацией программного интерфейса BLAS (Basic Linear Algebra Subprograms). Результаты, представленные далее, получены на видеокарте GeForce GTX 660Ti (табл. 1) и процессоре Intel Core2 Duo CPU E8500 (табл. 2).

Таблица 1. Основные характеристики GPU NVIDIA GeForce GTX 660Ti

Характеристика Значение

Количество мультипроцессоров, шт. 1228

Размер видеопамяти, Мб 2048

Максимальное число потоков в блоке, шт. 512

Максимальная размерность блока потоков (х, у, г), шт. 512x512x64

Максимальная размерность сетки блоков, шт. 65535x65535x1

Тактовая частота ядра, МГц 915

Тактовая частота памяти, МГц 1502

2

h

h

z

I

0

h

h

z

ÍJÍJ

0

h

h

y

h

Z

т, б

Таблица 2. Основные характеристики CPU Intel Core2 Duo CPU E8500

Характеристика Значение

Тактовая частота ядра, ГГц 3,16

Тактовая частота шины CPU, МГц 1333

Кеш L1, Кб 64x2

Кеш L2, Кб 6144

Тактовая частота ядра, ГГц 3,16

На рис. 2 представлено сравнение результатов, полученных с помощью первого алгоритма c «длинными» векторами, реализованного для двумерного уравнения теплопроводности, и с помощью коротко-векторного для сеточных областей размером N х N (N = 1..1000). Реализация алгоритмов происходила с использованием таких функций стандартной библиотеки CUBLAS, работающей с GPU, как cublasSaxpy(), cublasSgemv(). Время работы алгоритма с «короткими» векторами растёт линейно с увеличением N, т.к. происходит линейное увеличение количества вызовов функций библиотеки CUBLAS, а следовательно, и количества арифметических операций. Алгоритмы, оперирующие «длинными» векторами, показывают незначительную разницу по времени работы программы как для сеточной области с N = 100, так и для N = 1000. Выигрыш алгоритма, работающего с «длинными» векторами, обусловлен, во-первых, меньшим количеством арифметических операций (для первого алгоритма из п. 3), во-вторых, особенностью библиотеки CUBLAS для векторных операций, заключающейся в том, что на каждую векторную операцию запускается отдельное ядро (CUDA Kernel), что привносит фиксированную задержку, зависящую не от длины вектора (до некоторого достаточно большого значения), а только от количества таких векторов.

Т,нс х105_

8 6 4 2

250 500 750 1000

Рис. 2. Сравнение времени работы Т коротковекторного и длинновекторных алгоритмов для уравнения теплопроводности

Второй алгоритм из п. 1.3 работает медленнее первого, т.к. первый, обладая всеми плюсами «длинных» векторов, имеет меньшее количество арифметических операций. На рис. 3 приведено сравнение лучшего алгоритма с «длинными» векторами для уравнения теплопроводности с реализацией на ОрепСшггеШ, где видно, что первый алгоритм из п. 3, обозначенный пунктиром, работает быстрее - его ускорение составило 2.

5 4

3 2 1

250 500 750 1000

Рис. 3. Сравнение времени работы T длинновекторных алгоритма и пакета OpenCurrent

На рис. 4 приведены результаты для алгоритма, работающего с «длинными» векторами для FDTD-метода, и реализация этой же задачи при помощи пакета B-CALM. Алгоритм для «длинных» векторов даёт выигрыш по времени примерно в 1,4 раза. Т, не_

15

10 5

250 500 750 1000

Рис. 4. Сравнение времени работы T длинновекторного алгоритма и реализации двумерного FDTD-метода при помощи пакета B-CALM

Заключение

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

В результате экспериментальных исследований была показана эффективность алгоритмов с «длинными» векторами как для решения уравнения теплопроводности при помощи библиотеки CUBLAS, так и для двумерного варианта FDTD-метода, что доказывает актуальность предложенного подхода к решению исследуемой задачи.

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

Благодарности

Работа выполнена при поддержке Министерства образования и науки РФ и грантов РФФИ 14-07-00291А и 14-07-31178мол_а.

----Алгоритм с «короткими» векторами Алгоритм с «длинными» векторами - первый ---Алгоритм с «длиннцми» векторами - второй У У у у

у у

У У -У

У У .у

у ■а* ------

не х104

----Первый алгоритм с «длинными»

-OpenCurrent

VP""*

N

----Алгоритм с «длинными» векторами

-B-CALM

jf У

1 W V ^ / ^ у / .И

_____ У N

Литература

1. Глобальный прогноз развития кремниевых технологий [Электронный ресурс]. - http://wwwrussianelectronics.ru/le-ader-r/review/doc/64629/ (дата обращения 15.12.2014).

2. General-purpose computing on graphics processing units [Электронный ресурс]. - http://www.gpgpu.org/ (дата обращения 18.12.2014).

3. Flynn, M.J. Very High-Speed Computing System / M.J. Flynn // Proceedings IEEE - 1966. - Vol. 54(12). -P. 1901-1909.

4. Ortega, J. Introduction to Parallel and Vector Solution of Linear Systems / J. Ortega. - NY: Plenum Press, 1987. -313 p.

5. Golub, G.H. Matrix Computations / G.H. Golub, Ch.F. Van Loan. - JHU Press, 1996. - 747 p.

6. Murray, C.M. The Supermen: The Story of Seymour Cray and the Technical Wizards Behind the Supercomputer / Ch.J. Murray. - Wiley, 1997. - 232 p.

7. NVIDIA GeForce [Электронный ресурс]. - http://www.nvi-dia.ru/object/geforce-family-ru.html (дата обращения 18.12.2014).

8. NVIDIA Tesla [Электронный ресурс]. - http://www.nvi-dia.ru/object/tesla-high-performance-computing-ru.html (дата обращения 16.12.2014).

9. Самарский, А.А Введение в теорию разностных схем / А.А. Самарский - М.: Наука, 1971. - 552 с.

10. Волков, К.Н. Численное решение задач гидродинамики на графических процессорах общего назначения / К.Н. Волков, В.Н. Емельянов, А.Г. Карпенко, И.В. Ку-рова, А.Е. Серов // Вычислительные методы и программирование - 2013. - T. 14(1). - С. 82-90.

11. OpenCurrent is an open source C++ library for solving Partial Differential Equations (PDEs) over regular grids using the CUDA platform from NVIDIA [Электронный ресурс]. - https://code.google.com/p/opencurrent/ (дата обращения 18.12.2014).

12. Wahl, P. B-CALM: An open-source GPU-based 3D-FDTD with multi-pole dispersion for plasmonics / P. Wahl, D.S. Ly-Gagnon, Ch. Debaes, D.A.B. Miller, H. Thienpont // Optical and Quantum Electronics. - 2012. - Vol. 44. -P. 285-290.

13. Golovashkin, D.L. Solving finite-difference equations for diffractive optics problems using graphics processing units / D.L. Golovashkin, D. Vorotnikova, A. Kochurov, S. Ma-lysheva // Optical Engineering. - 2013. - Vol. 52(9). - Doi: 10.1117/1.0E.52.9.091719.

14. Головашкин, Д.Л. Решение сеточных уравнений на графических вычислительных устройствах. Метод пирамид / Д.Л. Головашкин, А.В. Кочуров // Вычислительные технологии. - 2012. - Т. 17, № 3. - C. 55-69.

15. Дифракционная нанофотоника / А.В. Гаврилов, Д.Л. Головашкин, Л.Л. Досколович, П.Н. Дьяченко, А.А. Ковалёв, В.В. Котляр, А.Г. Налимов, Д.В. Нестеренко, В.С. Па-вельев, Р.В. Скиданов, В.А. Сойфер, С.Н. Хонина, Я.О. Шуюпова. - Под ред. В.А. Сойфера. - М.: Физмат-лит, 2011. - 680 с.

References

1. The global forecast for the development of silicon technology [Electronical Resource]. - http://www.russianelectro-nics.ru/leader-r/review/doc/64629/ (request date 15.12.2014). -(In Russian).

2. General-purpose computing on graphics processing units [Electronical Resource]. - http://www.gpgpu.org/ (request date 18.12.2014).

3. Flynn, M.J. Very High-Speed Computing System / M.J. Flynn // Proceedings IEEE. - 1966. - Vol. 54(12). -P. 1901-1909.

4. Ortega, J. Introduction to Parallel and Vector Solution of Linear Systems / J. Ortega. - New York: Plenum Press, 1987. - 313 p.

5. Golub, G.H. Matrix Computations / G.H. Golub, Ch.F. Van Loan. - JHU Press, 1996. - 747 p.

6. Murray, C.M. The Supermen: The Story of Seymour Cray and the Technical Wizards Behind the Supercomputer / Ch.J. Murray. - Wiley, 1997. - 232 p.

7. NVIDIA GeForce [Electronical Resource]. - http://www.nvi-dia.ru/object/geforce-family-ru.html (request date 18.12.2014).

8. NVIDIA Tesla [Electronical Resource]. - http://www.nvi-dia.ru/object/tesla-high-performance-computing-ru.html (request date 16.12.2014).

9. Samarski, A.A. Introduction into the theory of difference schems / A.A. Samarski. - Moscow: "Nayka" Publisher, 1971. - 552 p. - (In Russian).

10. Volkov, K.N. Numerical solution of hydrodynamics problems by means of graphical processors / K.N. Volkov, V.N. Emelanov, A.G. Karpenko, I.V. Kyrova, A.E. Serov // Calculating methods and programming. - 2013. -Vol. 14(1). - P. 82-90. - (In Russian).

11. OpenCurrent is an open source C++ library for solving Partial Differential Equations (PDEs) over regular grids using the CUDA platform from NVIDIA [Electronical Resource]. -https://code.google.com/p/opencurrent/ (request date 18.12.2014).

12. Wahl, P. B-CALM: An open-source GPU-based 3D-FDTD with multi-pole dispersion for plasmonics / P. Wahl, D.-S. Ly-Gagnon, Ch. Debaes, D.A.B. Miller, H. Thienpont // Optical and Quantum Electronics. - 2012. - Vol. 44. - P. 285-290.

13. Golovashkin, D.L. Solving finite-difference equations for diffractive optics problems using graphics processing units / D.L. Golovashkin, D. Vorotnikova, A. Kochurov, S. Maly-sheva // Optical Engineering. - 2013. - Vol. 52(9). - Doi: 10.1117/1.0E.52.9.091719.

14. Golovashkin, D.L. Solution of difference equations in graphical computing devices. Method of pyramids / D.L. Golovashkin, A.V. Kochurov // Computing Technologies. - 2012. - Vol. 17(3). - P. 55-69. - (In Russian).

15. Diffractive nanophotonics / A.V. Gavrilov, D.L. Golovashkin, , L.L. Doskolovich, P.N. Dyachenko, A.A. Kovalev, V.V. Kot-lyar, A.G. Nalimov, D.V. Nesterenko, V.S. Pavelyev, R.V. Ski-danov, V.A. Soifer, S.N. Khonina, Ya.O. Shuyupova. - Ed. by V.A. Soifer. - Moscow, "Fizmatlit" Publisher, 2011. - 680 p. -(In Russian).

LONG VECTORS ALGORITHMS FOR SOLVING GRID EQUATIONS OF EXPLICIT DIFFERENCE SCHEMES

D.G. Vorotnikova, D.L. Golovashkin Image Processing Systems Institute of the Russian Academy of Sciences, Samara State Aerospace University

Abstract

We propose two variants of long vectors algorithms for solving grid equations of explicit difference schemes, allowing one to use the maximum number of CUDA cores even for a small dimension of the grid domain. The implementation of these algorithms is shown by the example

of the heat conduction equation and the solution of Maxwell's equations. The comparison between the proposed algorithms implemented by means of the CUBLAS library and free software packages B-CALM and OpenCurrent is done.

Keywords: heat equation, Maxwell's equations, difference scheme, PDE, CUBLAS.

Сведения об авторах

Воротникова Дарья Геннадьевна, 1989 года рождения. В 2012 году окончила Самарский государственный аэрокосмический университет по направлению «Прикладная математика и информатика». Работает стажёром-исследователем в Институте систем обработки изображений Российской академии наук и является аспирантом кафедры наноинженерии Самарского государственного аэрокосмического университета. Область научных интересов: векторные и параллельные вычисления, FDTD- и BPM-методы.

E-mail: davavorotnikova@gmail. com .

Daria Gennadievna Vorotnikova, (b. 1989). Graduated (2012) from S.P. Korolyov Samara State Aerospace University (SSAU). Since 2012 she is the graduate student of SSAU and the trainee-researcher of the Image Processing Systems Institute of the RAS. Scientific interests: BPM- and FDTD-method, vector and parallel algorithms for matrix computation.

Головашкин Димитрий Львович, доктор физико-математических наук, доцент, ведущий научный сотрудник Института систем обработки изображений РАН. Область научных интересов: разностное решение уравнений Максвелла (FDTD-метод), дифракционная оптика, векторные и параллельные матричные вычисления.

E-mail: dimitriv@smr.ru .

Dimitry Lvovich Golovashkin, Doctor of Physical and Mathematical Sciences, Associate Professor, Leader Researcher of the Image Processing Systems Institute of the Russian Academy of Sciences. Scientific interests: FDTD-method, sub wave optics, vector and parallel algorithms of matrix computation.

Поступила в редакцию 20 января 2015