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

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

Во многих областях науки и техники, при решении различных задач, возникает необходимость расчета систем линейных уравнений, которые иногда являются сильно разреженными. Поскольку обычно такие системы уравнений имеют большой порядок, то найти их решение (в сколь-нибудь обозримые сроки) возможно только с применением вычислительной техники. Для решения подобных задач можно использовать прямые методы (Гаусса, Крамера), и приближенные (Зейделя, градиентного спуска и др.). Прямые методы дают точное решение (ошибка может быть только из-за выхода величин при преобразованиях за пределы точности типов для используемой платформы), но обычно требуют большее количество итераций. Приближенные методы обычно работают быстрее, но могут требовать особых начальных условий или какой-либо подготовки исходных данных, и решение не всегда сходится. В данной статье рассматривается реализация метода Гаусса, модифицированная для решения разреженных систем линейных уравнений. Данный метод не является «математическим», а скорее «алгоритмическим» и показывает, как путем простых приемов можно ускорить расчет в сотни раз.

Условия применения

Рассматриваемый метод подходит для любых невырожденных (имеющих одно решение) систем линейных уравнений. В дальнейшем под словом матрица будет пониматься матрица коэффициентов линейных уравнений. Единственное условие — метод имеет смысл только для сильно разреженных матриц больших размерностей N>1000, где не менее двух третей коэффициентов равны 0. Для обычных матриц метод не даст выигрыша ни в скорости, ни в экономии памяти.

Описание метода

Метод Гаусса заключается в обнулении нижней (или верхней) полуматрицы a[i,k], где строка i=2..N, столбец k=1..i-1. Используя известное правило, что решение системы не изменится, если к одному из ее уравнений добавить другое, умноженное на произвольный коэффициент, строки матрицы поочередно, начиная с первой, домножаются на специальные коэффициенты и суммируются со всеми ниже лежащими строками. Коэффициенты подбираются таким образом, чтобы каждая строка обнуляла соответствующий ей (i=k) ниже лежащий столбец (только его поддиагональную часть). Так первая строка обнуляет первый столбец A[2..N,1], вторая – второй A[3..N,2] и т.д. до предпоследней строки. Коэффициенты рассчитываются так: например для того чтобы обнулить первый столбец нужно первую строку домножить на k=-A[2,1]/A[1,1] и прибавить ко второй, затем на k=-A[3,1]/A[1,1] и прибавить к третей и т.д. После операция повторяется со второй строкой (она домножается на соответствующие коэффициенты и складывается со всеми нижележащими строками) и т.д. до предпоследней. В результате получается матрица, под главной диагональю которой остаются одни нули. Теперь в последнем уравнении остался только один коэффициент – соответственно можно легко вычислить значение последней неизвестной: X[N]=B[N]/A[N,N], где B – вектор свободных членов уравнений, который кстати, также участвует в преобразованиях при обнулении матрицы A. В предпоследнем уравнении осталось два коэффициента, и последняя неизвестная уже известна, можно вычислить предпоследнюю: X[N-1]=(B[N-1]-A[N-1,N]*X[N])/A[N-1,N-1]. Таким образом, используя обратный ход, можно вычислить все неизвестные. В процессе преобразования матрицы не должны обнуляться коэффициенты главной диагонали. Если такое происходит, то матрица не была факторизована или неправильно составлены уравнения. Или же данная система имеет не одно решение. Более подробно о методе Гаусса можно узнать здесь.

Приведенный выше алгоритм легко реализуем. Количество итераций известно заранее и зависит только от порядка матрицы. Если его запрограммировать в классическом стиле, то он легко будет справляться с небольшими матрицами N

1000. Однако матрицы иногда бывают порядка 100 000 и более. На расчет такой матрицы классическим методом Гаусса уйдут годы, поскольку количество итераций будет огромным. Если данную матрицу хранить в обычном двумерном массиве, то при использовании типа double она будет занимать в памяти 100 000 х 100 000 х 8 = 800 000 000 байт = 800 Мб что совсем не мало. А если учесть что размер матрицы возрастает пропорционально квадрату порядка матрицы, то для матрицы порядка 400 000 это будет в 16 раз больше, т.е. 12.8 Гб!

