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

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла

Методы решения СЛАУ большойразмерностиМ.Ю.БаландинЭ.П.Шурина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, а соответствующие элементы вектора правой части приравниваются к значениям,взятым из краевых условий.Подобная процедура нарушает симметричность матрицы.

Характеристики

Тип файла PDF

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

Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.

Список файлов книги