Решение систем линейных алгебраических уравнений методом BiCGStab с предобусловливанием Текст научной статьи по специальности «Математика»

УДК 519.612.2

В. В. Пузикова

РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ МЕТОДОМ BICGSTAB С ПРЕДОБУСЛОВЛИВАНИЕМ

Рассмотрен стабилизированный метод бисопряженных градиентов (BiCGStab) и его модификации с предобусловливанием. Метод BiCGStab позволяет решать системы линейных алгебраических уравнений достаточно общего вида: симметрия и положительная определенность матрицы системы не имеют принципиального значения. Для анализа эффективности метода проведены тестовые расчеты, в которых решены линейные системы с разреженными матрицами. Также рассмотрены варианты метода с ILU и многосеточным предобусловливанием. Результаты расчетов подтверждают высокую эффективность метода BiCGStab с многосеточным предобусловливанием.

E-mail: [email protected]

Ключевые слова: разреженные линейные системы, стабилизированный метод бисопряженных градиентов, ILU-предобусловливание, многосеточные методы, многосеточное предобусловливание.

Метод BiCGStab. Стабилизированный метод бисопряженных градиентов (BiCGStab) был предложен в 1992 г. Х. ван дер Ворстом [1]. Метод характеризуется наиболее быстрой и гладкой сходимостью (под гладкостью понимается отсутствие сильных колебаний нормы невязки) из известных методов крыловского типа. Данный метод является развитием широко распространенного метода сопряженных градиентов (CG) [2], весьма эффективного для систем с симметричными положительно определенными матрицами. Преимуществом метода BiCGStab является возможность его использования для решения систем достаточно общего вида: симметрия и положительная определенность матрицы не имеют принципиального значения. Алгоритм метода BiCGStab для решения системы Ax = b, A е M (R)N xN, x, b e RN имеет следующий вид [2]:

Vx0, r0 = b - Ax0, p0 = r0, Vr°0 : (r0,r0) = 0

for n = 0,1,..., while (||rn+1||2 > e), do

(r0,rn) (Asn,sn)

\ * > / _n „n л „n... \ >'

an = т~л П—^; sn = rn - anApn; шп =

<

(Apn,r°)' n * , n )' (1)

an(rn+1,r0)

xn+i = xn + anpn + wnsn:

rn+1 = sn - wnAsn; ßn =

Jn

pn+1 = rn+1 + en(pn - wnApn). Здесь rn — вектор невязки, pn — вектор коррекции на n-м шаге.

^n(rn,r°)

•

Таблица 1

Число итераций при применении различных итерационных методов

Метод Число итераций

"Тест" 1 "Тест 2"

Метод Зейделя 3313 2987

Метод верхней релаксации (w =1,4) 2384 2168

BiCGStab 139 240

Для анализа эффективности данного метода было проведено несколько тестовых расчетов, в которых решались линейные системы с разреженными матрицами. ("Тест 1" — система с симметричной матрицей. "Тест 2" — система с несимметричной матрицей.) Размерности обеих систем N = 809; число ненулевых элементов в матрицах равно М = 5359 и М = 4988 соответственно. Данные системы были получены при решении стационарной задачи теплопроводности методом конечных элементов.

Во всех расчетах полагалось х0 = {0}^, г° = г0 = р0, е = 10-16, в качестве критерия останова использовалось условие ||гп+1||2 < е. Для сравнения системы также решались методом Зейделя и методом релаксации с параметром релаксации ш = 1,4 [2]. Эффективность указанных методов оценивалась не только по числу итераций (табл. 1, рис. 1), но и по числу выполненных арифметических операций умножения и деления (табл. 2).

Таким образом, метод БЮС81аЪ является весьма эффективным численным методом решения разреженных систем; в рассмотренных

Рис. 1. Убывание норм невязок при решении СЛАУ "Тест 2" (1 — метод Зейделя, 2 — метод релаксации, 3 — ВЮС81аЬ)

Таблица 2

Число операций при применении различных итерационных методов

Метод

Число операций

в 1 итерацию "Тест 1" "Тест 2"

Метод Зейделя

M 17 754 367 14 388 379

Метод верхней релаксации (w =1.4)

M

12 775 856 10 443 256

BiCGStab

2 M + 11N 2 726 763 4 447 920

примерах он позволил сократить число итераций в 10-20 раз по сравнению с методами Зейделя и релаксации, а число арифметических опреаций (и, соответственно, время счета) в 2,5-6 раз. Использование модификаций метода БЮС81аЪ с предобусловливанием позволяет повысить эффективность метода. Далее будут рассмотрены варианты с /Ьи и многосеточным предобусловливателями.

