Методы решения больших разреженных систем уравнений

Методы решения больших разреженных систем уравнений

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

Выводы.

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

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

Решение разреженных СЛАУ больших размерностей средствами ManagedCuda в .NET

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

В статье будет рассказано о том, как ваш покорный слуга значительно повысил эффективность компьютерной модели расчета нестационарных течений газа в крупных системах газоснабжения произвольной конфигурации, благодаря применения библиотеки ManagedCuda и nVidia CUDA 7.0. Однако изложение будет вестись без привязки к конкретной предметной области.

Постановка задачи

Рассматривается классическая задача решения СЛАУ:
Ax=b
Здесь матрица A размером nxn, x — вектор неизвестных размером n, b — вектор известных свободных членов размером n.

Автор статьи занимается разработкой специализированного программно-вычислительного комплекса (ПВК) для моделирования и оптимизации нестационарных течений газа в системах газоснабжения. В решаемых мною задачах A является положительно определенным якобианом с диагональным преобладанием некоторой системы нелинейных уравнений. A получается в результате матричных преобразований матрицы инциденций графа системы газоснабжения и других матриц. В моих практических расчетах A обычно заполнена на 6-10%, остальное нули. Размер же ее зависит от размера конкретной системы газоснабжения и варьируется n от 10 до 1000.

Наш ПВК разрабатывается под .NET 4.0 на C#. Расчетный модуль и вся математика также разрабатывается на C#. Изначально для решения СЛАУ я написал собственную, доморощенную, библиотеку, не использующую технологию разреженных матриц. Для решения СЛАУ применял метод LU декомпозиции. И до поры меня все устраивало, пока я не начал заниматься задачей оптимизацией нестационарных режимов течения газа, где приходится методом динамического программирования осуществлять большое количество переборов значений управляющих параметров и, соответственно, много раз решать СЛАУ. Стандартный профилировщик Visual Studio показал, что в процессе выполнения программы около 40% затрат приходится на решение СЛАУ.

Именно в этот момент я понял, что пора что-то менять.

Математическая библиотека Math.Net Numerics

Проведя анализ имеющихся математических библиотек под .Net я решил остановиться на библиотеке Math.Net Numerics. Подробно ознакомиться с ней вы можете здесь.

Меня заинтересовали ее ключевые возможности:

• Реализована технология разреженные матриц и векторов;
• Имеются встроенные все необходимые вектор-матричные операции;
• Присутствуют встроенные решатели;
• Интуитивно понятный API.

Приведу листинг с примером решения СЛАУ средствами Math.Net Numerics:

Как видно, все очень элегантно выглядит, но практически сразу меня настигло разочарование — встроенный решатель Math.Net Numerics работал значительно медленнее моего. Меня это совершенно не устроило, однако претензии к реализованным в библиотеке вектор-матричным операциям не столь существенные. Поэтому полностью отказываться от Math.Net Numerics я не стал, оставив код, где производится работа с векторами и матрицами.

Но тут мне вспомнился весьма удачный опыт коллег по аспирантуре, которые для решения задач подземной гидродинамики использовали графические процессоры. В распоряжении у меня имеется ноутбук с видеокартой nVidia GeeForce GT 540M с 2 GB памяти и поддержкой технологии CUDA. Было принято решение попробовать эту технологию на деле, о чем теперь не жалею.

Библиотека ManagedCuda

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

Меня заинтересовала библиотека cuSPARSE. При первом знакомстве с ней я наткнулся на проблемы:

• Непосредственная работа с библиотекой возможна только через C/C++. Конечно, проблема решится, если использовать качественный враппер, который очень не хотелось писать. По результатам поиска был найден CUDA-враппер для .Net под названием Cudafy.
• cuSPARCE позволяет решать СЛАУ с треугольными матрицами, а в моем случае матрицы таковыми не являются. Однако cuSPARSE содержит методы факторизации матриц, приводящие их к треугольной форме (метод LU факторизации, разложение Холецкого). Однако в API Cudafy отсутствовали соответствующие методы, поэтому от Cudafy пришлось отказаться.

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

