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

Методы решения СЛАУ большойразмерностиМ.Ю.БаландинЭ.П.Шурина1 августа 2000 г.АннотацияДокумент является электронной версией одноименногоучебного пособия, выпущенного в июне 2000 г. в издательствеНовосибирского Государственного Технического Университета(НГТУ). Исправлены замеченные в печатном издании опечаткии неточности.Документ может свободно распространяться при условииотсутствия внесения изменений и неизвлечения коммерческойвыгоды. При цитировании необходимо ссылаться на следующее издание:Баландин М.Ю., Шурина Э.П.Методы решения СЛАУ большой размерности. —Новосибирск: Изд-во НГТУ, 2000. — 70 стр.Связаться с авторами можно по следующим адресам электронной почты:dolphin@gts.nsk.su(М.Ю.Баландин)shurina@online.sinor.ru (Э.П.Шурина)Пособие подготовлено в системе LATEX1ВведениеРазвитие вычислительной техники и вызванный этим процессом переход к более сложным (трехмерным, в произвольных геометрических областях) моделям в виде систем дифференциальных уравнений в частных производных и их дискретным аналогам на неструктурированных сетках, привел к необходимости решения большихразреженных систем линейных алгебраических уравнений с матрицами нерегулярной структуры.Наиболее эффективными и устойчивыми среди итерационныхметодов решения таких систем уравнений являются так называемые проекционные методы, и особенно тот их класс, который связан с проектированием на подпространства Крылова.

Эти методыобладают целым рядом достоинств: они устойчивы, допускают эффективное распараллеливание, работу с различными строчными(столбцовыми) форматами и предобусловливателями разных типов.В данном учебном пособии рассматриваются форматы храненияразреженных матриц и реализация операций над ними (в частности, матрично-векторное умножение), предобусловливание (неполное LU-разложение), классические итерационные и проекционныеметоды на подпространствах Крылова. Приведены основные теоретические обоснования, рассмотрены системы с симметричнымии несимметричными матрицами, приведены конкретные примеры. Для желающих использовать при реализации методов готовоепрограммное обеспечение в приложении даны соответствующиеInternet-адреса.2Используемые обозначения исоглашенияОсновной рассматриваемой задачей будет являться решение СЛАУ(1)Ax = bс квадратной невырожденной матрицей A и ненулевым вектором bразмерности n, составленными из действительных коэффициентов.В силу невырожденности матрицы у системы (1) существуетнекоторое точное решение x∗ = A−1 b.Невязкой системы (1), соответствующей вектору x, будем называть вектор r = b − Ax = A (x∗ − x).

Здесь вектор ϑ = x∗ − x естьошибка решения.Линейный оператор A, заданный матрицей A, действует в пространстве Rn . Введем в нем скалярное произведение(x, y) =nX(2)xj yj .j=1Для действительного случая сопряженный оператор A ∗ определяется матрицей AT , так что (Ax, y) = (x, ATy). Иногда также оказывается удобным использование скалярного A-произведения(x, y)A = (Ax, y) = (x, ATy) = y TAx = xTATy .(3)На основе скалярных произведений (2) и (3) введем в R n векторные нормыkxk2 =kxkA =q2q2(x, x) ;(Ax, x) =q2(x, x)A(евклидова норма и A-норма соответственно).Если явно не оговорено обратное, то будет предполагаться, чтоматрица A не обладает свойством симметрии и положительной определенности (A-норма, очевидно, применима лишь к таким матрицам1 , иначе она вообще не является нормой).1Часто называемым SPD-матрицами (Symmetric and Positively Defined).3Для невырожденной A симметризованная матрица ATA является положительно определенной, и тогда для евклидовой и A-нормысправедливо соотношениеkx − x∗ k2ATA = (ATA(x − x∗ ), (x − x∗ )) == (A(x − x∗ ), A(x − x∗ )) = (rx , rx ) = krx k22 .(4)Собственные значения матрицы A (которые, вообще говоря, являются комплексными) будем обозначать λ j (A).

