Формат Skyline для хранения и обработки разреженных матриц нерегулярной структуры для численного решения СЛАУ Текст научной статьи по специальности «Математика»

84

УДК 519.612

С. В. Старовойтов

ФОРМАТ SKYLINE ХРАНЕНИЯ И ОБРАБОТКИ РАЗРЕЖЕННЫХ МАТРИЦ НЕРЕГУЛЯРНОЙ СТРУКТУРЫ ДЛЯ ЧИСЛЕННОГО РЕШЕНИЯ СЛАУ

Описан формат Skyline для хранения и обработки разреженных матриц нерегулярной структуры. Определен алгоритм ILUS-пре-добусловливания разреженных матриц в формате Skyline.

Skyline format for sparse linear systems and ILUS (an incomplete LU preconditioner for matrixes in sparse Skyline format) are examined.

Ключевые слова: формат Skyline, разреженная матрица, нерегулярная структура, алгоритм.

Введение

Современные требования к вычислительной математике и математическим моделям, основанным на решении сложных систем дифференциальных уравнений в частных производных на произвольных геометрических областях, привели к необходимости решения больших разреженных систем линейных алгебраических уравнений с матрицами нерегулярной структуры. Цель статьи — описание программной реализации решения СЛАУ с разреженной матрицей нерегулярной структуры большой размерности с условием хранения матрицы в формате Skyline, лучшим из мною исследованных форматов хранения разреженных матриц, так как он идеально подходит не только для выполнения всех стандартных операций над матрицами, но и для алгоритма ILUS-предобу-словливания, описанного в статье. Пример программной реализации (на Borland Pascal 7.0 for Windows) всех алгоритмов, приведенных в статье, доступен в Интернете: www.nwtorg.ru/skyline.pas.

1. Формат Skyline хранения и обработки разреженных матриц с симметричным портретом

l.l. Структура формата Skyline

Портретом разреженной матрицы A называется множество пар индексов (i, j) таких, что aj Ф 0, т.е. PA = {(i, j) | aj Ф 0}. Возможны три случая:

• матрицы с несимметричным портретом, в этом случае для них 3(i, j) є PA : (j, i) g PA (такие матрицы несимметричны A Ф A );

• несимметричные матрицы с симметричным портретом, для которых выполняется (в общем случае aj Ф a^):

(i, і)є pa » (j, i) є pa ; (1.1)

• симметричные матрицы, у которых a j = aц ((1.1) для них выполняется).

Вестник РГУ им. И. Канта. 2008. Вып. 10. Физико-математические науки. С. 84 — 94.

W

Рассмотрим случай несимметричных матриц с симметричным портретом и более строгий случай симметричных матриц при реализации метода CG — сопряженных градиентов. Перейдем к описанию структуры формата Skyline. Пусть матрица системы A является разреженной матрицей с симметричным портретом размерности N, тогда к ней может быть применена схема хранения Skyline. Матрица системы хранится в виде 5 динамических одномерных массивов:

1. adiag[] содержит N элементов главной диагонали;

2. altr[] содержит поддиагональные элементы (т. е. элементы нижнего левого треугольника) в строчном порядке;

3. autr[] содержит наддиагональные элементы (т. е. элементы правого верхнего треугольника), причем он хранит элементы не в строчном, а в столбцовом порядке (в силу симметрии портрета количество элементов в нем совпадает с altr[]);

4. jptr[] содержит столько же элементов, сколько altr[] (autr[]), и для каждого из них указывает столбец (строку);

5. iptr[] содержит N +1 элементов, и его i-й элемент указывает, с какой позиции в массивах altr[] и jptr[] начинается i-ая строка матрицы. Последний элемент iprt[ N +1 ] равен числу элементов в altr[] плюс 1 (в силу симметрии портрета для autr[] массива iptr[] не формируется).

Рассмотрим следующую разреженную матрицу:

A=

' 7 1 О 3 О О 14

2 1О О 2 1 О 2

О О 8 О О О О

1 3 О 12 1 О О

О 1 О 1 9 О 3

О О О О О 11 О

V 2 1 О О 2 О 9,

(1.2)

