Параллельно-последовательный алгоритм решения треугольных систем на процессорном кольце

Широкое применение вычислительного эксперимента в научной практике и развитие микропроцессорной архитектуры обуславливают постоянный интерес как к разработке новых численных методов, так и к совершенствованию хорошо известных. Не будет преувеличением отметить задачу решения систем линейных алгебраических уравнений (СЛАУ) как наиболее массовую [1]. Ей посвящены главы широко известных учебников, рассчитанных на читателей без математического образования [2] и на читателей -специалистов [3]. В компьютерной оптике решение СЛАУ является заключительной частью таких разных подходов к моделированию прохождения электромагнитного излучения через дифракционные оптические элементы, как модовые методы, метод связанных волн, разностного решения уравнений Максвелла, методов конечных и граничных элементов решения уравнения Гельмгольца [4].

Вместе с тем, треугольным системам, которым посвящена настоящая работа, в массовой учебной литературе не уделяется специального внимания. Объяснением тому служит очевидность методики нахождения решения таких систем в скалярном случае. При векторизации вычислений рассматриваемая задача несколько усложняется: разработчик оказывается перед выбором строчно-ориентированного алгоритма (основанного на скалярном произведении) либо столбцово-ориентированного (опирающегося на saxpy) [1]. Переход к параллельным вычислениям приводит к алгоритму [5], уже превышающему по сложности LU разложение [6]. В отличие от последнего, характеризуемого для СЛАУ размерности n объемом арифметических операций О(n3) и количеством коммуникаций О(n), ставший классическим алгоритм из [5] связан с О(n2) арифметическими операциями при том же количестве коммуникаций. Доля длительности арифметических операций в общем времени вычислений, а следовательно, и эффективность [7] у алгоритма решения треугольной системы оказывается много ниже по сравнению с параллельным алгоритмом LU разложения. При том, что и тот, и другой используются для решения одной СЛАУ. Указанное несоответствие обуславливает необходимость повышения эффективности параллельного алгоритма из [5], чему и посвящена предлагаемая работа.

Известное утверждение о неактуальности синтеза параллельных алгоритмов решения треугольных СЛАУ (цитируемое в [1]) в силу меньшего количества арифметических операций по сравнению с LU-разложением представляется авторам необоснованным. Например, в работе [8], посвященной разреженным СЛАУ большой размерности, особо указывается на параллельное решение треугольных систем, как на «камень преткновения» при реализации метода неполной факторизации в модификации Айзенштата.

Приведенный параллельно-последовательный алгоритм основан на комбинации векторного варианта решения треугольной СЛАУ из [1] и параллельного варианта из [5]. Поэтому первым будет рассмотрен векторный алгоритм, за ним следует параллельный, завершает изложение авторский вариант метода.

1.    Векторный столбцово- ориентированный алгоритм решения треугольных систем

Выбор в настоящей работе столбцовоориентированного алгоритма решения СЛАУ треугольного вида согласуется с общей ориентацией на столбцово-ориентированные алгоритмы, принятой в [1]. Такое предпочтение объясняется не только традиционным использованием в матричных вычислениях языка программирования Фортран, но и современной практикой. Например, в библиотеке CUBLAS версии 2.2 (появившейся в мае 2009 года) для варианта языка Си (обычно подразумевающего хранение матриц по строкам), используемого при кодировании векторных вычислений на видеокартах фирмы NVIDIA, предполагается хранени е двумерных массивов по столбцам [8].

Для определенности здесь и далее будем рассматривать систему

Lx = b, (1) где матрица L е Rnxn нижняя унитреугольная; x, b е Rnx1. Тогда традиционный столбцовоориентированный алгоритм из [1], записанный в векторной нотации, основанный на операции saxpy и использующий замещение вектора правой части b вектором неизвестных x, выглядит след ующим образом.

Алгоритм 1

for j=1:n – 1 {проход по столбцам матрицы L} {модификация правой части}

b(j+1:n)=b(j+1:n) – b(j)·L(j+1:n,j);

end

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

Алгоритм 2

j=0; z(1:n)=0;

while j

{модификация вектора z} z(j+1:n)=z(j+1:n)+L(j+1:n, j)·b(j);

end

В отличие от Алгоритма 1, в данном введен вектор z G R ^{n x 1}, модифицируемый посредством операции saxpy и накапливающий значения найденных неизвестных b(j), подставленных в уравнения системы (1).

2.    Параллельный алгоритм решения треугольных систем на кольце

