Баландин М.Ю., Шурина Э.П. - Методы решения СЛАУ большой размерности (1051191), страница 8
Текст из файла (страница 8)
в §6.2.Biconjugate Gradients. Иногда также употребляется аббревиатура BiCG.51Глава 6Симметричный случайВсе ранее описанные методы при построении не опирались на какиелибо предположения о симметричности матрицы A, а потому безспециальных изменений могут применены к симметричным СЛАУ.Однако свойство симметрии матрицы позволяет упростить вычислительную схему, что ведет к меньшим требованиям к вычислительным ресурсам при программной реализации.Как будет показано ниже, основные упрощения могут быть внесены в процедуру построения базиса пространства Крылова.6.1.Трехчленные соотношения для невязок.Метод ЛанцошаЗаметим, что при симметричности матрицы A пространства Крылова Km (v1 , A) и Km (w1 , AT ) совпадают, если v1 = w1 .
Следовательно, вэтом случае биортогонализация Ланцоша (§4.6) приведет к построеmнию одинаковых систем векторов {vi }mi=1 , {wi }i=1 . Условие биортогональности тогда сведется кi 6= j =⇒ (vi , vj ) = (wi , wj ) = 0 ,(6.1)т.е. будет построен ортогональный базис пространства Крылова.Сравнивая формулы ортогонализации Арнольди из §4.5 и формулы (4.22)–(4.25), нетрудно заметить, что они совпадают с тем лишьотличием что в (4.22)–(4.23) векторы не нормируются.
Таким образом, в симметричном случае ортогонализация Арнольди сводится ктрехчленному рекуррентному соотношениюβj+1 vj+1 = Avj − αj vj − βj vj−1 ,(6.2)называемому ортогонализацией Ланцоша. Коэффициенты α i , βi образуют трехдиагональную матрицу T m , для которой справедливо со52r0 := b − Ax0β := kr0 k2 ; v1 := r0 /ββ1 := 0; v0 := 0Для j = 1, mwj := Avj − βj vj−1αj := (wj , vj )βj+1 := kwj k2Если βj+1 = 0тоm := jПрервать цикл по jvj+1 := wj /βj+1увеличить jTm := tridiag(βi , αi , βi+1 )mi=1Найти y из трехдиагональной системы T m y = βe1xm := x0 +mPi=1yi viРис.
6.1. Метод Ланцошаотношение (аналог (4.20))TVmAVm = Tm .(6.3)Этот же факт можно получить и другим способом. ВоспользуемT AVT Tся соотношением (4.20). Если A = AT , то Vmm = Vm A Vm =TT= (Vm AVm ) . ТогдаTTHTm = (VmAVm )T = VmAVm = Hm .Очевидно, что матрица в форме Хессенберга может быть симметричной лишь тогда, когда она является трехдиагональной.Рассмотренный в §5.1 метод FOM может быть упрощен для симметричного случая, если процедуру Арнольди в нем заменить на(6.2); для решения трехдиагональной СЛАУ при этом наиболее естественно использовать метод прогонки [4, 3]. Такая вычислительнаясхема, называемая методом Ланцоша, приведена на рис. 6.1.6.2.CG: Метод сопряженных градиентовКак уже отмечалось, в случае симметричной матрицы пространстваКрылова Km (r0 , A) и Km (r0 , AT ) совпадают.
Это приводит к тому, что53Выбрать начальное приближение x0r0 := b − Ax0p0 := r0Для j = 1, 2, . . .αj := (rj , rj ) / (Apj , pj )xj+1 := xj + αj pjrj+1 := rj − αj Apjβj := (rj+1 , rj+1 ) / (rj , rj )Если krj+1 k2 < εто КОНЕЦ (решение достигнуто)pj+1 := rj+1 + βj pjувеличить jРис. 6.2. CGв методе бисопряженных градиентов (рис. 5.4) будут иметь место равенстваr̃i ≡ ri ,p̃i ≡ piпри выборе r̃0 = r0 .Очевидное упрощение вычислительной схемы сводится к тому,что отпадает необходимость в хранении и обработке векторов r̃ i и p̃i ;соответствующий алгоритм приведен на рис. 6.2.Соотношение (5.14) при этом превращается в условиеi 6= j =⇒ pTiApj = 0 .(6.4)Определение 6.2.1 Векторы {pi }mi=1 , удовлетворяющие условию(6.4), называются A-сопряженными или, если не возникает неоднозначности, просто сопряженными.Соответственно, метод рис.
6.2 называется методом сопряженных градиентов или CG1 .Если в дополнение к симметричности потребовать, чтобы матрица A была положительно определенной, то при p i 6= 0 всегда pTiApi 6== 0 и сходимость метода гарантирована.Присутствие в названии слова «градиент» обусловлено следующей причиной. Поправки к решению вычисляются вдоль векторовpi , которые генерируются из последовательности невязок процессомих приведения к сопряженной системе 2 .При этом CG относится к классу тех проекционных методов, длякоторых K = L (= Km (r0 , A)).
Теорема 4.3.1 утверждает, что задача12Conjugate Gradient Method.Т.е. A-ортогонализацией в скалярном произведении (3).54проектирования в этом случае эквивалентна минимизации функционала Φ1 (x) = kx − x∗ k22 , а в примере 1 из §4.2 было показано, чтонаправление его наискорейшего убывания −∇Φ 1 совпадает с направлением невязки.Метод бисопряженных градиентов исторически появился позднее, как обобщение CG на несимметричный случай, что и отразилось на его названии, так как вместо ортогонализации (6.2) в нем используется биортогонализация (4.22)–(4.25).Замечание.
В [14] показан пример подхода к выводу метода сопряженных градиентов, в котором не используется идеология LUразложения.6.3.О связи симметричного и несимметричногослучаяРассмотрение симметричных СЛАУ завершим табл. 6.1, краткообобщающей рассуждения двух предыдущих параграфов.Замечание. Метод Ланцоша и метод сопряженных градиентовэквивалентны алгебраически. Однако их реализации в конечнойарифметике будут давать, вообще говоря, разные результаты.Несимметричный случайОртогонализация АрнольдиМатрица в форме ХессенбергаБиортогонализация ЛанцошаСимметричный случайОртогонализация ЛанцошаТрехчленное соотношениеТрехдиагональная матрицаОртогонализация ЛанцошаLLT или LDLT -разложениеМетод ЛанцошаМетод CGМетод CGLU-разложениеМетод FOMМетод BCG с r̃0 = r0Таблица 6.1.
Связь симметричного и несимметричного случая55ЗаключениеПоскольку решение СЛАУ является стандартным этапом численного моделирования, для его реализации разработано большое количество программного обеспечения, большинство которого является бесплатным и, ведя свою историю еще от появления первыхэлектронно-вычислительных машин, стало стандартом de-facto.Для систем с плотными матрицами следует упомянуть такие пакеты, как BLAS, LINPACK, EISPACK и LAPACK; для СЛАУ с разреженными матрицами существует пакет SRARSKIT.
Перечисленные пакеты разработаны на языке FORTRAN и хотя не реализуют собственнометоды решения СЛАУ, но содержат богатый набор подпрограмм,описывающих всевозможные операции с матрицами и векторами,часто встречающиеся в схемах методов. Наличие таких библиотеки их высокая оптимизированность значительно упрощает и делаетболее эффективной разработку программ-решателей.Кроме того, авторами данного пособия разработан программныйпакет LinPar, содержащий готовые реализации ранее описанных методов CG, BCG, GMRES(m) и методов класса А.А.Абрамова (ориентированных на СЛАУ с плотными матрицами) для платформы IntelDOS/Windows.Internet-адреса, по которым доступно перечисленное программное обеспечение, приведены в табл.
6.2.ПакетыBLAS,LINPACK,EISPACK,LAPACKSPARSKITLinParInternet-адресаhttp://www.netlib.orghttp://www.cs.umn.edu/Research/arpa/SPARSKIT/http://www.ict.nsc.ru/linpar/Таблица 6.2. Программное обеспечение линейной алгебры56Приложение ALU-факторизациятрехдиагональных матрицРешение вспомогательных СЛАУ с трехдиагональными матрицамитребуется в тех методах, которые основаны на трехчленных соотношениях для построения базиса (см.
BCG, §5.4 и CG, §6.2). Как оказывается, LU-разложение таких матриц осуществляется очень простым способом, что позволяет обойтись без излишнего накопленияданных.Имеет место следующая теорема:Теорема A.0.1 Матрицы Lm и Um , образующие LU-разложениетрехдиагональной матрицыTma1 b2 c 2 a2=c30b3a3.........cm0bmam,являются нижнедвухдиагональной и верхнедвухдиагональной соответственно.Доказательство. Утверждение теоремы доказывается по индукции.Для матрицы третьего порядка T3 его легко проверить непосредственно.Предположим теперь, что для Tm утверждение справедливо иона представима в виде произведения T m = Lm Um нижне-двухдиагональной Lm и верхне-двухдиагональной Um . Заметим, что m-ыйстолбец Lm по предположению индукции имеет лишь один — последний — ненулевой элемент, и то же самое справедливо для m-ойстроки Um .57Необходимо показать, что существуют такие коэффициенты α, β,γ, для которыхTm em bmeTm cmam!Lm 0γeTm 1=!!Um βem0α.Выполнив умножение в правой части, получим блочно-матричноеуравнениеTm em bmeTm cmam!Lm Um βLm emγeTm Um α + γβ=!,которое сводится к трем скалярным: β[Lm ]mm = bmγ[U ]= cmm mm α + γβ = am−1имеющим очевидное решение β = bm [Lm ]−1mm , γ = cm [Um ]mm , α = am −− bm cm / ([Lm ]mm [Um ]mm ).
Тем самым теорема доказана.2Замечание. В [3] доказывается и более общий результат, заключающийся в том, что ленточный характер матриц сохраняется приLU-разложении для любой ширины ленты, в том числе для лент, имеющих разную ширину вверх и вниз относительно главной диагонали. Этот факт может оказываться полезным, например, для методаIOM (см. §5.1).Коэффициенты самих матриц Lm и Um с учетом простоты их видаможно найти непосредственно из выполнения умножения L m Um исравнения результата с матрицей Tm . Так, если отнести диагональ кматрице U, то разложение принимает видLm1 γ 2=001......γm1;UmВыполнив умножение, получимLm Umα1β2 α γ α +β γ22 2 1 2=α2 γ3=β3α3 + β 3 γ3...58α1β2α20....0......βmαm...βmαm−1 γm αm + βm γm.Сопоставление этой матрицы с Tm дает следующие рекуррентныеформулы:α1 = a 1 ;βi ≡ b i ;γi+1 = ci+1 /αi ;αi+1 = ai+1 − bi+1 γi+1 .При построении методов BCG и CG используются матрицы L −1m икоторые легко находится аналитическиU−1m ,[L−1m ]ij =[U−1m ]ij0, 1,j>ij=ii+j (−1)0, 1/α ,i=(−1)i+jαiiQk=j+1γk ,j<ij=ijQβk,αk=i+1 k59(A.1)j<ij>i(A.2)Приложение BМетоды, не использующиетранспонирование.
TFQMRВ §5.4 был описан метод бисопряженных градиентов BCG, имеющийпо сравнению с FOM и GMRES гораздо более простую вычислительную схему и меньшие требования к памяти.Однако BCG (как и любой другой метод, использующий операции с транспонированной матрицей A) плохо поддается реализациина многопроцессорных вычислительных системах с распределеннойпамятью.
Все более широкое применение таких систем привело кразработке целого класса методов, в которых операция транспонирования не используется1.Алгебраически это может быть достигнуто за счет изменения специальным образом полинома pm , которому (см. §4.4) удовлетворяетпоследовательность невязок в методах, использующих подпространства Крылова:rm = pm (A)r0(B.1)Так, на использовании вместо pm (A) полинома p2m (A) основан метод CGS [13]; метод BCGSTAB вместо (B.1) использует соотношениеrm = pm (A)qm (A)r0 , где qm — специальным образом строящийся полином, такой что произведение pm qm не содержит нечетных степеней. Построение подобных методов, как правило, оказывается значительно более сложным.Из числа методов, свободных от транспонирования, в настоящеевремя широко применяется TFQMR2 , алгоритм которого показан нарис. B.1.12Англ.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