/Ьи-предобусловливание метода В1СС81аЬ и его варианты. Представление матрицы в виде произведения нижнетреугольной и верхнетреугольной матриц (Ьи-факторизация) позволяет легко решать системы вида Ах = Ь путем двукратного решения систем Ьу = Ь и их = у с треугольными матрицами. Однако этот алгоритм непригоден для систем с разреженными матрицами, так как ведет к заполнению портрета, т.е. появлению большего числа ненулевых элементов в матрицах Ь и и по сравнению с исходной матрицей А.

Поэтому для разреженных матриц вместо задачи Ьи-факторизации формулируют задачу неполной Ьи-факторизации (/Ьи-факторизации): матрица А представляется в виде А = Ьи + Л, где Я — матрица ошибок разложения. При этом Ь и и — нижнетреугольная и верхнетреугольная матрицы соответственно; Рь С РА и Рц С РА;

V(i, j) G Pa : [LU]ij = ; PA П Pr = 0, PA, Pl, Pu, Pr - портреты,

т.е. множества позиций ненулевых элементов, матриц А, Ь, и и Я соответственно.

Приближенное представление А « Ьи, удовлетворяющее данным условиям, называется /Ьи-разложением матрицы А. Для его получения можно использовать следующие соотношения [3]:

г=1

Матрица М = Ьи (Ьи — /Ьи-факторизация матрицы А) удовлетворяет условиям, предъявляемым к матрице предобусловливателя для

(2)

(3)

системы вида Ax = b: она хорошо приближает матрицу A, поскольку на множестве индексов Pa точно воспроизводит ее, она легко вычисляется по приведенным формулам и она легко обратима, так как является произведением двух треугольных матриц. Метод BiCGStab с /LU-предобусловливанием по существу представляет собой решение системы M-1Ax = M-1b с помощью алгоритма (1); его алгоритм имеет вид

Ух0, r0 = b — Ax0, p0 = r0, Vr0 : (r0,r0) = 0 for n = 0,1,..., while (||rn+1||2 > e), do

(r° rn)

yn = M-1Pn• а = '_— • sn = rn — а Avn-

y M p • an (Ayn,r0); s ' anAy •

(Azn,Sn) (4)

zn = M-1 sn; шп = ^-Ц-Л;

' n (Azn,Azn)

xn+1 = xn + anyn + wnzn; rn+1 = sn — wnAzn; a (rn+1 r0)

