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

Коньшин, Игорь Николаевич

Параллельные методы решения систем линейных уравнений с симметричными положительно-определенными матрицами на основе аддитивного разложения с перекрытиями тема автореферата и диссертации по математике, 01.01.07 ВАК РФ

Коньшин, Игорь Николаевич АВТОР

кандидата физико-математических наук УЧЕНАЯ СТЕПЕНЬ

Москва МЕСТО ЗАЩИТЫ

2009 ГОД ЗАЩИТЫ

01.01.07 КОД ВАК РФ

Диссертация по математике на тему «Параллельные методы решения систем линейных уравнений с симметричными положительно-определенными матрицами на основе аддитивного разложения с перекрытиями»

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

На правах рукописи _

Коньшин Игорь Николаевич

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

01.01.07 — Вычислительная математика

Автореферат 1 ?

диссертации на соискание учёной степени кандидата физико-математических наук

Москва — 2009

003482890

Работа выполнена в Учреждении Российской академии наук Вычислительный ценр им. А. А. Дородницына РАН.

кандидат физико-математических наук Капорин И. Е.

доктор физико-математических наук,

профессор

Кобельков Г. М.

доктор физико-математических наук, доцент

Василевский Ю. В.

Научно-исследовательский вычислительный центр МГУ

Защита состоится" / &" ^-^¿^^о^ООЭ в ч. && мин. на заседании Диссертационного совета Д 002.017.01 при Учреждении Российской академии наук Вычислительный ценр им. А. А. Дородницына РАН по адресу: 119333, г. Москва, ул. Вавилова, дом. 40, тел: (499) 135-24-89.

С диссертацией можно ознакомиться в библиотеке ВЦ РАН.

Научный руководитель: Официальные оппоненты:

Ведущая организация:

Автореферат разослан

Г- 2009.

о о V

./ V)

Учёный секретарь Диссертационного совета V

доктор физико-математических наук, профессор Зубов В. И.

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Объект исследования и актуальность темы. Актуальность темы диссертационной работы обусловлена необходимостью эффективного численного решения систем линейных алгебраических уравнений с разреженными симметричными положительно-определенными матрицами большой размерности. Численное решение систем линейных уравнений является одной из наиболее часто встречающихся задач в научно-технических исследованиях и практических приложениях. Подобные задачи возникают, например, в математической физике при численном решении дифференциальных и интегральных уравнений. При этом прикладные задачи часто требуют решения систем линейных уравнений с большим числом неизвестных. К таким системам, например, приводит численное решение двумерных и особенно трехмерных задач математической физики, в которых условия физической и геометрической аппроксимации двумерной и трехмерной области диктуют использование достаточно мелкой расчетной сетки с большим числом расчетных узлов. Для возникающих систем линейных уравнений применение прямых методов решения оказывается неприемлемым как по причине большого времени решения, так и недостаточности обьема оперативной памяти для хранения данных задачи. Итерационные методы решения систем линейных уравнений намного экономичнее с точки зрения использования оперативной памяти, но для ускорения их сходимости требуется строить эффективные предобусловливания, особенно при решении систем с плохо обусловленными матрицами.

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

влечения информации о пространственном расположении узлов расчетной сетки.

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

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

В соответствии с целью исследования были поставлены следующие задачи:

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

2. Разработка эффективного способа балансировки вычислений для подзадач на различных процессорах.

3. Исследование эффективности метода в зависимости от параметров прсдобусловливапия и количества процессоров.

4. Подтверждение эффективности метода для различных типов решаемых задач и на различных вычислительных системах.

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

Научная новизна:

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

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

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

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

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

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

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

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

Положения, выносимые на защиту:

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

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

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

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

Апробация работы. Результаты работы докладывались и обсуждались па международной конференции "РаСТ99'! (Санкт-Петербург, Россия, 1999 г.), Всемирном 16-ом Конгрессе "IMACS 2000" по научным вычислениям, прикладной математике и моделированию (Лозанна, Швейцария, 2000), Голландско-российском симпозиуме NWO (МГУ, Москва, июнь 200Q), Симпозиуме NWO (Амстсрдам-Ниймсгсн, Голландия, ноябрь 2000),

Втором Международном научно-практическом Семинаре и Всероссийской молодежной школе "Высокопроизводительные параллельные вычисления на кластерных системах" (Нижний Новгород, Россия, 2002), Международной конференции "Parallel CFD 2003" (Москва, 2003), Международной конференции "VIII Забабахинские научные чтения" (РФЯЦ-ВНИИТФ, Сне-жинск, Россия, 2005), Международной конференции по численной геометрии, построении расчетных сеток и научным вычислениям "NUM-GRID2008" (ВЦ РАН, Москва, 2008), на научно-исследовательских семинарах Вычислительного центра РАН, Московского физико-технического института, Института автоматизации проектирования РАН, а также на семинарах Департамента информационных технологий Республики Иещия (С-DAC, Пуна, Индия), Institut Français du Pétrol (Париж, Франция), ExxonMobil Upstream Research Co. (Хьюстон, США).

Личный вклад соискателя. В работах соискателя соавтором является научный руководитель Капорин И. Е., ему принадлежит постановка задач и выбор математических методов исследований. Соискателем проведены теоретические исследования, разработаны параллельные предобус-ловливания, выполнены параллельные реализации предложенных алгоритмов и исследованы границы применения предложенных методов.

Публикации. Основные публикации по теме диссертации включают И работ, из них 3 работы [6, 7, 10] в ведущих отечественных и международных рецензированных научных изданиях, рекомендованных ВАК, 3 в других международных рецензируемых изданиях, 1 в российском рецензируемом издании, 1 в трудах международной конференции, 3 публикации в других научных изданиях.

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

Структура и объем работы. Диссертационная работа состоит из введения, 5 глав, заключения, списка литературы, содержит 18 таблиц и 23 рисунка. Объем работы 138 страницы. Список литературы включает в себя 112 наименований.

КРАТКОЕ СОДЕРЖАНИЕ РАБОТЫ

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

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

На примере решения системы линейных алгебраических уравнений

Ас = 6 (1)

с симметричной положительно-определенной (с.п.о.) разреженной п х п матрицей А рассматривается предобусловлепный метод сопряженных градиентов (МСГ), в стандартной двучленной форме:

г0 = Ь~Ах0, р0 = Нг0, г = 0,1,..., гТНп