Авторы [1,5] производят циклическую декомпозицию L и b, учитывая треугольную структуру матрицы системы и избавляя алгоритм от простоев. Перед началом вычислений в локальной памяти процессора с номерном µ (1≤µ≤p, где p – количество процессоров) хранится матрица L loc =L(1:n, µ: p: n) и столбец b loc =b(µ: p: n). Для простоты в [1] принимается n / p g N . Левый и правый соседи процессора по кольцу обозначаются как left и right; вектор, пересылаемый по кольцу, как w.

Вывод алгоритма начинает ся в [1,5] с представления j-ой компоненты вектора x, относящейся к ве- дению процессора µ, через скалярное произведение фрагмента j-ой строки матрицы L и уже найденны х компонент столбца x.

j ^{- 1}

x (j) =b (j)-E L (j,i)x (i)=

i = 1

min { p,j - 1 }

= b ( ^j )^- E ^L ( ^j,i:p:^{j -} 1 ) ^x ( ^i:p:^{j -} 1 )

i = 1

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

^L(j,i:p:j ^- Dx^pu ^- 1) =

/^(к) ( ^j ) ,i = ^1: ц ^{- 1} _ ^{z(k - 1)} ( ^j ) ,i = Ц :p

,(3)

где

z[k) = L(:, X :p: X + kp) x (X :p: X + kp)

вектор z, сформированный процессором λ на k–ом обороте по кольцу вектора w. Подставляя (3) в (2), получим:

x (j)=ь (j)-Ez(k)(j)-z?-1)( j)- :Ez(k-i,( j) .(5) i=i;

Пусть р-ая компонента вектора w (w G Rpx1), обращающегося по кольцу, формируется как w (u) = Ez(k)(j)+ Ebz‘k-1>( j).(6)

i=1!

Тогда расчет j-ой компоненты вектора неизвестных производится с учетом (5) и (6) по формуле

x(j)= b(j)-w(ц)-zJ,k-1)(j).

След уя [1], запишем параллельный алгоритм решения (1) на кольце.

Алгоритм 3

r=n/p; {число оборотов вектора w по кольцу} col= µ:p:n; {вектор, ставящий в соответствие локальные и глобальные индексы} k=0; {счетчик оборотов по кольцу}

z(1:n)=0; w(1:p)=0; {задание векторов z и w} while k < r {в ектор w сделает r оборотов по кольцу}

If k ≠ 0 or µ >1

recv(w, left); {первый оборот для первого процессора не сопровождается получени ем вектора w от левого соседа} end k=k+1; j=col(k);

b loc (k)=b loc (k) – z(j) – w(µ); {реализация ф. (7)}

If j < n {если вычисления не завершились}

{модификация части z, необходимой для формирования w} q=min(n,j + p – 1);

z(j+1: q) = z(j+1: q)+ L loc (j+1: q, k)·b loc (k); { реализация ф. (4)} {модификация циркулирующего вектора w в соответствии с ф. (6)} if k < r {если оборот не последний и в ерхняя часть w еще понадобится} w(1: µ – 1)= w(1: µ – 1)+ z(j + p – µ + 1 : j + p – 1);

w(µ) = 0;

end

w(µ + 1: p)= w(µ + 1: p)+z(j + 1 : j + p – µ);

send(w, right); {пересылка w правому сосед у по кольцу}

{модификация остальной части z по ф. (4)} if q < n

z(q +1: n)= z(q +1 : n)+ end end end

Отметим, что каждый процессор хранит и модифицирует свой вектор z, в то время как вектор w формируется сообща.

Авторы настоящей работы реализовали приведенный алгоритм на языке Фортран 77 с применением библиотеки MPI, разворачивая векторные участки алгоритма в циклические конструкции. Расчеты,

L loc (q + 1: n, k)·b loc (k)

результаты которых представлены в табл. 1, производились на учебном кластере СГАУ, объединяющем ЭВМ с процессорами AMD Sempron (тактовая частота 1000 MHz), ОЗУ 1 Gb, работающими под управлением операционной системы Linux. Коммуникации межд у компьютерами осуществлялись посредством сети Ethernet (100 Mb).

Таблица 1. Ускорение (S) двухпроцессорного параллельного алгоритма

n/1000	10	11	12	13	14	15	16	17	18	19	20
S	0,94	0,83	1,04	1,06	1,15	1,21	1,166	1,22	1,28	1,37	1,34