Pn = nV( n '0:;; pn+1 = r^1 + ¡3n(pn — UnAyn).

Отметим, что предобусловленный алгоритм не требует в явном виде обращения матриц и их умножения на векторы: необходимо однократно вычислить матрицы L и U, а выполнение операций алгоритма yn = M-1 pn и zn = M-1 sn сводится к решению двух систем с треугольными матрицами.

Автор метода BiCGStab, Х. ван дер Ворст, обращает внимание на то, что при таком вычислении шп предобусловливание несколько снижает скорость сходимости и норма невязки для системы M-1Ax = M-1b убывает медленнее. Для исходной системы Ax = b [4]. Можно показать, что для системы M-1Ax = M-1b норма невязки будет убывать быстрее, чем для исходной, если вычислять коэффициенты un по следующей формуле:

(L-1Az n,L-1sn) Шп = (L-1Azn,L-1Azn)'

Данную модификацию /LU-предобусловливания будем называть модификацией ван дер Ворста. Еще одна модификация /LU-предобусловливания — S/P(а)-предобусловливание [5]. Она получается при использовании в /LU-разложении параметра а, при этом S/P(а)-предобусловливание эквивалентно /LU-предобусловливанию. Однако такой подход применим только для разреженных матриц из трех, пяти, семи или девяти диагоналей. На разреженные матрицы общего вида алгоритм построения такого разложения не обобщен. В настоящей работе предпринята попытка ввести в /LU-разложение для

матриц общего вида параметр а в целях ускорения сходимости. При а = 0 оно, как и SIP (а), превращается в ILU-предобусловливание.

Поскольку ILU-разложение неединственно, можно допустить любую модификацию формул для вычисления элементов матрицы L или U (при вычислении элементов второй матрицы по формулам (2, 3) она будет скорректирована так, чтобы произведение по-прежнему оставалось неполным ILU-разложением). После некоторых методических экспериментов для вычисления элементов матрицы U была выбрана формула

к-1

Uk-j = a,к - (1 + а) l

kj Uij

i=l

Элементы матрицы Ь вычисляются по формуле (2). Такое разложение для краткости будем называть а1Ьи-разложением. Для сравнения описанных алгоритмов были рассмотрены те же тестовые системы ("Тест 1" и "Тест 2") и те же начальные значения параметров алгоритмов, что и ранее. При использовании метода с а1Ьи-предобусловливанием оптимальное значение а определялось экспериментально. Рис. 2 иллюстрирует убывание нормы невязки для системы "Тест 2" при использовании метода БЮС81аЪ с различными вариантами предобусло-вливания. В табл. 3 и 4 приведены сведения о числе итераций, затраченных на решение систем, и числе выполненных арифметических операций; при использовании а1Ьи-предобусловливания в скобках указано оптимальное значение параметра а.

Видно, что использование вариантов метода БЮС81аЪ с предобусловливанием увеличило скорость сходимости в среднем в 3,5 раза по сравнению с БЮС81аЪ без предобусловливания. По числу вы-

Рис. 2. Убывание норм невязок при решении СЛАУ "Тест 2" (1 — Б1СС81аЬ, 2 — варианты метода БЮС81аЬ с предобусловливанием)

Таблица 3

Число итераций при применении различных итерационных методов

Метод Число итераций

"Тест 1" "Тест 2"

BiCGStab 139 240

BiCGStab+ILU 41 30

BiCGStab+ILU (ван дер Ворст) 40 30

BiCGStab+oJXU 37 27

BiCGStab+oJXU (ван дер Ворст) 37 27

Таблица 4

Число операций при применении различных итерационных методов

Число операций

Метод в 1 итерацию "Тест 1" "Тест 2"

BiCGStab 2 M + 11N 2 726 763 4 447 920

BiCGStab+ILU 4M + 9N 1 177 397 796 470

BiCGStab+JLU (ван дер Ворст) 5 M + 7N 1 298 320 892 440

BiCGStab+alLU 4M + 9N 1 062 529 716 823

BiCGStab+alLU (ван дер Ворст) 5 M + 7N 1 200 946 803 196

полненных умножении предпочтительными также являются варианты BiCGStab с предобусловливанием.

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

Ошибка zn решения системы уравнений на n-й итерации может быть представлена в виде разложения с некоторыми коэффициентами СП £ (R) по собственным векторам соответствующей разностной задачи Штурма-Лиувилля, которые в случае однородных граничных условий имеют вид

N — 1

sin N-a £ rn.

) i=0

Показано, что при решении системы различными итерационными методами, такими как методы Якоби и Зейделя, коэффициенты сП, соответствующие большим номерам k (k > 0,5N), уменьшаются достаточно быстро, в то время как коэффициенты при низких гармониках (k < 0,5N) убывают сравнительно медленно [6]. Поэтому после

Рис.3. Гармоника , к = 2 на подробной = 7), (а) и грубой = 3), (б) сетках

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

Шаг 1. Несколько сглаживающих итераций. В настоящей работе используется ЛБЬ//-метод (метод релаксированной переменной линии Якоби) [7].

Шаг 2. Несколько итераций для подавления низкочастотных компонент в ошибке.

Шаг 3. Если требуемая точность не достигнута, перейти к шагу 1.

Рассмотрим алгоритм подавления низкочастотных компонент в ошибке. Заметим, что половина низкочастотных гармоник отображается в высокочастотные гармоники на более грубой сетке. Это показано на рис. 3. Слева приведена гармоника, соответствующая к = 2 на подробной сетке (Ы = 7). На ней эта гармоника является низкочастотной (к < 0,5Ы). Справа представлена проекция этой же гармоники на сетку, которая в 2 раза грубее исходной (Ы = 3), и на ней гармоника является уже высокочастотной (к > 0,5Ы).

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

Таким образом, если ошибка "гладкая", то ее подавление будет эффективным на более грубой сетке (в теории многосеточных методов это называется коррекцией с грубой сетки). При этом на сетке, которая в 2 раза грубее исходной, вычислительные затраты на одну итерацию

1 1 будут составлять ^ в 2D-случае и ^ в 3D-случае от вычислительных

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

Зейделя при решении двумерных задач сходится в 4 раза быстрее на сетке, являющейся в 2 раза более грубой).

Это указывается на то, что большая часть работы может быть проделана на более грубых сетках. Для этого требуется определить соотношения между двумя сетками, оператор на грубой сетке, метод сужения невязки с подробной сетки на грубую и метод пролонгации коррекции с грубой сетки на подробную. Если число ячеек подробной сетки в каждом направлении четно, то грубая сетка состоит из каждой второй линии подробной сетки, а в качестве метода пролонгации используется обычная линейная интерполяция.

