Баландин М.Ю., Шурина Э.П. - Методы решения СЛАУ большой размерности, страница 2

Избежать этого можно за счет вычеркивания из СЛАУ соответствующейстроки и столбца (поскольку значение решения в узле известно изкраевого условия, ненулевые элементы вычеркиваемого столбцаможно учесть в виде поправок к вектору правой части). Однако соответствующий учет изменяет размерность и портрет матрицы иреализуется достаточно сложно.Сохранить симметричный портрет матрицы можно более простым способом, если ввести в него некоторое количество нулевыхэлементов, дополняющих портрет до симметричного. Легко видеть, что все позиции таких элементов принадлежат портрету P Aисходной матрицы до учета краевых условий.Так, обнуление внедиагональных элементов третьей строки матрицы (1.2) приведет ее к несимметричной матрице Â, в которой позиции элементов, дополняющих портрет до симметричного, обведены:9 0 0 31 0 1 0 11 2 10 0 2  0 0 1 0000. = 2129100 1 0 0 1 12 0 1  0 0 0 00 8 0 2 2 0 03 0 8Для формата CSR обнуление внедиагональных элементов k-йстроки матрицы описывается следующим простым циклом:для i от iptr[k] до iptr[k+1]-1 {если jptr[i]=kтоaelem[jptr[k]]:=1иначе aelem[jptr[k]]:=0}Для случая раздельного хранения верхнего треугольника, нижнего треугольника и диагонали процедура несколько сложнее:9adiag[k]:=1для i от iptr[k] до iptr[k+1]-1 {altr[jptr[k]]:=0}для i от 1 до n {найти такое j (iptr[i]<j<iptr[i+1]) что jptr[j]=kесли найденотоautr[j]:=0}Очевидно, если элементы матрицы внутри строк упорядочены постолбцовым индексам (см.

замечание 2, §1.1), операция поиска может быть реализована более оптимально.Поскольку на практике известно множество D индексов тех узлов сетки, которые лежат на границе области с заданными условиями Дирихле, эту процедуру можно обобщить:для всех k из множества D {adiag[k]:=1для i от iptr[k] до iptr[k+1]-1 {altr[i]:=0если jptr[i]∈ Dтоautr[i]:=0}}1.4.Прямой и обратный ход по разреженнымтреугольным матрицамХранение матрицы СЛАУ в формате CSlR оказывается удобным длятех случаев, когда используется предобусловливание с разложениемна множители треугольной структуры, например ILU-факторизация 5(см.

§2.2).Поскольку после выполнения факторизации 6 треугольные матрицы L и U имеют те же портреты, что нижний и верхний треугольники матрицы A, для них достаточно хранить только самиэлементы матриц и возможно (в зависимости от модификации ILU)диагональные элементы; портреты будут полностью определятьсямассивами iptr и jptr.56Incomplete Lower-Upper triangles decompozitionО программной реализации ILU-декомпозиции см. также §2.4.10Пусть для матрицы A, заданной массивами adiag, altr, autr, jptr иiptr, найдена факторизация A ≈ LU, в которой L задается массивамиlltr и ldiag, а матрица U — массивами uutr и udiag.Тогда, если дан некоторый вектор f , прямой ход для системыLz = f может быть выполнен следующим образом:Входные данные: f; ldiag, lltr, jptr, iptr; nРезультат: z=L−1 fПобочные эффекты: изменяется массив fдля i от 1 до n {для j от iptr[i] до iptr[i+1]-1 {f[i]:=f[i]-z[jptr[j]]*lltr[j]}z[i]:=f[i]/ldiag[i]}Как видно, классическая процедура прямого хода реализуется без каких-либо существенных изменений, лишь с поправкой наспециальный способ хранения матрицы.Массив uutr, в отличие от lltr, хранит элементы матрицы U не встрочном, а в столбцовом порядке.