Для матрицы (1.2) указанные массивы должны выглядеть так: аШ = [2,1, 3,1,1, 2,1, 2], аыЬг = [1, 3, 2,1,1,1, 2, 3], adiag= [7,10,8,12,9,11,9],

]рЬг = [1,1,2,2,4,1, 2,5], iptr = [1,1, 2, 2, 4,6, 6, 9].

Разберем массив 1рЬ"[], рассмотрим поддиагональные элементы:

85

1 2 1

0 2

1 2 О

4

О

б

21

V б

О

1 О 1 •.

ОООО ОО2

iptr[1] = 1 iptr[2] = 1 iptr[3] = 2

iptr[4] = 2 Последний элемент iptr[N +1] = 9 iptr[5] = 4 iptr[6] = 6 iptr [7] = 6

О

Рис. 1. Схема построения массива iptr[]

Ненулевой элемент появляется во 2-й строке матрицы, поэтому iptrp] = 2. В 4-й строке есть два ненулевых элемента: 1, З и соответственно iptr[5] = 2 + 2 = 4. Таким образом, iptr[k] = iptr[k-1] + {кол-во элементов из PA в строке k-1}, причем всегда iptr[1]=iptr[2]=1 (конечно, хранить 2 элемента в памяти компьютера необязательно, но на скорость вычислений это не влияет, и они остаются в соответствии с форматом).

В случае, когда матрица симметрична, массивы altr[] и autr[] для нее совпадают, поэтому достаточно хранить только один из них.

1.2. Матрично-векторное умножение

Умножение разреженной матрицы на вектор z := Ax в случае использования Skyline формата реализуется так: поскольку нижний треугольник хранится по строкам, производится последовательный перебор элементов массива altr, которые умножаются на коэффициент вектора x с индексом, взятым из соответствующей позиции массива jptr; результат добавляется к коэффициенту вектора z, соответствующему текущей сканируемой строке. Для верхнего треугольника строчные и столбцовые индексы меняются местами, а элементы главной диагонали учитываются отдельно:

Листинг 1

Схема выполнения матрично-векторного умножения

Входные данные: 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, а какого-нибудь одного из них.

1.3. Прямой и обратный ход по матрицам в Skyline формате

Хранение матрицы СЛАУ в Skyline формате удобно для тех случаев, когда используется предобусловливание с разложением на множители треугольной структуры, например ILU-факторизация (например, ILUS: an incomplete LU preconditioner in sparse Skyline format [1]).

Поскольку после выполнения факторизации треугольные матрицы L и U имеют те же портреты, что нижний и верхний треугольники матрицы A, для них достаточно хранить только сами элементы матриц

и, возможно (в зависимости от модификации 1ЬИ), диагональные элементы; портреты будут полностью определяться массивами ірії и ]'рЬг.

Пусть для матрицы А, заданной массивами adiag, а1іг, au.tr, ]'р1г и іріг, найдена факторизация А и Ьи, где Ь задается массивами 11іг и 1diag, матрица и — массивами uu.tr и udiag. Тогда, если дан некоторый вектор /, имеем прямой ход для системы Ьг = / :

Листинг 2

Схема выполнения прямого хода для системы Lz=f

Входные данные: ^ ldiag, lltr, jptr, iptr,• п Результат: z=L-1f

Побочные эффекты: изменяется массив f для ± от 1 до п {

для j от iptr[i] до iptr[i+1]-1 { f [±] ^[±]^^гШ]*ЦЛ;гШ

}

z[i]:=f[i]/ldiag[i]

}

S7

Итак, классическая процедура прямого хода реализуется без существенных изменений, лишь с поправкой на специальный способ хранения матрицы. Массив uutr, в отличие от Шг, хранит элементы матрицы и не в строчном, а в столбцовом порядке. Поэтому обратный ход для системы иг = / должен быть организован по-другому; соответствующая процедура имеет вид:

Листинг 3

Схема выполнения обратного хода системы Uz=f

Входные данные: ^ ud±ag, uutr, jptr, iptr,• п Результат: 2 = и 1/