Отметим, что ускорение достигается лишь для значений n≥12000, в то время как для получения ускорения при реализации LU разложения на данном кластере достаточно систем размерностью на два порядка меньше. Даже в предельном рассмотренном случае n=20000 (матрицы большего размера не умещались в ОЗУ) получена весьма скромная для двухпроцессорной вычислительной системы величина S=1,34 (табл. 1). Причиной тому является отмеченная во введении невысокая эффективность Алгоритма 3. С увеличением числа задействованных в расчетах процессоров картина лишь усуг убляется: при n=20000 трехпроцессорный алгоритм показал ускорение S=1,11, четырехпроцессорный – S=1,01. Сокращение длительности арифметических операций не компенсирует рост коммуникационных издержек.

3.    Параллельно-последовательный алгоритм

Анализируя структуру Алгоритма 3, ук ажем, что по мере заполнения вектора неизвестных количество арифметических операций, связанны х с модификацией вектора z, сокращается (линейно), а объем пересылаемых данны х остается постоянным (вектор из р чисел). В силу этого доля длительности арифметических операций в общем времени вычислений на одном обороте вектора w по кольцу падает с уве- личением числа оборотов. Следовательно, на заключительных оборотах по кольцу падает и эффективность распараллеливания.

Предположим, что повысить ускорение можно, разделив алгоритм на две части: параллельную (Алгоритм 3 вплоть до некоторого оборота r 0 включительно) и последовательную (Алгоритм 2, начиная с j=r 0 p+1). При этом для недостаточных значений r 0 ув еличение ускорения не будет достигаться (за счет существенной доли последовательной части в общем времени вычислений), а для избыточных значений даст знать о себе отмеченный выше недостаток Алгоритма 3.

Переход от параллельного участка алгоритма к последовательному сопровождается пересылкой векторов z от всех процессоров выделенному с их суммированием. Получившийся вектор б удет соответствовать z из Алгоритма 2. Вектор правой части b(r 0 p+1:n) собирается из соответствующих компонент, хранящихся в локальной памяти всех процессоров вычислительной системы. То же касается матрицы системы.

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

Таблица 2. Ускорение двухпроцессорного параллельно-последовательного алгоритма для n=10000

r 0 /100	5	10	15	20	25	30	35	4	45
S	0,92	1,02	1,04	1,05	1,06	1,17	0,99	0,95	0,95

В случае n=10000 параллельный алгоритм характеризовался отсутствием приемлемого ускорения (S=0,94, табл. 1), применение параллельнопоследовательного с параметром r 0 =3000 позволило ув еличить ускорение на 25%, (до S=1,17, табл. 2).

Таблица 3. Ускорение двухпроцессорного параллельнопоследовательного алгоритма для n=12000

r 0 /100	10	20	30	40	50
S	1,10	1,22	1,27	1,23	1,12

Для n=12000 использование параллельного алгоритма только начало становиться обоснованным (S=1,04, табл. 1), для параллельно-последовательно- го с параметром r0=3000 ускорение увеличилось на 22%, (до S=1,27, табл. 3).

Таблица 4. Ускорение двухпроцессорного параллельнопоследовательного алгоритма для n=17000

r 0 /100	40	50	60	70	80
S	1,37	1,43	1,43	1,38	1,26

При n=17000 для параллельного алгоритма достигается уже ощутимый результат (S=1,22, табл. 1), использование параллельно-последовательного с параметром r0=5000 характеризуется улучшением этого результата на 17%, (до S=1,43, табл. 3). Приведенное значение ускорения параллельно- последовательного алгоритма представляется достаточным для организации расчетов, как дающее ощутимый выигрыш от использования двухпроцессорной вычислительной системы.

Отметим также соответствие результатов вычислительных экспериментов предположению о наличии оптимального значения r 0 , отклонение от которого в большую либо меньшую сторону приводит к снижению ускорения последовательно-параллельного алгоритма.

Выводы

Предложенный в работе параллельно-последовательный алгоритм решения треугольных систем на кольце в силу структуры решаемой задачи характеризуется меньшей длительностью вычислений по сравнению с параллельным. Для приведенных параметров вычислительных экспериментов наблюдался рост ускорения вычислений от 17% до 24%. Предлагаемый подход к синтезу алгоритмов может быть распространен на другие задачи, связанные с матрицами треугольной структуры (например, для решения разреженных СЛАУ большой размерности [8]).

Работа выполнена при поддержке российско-американской программы «Фундаментальные исследования и высшее образование» ("BRHE"), гранта Президента РФ поддержки ведущих научных школ (НШ-3086.2008.9) и гранта РФФИ 07-07-00210а.м.