Поэтому обратный ход для системы Uz = f должен быть организован по-другому; соответствующаяпроцедура имеет видВходные данные: f; udiag, uutr, jptr, iptr; nРезультат: z=U−1 fПобочные эффекты: изменяется массив fдля i от n до 1 с шагом -1 {z[i]:=f[i]/udiag[i]для j от iptr[i] до iptr[i+1]-1 {f[jptr[j]]:=f[jptr[j]]-z[i]*uutr[j]}}11Глава 2Предобусловливание.Неполное LU-разложение2.1.ПредобусловливаниеПусть M — некоторая невырожденная матрица размерности n. Домножив (1) на M−1 , получим системуM−1Ax = M−1b ,(2.1)которая в силу невырожденности M имеет то же точное решение x ∗ .Введя обозначения  = M−1A и b̂ = M−1b, запишем (2.1) в видеÂx = b̂ .(2.2)Хотя (2.2) алгебраически эквивалентна (1), спектральные характеристики матрицы  отличаются от характеристик матрицы A,что, вообще говоря, ведет к изменению скорости сходимости численных методов для (2.2) по отношению к (1) в конечной арифметике.Процесс перехода от (1) к (2.2) с целью улучшения характеристикматрицы для ускорения сходимости к решению называется предобусловливанием1 , а матрица M — матрицей предобусловливателя.Из (2.1) сразу же вытекает важное требование: матрица M должна быть близка к матрице A.

(Выбор M = A сразу же приводит (1) квиду Ix = A−1b, однако не имеет практического смысла, так как требует нахождения A−1 , что, по существу и сводится к решению (1)).Вторым естественным требованием является требование легкой вычислимости матрицы M.Нахождение произведения M−1A для вычисления Â в общем случае ведет к изменению портрета (PA 6= PÂ ). Поэтому на практике используется другой подход.1Иногда переобусловливанием.12Невязка r̂ системы (2.2) связана с невязкой r системы (1) очевидным соотношениемMr̂ = r ,(2.3)которое справедливо и для матрично-векторных произведений z == Aq и ẑ = Âq:ẑ = Âq = M−1Aq =⇒ Mẑ = Aq = z .(2.4)Это позволяет вместо явного перехода от (1) к (2.2) вводить в схемы методов корректирующие шаги для учета влияния предобусловливающей матрицы2 .

Пример такого подхода будет приведен в §5.2.Из (2.4) следует еще одно условие: структура матрицы предобусловливателя должна допускать легкое и быстрое решение «обратных к предобусловливателю» систем вида Mẑ = z.Таким образом:• M должна быть по возможности близка к A;• M должна быть легко вычислима;• M должна быть легко обратима.Замечание. Выбор в качестве M некоторой диагональной матрицы удовлетворяет двум последним требованиям из трех перечисленных и, по существу, сводит предобусловливание к масштабированию СЛАУ. Как известно, в ряде случаев масштабирование способнозначительно ускорить процесс решения.Описанное выше предобусловливание иногда называется левым,так как домножение матрицы СЛАУ на матрицу предобусловливателя производится слева. Другой возможный подход основан на переходе от (1) к системеAM−1 y = b ,(2.5)у которой точное решение y∗ связано с точным решением x∗ исходнойСЛАУ соотношениемx∗ = M−1 y∗ .(2.6)Предобусловливание (2.5) реализуется путем выполнения вместоматрично-векторных умножений z := Aq двойных умножений z :== A(M−1 q); кроме того, при достижении требуемой точности осуществляется пересчет решения в соответствии с (2.6).Подобная схема предобусловливания называется правой.2В некоторых случаях, когда матрица M имеет простую структуру (например,диагональную — см.

ниже замечание), переход от A и b к Â и b̂ может выполнятьсяи явно. Подобное предобусловливание, чтобы подчеркнуть что оно выполняется внеметода, иногда называют внешним.132.2.ILU-факторизацияОдним из широко известных способов разложения матриц на множители является LU-факторизация [15], позволяющая представитьматрицу A в видеA = L A UA ,(2.7)где LA и UA — нижне- и верхнетреугольные матрицы соответственно.Подобное представление позволяет легко решать СЛАУ видаAx = b путем выполнения прямого хода для нижнетреугольнойсистемы LA y = b и затем обратного хода для верхнетреугольнойсистемы UA x = y. Однако алгоритм факторизации непригоден дляСЛАУ с разреженными матрицами, так как ведет к заполнению портрета, т.е. появлению в матрицах LA и UA ненулевых элементовв тех позициях, для которых aij = 0, — и как следствие, резкомуувеличению объема памяти, требуемой для хранения матриц.Сказанное можно проиллюстрировать на примере матрицы (1.2).Применение к ней как к плотной матрице алгоритма LU-факторизации дает (с точностью три знака) следующее разложение:LA =UA =10100.09110.222 0.091 0.18510.111000.08510000010.222 0.182 −0.037 −0.099 0.241 0 1901103219.818 1.9097.8891000.77811.82301020 −0.182 0 −0.37 .00.92 807.149;(2.8)(2.9)Сравнение (2.8) и (2.9) с (1.2) показывает, что коэффициенты l 73 ,l74 , u37 и u47 не равны нулю, тогда как соответствующие позиции матрицы A содержат нули.Вместо задачи нахождения факторизации (2.7) сформулируемдругую задачу.

