CC BY
28
Известия Института математики и информатики УдГУ
2015. Вып. 2 (46)
УДК 519.612
© Н. С. Недожогин, С. П. Копысов, А. К. Новиков
ВАРИАНТ ПАРАЛЛЕЛЬНОГО РАЗЛОЖЕНИЯ В ПРЕДОБУСЛАВЛИВАТЕЛЕ AISM1
Предлагается вариапт параллельного разложения при формировании предобуславливателя. основанного па приближенном обращении Шермапа Моррисопа. Проведены численные эксперименты по решению тестовых систем уравнений па графических ускорителях.
Ключевые слова: решение систем липейпых алгебраических уравнений, предобуславливапие. параллельные алгоритмы. графические ускорители вычислений.
Введение
Разработка и реализация новых подходов к решению больших систем .линейных алгебраический уравнений пред обусловленными итерационными методами являются достаточно трудоемкой задачей. Матрица иредобуславливателя не только должна быть близка к обратной матрице коэффициентов системы уравнений, но и должна допускать эффективно распараллеливаемый алгоритм формирования. Наиболее используемыми на сегодняшний день методами, ориентированными на разреженные матрицы, можно считать метод неполного LU-разложения и метод неполного разложения Холецкого [1]. Имея высокую эффективность, данные методы сталкиваются с известными проблемами при их параллельной реализации.
Особенности архитектуры GPU (ставшей широко используемой в вычислениях) накладывают определенные ограничения на подходы к реализации параллельных алгоритмов формирования предобуелавливателей. Исходя из этого, высоким потенциалом обладают предобу-елавливатели на основе аппроксимации обратной матрицы: полиномиальные (TNS [1] и др.), разреженные аппроксимации обратной матрицы (например, AINV), аппроксимации обратной матрицы в факторизованной форме (такие как FSAI, SPAI и др.) [2 4|, а также методы, основанные на приближенном обращении матриц по формулам Шермана Моррисона (AISM Approximate Inverse Sherman Morrison formula) [5].
Эффективное распараллеливание алгоритмов построения предобуелавливателей связано с возможностью использования блочных вариантов формирования иредобуславливателя, что дает сокращение перемещения данных по уровням иерархии памяти GPU, с операциями линейной алгебры, таких как матрично-векторные и матричные произведения, обладающих большим потенциалом распараллеливания.
В работе будет рассмотрен один из подходов к параллельному разложению при формировании обратного иредобуславливателя AISM [6].
§ 1. Последовательный алгоритм
Рассмотрим оригинальный вариант алгоритма построения иредобуславливателя (предложенного в [7]), так как именно он наиболее близок к предлагаемой в статье параллельной форме.
За основу построения иредобуславливателя вида P = Л-1 возьмем матрицу B той же размерности, что и Л, но с известной обратной матрицей.
Теорема 1.1 ([5]). Пусть B ^ невырожденная матрица, векторы и и v т,акие, что
r = 1 + vT B-1и, r = 0.
Тогда матрица Л = B + uvT является обратимой и ее обращение находится как
Л-1 = B-1 - r-1B-1uvTB-1. (1.1)
Работа поддержана РФФИ (гранты 14 01 00055^а, 14 01 31066-мол_а).
Обозначим Ао — невырожденную матрицу, обращение которой легко вычисляется, например диагональную или единичную. Тогда Ак = Ао + ^к=1 иV2, где к = 1,..., п и А = Ага. Если А^, и^ ^ удовлетворяют выражению (1.1), тогда обращение матрицы А может быть вычислено при п-кратном использовании (1.1):
A-1 = A-1 - £ r-U-V^ A-,1!. (1.2)
k=1
Представим (1.2) в матричной форме:
A-1 - A-1 = ФО-1Ф2, (1.3)
Ф = [A0 1U1, A- 1u2,..., A--11ura], Ф = К A0-1,v2 A-1 A-- ] Q-1 =diag[r-V-1,...,r-1].
Покажем, что факторизация (1.3) записывается без явного вычисления A° через векторы uk, Vfc как
k-1 t2 4-1^ k-1 v2
efc = Ufc_£li4О*,., /, r, /, 1.....„. (1.4)
• n • n
¿=0 ¿=0
Тохда выполняются соотношения
A-Vk = A-1Sfc, u2 A-_11 = $ A-1, (1.5)
rfc = 1+ v2 A-1sfc = 1 + A-1 ufc. (1.6)
С учетом (1.5) соотношение (1.3) запишем в виде
A-1 - A-1 = A-1SQ-1T2A-1, (1.7)
где матрицы S = [s1, s2,..., sn], T = [t1, t2,..., tn], столбцы которых вычисляются по uk,Vfc. Выбор A0, uk, Vfc определим следующим образом [7]:
Ao = gin, Ufc = efc, vfc = (ak - a^22, k = 1,...,n,
где /n, ek — единичная матрица и ее k-столбец; вектора ak, ak — k-ая строка матриц A Ao; g — скаляр. Тогда аппроксимация обратной матрицы и соответствующий предобуславливатель примут вид
P = A-1 = g/„ - g-2UQ-1V2.
Рассмотрим еще один вариант разложения обращаемой матрицы A = W - Z, где W -
n
обратимая матрица, Z = UV2 = ^ ufcv)2, а vfc, ufc такие, что dk = 1 - v)2W-^uk = 0, где
k=1
Wfc = Wo - Ek=1 u^v2.
Задавая выбор матриц W = в ■ diag (A) в > 0 U = I, V = Z2 и следуя соотношениям (1.4), (1.5) и (1.6), получим выражения для вычисления столбцов матриц S и T в виде
■ч = уф - }_j —1-fk = г'к - 2_j —Ь
¿=1 ¿=1
и обратной матрицы в виде
A-1 = W-1 - W-1SD-1T2W-1. (1.8)
Алгоритм 1: Формирование нредобуелавливателя AISM
A = W - Z; W = в ■ diag (A); Z = W - A; U = I; V = ZT;
l for к = 1 to n do
Sk = Uk ; tk = vk ;
for г = 1 to fc — 1 do
5 = (i^V-1,^) ;
if >rs then
_ Sfc = uk - j- ■ Si
8 = {vTkW-\Sl) ;
if ||| >tt then
Ll ^
tk = Vk — -Щ ■ U
for j = 1 to n do
if \(sk)j| < Ts then
L (Sfc)j = о
if |(iA-)il < тт then
L (ifc)j = о
dk = 1 - (tTW-1,Uk) ;
17 S = {si, S2, . . . , Sn},T = {ti,t2, ...,tn} ;
18 D = {di ,d2,... ,dn} ;
19 P2 = w-1 - w-1,SD-1:fTw-1;
Используя стратсх'ию фильтрации но значениям элементов (оставляем элементы, значения которых больше некоторой величины т) при вычислении матриц S, T и основываясь на (1.8), выпишем предобуславливатель, аппроксимирующий обратную матрицу в виде
P2 = W-1 - W-1«S'D-1 TTW-1. (1.9)
Последовательный процесс формирования рассматриваемох'о нредобуелавливателя представлен в виде алгоритма 1.
Основные затраты последовательнохх) а;п'оритма составляет вложенный цикл (строка 4 алгоритма 1), в котором определяются матрицы разложения S, Т. Отметим, что матрицы формируются последовательно так, что существует зависимость по данным.
§ 2. Параллельный алгоритм
Первоначальным вариантом распараллеливания данного алгоритма является вычисление скалярных произведений на GPU. Такой подход показал, что достижение эффективного ускорения существенно охраничено доступом к памяти, и не дал существенного увеличения производительности .
Рассмотрим параллельный алгоритм 2 формирования предобуславливателя P2 (1-9) на GPU, в котором скалярные произведения в строках 5 и 8 алгоритма 1 заменены матрично-векторными произведениями.
Для этого введены матрицы Sk = {s1,..., Sk-1} и Tk = {t1,..., tk-1}, состоящие из столбцов, вычисленных на шаге k - 1. В этом случае k - 1 скалярное произведение выполняется в виде матрично-векторных произведений (строка 4 алгоритма 2), результаты которых обозначены Xdi и xd2) а чеРез (xdi )и (xd2 )i — соответствующие компоненты векторов. Также стоит отметить, что в памяти GPU1 и GPU2 хранятся локадьпые копии векторов Xd1 и Xd2) в то время как Xh1 и х^2 — глобальные векторы в оперативной памяти CPU.
Алгоритм 2: Формирование предобуелавливателя AISM на двух GPU Л = W — Z W = в ■ diag (Л) U = Ц V = ZT
1 for к = 1 to п do
// GPU1 // GPU2 tk = vk; = uk;
yi = W-luk] У2 = W -V;
xdi = Tkyi,Tk = (tb . . . ,tk-1) xd2 = Sky2,Sk = (sb . . . ,sk-l);
xdl — xhl II Пересылка GPU1 — CPU xd2 — xh2 II Пересылка GPU2 — CPU xh2 — xd2 // Пересылка CPU — GPU 1 xhl — xdl // Пересылка CPU — GPU2 for i = 1 to к — 1 do for i = 1 to к — 1 do
if |
(xd2 )i
_ tk
= tk —
> ts then
(Xd2)
ti)
I (xd1 )i
if I ; I
_ Sfr = Sfr -
> tt then
(•fdi)
for j = 1 to n do
if \(tk)j\ < Ts then
L (ifc)j = 0;
dk = l-{tlw~1,uky,
14 = {S1,S2, . . . ,Sn]jT = {ti,t2, .
is P2 = W-1 — W-1,SD-1TTW-1
for j = 1 to ?? do
if |(sa.)j| < tt then
_ (sk)j = 0;
, tn} D = {di, d2j... j dn};
d
Sk tk
не возникает ситуации блокировки памяти, что позволило выполнить операции матрично-векторжнх) и скалярных произведений (строки 4, 15) независимо в параллельных нитях. Такой подход применяется при вычислениях на нескольких GPU. В этом случае каждая последующая итерация цикла зависит от данных, полученных на предыдущем uiai'e, и требуется выполнение Sk tk
Шаги алгоритма, содержащие матричные операции (строка 15), были распараллелены в рамках технологии CUDA с разработкой ядра, в котором произведение матриц вычисляется в виде последовательности матрично-векторных произведений. Эффективность распараллеливания данного этана очень высокая, и затраты существенно меньше, чем на предыдущих шагах алгоритма.
При построении предобуелавливателя Paism разреженность матриц неизвестна. Промежу-
Sk tk
столбцами матриц S, Т, хранящихся по строкам. Расходы по памяти увеличивались, но сокращалось время обращения к элементам векторов. Для преобразования из сжатого формата хранения матриц CSR в формат хранения полных строк (этап построения предобуелавливателя) и обратное преобразование (матрично-векторное произведение при решении СЛАУ) были разработаны эффективные параллельные алгоритмы, позволяющие пренебречь затратами на преобразование матриц [6].
§ 3. Численные эксперименты
В численных экспериментах были использованы матрицы из коллекции The University of Florida Sparse Matrix Collection. Решались системы уравнений Лу = / с известным точным решением у = [1,1,... , 1], матрицы которых хранились в сжатом строчном формате (CSR). В качестве начального приближения выбиралось уо = [0, 0,... , 0], а критерий сходимости -11ri|| ^ 10-6||го||, где п = / — Луг.
В расчетах рассматривался один из вариантов представления матрицы W = в-diag (Л), хотя возможны и другие, например W = в1- Выбор стратегии фильтрации по значениям элементов и предельных значений т основывался на рекомендациях, приведенных в работе [7]. Отметим, что затраты на формирование возрастают несущественно при уменьшении порогового значе-
ния. Выбор оптимального значения параметра т для рассматриваемых матриц заключался в минимизации числа итераций и времени решения. Точность фильтрации для матриц S и T выбиралась следующая: ти = 0.0001 в = 100.
Сравним сначала результаты, полученные при решении несимметричных систем. В таблицах приведены затраты на построение предобуславливателей неполного разложения с контролем заполнения ILU(p) (рассматривались варианты p = 0 и p = 1), явном вычислении приближенной обратной матрицы AINV, TNS в реализации пакета PARALUTION (http: //www.paralution.com/) на центральных процессорах и графических ускорителях.
Предобуславливатели ILU(p) формировались на центральном процессоре, а итерационный процесс решения BiCGStab выполнялся на ускорителе вычислений. В таблице 1 приведены времена формирования предобуславливателя (tp), итерационного процесса (tits) и число итераций (its), необходимых для решения систем линейных уравнений с использованием предобуславливателей LU и AISM. Временные затраты на пересылку предобуславливателя, сформированного на CPU, на графический ускоритель в приведенных табличных данных не учитывались.
Таблица 1. Вычислительные затраты предобуславливателей при решении систем с несимметричными матрицами, tp/tits(its)
Матрица PILU (0) PILU (1) Paism (t = 0.0001)
A CPU/GPU CPU/GPU GPU/GPU 2 x GPU/GPU
eddeS 0.0002/0.13 (140) 0.002/0.12 (106) 0.57/0.13(167) 0.76/0.13 (165)
ех.37 0.003/0.03 (3) 1.4/0.03 (2) 15.44/0.004 (4) 8.76/0.004 (4)
raj at 03 169/0.11 (135) 79.08/0.28 (333)
flowmeters 0.002/0.36 (63) 0.04/0.34 (32) 357.9/0.15 (112) 165.41/0.14 (105)
ex 19 0.04/ 699/0.5 (279) 703/0.39 (128) 325.96/0.44 (127)
sme3Da 0.15/13.23(1700) 495/16.61(158) 843/13.5(1338) 389.113/11.5 (1005)
poisson.3Da 0.04/0.13 (24) 56.3/0.58 (11) 1066/0.24 (30) 495.475/0.24 (30)
При рассмотрении результатов экспериментов основное внимание будет уделяться вычислительным затратам. Прежде всего остановимся на сравнении вариантов распараллеливания при формированиии предобуславливателя Paism- Для алгоритма 2 при использовании двух
ST
случае коммуникационные затраты при обмене векторами Xdi ,Xhi между GPU через общую память в рамках ОрепМР существенны только для небольших матриц до N = 3665 (cdde5, ех.37).
По времени работы алгоритмов решения систем методом BiCGStab видны преимущества предобуславливателя Paism• Затраты па одну итерацию в этом случае существенно ниже, чем при ILU(p). Достигнутое ускорение при формировании явного предобуславливателя Paism все-таки сохраняет достаточно большие вычислительные затраты и требует дальнейших исследований.
При рассмотрении плохо обусловленных матриц больших) размера (ех19, sme3Da) затраты па формирование предобуславливателя Paism сравнимы с вариантом Pilu(1)- Как видно из таблицы 1, в ряде случаев предобуславливатели на основе неполного разложения не обеспечили сходимость к решению. Так, использование ILU(O) в случае матрицы ех19 не привело к решению системы, а при решении системы с матрицей rajat03 из приведенных предобуславливателей только AISM обеспечил сходимость к точному решению.
Исключительно для целей тестирования в таблице 2 приведены результаты для явных предобуславливателей AINV, TNS при решении симметричных систем уравнений методом сопряженных градиентов с формированием матрицы предобуславливателя на центральных процес-
еорах и графических ускорителях. В этом случае симметричность матриц при формировании предобуславливателя Paism не учитывалась. Потенциально сокращение затрат при формировании раеематриваемш'о предобуславливателя в случае симметричных матриц представляется возможным и перспективным.
По скорости сходимости результаты многих тестов для различных предобуславливателей сопоставимы, а для матриц bcsstklS, vibrobox применение Paism позволило получить меньшее число итераций. Затраты на решение систем уравнений с предобуславливателями Paism и Painv примерно одинаковы.
Дальнейшее развитие предобуславливателей, основанных на приближенном обращении Шермана Морриеона, должно быть связано блочным представлением разложения, выполняемом на нескольких GPU.
Таблица 2. Вычислительные затраты явных предобуславливателей при решении систем с симметричными матрицами, tp/tits(its)
Матрица PAINV PTNS PAISM
A CPU/GPU GPU/GPU GPU/GPU 2x GPU/GPU
nasa2910 2.5/0.65 (314) 0.002/0.32 (962) 8.78/0.62 (387) 5.62/0.57 (339)
bcsstklS 0.15/0.6 (293) 0.002/0.08 (259) 20.37/0.17 (81) 11.23/0.16 (70)
Kuu 4.03/0.16 (75) 0.004/0.1 (241) 142.03/0.18 (103) 66.7/0.18 (102)
mse.10848 1.48/2.68 (1190) 0.006/14.7 (21871) 505.5/5.12 (846) 230.715/7.03 (757)
vibrobox 1.29/1.42 (683) 0.003/1.98 (4682) 814/0.81 (52) 391.8/0.8 (52)
Список литературы
1. Saad Y. Iterative methods for sparse linear systems. SIAM. 2003. xviii — 528 p. ISBN-10: 0-89871-534-2
2. Berizi M. Preconditioning techniques for large linear systems: a survey // Journal of Computational Physics. 2002. Vol. 182. № 2. P. 418 477.
3. Kolotilina L.Yu., Yeremin A.Yu. Factorized sparse approximate inverse preconditionings. I: Theory // SIAM Journal on Matrix Analysis and Applications. 1993. Vol. 14. №1. P. 45 58.
4. Grote M.. Hnckle T. Parallel preconditioning with sparse approximate inverses // SIAM Journal on Scientific Computing. 1997. Vol. 18. №3. P. 838 853. DOI: 10.1137/S1064827594276552
5. Sherman J.. Morrison W.J. Adjustment of an inverse matrix corresponding to a change in one element of a given matrix // The Annals of Mathematical Statistics. 1950. Vol. 21. № 1. P. 124 127.
6. Недожогин H.C.. Копысов С.П.. Новиков А.К. Параллельное формирование предобуславливателя. основанного на аппроксимации обращения Шермана Морриеона /'/' Вычислительные методы и программирование. 2015. Т. 16. С. 86 93.
7. Brn М.. Cerdan J.. Marin J.. Mas J. Preconditioning sparse nonsymmetric linear systems with Sherman Morrison formula /'/' SIAM Journal on Scientific Computing. 2003. Vol. 25. №2. P. 701 715.
DOI: 10.1137/S1064827502407524
Поступила в редакцию 30.09.2015
Недожогин Никита Сергеевич, младший научный сотрудник. Институт механики УрО РАН. 426067. Россия, г. Ижевск, ул. Т. Барамзиной. 34. E-mail: nedozhoginiiinbox.rn
Копысов Сергей Петрович, д. ф.-м. н.. профессор, кафедра вычислительной механики. Удмуртский государственный университет. 426034. Россия, г. Ижевск, ул. Университетская. 1: заведующий лабораторией. Институт механики УрО РАН. 426067. Россия, г. Ижевск, ул. Т. Барамзиной. 34. E-mail: s.коруsov¡йgmail.com
Новиков Александр Константинович, к. ф.-м. н.. старший научный сотрудник. Институт механики УрО РАН. 426067. Россия, г. Ижевск, ул. Т. Барамзиной. 34. E-mail: sc_work:<imail.ru
N. S. Nedozhogin, S.P. Kopysov, A. K. Novikov Parallel decomposition version in the AISM preconditioner
Keywords: solving systems of linear algebraic equations, preconditioners, parallel algorithms, graphics processing.
MSC: 15A23, 65F08
We propose a variant of the parallel decomposition in the formation of preconditioner based on the approximate inversion of the Sherman-Morrison. We conduct numerical experiments for solving test systems of equations using GPUs.
REFERENCES
1. Saad Y. Iterative methods for sparse linear systems, SIAM, 2003. xviii + 528 p. ISBN-10: 0-89871-534-2
2. Benzi M. Preconditioning techniques for large linear systems: a survey, Journal of Computational Physics, 2002, vol. 182, no. 2, pp. 418-477.
3. Kolotilina L.Yu., Yeremin A.Yu. Factorized sparse approximate inverse preconditionings. I: Theory, SIAM Journal on Matrix Analysis and Applications, 1993, vol. 14, no. 1, pp. 45-58.
4. Grote M., Huckle T. Parallel preconditioning with sparse approximate inverses, SIAM Journal on Scientific Computing, 1997, vol. 18, no. 3, pp. 838-853. DOI: 10.1137/S1064827594276552
5. Sherman J., Morrison W.J. Adjustment of an inverse matrix corresponding to a change in one element of a given matrix, The Annals of Mathematical Statistics, 1950, vol. 21, no 1, pp. 124-127.
6. Nedozhogin N.S., Kopysov S.P., Novikov A.K. Parallel forming preconditioner based on the approximate Sherman-Morrison inversion, Vychislitel'nye metody i programmirovanie, 2015, vol. 16, pp. 86-93 (in Russian).
7. Bru M., Cerdan J., Marin J., Mas J. Preconditioning sparse nonsymmetric linear systems with Sherman-Morrison formula, SIAM Journal on Scientific Computing, 2003, vol. 25, no. 2, pp. 701-715.
DOI: 10.1137/S1064827502407524
Received 30.09.2015
Nedozhogin Nikita Sergeevich, Junior Researcher, Institute of Mechanics, Ural Branch of RAS, ul. T. Baramzi-noi, 34, Izhevsk, 426067, Russia. E-mail: nedozhogin@inbox.ru
Kopysov Sergei Petrovich, Doctor of Physics and Mathematics, Professor, Department of Computational Mechanics, Udmurt State University, ul. Universitetskaya, 1, Izhevsk, 426034, Russia; Head of the Laboratory, Institute of Mechanics, Ural Branch of RAS, ul. T. Baramzinoi, 34, Izhevsk, 426067, Russia. E-mail: s.kopysov@gmail.com
Novikov Aleksandr Konstantinovich, Candidate of Physics and Mathematics, Associate Professor, Department of Computational Mechanics, Udmurt State University, ul. Universitetskaya, 1, Izhevsk, 426034, Russia; Senior Researcher, Institute of Mechanics, Ural Branch of RAS, ul. T. Baramzinoi, 34, Izhevsk, 426067, Russia.
E-mail: sc_work@mail.ru
