Блочные алгоритмы решения сеточных уравнений Zheng/Chen/Zhang

Развитие вычислительных оптики и электродинамики в XXI веке неразрывно связано с FDTD-методом (Finite-Difference Time-Domain), широко применяемым [1, 2] для моделирования разнообразных оптических и электродинамических явлений в рамках строгой теории дифракции. С точки зрения теории разностных схем принято различать два подхода внутри FDTD-метода, основанных на разработке явных и неявных схем. К достоинствам первого следует отнести простоту понимания и реализации, важную для исследователей, далеких от вычислительной математики (ведь метод предназначен для физиков); второго – абсолютную устойчивость, столь ценимую специалистами-математиками. Длительное время упомянутое разделение по специализации сдерживало развитие неявных разностных схем; так, между публикациями первых работоспособных явной [3] и неявной [4] схем прошло целых 33 года. Со временем необходимость проведения исследований в тех областях оптики и электродинамики, где ограничения, налагаемые критерием Куранта, представляются чрезмерными [1], обусловила развитие математического аппарата неявных схем для моделирования распространения излучения: в средах с частотной дисперсией [5], плазме [6], тонких ферромагнитных пленках [7] и др.

Суждение о высокой вычислительной сложности FDTD-метода в полной мере относится и к его неявным схемам. Практически сразу с момента их появления основным способом преодоления данной проблемы являлась организация вычислений на суперкомпьютерах: от легендарных Cray в начале века [8] до использующих графические процессоры в настоящее время [9]. Вместе с тем получают актуальность другие, более доступные широкому кругу исследова-

телей способы сокращения длительности вычислений по FDTD-методу, в частности, разработка и реализация блочных алгоритмов. Для явных схем FDTD-метода такие алгоритмы появились относительно недавно [10, 11], хотя сам прием составления блочных алгоритмов для снижения вычислительной сложности имеет длительную историю, пусть и в другой предметной области – матричных вычислениях [12].

Задачей авторов настоящей работы была демонстрация возможности синтеза различных блочных алгоритмов для неявной разностной схемы FDTD-метода, чего ранее в вычислительной практике не делалось. В начале нового направления исследований принято решение ограничиться одномерным случаем классической схемы Zheng/Chen/Zhang [13], который при успехе далее можно развивать до многомерного.

1. Векторные алгоритмы решения сеточных уравнений Zheng/Chen/Zhang

Ориентируясь на материал фундаментальной монографии [1], разностную схему Zheng/Chen/Zhang в одномерном случае можно представить следующим образом. Этап I – Переход со слоя n на слой n + 1 /2:

S o

77 n + 1/2 Z7 n

EXk EXk ht/2

H ⁿ

У к + 1/2

У к - 1/2

h z

Ц о

Hr , + 1/2 _ H n

У к + 1/2 ___________ У к + 1/2

h t /2

h z

Этап II – переход со слоя n + 1 /2 на слой n + 1:

S o

77 n +1    77 n +1/2

Exk    Exk ht/2

H n + 1 _ H + 1

У к + 1/2 __________ У к - 1/2

h z

Ц о

H n + 1 _ H n + 1/2

У к + 1/2 __________ У к + 1/2

h t /2

E n + 1_ E n + 1 ^х к + 1         ^х к

Диэлектрическая и магнитная проницаемости для простоты приняты равными 1 и в выражениях (1 –4) не фигурируют, ε 0 , μ 0 – электрическая и магнитная постоянные. Сеточная проекция E _x ⁿ _k напряженности электрического поля на ось X определена в узлах

{( t n , z k ): t n = n ⋅ h t , n =0, 1/2, 1, 3/2,.., N ; h t = T ∕ N ;

zk =(k –1)⋅hz, k =1, 2,..., K; hz = L ∕ (K –1)}, а сеточная проекция Hnt+1/2 напряженности магнитного поля на ось Y определена в узлах

{( t n , z k +1∕2 ): t n = n ⋅ h t , n =0, 1 /2, 1, 3 / 2,..., N ;

z k +1/2 = ( k – 1 ∕ 2)⋅ h z , k = 1, 2,.., K – 1}.

