Предобусловливание итерационных методов решения систем линейных алгебраических уравнений

Капорин, Игорь Евгеньевич

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

Капорин, Игорь Евгеньевич АВТОР

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

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

2011 ГОД ЗАЩИТЫ

01.01.07 КОД ВАК РФ

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

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

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

КАПОРИН Игорь Евгеньевич

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

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

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

1 з ОНТ 2011

Москва — 2011

4856837

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

Официальные оппоненты: доктор физико-математических наук,

профессор

Ильин Валерий Павлович доктор физико-математических наук, Жуков Виктор Тимофеевич доктор физико-математических наук, Нечепуренко Юрий Михайлович

Ведущая организация: Казанский (Приволжский) федеральный

университет

Защита состоится ок-т&^гШ г. в ч. 00_ мин

на заседании Диссертационного совета Д 002.045.01 при Учреждении Российской академии наук Институте вычислительной математики РАН по адресу: 119333, г. Москва, ул. Губкина, 8.

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

Автореферат разослан СЛпТА^А 2011 г.

Учёный секретарь Диссертационного совета доктор физико-математических наук Бочаров Г.А.

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

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

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

(а) разработка Хестенсом и Штифелем в 1952 году1 метода сопряженных градиентов, который, в свою очередь, тесно связан с алгоритмом тридиагонализации, опубликованным Ланцошем в 1950 году;2

(б) разработка техники предобусловливания для итерационных методов (которая эквивалентна предварительному линейному преобразованию матрицы системы), сформулированной в связи с ускорением метода сопряженных градиентов в работах Даниэля 1967 года,3 Г. И. Марчука и Ю. А. Кузнецова 1968 года,4 О. Аксельсона5 и др.

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

гМ. R. Hestenes and Е. Stiefel. Methods of conjugate gradients for solving linear systems, J. Research Nat. Bur. Standards, 49 (1952) 409-436.

2C. Lanczos. An iteration method for the solution of the eigenvalue problem of linear differential integral operators. J. Research Nat. Bur. Standards, 45 (1950) 255-282.

3J. W. Daniel. The conjugate gradient method for linear and nonlinear operator equations. SIAM J. Nu-mer. Analysis, 4 (1967) 10-26.

4Г. И. Марчук, Ю. А. Кузнецов. Об оптимальных итерационных процессах. Докл. Акад. Наук СССР, т.181, №6 (1968) 1331-1334.

50. Axelsson. A class of iterative methods for finite element equations. Computer Meth. Appl. Mech. Engrg., 9 (1976) 123-137.

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

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

Отметим первые попытки разработки нестандартных подходов к

предобусловливанию метода сопряжённых градиентов, см., напр., [1, 2],

а также недавний обзор,7 которые опирались на критерий' качества

предобусловливания, выбранный в виде евклидова расстояния ||/„ — M\\f 2

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

бО. Axelsson and G. Lindskog. On the rate of convergence of the preconditioned conjugate gradient method, Numer. Math. 48, (1986) 499-523.

7T. Huckle. Factorized Sparse Approximate Inverses for Preconditioning. The Journal of Supercomputing v.25, no. 2 (2003) 109-1 IT.

во многих важных ситуациях, см. напр., обсуждение в работе [26]. В последней работе содержится законченная теория сходимости метода обобщенных минимальных невязок в терминах величины ||/„ — М\\р и дополнительного параметра, характеризующего отделение спектра М от нуля. Однако, в контексте разреженных матриц большого размера, не представляется реалистичной постановка задачи о минимизации ||/„ — М\\р при условии, например, ограниченности ||А/-1||.

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

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

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

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

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

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

зрения нового числа обусловленности.

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

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

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

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

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

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

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

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

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

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

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

s0. Axelsson. Iterative solution methods. Cambridge University Press, Cambridge, 1994.

как с точки зрения оптимизации как спектрального числа обусловленности, так и К-числа обусловленности.

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

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

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

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

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

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

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

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

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

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

3. Теория предобусловливания симметричных положительно определенных матриц посредством обратных приближенных

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

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

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

Обоснованность и достоверность результатов.

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

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

Россия, 2002), международной конференции "Parallel CFD 2003" (Москва, 2003), Международной конференция "VIII Забабахинские научные чтения" (РФЯЦ-ВНИИТФ, Снежинск, Россия, 2005), международной конференции по вычислительной геометрии, генерации сеток и научным вычислениям "NUMGRID2008" (ВЦ РАН, Москва, 2008), международной конференции по вычислительной геометрии, генерации сеток и высокопроизводительным вычислениям "NUMGRID2010" (ВЦ РАН, Москва, 2010), международной конференции по матричным методам в математике и приложениях "МММА-2011" (ИВМ РАН, Москва, 2011), научно-исследовательских семинарах Вычислительного центра РАН и Института вычислительной математики РАН, рабочих семинарах ExxonMobil Upstream Research Со. (Хьюстон, США).

Публикации. По теме диссертации опубликовано 28 работ, из них 5 в отечественных рецензированных изданиях, рекомендованных ВАК, 8 в международных рецензируемых изданиях из рекомендованного ВАК списка "Web of Science: Science Citation Index Expanded" (база по естественным наукам), 2 в других международных рецензируемых изданиях, 1 в российском рецензируемом издании, 4 в трудах всероссийских конференций, 5 в трудах международных конференций, 3 публикации в других научных изданиях, (список работ приводится в конце автореферата).

В совместных работах [15, 19, 22, 24, 25, 27, 28] И.Н.Коныпину принадлежит общая схема параллельной реализации метода и описание численных экспериментов, а автору диссертации - разработка и исследование предобусловливания.

Благодарности. Автор выражает благодарность А. Ю. Еремину, привлекшему автора к работе по теме диссертации, J1. Ю. Колотилиной за внимание к работе и ценные обсуждения, профессору О. Аксельсону за интерес к работе и всестороннюю ее поддержку на протяжении многих лет, И. Н. Коньшину за длительное сотрудничество в области численного тестирования и параллельной реализации разработанных автором методов, В. А. Гаранже за полезные обсуждения, замечания и поддержку работы, В. Ю. Правильникову за сотрудничество в области

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

Структура и объем работы. Диссертационная работа состоит из введения, 6 глав, заключения и списка литературы. Объем работы 216 страниц. Список литературы включает 192 наименования.

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

Первая глава содержит изложение теории нового числа обусловленности

далее называемого К-числом обусловленности. Указанный матричный функционал определяется для вещественной симметризуемой матрицы М (т.е., существует невырожденная матрица Т такая, что М = ТБТ-1, где 5 - симметричная положительно определенная матрица). Проводится детальное сопоставление К-числа обусловленности со стандартным спектральным числом обусловленности

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

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

(1)

с = С(М) = Атах(М)/Ага;п(М).

(2)

положительно определенной (с.п.о.) разреженной п X гс-матрицей М. Рассматривается стандартный метод сопряженных градиентов (м.с.г.):

г0 = Ь-М:с0, ро = г0; ¿ = 0,1,...: гТг-

