Методы решения больших СЛАУ на разреженных матрицах общего вида

Авторы: Дяченко Т.Ф., Михайлова Т.В.
Источник: Тезизы доклада на I международной научно-технической конференции студентов, аспирантов и молодых ученых «Информационные управляющие системы и компьютерный мониторинг – 2010», 19-21 мая 2010г., ДонНТУ, Донецк

Аннотация

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

Общая постановка проблемы.

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

В научной литературе широко представлены методы решения СЛАУ для матриц особого вида или сводящихся к ним [1,2].

Постановка задач исследования.

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

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

Однако до разработки подобных параллельных алгоритмов следует решить две подзадачи:
1) Выбрать эффективный метод хранения ненулевых элементов
2) Выбрать аналитический метод решения СЛАУ.

Выбор метода решения СЛАУ

Методы решения большинства задач численных методов (в том числе задач решения линейных алгебраических систем) можно разделить на два типа: прямые и итерационные.

Метод решения задачи называется прямым, если он дает ее точное решение за конечное число действий.

Метод решения задачи называется итерационным, если он дает точное решение как предел последовательности приближений , вычисляемых по единообразной схеме, то есть где x - точное решение задачи, x(k) - k-тое приближение (итерация) к решению.

При реализации итерационного метода на ЭВМ его приходится обрывать на каком-то приближении с номером K и полагать x ≈ x^(k) . Обычно указывают погрешность (точность) ε, с которой требуется найти решение x и условие окончания итераций задают таким образом, чтобы норма разности точного решения и очередного итерационного приближения не превышала ||x - x(k)|| ≤ ε.

К прямым методам решения СЛАУ относятся известный из алгебры метод Крамера, метод Гаусса и его модификации и другие, к итерационным - метод простой итерации, метод Зейделя и другие.

Стандартные прямые методы, такие как поиск обратной матрицы (численно неустойчив) или метод Гаусса не учитывают разреженность матрицы А и мало подходят для программирования, особенно в нашем конкретном случае – при преобразованиях исходной матрицы А на каждом шаге возникает вычислительная погрешность, связанная с конечным числом разрядов ЭВМ для представления вещественных чисел.

Эффективность применения приближенных методов зависят от выбора начального вектора и быстроты сходимости процесса.

Рассмотрим метод простой итерации (метод последовательных приближений). Пусть дана система n линейных уравнений с n неизвестными:

Обозначим ; , где i = 1, 2, ...,n; j = 1,2,..., n. Тогда система (2) запишется таким образом в матричной форме

Решим систему (2) методом последовательных приближений. За нулевое приближение примем столбец свободных членов. Любое (k+1)-е приближение вычисляют по формуле

Если последовательность приближений x(0),...,x(k) имеет предел , то этот предел является решением системы (1), поскольку в силу свойства предела , т.е. [2].

Метод Зейделя [3] представляет собой модификацию метода последовательных приближений. В методе Зейделя при вычислении (k+1)-го приближения неизвестного xi (i>1) учитываются уже найденные ранее (k+1)-е приближения неизвестных xi-1.

Пусть дана линейная система, приведенная к нормальному виду:

(3)

Выбираем произвольно начальные приближения неизвестных и подставляем в первое уравнение системы (3). Полученное первое приближение подставляем во второе уравнение системы и так далее до последнего уравнения. Аналогично строим вторые, третьи и т.д. итерации.

Таким образом, предполагая, что k-е приближения известны, методом Зейделя строим (k+1)-е приближение по следующим формулам: , где k=0,1,...,n

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

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

Схемы хранения разреженных матриц

Основной целью является хранение только ненулевых элементов матрицей и возможность производить над ней простейшие матричные операции.

Пусть Nz – количество ненулевых элементов матрицы A. Рассмотрим наиболее распространенные схемы хранения, более подробная информация приводится в [4].

Самой простой схемой хранения разреженных матриц является так называемый «координатный» формат. Структура исходной матрицы отображается на три одномерных массива: (1) – массив, содержащий значения ненулевых элементов матрицы А в произвольном порядке; (2) целочисленный массив, содержащий их строковые индексы; и (3) целочисленный массив, содержащий их столбцовые индексы. Все три вектора имеют длину Nz.

Пример 1. Матрица

Может быть представлена как:

AA	12	9	7	5	1	2	11	3	6	4	8	10
JR	5	3	3	2	1	1	4	2	3	2	3	4
JC	5	5	3	4	1	4	4	1	1	2	4	3

В приведенном примере элементы перечислены в произвольном порядке, однако обычно они перечисляются по строкам или столбцам. Если элементы перечисляются по строкам, то вектор JC, который содержит избыточную информацию, может быть заменен массивом указателей на начало каждой строки. Это даст заметную экономию при хранении. Новая структура данных будет иметь три массива, которые будут выполнять следующие функции:

- массив АА содержит ненулевые значения aij, сохраняемые построчно, от строки 1 до n. Длина массива АА есть Nz.

- целочисленный массив JA содержит столбцовые индексы элементов aij в матрице А. Длина массива JА есть Nz.

- целочисленный массив содержит указатели входа для каждой строки в массивах АА и JА. Таким образом, значение IA(i) показывается позиция, где начинается i-тая строка в массивах АА и JА. Длина IA есть n+1, где IA(n+1) содержит число IA(1)+ Nz, то есть адрес в А и JА начала фиктивной строки с номером n+1.

Таким образом, представленная выше матрица А может быть сохранена как:

AA	1	2	3	4	5	6	7	8	9	10	11	12
JА	1	4	1	2	4	1	3	4	5	3	4	5
IA	1	3	6	10	12	13

Такой формат является наиболее популярным для хранения разреженных матриц общего вида. Он носит название разреженного строчного формата или Compressed Sparse Row (CSR). Эта схема более предпочтительна, чем координатная, так как часто является более удобной для нескольких важных операций над разреженными матрицами: сложения, умножения, перестановок строк и столбцов, транспонирования, решения ли-нейных систем с разреженными матрицами коэффициентов как прямыми, так и итерационными методами. С другой стороны, координатная схема более проста, наглядна и гибка, она часто используется в качестве «входного» формата для вычислительных библиотек, предназначенных для работы с разреженными матрицами.

Существует несколько модификаций данной схемы хранения, наиболее очевидной является разреженная столбцовая схема или Compressed Sparse Column (CSC).

Другая распространенная схема учитывает тот факт, что диагональные элементы большинства матриц ненулевые и/или используются чаще других. Модифицированная строчная схема (Modified Sparse Row, MSR) использует только два массива: массив значении АА и целочисленный массив JА. В первых n позициях АА содержит диагональные элементы исходной матрицы по порядку. Элемент АА(n+1) не заполняется (или же несет дополнительную информацию о матрице). Начиная с позиции n+2 записываются ненулевые элементы исходной матрицы по строкам, исключая диагональные. Для каждого элемента АА(k) элемент JA(k) показывает столбцовый индекс в исходной матрице. На n+1 позициях матрицы JA размещаются указатели входа для каждой строки матрицы АА и JA.

Поэтому для матрицы из примера выше эти два массива будут иметь вид:

AA	1	4	7	11	12	*	2	3	5	6	8	9	10
JА	7	8	10	13	14	14	4	1	4	1	4	5	3

Звездочка указывает на неиспользуемый элемент. Заметим, что JA(n)=JA(n+1)=14, показывая, что последняя строка является нулевой, так как диагональный элемент уже описан ранее.

Диагонально структурированные матрицы – это матрицы, чьи ненулевые элементы расположены вдоль небольшого числа диагоналей. Эти диагонали могут храниться в прямоугольном массиве DIAG(1:n, 1:Nd), где Nd–число диагоналей. Смещение каждой диагонали вычисляется относительно главной диагонали, которая должна быть известна. Смещения хранятся в массиве IOFF(1:Nd). Таким образом, элемент ai,i+ioff(j) исходной матрицы будет храниться в позиции (i,j) в массиве DIAG. Порядок, в каком хранятся диагонали, в общем случае, неважен, особенно если большинство операций сосредоточенно около главной диагонали, ее имеет смысл хранить в первом столбце. Также следует отметить, что все диагонали, кроме главной, имеют менее n элементов, так что не все элементы матрицы DIAG будут использоваться.

Пример 2. Трехдиагональная матрица

Более общей схемой, популярной для векторных машин является так называемый формат. Эта схема предполагает, что количество диагоналей Nd невелико. Для хранения требуются два прямоугольных массива размером n*Nd каждый.

В первом, COEF, подобном DIAG, хранятся ненулевые элементы матрицы А. Ненулевые элементы каждой строки могут сохраняться в строке массива COEF(1:n,1:Nd), дополняя эту строку нулями, если необходимо. Вместе с COEF, целочисленный массив JCOEF(1:n,1:Nd) должен хранить индекс столбца каждого элемента COEF.

Для матрицы А из примера 2, схема хранения Ellpack-Itpack примет вид:
        COEF =

1	2	0
3	4	5
6	7	8
9	10	0
11	12	0

                      JCOEF =

1	3	1
1	2	4
2	3	5
3	4	4
4	5	5

Выводы.

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

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

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

1. Тьюарсон Р. Разреженные матрицы. М.: Мир, 1977.
2. Джордж А., Лю Д. Численные методы решения больших разреженных систем уравнений. М.: Мир, 1984.
3. Фельдман Л.П., Петренко А.И., Дмитриева О.А. Численные методы в информатике. BHV, 2005.
4. Писсанецки С. Технология разреженных матриц. М.: Мир, 1988.
5. Saad Y. Iterative Methods for Sparse Linear Systems, 2000.