На основании выше приведенных рассуждений можно выяснить два направления оптимизации расчета:

• Более компактное хранение информации.
• Уменьшение количества итераций.

Для разреженных матриц возможно и то, и другое:

Более компактное хранение информации — разреженные матрицы имеет смысл хранить в памяти не в виде двумерного массива, а в виде некоторой структуры, в которой сохранены только ненулевые значения коэффициентов, а также их координаты в матрице. Уменьшение количества итераций — если выше упомянутую структуру организовать особым образом, чтобы из нее можно было быстро получить нужную информацию в нужном для метода Гаусса виде, то можно сократить кол-во итераций до минимального (ведь отпадает необходимость перебирать нулевые коэффициенты, которых в матрице может быть достаточно много).

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

• для главной диагонали это одна координата d, которая равна соответствующему индексу элемента главной диагонали
• для верхней и нижней полуматрицы координаты d,j, где j – смещение от главной диагонали вправо или вниз (как показано на рисунке)

Диаграмма классов

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

Реализация метода

Все исходные коды приведены на C#, и представляют полный набор классов, необходимых для функционирования метода. Код написан максимально просто для понимания принципа работы метода, и поэтому не претендует на максимальную производительность и не является оптимальным. О необходимых доработках сказано в заключении. Код снабжен комментариями и не нуждается в дополнительном разъяснении.

Реализация интерфейса строк и столбцов:

Реализация интерфейса матрицы:

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

На данном этапе нужно решить в каком формате хранить значения ячеек матрицы. Для очень сильно разреженных матриц, когда в каждой строке/столбце матрицы ненулевых значений совсем мало, данные лучше хранить в обычных массивах, так как поиск в небольшом массиве работает быстро. Если же матрица не сильно разрежена, то лучше использовать хэш-таблицу. Ниже приведены реализации интерфейса IRow для обоих случаев:

Реализация интерфейса строки c использованием массива для хранения значений

IРеализация интерфейса строки c использованием хэш-таблицы для хранения значений:

Полная реализация класса «матрица»:

И, наконец, реализация метода Гаусса:

Заключение

В статье рассмотрен метод решения систем линейных уравнений по алгоритму Гаусса, и приведен код реализующий этот алгоритм для C#. Хотя данный код и является полным, для использования в пользовательских приложениях он нуждается в следующей доработке:

• Как минимум нужно сделать обработку исключительных ситуаций, например таких как деление на ноль или выход индекса за пределы массива
• Реализовать расчетный алгоритм в отдельном потоке команд. Если вызывать метод calculate из основного UI потока, это остановит прием сообщений окружения («повесит программу») до момента завершения расчета
• Если необходимо рассчитать несколько систем линейных уравнений, которые отличаются только правой частью (свободными членами), то естественно нет необходимости каждый раз обнулять нижнюю полуматрицу. Для этого нужно доработать класс Gauss таким образом, чтобы он, обнуляя матрицу, сохранял историю изменений. После этого для расчета неизвестных нужно будет преобразовать только вектор свободных членов, передаваемых в метод calculate в соответствии с историей. Историю можно хранить в отдельной матрице, или в нижней полуматрице основной матрицы, которая все равно не используется при обратном ходе
• Если данный алгоритм планируется реализовать на платформе, где возможна прямая работа с указателями (Delphi, C++), то лучше по возможности использовать их, это позволит в некоторых местах выполнять меньшее количество преобразований
• Для некоторых матриц необходима доработка алгоритма факторизации, для исключения попадания нулей на главную диагональ
• В реализациях интерфейса IRow, при создании экземпляров массивов и хэш-таблиц, в явном виде не задавалось свойство capacity. На практике, чтобы оптимально использовать быстродействие и экономить память, нужно подобрать компромисное значение capacity. Так с увеличением capacity увеличивается быстродействие и расход памяти, с уменьшением – наоборот. Оптимальное значение capacity зависит от структуры матрицы и не превышает максимального количества ненулевых элементов в строке или столбце матрицы

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

Литература

1. Эндрю Троелсен. С# и платформа .NET. Библиотека программиста. — СПб.: Питер, 2004.

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