°' ~ рТМр•' х<+х = Х{ +1на" г,+1 = г'' ~ (3)

Р% = ^Г-> Р«+1 = г«+1 + РФ*-

г- г,-

Здесь ж,- и г,- — Ъ — Мх{ есть приближение к решение и невязка на г-й итерации, соответственно. Хорошо известная оценка убывания М-1-нормы невязки, где

1И1 = у/гтМ~хг. имеет вид

2 / ((ус(М) + Л(уЩм)-1\ \ЫЫ-г - / \^ус(м)+ 1у )' [4)

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

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

Новая теория сходимости м.с.г., использующая предложенное в работах [3, 4] К-число обусловленности (1), основывается на верхней границе евклидовой нормы невязки ||г,|| = \/ггг вида

Ы<(К(М).'<(5)

Окончательный вид этой оценки был установлен в работах [8, 11], где была также доказана ее неулучшаемость в том же смысле, в котором это справедливо для стандартной оценки (4).

Вторая глава продолжает изложение теории метода сопряженных градиентов (м.с.г.), для решения с.л.а.у. Ах — Ь. а также практики его применения. Особое внимание уделяется предобусловленному варианту м.с.г., позволяющему значительно повысить его вычислительную

эффективность. Если с.п.о. матрица Н ~ А-1 выбрана в качестве предобусловливателя для матрицы А, то этот метод описывается следующими формулами:

г0 = Ъ-Ах о, ро = Нг0\ г' = 0,1,.

rjHn

(6)

Важным новым результатом, содержащимся в этой главе, является (упрощенная) оценка числа итераций метода сопряженных градиентов, достаточного для уменьшения Я-нормы начальной невязки в раз через К-число обусловленности

которую можно соотнести со стандартной оценкой числа итераций

Основными недостатками критерия качества предобусловливания, основанного на оценке (8) и связанного с минимизацией спектрального числа обусловленности С {НА), является как проблематичность его непосредственного практического использования в сколько-нибудь общем случае (например, при построении предобусловливаний для произвольной с.п.о. матрицы), так и его недостаточная согласованность с реальной скоростью сходимости м.с.г., особенно в условиях действия вычислительной погрешности9. Поэтому в заключение главы формулируется основной тезис диссертации, выносимый на защиту, а именно, вывод о целесообразности построения предобусловливателя Я исходя из условия минимизации K-числа обусловленности К (НА) (или его верхней оценки) при ограничениях, накладываемых структурой предобусловливателя Я.

9Y. Notay. On the convergence rate of the conjugate gradients in presence of rounding errors. Numer. Math, v.65 (1993) 301-317.

*к(е) = ["log2K(#A) + log2(e-1)] .

(7)

ic(e)= -\/С(ЯА) log (2e-1) .

(8)

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

Теорема 1 [12]. Пусть А = АТ > 0 и Ю1ад{А) ■= 1„, и пусть при некотором 0 < 5 -С 1 выполнены условия структурной устойчивости

Ъже(иТи - А) < п5, йе^ити) > с^ А. Тогда из условия отсечения

кад > Л з>

следует оценка заполнения

пг{и) < п I 1 +----I ,

являющаяся асимптотически неулучшаемой при малых г, 6 и больших п,К(А).

Анализ последних достижений в области построения приближенных треугольных разложений, прежде всего алгоритма Тисменецкого,10 и его более практичной приближенной версии,11 привел к следующей матричной формулировке [12], которая позволила не только дать теоретическое обоснование соответствующему предобусловливанию, но и сформулировать понятие приближенного треугольного разложения второго порядка:

А = иТи + итЕ1 + - Е2 (9)

10M. Tismenetsky. A new preconditioning technique for solving large sparse linear systems. Linear Algebra Appls. v.154-156 (1991) 331-353.

UM. Suarjana and K. H. Law. A robust incomplete factorization based on value and space constraints. Int. J. Numer. Methods Engrg. v.38 (1995) 1703-1719.

Здесь U - невырожденная верхняя треугольная матрица (используемая как приближение правого множителя Uq точного треугольного разложения А = U{fUo для матрицы А), Е\ - строго верхняя треугольная матрица погрешности с относительно малыми по модулю элементами, и Еч - симметричная неотрицательно определенная матрица погрешности с элементами, существенно меньшими по сравнению с Е\.

Принципиальным преимуществом такого разложения (9) (при выборе Е\ = О(т), Е'2 = 0(т2)) является тот факт, что при той же нижней границе на модули ненулевых элементов U (а тогда и сопоставимой верхней границе ее заполнения nz(?7), см. выше Теорему 1), точность приближения оказывается на порядок лучше в смысле зависимости от числа обусловленности матрицы А:

С(U-TAU~l) = 1 + 0[ту/сЩ + r2C(A)), logЩи^Аи-1) < с2{А)т\

тогда как для схемы типа Дженнингса-Малика,12 см. также13 (где Ei =0, E'i = О(г)) имеем гораздо более слабую оценку

C{U-TAU~l) = 1 + О(гОД), logK(и~тАи~1) < а(А)т.

Рекуррентные соотношения для вычисления приближенного треугольного разложения могут быть получены непосредственно из (9), например, если потребовать, чтобы в структуре разреженности матриц U и Ei не было совпадающих позиций ненулевых элементов и чтобы значения этих ненулевых элементов удовлетворяли соотношениям |(С0у| > т и |(Z?i),-j| < г соответственно. Здесь 0 < г <С 1 представляет собой малый параметр, определяющий качество разложения (будем также называть г параметром отсечения малых элементов в приближенном треугольном разложении).

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

12А. Jennings and G. М. Malik. Partial elimination. J. Inst. Math. Appl., 20, 307-316, 1977.

13M. A. Ajiz and A. Jennings. A robust incomplete Choleski-conjugate gradient algorithm. Int. J. Numer. Methods Engrg., 20, 949-966,1984.

Теорема 2 [12]. Пусть в тождестве (9) матрицы удовлетворяют следующим условиям: V - верхняя треугольная с положительной диагональю, Е\ - строго верхняя треугольная,

А = АТ> О, Е2 = Е% > О, Е^Ех <

ЦЕгА-'Е^ < а, \\А-^Е2А~^\\ < /?, Тогда справедливы оценки

АтЬ(СГ1Л£Гт) > (уЪ + + а + ~2, Хтии~1Аи-т) <1 + 7-С{и~1Аи-т) < (1 + 7) (у^ + +а + (3)2-

Результаты, полученные в этом случае при оценке К-числа обусловленности, дает

Теорема 3 [12]. Пусть в тождестве

А = иТи + итЕ1 + Е'[и - Е2,

матрицы удовлетворяют следующим условиям: II - верхняя треугольная с положительной диагональю, Е\ - строго верхняя треугольная,

А = АТ> О, Е2 = > О, Тогда справедливы оценки

К(и~1Аи~т) < Ле^Л + Е(Е! + Е^/А^А,

К(1/~1Аи~т) < пип ^^гасе(£1г£1 + Е2) + £ >

где 1>а(Ь) представляет собой кусочно-постоянную неубывающую функцию, равную количеству собственных значений матрицы А, не превосходящих

Заметим, однако, что как вычисление (приближенная факторизация), так и применение (решение систем с матрицами

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

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

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

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

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

GAGT = 1 + Е,

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

ненулевых элементов нижнего треугольника матрицы Aq, g = 1,2,3,____

Значения определяются из условия минимума K(GAGT). Для этого достаточно выполнить п треугольных разложений подматриц А, построенных на индексных множествах ненулевых элементов каждой строки G. Такое предобусловливание будем обозначать IIC(q) (англ. Incomplete Inverse Cholesky).

Более совершенный алгоритм предусматривает построение

усеченной структуры разреженности {(г, _;')} согласно критерию ^ где 0 < т << 1, с последующим перевычислением

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

Однако, несмотря на "идеальную" параллелизуемость алгоритмов н.о.т.-разложения, получаемый предобусловленный м.с.г. может иногда оказаться лишь ненамного более эффективным, чем в случае простейшего диагонального предобусловливания Я = (01ад(Л))-1.

Проблема в том, что для указанных "очевидных" реализаций неполного обратного треугольного разложения при увеличении количества ненулевых элементов матрицы 6? улучшение свойств матрицы (ЗАС?7 происходит гораздо медленнее, чем при сопоставимой степени заполнения в обычном приближенном треугольном разложении (и тем более в разложении 2-го порядка), см. [9]. Кроме того, слишком быстро растут затраты на вычисление и применение предобусловливателя (3, см. далее Фиг.4.

Более эффективный подход описан в работах [4, 11], где был предложен блочный вариант метода неполного обратного треугольного разложения, позволивший в значительной мере преодолеть упомянутые недостатки простейшего н.о.т.-разложения.

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

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

исходной матрицы Ая, ^ = 1,2,3,____Тогда, в силу совпадения блочно-

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

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

И, наконец, упомянутые факторизации блоков вовсе не обязаны быть точными, вместо них можно использовать их приближенные треугольные факторизации 2-го порядка с порогом отсечения т, как это было предложено и реализовано в работах [15,19, 22, 24, 25, 27, 28]. Такое предобусловливание будем обозначать BIIC(p,q)-IC2(r) (Block Incomplete Inverse Cholesky - Incomplete Cholesky of the 2nd order).

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

Этот факт был проверен экспериментально в [15, 19] и обоснован теоретически в работе [22]. Дальнейшее развитие этот метод получил в работах [27, 28], где на основе специальных свойств треугольного разложения 2-го порядка был построен, обоснован и испытан метод улучшения балансировки вычислений по блокам такого предобусловливателя.

Приведем более детальное описание блочной версии неполного обратного треугольного разложения (BIIC) [4, 9, 11, 22]. Пусть матрица А переупорядочена и разбита на блоки тем же образом, что и в для известного блочно-диагонального предобусловливания (Block Jacobi, BJ), т.е., i-й диагональный блок симметрично переупорядоченной

матрицы имеет размер щ, где щ Л-----\-пр = п. Здесь f = 1,2,...,р и

р - блочная размерность матрицы А. Для t-го диагонального блока определим "базисное" множество индексов как +1,... ,kt}, где

kt-i = щ -)-----hnt-i (ко = 0, кр = п), и введем "перекрывающиеся"

множества индексов {jt(l),... ,ji(m{ — n*)}, jt{p)<kt-v Для каждого t такое индексное множество, как правило, включает те индексы, не превосходящие kt, которые оказываются связанными с i-м базисным множеством индексов, например, относительно структуры разреженности матрицы Ая (или q-й степени матрицы, полученной из

А заменой относительно малых элементов на нули). Здесь тп* > тц и гп\ = щ, т.е. первое такое множество является пустым.

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

игтУ?, (10)

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

^ = [ел(1)1'"1ел(ш«-п|)К.1+1|"-К], ¿ = 1,2,...,р, (11)

а каждая верхнетреугольная матрица С/« является правым множителем Холецкого ¿-й "расширенной" диагональной (ш4 х т()-подматрицы У?АУХ, т.е.

ЦтАЦ = и?ии ¿=1,2,...,р. (12)

В работах [4, 9, 10] показано, что ВПС-предобусловливатель Я обладает свойством /^-оптимальности в следующем смысле. Обозначим через С множество разреженных нижнетреугольных матриц, которые имеют ненулевые элементы только в позициях (г,^), где

€ 0'((1), • • ■ — + 1,..., г}, к^ + 1 < г < ки (13)

при 1 = 1,... ,р. Тогда для построенного предобусловливателя справедливо представление Н = где

С = а^ттК (<ЗЖ?Г)

(напомним, что имеет место тождество К(С?А£?Т) = К(НА)).

Далее рассматривается приближенный вариант ВНС-предобусловливания, получаемый заменой для каждого £ = 1,2,... ,р точного ит[/-разложения (12) на соответствующее приближенное 1С2-разложение, аналогичное (9):

= + \jjRt + ЩЩ. (14)

Я = ¿WГ

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

я-Ew1

о о О /„,

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

Теорема 4 [22]. Пусть А является с.п.о. матрицей; тогда из свойств разложений 1С2(т) и ВПС(р) вытекает следующая оценка отношения К-чисел обусловленности матриц НА и НА:

= П ** С- + ВД^Г'яГ) < ехр(Сог2),

где 0 < со = со(А).

Таким образом, из (7) следует оценка числа итераций

(15)

*к(е) =

log2K(#A) + log2 (е~1) +

-и , со(Л1

(16)

log 2

т. е., для некоторого достаточно малого значения параметра отсечения г в 1С2-разложениях блоков построенный ВИС-1С2-предобусловливатель Н имеет почти такую же оценку числа итераций м.с.г., что и К-оптимальное ВПС-предобусловливание (10). В то же время, BIIC-IC2-разложение требует значительно меньших затрат для вычисления, хранения и применения предобусловливателя по сравнению с "точным" BIIC-предобусловливанием (10), (см. [15]).

Численное тестирование описанных выше методов предобусловливания IIC(q), IIC(q;r), BIIC(p;q)-IC2(r) проводилось с использованием вещественной симметричной положительно определенной матрицы apacheB, взятой из коллекции разреженных матриц университета Флориды14. Эта матрица, возникающая в

14Т. А. Davis, Y .F .Hu Y. University of Florida sparse matrix collection. To appear in: ACM Trans, on Math. Software, 2011. V. 38; http://www.cise.ufl.edu/research/sparse/matrices

некоторой трехмерной конечно-разностной задаче, имеет достаточно большой размер и довольно плохо обусловлена: п = 715176, nz(A) = 4817870, Лтах(А5) = 2, Amax(As)/Amin(^) = 1.2 • 106, где As = (Diag^))~1/2A(Diag(yl))~1/'2 представляет собой матрицу коэффициентов, симметрично масштабированную к единичной диагонали.

Итерации предобусловленного м.с.г. начинались с нулевого начального приближения xq = 0 и критерием их остановки служило условие ||i — Ахк\\н < 10_12||b||#, где правая часть b = Ах, вычислялась посредством умножения матрицы А на вектор тестового решения х,(г) = 1, i = 1,..., п.

Параметры предобусловливаний задавались следующим образом:

(а) для IIC(q) и IIC(q;r)) использовалась перестановка исходной матрицы коэффициентов, уменьшающая среднюю ширину ленты матрицы; параметр отсечения выбирался равным г = 0.01;

(б) для BIIC(p;q)-IC2(r)) использовалась перестановка исходной матрицы коэффициентов, обеспечивающая достаточно заполненную блочно-диагональную часть с почти совпадающими размерами блоков; для этого использовалась подходящая версия алгоритма PABLO15; значения параметров в методе приближенного треугольного разложения IC2 [12] выбирались равными г = 0.01, г2 = 0.0001, а = 0, 9 = 0.1;

Расчеты проводились на настольном ПК AMD Athlon 64x2 Dual Core Processor 4200+ (2.21 GHz, память 2 Gbytes) под ОС MS Windows XP v.2002 Service Pack 2 с компилятором Intel Fortran Compiler 6.0 for Windows и опциями компиляции ifl /Ох /G6 /Qxi /Qip *.for. Время счета приводится в секундах.

Фиг. 1 и 2 иллюстрируют результаты, полученные для числа блоков р = 1, 2, 4, ... , 256, 512 при использовании предобусловливаний BJ(p)-IC2(0.01) и ВНС(р;4)-1С2(0.01). Заметим, что для обоих методов случаю р = 1 отвечает один и тот же метод IC2(0.01)-CG, описанный в предыдущей главе.

15D. Fritzsche, A. Frommer, and D. В. Szyld, Extensions of certain graph-based algorithms for preconditioning. SIAM Journal on Scientific Computing, v.29, no.5, (2007) 2144-2161.

Фиг.1. Число итераций м.с.г. в зависимости от числа блоков для предобусловливаний BJ(p)-IC2(0.01) и BIIC(p;4)-IC2(0.01) и задачи apacheZ.

ЖЗЕЕЯЗЗЕЖЗЯСЗЭИ

200

175

150

„'И

5.00 я

50 25

— apachs2. BJ-IC2 - apactw2. BHC-IC2

Nblocks

ю5

Фиг.2. Общее время счета м.с.г. в зависимости от числа блоков р для предобусловливаний ВЛ(р)-1С2(0.01) и ВНС(р;4)-Ю2(0.01) и задачи араскеЗ.

Фиг.З. Число итераций м.с.г. в зависимости от степени перекрытия д для методов ПС(д), ПС(д;0.01) и ВНС(512;?)-1С2(0.01) и задачи араМ.

Фиг.4. Общее время счета м.с.г. в зависимости от степени перекрытия д для методов ПОД, ПС(д;0.01) и ВПС(512;д)-1С2(0.01) и задачи арасЬе2.

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

Замечание. Наблюдаемый эффект независимости числа итераций м.с.г. при числе блоков до р = 512 объясняется тем, что возмущение, вносимое разделением на блоки, не только сравнительно невелико (так как является, в смысле К-оптимальности, минимально возможным), но и остается незаметным на фоне погрешности, вносимой приближенной факторизацией блоков. При этом с ростом р размеры блоков уменьшаются, их факторизация становится точнее, и при достаточно большом р уже начинается рост числа итераций за счет выявления возмущений, вносимых разделением на блоки. Если факторизацию блоков осуществлять точно, то, несомненно, рост числа итераций м.с.г., будет наблюдаться уже начиная с малых значений р.

Фиг. 3 и 4 иллюстрируют результаты, полученные для разных значений степени перекрытия д от 0 до 6 при использовании предобусловливаний НС(я), ИС(ч;0.01) и ВНС(512;я)-1С2(0.01) соответственно. Заметим, что при д = О методы НС(0) и НС(0;0.01) сводятся к точечному предобусловливанию Якоби, а метод ВПС(512;0)-1С2(0.01) переходит в ВЛ(512)-1С2(0.01). При возрастании д происходит существенное сокращение числа итераций, но также и удорожание затрат на предобусловливание.

Согласно показанной на Фиг.З. зависимости числа итераций от степени перекрытия, при равных значениях q блочное предобусловливание требует существенно меньше итераций, чем точечные. С другой стороны, усовершенствованный алгоритм ПС(я;0.01) требует почти столько же или даже меньше итераций по сравнению с простейшим методом ПС^) (который, как правило, использует

избыточное количество ненулевых элементов для формирования матрицы (3).

Зависимость времени счета от степени перекрытия качественно примерно одинакова для всех методов, см. Фиг.4. Здесь во всех случаях наблюдается выраженный оптимум по д (для методов ПС это 5 = 2 и для ВИС это д = 4). Однако количественно можно видеть заметное преимущество метода ПС^;0.01) над ПС(ц) (в полтора раза быстрее для 5 = 3 за счет существенного уменьшения заполненности предобусловливателя). В свою очередь, метод ВНС(512;3)-1С2(0.01) уменьшает время счета еще почти в два раза (по сравнению с ПС(3;0.01)).

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

Пятая глава представляет весьма наглядную иллюстрацию целесообразности применения К-оптимальных предобусловливаний метода сопряженных градиентов. В ней проведен теоретический анализ методов полиномиального предобусловливания Н=-р,1~1{А) симметричных положительно определенных матриц А общего вида.

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

где Тл(х) = ((х + у/х2 — 1)^ + (х — \/х2 — 1)^)/2 - многочлен Чебышева первого рода й-й степени,

0< (А) < Д,

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

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

Н = 1 + СБСТ.

Здесь С - фиксированная п х т-матрица полного столбцового ранга, причем т « п, а Я ■ симметричная тп х т-матрица, определяемая из условия минимума К (ЯЛ). Указанная конструкция была предложена в [4] и изучена в [8], где было установлено, что К-оптимизация такого предобусловливания приводит к выбору

5 = -(С^С)-1 + о{СтАС)~\ (17)

при некотором значении <т = о(А,С). Для практического использования можно положить а = 1.

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

Простая подстановка непосредственно приводит к обобщению полученной формулы в виде

Н= В'1+ У (~(УТВУУХ + (УТАУ)~1) Ут,

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

Заключение содержит следующие основные результаты диссертации.

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

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

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

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

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

Основные публикации автора по теме диссертации:

[1] Еремин А.Ю., Капорин И.Е. Спектральная оптимизация явных итерационных методов. I // В кн.: Численные методы и вопросы организации вычислений, вып. 7, Зап. научн. семин. ЛОМИ АН СССР, 1984, т.139, С.51-60.

(перевод: A. Yu. Eremin and I. E. Kaporin, Spectral optimization of explicit iterative methods. I. Journal of Mathematical Sciences, 1987, V.36, no.2, P.207-214)

[2] Еремин А.Ю., Капорин И.Е., Марьяшкин Н.Я. Решение линейных систем большой размерности на векторных и конвейерных суперЭВМ: новый подход к оптимизации явных итерационных методов //

В кн.: Вопросы кибернетики. Конвейеризация вычислений, (под ред. О.М.Белоцерковского и В.В.Щенникова), Научный совет АН СССР по комплексной проблеме "Кибернетика", Москва, 1986, С.114-124.

[3] Капорин И.Е. Альтернативный подход к оценке числа итераций метода сопряженных градиентов //В кн.: Численные методы и программное обеспечение (под ред. Ю. А. Кузнецова) - М.: ОВМ АН СССР, 1990, С.55-72.

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

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

[6] Капорин И.Е. Итерационное решение систем линейных уравнений с использованием неполной обратной треугольной факторизации //В кн: Прямые и обратные задачи математической физики. М.: Изд-во МГУ, 1991. С. 71-77.

(перевод: I. Е. Kaporin. Iterative solution of systems of linear equations using incomplete inverse triangular factorization. Computational Mathematics and Modeling vol.4, no.l (1993) P.28-32)

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

[8] Kaporin, I.E. Explicitly preconditioned conjugate gradient method for the solution of nonsymmetric linear systems // Int. J. Computer Math., 1992, V.40, P.169-187.

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

[10] Капорин И.Е. Оптимизация алгоритмов метода сопряженных градиентов. // В кн: Математические модели и оптимизация вычислительных алгоритмов. М.: Изд-во МГУ, 1994. С. 50-62. (перевод: I. Е. Kaporin, Optimization of conjugate gradient algorithms, Computational Mathematics and Modeling 1994, V.5, no.2, P.139-147)

[11] Kaporin, I.E. New convergence results and preconditioning strategies for the conjugate gradient method // Numer. Linear Algebra with Appls., V.l, no.2, 1994, P.179-210.

[12] Kaporin, I.E. High quality preconditioning of a general symmetric positive matrix based on its UTU+UTR+^¿/-decomposition // Numerical Linear Algebra Appl., 1998, V.5, P.484-509.

[13] Kaporin, I.E. Second order incomplete Cholesky preconditionings for the CG method // in: Proc. of IMMB'98, Iterative solution methods for the elasticity equations as arising in Mechanics and Biomechanics, Nijmengen, The Netherlands, Sept. 28-30, 1998, P.47-49.

[14] Еремин А.Ю., Капорин И.Е. Влияние наибольших собственных значений на численную сходимость метода сопряженных градиентов, Численные методы и вопросы организации вычислений. XIII, Зап. научн. сем. ПОМИ, 248, ПОМИ, СПб., 1998, С.5-16. (перевод: A. Yu. Yeremin and I. E. Kaporin. The influence of isolated largest eigenvalues on the numerical convergence of the CG method. Journal of Mathematical Sciences 2000, V.101, no.4, P.3231-3236)

[15] Kaporin,I.E. and Konshin,I.N. Parallel solution of large sparse SPD linear systems based on overlapping domain decomposition // in: Parallel Computing Technologies (Ed. V.Malyshkin), Proceedings of the 5th International Parallel Computing Technologies Conference (PaCT-99), St.-Petersburg, Russia, September 6-10, 1999. Lecture Notes in Computer Science, Vol. 1662, Springer-Verlag, Berlin - Heidelberg - New-York, 1999, P.436-445.

[16] 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.

[17] Axelsson O., Kaporin I., Konshin I., Kucherov A., Neytcheva M., Polman В., Yeremin A. 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, 89p.

[18] Axelsson O., Kaporin I. On the sublinear and superlinear rate of convergence of conjugate gradient methods // Numerical Algorithms, 2000, V.25, P.l-22.

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

[20] Axellson О., Kaporin I. Optimizing two-level preconditionings for the conjugate gradient method // Rep.0116, August 2001, Department of Mathematics, Catholic University of Nijmegen, Nijmegen, The Netherlands, 21p.

[21] Kaporin I.E. Using the Modified 2nd Order 1С Decomposition as the Conjugate Gradient Preconditioning // In: Proc. of PRISM Conference, 20-23 May. 2001, Nijmegen, The Netherlands.

[22] Kaporin I.E., Konshin I.N. A parallel block overlap preconditioning with inexact submatrix inversion for linear elasticity problems // Numerical Linear Algebra with Applications, 2002, V.9, no.2, P.141-162.

[23] Kaporin I.E. Using the Modified 2nd Order Incomplete Cholesky Decomposition as the Conjugate Gradient Preconditioning // Numerical Linear Algebra with Applications, 2002, V.9, P.401-408.

[24] Kaporin I., Konshin I.N. Parallel Conjugate Gradient Preconditioning via Incomplete Cholesky of overlapping Submatrices // In: Book of Abstracts, Parallel Computational Fluid Dynamics, May 13-15, 2003, Moscow, Russia, P.52-54.

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

[26] Kaporin I. Superlinear convergence in minimum residual iterations // Numerical Linear Algebra with Applications, 2005, V.12, P.453-470.

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

[28] 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.

Подписано в печать 26.08.2011 Формат 60x84 1/16

Объем 1,33 п. л. Тираж 100 экз. Заказ № 785

Типография издательства «Фолиум» 127238, Москва, Дмитровское ш., 58

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

Введение.

1 Теория К-числа обусловленности

1.1 Основные свойства функционалов С и К.

1.1.1 Определения и элементарные свойства.

1.1.2 Выпуклость и ортогональные преобразования.

1.1.3 Оценка сходимости метода сопряженных градиентов через К

1.1.4 Различия в свойствах функционалов С и К.

1.1.5 Применение функционалов С и К к несимметризуемым матрицам

1.2 Свойства С и К как функций от спектра собственных значений

1.2.1 Выпуклость и экстремальные значения.

1.2.2 Верхняя оценка К через С.

1.2.3 Связь между значениями К от разных аргументов.

1.3 Оценки собственных значений М через К (М).

1.3.1 Оценки для наименьших собственных значений.

1.3.2 Оценки для наибольших собственных значений.

1.4 Выводы к главе

2 Алгоритмы метода сопряженных градиентов

2.1 Определение и свойства пространств Крылова

2.2 Оценка сходимости через чебышевские многочлены

2.3 Проение методапряженных градиентов (мг.).

2.3.1 Определение итерационных параметров мг.

2.3.2 Ретные формулы мг.

2.3.3 Стандартный алгоритм метода.

2.4 Сходимь итераций мг.

2.4.1 Оценка убывания нормы невязки мг. через С.

2.4.2 Оценка ча итераций мг. через С

2.4.3 Оценка ча итераций мг. через К.

2.5 Блочный метод сопряженных градиентов.

2.6 Проение предобовленного мг.

2.7 Критерии качества предобусловливания

2.7.1 Обобщенная симметричная верхняя релаксация.

2.7.2 Блочная симметричная верхняя релаксация.

2.8 Параллельная реализация крыловских итераций.

2.8.1 Ускорение и эффективность вычислений.

2.8.2 Масштабируемость алгоритмов.

2.8.3 Параллельные архитектуры и модели программирования

2.8.4 Параллельная реализация базовых процедур.

2.9 Выводы к главе 2.

3 Приближенное Т1ти-разложение

3.1 Структура предобусловливадия и оценивание его качества'.

3.2 Предобовливание мг. при помощи приближенного ит11-разложения

3.3 Спектральная устойчивость приближенного ит1Т-разложения

3.4 Алгоритмы приближенного ит11-разложения.

3.5 Оценка заполнения треугольных множителей через параметр отсечения

3.5.1 Структурная устойчивость приближенного 11ти-разложенпя

3.5.2 Оценка числа ненулевых элементов II.

3.5.3 Неулучшаемость оценки заполнения множителя II.

3.6 Верхние оценки обусловленности через параметр отсечения.

3.6.1 Простая оценка нижней границы спектра

3.6.2 Верхняя оценка общих затрат на итерации.

3.7 Базовый алгоритм приближенного ити-разложения

3.8 Безотказные алгоритмы приближенного Х1ти-разложения.

3.8.1 Приближенное ит11-разложение с пробными сдвигами.

3.8.2 Приближенное ити-разложение с пошаговой коррекцией.

3.8.3 Приближенное ит11-разложение с отложенным отсечением

3.8.4 Стабилизированный метод с отложенным отсечением.

3.9 Общие затраты на итерации для методов 2-го порядка.

3.10 Оценка К-числа обусловленности для методов 2-го порядка.

3.11 Численный пример.

3.12 Выводы к главе

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

4.1 Общая структура предобусловливания.

4.1.1 Предварительное масштабирование.

4.1.2 Предварительное переупорядочение.

4.1.3 Предобусловливание при помощи н.о.т.-факторизации.

4.2 Явное предобовливание в алгоритме мг.

4.3 Выбор значений ненулевых элементов матрицы С?.

4.3.1 Доказательство К-оптимальности предобусловливателя.

4.3.2 Оценки достигаемого уменьшения К-числа обусловленности

4.3.3 Многоуровневая структура н.о.т.-разложения.

4.4 Приближенное блочное н.о.т.-предобусловливание.

4.4.1 Структура предобусловливания.

4.4.2 Рекурсивное описание предобусловливания.'

4.4.3 Пример построения для трех блоков.

4.4.4 Связь аддитивной формы с н.о.т.-разложением.

4.4.5 Приближенное блочное н.о.т.-разложение.

4.5 Анализ блочного метода Якоби для модельных задач.

4.6 Численные результаты для тестовой задачи.

4.6.1 Тестовая матрица.

4.6.2 Методика тестирования.

4.6.3 Исследование зависимости от числа блоков р

4.6.4 Исследование зависимости от порога отсечения.

4.6.5 Исследование зависимости от степени перекрытия.

4.7 Результаты для набора тестовых задач.

4.8 Выводы к главе 4.

5 Полиномиальное предобусловливание

5.1 Анализ полиномиального предобусловливания.

5.1.1 Структура полиномиального предобусловливания.

5.1.2 Построение многочленов Чебышева-наименыних квадратов

5.1.3 Свойства многочленов Чебышева-наименыних квадратов

5.1.4 Вывод оценки для K(Mpd-i(M))

5.1.5 Оценка для многочленов Чебышева-наименыних квадратов

5.2 Оценкаодими мг. через верхнюю границуектра и изолированные наименьшиевенные значения.

5.3 Полиномиальное предобовливание мг.

5.3.1 Предобусловливание, использующее полиномы Чебышева.

5.3.2 Алгоритм полиномиального предобусловливания.

5.3.3 Сочетание полиномиального предобусловливания с предварительным.

5.4 Оценка погрешни мг. для полиномиального предобовливания

5.4.1 Полиномиальное предобусловливание и изолированные наименьшие собственные значения

5.4.2 Сравнение полиномиальных предобусловливаний на модельной задаче.

5.4.3 Оценка евклидовой нормы погрешности через Н-норму невязки

5.5 Численные резуль гаты для тестовых за дач

5.5.1 Набор тестовых матриц.

5.5.2 Методика тестирования.

5.5.3 Результаты расчетов тестовых задач.

5.6 Выводы к главе 5.

6 Предобусловливание малоранговой модификацией

6.1 Структура предобусловливания.

6.2 Анализ K-числа обусловленности.

6.2.1 Доказательство К-оптимальности.

6.2.2 Случай (7 = 1 и выбор матрицы С.

6.3 Оценки спектральной обусловленности.

6.3.1 Оценки границ спектра предобусловленной матрицы.

6.3.2 Оценки спектрального числа обусловленности.

6.4 Связь с многоуровневыми и многосеточными методами

6.5 Анализ предобусловливания для модельных задач.

6.6 Расчеты тестовых задач

6.7 Выводы к главе

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

Двумя существенными шагами в развитии итерационных алгоритмов решения систем линейных уравнений Ах = Ь с симметричной положительно определенной матрицей А стали а) разработка Хестенсом и Штифелем в 1952 году [103] метода сопряженных градиентов, который, в свою очередь, тесно связан с алгоритмом тридиагонализации, опубликованным Ланцошем в 1950 году [136]; б) разработка техники предобусловливания для итерационных методов (эквивалентной предварительному линейному преобразованию матрицы системы, например, А —> М, где М = НА, и преобразование задается симметричной положительно определенной матрицей-предобусловлпвателем Н). сформулированной в связи с ускорением метода сопряженных градиентов в работах Даниэля 1967 года [69], Г. И. Марчука и Ю. А. Кузнецова 1968 года [26], О. Аксельсона [39] и др.

Как указывает, в частности, фу н д а мел I т ал ьн а я обзорная работа [94], теория предобусловленного метода сопряженных градиентов в основном развивалась вокруг оценок отношения границ спектра Атах(М)/Ат„(]1/) (т.наз. спектрального числа обусловленности) предобусловленной матрицы М, что порождало соответствующие критерии качества предобусловливания.

Это означало, что не только в теорию, но и в построение методов решения систем линейных уравнений в той или иной мере вовлекалась проблема, еще более сложная по сравнению с исходной, а именно, обобщенная задача на собственные значения Ау = АН~1у. Причем эту задачу требуется решить не в прямой, а в обратной постановке, выбирая параметры предобусловливающей матрицы Н из условия минимума отношения крайних собственных значений матрицы М — НА.

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

Таким образом, в идеале следовало бы найти некоторый матричный функционал, который

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

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

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

Отметим первые попытки разработки нестандартных подходов к предобусловливанию метода сопряженных градиентов, см., напр., [4, 6], а также более поздний обзор Т.Хукле [105], которые опирались на критерий качества предобусловлпвания, выбранный в виде евклидова расстояния \\1п — М\\р в К"2 от предобусловленной матрицы М — НА до единичной 1п. Однако, такому подходу присущи серьезные ограничения, делающие его непрактичным во многих важных ситуациях, см. напр., обсуждение в работе [122]. Прежде всего, если матрица А симметрична, а для матрицы Н фиксирована структура разреженности и ищутся ее элементы, то в результате такой оптимизации получается несимметричная матрица Н, невырожденность которой, как правило, ничем не гарантируется. Так, в [122] указывалось, что для того, чтобы быть уверенным хотя бы в невырожденности М, необходимо потребовать ||/„ — М||,р < 1, а это является чрезмерно жестким условием, непригодным для практического использования. Соответственно, не существует никакой содержательной оценки числа итераций, которая бы зависела только от величины ||/п — М\\р. Следовательно, минимизация указанного евклидова расстояния (до значения, большего единицы) является не более чем эвристическим приемом, за которым не стоит никаких теоретических гарантий. Заметим, что в работе [122] содержится законченная теория сходимости метода обобщенных минимальных невязок в терминах двух величин: во-первых, расстояния ||/п~ Л/1j/г, и, во-вторых, дополнительного параметра, характеризующего отделение спектра М от нуля. Однако, в контексте разреженных матриц большого размера, не представляется реалистичной постановка задачи о минимизации ||/„ — M\\f при условии, например, ограниченности ||А/-1||.

Таким образом, критерии малости евклидова расстояния ||/п — НА\\р сам по себе плохо приспособлен к случаю оптимизации предобусловливаний симметричных положительно определенных матриц. Поэтому в [132] (вслед за работами автора [11, 12, 110], см. также [17], где исследовалось улучшение спектрального числа обусловленности) был рассмотрен случай факторизованного предобусловливателя Н = GTG, где G - нижняя треугольная матрица. К настоящему времени этот класс FSAI-предобусловливателей (от англ. Factorized Sparse Approximate Inverse) достаточно широко распространен в зарубежных работах; соответствующие методы используется, в основном, в высокопараллельных вычислениях. В отличие от оптимизации К-числа обусловленности, лежащей в основе работ автора, в [132, 16, 66] построение предобусловливания снова использует минимизацию евклидова расстояния, но теперь уже в виде ||/n — GL\\f, где А = LLT представляет собой точное симметрично-треугольное разложение исходной матрицы А (англ. Cholesky factorization). Для точного вычисления такой матрицы G оказывается необходимым знать все диагональные элементы которые, естественно, на практике недоступны (в противном случае, можно было бы решить заданную систему быстро и точно, не прибегая к итерационным методам). Поэтому недостающая информация извлекается из естественного требования равенства единице всех диагональных элементов симметрично предобусловленной матрицы GAG1. Получаемое в итоге предобусловливание оказывается, таким образом, субоптимальным с точки зрения минимума \\Iu — GL\\f, однако, как показано в [16], в точности совпадает с К-оптимальным выбором матрицы G, выносимым на защиту в настоящей диссертации. Более того, в той же работе [16] построен пример параметризованной матрицы 2-го порядка, показывающий, что К-оптимизация может давать сколь угодно лучшие результаты, чем безусловная оптимизация евклидовой нормы || In — GL\\f, (в смысле уменьшения спектрального числа обусловленности). В связи с этим, во вводной части упомянутой работы ее авторы справедливо указывают на "потенциальную опасность построения FSAI-предобусловливателей на основе безусловной минимизации фробениусовой нормы".

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

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

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

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

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

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

Таким образом, в условиях взаимообусловленного роста вычислительных мощностей и усложнения постановок прикладных задач, общая проблема построения "наилучшего" метода решения больших разреженных систем линейных уравнений все еще остается далекой от окончательного решения (даже если ограничиться задачами с симметричной положительно определенной матрицей). Заметим, что с увеличением размеров задач для многих важных классов разреженных матриц (возникающих, например, при дискретизации пространственных объектов) методы "точной" треугольной факторизации [82, 177] неизбежно начинают требовать нелинейно растущего объема памяти [139], что и является причиной постоянного интереса к усовершенствованию итерационных методов.

Кроме того, для итерационного решения сложных нелинейных задач могут представлять интерес методы предобусловливания (т.е., построения приближенных обратных матриц) как таковые. Например, эффективность применения методов нелинейной минимизации типа сопряженных градиентов (см., напр., [163]) к большим разреженным задачам на собственные значения [171, 192, 187, 81] или к нелинейным задачам наименьших квадратов [19, 112] сильно улучшается при использовании подходящих предобусловливателей (хотя сколько-нибудь точного решения какой-либо последовательности с.л.а.у. при этом не происходит).

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

1. Предложен и исследован новый критерий эффективности предобусловливания, сформулированный в терминах минимизации матричного функционала К(у1/) = (;г-1 trace M)nf del М, определенного для любой симметризуемой п х n-матрицы М и пропорционального отношению n-й степени следа матрицы к ее определителю. Этот функционал мы называем далее К-числом обусловленности (термин был введен О.Акссльсоном и использован в его известной монографии [41], где обширно цитированы некоторые результаты настоящей диссертации.)

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

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

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

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

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

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

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

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

2. Теория предобусловливания симметричных положительно определенных матриц посредством приближенных треугольных разложений 2-го порядка, с обоснованием как через К-число обусловленности (п. 3.10), так и в терминах спектрального числа обусловленности (п. 3.8.4).

3. Теория предобусловливания симметричных положительно определенных матриц посредством неполных обратных треугольных разложений (п. 4.3), с использованием блочности (п. 4.4.4), а также внутриблочной аппроксимации, с обоснованием в терминах К-числа обусловленности (п. 4.4.5).

4. Теоретическое объяснение несоответствия точной оптимизации спектрального числа обусловленности и получаемого качества полиномиального предобусловливания, полученное на основе анализа К-числа обусловленности; отыскание параметризации чебышевского многочлена, обеспечивающего близкое к оптимальному качество предобусловливания (глава 5).

5. Теория предобусловливания симметричных положительно определенных матриц посредством К-оптимальной малоранговой модификации, с обоснованием как через К-число обусловленности (п. 6.2.1), так и в терминах спектрального числа обусловленности (п. 6.3.2).

Обоснованность и достоверность результатов. Представленные в диссертации результаты имеют строгое математическое обоснование. С другой стороны, достоверность эффективности построенных предобусловливаний основывается на прямом сравнении результатов расчетов с результатами, полученными при использовании других предобусловливаний: точечного метода Якоби, приближенного блочного метода Якоби. Кроме того, для ряда задач из коллекции университета Флориды, использовавшихся другими* авторами для проверки разработанных ими алгоритмов;, были получены, существенно. лучшие результаты. ;. ';,

Апробация, работы. Результаты работы докладывались и обсуждались на 7 всероссийских и 16 международных конференциях, в частности: международной конференции "РаСТ99" (Санкт-Петербург, . Россия, 1999г; г.), всемирном 16-ом Конгрессе "IMACS 2000". по научным вычислениям, прикладной математике и моделированию' (Лозанна, Швейцария, 2000), Голландско-российском, симпозиуме NWO (МГУ, Москва, июнь, 2000),. Симпозиуме NYVO (Амстердам-Ниймеген, Голландия, ноябрь, 2000); .Втором. Международном научно-практическом^ Семинаре* и Всероссийской молодежной школе "Высокопроизводительные параллельные; вычисления' на кластерных системах" (Нижний Новгород, Россия,. 2002), международной, конференции "Parallel CFD 2003"" (Москва-, 2003), Международной конференция; "VIII Забабахинские научные чтения" (РФЯЦ-ВНИИТФ, Снежинск, Россия, 2005), международной конференции по- вычислительной геометрии, генерации сеток .-и научным вычислениям "NUMGRID2008" (ВЦ РАН, Москва, 2008), международной конференции по-вычислительной-.геометрии, генерации сеток и высокопроизводительным вычислениям "NUMGRID2010" (ВЦ РАН, Москва, 2010), международной; конференции по прикладной, математике и информатике; посвященной 100-летию со дня рождения академика А.А.Дородницына; (ВЦ РАН, Москва, 2010); международной, конференции по матричным методам в математике и приложениях "МММА-2011" (ИВМ РАН; Москва, 2011),. научно-исследовательских семинарах Вычислптельного цетттра РАН и Института вычислительной, математики РАН (Москва, 2009-2011), рабочих семинарах ExxonMobil Upstream Research Co. (Хьюстон, США, 2006-2010); '

Публикации. г

По теме диссертации опубликовано 32 работы,, из них 5 в отечественных рецензированных изданиях,, рекомендованных ВАК, 9 в международных рецензируемых изданиях из рекомендованного ВАК списка "Web of Science: Science Citation Index Expanded" (база, по естественным наукам), 2 в других международных рецензируемых изданиях, 1 в российском рецензируемом издании, 4 в? трудах всероссийских конференций, 7 в труд ах. между народных конференций, 4 публикации1 в других научных изданиях. '

В совместных работах [115, 20, 119; 121, 22, 23, 24] И.Н.Коныиину принадлежит разработками реализация параллельного алгоритма метода и описание численных экспериментов, а автору диссертации - разработка и теоретическое исследование предобусловливания.

Благодарности. Автор выражает благодарность А.Ю. Еремину, привлекшему автора к работе по теме диссертации, Л. Ю. Колотилиной за внимание к работе и ценные обсуждения, профессору О. Аксельсону за интерес к работе и всестороннюю ее поддержку на протяжении многих лет, И. Н. Коньшину за длительное сотрудничество в области численного тестирования и параллельной реализации разработанных автором методов, В. А. Гаранже за полезные обсуждения, замечания и поддержку работы, В. 10. Правильникову за сотрудничество в области практической реализации и внедрения разработок по теме диссертации, академику РАН Ю. Г. Евтушенко и член-корреспонденту РАН Е. Е. Тыртышникову за внимание к работе и ее поддержку. Особенную ценность имели обсуждения материала главы 4 с 0. Ю. Милюковой и Ю. В. Василевским, касающиеся оценок вычислительной эффективности предложенных автором методов.

Структура и объем работы. Диссертационная работа состоит из введения, б глав, заключения и списка литературы. Объем работы 216 страниц. Список литературы включает в себя 192 наименования.

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

6.7 Выводы к главе 6

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

Заключение

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

К(М) = (тГ1 trace М)" / det М, для симметризуемой вещественной матрицы А/. В основном изучался предобусловленный метод сопряженных градиентов (м.с.г.) для решения систем линейных алгебраических уравнений (с.л.а.у.) Ах = Ъ с симметричной положительно определенной матрицей коэффициентов А.

С одной стороны, были получены окончательные результаты, описывающие зависимость убывания нормы погрешности м.с.г. с ростом числа итераций от К-числа обусловленности предобусловленной матрицы М = НА, а именно

Ь - Ахк\\н < (K(tf А)1/* - if'2 ||6 - АхоНя

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

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

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

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

1. Богачев К.Ю., Жабицкий Я.В. Метод Капорина-Конышша параллельной реализации блочных предобусловливателей для несимметричных матриц в задачах фильтрации многокомпонентной смеси в пористой среде // Вестник Моск. ун-та. Матем. Механ. 2010, №1, С.46-52.

2. Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. М., Наука, 1984, 320 с.

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

4. Еремин А.Ю., Капорин И.Е. Реализация явных чебышевскпх методов прп решении задач большой размерности. //Многопроцессорные вычислительные структуры: Междуведомственный тематический научный сборник. Таганрог: ТРТИ, 1985, вып. 7 (XVI), С.43-46.

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

6. Еремин А.Ю., Капорин И.Е. Влияние наибольших собственных значений на численную сходимость метода сопряженных градиентов, //Численные методы и вопросы организации вычислений. XIII, Зап. научн. сем. ПОМИ, Т.248,

7. ПОМИ, СПб., 1998, 5-16. (перевод: A. Yu. Yeremin and I. E. Kaporin. The influence of isolated largest eigenvalues on the numerical convergence of the CG method. Journal of Mathematical Sciences, 2000. V.101. N 4. P.3231-3236)

8. Ильин В.П. Методы и технологии конечных элементов. Новосибирск: Изд. ИВМиМГ СО РАН, 2007, 371 с.

9. Капорин И.Е. Альтернативный подход к оценке числа итераций метода сопряженных градиентов. //Численные методы и программное обеспечение (под ред. Ю. А. Кузнецова) М.: ОВМ АН СССР, 1990, С.55-72.

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

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

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

13. Еремин АЛО., Колотилина Л.Ю., Никишин А.А. Переобуславливание с помощью факторизованных разреженных приближенных обратных III. Итерационное построение переобуславливателей. //Зап. научн. семин. ПОМИ, 1998, Т.248, С.17-48.

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

15. Капорин И.Е. Использование предобусловленных подпространств Крылова в методах типа сопряженных градиентов для решения нелинейной задачи наименьших квадратов //Вестн. МГУ, сер.15. Вычисл. матем. и кибернетика. 1995. №3, С.26-31.

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

17. Капорин И.Е. Использование внутренних итераций метода сопряженных градиентов при решении больших разреженных нелинейных задач оптимизации //Журн. вычисл. матем. и матем. физ. 2003, Т.43. С.802-807.

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

19. Kaporin I.E., Ivonshin 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.

20. Лебедев В.И., Финогенов С.А. О порядке выбора итерационных параметров в' чебышевском циклическом итерационном методе. //Журн. Вычисл. Матем. и Матем. Физ., 1971, T.ll, N 2, С.425-438.

21. Марчук Г.И., Кузнецов Ю.А. Об оптимальных итерационных процессах. //Докл. Акад. Наук СССР, 1968, Т.181, N 6, С.1331-1334.

22. Ортега Дж. Введение в параллельные и векторные методы решения линейных систем. М.: Мир, 1991, 367 с.

23. Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978, 592 с.

24. Тыртышников Е.Е. Краткий курс численного анализа. М., ВИНИТИ, 1994, 220 с.

25. Фаддеев Д.К., Фаддеева В.Н. Вычислительные методы линейной алгебры. СПб., изд-во "Лань", 2002, 736 с.

26. Хейгеман JL, Янг Д. Прикладные итерационные методы. М., Мир, 1986, 448 с.

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

28. Armstrong B.A. Near-minimal matrix profiles and wavefronts for testing nodal re-sequencing algorithms. //Int. J. Numer. Meth. Engng. 1985. V.21. P.1785-1790.

29. Ashby S.F., T.A.Manteuffel T.A., Otto.J.S. A Comparison of Adaptive Chebyshev and Least Squares Polynomial Preconditioning for Hermitian Positive Definite Linear Systems. //SIAM J. Sci. Statist. Comput. 1992. V.13. P. 1-29.

30. Ashby S.F. Minimax Polynomial Preconditioning for Hermitian Linear Systems //SIAM. J. Matrix Anal. Appl. 1991, V.12, N.4, P.766-789.

31. Auchmann B. The Coupling of Discrete Electromagnetism and the Boundary Element Method for the Simulation of Accelerator Magnets //CERN-THESIS-2009-017, 02 December 2004, 114 p.

32. Avron H., Gupta A., Toledo S. New Krylov-subspace solvers for Hermitian positive definite matrices with indefinite preconditioners //Tech. Rep. RC 24698 (W0812-001), IBM T.J.Watson Research Center, Yorktown Heights, NY, December 1, 2008, 13 p.

33. Axelsson 0. A generalized SSOR method. //BIT, 1972. V.13. N 4, P.443-467.

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

35. Axelsson 0. Bounds of eigenvalues of preconditioned matrices. //SIAM J. Matrix. Anal. Appl. 1992. V.13. P.847-862.

36. Axelsson O. Iterative solution methods. Cambridge University Press, Cambridge, 1994.

37. Axelsson 0., Neytcheva M., Polman B. The bordering method as a preconditioning method. //Vestnik Moscow Univ., Ser. 15: Vychisl. Mat. Cybern. 1995, P.3-24.

38. Axelsson 0., Kaporin I. On the sublinear and superlinear rate of convergence of conjugate gradient methods. //Numerical Algorithms, 2000. V.25, P.l-22.

39. Axelsson 0., Kaporin I, Optimizing two-level preconditionings for the conjugate gradient method. //Rep.0116, August 2001, Department of Mathematics, Catholic University of Nijmegen, Nijmegen, The Netherlands, 21 p.

40. Axelsson 0., Kaporin I. Error norm estimation and stopping criteria in preconditioned conjugate gradient iterations. //Numerical Linear Algebra with Applications, 2001. V.8. P.265-286.

41. Axelsson 0., Lindskog G. On the rate of convergence of the preconditioned conjugate gradient method //Numer. Math. 1986. V.48, N 5. P.499-523.

42. Axelsson O., Carey G., and Lindskog G. On a class of preconditioned iterative methods on parallel computers. //Int. J. Numer. Meth. Engrg. 1989. V.27. P.637-654.

43. Barnard S.T., Pothen A., and Simon H. A spectral algorithm for envelope reduction of sparse matrices. //Numerical Linear Algebra with Applications, 1995. V.2. P.317-334.

44. Basermann A., Reichel B., Schelthoff C., Preconditioned CG methods for sparse matrices on massively parallel machines. //Parallel Computing. 1997. V.23, N 3, P. 381-398.

45. Bateman H., Erdelyi A., Higher transcendental functions, Vol.2, McGraw Hill, Inc., N.Y.etc., 1953.

46. Beck E. Uber die differenz zwischen dem arithmetishen und dem geometrishen miltel von zallien, //Math.-Phys.Semesterber, 1969. V.16. P.117-123.

47. Beckermann B., Kuijlaars A.B.J. Superlinear convergence of conjugate gradients. //SIAM J. Numer. Anal. 2001. V.39, N 1, P.300-329.

48. Behie A., Forsyth P.A. Comparison of fast iterative methods for symmetric systems, //Inst. Math. Applies. J. Num. Anal. 1983. V.3. P.41-63.

49. Benzi M., Haws J.C., Tuma M. Preconditioning highly indefinite and nonsymmetric matrices. //SIAM J. Sci. Comput. 2000. V.22, P. 1333-1353.

50. Benzi M., Tuma T. A robust incomplete factorization preconditioner for positive definite matrices. //Numer. Linear Algebra Appl. 2003. V.10. P.385-400.

51. Bollhofer M., Saad Y. Multilevel preconditioners constructed from inverse-based ILUs. //SIAM J. Sci. Comput. 2006. V.27. P.1627-1650.

52. Bridson R., Tang W.-P. A structural diagnosis of some IC orderings. //SIAM J. Sci. Comput. 2000. V.22. P.1527-1532.

53. Bunch J.R. and Parlett B.N., Direct methods for solving symmetric indefinite systems of linear equations. //SIAM J. Numer. Anal. 1971. V.8. P.639-655.

54. Cai X.-C., Dryja M., Sarkis M. Restricted additive Schwarz preconditioners with harmonic overlap for symmetric positive definite linear systems //SIAM J. Numer. Anal. 2003, V.41, N 4, P.1209-1231.

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

56. Canga M.E., Becker E.B. An iterative technique for the finite element analysis"of near-incompressible materials, //Comm. Meths Appl. Mech. Engrg. 1999, V.170, P. 79-101.

57. Chin P., Forsyth P.A. Comparison of GMRES and CGSTAB accelerations for incompressible viscous flow, //J. Comp. Appl. Math. 1993. V.46. P.415-426.

58. Chow E. A priori sparsity patterns for parallel sparse approximate inverse preconditioners. //SIAM J. Sci. Comput. 2000. V.21, N 5, P.1804-1822.

59. Chow E. Parallel implementation and practical use of sparse approximate inverse preconditioners with a priori sparsity patterns. //International Journal of High Performance Computing Applications, 2001, V.15, N 1, P.56-74.

60. Cuthill E., McKee J. Reducing the bandwidth of sparse symmetric matrices. //Proceedings 24th National Conference of the Association for Computing Machinery, Brandon Press, New Jersey, 1969, P. 157-172.

61. Daniel J.W. The conjugate gradient method for linear and nonlinear operator equations. //SIAM J. Numer. Analysis, 1967. V.4. P. 10-26.

62. Davis T.A., Hu Y.F. University of Florida sparse matrix collection. To appear in: ACM Trans, on Math. Software, 2011. V. 38 http://www.cise.ufl.edu/research/sparse/matrices

63. Dennis J.E., Wolkowicz H. Sizing and least change secant methods //SIAM J. Nu-mer. Anal. 1993. V.10. P.1291-1314.

64. Dewilde P., Deprettere E.F. Approximate inversion of positive matrices with applications to modelling //Modelling, Robustness, and Sensitivity Reduction in Control Systems (R.F. Curtain, ed.) NATO ASI Series, 1984, V.F 34. P.211-238.

65. Dickinson J.K., Forsyth P.A., Preconditioned conjugate gradient methods for three-dimensional elasticity //Int. J. Numer. Methods Eng. 1994. V.37. P.2211-2234.

66. Diyankov O.V., Koshelev S.V., Kotegov S.S., Krasnogorov I.V., Kuznetsova N.N., Pravilnikov V.Y., Beckner B.L., Usadi A.K. SparSol: Sparse linear systems solver // Russian J. Numer. Analysis and Math. Modelling, 2007. V.22, N 4, P.325-339.

67. Dubrulle A.A. Retooling the method of block conjugate gradients. //Electronic Transactions on Numerical Analysis. 2001. V.12. P.216-233.

68. Duff I.S. User's guide for the Harwell-Boeing sparse matrix collection (Release I). //Technical Report TR/PA/92/86, CERFACS, October 1992.

69. Duff I.S., Scott J.A. Stabilized block diagonal forms for parallel sparse solvers. //Rutherford Appleton Laboratory Technical Report RAL-TR-2004-006 CCLRC, February 2004 (revised November 2004).

70. Duff I.S., Koster J. The design and use of algorithms for permuting large entries to the diagonal of sparse matrices. //Rutherford Appleton Laboratory Technical Report RAL-TR-97-059 October 24, 1997

71. Duff I.S., Koster J. On algorithms for permuting large entries to the diagonal of a sparse matrix. //Rutherford Appleton Laboratory Technical Report RAL-TR-1999-030 April 19, 1999

72. D'yakonov E.G. Optimization in solving elliptic problems. CRC Press, Boca Raton, FL, 1996. Translated from the 1989 Russian original, Translation edited and with a preface by Steve McCormick.

73. 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. V.18. P. 11451151.

74. Everstine G.C. A comparison of three resequencing algorithms for the reduction of matrix profile and wavefront. //Int. J. Numer. Meth. Engng. 1979. V.14. P.837-853.

75. Feng Y.Y., Owen D.R.J., Peric D. A block conjugate gradient method applied to linear systems with multiple right-hand sides //Computer Meth. in Applied Mech. Engrg. 1995. V.127. P.203-215.

76. Fiedler M. A property of eigenvectors of nonnegative symmetric matrices and its application to graph theory. //Czecheslovak Mathematical Journal. 1975. V.25. P.619-633.

77. Field M.R. Adaptive polynomial preconditioning for the conjugate gradient algorithm. //Applied Parallel Computing Computations in Physics, Chemistry and Engineering Science. Lecture Notes in Computer Science, 1996, V.1041, P.189-198.

78. Frank J., Vuik C. On the construction of deflation-based preconditioners. //SLAM J. Sei. Comput. 2001, V.23. N 2. P.442-462.

79. Fritzsche D., Frommer A., Szyld D.B. Extensions of certain graph-based algorithms for preconditioning. //SIAM Journal on Scientific Computing. 2007. V.29, N 5, P.2144-2161.

80. George A., Liu J. W. II. An implementation of a pseudoperipheral node finder. //ACM Trans. Math. Softw. 1979. V.5. P.284-295.

81. George T., Gupta A., Sarin V. An experimental evaluation of iterative solvers for large SPD systems of linear equations //IBM Research Report RC 24479, Jan.25, 2008.

82. Gibbs, N. E. A hybrid profile reduction algorithm. //ACM Trans. Math. Softw. 1976. V.2. P.378-3S7.

83. Gibbs N.E., Poole W.G.,Jr., Stockmeyer P.K. An algorithm for reducing the bandwidth and profile of a sparse matrix //SIAM J. Numer. Anal. 1976. V.13. P.236-250.

84. Giraud L., Tuminaro R.S. Algebraic Domain Decomposition Preconditioners //Technical Report ENSEEIHT-IRIT RT/APO/06/07, Toulouse, France, October 2006, 24 p.

85. Golub G.H., O'Leary D.P. Some history of the conjugate gradient and Lanczos algorithms: 1948-1976. //SIAM Review, 1989. V.31, P.50-102.

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

87. Greenbaum A., Strakos Z. Predicting the behavior of finite precision Lanczos and conjugate gradient computations //SIAM J. Matrix Anal. Appl. 1992. V.13. P.121-137.

88. Greenbaum A., Li C., Chao H.Z. Parallelizing preconditioned conjugate gradient algorithms. //Comput. Phys. Comm., 1989. V.53, N 1-3. P.295-309.

89. Gonzalez L. Orthogonal projections of the identity: spectral analysis and applications to approximate inverse preconditioning //SIAM Review, 2006. V.48, N 1, P. 66-75.

90. Greenbaum A. Comparison of splittings used with the conjugate giadient algorithm //Numer Math. 1979. V.33. P.181-194.

91. Gupta A., Ying L. A fast maximum-weight-bipartite-matching algorithm for reducing pivoting in sparse Gaussian elimination. //IBM T. J. Watson Research Center., Yorktown Heights, NY, Tech. Rep. RC-21576 (97320), October 1999.

92. Gustafson K.E. Operator trigonometry of iterative methods, //Numer. Linear Algebra Appls. 1997. V.4, P.333-347.

93. Gustafson K.E. A computational trigonometry, and related contributions by russians Kantorovicli, Krein, Kaporin //Вычислительные технологии. 1999. Т.4. №3. С.73-83.

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

95. Hladik I., Reed M.B., Swoboda G. Robust preconditioners for linear elasticity FEM analyses, //Int. J. Numer. Methods Eng., 1997. V.40. P.2109-2127.

96. Huckle T. Factorized Sparse Approximate Inverses for Preconditioning //The Journal of Supercomputing. 2003. V.25. N 2. P.109-117.

97. Jennings A. Influence of the eigenvalue spectrum on the convergence rate of the conjugate gradient method, //Journal of the Institute of Mathematics and Its Applications. 1977. V.20. P.61-72.

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

99. Johnson O.G., Micclielli C.A., and Paul G. Polynomial Preconditioning for Conjugate Gradient Calculations, //SIAM J. Numer. Anal. 1983. V.20. P.362-376.

100. Jones M.T., Plassman P.E. An improved incomplete Cholesky factorization. //АСМ Trans. Math. Software, 1995. V.21. P.5-17.

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

102. Kaporin I.E. New convergence results and pieconditioning strategies for the conjugate gradient method. //Numer. Linear Algebra Appls. 1994. V.l. N 2. P.179-210.

103. Kaporin I.E., Axelsson 0. On a class of nonlinear equation solvers based on the residual norm reduction over a sequence of affine subspaces, //SIAM J. Sci. Com-put.1994. V.16. N 1. P.228-249.

104. Kaporin I.E. High quality preconditionings of a general symmetric positive definite matrix based on its UTU + UTR-\- Rri7-decomposition. //Numer. Lin. Alg. Appl. 1998. V.5. P.4S3-509.

105. Kaporin I.E. Second order incomplete Cholesky preconditionings for the CG method, //Proc. of IMMB'98, Iterative solution methods for the elasticity equations as arising in Mechanics and Biomechanics, Nijmegen, Netherlands, Sept. 28-30, 1998, P.47-49.

106. Kaporin I.E., Konshin I.N. Benchmark problems in linear elasticity: parallel solution // Rep.0028, December 2000, Department of Mathematics, Catholic University of Nijmegen, Nijmegen, The Netherlands, 23 p.

107. Kaporin I.E., Konshin I.N. A parallel block overlap preconditioning with inexact submatrix inversion for linear elasticity problems. //Numerical Linear Algebra with Applications, 2002. V.9. P.141-162.

108. Kaporin I.E. Using the Modified 2nd Order Incomplete Cholesky Decomposition as the Conjugate Gradient Preconditioning. //Numerical Linear Algebra with Applications, 2002, V.9. P.401-408.

109. Kaporin,I.E., Konshin,I.N. Parallel conjugate gradient preconditioning via incomplete Cholesky of overlapping submatrices //Book of Abstracts, Parallel Computational Fluid Dynamics, May 13-15, 2003, Moscow, Russia, P.52-54.

110. Kaporin I.E. Second order ILU preconditioning for nonsymmetric matrices. //International Conference on Matrix Methods and Operator Equations (June 20 25, 2005, Moscow) - Abstracts, P.14-15.

111. Kaporin I. Scaling, Reordering, and Diagonal Pivoting in ILU Preconditionings. //Russian Journal of Numerical Analysis and Mathematical Modelling. 2007. V.22. N 4, P.341-375.

112. Kaporin I.E. Incomplete ILU factorization of general sparse matrices. //2nd International Conference on Matrix Methods and Operator Equations (July 23 27, 2007, Moscow) - Abstracts, P.31-33.

113. I. E. Kaporin, Scaling, Preconditioning, and Superlinear Convergence in GMRES-type iterations. //Matrix Methods: Theory, Algorithms, Applications (V.Olshevsky, E.Tyrtyslmikov, eds.), World Scientific Publ., 2010, P.273-295.

114. Karypis G., Kumar V. Multilevel k-way hypergraph partitioning // Techn. Rept. 98-036. Minnesota: Dept. Comput. Sci. Engrgs. Army HPC Res. Center, Univ. Minnesota, 1998.

115. Kershaw D.S. The Incomplete Cholesky-Conjugate Gradient Method for the Iterative Solution of Systems of Linear Equations. //J. Comput. Physics. 1978. V.26, P. 43-65.

116. Kershaw D.S. On the problem of unstable pivots in the incomplete LU-conjugate gradient method. //J.Comput. Physics. 1980. V.38. P.114-123.

117. Kolotilina L.Yu., Yeremin A.Yu. Factorized sparse approximate inverse preconditionings I. Theory //SIAM J. Matrix Anal. Appl. 1993. V.14. P.45-58.

118. Kolotilina L.Yu., Nikishin A.A., Yeremin A.Yu. Factorized sparse approximate inverse preconditionings IV. Simple approaches to rising effciency //Numer. Linear Alg. Appl. 1999. V.6. P.515-531.

119. Kuijlaars A.B.J. Convergence analysis of Krylov subspace iterations with methods from potential theory. //SIAM Review, 2006. V.48, N 1. P.3-40.

120. Kumfert G., Pothen A. Two improved algorithms for envelope and wavefront reduction // BIT 1977. V.18. P.559-590.

121. Lanczos C. An iteration method for the solution of the eigenvalue problem of linear differential integral operators //J. Research Nat. Bur. Standards, 1950. V.45. P.255-282.

122. Lewis J. G. Implementation of the Gibbs-Poole-Stocknieyer and Gibbs-King algorithms // ACM Trans. Math. Softw. 1982, V.8, P.180-194.

123. Lindskog G., Gustafsson I. Pieconditioning by incomplete factorization: an overview and an application to a linear elasticity problem //Report 2001:78, Department of Mathematics, Chalmers University of Technology, Goteborg, Sweden, 2001

124. Lipton R.J., Rose D.J., Tarjan R.E. Generalized nested dissection // SIAM J. Nu-mer. Anal. 1979. V.16. P.346-358.

125. Li N., Saad Y., Chow E. Crout versions of ILU for general sparse matrices. // SIAM J. Sci. Comput. 2003. V.25. P.716-728.

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

127. Mandel J. Balancing domain decomposition // Communications in Applied and Numerical Methods, 1993. V.9. P.233-241.

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

129. Marcus M., Mine H. A Survey of Matrix Theory and Matrix Inequalities. Allyn and Bacon, Inc., Boston. 1964.

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

131. Meinardus G. Uber Eine Verallgemeinerung Einer Ungleichung von L.V.Kantorowitsch. // Numer. Math. 1963. V.5. P.14-23.

132. Migallon V., Penades J., Szyld D.B. Nonstationary multisplitting with general weighting matrices, // SIAM J. Matrix Anal. Appl. 2001. V.22. N 4. P.1089-1094.

133. Milyukova O.Yu. Parallel approximate factorization method for solving discrete elliptic equations // Parallel Comput. 2001. V.27. P.1365-1379.

134. Milyukova O.Yu. Parallel version of approximate factorization method for solving 2D and 3D elliptic equations //J. of Comput. Methods in Science and Engineering. 2002. V.2. N 1-2, P. 195-200.

135. Mroueh H., Shahrour I. Use of sparse iterative methods for the resolution of three-dimensional soil/structure interaction problems, //Int. J. Numer. Anal. Meths. En-grg. 1999. V.23, P.1961-1975.

136. Munksgaard N. Solving sparse symmetric sets of linear equations by preconditioned conjugate gradients. //ACM Trans. Math. Software. 1980. V.6. P.206-219.

137. Nataf F., Xiang H., Dolean V. A two level domain decomposition preconditioner based on local Dirichlet to Neumann maps // C.R.Mathematique Acad.Sci.Paris, 2010. V.348, N 21-22. P.1163-1167.

138. Nicolaides R.A. Deflation of conjugate gradients with applications to boundary value problems. // SIAM J. Numer. Anal. 1987. V.10. N 2. P.355-365.

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

140. Notay Y. An efficient parallel discrete PDE solver // Parallel Comput. 1995. V.21. P. 1725-1748.

141. O'Leary D.P. The block conjugate gradient algorithm and related methods. // Linear Algebra Appl. 1980. V.29. P.293-322.

142. Olschowska M., Neumaier A. A new pivoting strategy for Gaussian elimination. // Linear Algebra Appl. 1996. V.220. P.131-151.

143. O'Neil J., Szyld D.B. A block oidering method for sparse matrices // SIAM J. Sci. Statist. Comput. 1990. V.ll. P.811-823.

144. Padiy A., Axelsson 0., Polman B. Generalized augmented matrix preconditioning approach and its application to iterative solution of ill-conditioned algebraic systems. //SIAM J. Matrix Anal. Appl. 2000. V.22. N 3. P.793-818.

145. Paulino G.H., Menezes I.V.M., Gattass,M., Mukherjee,S. A new algorithm for finding a pseudoperipheral vertex or the endpoints of a pseudodiameter in a graph // Communications in Numer. Meth. Engng. 1994. V.10. P.913-926.

146. Polya G., Szego G. Problems and Theorems in Analysis, Vol.1, Springer-Verlag, Berlin, New York, 1972.

147. Pothen A., Simon H., Liou K.-P. Partitioning sparse matrices with eigenvectors of graphs. // SI AM J. Matrix Anal. Appl. 1990. V.ll. P.430-452.

148. Qian II., Sapatnekar S. Stochastic preconditioning for diagonally dominant matrices // SIAM J. Sci. Stat. Comput. 1988. V.9. N 3. P.152-163.

149. Radicati di Brozolo G., Robert Y. Parallel conjugate gradient-like algorithms for solving sparse nonsymmetric linear systems on a vector multiprocessor, // Parallel Comput. 1989. V.ll. P.223-239.

150. Rao V.N., Kumar V. Parallel depth-first search, part II: Analysis // International Journal of Parallel Programming, 1987. V.16. N 6. P.501-519.

151. Reid J.K., Scott J.A. Ordering symmetric sparse matrices for small profile and wavefront // Rep. no. RAL-TR-98-016. Computing and Information Systems Department, Rutherford Appleton Laboratory, Chilton, Feb. 1998.

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

153. Saad Y. Practical use of polynomial preconditionings for the conjugate gradient method // SIAM Journal of Scientific and Statistical Computing. 1985. V.6. N 4. P.865-881.

154. Saad Y. Numerical Methods for Large Eigenvalue Problems, Halsted Press, New York, 1992.

155. Saad Y. ILUT: a Dual Threshold Incomplete LU Factorization. // Numer. Linear . Algebra Appls. 1994. V.l. P.387-402.

156. Saad Y. Iterative Methods for Sparse Linear Systems. PWS Publishing Co., Boston-New York, 1996.

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

158. Saad Y. Multilevel ILU with reorderings for diagonal dominance. // SIAM J. Sci. Comput., 2005. V.27, P.1032-1057.

159. Saint-Georges P., Warzee G., Notay Y., Beauwens R. Problem-dependent preconditioned for iterative solvers in FE elastostatics // Comp. Struct. 1999. V.73. P.33-43.

160. Schenk 0., Garther K. Solving unsymmetric sparse systems of linear equations with PARDISO // Future generation computer systems, 2003.

161. Sloan S.W. An algorithm for profile and wavefront reduction of sparse matrices// Int. J. Numer. Meth. Engrg. 1986. V.23. P.239-251.

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

163. Tang W.-P. Generalized Schwarz splittings. //SIAM J. Sci. Statist. Comput. 1992. V.13. P.573-595.

164. Vassilevski Yu. A hybrid domain decomposition method based on aggregation. //Numer. Linear Algebra Appls. 2004. V.ll, P.327-341.

165. Vollan A., Kaporin I., Babikov P. Practical Experience with Different Iterative Solvers for Linear Static and Modal Analysis of Large Finite Element Models. //Proc. of the 21st MSC European Users' Conf., Italian Session, September, 1994.

166. Winther R. Some superlinear convergence results for the conjugate gradient method. //SIAM Journal on Numerical Analysis, 1980. V.17. P.14-17.

167. Wolkowicz H., Styan G.P.H. Bounds for eigenvalues using traces, //Lin. Alg. Appls. 1980. V.29. P.471-506.

168. Yamazaki I., Bai Z., Chen W., Scalettar R. A high-quality preconditioning technique for multi-length-scale symmetric positive definite linear systems. //Numer. Math. Theor. Meth. Appl. 2009. V.2, N 4, P.469-484.

169. Yang H. Conjugate gradient methods for the Rayleigh quotient minimization of generalized eigenvalue problems //Computing, 1993. .V.51. N 1. P.79-94.