Значения T и L задают протяженность вычислительной области по времени и пространству, h t и h z – соответствующие шаги дискретизации, n и k – индексы дискретизации.

В качестве начального условия примем отсутствие поля при t =0 соответственно, E ⁰ и H ⁰ положим ,      x k             У к + 1/2

равными 0 при любых k. На правой границе установим электрическую стенку EnK = 0 ; на левой - «жесткий» источник [1] En = sin(2nchtn), где c - скорость света в свободном пространстве. При выбранном способе наложения сеточной области необходимости в задании магнитного поля на ее краях не возникает.

Отметим сходства и отличия схемы Zheng/Chen/Zhang из [1] от классической схемы Yee [3]. В обоих случаях для определения каждой сеточной проекции напряженности полей используются свои узлы и их подмножества не пересекаются. Вместе с тем таких узлов в схеме Zheng/Chen/Zhang задается вдвое больше – в отличие от Yee здесь обе сеточные проекции вычисляются как на слое n , так и на слое n + 1 /2. Указанное обстоятельство определяет двухэтапный переход от слоя n к слою n + 1. Сеточные уравнения первого этапа (1), (2) следует признать явными (как в Yee), уравнения (3), (4) второго этапа очевидно следует отнести к неявным (в отличие от Yee). Последнее обуславливает как преимущество схемы Zheng/Chen/Zhang – абсолютную устойчивость, так и ее недостаток – необходимость решения систем линейных алгебраических уравнений в ходе производства вычислений.

Определяя вид системы, появляющейся на этапе II, дважды подставим (4) в (3), выразив из (4) H " ^{+ 1} и H " ^{+ 1} , получив в итоге сеточное уравнение У к + 1/2 У к - 1/2 ' v .гл.

t    E« + 1

4 h _z ²So P o x ^- 1

+ 1 + —h   E n ^{+ 1} h t E n ^{+ 1} = ( 2 h |s o ^ o ) x^k 4 h _z ²S o Ц o x^{k +} 1

= E ⁿ + 1/2

h t

2 h _z S o