Для заданной матрицы A потребуем представить еев видеA = LU + R ,(2.10)где матрицы в правой части удовлетворяют следующим свойствам:14• матрицы L и U являются нижнетреугольной и верхнетреугольной соответственно;• PL ⊂ PA и PU ⊂ PA ;• ∀(i, j) ∈ PA : [LU]ij = [A]ij ;• PA ∩ PR = ∅.Тогда приближенное представление A ≈ LU называется неполнойLU-факторизацией матрицы A или коротко ее ILU-разложением 3 .(Последние два требования, по существу, означают, что на множестве индексов PA матричное произведение LU должно точно воспроизводить элементы A.)Для нахождения матриц L и U будем генерировать их построчно.

Предположим, что первые (k − 1) строк уже найдены и необходимо найти k-ю. Запишем в блочном виде первые k строк разложения(2.10):A11 A12aT21 aT22!=L11 0lT21 1!U11 U120 uT22!R11 R12rT21 rT22+!, (2.11)где l21 , u22 , r21 и r22 — некоторые векторы. Выполнив действия надматрицами в правой части (2.11), получимA11 A12aT21 aT22!=L11 U11 + R11 L11 U12 + R12lT21 U11 + rT21 lT21 U12 + uT22 + rT22!.Из равенства матриц следует, что искомые векторы l 21 и u22должны удовлетворять условиям:lT21 U11 + rT21 = aT21 ;uT22+rT22=aT22−(2.12)lT21 U12.(2.13)Решив эти системы, можно найти коэффициенты k-х строк матрицразложения lk1 , .

. . , lk,k−1 , ukk , . . . , ukn 4 .Определим lkj из (2.12) в предположении, что lk1 . . . lk,j−1 уже найдены5 . Согласно ранее сформулированным условиям, если a kj = 0,3Чтобы подчеркнуть тот факт, что в данной постановке задачи в портрет не вводятся новые ненулевые элементы, такую факторизацию иногда называют факторизацией с нулевым заполнением.4Запись (2.11) подразумевает, что lkk ≡ 1.5Поскольку матрица L нижнетреугольна, для ее искомых коэффициентов всегдаj < k.15то сразу lkj = 0. В противном же случае rkj = 0 и (2.12) можно записать в виде6jXlki uij =j−1X(2.14)lki uij + lkj ujj = akj .i=1i=1Это позволяет вычислить lkj следующим образом7 :lkj =j−1X1 akj −lki uij  .ujji=1(2.15)Аналогичными рассуждениями с учетом l jj ≡ 1 из (2.13) можно получить выражение для ukj 8 :ukj = akj −k−1X(2.16)lki uiji=1(для тех случаев, когда akj 6= 0, иначе сразу ukj = 0).На основании формул (2.15) и (2.16) можно записать следующийалгоритм ILU-разложения:Для k = 1, nДля j = 1, k − 1Если (k, j) ∈/ PAто lkj := 0иначе lkj :=увеличить jlkk := 1Для j = k, nЕсли (k, j) ∈/ PAто ukj := 061 akj −lki uij ujji=1иначе ukj := akj −увеличить jувеличить kj−1Xk−1Pi=1lki uijС учетом того, что для i > j элементы верхнетреугольной матрицы U равнынулю.7Напомним, что строки матриц L и U с 1-й по k-ю предполагаются уже построенными.