Максимальный извсех модулей |λj (A)| есть спектральный радиус матрицы.Вектор ek будет обозначать k-й единичный орт:ek = (0, . . . , 0, 1k , 0, . . . , 0)T .4Глава 1Хранение и обработкаразреженных матриц1.1. Форматы храненияОпределение 1.1.1 Портретом разреженной матрицы A называется множество пар индексов (i, j), таких что a ij 6= 0:defPA = {(i, j) | aij 6= 0} .Можно выделить следующие три случая:1. матрицы с несимметричным портретом, т.е. матрицы, у которых∃(i, j) ∈ PA : (j, i) ∈/ PA .Такие матрицы несимметричны и в обычном смысле A 6= AT ;2.

несимметричные матрицы с симметричным портретом, у которых(1.1)(i, j) ∈ PA ⇐⇒ (j, i) ∈ PA ,хотя в общем случае aij 6= aji ;3. симметричные матрицы, у которых a ij = aji (очевидно, (1.1)для них выполняется).Разреженную матрицу можно представить графом со множеством вершин {1, 2, . . . , n}, тогда PA будет являться множеством ребер этого графа. В случае симметричного портрета такой граф будетнеориентированным, в случае несимметричного — ориентированным.5Очевидно, что для симметричных матриц достаточно хранитьтолько один из треугольников (верхний или нижний).

Для несимметричных матриц с симметричным портретом необходимо хранение всех ненулевых элементов, однако схему их размещения опятьже достаточно задать лишь для одного из треугольников. Следуеттакже принимать во внимание тот факт, что генерируемые в МКО иМКЭ глобальные матрицы СЛАУ всегда имеют ненулевую главнуюдиагональ.Распространенным способом хранения несимметричных матрицпроизвольной структуры является CSR 1 [6].

В нем разреженная матрица A хранится с использованием следующих массивов:• aelem, который содержит все ненулевые элементы матрицы A,перечисленные в строчном порядке;• jptr, который содержит столько же элементов, сколько aelem идля каждого из них указывает, в каком столбце находится данный элемент;• iptr, который хранит число элементов, равное увеличеннойна единицу размерности СЛАУ.

Его i-й элемент указывает, скакой позиции в массивах aelem и jptr начинается i-я строкаматрицы. Соответственно iptr[i+1]-iptr[i] равно числу ненулевых элементов в i-й строке. Последний элемент iptr[n+1] равенчислу элементов в aelem, увеличенному на единицу. 2П РИМЕР. Проиллюстрируем сказанное на следующей разреженной матрице:9 0 0 3 1 0 1 0 11 2 1 0 0 2  0 1 10 2 0 0 0 A=2 9 1 0 0 (1.2) 2 1. 1 00112010 0 0 8 0  0 02 2 0 0 3 0 8Для нее n = 7 и массивы имеют вид:aelem = [9, 3, 1, 1 , 11, 2, 1, 2 , 1, 10, 2 , 2, 1, 2, 9, 1 , 1, 1, 12, 1 , |{z}8 , 2, 2, 3, 8 ]| {z } |4(1){z4(5)} | {z } |3(9){z5(12)} |{z4(17)}1(21)| {z }4(22)jptr = [1, 4, 5, 7 , 2, 3, 4, 7 , 2, 3, 4, 1, 2, 3, 4, 5 , 1, 4, 5, 7 , |{z}6 , 1, 2, 5, 7 ]| {z } | {z } | {z } |4(1)4(5)3(9)iptr = [1, 5, 9, 12, 17, 21, 22, 26]1{z5(12)} | {z }4(17)1(21)| {z }4(22)Compressed Sparse Row.Поскольку всегда iptr[1]=1, его можно и не хранить.

В этом случае iptr[i] указывает на начало не i-й, а (i + 1)-й строки.26(фигурными скобками сгруппированы элементы, принадлежащиеодной и той же строке матрицы. Под скобками указано число сгруппированных элементов и порядковый номер первого из них в массиве).Замечание 1. Возможно использование модификации CSR, в которой данные хранятся тем же образом и в тех же массивах, однако ненулевые элементы матрицы перечисляются не по строкам, а постолбцам. Данный вариант известен как CSC. 3Замечание 2.