Пусть P — матрица оператора пролонгации сеточной функции с грубой сетки на подробную, R — матрица оператора сужения сеточной функции с подробной сетки на грубую, индекс h указывает на принадлежность подробной сетке, а индекс H — на принадлежность грубой сетке.

Пусть на подробной сетке нужно решить систему Ahuh = fh. Тогда на грубой сетке ей будет соответствовать система AHuh = fH, где uh = Ruh, fH = Rfh. Матрица системы на грубой сетке определяется с помощью подхода Галеркина [8]: AH = RAhP.

Алгоритм метода BiCGStab с многосеточным предобусловливате-лем аналогичен алгоритму метода BiCGStab с ILU-предобусловлива-телем (4) с той лишь разницей, что две системы уравнений, которые необходимо решать на каждой итерации этого метода, решаются не с помощью ILU-разложения матрицы системы, а многосеточным методом (запись MG(Ax = b) в приведенном далее алгоритме означает, что система Ax = b решается многосеточным методом):

Vx0, r0 = b - Ax0, p0 = r0, Vr0 : (r°°,r0) = 0 for n = 0,1,..., while (\\rn+1\\2 > e), do

(r0 r'n)

MG(Ayn = pn); an = Д*' 0ч; sn = rn - anAyn; MG(Azn = sn);

(Ay ,r*)

(Azn,sn)

Wn —

(Лгп,Лгп)'

хп+1 = хп + апуп + Шпхп; тп+1 = вп — шпЛгп;

а (тп+1 то\

вп = п\ п ; рп+1 = тп+ + вп(рп — ШпЛуп). Мп[.т ,т* )

На рис. 4 и 5 показаны результаты подсчета числа итераций, выполненных при решении разностного аналога уравнения Пуассона в задаче об обтекании квадратного препятствия методом БЮС81аЪ с 1Ы7-предобусловливанием и многосеточным предобусловливанием соот-

Рис. 4. Решение разностного аналога уравнения Пуассона методом BiCGStab с 1Ы1-предобусловливанием

Рис. 5. Решение разностного аналога уравнения Пуассона методом BiCGStab с многосеточным предобусловли-ванием

ветственно. Задачи решались на неравномерной сетке 120 х 148 с шагом по времени 0,01, ограничение на норму невязки е = 10-6, системы линейных уравнений получались путем LS-STAG-дискретизации уравнений гидродинамики [9]. Отметим, что при использовании многосеточного предобусловливания время расчета уменьшилось в 2 раза. Этот эффект планируется усилить, реализовав решатель в параллельной среде.

Заключение. Метод BiCGStab позволяет эффективно решать большие системы линейных алгебраических уравнений с разреженными матрицами, не требуя их симметрии и положительной определенности. Использование вариантов метода с 1Ьи и многосеточным пре-добусловливанием дает возможность дополнительно снизить затраты вычислительных ресурсов в несколько раз. Данные методы представляются весьма эффективными при решении задач вычислительной гидродинамики. В дальнейшем планируется разработать параллельную реализацию рассмотренных алгоритмов и использовать полученные результаты при решении практических задач.

СПИСОК ЛИТЕРАТУРЫ

1.VanderVorstH. A. Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for solution of non-symmetric linear systems// SIAM J. Sci. Stat. Comp. -1992. - No. 2. - P. 631-644.

2. S a a d Y. Iterative Methods for Sparse Linear Systems. - N.Y.: PWS Publ., 1996. -547 p.

3.Баландин М. Ю., Шурина Э. П. Методы решения СЛАУ большой размерности. - Новосибирск: Изд-во НГТУ, 2000. - 70 с.

4. VanderVorst H. A. High performance preconditioning // SIAM J. Sci. Stat. Comp. - 1989. -No. 6. - P. 1174-1185.

5. Stone H. L. Iterative solution of implicit approximations of multidimensional partial differential equations // SIAM J. Sci. Numer. Anal. - 1968. - No. 3. - P. 530558.

6. Ольшанский М. А. Лекции и упражнения по многосеточным методам. - М.: Изд-во ЦПИ при механико-математическом факультете МГУ им. М.В. Ломоносова, 2003. - 163 с.

7. WesselingP. An introduction to multigrid methods. - Chichester: John Willey & Sons Ltd., 1991.-284 p.

8. V a n K a n J., Vu i k C., We s s e l i n g P. Fast pressure calculation for 2D and 3D time dependent incompressible flow// Numer. Linear Algebra Appl. - 2000. -No. 7. - P. 429-447.

9. C h e n y Y., B o t e 11 a O. The LS-STAG method: A new immersed boundary/level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties // J. Comput. Phys. - 2010. -No. 229. - P. 1043-1076.

Статья поступила в редакцию 25.10.2011