Решение разреженных систем линейных уравнений

Решение разреженных систем линейных уравнений

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

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

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

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

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

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

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

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

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

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

В приведенном примере элементы перечислены в произвольном порядке, однако обычно они перечисляются по строкам или столбцам. Если элементы перечисляются по строкам, то вектор 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.

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

Такой формат является наиболее популярным для хранения разреженных матриц общего вида. Он носит название разреженного строчного формата или 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.

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

Звездочка указывает на неиспользуемый элемент. Заметим, что 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. Трехдиагональная матрица

может быть сохранена в двух массивах согласно вышеприведенной схеме:
DIAG =

IOFF =

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

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

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

JCOEF =

Выводы.

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

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

Научная статья на тему «Решение плохо обусловленных разреженных систем линейных алгебраических уравнений с помощью крыловского подпространства»

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

Когда большая доля коэффициентов матрицы состоит из нулей, вполне очевидно, что нам не хотелось бы хранить в памяти компьютера все эти нули. Поэтому матричные алгоритмы должны проектироваться таким образом, чтобы обрабатывались только ненулевые элементы и чтобы на основании предварительного знания о расположении ненулевых элементов избегались операции типа сложения с нулем или умножения на нуль. Таким образом число операций, производимых машиной при исполнении алгоритма, пропорционально числу ненулевых элементов, а не числу всех элементов матрицы. Серьезную проблему при работе с разреженными матрицами представляет численная устойчивость [1, с. 195].

Когда методы типа гауссова исключения требуют слишком много времени или памяти для решения систем уравнений используются итерационные методы. При решении плохо обусловленных разрежен­ных СЛАУ возникает необходимость выбора метода, который позволит получить при решении точный результат и наименьшее заполнение (возникновение новых ненулевых элементов) [3, с. 29]. Наиболее эффективными и устойчивыми среди итерационных методов решения таких систем уравнений являются так называемые проекционные методы, и особенно тот их класс, который связан с проектированием на подпространства Крылова. Эти методы обладают целым рядом достоинств: они устойчивы, допускают эффективное распараллеливание, работу с различными строчными (столбцовыми) форматами и предобуславливателями разных типов.

Рассмотрим систему линейных алгебраических уравнений

, (1)

где: — плохо обусловленная разреженная матрица,

.

В данной работе проводится сравнительный анализ итерацион­ных методов для решения плохо обусловленных разреженных СЛАУ. В качестве исследуемых методов выбраны: метод сопряженных градиентов (CG) , метод минимальных невязок (MinRes), сдвоенный метод сопряженных направлений (CGS), квазиминимальных невязок (QMR) [6, с. 2].

В вопросах выбора того или иного способа решения СЛАУ важно учитывать структуру матрицы A [4, с. 2; 5, с. 47]. Это связано с тем, что не каждый метод дает возможность получить гарантированный результат для определенной системы линейных уравнений.

Таким образом, критерием сравнения итерационных методов решения СЛАУ будут: погрешность результатов, скорость сходимости, структура матрицы.

Результаты численных исследований показали, что для решения СЛАУ с матрицей A являющейся симметричной/несимметричной и хорошо обусловленной к нормальным уравнениям лучше применять метод CG. Если матрица А- симметричная и плохо обусловленная, то лучшую сходимость показал метод MinRes. Для А- несимметрич­ной, плохо обусловленной — метод квазиминимальных невязок [7].

Нужна помощь в написании статьи?

Мы — биржа профессиональных авторов (преподавателей и доцентов вузов). Пишем статьи РИНЦ, ВАК, Scopus. Помогаем в публикации. Правки вносим бесплатно.

Для улучшения скорости сходимости итерационных методов используют предобуславливание матрицы системы. Оно заключается в том, что подбирается такая матрица предообуславливания, что при этом процедура решения СЛАУ является не слишком трудоемкой и численно устойчивой [2, с. 313]. Правильный выбор предобуславливателя, зависящий от конкретной задачи, может в огромной степени ускорить сходимость [9, с. 211]. В действитель­ности, хороший предобуславливатель часто необходим для того, чтобы итерационный метод вообще сходился.

В данной работе было рассмотрено несколько видов предобус­лавливания для метода квазиминимальных невязок с разреженными плохо обусловленными матрицами: левое и правое предобус­лавливание с использованием QR — разложения, левое и правое предобуславливание с использованием LU — разложения, а также с использованием модификации LU — разложения [8, с. 315].

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

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

Нужна помощь в написании статьи?