Аннотация научной статьи по математике, автор научной работы — Чадов Сергей Николаевич

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

Похожие темы научных работ по математике , автор научной работы — Чадов Сергей Николаевич

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

О РЕШЕНИИ РАЗРЕЖЕННЫХ СИСТЕМ ЛИНЕЙНЫХ УРАВНЕНИЙ ПРИ ПОМОЩИ СТАБИЛИЗИРОВАННОГО МЕТОДА БИСОПРЯЖЕННЫХ ГРАДИЕНТОВ

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

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

THE SOLUTION OF LOW-DENSITY LINEAR EQUATION SYSTEMS WITH THE HELP OF STABILIZED BI-CONJUGATE GRADIENT METHOD

CHADOV S.N., postgraduate.

The article concerns one of the iterative methods of linear equation systems solution. The article contains the description of the method, some ways of its paralleling and the advantages of this method against more common ones.

Key words: linear equation system, high dimensionality systems, fore-stipulation matrix, parallel realization, biconjugate gradient method.

Введение. Решение систем линейных уравнений является задачей, имеющей применение как в чистом виде, так и в качестве примитива в более сложных численных алгоритмах. В качестве примера можно привести решение систем дифференциальных уравнений с использованием неявных схем [1, 2].

Существует множество методов решения систем линейных уравнений, каждый из которых имеет свои достоинства и недостатки. Наиболее часто на практике используется один из вариантов метода LU-разложения (см., например, [1]). Однако в реальных задачах при моделировании сложных систем (например, в энергетике) приходится иметь дело с системами очень большой размерности (десятки и сотни тысяч элементов и более), но при этом часто эти системы являются достаточно разреженными. Для таких систем распространенные методы не всегда являются самыми оптимальными, поэтому часто используют итерационные методы [3]. Рассмотрим один из таких методов, известный под названием BiCG-Stab [3]. Этот метод был предложен H. van der Vorst в 1992 г. и является модификацией метода бисопряженных градиентов, обеспечивающей лучшую сходимость [3, 4] .

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

Постановка задачи и описание метода. Итак, пусть требуется решить систему уравнений Ax = b, (1)

где х — искомый вектор; A — некоторая матрица.

Общая идея метода бисопряженных градиентов для решения линейной системы уравнений (1) заключается в итеративной генерации пары остатков r и r, вычисляемых как r(i) = r(i-1) -a^p(i)

и r(i> = r(i-1) -aiATp(i»,

где p — векторы направлений: p(i) = r(i-1) +pi-1p(i-1),

p(i) = p(i-1) +pi-1p(i-1); a и p выбираются таким образом, чтобы выполнялись отношения ортогональности (r(i))Tr(j) = (p(i))TAp(j) = 0 для i t j:

p(i-1)T г(і-1) f(i)T г(і)

ai= p(i)T Ap(i) , e = f(i-1)T г(і-1).

Можно заметить, что r(j) и f(j) выразимо в виде г() = P^A)^, г() = P(j)(AT )f(0)

Из отношения ортогональности получаем

где (■) — оператор скалярного произведения; P — некоторый полином.

Таким образом, используя последнее выражение, возможно полностью избежать транспонирования

матрицы А, а также вычисления вектора r(j).

H. van der Vorst в 1992 г. развил эту идею, предложив представить r(j) в виде

где Q(j)(x) = (1 -ro1x)(1 —

Следует отметить использование матрицы M в качестве матрицы предобусловливания (preconditioner). Использование предобусловливания во многих случаях позволяет на порядки сократить время вычислений. Причина заключается в том, что скорость сходимости итеративных методов зависит от спектральных характеристик матрицы коэффициентов. Соответственно, можно попытаться так преобразовать систему, чтобы она имела те же решения, что и исходная, но при этом имела бы лучшую (относительно метода решения) спектральную характеристику. Выбор матрицы предобусловливания часто является достаточно сложной творческой задачей, поскольку на настоящий момент неизвестно удобных на практике формальных алгоритмов для такого выбора (хотя существует несколько наиболее распространенных вариантов [3]). Одним из наиболее распространенных является реализация предобусловливателя на основе алгоритма неполного LU-разложения [3, б]

Вычислить г(0) = Ь- Ах(0) для некоторого начального значения х(0) г = г(0)

если рі-1 == 0 -> вернуть ошибку

если і==1, р(і) = г(і-1)

Рі-1 = (Рі-1/ Рі-2)(“і-1/ “і-1) p(i) = г(і-1) + Рм(р(і-1) -юь/-1))

Решить Мр’ = р(і) у(і) = Ap’

аі = Рі-1 / г V0 s = г(і-1) — а іv(i )

если (норма б достаточно мала), Xі = хі-1 + аір’ и остановиться

x(i) = x(i-1) + а ip + ю is’

Проверить на сходимость и, если она не достигнута и юі £0, выполнить следующую итерацию

Алгоритм метода Bicg-STAB

Предобусловливание может не использоваться вообще, тогда в алгоритме следует заменить шаги «Решить Мр’ = р(|)» и «Решить Мб’ = б» на шаги «р’=р» и «б’=б», соответственно.

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

1) не требует серьезной модификации существующего программного кода, при условии, что этот код грамотно спроектирован;

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