Идея заключается в использовании SparseMatrix из Math.Net Numerics для формирования и хранения разреженной матрицы в формате CSR, который принимают на вход cuSPARSE и ManagedCuda. Далее приводится листинг соответствующей программы:

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

Дабы не утомлять читателя сухим текстом и листингами программ, рассмотрим результаты расчетов СЛАУ размерностей от 50 до 500 с помощью различных технологий:

• Доморощенный решатель, написанный автором статьи. Назовем его Mani.Net.
• Стандартный решатель библиотеки Math.Net Numerics.
• Метод LU декомпозиции библиотеки ManagedCuda.

Матрица A заполняется на 10% случайным образом.

Рисунок наглядно демонстрирует, почему мне пришлось отказаться от Math.Net Numerics для решения СЛАУ.
Отметим, при размерности матрицы равной 500 моему собственному решателю (Mani.Net) потребовалось 1170 мс, Math.Net Numerics — 12968 мс, ManagedCuda — 70 мс.

Заключение

Следует ожидать в комментариях замечания, мол де не на всех машинах есть GPU nVidia с поддержкой CUDA. Действительно, это так. Поэтому наше приложение настроено на две конфигурации компиляции: с собственной библиотекой решения СЛАУ и с ManagedCuda.

Насчет Mani.Net отмечу, что это не реклама моей собственной библиотеки. Найти ее нигде невозможно и никому её я передавать не буду. Нет, я не жадный. Мне стыдно за код.

Спасибо за прочтение статьи! Буду рад узнать о вашем мнении и замечаниях в комментариях.

Численное решение больших разреженных систем уравнений, Джордж А., Лю Дж., 1984

Численное решение больших разреженных систем уравнений, Джордж А., Лю Дж., 1984.

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

Положительно определенные и неопределенные матричные задачи.

В этой книге мы будем иметь дело исключительно со случаем, когда А симметрична и положительно определена. Как уже было отмечено, существенная часть линейных систем, возникающих в научных и инженерных расчетах, обладает этим свойством, и проблема упорядочения для них решается иначе и проще, чем для разреженной матрицы А общего вида. В последнем случае необходима для обеспечения численной устойчивости та или иная форма выбора главного элемента, т. е. перестановки строк и/или столбцов iForsythe 1967). Таким образом, при заданной А обычно получают разложение для РА или PAQ, где Р н Q —матрицы перестановок соответствующих размеров. (Заметим, что умножение на Р слева переставляет строки 4, а умножение на Q справа переставляет столбцы А.)
Эти перестановки определяются в процессе разложения путем компромисса между (обычно конкурирующими) требованиями численной устойчивости и разреженности (Duff 1974). Различные матрицы, хотя бы они и имели одинаковую структуру нулей-ненулей, обычно приводят к различным Р и Q и, следовательно, имеют множители с различной структурой разреженности. Другими словами, для разреженных матриц общего вида, как правило, нельзя предсказать, 1де произойдет заполнение, пока не начались собственно вычисления. Тем самым мы вынуждены пользоваться какой-либо схемой динамического хранения, в которой память для заполнения выделяется в ходе вычислений.

Содержание.

От переводчика
Предисловие
Глава 1. Введение.
Глава 2. Вводные сведения.
Глава 3. Некоторые сведения из теории графов и ее применение к исследованию разреженных симметричных матриц.
Глава 4. Ленточные и профильные методы.
Глава 5. Универсальные разреженные методы.
Глава 6. Методы фактор-деревьев для конечноэлементных и конечноразностных задач.
Глава 7. Методы параллельных сечений для конечноэлементных задач
Глава 8. Методы вложенных сечений.
Глава 9. Численные эксперименты.
Приложение А. Некоторые указания к использованию подпрограмм.
Приложение В. SPARSPAK: Пакет для разреженных матриц.
Литература.

Бесплатно скачать электронную книгу в удобном формате, смотреть и читать:
Скачать книгу Численное решение больших разреженных систем уравнений, Джордж А., Лю Дж., 1984 — fileskachat.com, быстрое и бесплатное скачивание.

Скачать djvu
Ниже можно купить эту книгу по лучшей цене со скидкой с доставкой по всей России. Купить эту книгу