Р(АР1 '

хг+1 = х{ + ргаг, 74+1 = п ~ ар^щ, (2)

"Т+гНп+х

тТНп

А = -т ,,.. , Рг+1 = Нп+1 +

Здесь Н ~ А'1 - предобусловливатель для матрицы А, х* и г^ = Ь — Ах1 есть приближение к решению и повязка на г-й итерации, соответственно. Хорошо известная верхняя граница количества итераций, необходимых для уменьшения А_1-пормы начальной невязки в в раз, т.е. у/г[А~^Г{ < гугд Л_1го, выражается приближенной формулой

*с(е) =

\Vc\og-2 £

где

С = С(НА) = Атах(ЯЛ)/ЛтЬ(ЯА)

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

С целью нахождения более конструктивного критерия качества пре-добусловливапия была разработана альтернативная теория сходимости МСГ. И. Е.Капориным показано, что справедлива неулучшаемая оценка, измеряемая Я-нормой невязки, вида

^ (.К^ ~ 1) (3)

где

К = К (НА) = /¿еЬ(НА)

есть так называемое /Г-число обусловленности предобусловлеиной матрицы НА. Из неравенства (3) можно получить упрощенную, но наглядную верхнюю оценку количества итераций, необходимых для уменьшения Я-нормы начальной невязки в е-1 раз:

1к(е)= [к^ЯЧ-к^е"1] . (4)

В заключение формулируется основной вывод о возможности построить прсдобусловливатсль Я путем минимизации К(НА) (или соответственной верхней оценки) при ограничениях, накладываемых структурой предобусловливателя Н.

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

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

Далее для удобства предполагается, что матрица А симметрично от-маештабировапа таким образом, что па главной диагонали стоят единичные элементы. Приводится матричное соотношение, лежащее в основе пре-добусловливапия с использованием приближенного треугольного разложения второго порядка (Incomplete Cholesky 2-nd order, IC2), которое имеет вид

А = UTU + UTR + RTU, (5)

где U - невырожденная всрхпстрсугольная матрица (используемая как приближение правого множителя точного разложения Холецкого для матрицы А) и R - строго всрхпстрсугольная матрица ошибок с ''малыми" элементами.

Важными свойствами матрицы

М = JJ~T AU~l = 1 + RU~l + U~TRT,

полученной в результате двустороннего 1С2-предобусловливапия, являются следующие:

diagM = I, К{М) = det(7 + RA^R7), (6)

К(НА) = K{U~T AU~l) = К(М),

где Я = {UTUyl.

Рекуррентные соотношения для вычисления 1С2-разложепия могут быть легко получены из (5), например, если потребовать, чтобы в структуре разреженности матриц U и R не было совпадающих позиций ненулевых элементов, и чтобы их ненулевые элементы удовлетворяли соотношениям \uij\ ^ т и |г,-7'| < т соответственно. Здесь 0 < т -С 1 представляет собой малый параметр, определяющий качество разложения (будем называть также этот параметр параметром отсечения малых элементов в IC2 разложении). В частности, из второго уравнения (6) легко получается оценка

\ogK(M)^c0(A)r\ (7)

сопоставление которой с (4) объясняет, почему это разложение называется разложением 2-го порядка.

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

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

Приведено описание блочной версии алгоритма неполного обратного треугольного разложения (Block Incomplete Inverse Cholesky, ВИС). Пусть матрица А переупорядочена и разбита на блоки тем же образом, что и в случае применения блочио-диагонального предобусловливания, т.е. 2-й диагональный блок симметрично переупорядоченной матрицы имеет размер nt, где щ + ■■■ + rip = п. Здесь t = 1,2,... ,р и р - блочная размерность матрицы А. Для i-ro диагонального блока определим 'базисное" множество индексов как {A:(_i + 1,..., kt}, где kt~\ = щ + ■ ■ ■ + (ко = 0, кр = п), и введем "перекрывающиеся" множества индексов {jt{ 1),..., jt(mt - пг)}, jt(p) < kt-1- Для каждого t такое индексное множество, как правило, включает те индексы, по превосходящие kt, которые оказываются наиболее существенно связанными с i-м базисным множеством индексов, например, в смысле графа связности разреженной матрицы А (или матрицы, полученной из А заменой относительно малых элементов на нули). Здесь mt ^ щ и mi = щ, т.е. первое такое множество является пустым.

В этих обозначениях BIIC-предобусловливатсль может быть записан в следующем аддитивном виде:

где У( выбираются в виде прямоугольных матриц, столбцы которых являются единичными п-вскторами с определенными выше индексами,

vt = 1ел(1)I'"" \ej,(mt-nt)\eh^+i\ ■ ■ ■ Ы, t = 1,2,... ,р, (9)

а каждая всрхнстрсугольпая матрица Ut является правым множителем Хо-лецкого t-й "расширенной"' диагональной (mt х ?п,()-подматрицы V^AVt} т. с.

Далее рассматривается использование 1С2-разложснпя в ВПС-иредобусловливапии. Для каждого £ = 1,2.... ,р заменим точное разложение Холецкого (10) соответствующим приближенным 1С2-разложсписм, аналогичным (5):

где Rt ~ локальные матрицы ошибок с элементами < т. Обозначим

и будем для простоты предполагать, что стратегия отсечения малых элементов выбрана таким образом, что первые mt — щ диагональных элементов Ut совпадают с соответствующими элементами [/¡. Используя это предположение и изложенные выше свойства IC2- и ВНС-разложепий, доказывается следующая теорема.

Теорема 1. Пусть А - с.п.о. матрица, тогда из свойств разложений IC2 с порогом т отсечения малых элементов и ВИС с выделением р блоков следует оценка отношения К-чисел обусловленности матриц НА

VtTAVt = U?Uty t = 1,2,...,p.

(10)

VtTAVt = UjUt + UjRt + RTtUt.

(H)

и HA:

= ndet(Jmi + R^AVr'Rf) < e-2. (12)

К (HA) К (HA)

Используя (7) и (12), получаем

log к (ha) ^ log к (ha) + c0t2, и из (4) непосредственно следует оценка

Ч'М) ~ *я(е) + с0т2.

Результаты численных экспериментов подтверждают эту зависимость от т для набора из 5-и тестовых задач, подробно описанных далее. Они также показывают, что при т = 1СП3 количество итераций незначительно отличается от количества итераций для "точного" BIIC-предобусловливаиия (8), отвечающего значению т — 0, и остается практически неизменным при т < 1СГ4.

В заключение делается вывод, 'что для некоторого достаточно малого значения параметра отсечения т такое комбинированное BIIC-IC2-разложепие будет давать предобусловливатель Я, позволяющий достичь почти такой же скорости сходимости итераций МСГ, как и /^-оптимальное ВНС-предобусловливашге (8). В то же время ВПС-Ю2-разложсние требует значительно меньших затрат для вычисления, хранения и применения пре-добусловливателя по сравнению суточным" ВИС-предобусловливанисм (8).

Четвертая глава содержит описание параллельной реализации и методик балансировки вычислений. В начале главы дан обзор современных параллельных ЭВМ, моделей параллельных вычислений и способов оценки эффективности вычислений на таких ЭВМ. Для повышения параллельной эффективности применяется балансировка вычислений на различных этапах решения. Приведены различные варианты проведения балансировки на этапе выполнения итераций МСГ.

Далее описана параллельная реализация предобусловлениого МСГ с использованием аддитивного приближенного разложения с перекрытиями. Данный способ построения эффективных предобусловливаний был реализован с использованием обменов данными на основе стандарта MPI. Так была достигнута переносимость программного обеспечения на широкий класс ЭВМ с параллельной архитектурой.

Предположим, что система линейных уравнений (1) решается на параллельной ЭВМ с распределенной памятью, имеющей р процессорных элементов. Будем также предполагать, что количество блоков совпадает с количеством процессоров, на которые разбита исходная матрица A, a i-il блок соответствует t-му процессору, t = 1,2 ,...,р.

Предложенный алгоритм может быть реализован следующим образом. На этапе построения BIIC-IC2 предобусловливания выполняются IC2-разложения (11) локальных подматриц VfAVt на локальном t-м процессоре. При этом никаких обменов данными не происходит.

На этапе итераций МСГ (2) с выполнением обменов данными происходят следующие операции: а) умножение матрицы А на вектор, б) умножение прсдобусловливатсля Я на вектор, в) скалярные произведения.

Выполнение умножения матрицы А на вектор является часто используемой и хорошо исследованной операцией. Ее параллельная эффективность зависит от графа связности матрицы и топологии разбиения па подобласти. Эффективность умножения предобусловливателя Я па вектор зависит кроме этого еще и от геометрии перекрытия подобластей. Параллельная эффективность выполнения скалярных произведений определяется только количеством неизвестных в подобластях. Балансировки вычислений па всех описанных этапах взаимосвязаны, причем наиболее сложной является балансировка заполненности предобусловливателя. Основное внимание уделено именно этой проблеме.

Из вышеизложенного следуют свойства параллельной реализации метода ВПС-Ю2-МСГ:

1. Общее количество инициализаций обменов не зависит от количества блоков и размерности системы линейных уравнений и равно 5 па каждой итерации МСГ.

2. Коммуникационные затраты невелики по сравнению с арифметическими затратами.

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

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

Перекрытие подобластей строится с использованием структуры разреженности матрицы А9, где q - целое число, q ^ О, называемое в дальнейшем параметром перекрытия.

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

ки ого заполненности. Отдельно рассматриваются отсечения элементов в традиционном Ю-разложенни первого порядка и два варианта отсечения элементов в 1С2-разложепнп: по треугольным сомножителям л по матрица ошибок 1С2-разложспня.

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

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

Рассматривается постфильтрация множителей [/< по традиционно применяемой стратегии. Главным недостатком традиционной фильтрации является то, что при се применении сходимость итераций может ухудшиться. Это может произойти, даже если используется ставил тированное 1С2-разложенне, которое запишем в виде:

А = иТи + 1/тЛ + Нти - 5,

где Я - строго верхнетреугольпая матрица ошибок с элементами г^, |гу| < г, 5 - симметричная матрица ошибок 5 = Б7 > 0 с элементами вц, |ву| < т2, ! ^ |5г,| < Схт2 для некоторого С\ > 1, и т - параметр отсечения, 0 < т « 1. Известно, что для предобусловлсшюй матрицы М = и~тАи~1 при выполнении неравенства Н1 II ^ 76" для всех собственных значений матрицы М верпа оценка

А,-(Л/)

В случае, когда 7 является величиной порядка 1, можно расчитывать па хорошую устойчивость и быструю сходимость МСГ.

В связи с этим рассматривается расщепление матрицы {/, используемое в традиционной стратегии фильтрации:

и = и + й,

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

Теорема 2. Для всех собственных значений пре.Ообусловленной матрицы М = и~ТА0~1, при выполнении условий ^ 75 и

йтй ^ 7Б, где 7 < (^/7 + у/1 + гу)~2 верна оценка

А,(М) < г = 1 ,...,п. (13)

1 - 7 - 2ут7

Таким образом, порчены условия, показывающие, что при использовании традиционной фильтрации сходимость МСГ может ухудшиться. Результаты приведенных иижс численных экспериментов подтверждают, что примерно в 1/4 случаев происходит резкое ухудшение сходимости более чем в 2 раза, вызванного именно возрастанием Атах(М).

Рассматривается постфильтрация, выполняемая по принципиально новой стратегии. Для приближенного треугольного разложения второго порядка Ю2 предложено переносить наибольшие по величине элементы из матрицы ошибок в вычисленные треугольные множители. Это не только обеспечивает сбалансированность вычислений, но и но ухудшает скорость сходимости МСГ.

Для этого Ю2-разложеиие переписывается в эквивалентном виде:

А + + Б = (и + Я)т(и + Д),

которое позволяет рассмотреть расщепление матрицы ошибок Ю2-разложения

д = д + д,

где Д включает некоторое количество относительно больших по величине элементов матрицы Д, и оценить качество предобусловливания, полученного с использованием восполнения треугольного сомножителя 1С2-разложения

и = и + А.

Теорема 3. Для предобусловленной .матрицы М = U~TAU~l при выполнении условия RTR ^ S + RTR верна оценка

К{М) ^ (det U)2/ det А.

Таким образом, если норма матрицы R достаточно мала, то К-число обусловленности иредобусловлсшюй матрицы М имеет верхнюю границу, которая не зависит от возмущения, вносимого описанной постфильтрацисй. Это дает основания полагать, что сходимость МСГ при этом, по крайней мерс, не ухудшится.

Результаты приведенных численных экспериментов показывают, что для рассмотренных задач более чем в 70% случаев было зафиксировано уменьшение времени счета, обусловленное улучшением сходимости МСГ, а в остальных случаях результаты оказались сравнимы.

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

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

Для сравнения предложенного метода с некоторыми другими стратегиями иредобусловливания, также требующими малого количества инициализаций обменов на этапе итераций, приводятся данные численных экспериментов для точечного метода Якоби (Point Jacobi, PJ), приближенного блочного метода Якоби (Block Jacobi, BJ), а также обобщения метода Якоби, использующее взвешенное перекрытие блоков (Overlapping Block Jacobi, OBJ).

Необходимо заметить, что в последнем случае перекрытие подобластей (или блоков) строится "во все стороны" (OBJ-тии перекрытия), а не только в сторону меньших индексов, как это происходит в случае BIIC-предобусловливания. Второй особенностью метода OBJ является то, что до и после обращения диагональных подматриц применяется масштабирование весовой диагональной матрицей где j-ii элемент D равен числу взаимных перекрытий подобластей для j-го узла. Построение OBJ-предобусловливателя может быть представлено в виде

Ною = 1Г1/2 Vt(VtTAVt)-lVtT^j D-V2, (14)

где Vt определены выше в (9).

Заметим, что блочно-диагональный предобусловливатсль, эквивалентный применению блочного метода Якоби без перекрытия (где mt =

nt, t — 1,2,... ,р), может быть записан кал

#bj = Et(E?AEt)'1 Е?, (15)

где Et = [efct_1+1| • ■ • |ekt], t = 1,2,... ,p, так что EfAEt в точности является 1-й диагональным блоком матрицы А. Рассмотренные методы BJ и OBJ называют также методами типа декомпозиции области (Domain Decomposition, DD).

Предобусловливанию по точечному методу Якоби (PJ) соответствует матрица-предобусловливатель

#pj = Щ1,

где D'aJ1 - диагональ матрицы А. Если же исходная матрица А уже отмас-штабироваиа по Якоби, то НР] = I.

Во всех рассматриваемых метода.х (кроме точечного метода Якоби) использовался одинаковый алгоритм приближенного 1С2-разложения для диагональных блоков.

Приведены результаты численных экспериментов для задач из коллекции университета Флориды1. Численные эксперименты выполнялись на высокопроизводительном вычислительном комплексе MVS6000IM, использующем процессоры Dual Core Intel Itanium 2 с частотой 1.6 ГГц.

1http://ww.<^.ufl.edu/tesearch/sparse/matrices

В качество основной модельной задачи (обозначаемой далее DCD) была рассмотрена конечно-разностная дискретизация задачи минимизации квадратичного функционала min J с|Дк|2с1а;^Х2 на квадрате [0,1] х [0,1], с кусочно-постоянными коэффициентами в виде

и начальным приближением щ = 0. Было выбрано значение с* = 103 и использовалось 300 точек дискретизации вдоль каждого из направлений.

Также были рассмотрены несколько задач с плохо обусловленными разреженными матрицами из коллекции, собранной в университете Флориды. Характеристики тестовых матриц представлены в табл.1, где использованы следующие обозначения: п - размерность матрицы А: ^(Л) - количество ненулевых элементов матрицы Л; СопсЦЛэ) - нижняя оценка спектрального числа обусловленности матрицы Лэ = Б^2 АО^2, т.е. исходной матрицы, симметрично отмасштабированпой к единичной главной диагонали. Следует заметить, что только одна из рассмотренных матриц является М-матрицей. Необходимо также обратить внимание на широкий круг областей потенциального применения рассматриваемого метода решения систем линейных уравнений.

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

Во всех экспериментах правая часть выбиралась равной Ъ = Ах*, используя х* = [1 1 • • ■ 1]т в качестве точного решения и хо = [0 0 • • • 0]т в качестве начального приближения. Итерации МСГ проводились до достижения критерия остановки

В табл.2 использованы следующие обозначения: р - количество процессоров (блоков); Dens - плотность предобусловливатсля NZ({7)/NZ(i/J4)

1, иначе,

с граничными условиями и = 0и ди/дп = 0, с точным решением и* = u(xi, х2) = х\еТлХ2 Sm(7TXi) sin(7TX2)

ЫЮ1М1, £ = 10"8.

(16)

Таблица 1. Характеристики матриц из коллекции университета Флориды.

Матрица n NZ(A) Cond(As) Описание

DCD 90000 1164004 9.61xl010 расчет конструкций (пластина с утолщением)

bcsstk25 msc23052 gridgena 15439 23052 48962 252241 1154814 512084 3.55x106 7.57xl07 1.35x105 расчет конструкций (небоскреб) расчет конструкций оптимизация вычислительных сеток

evxbqpl 50000 349 963 1.31x10s выпуклое квадратичное программирование

oilpaii s3dkq4m2 s3dkt3m2 xl04 73 752 90119 90449 108384 3597188 1820891 3753461 10167624 1.03x10s 1.74X1010 3.12xl010 2.12xl09 расчет конструкций МКЭ для цилиндрической оболочки МКЭ для цилиндрической оболочки расчет конструкций (SESAM)

bmw7st_l g2_circuit pvvtk hood 141347 150102 217918 220542 7339667 726674 11634424 10768436 1.03xl07 2.34x10s 4.25xl07 5.42xl05 растет конструкций (матрица жесткости) моделирование СБИС (AMD, М-матрица) расчет конструкций (воздуховод) расчет конструкций (автомобилестроение)

BenElechil 245874 msdoor 415863 af_l_kl01 503625 parabolic_fcm 525 825 13150496 20240935 17550675 ■3674625 1.84x10s 1.95xl08 4.43x10" 2.10xl05 пространственно-плоскостная задача расчет конструкций (автомобилестроение) расчет конструкций (листовая штамповка) МКЭ для конвекции-диффузии

по отношению к плотности верхнего треугольника исходной матрицы А; Iter - количество итераций, необходимых для достижения требуемой точности s (16); Тил - общее время решения задачи в секундах. Приведенные здесь данные свидетельствуют, что рассматриваемая задача DCD является достаточно сложной для решения. Это видно как из нижней оценки ее числа обусловленности (см. табл. 1), так и из приведенных данных о сходимости традиционного применяемого предобусловливания первого порядка 1С(т) (эквивалентного использованию 1С2(т, т)), т = Ю-3, Ю-4, Ю-5,10~6, а также (ходимости применяемого предобусловливания второго порядка 1С2(т, т2). Количество итераций для 1С(т), при использовании разумных значений т, достаточно велико. Результаты, приведенные в табл. 2, показывают значительно бблыную эффективность предобусловливания второго порядка точности IC2, сочетающего скорость сходимости 1С(т2) с заполненностью 1С(т).

В табл. 2 также представлены результаты сравнения стратегий предобусловливания: 1С2-разложение с описанным выше ВИС предобусловлива-нисм и с предобусловливанием по блочному методу Якоби (BJ) (который эквивалентен BIIC с параметром q = 0). Можно видеть, что сходимость

Таблица 2. Сравнение предобусловлнваний для задачи DCD (г = 10 3).

Метод Р Dens Iter Ttat

10(10" -3) 1 2.44 9013 513.80

10(10" "4) 1 21.34 5014 1273.97

1С(10" -Ь) 1 45.65 2113 1058.28

1С(10- -ь) 1 51.79 С01 377.20

Ю2(т, Т*) 1 4.47 62G 94.38

1С2(т, т2)-ВНС 8 5.37 528 8.17

1С2(т, t2)-BJ 8 4.02 7033 55.13

для BJ достаточно медленная по сравнению с предобусловлнваписм ВИС, где применяются специальным образом построенные перекрытия.

Еще более наглядно сравнение с другими методами предобусловли-ваиия представлено на рис. 1а. Приведены данные времени счета для предобусловлнваний по методам PJ, BJ, OBJ, а также используемого метода ВИС. Для аппроксимации задач в подобластях во всех случаях использовалось разложение второго порядка IC2. Видно, что время решения для точечного метода Якоби PJ неприемлемо велико во всех случаях. При использовании предобуеловливания OBJ обычно происходит несущественное изменение количества итераций по сравнению с предобусловлнваписм без перекрытия BJ, но из-за значительно более плотного предо-бусловливателя выигрыша во времени решения не происходит. Результаты проведенных расчетов показывают, что для предложенного BIIC-IC2-предобусловливапия получаются наилучшие результаты как по количеству итераций, так по общему времени решения.

Будем в дальнейшем использовать следующие обозначения: STD -описанное выше ВНС-1С2-прсдобуеловливапие, без использования какой-либо постфильтрации и балансировки; MIN - традиционно применяемая постфильтрация каждой матрицы Üt, которая уменьшает все NZ(f/i) до их минимального значения среди всех р блоков; МАХ - предложенный способ дополнения в каждую Ut максимальных по модулю элементов из матрицы ошибок lit. чтобы увеличить все NZ(f/t) до их максимального значения.

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

Time

gSTD

(в) (г)

Рис. 1. Задачи из коллекции университета Флориды: (а) сравнение различных методов для р = 8; (б) ускорение без балансировки для ВНС-1С2-МСГ; (в) ускорение для балансировки по традиционной, стратегии MIN; (г) ускорение для балансировки по предложенной стратегии МАХ.

1 Ttot(p)'

где Ttot(p) - общее время решения задачи на р процессорах, р = 1,2,4,8.

Данные, приведенные па рис. 16, показывают, что предложенное BIIC-IC2 предобусловливание даже без применения какой-либо балансировки показывает достаточно устойчивые высокие результаты ускорения. Если исключить из рассмотрения задачи небольшой размерности, то достигнутое ускорение составит более 5 на 8 процессорах.

На рис. 1в представлены результаты, показывающие, что наиболее очевидная и традиционно применяемая балансировка, основанная на стратегии MIN, даст медленную скорость сходимости. Так, резкое возрастание количества итераций для задачи DCD для BIIC-IC2 прсдобусловливаиия с 8 процессорами и традиционно используемой стратегией балансировки MIN объясняется чрезмерно большим значением Amax(Ai), где М = НА - пре-добусловлсппая матрица (см. (13)). В этом случае оцененная граница для этой величины получилась равной 101.55, в то время как типичное значение, которое гарантирует численную устойчивость итераций МСГ, ниже 2. По этой же причине для задачи ''1118023052'' сходимость стала хуже почти в 100 раз. Здесь соответствующие спектральные характеристики оценены как Ащах(М) > 5 ■ 104 и Cond(iTi) > 6 • 10®. Примерно в половине рассмотренных случаев количество итераций при применении стратегии MIN возросло более чем в 2 раза, а в 17% случаев более чем в 10 раз. Наблюдаемые отрицательные эффекты могут быть связаны с нарушением сформулированного в теореме 2 условия па верхнюю границу нормы отбрасываемой части Ut-

С другой стороны, рис. 1г показывает, что предложенная балансировка на основе стратегии МАХ достаточно работоспособна и эффективна. Улучшение достигнуто более чем в 70% случаев, в остальных случаях время счета отличалось не более чем па 5%, что можно считать лежащим в области ошибок измерения.

Отметим, что применение стратегий STD и МАХ позволило получить свсрхлинсйпос ускорение в 14% случаев. Причем для стратегии! МАХ выигрыш был более значителен. Помимо возможных эффектов, зависящих от оборудования (улучшение работы кэш-памяти при уменьшении размерности подзадачи на каждом из процессоров) указанное явление может иметь

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

Проводилось также сравнение с опубликованными результатами других исследователей, которыми рассматривался алгебраический многосеточный метод (AMG) и алгебраический рекурсивный многоуровневый метод (ARMS) на примере решения 5 задач из рассмотренной коллекции. Исходя из опубликованных ими результатов, сравнения проводились для решения на одном процессоре. Практически во всех случаях предложенный метод дал лучшие результаты (иногда до 100 раз по времени счета), что свидетельствует о его более широкой области применимости при решении жестких задач.

Следует также заметить, что исследования не были ограничены только описанной коллекцией матриц. Были рассмотрены и другие области применения предложенного метода, а именно, задачи упругости тонкостенных оболочек и задачи линейной теории упругости. Кроме этого исследовалась эффективность предложенного метода при использовании до 64 блоков в ВИС-1С2-предобусловливании. Результаты этих расчетов, также приведенные в диссертационной работе, еще раз продемонстрировали надежность и эффективность предложенной методики построения параллельных предо-бусловливапий.

Выводы

Как теоретические исследования, так и экспериментальные результаты параллельной реализации ВПС-1С2-предобусловливания для МСГ показали его высокую эффективность для решения систем линейных алгебраических уравнений с разреженными с.п.о. плохо обусловленными матрицами. Результаты проведенных численных экспериментов свидетельствуют, что чем "атожнее" для решения выбирается задача, чем больше число обусловленности матрицы системы линейных уравнений, тем больше оказывается эффективность предложенного метода предобусловливаиия по сравнению с традиционно используемыми параллельными методами.

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

размерности берется решаемая задача, тем большего ускорения удастся достигнуть. Отсутствие замедления даже для задач малой размерности показывает высокую эффективность использования ресурсов параллельной ЭВМ. Повышение эффективности использования каждого процессора было достигнуто за счет повой методики балансировки вычислений. Таким образом, параллельные свойства предложенного алгоритма ВНС-Ю2-МСГ при использовании умеренного количества процессоров оказались достаточно хорошими, как и ожидалось из теоретических построений.

Основные результаты

1. Разработан, исследован и реализован новый параллельный метод ВПС-1С2-МСГ решения жестких систем линейных уравнений с разреженными с.п.о. матрицами.

2. Разработана, исследована и реализована новая методика балансировки вычислений для подзадач на различных процессорах.

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

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

Основное содержание диссертационной работы изложено в следующих публикациях:

[1] Коньшин H.H. О реализации методов разбиения области на вычислительных системах с матричными процессорами /'/ Методы математической физики. - М.: ОВМ АН СССР. 1988. 23 с.

[2] Коньшин H.H. Оптимизация многосеточных методов декомпозиции области //. Числсипыс методы и программное обеспечение. - М.: ОВМ АН СССР. 1990. 22 с.

[3] Ibragimov I.V., Kolesnikov P.A., Konshin I.N., et al. Numerical experiments with industrial iterative solvers on massively parallel computers. I: Numerical experiments with the A_SPARSE solver on CRAY T3D /,/

In: Proc. of the High Performance Computing Symposium'95. Phoenix, Arizona, USA. 1995. P. 283-289.

[4] Kaporin I.E.. Konshin I.N. Parallel solution of large sparse SPD linear systems based on overlapping domain decomposition //' Lecture Notes in Computer Science, Vol.1662. Springer-Verlag, Berlin-Heidelberg-New-York. 1999. P. 436-445.

[5] Axclsson 0., Kaporin I., Konshin I., et al. Comparison of algebraic solution methods on a set of benchmark problems in linear elasticity // Tech. Report of Department of Mathematics, University of Nijmegen, The Netherlands. 2000. 89 p.

[6] Капорин И.Е.. Коныиин И.Н. Параллельное решение симметричных положительно-определенных систем на основе перекрывающегося разбиения на блоки //' Ж. Вычисл. Матем. и Матем. Физ. 2001. Т. 41. №4. С. 515-528.

[7] Kaporin I.E., Konshin I.N. A parallel block overlap preconditioning with inexact submatrix inversion for linear elasticity problems // J. Numer. Lin. Alg. Appl. 2002. V.9. №2. P. 141-162.

[8] Капорин И.Е., Коньшин И.Н. Параллельное решение линейных систем при использовании приближенной факторизации перекрывающихся блоков // В кн.: "Математическое моделирование. Проблемы и результаты." (Ред. О.М.Белоцерковский, В.А.Гущин) - М.: Наука. 2003. С. 308-319.

[9] Kaporin I.E., Konshin I.N. Post-filtering of 1С factors for load balancing in parallel preconditioned CG solvers // Numerical geometry, grid generation and high performance computing (NUMGRID2008). (Eds. V.A.Garanzha, Yu.G.Evtushenko, B.K.Soni, N.P.Weatherill) - M.: Folium, 2008. P. 158164.

[10] Капорин И.E., Коньшин И.Н. Постфильтрация множителей Юг-разложения для балансировки параллельного предобусловливания // Ж. Вычисл. Матем. и Матем. Физ. 2009. Т. 49. №6. С. 940-957.

[11] Kaporin I.E.. Konshin I.N. Load balancing of parallel block overlapped incomplete Cholesky preconditioning // Lecture Notes in Computer Science. Vol. 5698. Springer-Verlag, Berlin-Heidelberg. 2009. P. 304-315.

Заказ № 305/ 11 /09 Подписано а печать 02.11.2009 Тираж 100 экз. Усл. п. л. 2

ООО "Дизайнбюро" 119333, г. Москва, ул. Вавилова, д. 46 (495)649-93-76, (495)518-73-78

e-mail: designbr@list.ru designbr@yandex.ru www.designbr.ru

Содержание диссертации автор исследовательской работы: кандидата физико-математических наук, Коньшин, Игорь Николаевич

Введение

1 Метод сопряженных градиентов (МСГ)

1.1 Прямые и итерационные методы решения линейных систем

1.2 Описание предобусловленного метода сопряженных градиентов

1.3 Оценка сходимости МСГ через спектральное число обусловленности

1.4 Оценка сходимости МСГ через К-число обусловленности

1.5 Устойчивость предобусловливаний МСГ.

2 Методы приближенных треугольных разложений

2.1 Предобусловливания, основанные на треугольном разложении.

2.2 Неполные и приближенные треугольные разложения

2.3 Предварительное масштабирование как этап предобусловливания.

2.4 Теория приближенного треугольного разложения

2-го порядка.

2.4.1 Приближенное треугольное разложение 2-го порядка

2.4.2 Улучшение обусловленности, достигаемое применением приближенных треугольных разложений.

2.4.3 Устойчивость приближенных треугольных разложений

2.5 Алгоритмы безотказного приближенного треугольного разложения.

2.6 Трудности распараллеливания IC2 разложения.

3 Параллелизуемое аддитивное предобусловливание

3.1 Методы построения параллельных предобусловливаний

3.1.1 Использование окаймленной блочно-диагональной структуры.

3.1.2 Использование приближенных обратных матриц

3.2 Блочное неполное обратное треугольное разложение.

3.2.1 Построение BIIC-предобусловливания.

3.2.2 Оценка качества предобусловливания по методу BIIC-IC2.

3.3 Диагональное и блочно-диагональное предобусловливание

4 Параллельная реализация и балансировка вычислений

4.1 Параллельные ЭВМ и параллельные вычисления

4.2 Параллельная реализация итерационных методов.

4.3 Описание параллельной реализации.

4.4 Способы балансировки вычислений

4.5 Теоретический анализ стратегий постфильтрации построенного предобусловливателя.

4.6 Описание реализации балансировки и применения параллельного ВПС-1С2-предобусловливания.

5 Численные эксперименты

5.1 Тестовые задачи и методика проведения численных экспериментов

5.2 Численные эксперименты для задач упругости тонкостенных оболочек.

5.3 Численные эксперименты для задач теории линейной упругости в механике упругого тела.

5.4 Численные эксперименты по балансировке вычислений для задач из коллекции университета Флориды.

5.5 Сравнение метода ВНС-1С2-МСГ с другими методами

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

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Объект исследования и актуальность темы. Актуальность темы диссертационной работы обусловлена необходимостью эффективного численного решения систем линейных алгебраических уравнений с разреженными симметричными положительно-определенными матрицами большой размерности. Численное решение систем линейных уравнений является одной из наиболее часто встречающихся задач в научно-технических исследованиях и практических приложениях. Подобные задачи возникают, например, в математической физике при численном решении дифференциальных и интегральных уравнений. При этом прикладные задачи часто требуют решения систем линейных уравнений с большим числом неизвестных. К таким системам, например, приводит численное решение двумерных и особенно трехмерных задач математической физики, в которых условия физической и геометрической аппроксимации двумерной и трехмерной области диктуют использование достаточно мелкой расчетной сетки с большим числом расчетных узлов. Для возникающих систем линейных уравнений применение прямых методов решения оказывается неприемлемым как по причине большого времени решения, так и недостаточности объема оперативной памяти для хранения данных задачи. Итерационные методы решения систем линейных уравнений намного экономичнее с точки зрения использования оперативной памяти, но для ускорения их сходимости требуется строить эффективные предобусловливания, особенно при решении систем с плохо обусловленными матрицами.

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

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

В соответствии с целью исследования были поставлены следующие задачи:

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

2. Разработка эффективного способа балансировки вычислений для подзадач на различных процессорах.

3. Исследование эффективности метода в зависимости от параметров предобусловливания и количества процессоров.

Научная новизна:

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

Положения, выносимые на защиту:

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

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

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

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

Апробация работы. Результаты работы докладывались и обсуждались на международной конференции "РаСТ99" (Санкт-Петербург, Россия, 1999 г.), Всемирном 16-ом Конгрессе "IMACS 2000" по научным вычислениям, прикладной математике и моделированию (Лозанна, Швейцария, 2000), Голландско-российском симпозиуме NWO (МГУ, Москва, июнь 2000), Симпозиуме NWO (Амстердам-Ниймеген, Голландия, ноябрь 2000), Втором Международном научно-практическом Семинаре и Всероссийской молодежной школе "Высокопроизводительные параллельные вычисления на кластерных системах" (Нижний Новгород, Россия, 2002), Международной конференции "Parallel CFD 2003" (Москва, 2003), Международной конференции "VIII Забабахинские научные чтения" (РФЯЦ-ВНИИТФ, Сне-жинск, Россия, 2005), Международной конференции по численной геометрии, построении расчетных сеток и научным вычислениям "NUM-GRID2008" (ВЦ РАН, Москва, 2008), на научно-исследовательских семинарах Вычислительного центра РАН, Московского физико-технического института, Института автоматизации проектирования РАН, а также на семинарах Департамента информационных технологий Республики Индия (С-DAC, Пуна, Индия), Institut Frangais du Petrol (Париж, Франция), ExxonMobil Upstream Research Co. (Хьюстон, США).

Личный вклад соискателя. В работах соискателя соавтором является научный руководитель Капорин И. Е., ему принадлежит постановка задач и выбор математических методов исследований. Соискателем проведены теоретические исследования, разработаны параллельные предобусловливания, выполнены параллельные реализации предложенных алгоритмов и исследованы границы применения предложенных методов.

Публикации. Основные публикации по теме диссертации включают 11 работ, из них 3 работы [27, 28, 37] в ведущих отечественных и международных рецензированных научных изданиях, рекомендованных ВАК, 3 в других международных рецензируемых изданиях, 1 в российском рецензируемом издании, 1 в трудах международной конференции, 3 публикации в других научных изданиях.

Структура и объем работы. Диссертационная работа состоит из введения, 5 глав, заключения, списка литературы, содержит 18 таблиц и 23 рисунка. Объем работы 138 страниц. Список литературы включает в себя 112 наименований.

Заключение диссертации по теме "Вычислительная математика"

Основные результаты

1. Разработан, исследован и реализован новый параллельный метод ВНС-1С2-МСГ решения жестких линейных систем с разреженными с.п.о. матрицами.

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

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

Заключение

Список источников диссертации и автореферата по математике, кандидата физико-математических наук, Коньшин, Игорь Николаевич, Москва

1. Коиъшин И.Н. Особенности решения линейных систем на вычислительных системах с матричными процессорами // 5-я Московская конференция по проблемам кибернетики и вычислений. Научный совет "Кибернетика" АН СССР, Статистический отдел "Москва". Москва. 1986.

2. Коньшин И.Н. Блочные методы исключения на параллельных ЭВМ для ленточных СЛАУ большой размерности // Организация вычислений на супер-ЭВМ. МФТИ. Москва. 1987. Деп. ВИНИТИ N.2184-B87, 7с.

3. Коньшин И.Н., Малясов С.А. О применении метода разбиения области в одной задаче квантовой механики // Сопряженные уравнения и алгоритмы возмущения. ОВМ АН СССР. Москва. 1988. С.92-98.

4. Коньшин И.Н. О реализации методов разбиения области на вычислительных системах с матричными процессорами j j Методы математической физики. ОВМ АН СССР. Москва. 1988. 23с.

5. Коньшин И.Н., Еремин А.Ю. Проблемы отображения многоуровневых алгоритмов разбиения области на архитектуру параллельных ЭВМ // 2-й Семинар соц. стран "Вычислительная механика и автоматизация проектирования". Москва-Ташкент. 1988.

6. Коньшин И.Н. Оптимизация многосеточных методов декомпозиции области // Численные методы и программное обеспечение. ОВМ АН СССР. Москва. 1990. 22с.

7. Gushchin V.A., Kononov I.N., Konshin I.N., Konshin V.N. Parallel numerical modeling of clean room air flows // The third Russian-Japan joint symposium on CFD. Book of abstracts. Part II. August 25-30, 1992. Vladivostok. Russia. P.121-122.

8. Кононов И.Н., Коньшин B.H., Коньшин И.Н. Численный анализ потоков воздуха в чистом производственном помещении // Сб. докл. 2-й конференции ассоциации инженеров по контролю микро-загрязнений. 12-16 октября, 1992. АСИНКОМ-92. Суздаль. 1992. С.136-140.

9. Kaporin I.E., Konshin I.N. Parallel solution of SPD linear systems based on algebraic DD // In: Proc. of European Conference on Parallel Computing, Euro-Par 2000. 29.8.-1.9.2000, Munich, Germany. 2000. 6p.

10. Kaporin I., Konshin I. Accuracy and efficiency of parallel implementation of matrix multiplication. // In: Proc. of 6th International IMACS Conference on Applications of Computer Algebra. St. Petersburg, Russia. June 25-28, 2000. P.ll.

11. Konshin I.E., Kaporin I.E. Parallel solution of STW linear systems based on BIIC2 preconditioning // Dutch-Russian NWO Workshop. Moscow State University, Moscow. June 23-24, 2000.

12. Konshin I., Kaporin I. Numerical testing of parallel solvers for a set of benchmark problems // NWO Workshop 2000. 16-17 November 2000, CWI, Amsterdam. 21 November 2000, KUN, Nijmegen.

13. Kaporin I.E., Konshin IN. Benchmark problems in linear elasticity: parallel solution // Tech. Report No.0028. Department of Mathematics, University of Nijmegen, The Netherlands, December 2000. 22p.

14. Капорин И.Е., Коньшин И.Н.: Параллельное решение симметричных положительно-определенных систем на основе перекрывающегося разбиения на блоки // Ж. Вычисл. Матем. и Матем. Физ. 2001. Т.41, N.4, Р.515-528.

15. Kaporin I.E., Konshin I.N. A parallel block overlap preconditioning with inexact submatrix inversion for linear elasticity problems //J. Numer. Lin. Alg. Appl. 2002. Vol.9. No.2. P.141-162.

16. Garanzha V.A., Kaporin I.E., Konshin I.N. Truncated Newton type solver with application to grid untangling problem //J. Numer. Lin. Alg. Appl. 2004. Vol.11. No.5-6. P.525-533.

17. Kaporin I.E., Konshin I.N. Recursive scaling, permutation and 2x2-block splitting in ILU preconditionings // In: Proc. 2-nd International Conference on Matrix Methods and Operator Equations. July 23-27, 2007. Moscow. P.33-34.

18. Капорин И.Е., Коньшин И.Н. Постфильтрация множителей IC2-разложения для балансировки параллельного предобусловливания // Ж. Вычисл. Матем. и Матем. Физ. 2009. Т.49. N.6. С.940-957.

19. Вогачев К.Ю. Основы параллельного программирования М.: Бином. 2003.

20. Еремин А.Ю., Капорин И.Е. Частичное исключение как способ переобусловливания // Препр. №216. М.: Отд. Вычисл. Матем. АН СССР. 1988. 11 с.

21. Капорин И.Е. О предобусловливании и распараллеливании метода сопряженных градиентов //В (Ортега Дж.) Введение в параллельныеи векторные методы решения линейных систем. М.: Мир. 1990. С.343-355.

22. Капорин И.Е. О предобусловливании метода сопряженных градиентов при решении дискретных аналогов дифференциальных задач // Дифферент ур-ния. 1990. Т.26. N.7. С.1225-1236.

23. Капорин И.Е. Двухуровневые явные предобусловливания метода сопряженных градиентов // Дифференц. ур-ния. 1992. Т.28. N.2. С.329-339.

24. I. Е. Kaporin. Explicitly preconditioned conjugate gradient method for the solution of nonsymmetric linear systems // Int. J. Computer Math. 1992. Vol.40. P.169-187.

25. Капорин И.Е. Оценки границ спектра двусторонних явных предобус-ловливаний // Вестн. МГУ. Выц.15. Вычисл. матем. и кибернетика 1993. Т.2. С.28-42.

26. Милюкова О.Ю. Некоторые параллельные итерационные методы с факторизованными матрицами предобусловливания для решения эллиптических уравнений на треугольных сетках // Ж. Вычисл. Матем. и Матем. Физ. 2006. Т.46. N.6. С.1096-1113.

27. Ортега Дою. // Введение в параллельные и векторные методы решения линейных систем. М.: Мир. 1990.

28. Посыпкин М.А., Хританков А.С. О понятии производительности в распределенных вычислительных системах // Проблемы вычислений в распределенной среде. М.: Труды Института Системного Анализа РАН. 2008. Т.32. С.26-32.

29. Хейгеман Л., Янг Д. // Итерационные методы. М.: Мир. 1990.

30. Ajiz М.А., Jennings A. A robust incomplete Choleski-conjugate gradient algorithm // Int. J. Numer. Methods Engrg. 1984. Vol.20. P.949-966.

31. Axelsson 0. A class of iterative methods for finite element equations // Computer Meth. Appl. Mech. Engrg. 1976. Vol.9. No. 2. P.123-137.

32. Axelsson 0., Lindskog G. On the rate of convergence of the preconditioned gradient methods // Numer. Math. 1986. Vol.48. No. 5. P.499-523.

33. Axelsson 0. // Iterative Solution Methods. Cambridge University Press, Cambridge, New York. 1994.

34. Axelsson O., Kaporin I., Kucherov A., Neytcheva M. Benchmark problems in linear elasticity. Part I: Direct and robust second order incomplete factorization methods, (submitted to Int. J. Numer. Methods Engrg.)

35. Benzi M., Kouhia R., Тита M. An assessment of some preconditioning techniques in shell problems // Techn. Rept LA-UR-97-3892. Los Alamos Nat. Lab., Los Alamos, NM. 1997.

36. Benzi M. Preconditioning Techniques for Large Linear Systems: A Survey. // Journal of Computational Physics. 2002. Vol.182. No.2. P.418-477.

37. Benzi M., Тита T. A robust incomplete factorization preconditioner for positive definite matrices // Numer. Linear Algebra Appl. 2003. V.10. P.385-400.

38. Campos F.F., Birkett N.R.C. An efficient solver for multi-right-hand-side linear systems based on the CCCG(?7) method with applications to implicit time-dependent partial differential equations // SIAM J. Sci. Com-put. 1998. Vol.19. P.126-138.

39. Chow E., Saad Y. Approximate inverse preconditioners via sparse-sparse iterations. // SIAM Journal on Scientific Computing. 1998. Vol.19. No. 3. P. 995-1023.

40. Chow E. A priori sparsity patterns for parallel sparse approximate inverse preconditioners // SIAM J. on Scientific Computing. 2000. Vol.21. No.5. P. 1804-1822

41. Chow E. Parallel implementation and practical use of sparse approximate inverse preconditioners with a priori sparsity patterns // Int. Journal of High Performance Computing Applications. 2001. Vol.15. No.l. P.56-74.

42. Concus P., Golub GO'Leary D. Generalized conjugate gradient method for the numerical solution of elliptic partial differential equations // In: Sparse matrix computations. (Eds. J.Bunch, D.Rose) 1976. P.309-332.

43. Cuthill E., McKee J. Reducing the bandwidth of sparse symmetric matrices. // In: Proceedings of the 24th National Conference of the ACM. 1969. P.152-172.

44. Davis T.A. // University of Florida Sparse Matrix Collection. 2009 http: / / www.cise.ufi.edu/research / sparse / matrices

45. Dennis J.E., Wolkowicz H. Sizing and least change secant methods // SIAM J. Numer. Anal. 1993. Vol.10. P.1291-1314.

46. Dubrulle A. A. Retooling the method of block conjugate gradients // Electr. Transact. Numer. Anal. 2001. Vol.12. P.216-233.

47. Duff I.S., Reid J.K. Performance evaluation of codes for sparse matrix problems //In: Performance Evaluation of Numerical Software (Ed. L.D.Fosdick). Elsevier, North-Holland, New York. 1979. P.121-135.

48. Duff I.S., Grimes R.G., Lewis J.G. Sparse matrix test problems // ACM Trans. Math. Softw. 1989. Vol.15. No.l. P.l-14.

49. Duff I.S., Grimes R.G., Lewis J.G. The Rutherford-Boeing sparse matrix collection // Tech. Report TR/PA/97/36. CERFACS, Toulouse, France.

50. Tech. Report RAL-TR-97-031. Rutherford Appleton Lab. Tech. Report ISSTECH-97-017. Boeing Information and Support Service. 1997.

51. Duff I.S., Riyavong S., Van Gijzen B. Parallel preconditioners based on partitioning sparse matrices // Report RAL-TR-2004-040. Сотр. Sci. and Engrg. Dept., Atlas Centre, Rutherford Appleton Lab., Oxon, 2004. 19p.

52. Eisenstat S.C., Gursky M.C., Schultz M.H., Sherman A.H. Yale sparse matrix package I: The symmetric codes // Int. J. Numer. Methods Engrg. 1982. Vol.18. P.1145-1151.

53. Feng Y.Y., Owen D.R.J., Peric D. A block conjugate gradient method applied to linear systems with multiple right-hand sides // Comput. Meth. Appl. Mech. Engrg. 1995. Vol.127. P.203-215.

54. Graham, E., Forsyth P.A. Preconditioning methods for very ill-conditioned three-dimensional linear elasticity problems // Int. J. Numer. Meths. Engrg. 1999. Vol. 44. P.77-99.

55. Hestenes M.R., Stiefel E. Methods of conjugate gradients for solving linear systems // J. Research Nat. Bur. Standards. 1952. Vol.49. P.409-436.

56. Hypre, High Performance Preconditioners: Users Manual. // Technical Report. Lawrence Livermore National Laboratory. 2006.

57. Jennings A., Malik G.M. Partial elimination // J. Inst. Math. Appl. 1977. Vol.20. P.307-316.

58. Jimack P.K. Domain decomposition preconditioning for parallel PDE software // Engineering computational technology. Civil-Comp press. Edinburgh. UK. 2002.

59. Jones M. Т., Plassman P.E. An improved incomplete Cholesky factorization // ACM Trans. Math. Software. 1995. Vol.21. P.5-17.

60. И.Е. Капорин. О предобусловливаиии метода сопряженных градиентов при решении дискретных аналогов дифференциальных задач // Дифферент ур-ния 1990. Vol.26. No. 7. Р.1225-1236.

61. Kaporin I. Explicitly preconditioned conjugate gradient method for the solution of unsymmetric linear systems // Internat. J. Comput. Math. 1992. Vol.40. P.169-187.

62. Kaporin I.E. New convergence results and preconditioning strategies for the conjugate gradient method // Numer. Linear Algebra and Appl. 1994. Vol.1. No.,2. P.179-210.

63. Kaporin I.E. High quality preconditioning of a general symmetric positive definite matrix based on its UTU + UTR + RT/-decomposition. // Numer. Linear Algebra and Appl. 1998. Vol.5. No.5. P.483-509.

64. Kaporin I.E. Optimizing the UTU + UTR + Дг /-decomposition based Conjugate Gradient preconditionings // Rep. 0030. November 2000. Department of Mathematics, Catholic University of Nijmegen, Nijmegen, The Netherlands. 16p.

65. Karypis G., Kumar V. Multilevel k-way hypergraph partitioning // Tech. Report 98-036. Dept. Сотр. Sci. Engrg., Army HPC Research Center, Univ. of Minnesota, MN. 1998.

66. Kershaw D.S. The incomplete Cholesky-conjugate gradient method for the iterative solution of systems of linear equations //J. Comput. Physics. 1978. Vol.26. P.43-65.

67. Kolotilina L.Yu., Yeremin A.Yu. Factorized sparse approximate inverse preconditionings, I: theory // SIAM Journal on Matrix Analysis and Applications. 1993. Vol.14. No.l. P.45-58.

68. Lipton R.J., Rose D.J., Tarjan R.E. Generalized nested dissection // SIAM J. Numer. Anal. 1979. Vol.16. P.346-358.

69. Maclahlan S., Saad Y. A greedy strategy for coarse-grid selection // SIAM J. Sci. Comput. 2007. V. 29. No. 5. P. 1825-1853.

70. Manteuffel T.A. An incomplete factorization technique for positive definite linear systems // Math. Comput. 1980. Vol.34. P.473-497.

71. Marshall A.W., Olkin I. /j Inequalities: Theory of Majorization and its Applications. Academic Press, New York. 1979.

72. MatrixMarket collection of sparse matrices. 2009 http://math.nist.gov/MatrixMarket.

73. Message Passing Interface (МР1'). http://www.mpi-forum.org/docs/docs.html

74. MPICH. http://info.msc.anl.gov/pub/mpi

75. Notay Y. On the convergence rate of the conjugate gradients in presence of rounding errors // Numer. Math. 1993. Vol.65. P.301-317.

76. O'Leary D.P. The block conjugate gradient algorithm and related methods // Linear Algebra and Appls. 1980. Vol.29. P.293-322.99. http://www.parallel.ru

77. Rao V.N., Kumar V. Parallel depth-first search, part II: Analysis j I Int. J. Par all. Programming. 1987. Vol.16, No.6. P.501-519.

78. Ruge J. W., Stiiben K. Algebraic multigrid (AMG) // Multigrid Methods (Ed. McCormick S. F.) Frontiers Appl. Math. Vol. 3. Philadelphia: SIAM. 1987. P. 73-130.

79. Saad Y. // Iterative Methods for Sparse Linear Systems. Society for Industrial and Applied Mathematics. Philadelphia. PA. 2003.

80. Saad Y., Suchomel B. ARMS: An algebraic recursive multilevel solver for general sparse linear systems // Numer. Linear Algebra Appl. 2002. Vol.9. P.359-378.

81. Suarjana M., Law K.H. A robust incomplete factorization based on value and space constraints // Int. J. Numer. Methods Engrg. 1995. Vol.38, P.1703-1719.

82. The TOP500 Supercomputer Sites. http://www.top500.org, 2009.

83. Tismenetsky M. A new preconditioning technique for solving large sparse linear systems // Linear Algebra and Appl. 1991. Vol.154-156. P.331-353.

84. Tuff A. D., Jennings A. An iterative method for large systems of linear structural equations // Int. J. Numer. Methods Engrg. 1973. Vol.7. P.175-183.

85. Ugar В., Aykanat C. Partitioning sparse matrices for parallel preconditioned iterative methods // SIAM J. Sci. Comput. 2007. Vol.29. No.4. P.1683-1709.

86. Ugar В., Aykanat C. Revisiting hypergraph models for sparse matrix partitioning // SIAM Review. 2007. Vol.49. No.4. P.595-603.

87. Vastenhouw В., Bisseling R.H. A two-dimensional data distribution method for parallel sparse matrix-vector multiplication // SIAM Review. 2005. Vol.47. No.l. P.67-95.

88. Wilkinson J.H., Reinsch C. // Handbook for automatic computation. Linear algebra, Part II. Springer-Verlag, Berlin. 1971.