3) масштабируется практически на любое количество потоков.

Описанный подход, ввиду своей простоты, находит достаточно широкое применение (см., например, [2]).

Перспективной представляется реализация такого алгоритма на машине с общей памятью и использованием векторных команд процессора (например, современные процессоры персональных компьютеров поддерживают наборы Э!МО-инструкций БЭБ, Э8е2 и так далее, позволяющие выполнять две операции с числами с двойной точностью или 4 команды с числами с одинарной точностью одной командой).

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

V0) = Ap = A(r(M) + Рм(р0-1) — Ю|_УМ))) =

где u(i-1) = Ar(i-1); q(i-1) = Au(i-1); для остальных используемых векторов — аналогично.

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

Отдельную задачу представляет собой параллельная реализация предобусловливателя. Эксперименты показывают, что на решение предобу-словливающей системы уравнений уходит не менее 30 % времени вычислений. Разработать какой-либо универсальный алгоритм распараллеливания в данном случае невозможно, поскольку шаг вычисления предобусловливателя может выполняться совершенно произвольным образом, зависящим от конкретной решаемой задачи. Однако приведем вариант параллельной реализации наиболее часто используемого метода неполного LU-разложения.

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

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

Рассмотренный метод поддается распараллеливанию и соответственно может быть использован в пакетах параллельных программ. Возможным слабым местом в реализации метода может являть-

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

Данный вариант метода BICG-Stab реализован как часть комплекса программ для решения систем жестких дифференциальных уравнений (применяется для вычисления корректора неявной схемы интегрирования) как в последовательном, так и в параллельном вариантах. Параллельные варианты были разработаны для многопроцессорных систем с общей памятью (с использованием OpenMP) и на основе SIMD-инструкций процессоров x86.

1.William H. Press . [et al]. Numerical recipes in C: the art of scientific computing . 2nd ed ISBN 0-521-43108-5.

2.Hindmarsh A.C., Serban R. User Documentation for cvode v2.4.0. Center for Applied Scientific Computing Lawrence Livermore National Laboratory, 2006.

3. Barrett R., Berry M., Chan T.F., et al. Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, 2nd Edition.

4.Van der Vorst H. Linear System Solvers: Sparse Iterative methods.

5.Chan T.F., Van der Vorst H. Approximate and incomplete factorizations.

6.Ouarraui С., Kaeli D. Developing object-oriented parallel iterative methods.Department of Electrical and Computer Engineering,Northeastern University, Boston, MA 02115.

7.Gonzalez P., Cabalero J., Pena T. Parallel Incomplete LU factorization as a Preconditioner for Krylov Subspace methods.

Огарёв-online

Электронное периодическое издание для студентов и аспирантов

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

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

ALGORITHM OF SOLVING A SPARSE SYSTEM OF LINEAR ALGEBRAIC EQUATIONS OF LARGE DIMENSION BY CONJUGATE GRADIENT METHOD

The algorithm of solving sparse systems of linear algebraic equations of large dimension of a special type by using of the conjugate gradient method is described. The results of a computational experiment for different values of the input data are presented.