Мы — биржа профессиональных авторов (преподавателей и доцентов вузов). Пишем статьи РИНЦ, ВАК, Scopus. Помогаем в публикации. Правки вносим бесплатно.

1.Голуб Дж., Ван Лоун Ч. Матричные вычисления/ Под ред. В.В. Воеводина. — М.: «Мир», 1999. — 548 с.

2.Деммель Дж. Вычислительная линейная алгебра. Теория и приложения / Пер. с англ.Х.Д. Икрамова. — М.: «Мир», 2001. — 430 c.

3.Писсанецки С. Технология разреженных матриц/ Под ред. Х.Д. Икрамова — М.: «Мир», 1988. — 410 с.

4.Станкевич, И.В. Хранение и использование разреженных матриц в конечно- элементной технологии. Журнал «Наука и Образование». — 2005. — 10 октября.

5.Тьюарсон Р. Разреженные матрицы/ Под ред. Х.Д. Икрамова. — М.: «Мир», 1977. — 172 с.

6.Bucker Martin, Basermann Achim. A comparison of QMR, CGS and TFQMR on a distributed memory machine / Bucker Martin //Mathematics of computation. — 1994 — 31 may

7.Harwell-Boeing Collection — [Электронный ресурс] – Режим доступа. — URL: http://math.nist.gov/MatrixMarket/data/Harwell-Boeing/ (дата обращения: 15.12.2012)

8.Roland W. Freund, Noel M. Nachtigal. QMR: a Quasi-Minimal Residual Method for Non—Hermitian Linear Systems / Roland W. Freund, Noel M. Nachtigal // Journal Math. — 1991. — №60. — p. 315—339.

9.Saad, Y. Iterative methods for sparse linear systems / Y. Saad. // SIAM. — 2003. — 447 p.

Книга: Решение систем линейных алгебраических уравнений для больших разреженных матриц

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ ОБРАЗОВАТЕЛЬНОЕ

УЧРЕЖДЕНИЕ ВЫСШЕГО ПРОФЕССИОНАЛЬНОГO

“ЮЖНЫЙ ФЕДЕРАЛЬНЫЙ УНИВЕРСИТЕТ”

РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ

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

В.А.Еремеев. Решение систем линейных алгебраических уравнений для больших разреженных матриц: Учебно-методическое пособие. – г.Ростов-на-Дону, 2008 г. – 39 с.

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

1. задачи линейной алгебры с разреженными матрицами на примере дискретизации уравнения Пуассона;

2. векторные и матричные нормы;

3. итерационные методы;

4. методы подпространств Крылова.

Пособие предназначено для преподавания дисциплин в рамках магистерской образовательной программы “Математическое моделирование и компьютерная механика”.

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

Пособие подготовлено в рамках гранта ЮФУ K-07-T-41.

1 Задачи линейной алгебры с разреженными матрицами на примере дискретизации уравнения Пуассона 4

1.1 Решение уравнения −x 00 = f (t ) . . . . . . . . . . . . . . . 4

1.2 Решение уравнения Пуассона . . . . . . . . . . . . . . . . 6

2 Векторные и матричные нормы 11

2.1 Векторные нормы . . . . . . . . . . . . . . . . . . . . . . 11

2.2 Матричные нормы . . . . . . . . . . . . . . . . . . . . . . 13

2.3 Связь векторных и матричных норм . . . . . . . . . . . 15

3 Итерационные методы 19

3.1 Определения и условия сходимости итерационных методов 19

3.2 Метод простой итерации . . . . . . . . . . . . . . . . . . . 23

3.3 Метод Якоби . . . . . . . . . . . . . . . . . . . . . . . . . 24

3.4 Метод Гаусса-Зейделя . . . . . . . . . . . . . . . . . . . . 25

3.5 Метод последовательной верхней релаксации (SOR) . . . 26

3.6 Метод симметричной последовательной верхней релак-

сации (SSOR) . . . . . . . . . . . . . . . . . . . . . . . . . 27

4 Методы подпространств Крылова 30

4.1 Инвариантные подпространства . . . . . . . . . . . . . . 30

4.2 Степенная последовательность и подпространства Кры-

4.3 Итерационные методы подпространств Крылова . . . . . 33

Дополнительная литература 39

1 Задачи линейной алгебры с разреженными матрицами на примере дискретизации уравнения Пуассона

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

1.1 Решение уравнения −x 00 = f (t )

Рассмотрим простейшую краевую задачу

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

так, чтобы t ₀ = 0, t_{n +1} = 1. Используем обозначения

Для дискретизации второй производной воспользуемся центральной конечной разностью

