3_19

отличием что в (4.22)–(4.23) векторы не нормируются. Таким образом, в симметричном случае ортогонализация Арнольди сводится к трехчленному рекуррентному соотношению

называемомуортогонализациейЛанцоша. Коэффициенты i, i образуют трехдиагональную матрицу Tm, для которой справедливо соотношение (аналог (4.20))

Этот же факт можно получить и другим способом. Воспользуемся соотношением (4.20). Если A = AT , òî VmT AVm = = (VmT AVm)T . Тогда

HTm = (VmT AVm)T = VmT AVm = Hm :

Очевидно, что матрица в форме Хессенберга может быть симметрич- ной лишь тогда, когда она является трехдиагональной.

Рассмотренный в §5.1 метод FOM может быть упрощен для симметричного случая, если процедуру Арнольди в нем заменить на (6.2); для решения трехдиагональной СЛАУ при этом наиболее естественно использовать метод прогонки [4, 3]. Такая вычислительная схема, называемаяметодом Ланцоша, приведена на рис. 6.1.

6.2. CG: Метод сопряженных градиентов

Как уже отмечалось, в случае симметричной матрицы пространства Крылова Km(r0; A) è Km(r0; AT ) совпадают. Это приводит к тому, что в методе бисопряженных градиентов (рис. 5.4) будут иметь место равенства

r~i ri ; p~i pi

при выборе r~0 = r0.

61

r0 := b Ax0

:= kr0k2; v1 := r0=1 := 0; v0 := 0

Äëÿ j = 1; m

wj := Avj j vj 1

j := (wj ; vj )

j+1 := kwj k2

Åñëè j+1 = 0

vj+1 := wj = j+1

увеличить j

Tm := tridiag( i; i; i+1)mi=1

Найти y из трехдиагональной системы Tmy = e1

m

P

xm := x0 + yivi

i=1

Рис. 6.1. Метод Ланцоша

Очевидное упрощение вычислительной схемы сводится к тому, что отпадает необходимость в хранении и обработке векторов r~i è p~i; соответствующий алгоритм приведен на рис. 6.2.

Соотношение (5.14) при этом превращается в условие

Определение 6.2.1 Векторы fpigmi=1,удовлетворяющиеусловию(6.4),

называются A-сопряженными или, если не возникает неоднозначности, просто сопряженными.

Соответственно, метод рис. 6.2 называется методом сопряженных градиентов èëè CG1.

Если в дополнение к симметричности потребовать, чтобы матрица A была положительно определенной, то при pi 6= 0 всегда pTi Api =6 0 и сходимость метода гарантирована.

1Conjugate Gradient Method.

62

Выбрать начальное приближение x0 r0 := b Ax0

p0 := r0

Äëÿ j = 1; 2; : : :

rj+1 := rj j Apj

j := (rj+1; rj+1) (rj ; rj )

Åñëè krj+1k2 < "

òî КОНЕЦ (решение достигнуто)

pj+1 := rj+1 + j pj

увеличить j

Ðèñ. 6.2. CG

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

При этом CG относится к классу тех проекционных методов, для которых K = L (= Km(r0; A)). Теорема 4.3.1 утверждает, что задача проектирования в этом случае эквивалентна минимизации функционала 1(x) = kx x k22, а в примере 1 из §4.2 было показано, что направление его наискорейшего убывания r 1 совпадает с направлением невязки.

Метод бисопряженных градиентов исторически появился позднее, как обобщение CG на несимметричный случай, что и отразилось на его названии, так как вместо ортогонализации (6.2) в нем используется биортогонализация (4.22)–(4.25).

Замечание. В [14] показан пример подхода к выводу метода сопряженных градиентов, не использующего идеологию LU-разложения.

2Ò.å. A-ортогонализацией в скалярном произведении (3).

63

6.3.О связи симметричного и несимметрич- ного случая

Рассмотрение симметричных СЛАУ завершим табл. 6.1, кратко обобщающей рассуждения двух предыдущих параграфов.

Замечание. Метод Ланцоша и метод сопряженных градиентов эквивалентны алгебраически. Однако их реализации в конечной арифметике будут давать, вообще говоря, разные результаты.

Таблица 6.1. Связь симметричного и несимметричного случая

64

Заключение