( Hⁿ + 1/2

^{( 1} Ук + 1/2

TJ n + 1/2 \

^НУ к - 1/2 )

Варьируя k от 2 до K – 1, получаем систему трехдиагонального вида, после решения которой этап II завершается подстановкой найденных значений се- точной напряженности электрического поля на слое n + 1 в формулу (4).

Признавая вслед за авторами работы [9] уместность метода Томаса (известен как метод прогонки в отечественной литературе) для решения набора таких систем в многомерном случае, отметим следующую особенность рассматриваемого одномерного. На этапе II возникает лишь одна система и векторизация вычислений по ней при использовании метода Томаса невозможна, как впрочем и для любого другого прямого метода [14]. Следовательно, решение (5) необходимо проводить итерационными методами, например Якоби или Гаусса–Зейделя. И хотя в вычислительной математике их считают «простыми», однако ранее для синтеза векторных [14], а теперь и блочных [15] алгоритмов принято использовать именно их; возможно, именно в силу упомянутой «простоты», допускающей их эффективную векторизацию и удачный переход к блочности. Сходимость обоих методов в данном случае обеспечивается в соответствии с теоремами 3.1.1 и 3.2.1 из [14], строгим диагональным преобладанием матрицы системы, образованной левой частью (5).

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

Алгоритм 1.

01 do n=1,N

02    Exx=Ex;

03    Ex(2:Nz-1)=Ex(2:Nz-1)-c4*(Hy(2:Nz-1)-

Hy(1:Nz-2));

04    Hy=Hy-c1*(Exx(2:Nz)-Exx(1:Nz-1));

05    Ex(1)=sin(c5*n);

06 b=c6*(c4*(Hy(1:Nz-2)-Hy(2:Nz-1))+

Ex(2:Nz-1));

07    do g=1,M

08      Ex(2:Nz-1)=c2*(Ex(1:Nz-2)+Ex(3:Nz))+b;

09    end do

10    Hy=Hy-c1*(Ex(2:Nz)-Ex(1:Nz-1));

• 11    end do

В теле цикла 01 – 11 производится переход со слоя n на слой n + 1 в два этапа: I, строчки 02–04, и II, строчки 06– 10. При этом вычисления по формуле (1) производятся в строчке 03; по (2) – в строчке 04; по (5) – в 07–09, где в цикле реализуются M итераций метода Якоби для решения системы; по (4) – в строчке 10. В строчке 02 производится дублирование вектора Ex в Exx с целью предотвращения преждевременного затирания значений E _x ⁿ _K , необходимых при реализации (2). В 05 задается «жесткий источник», в 06 формируется правая часть системы. Константы c1, c4, c3, с2, с6 и c5 определяются до алгоритма как h t /(2 h z µ 0 ), h t /(2 h z ε 0 ), c1*c4, c3/(1+2*c3); 1 /(1 + 2*c3) и 2π ch t соответственно.

На основе Алгоритма 1 предложим Алгоритм 2, также реализующий метод Якоби, однако несколько иначе. Пусть вместо цикла 07–09 будет записано:

07    do g=1,M/2

0 8      Exx(2:Nz-1)=c2*(Ex(1:Nz-2)+Ex(3:Nz))+b;

0 9      Ex(2:Nz-1)=c2*(Exx(1:Nz-2)+Exx(3:Nz))+b;

• 10    end do

Тогда строчки 10– 11 Алгоритма 1 получат номера 11 – 12 в Алгоритме 2. В ходе нечетных итераций по g , строка 08, приближенное решение системы будет записываться в вектор Exx (который ранее уже был введен, пусть и с другой целью), в ходе четных, строка 09, – в Ex. Предложенная модификация, как будет показано в ходе экспериментального исследования, позволит сэкономить память и ускорить вычисления.

При векторизации метода Гаусса–Зейделя принято [14] использовать его вариант, основанный на красночерном упорядочивании. Тогда для обсуждаемой разностной схемы можно составить третий векторный алгоритм (Алгоритм 3), отличающийся от первого небольшим фрагментом. Вместо строчки 08, в нем будет записано две следующих:

08 Ex(3:Nz-2:2)=c2*(Ex(2:Nz-3:2)+Ex(4:Nz-1:2))+

b(2:Nz-3:2);

09 Ex(2:Nz-1:2)=c2*(Ex(1:Nz-2:2)+Ex(3:Nz:2))+

b(1:Nz-2:2);

В первой реализуется расчет итерационного приближения в «красных» (нечетных) позициях массива через использование значений на предыдущей итерации, во второй – в «черных» (четных) позициях через использование только что найденных значений на текущей итерации. Соответственно, строчки 09– 11 Алгоритма 1 получат номера 10– 12 в Алгоритме 3. Следует различать итерационные приближения (в методах Якоби или Гаусса–Зайделя) для решения (5), реализуемого в цикле по s , и переход на следующий временной слой по формулам (1, 2, 5) и (4), в цикле по n .

• 2.    Блочные алгоритмы решения сеточных уравнений Zheng/Chen/Zhang с использованием метода Якоби

Поясняя смысл и особенности блочных алгоритмов, предостережем читателя от путаницы в терминологии. Далее речь пойдет не о «блочных методах» из [14], которые предназначены для работы с блочными матрицами (таких в нашем изложении нет), а о приеме организации вычислительного процесса, позволяющем сократить длительность расчетов за счет снижения коммуникационных издержек. В классических работах по теории разностных схем, вышедших как в период становления этой теории [16], так и значительно позже [17], вычислительная сложность численного метода определяется как количество скалярных арифметических операций по нему и только. Однако в соседней предметной области, вычислитель- ной линейной алгебре, самое серьезное внимание (и это отражено в фундаментальных монографиях [12, 18]) всегда уделялось учету коммуникаций внутри иерархической структуры памяти ЭВМ и синтезу алгоритмов (их и называют блочными), предназначенных для минимизации таких коммуникаций. Основной идеей блочности является максимально полное использование данных, загруженных в верхний уровень иерархической памяти. Один из смыслов настоящей работы ее авторы видят в распространение практики составления блочных алгоритмов на теорию разностных схем.

Обращаясь к решению сеточных уравнений (1, 2, 5, 4), зададимся целью перейти от Алгоритма 2 к новому путем фрагментации сеточной области (рис. 1), допускающей раздельное решение (5) по методу Якоби для каждого фрагмента (блока). В этом случае значения E _x ⁿ _k ^{+ 1} формируются не одновременно для всех k (как в строчках 07– 10 Алгоритма 2), а лишь для определенного интервала. Такого, что все необходимые данные помещаются в верхний уровень иерархической памяти.

слой n+1

• - ф--G—6--G—G—ф- итерация 3

• - ф--G—G--G—G— Щ-итерация 2

• - Ф—G—G—G—G—Ф- итерация 1

слой п+1/2 ФПОПОПОПОФ итерация 0 слойп ФПОПОПОПОЯ

Рис. 1. Блок сеточной области при использовании метода Якоби

Так, на рис. 1 изображены узлы сеточной области и значения итерационных приближений, составляющие определенный блок. Левее – их узлы и значения, относящиеся к предыдущему блоку, правее – к следующему. Представлен частный случай для ширины блока dl =4 и M =4. Черные окружности обозначают узлы, на которых определены сеточные проекции электрического поля, границы квадратов – магнитного. Серым цветом выделены значения (кроме нулевого и последнего) итерационных приближений. Круги и квадраты соответствуют узлам и значениям из соседних блоков. Отметим, что нулевое итерационное приближение составляется из значений E _x ⁿ _k ^{+ 1/2} , а последнее задает E _x ⁿ _k ^{+ 1} .

Смещение границ блока на одну позицию влево при переходе на следующую итерацию обусловлено необходимостью соблюдения информационных зависимостей на пространстве итерации циклической конструкции 07– 10 Алгоритма 2. Хранение значений на четных и нечетных итерациях в различных массивах (Ex и Exx) также необходимо для удовлетворения упомянутых информационных зависимостей. В работе [15] предлагается соблюдать зависимости за счет хранения каждого итерационного приближения в отдельном массиве, что приводит к значительному увеличению объема занимаемой памяти. Следующий блочный алгоритм лишен этого недостатка.

Алгоритм 4.

01 do n=1,N

02     { ЛЕВАЯ ТРАПЕЦИЯ }

03    do q=2,p-1

04      wl=(q-1)*dl+2; wr=q*dl+1;

05      Exx(wl:wr)=Ex(wl:wr);

0 6      Ex(wl:wr)=Ex(wl:wr)-c4*(Hy(wl:wr)-

Hy(wl-1:wr-1));

07      Hy(wl-1:wr-1)=Hy(wl-1:wr-1)- c1*(Exx(wl:wr)-Exx(wl-1:wr-1));

08      b(wl-2:wr-2)=c6*(c4*(Hy(wl-2:wr-2)-

Hy(wl-1:wr-1))+Ex(wl-1:wr-1));

09      do g=1,M/2

• 10          wl=wl-1; wr=wr-1; Exx(wl:wr)=

c2*(Ex(wl-1:wr-1)+Ex(wl+1:wr+1))+b(wl-1:wr-1);

• 11          wl=wl-1; wr=wr-1; Ex(wl:wr)=

c2*(Exx(wl-1:wr-1)+Exx(wl+1:wr+1))+b(wl-1:wr-1);

• 12      end do

13      Hy(wl-1:wr-1)=Hy(wl-1:wr-1)- c1*(Ex(wl:wr)-Ex(wl-1:wr-1));

• 14    end do

• 15    { ПРАВАЯ ТРАПЕЦИЯ }

• 16    end do

• 3.    Блочные алгоритмы решения сеточных уравнений Zheng/Chen/Zhang с использованием метода Гаусса–Зейделя

Вычисления в центральных блоках, имеющих форму, близкую к параллелепипеду, расписаны подробно (строки 03– 14); первый и последний блок по форме напоминают трапецию, о вычислениях в них для краткости лишь упоминается (строки 02 и 15). Сравнивая Алгоритмы 2 и 4, отметим основные отли- чия. Векторные операции по формулам (1, 2, 5) и (4) теперь проводятся не для всей сеточной области, а для заданного фрагмента с границами wl (левая граница) и wr (правая). Границы фрагмента (блока) смещаются влево с увеличением номера итерации (начало строк 10 и 11). Сами блоки перебираются слева направо в цикле по q (строки 03 – 14), где q– номер блока. Общее количество блоков определяется до выполнения обсуждаемых алгоритмов по правилу p = K / dl.

Можно избежать хранения массивов Exx и b на всей сеточной области, ограничившись лишь той их частью, что необходима при работе с текущим блоком. Ведь для начала вычислений в каждом блоке необходимы лишь массивы Ex и Hy, а Exx и b формируются далее с их использованием. Оформим это в Алгоритме 5, характеризующемся в силу указанного обстоятельства вдвое меньшими (ведь dl<< K ) затратами памяти по сравнению с Алгоритмами 2 и 4.

Метод Гаусса–Зейделя не только выгодно отличается от метода Якоби лучшей сходимостью (в большинстве случаев), но и служит основой для построения других, более развитых итерационных методов [14]. Декомпозиция сеточной области на блоки при его применении несколько изменится (рис. 2).

слой n+1 ♦□ODODODO H-- итерация 4

•—e—«—e— q—• итерация 3

e—e—о—e— q—• итерация 2

•—e—e—e—e—• итерация 1 слой n+1/2 —® ПОПОПОПО1 итерация О слой п —• □ODODODO1

Рис. 2. Блок сеточной области при использовании метода Гаусса–Зейделя

На рис. 2 черные окружности обозначают узлы, на которых определены сеточные проекции электрического поля, границы квадратов – магнитного. Серым цветом выделены значения (кроме нулевого) итерационных приближений. Круги и квадраты соответствуют узлам и значениям из соседних блоков. Закрашенные сверху окружности – «черные» позиции массивов, снизу – «красные».

Действительно, при переходе к новому итерационному приближению (за исключением первого) границы блока теперь сдвигаются на две позиции влево (строчка 12 Алгоритма 6), что вызвано необходимостью согласования позиций: «красные» должны быть под «красными», «черные» – под «черными». Начинается блок всегда с «черной» позиции, а заканчивается «красной»; это обусловлено состоянием вычислительного процесса, при котором значения на текущей итерации слева от блока уже вычислены, а справа еще нет. При работе с «красными» позициями используются значения, найденные только на предыдущей итерации, с «черными» – на текущей.

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

Алгоритм 6.

01 do n=1,N

02     { ЛЕВАЯ ТРАПЕЦИЯ }

03    do q=2,p-1

04      wl=(q-1)*dl+2; wr=q*dl+1;

05      Exx(wl+1:wr+1)=Ex(wl+1:wr+1);

06      Ex(wl+1:wr+1)=Ex(wl+1:wr+1)- c4*(Hy(wl+1:wr+1)-Hy(wl:wr));

07      Hy(wl:wr)=Hy(wl:wr)-c1*(Exx(wl+1:wr+1)-

Exx(wl:wr));

08      b(wl-1:wr-1)=c6*(c4*(Hy(wl-1:wr-1)-

Hy(wl:wr))+Ex(wl:wr));

09      do g=1,M

10         Ex(wl+1:wr:2)=c2*(Ex(wl:wr-1:2)+

Ex(wl+2:wr+1:2))+b(wl:wr-1:2);

11         Ex(wl:wr-1:2)=c2*(Ex(wl-1:wr-2:2)+

Ex(wl+1:wr:2))+b(wl-1:wr-2:2);

• 12         wl=wl-2; wr=wr-2;

• 13      end do

• 14      Hy(wl+1:wr+1)=Hy(wl+1:wr+1)-

• c1*(Ex(wl+2:wr+2)-Ex(wl+1:wr+1));

• 15    end do

• 16    { ПРАВАЯ ТРАПЕЦИЯ }

• 17    end do

4.    Экспериментальное исследование построенных алгоритмов

Как и в предыдущем параграфе, здесь тоже уместен разговор о построении нового алгоритма (Алгоритма 7), отличающегося от только что представленного двукратным сокращением затрат памяти за счет хранения только тех частей массивов Exx и b, что необходимы в ходе вычислений внутри текущего блока.

Как известно, единственным научным критерием истины является практика. Эксперименты проводились на вычислительной системе, оснащенной процессором Intel Core i3-4150 3,5 ГГц (кэш L1 =64 Кб, L2 =512 Кб, L3 =3 Мб), материнской платой Gigabyte GA-Z97M-DS3H (частота системной шины DMI 5000 МГц), ОЗУ с 12 Гб памяти (DIMM DDR3-1333 МГц) и работающей под управлением операционной системы Ubuntu 18.04.4. Использовался компилятор GCC 7.5.0, ключ компиляции -O3.

Предваряя представление экспериментальных данных, отметим, что рассмотрению подлежала исключительно длительность вычислений по алгоритмам при равных параметрах сеточной области. Проблемы устойчивости, аппроксимации, корректности и сходимости не рассматривались, как выходящие за границы предлагаемой работы (для схемы Zheng/Chen/Zhang все это подробно изложено в [1]). Новизной здесь характеризуются лишь авторские блочные алгоритмы, перечисленные характеристики которых полностью совпадают с соответствующими характеристиками не блочных. В ходе программной реализации каждого блочного алгоритма проводилась проверка совпадения разностных решений в блочном и неблочном случаях. Правильно реализованным признавался лишь тот алгоритм, где такое совпадение было полным.

Целью первой серии экспериментов был подбор сеточной области, обеспечивающей, с одной стороны, достаточно длительные вычисления по неблочным алгоритмам – случайные системные события не должны значимо влиять на результаты; с другой стороны, нет смысла в избыточно большом времени расчетов – ускорение блочных алгоритмов не зависит от числа временных слоев сеточной области. Выбор был остановлен на следующих значениях: N = 50, K = 1e8.

Количество итераций M = 16 в методах Якоби и Гаусса– Зейделя обеспечивало приемлемую точность разностного решения на заданной сеточной области. Результаты первой серии экспериментов представлены в табл. 1.

Табл. 1. Результаты первой серии экспериментов

Алгоритм	Длительность вычислений, с	Занимаемая память, Гб
1	268,06	от 1,5 до 1,9
2	168,72	1,5
3	232,39	1,5

Как видно, Алгоритм 2 характеризуется меньшим временем расчетов по сравнению с Алгоритмом 1 за счет лучшей работы с памятью. При вычислениях по Алгоритму 1 объем выделяемой памяти непрерывно варьировался в указанных пределах (табл. 1), превышая ожидаемое значение в 1,5 Гб (четыре массива длиной K действительных чисел одинарной точности). Очевидно, это было предусмотрено при компилировании для предотвращения возможности преждевременного затирания значений массива Ex в строке 08 Алгоритма 1, для чего массив автоматически дублировался. Во втором Алгоритме необходимости в таком дублировании нет. Поэтому блочные Алгоритмы 4 и 5 строились на его основе.

Отметим, что расчеты по Алгоритму 3 производились заметно дольше, чем по Алгоритму 2 при тех же затратах памяти и одинаковом количестве скалярных арифметических операций. Причинами тому могут служить особенности производства векторных операций в цикле по g третьего Алгоритма: таких операций две на одну итерацию метода Гаусса–Зейделя (вместо одной на одну итерацию Якоби) и при составлении участвующих в них векторов шаг выборки 2 (вместо шага 1 у Якоби).

В ходе второй серии экспериментов исследуем зависимость длительности расчетов по двум блочным алгоритмам, основанным на методе Якоби, от ширины блока dl (табл. 2). Остальные параметры остаются прежними.

Табл. 2. Результаты второй серии экспериментов

Алгоритм	Ширина блока	Объем блока	Длительность вычислений, с	Занимаемая память, Мб	Ускорение
4	4 e 1	896 байт	53,79	1526	3,14
	4 e 2	6,5 Кб	38,07		4,43
	4 e 3	62,75 Кб	40,73		4,14
	4 e 4	625,25 Кб	46,65		3,62
	4 e 5	6,12 Мб	135,80		1,25
	4 e 6	61,04 Мб	167,14		1,01
5	4 e 1	896 байт	48,41	765	3,49
	4 e 2	6,5 Кб	29,32	765	5,75
	4 e 3	62,75 Кб	31,81	765	5,30
	4 e 4	625,25 Кб	38,59	766	4,37
	4 e 5	6,12 Мб	130,46	775	1,29
	4 e 6	61,04 Мб	164,89	859	1,02

Под ускорением здесь понимается отношение времени вычислений по Алгоритму 2 к длительности расчетов по блочным Алгоритмам 4 и 5. Если объем занимаемой памяти в табл. 1 и 2 измерялся по показаниям системного монитора, то объем блока оценивался теоретически с учетом его ширины и количества задействованных массивов. Ширина блока выбиралась из следующих соображений: вычисления в правой и левой трапеции (строки 02 и 16 Алгоритма 4) удобно, хотя и не обязательно, проводить при соблюдении условия dl> M, когда используется метод Якоби, и dl>2M, когда метод Гаусса–Зейделя. Для корректности будущего сравнения ориентируемся на последнее неравенство, выбирая минимальную ширину блока равной 40.

Обращаясь к четвертому столбцу табл. 2, наблюдаем U-образную зависимость длительности вычислений от ширины блока для обоих рассматриваемых алгоритмов. Действительно, при росте значений dl время расчетов сначала сокращается в силу уменьшения числа блоков и, как следствие, коммуникационных издержек между оперативной памятью и кэшпамятью процессора. С дальнейшим увеличением dl объем блока выходит за размеры кэш-памяти и эффективность алгоритмов резко падает. Например, при dl=4e4 (блок еще умещается в кэш) ускорения Алгоритмов 4 и 5 составляют 3,62 и 4,37 соответственно; при dl=4 e 5 (объем блока вдвое превосходит емкость уровня L 3) ускорения резко падают до 1,25 и 1,29; для dl=4 e 5 (емкость L 3 превышена в 20 раз) эффект от блочности исчезает совсем. Достаточно высокие ускорения при минимальном значении dl (3,14 и 3,49) объясняются существенной величиной этого значения. Выбрав dl небольшим (что возможно, но приведет к некоторому усложнению алгоритмов), будем наблюдать и небольшие ускорения, как это демонстрируется в [19] на примере явных разностных схем. Проведя дополнительные эксперименты, авторы установили, что при dl=2 e 5 (объем блока совпадает с емкостью L 3) величины ускорений достаточно велики и составляют 3,47 и 3,94 раза. Учитывая все изложенное, причиной ускорения вычислений однозначно следует считать блочность исследуемых алгоритмов.

Сравнивая блочные алгоритмы между собой, отметим преимущество пятого, наилучшее ускорение которого 5,75; в то время как максимальное ускорение четвертого всего 4,43 раза. Различие, очевидно, вызвано сокращением занимаемой памяти вдвое при переходе от четвертого Алгоритма к пятому. Отметим также увеличение занимаемой памяти у пятого Алгоритма при больших блоках. Испытав алгоритмы для других значениях dl, авторы нашли наибольшее ускорение Алгоритма 5 равным 6,12 при dl=2000. Возвращаясь к общему взгляду на теорию разностных схем, отметим несостоятельность оценки длительности вычислений исключительно по количеству арифметических операций. Количество таких операций в Алгоритмах 2 и 5 одинаково, однако время расчетов отличается в 6 раз. Наибольший вклад в длитель- ность, как видно из анализа экспериментальных результатов, оказывают не арифметические операции, а работа с памятью, и при построении новых разностных схем важно учитывать это обстоятельство.

Третья серия экспериментов (табл. 3) посвящена исследованию эффективности Алгоритмов 6 и 7, использующих метод Гаусса–Зейделя для решения сеточных уравнений Zheng/Chen/Zhang.

Табл. 3. Результаты третей серии экспериментов

Алгоритм	Ширина блока	Объем блока	Длительность вычислений, с	Занимаемая память, Мб	Ускорение
6	4 e 1	896 байт	76,65	1526	3,03
	4 e 2	6,5 Кб	76,48		3,04
	4 e 3	62,75 Кб	79,31		2,93
	4 e 4	625,25 Кб	83,15		2,8
	4 e 5	6,12 Мб	100,03		2,33
	4 e 6	61,04 Мб	231,70		1
7	4 e 1	896 байт	76,45	765	3,04
	4 e 2	6,5 Кб	72,34	765	3,21
	4 e 3	62,75 Кб	72,90	765	3,19
	4 e 4	625,25 Кб	74,57	766	3,12
	4 e 5	6,12 Мб	98,98	773	2,35
	4 e 6	61,04 Мб	232,34	841	1

Здесь под ускорением понимается отношение времени вычислений по Алгоритму 3 к длительности расчетов по блочным Алгоритмам 6 и 7.

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

• 1.    Вид зависимости времени расчетов от ширины блока остался прежним, однако левая ветка U-образного графика выражена слабее. Это связано либо с уменьшением максимально достижимого ускорения (ведь сами значения ускорений при dl=4 e 1 в табл. 2 и 3 достаточно близки); либо с особенностями другой зависимости – длительности вычислений от M , которая в данной работе не исследуется. Безусловно, M является блочным параметром, задающим высоту блока и влияющим на ускорение; однако любые его изменения скажутся и на точности разностного решения. Поэтому авторы настоящей работы не рекомендуют играть величиной M , как блочным параметром. Такова отличительная особенность предложенных блочных алгоритмов для неявных разностных схем от алгоритмов для явных (например, из [10, 11]), где можно произвольно варьировать любые параметры блока.

• 2.    При незначительном выходе за емкость кэшпамяти процессора (dl=4 e 5) происходит снижение ускорения, однако не такое резкое, как в табл. 2. Возможно, это объясняется меньшим (2 из 4) по сравнению с Алгоритмами 4 и 5 (там 3 из 4) количеством массивов, участвующих в итерационном процессе – к этому процессу относится большинство арифметических операций внутри блока. Ес-

• ли это так, то теоретическую оценку объема блока в табл. 2 следует умножить на 3/4 а в табл. 3 – на 1 /2. К сожалению, прямое обращение к кэшпамяти процессора при выполнении программы на языке высокого уровня (впрочем, и на низком ситуация немногим лучше) невозможно, и подтвердить высказанное предположение в ходе прямого эксперимента затруднительно.

• 3.    Значения ускорений для Алгоритмов 6 и 7 при одинаковой ширине блоков различаются не столь существенно (максимальное отличие в 1,15 раз для dl=4 e 4), как для Алгоритмов 4 и 5 (максимальное отличие в 1,3 раза для dl=4 e 2), хотя утверждение о двукратной экономии памяти при переходе к алгоритму с большим номером верно для обеих пар алгоритмов. Впрочем, двукратное сокращение объема занимаемой памяти имеет самостоятельную ценность.

• 4.    Наконец, максимально достигнутое ускорение в 3,21 раза (Алгоритм 7 при dl=4 e 2), хотя и существенно (например, в [10, 11] трехкратное ускорение для блочного алгоритма, связанного со схемой Yee, является наибольшим при задействовании одного вычислительного процесса), однако значительно ниже аналогичной характеристики Алгоритма 5.

Причинами последних двух обстоятельств по-видимости являются уже отмеченные при анализе первой серии экспериментов различия Алгоритмов 2 и 3.

Обсуждая развитие предложенного подхода к составлению блочных алгоритмов решения сеточных уравнений Zheng/Chen/Zhang на многомерный случай, следует упомянуть возможность синтеза волновых алгоритмов [20], хорошо зарекомендовавших себя для двумерных явных разностных схем [19]. Кроме того, безусловный интерес представляет реализация блочных алгоритмов решения упомянутых сеточных уравнений на графических процессорах, как это сделано в работе [21] для явной схемы Yee. Задачами такой реализации будет не только ускорение вычислений, но и преодоление ограничения на размер видеопамяти графического ускорителя.

Заключение

Продемонстрирована возможность построения блочных алгоритмов FDTD-метода для организации вычислений по неявной разностной схеме Zheng/Chen/Zhang, характеризующейся абсолютной устойчивостью. Экспериментально подтверждена эффективность предложенных блочных алгоритмов, обеспечивающих ускорение расчетов до шести раз на той же аппаратной и программной базе. Изложены соображения по развитию исследований в выбранной предметной области.

Работа выполнена при поддержке гранта РФФИ 19-07-00423 А.