Побочные эффекты: изменяется массив f для ± от п до 1 с шагом -1 { z [±]:=f[±]/udiag[i] для j от iptr[i] до iptr[i+1]-1 {

f[jptr[j]]:=f[jptr[j]]-z[±]*uutr[j]

}

}

2. Решение СЛАУ методом CG с ILUS-предобусловливанием

2.1.1ШБ-факторизация с сохранением портрета.

Сформулируем задачу 1ЬиБ-факторизации. Для заданной матрицы А потребуем представить ее в виде:

А = Ш + Я, (2.1)

где матрицы в правой части удовлетворяют следующим свойствам:

88

• матрицы Ь и и являются нижнетреугольной и верхнетреугольной соответственно;

• Рь с Ра и ри с Ра ;

• V(i',;) є Ра : [Ьи] = [А];

• Ра ° Рк = 0-

Тогда приближенное представление А и Ьи называется неполной Ьи -факторизацией матрицы А, или коротко, ее ІЬиБ-разложением. (Последние два требования, по существу, означают, что на множестве индексов Ра матричное произведение Ьи должно точно воспроизводить элементы А.)

Для нахождения матриц Ь и и будем генерировать их построчно. Предположим, что первые (к -1) строк уже найдены и необходимо найти к -ю. Запишем в блочном виде первые к строк разложения (2.1):

А“ А21=\^ 0їи 11 иі21+\Кт21, (2.2)

а21 а22 ) V 121 1)1 0 и22 ) V Г21 Г22 )

где І21, и22, ?21 и г22 — некоторые векторы. Выполнив действия над матрицами в правой части (2.2), получим

Ац А12 | \ Ьциц + Ьци12 + -^12

0^21 0-22 ) і Іц^^11 + Г*21 І21^^12 + и22 + Г22

Матрицы равны, поэтому искомые векторы І21 и и22 удовлетворяют

І21и 11 + Г21 = °21, (2.3)

и12 + г222 = Й22 - І21и12 . (2.4)

Решив эти системы, найдем коэффициенты к-х строк матриц разложения Ік1,..., Ік к-1, икк,..., икп . Определим ІЩ из (2.3) в предположении, что Ік 1,...,Ік,]-1 уже найдены. Согласно нашим условиям, если а. = 0, то 1к. = 0. В противном случае г. = 0, и (2.3) запишем в виде

І І-1

X Ікіиіі = X Ікіиіі + ІкІиІІ = акІ . (2.5)

і=1 і=1

Это позволяет вычислить іці следующим образом:

\ і- 1 1

І =-1-

Ікі =

І иІІ

/

акІ - X Ікіи

V і=1

Аналогично с учетом і.. = 1 из (2.4) получаем выражение для и^ :

Ікі = а1д

(2.6)

ІІ к-1

икі = акі - X Ікіиі/' (2.7)

1=1

(для тех случаев, когда ац Ф 0, иначе сразу и^ = 0).

На основе формул (2.6) и (2.7) можно записать следующий алгоритм ТЬИБ-разложения:

Листинг 4

Общая схема алгоритма ГШБ-факторизации

Для k = 1, n

Для ] = 1, к -1

Если (к, ]) £ РА то 1к] := 0

иначе 1ц :=

увеличить ]

1кк := 1 Для ] = к, и

Если (к, ]) £ РА

то ик := 0

;_1

akj _ X lkiuij i=1

s9

k_1

иначе и

= akj _ X lkiuij

i=1

увеличить ] увеличить к

Пример.

Для матрицы (1.2) приведенный алгоритм дает следующие матрицы Ь и и (с точностью три знака):

L=

Г 1 О,286 О

О,143

О

О

О,286

Г 7

U=

1

О

О,294

О,1О3

О

0,074

О

1

О,О79 О О

1

О

9,714 О 1,143

8

О

О, 218 0 -1J

О О 1 1

1 О 1,714

О О 0

О,7О6 О 0

8,842 О 2,824

11 0

7,973

Соответствующая им матрица ошибки разложения И равна

R=A_LU=

Г о О О О О О О

0 О О О О О О

0 О О О О О О

0 О О О О О _ О,647

0 О О О О О О

0 О О О О О О