Поскольку решение СЛАУ является стандартным этапом численного моделирования, для его реализации разработано большое количе- ство программного обеспечения, большинство которого является бесплатным и, ведя свою историю еще от появления первых электронновычислительных машин, стало стандартом de-facto.

Для систем с плотными матрицами следует упомянуть такие пакеты, как BLAS, LINPACK, EISPACK и LAPACK; для СЛАУ с разреженными матрицами существует пакет SRARSKIT. Перечисленные пакеты разработаны на языке FORTRAN и хотя не реализуют собственно методы решения СЛАУ, но содержат богатый набор подпрограмм, описывающих всевозможные операции с матрицами и векторами, часто встречающиеся в схемах методов. Наличие таких библиотек и их высокая оптимизированность значительно упрощает и делает более эффективной разработку программ-решателей.

Кроме того, авторами данного пособия разработан программный пакет LinPar, содержащий готовые реализации ранее описанных методов CG, BCG, GMRES(m) и методов класса А.А.Абрамова (ориентированных на СЛАУ с плотными матрицами) для платформы IntelDOS/Windows.

Internet-адреса, по которым доступно перечисленное программное обеспечение, приведены в табл. 6.2.

65

Таблица 6.2. Программное обеспечение линейной алгебры

66

Приложение A

LU-факторизация

трехдиагональных матриц

Решение вспомогательных СЛАУ с трехдиагональными матрицами требуется в тех методах, которые основаны на трехчленных соотношениях для построения базиса (см. BCG, §5.4 и CG, §6.2). Как оказывается, LU-разложение таких матриц осуществляется очень простым способом, что позволяет обойтись без излишнего накопления данных.

Имеет место следующая теорема:

Теорема A.0.1 Матрицы Lm è Um,образующие LU-разложениетрех-

диагональной матрицы

являются нижнедвухдиагональной и верхнедвухдиагональной соответственно.

Доказательство. Утверждение теоремы доказывается по индукции.

67

Для матрицы третьего порядка T3 его легко проверить непосредственно.

Предположим теперь, что для Tm утверждение справедливо и она представима в виде произведения Tm = LmUm нижнедвухдиагональной Lm и верхне-двухдиагональной Um. Заметим, что m-ый столбец Lm по предположению индукции имеет лишь один — последний — ненулевой элемент, и то же самое справедливо для m-ой

строки Um.

Необходимо показать, что существуют такие коэффициенты , ,

Выполнив умножение в правой части, получим блочно-матричное уравнение

которое сводится к трем скалярным:

8

<Lm mm = bm

Um mm = cm

: + = am

Замечание. В [3] доказывается и более общий результат, заключа- ющийся в том, что ленточный характер матриц сохраняется при LUразложении для любой ширины ленты, в том числе для лент, имеющих разную ширину вверх и вниз относительно главной диагонали. Этот факт может оказываться полезным, например, для метода IOM (см. §5.1).

Коэффициенты самих матриц Lm è Um с учетом простоты их вида можно найти непосредственно из выполнения умножения LmUm è

68

сравнения результата с матрицей Tm. Так, если отнести диагональ к матрице U, то разложение принимает вид

Сопоставление этой матрицы с Tm дает следующие рекуррентные формулы:

1 = a1 ;i bi ;

i+1 = ci+1=i ;

i+1 = ai+1 bi+1 i+1 :

При построении методов BCG и CG используются матрицы Lm1 è Um1, которые легко находится аналитически

69

Приложение B

Методы, не использующие транспонирование. TFQMR

В §5.4 был описан метод бисопряженных градиентов BCG, имеющий по сравнению с FOM и GMRES гораздо более простую вычислительную схему и меньшие требования к памяти.

Однако BCG (как и любой другой метод, использующий операции с транспонированной матрицей A) плохо поддается реализации на многопроцессорных вычислительных системах с распределенной памятью. Все более широкое применение таких систем привело к разработке целого класса методов, в которых операция транспонирования не используется1.

Алгебраически это может быть достигнуто за счет изменения специальным образом полинома pm, которому (см. §4.4) удовлетворяет последовательность невязок в методах, использующих подпространства Крылова:

Так, на использовании вместо pm(A) полинома p2m(A) основан метод CGS [13]; метод BCGSTAB вместо (B.1) использует соотношение rm = pm(A)qm (A)r0, ãäå qm — специальным образом строящийся полином, такой что произведение pmqm не содержит нечетных степеней.

1Àíãë. transpose-free methods.

70