Элементы, принадлежащие одной строке, могутхраниться как с упорядочиванием по столбцовому индексу, так и безупорядочивания. Это не влияет на алгоритм матрично-векторногоумножения из §1.2, однако как будет показано уже в §1.3, в определенных условиях упорядочивание способно дать выигрыш.Для хранения матриц с симметричным портретом и ненулевойглавной диагональю может быть применена та же CSR-схема, со следующими изменениями4 :• элементы главной диагонали хранятся отдельно в массивеadiag;• вместо массива aelem используется массив altr, в котором таким же образом хранятся только поддиагональные элементыматрицы. Для него формируются массивы jptr и iptr;• для наддиагональных элементов матрицы формируется массивautr, причем он хранит элементы не в строчном, а в столбцовомпорядке.

Для autr в силу симметрии портрета массивы jptr и iptrне формируются.Для матрицы (1.2) указанные массивы должны выглядеть следующим образом:adiag = [9, 11, 10, 9, 12, 8, 8]altr = [1, 2, 1, 2 , 1, 1 , 2, 2, 3]| {z } |{z} | {z }autr = [2, 3, 1, 2 , 1, 1 , 1, 2, 1]| {z } |{z} | {z }jptr = [2, 1, 2, 3 , 1, 4 , 1, 2, 5]| {z } |{z} | {z }iptr = [1, 1, 1, 2, 5, 7, 7, 10](Как и прежде, фигурными скобками сгруппированы элементы одной строки.)3Compressed Sparse Column.Такая схема иногда обозначается CSlR — Compressed Sparse (lower triangle)Row.

Встречается также термин «Skyline Format» [12].47В случае, когда матрица симметрична, массивы altr и autr для неесовпадают, поэтому достаточно хранить только один из них.1.2.Матрично-векторное умножениеУмножение разреженной матрицы на вектор z := Ax наиболее просто реализуется для формата CSR. Поскольку матрица хранится построкам, производится последовательный перебор элементов массива aelem, которые умножаются на коэффициент вектора x с индексом, взятым из соответствующей позиции массива jptr; результатдобавляется к коэффициенту вектора z, соответствующему текущейсканируемой строке:Входные данные: x; aelem, jptr, iptr; nРезультат: z=Axдля i от 1 до n {z[i]:=0для j от iptr[i] до iptr[i+1]-1 {z[i]:=z[i]+x[jptr[j]]*aelem[j]}}При использовании раздельного хранения главной диагонали итреугольников схема умножения для нижнего треугольника остается той же.

Для верхнего треугольника строчные и столбцовые индексы меняются местами, а элементы главной диагонали учитываютсяотдельно:Входные данные: x; adiag, altr, autr, jptr, iptr; nРезультат: z=Axдля i от 1 до n { z[i]:=x[i]*adiag[i] }для i от 1 до n {для j от iptr[i] до iptr[i+1]-1 {z[i]:=z[i]+x[jptr[j]]*altr[j]z[jptr[j]]:=z[jptr[j]]+x[i]*autr[j]}}Если поменять местами обращения к массивам altr и autr, то вместо z := Ax будет вычислено произведение z := ATx.Этот же алгоритм может быть применен и для случая симметричных матриц. Единственное необходимое изменение заключается виспользовании не двух массивов altr и autr, а какого-нибудь одногоиз них.81.3.Симметричность портрета и учеткраевых условийИзвестно, что при МКО и МКЭ-моделировании физических процессов краевые условия первого рода (условия Дирихле) учитываютсяпутем обнуления внедиагональных элементов тех строк глобальной матрицы, которые соответствуют узлам сетки на границе, приэтом диагональные элементы приравниваются к 1, а соответствующие элементы вектора правой части приравниваются к значениям,взятым из краевых условий.Подобная процедура нарушает симметричность матрицы.