,

которая аппроксимирует x 00 (t_i ) c точностью O (h 2 ). Тогда исходная краевая задача сводится к системе линейных алгебраических уравнений (СЛАУ) относительно

,

.

Эту СЛАУ можно записать также в матричном виде

Обратим внимание на следующие свойства матрицы A .

• A – разреженная матрица, она трехдиагональная;

• A – симметричная и, более того, положительно определенная. Это свойство она унаследовала от свойств оператора исходной краевой задачи;

• Структура A достаточно простая, ее элементы вычисляются по простым формулам. Следовательно, для ее хранения в памяти компьютера нет необходимости использовать типы данных, предназначенные для хранения заполненных матриц.

Обсудим размерность полученной СЛАУ. Предположим, что требуемая точность аппроксимации равна ² = 10 −6 . Поскольку ² ∼ h 2 , то

−3 , т.к. h ∼ √² . Т.е. требуемый шаг сетки нужно выбрать порядка 10 мы имеем матрицу размерности 10 3 × 10 3 . На практике, размерность должна быть выше, чтобы удовлетворить большей точности.

Рис. 1: Конечно-разностная сетка. Естественная нумерация узлов.

1.2 Решение уравнения Пуассона Рассмотрим более сложную краевую задачу

где u = u (x,y ) – неизвестная функция, f (x,y ) – известная, заданные на единичном квадрате (x,y ) ∈ Ω = [0, 1]×[0, 1]>, Γ = ∂ Ω – граница Ω. Введем на Ω равномерную сетку, так что координаты узлов даются формулами x_i = hi, y_i = hi, i = 0. N + 1, h = 1/ (N + 1).

Для дискретизации вторых производных в операторе Лапласа ∆ воспользуемся формулами для центральной конечной разности

,

Таким образом, исходная краевая задача сводится к СЛАУ

−u i −1,j + 4u i,j − u i +1,j − u i,j −1 − u i,j +1 = h 2f i,j (1) относительно u_i,j . Для того, чтобы свести ее к стандартному виду A x = b, необходимо преобразовать u_i,j к выражению с одним индексом, т.е. перенумеровать каким-то образом. Выбор нумерации влияет, вообще говоря, на структуру матрицы A .

Рассмотрим случай N = 3. Соответствующая сетка показана на рис. 1, где использована естественная нумерация внутренних узлов. Полые кружки соответствуют граничным узлам, в которых известно значение функции, а сплошные – внутренним.

Если преобразовать конечно-разностное уравнение (1) в матричное уравнение A x = b, вводя вектор неизвестных по правилу u 1, 1 = x 1, u 1, 2 = x 2, u 1, 3 = x 3, u 2, 1 = x 4. u 3, 3 = x 9,

Видно, что A является симметричной ленточной разреженной матрицей с диагональным преобладанием.

В результате нужно отметить, что в целом, матрица A обладает теми же свойствами, что и в случае одномерной, задачи, хотя ее структура может не быть ленточной, что зависит от способа нумерации узлов сетки

Обсудим размерность полученной СЛАУ. Предположим как и ранее, что требуемая точность аппроксимации равна ² = 10 −6 . Поскольку ² ∼ h 2 , то требуемый шаг сетки нужно выбрать порядка 10 −3 , т.к.

√ 3 × 10 3 = 10 6 . Таh ∼ ² . Число узлов здесь оказывается равным 10 кому же числу равно число неизвестных n . Таким образом, мы имеем матрицу размерности 10 6 × 10 6 .

Нетрудно догадаться, что если рассмотреть, не квадрат, а куб, то число неизвестных будет примерно равно n = 10 9 . В расчетной практике встречаются задачи, где нужно определять не одну функцию, а несколько. Например, в гидродинамике при учете теплопереноса имеем четыре скалярных неизвестных (три компоненты поля скорости и поле температуры). Если рассматривать конечно-разностную аппроксимацию этой задачи с точностью ² = 10 −6 , то получим размерность n = 4 · 10 9 .

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

Рассмотрим другую нумерацию узлов, показанную на рис. 2. Здесь использовано так называемое красно-черное разбиение (или чернобелое). Вначале нумеруются узлы, сумма индексов которых – четное число, а потом узлы, сумма индексов которых дает нечетное число. Таким образом, вектор переменных вводится по правилу вектор неизвестных по правилу

– это красные неизвестные, а

Рис. 2: Конечно-разностная сетка. Красно-черное разбиение.

Соответствующая красно-черному разбиению матрица дается формулой