10 О О _ О,943 О О О

1

u

1

9д

Очевидно, что матрица LU удовлетворяет всем трем требованиям к матрице предобусловливателя: она приближает матрицу A, так как на множестве индексов Pa точно воспроизводит ее; она легко вычисляется по приведенному выше несложному алгоритму; наконец, она легко обратима, так как является произведением двух треугольных матриц.

Таким образом, описанный выше алгоритм является достаточно хорошим способом предобусловливания и идеально подходит для Skyline-формата хранения разреженных матриц, так как порядок доступа совпадает с порядком хранения для обеих матриц. Полученная схема нахождения ILUS-факторизации выглядит так:

Листинг Б

Схема алгоритма ILUS-факторизации для матриц в формате Skyline

для k от 1 до n {

для j от iptr[k] до iptr[k+1]-1 {

вычислить ltr[j] по формуле для lk, iptr[ j]

}

вычислить ludiag[k] по формуле для Ukk для j от k+1 до n{

найти такое q (iptr[j]<=q<iptr[j+1] что jptr[q]=k если найдено

то вычислить utr[q] по формуле для Uk j

}

}

2.2. Метод CG с ILUS-предобусловливанием

До сих пор мы требовали от матрицы системы лишь симметричность портрета, однако на этот раз матрица системы должна быть симметричной. Мы не будем описывать вывод метода CG (его можно найти в любом учебнике по вычислительной математике), а приведу лишь общую схему метода, которая также подходит для хранения матрицы системы в Skyline-формате.

Пусть A — матрица (1), B — матрица предобусловливателя: B и A и B-1 легко вычисляется. Б нашем случае B = LU, где L и U — матрицы, полученные при ILUS-факторизации, описанной ранее. Тогда имеем схему метода CG:

Листинг 6

Схема метода CG с предобусловливанием

Выбрать начальное приближение Uq

Начальная невязка: Tq := b _ AUo , Mq := B_1ro Вычисляем параметр: Po := (mq,Jo)

Начинаем цикл: m := 1, N

am :_ A®m_1

Pm_1

(am r ®m_l)

Im := Um_1 +TmMm_1

m m_1

qm := B rm

pm := (qm, rm ) Pm

®ш • + ®т-1

рт-1

увеличить т

Это прямой метод с теоретической точки зрения, но как прямой он не применяется. Метод сопряженных градиентов с предобусловлива-нием иногда также называют РСС.

9l

3. Результаты

3.1. Общие сведения

Тестирование скорости работы отлаженной программы проводилось на алгоритмах скорейшего спуска и метода СЄ (признанного одним из самых быстрых и устойчивых медов решения СЛАУ). При практических экспериментах рассматривалась модельная задача (задача Пуассона на прямоугольной области) и задача Пуассона для области, представляющей прямоугольный равнобедренный треугольник.

Напомним, задача Пуассона состоит в нахождении такой функции и : Яп з В ^ Я, чтобы V* є В с Яп : УЫ = /(х).

Такие области были выбраны не случайно. При построении матрицы системы для прямоугольной области мы получим регулярную матрицу — блочно-трехдиагональную с размером и количеством блоков соответствующих количеству узлов на горизонтали и вертикали. Матрица прямоугольной области 3 х 3 представлена на рисунке 2.

( 4 1 О

1 4 1

0 1 4

1 О О

0 1 О

0 О 1

0 О О

О О О

10 О О

1 О О

О 1 О

О О 1

4 1 О

1 4 1

О 1 4

1 О О

О 1 О

О О 1

О О 01

О О 0

О О 0

1 О 0

О 1 0

О О 1

4 1 0

1 4 1

О 1 4 J

Рис. 2. Блочно-трехдиагональная матрица СЛАУ для задачи Пуассона на прямоугольной области 3 х 3

Для СЛАУ с такими матрицами давно известны результаты, и они учитывались при отладке программы.

В случае же треугольной области ситуация сложнее. Попробуем построить матрицу для небольшой области такого типа. Занумеруем узлы:

92

1

2 3 4 5 6 7 8 9 10

Рис. 3. Последовательность обхода по узлам

Двигаясь по узлам сетки в соответствии с последовательностью обхода (рис. 3), последовательно получаем симметричную разреженную матрицу системы. Эта матрица имеет относительно нерегулярную структуру. Для решения СЛАУ с матрицами подобной (нерегулярной) структуры написаны данная статья и программа:

Г 4 1 О О О О О О О

1 4 1 1 О О О О О

О 1 4 О 1 О О О О

О 1 О 4 1 О 1 О О

О О 1 1 4 1 О 1 О

О О О О 1 4 О О 1

О О О 1 О О 4 1 О

О О О О 1 О 1 4 1

О О О О О 1 О 1 4

V

Рис. 4. Матрица СЛАУ для 9 узлов треугольной области 3.2. Численные эксперименты

Тестирование программы проводилось на 32-разрядном компьютере с использованием компилятора Borland Pascal 7.0 for Windows.

Для оценки скорости вычислений сравним скорости работы программы для СЛАУ, получающейся при решении задачи Пуассона на прямоугольной и на треугольной области методом CG. Также рассмотрим решение задачи методом скорейшего спуска. Условием остановки методов будет максимальное отклонение от точного значения вектора решения < 10-10.

Таблица 1

Тестирование метода CG на прямоугольной области

Размер матрицы Бремя предобусловливания Число итераций Общее время, с

4S4 0,00 39 0,06

1600 0,05 70 0,44

3249 0,33 95 1,37

5041 0,77 126 3,02

Формат Skyline для хранения и обработки разреженных матриц у.—^

= ^

Таблица 2

Тестирование метода CG на треугольной области

Размер матрицы Бремя предобусловливания Число итераций Общее время, с

49б 0,00 49 0,0б

159б 0,05 S7 0,55

3240 0,33 134 1,81

5050 0,77 177 3,90

Для сравнения приведем результатні вычислений методом скорейшего спуска (матрицу СЛАУ будем также хранить в формате Skyline и предобусловливание будем выполнять по тому же алгоритму).

Таблица З

Тестирование метода скорейшего спуска на прямоугольной области

93

Размер матрицы Бремя предобусловливания Число итераций Общее время, с

4S4 0,00 25б 0,38

1б00 0,11 S24 4,34

3249 0,33 1б51 17,бЗ

5041 0,77 2545 42,35

Таблица 4

Тестирование метода скорейшего спуска на треугольной области

Размер матрицы Бремя предобусловливания Число итераций Общее время, с

49б 0,00 214 0,33

159б 0,05 б59 3,41

3240 0,33 1317 14,01

5050 0,77 2043 34,00

Из результатов видно, что СЛАУ методом CG решается на порядок быстрее, чем методом скорейшего спуска. Кроме того, число итераций растет как O(n) в отличие от метода спуска, у которого это число растет

как O(n2) (размер матрицы СЛАУ n х n).

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

Заключение

В статье показано практическое применение всего потенциала итерационных методов для решения СЛАУ и построен эффективный

94

«программный решатель» для систем с разреженными матрицами нерегулярной структуры. В случае более узкого класса задач по решению больших разреженных СЛАУ возможны улучшения, но они будут лишь этапными: более удобная схема хранения, более эффективное предобусловливание, более подходящий итерационный метод, кроме того, каждый этап, в случае решения на многопроцессорной ЭВМ, может быть распараллелен — общий же подход останется без изменений.

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

Данная работа выполнена при поддержке Российского фонда фундаментальных исследований (РФФИ) по проекту №08-01-00431.

Список литературы

1. Chow E., Saad Y. ILUS: an incomplete LU preconditioner in sparse Skyline format // Int. J. for Num. Meth. in Fluids. 1997. Vol. 25. P. 739 - 748.

2. Saad Y. Iterative methods for sparse linear systems. PWS Publishing Company, 1996.

3. Писсанецки С. Технология разреженных матриц. М., 1988.

4. Баландин М. Ю., Шурина Э. П. Методы решения СЛАУ большой размерности: Учеб. пособие. Новосибирск, 2000.

Об авторе

Старовойтов С. В. — асп., РГУ им. И. Канта.