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

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

• 3 од

■ > . На правах рукописи

ЛАРИН Максим Рудольфович

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

01.01.07- вычислительная математика

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

Новосибирск 1995

Работа выполнена в Вычислительном Центре СО РАН.

Научные руководители:

доктор физико-математических наук профессор Ильин В.Г1. кандидат физпко-математичеашх наук с.н.с.

Карначук В.И.

Официальные оппоненты:

доктор физико-математических наук профессор Воеводин А.Ф. доктор физико-математических наук профессор Мацокнн А.М.

Ведущая организация: Институт математики СО РАН.

Защита состоится _ 1996 г.

в /5№ часов на заседании Специализированного совета К 002.10.01 в Вычислительном Центре СО РАН по адресу: Новосибирск, 030090, Проспект Лаврентьева, 0, Вычислительный Центр.

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

Автореферат разослан " '12." _ 1996 г.

Ученый секретарь Специализированного совета доктор физико-математических наук

I . Ю.И. Кузнецов

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

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

Ах = Ь, (1)

где А - разреженная симметричная положительно определенная матрица, прямыми и итерационными методами.

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

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

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

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

Практическая ценность. Результаты диссертации находят приме-

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

Апробация работы. Основные результаты диссертации докладывались на Всероссийской конференции по комплексам программ математической физики (Новосибирск, 1992), 8-ой Всесибир-ской школе по методам вычислительной математики (Шушенское, 1993), международном семинаре "Numerical Linear Algebra with Applications" ("Численные методы линейной алгебры и приложения") (Неймеген, Нидерланды, 1993), конференции молодых ученых ВЦ СО РАН (Новосибирск, Россия, 1995), международной конференции "Advanced Mathematics, Computations and Applications" ("Современные проблемы прикладной и вычислительной математики") (Новосибирск, Россия, 1995), на научных семинарах ВЦ СО РАН.

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

Структура и объем работы. Диссертационная работа состоит из введения, трех глав и списка литературы. Объем работы - 168 машинописных страниц. Список литературы содержит 90 наименований.

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

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

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

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

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

Ах=Ь, (1.1.1)

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

Теорема 1.1. Если А - N х ./V, симметричная положительно определенная матрица, то существует и единственно треугольное разложение ЬЬТ, где Ь - нижняя треугольная матрица с положительными диагональными элементами.

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

Теорема 1.2. Число операций, необходимых для вычисления треугольного множителя Ь матрицы А, равно

1 ¿=1

где г}(Ь„{) - число ненулевых элементов г-го столбца множителя Холесского.

В §1.2 формулируется основная проблема применения метода Холесского - заполнение, т.е. появление новых ненулевых элементов,

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

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

Теорема 1.3. Пусть А - несингулярная матрица, имеющая треугольное разложение А = LU, тогда всякий ненулевой элемент матрицы (L + U)ij принадлежит оболочке Env(A).

Для любой матрицы перенумерации Р матрица РАРТ также является положительно определенной и метод Холесского по-прежнему применим. Таким образом, вместо системы (1) можно решать эквивалентную систему

(.РАРт){Рх) = РЬ, (1.2.1)

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

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

Одним из лучших методов построения матрицы Р является метод вложенных сечений для матрицы, ассоциированной с конечно-элементной (п х п)-сеткой.

Теорема 1.10. Число ненулевых элементов треугольного множителя L матрицы А, ассоциированной с регулярной п х п-сеткой,

которая упорядочена посредством метода вложенных сечений, выражается формулой

^n2log2n + 0(п2).

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

829 з 2,

Эти оценки показывают, что как с точки зрения динамического заполнения, так и количества арифметических операций, необходимых для вычисления множителя Холесского, метод вложенных сечений превосходит другие известные способы, включая и те из них, которые основаны на ленточных матрицах, поскольку заполнение уменьшается с обычных 0(п3) для ленточных матриц до 0(п2 log2 п) и, соответственно, количество арифметических операций с 0(п4) до 0(п3).

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

В §1.5 рассматриваются основные особенности блочного метода Холесского в отличие от оригинального метода, его различные алгоритмы и соответствующие способы решения блочных треугольных систем.

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

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

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

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

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

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

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

уровень = 1,.теквер = корень дерева, предвер = корень дерева.

2. Опускаемся из теквер в его первого левого потомка, ранее не рассматривавшегося.

предвер(теквер) = теквер, теквер = выбраной вершине,

уровень = уровень + 1.

Если уровень < //, то шаг 2, иначе шаг 3.

3. Обозначим количество узлов сетки из теквер через к, а их номера - через р1,...,рь- Разлагаем столбцы с номерами р1, ...,рк и вносим соответствующие модификации.

4. Выбираем новую текущую вершину:

Если рассмотрены не все потомки предеср(текаер), то шаг 2, иначе теквер = предвер(теквер), предвер — предвер(предвер), уровень = уровень — 1.

5. Для всех узлов из теквер выполняется разложение соответствующих столбцов матрицы

6. Если уровень > 1, то шаг 2, иначе конец.

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

В §2.3 предлагается квадросхема хранения разреженного множителя Холесского, даются оценки памяти и числа арифметических операций, проводится анализ численных экспериментов.

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

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

к а а и)

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

к - параметр, определяющий расстояние от данной вершины до корня дерева (для корневой вершины к = 1);

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

сг - указатель на численные значения блока 5,, привязанного к данной вершине;

и> - указатель на численные значения соответствующих блоков

Сами служебные коды, а также блоки 5 и V, хранятся в отдельных массивах АН, 5Л и У Л таким образом, что данные для потомков одной и той же вершины расположены последовательно (рис. 2). В предлагаемой схеме хранения квадратный блок 5 и прямоугольные блоки VI,...,VI представлены в виде регулярных (гп х т) и (п х га) массивов соответственно (при этом учтена симметрия блока 5). Для блоков ...,А\ используется аналогичный (структурный) способ хранения, как и для исходной матрицы. Таким образом, получается иерархическая схема хранения, совпадающая с построенным выше квадродере-

вом. Это дает значительное преимущество при реализации древовидного алгоритма Холесского.

АЯ

00

10

20

30

40

11

12

13

44

500 ^20 ¿30 5-40 5и !?4 4

^00 Ко у20 Узо у40

Рис.2. Древовидная схема хранения

Практические результаты показывают, что предложенная схема хранения разреженной матрицы оказывается более эффективной по сравнению с другими схемами, применяемыми на практике. Однако, предлагаемая схема хранения требует 0(п3) памяти, при реальном заполнении 0(п21о§2 п).

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

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

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

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

Проведенные теоретические и экспериментальные исследования показали, что наилучшими характеристиками обладает асимметричный вариант кросс-метода, который позволяет сократить как объем основной памяти до 0(n2log2n), так и количество вычислений в 3.7 раза.

Теорема 2.3. Объем памяти, требуемый для хранения ненулевых элементов треугольного множителя L, полученного с помощью аси-метричного варианта кросс-метода, матрицы А, ассоциированной с регулярной (п х п)-сеткой, выражается формулой

S(n) = 2п2 log2 п + 0(п2).

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

-n3 + 0(n2 log,n).

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

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

Предобуславливающая матрица К метода неполной факторизации для Стильтьесовых матриц вида

А = D — L - U, (3.1.6)

где D — {£>,} - блочно-диагональная матрица, L — {¿¿} - блочная нижняя треугольная матрица и U — {[/,} - блочная верхняя треугольная матрица, может быть записана в следующей форме:

К ={G-L)G~\G-U) (3.1.8)

где G = {С?,-} - блочно-диагональная матрица той же структуры, что и матрица D. Матрицу G можно искать в форме

G\ = —D\ — tSi, S\e = — ——-.Diei,

и UJ

Gi = -Д - rSi - ХДГ1 £/,._,,

ш (3.1.15)

S^ = {LiGj}xUi-\ - LiG^U^ -

u)

i = 2,..., I,

где 5 - диагональная матрица; С - некоторая "аппроксимация" матрицы С, т.е. матрица, имеющая ненулевыми элементами (совпадающими с соответствующими элементами б-1) только те, которые расположены на тех же позициях, что и ненулевые элементы в С; е - вектор-столбец с единичными элементами, а ы и т - числовые итерационные параметры, удовлетворяющие условиям 0 < ш < 2, О < т < 1.

Теорема 3.3. Для метода сопряженных градиентов с предобусла-вливающей матрицей К, определяемой по формулам (3.1.8) «(3.1.15) при и> — т = 1, необходимое число итераций п(е) оценивается величиной

п{е)< 2^1^+1. (3.1.28)

Анализ полученных теоретических и экспериментальных результатов позволяет сделать следующий главный вывод:

Из теоретических оценок (см.формулу (3.1.28)) следует, что процесс сходится со скоростью 0{Д0-1). Однако экспериментальные данные показывают более высокую скорость сходимости за одну-две итерации.

В §3.2 предлагается многоуровневый итерационный метод построения" предобуславливающих матриц к симметричным положительно определенным матрицам, возникающим при конечно-элементной аппроксимации двумерных задач математической физики. Суть предложенного метода заключается в специальном способе аппроксимации последовательности матриц {Л'1'}, к — 0,1,...,£, основанном на сеточной природе матриц, с тем чтобы получающиеся дополнения Шура имели заранее спрогнозированное заполнение, структура которого соответствует некоторой новой треугольной сетке. В результате анализа качества полученной предобуславливающей матрицы были

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

Для создания предобуславливающей матрицы М построим последовательность матриц Л'*), к = 0,1, ...,£— возрастающего порядка Пк, причем А™ = А.

Представим матрицу к > 0 в следующем блочном виде:

А(Ш) _

¿п+1) Л^] }ХШ\Х»

(3.2.5)

}Хк

41 л22

Аппроксимируем симметричную положительно определенную матрицу А'*+1) посредством положительной диагональной матрицы такой, что

4*+1,е = А^+1)е, (3.2.6)

где е - вектор-столбец порядка пи+1 — тц, состоящий из единиц. Определим матрицу равенством

АМ= А£+1) - АГ1,Л(п+1)"^(12+1>. (3-2.8)

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

Предобуславливающая матрица М определяется рекурсивно как последовательность предобуславливающих матриц {М^} к последовательности матриц {АМ} следующим образом:

М(°> = А«».

Для к = 1,2,...,Ь-1,Ь

м(*+1> =

ИГ" О'

А&+1) I

Т Т<*:+1)-1 ЛА+1)-, 1 Л11 Л12

(3.2.9)

где - определенная выше матрица и

(3.2.10)

где Рщ{х) - полином степени щ, Р„к(0) = 1, являющийся минимальным на промежутке = [а/., 6*], содержащем все собственные числа матрицы и определяемый как

V Ьк-ак ) + 1

71 (Ь±2к.\ 4. 1 ' \Ьк-ак) + 1

(3.2.24)

где - чебышевские полиномы 1-го рода степени и:

То — 1, Т1 = I, Т„+1 = 2ЬТ„ — Т„_1.

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

Теорема 3.6. Двухуровневый метод разложения для конечно-элементных матриц, базирующийся на последовательности матриц определяемых по формулам (3.2.5), (3.2.6) и (3.2.8), а также последовательности предобуславливаюгцих матриц {Мрекурсивно определяемых по (3.2.9), (3.2.10) и (3.2.24), имеет, оптимальную скорость сходимости и оптимальный порядок вычислительных затрат, если

О

где Оз и - определенные выше константы, р - нижняя граница отношения числа узлов на к + 1 -ом уровне к числу узлов на к-ом уровне, и - степень матричного полинома, используемого в (3.2.10), и и > 1 для каждого (ц + 1 )-го шага.

Таким образом, выбирая соответствующие степени матричных полиномов 1>ь, получаем оптимальную скорость сходимости, т.е. число обусловленности матрицы М~1А есть величина порядка 0(1), и оптимальный порядок общих вычислительных затрат, т.е. число операций пропорционально числу узлов исходной сеткн.

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

Пусть = А - симметричная положительно определенная М-матрица. Для определения последовательности матриц {А^} представим матрицу > 0, в следующем блочном виде:

Лп Л у2

Л") Л)

1Л21 Л 22 J

(3.3.5)

Введем в рассмотрение дополнение Шура матрицы А^

5(*> = 4« - 4М,г142) = - аЯ^Г1^ -вм = лВД^-'сВД^ - СК,Г14*2).

(3.3.8)

где С]]' - диагональная часть матрицы А[\\ а Си' - ее внедиагональ-ная часть. Определим матрицу Л(А'+1> как аппроксимацию матрицы

А{м) = 42 - Г14? - (3-3.9)

где в - числовой итерационный параметр, 0 < в < 1, а -

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

А(*+1)е = 5Ме прп #= 1.

(3.3.10)

Далее, применяем аналогичную процедуру к матрице А^к\ и повторяем этот процесс до тех пор, пока матрица Л'^' не станет трех-диагональной. Отметим, что полученная последовательность матриц {Л^} есть последовательность матриц Стильтьесовского типа.

Предобуславливаюхцая матрица М определяется по рекурсии сверху-вниз следующим образом:

л (к) Л11

М<*> =

м Л21

Ап Л12

0 М(*+1>

(3.3.12)

При исследовании скорости сходимости метода получено следующее равенство для отношения Рэлея

хтАх

__ 1

хтМх 1 -г'

(3.3.22)

где

I

т,- = х-('+1)Г(5,') - Л(,+1))х(,+1'. (3.3.23)

х Л.Г ,_о

В силу свойств матриц С^41' — £(*'+') и (¿М, получаем следующую оценку для числа обусловленности матрицы Л/"1 Л:

«<-1+1(1_"б)а, «>0, Р > 0, (3.3.27)

где

^ Л .Г

(3.3.28)

К сожалению, т полученной оценки не видно, что скорость сходимости многоуровневого метода неполной факторизации для двумерных сеточных задач есть величина 0(/г-1), где И - характерный размер сетки. Наоборот, наличие в знаменателе множителя хтАх при независимости от Л норм матриц и С)^ — В'1' указывает на то, что а и /3 есть величины 0(Л-2), и следовательно, скорость сходимости многоуровневого метода неполной факторизации также оценивается величиной 0(к~2). Последнее противоречит полученным экспериментальным результатам, которые приводятся в заключительном разделе параграфа. Таким образом, для получения соответствия между результатами теории и практики необходимо уметь более точно оценивать значения т, или, что, то же самое, нормы матриц и

<2« - В«К

Тем не менее, данная оценка для скорости сходимости методов неполной факторизации имеет несомненное преимущество по сравнению с известной оценкой Боуэнса-Нотея в смысле ее замыкания до прямых методов, т.е. при стремлении аппроксимации дополнения Шура матрицы Л(А+1) к точному дополнению Шура матрице В этом случае оценка Боуэнса-Потея для числа обусловленности матрицы

М не стремится к единице, поскольку

г = max Н^+^'Л^Н -> тах||5(*,",Л(1*,|| = const > О, в нашем йсе случае из (3.3.24) имеем

Q(i) _ о при А(<> -¥ S{i~l)

и, следовательно,

г = -Д- • £ z{i]T(Q(i) - B{i))xM 0.

X1 Ах

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

Теорема 3.7. Пусть А = D — L — U - блочно-трехдиагоналъная матрица Стильтьесовского типа, а матрицы М, G и G определяются соотношениями

М = (G — L)G~l(G - U),

G = D — 9S — LG^U, Se = (LG~lU - LG^Uje,

G = D-LG~lU.

Тогда для собственных чисел матрицы М~1А справедливо неравенство

1 < Х(М"1А) < ||Gf-1||oo!|Gf||oo(l + гЦДуЦоо + Иад/Иоо), (3.3.38)

где

Rl = (I-LG-1)-1L(G-1-G~1),

Rv = (G_1 - G~l)U(I - G~lU)~l.

При условиях Теоремы 3.7 в случае увеличения "степени неявности" метода, или точности аппроксимации обратной матрицы ¿М, оценка сверху стремится к 1. Таким образом, оценка (3.3.38) отражает замкнутость методов неполной факторизации относительно полного разложения.

ОСНОВНЫЕ РЕЗУЛЬТАТЫ РАБОТЫ

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

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

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

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

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

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

1. Карначук В.И., Ларнн М.Р. Программная реализация метода вложенных сечений. Сб."Вычислительный эксперимент в задачах математической физики", ВЦ СО РАН, Новосибирск, - 1991. - стр. 123-132.

2. Ильин В.П., Ларин М.Р. Об итерационном решеннн разностных уравнений Гельмгольца на сфере. - Препринт ВЦ СО РАН, №966, Новосибирск, - 1992. - 15с.

3. Карначук В.И., Ларин М.Р. Реализация прямых методов для решения систем линейных алгебраических уравнений с разреженной матрицей. - Препринт ВЦ СО РАН, N.971, Новосибирск, -1993. - 46с.

4. Карначук В.И., Ларин М.Р. Кросс-алгоритм для решения систем сеточных уравнений. Сб."Вычислительные технологии", ИВТ СО РАН, Новосибирск, 1993, Том 2, №6, стр. 189-200.

5. Гурьева Я.Л., Карначук В.И., Ларин М.Р. Конвейерный вариант кросс-алгоритма. Сб. "Вычислительные технологии", ИВТ СО РАН, Новосибирск, 1993, Том 2, №6, стр. 201-207.

6. V.P. Il'in and M.R. Larin, One iterative method of solving finite differential equations Helmgoltz on sphere, Rus. J. Numer. Anal. Math. Modelling, 8 (1993), N.4, pp. 297-309.

7. Ильин В.П., Карначук В.И., Ларин М.Р. Древовидный подход к организации структуры данных для разложения Холесского.

ЖВМ II МФ, 1994, Том 34, №12, стр. 1747-1757.

8. Ларин М.Р., Ширин А.В. Параллельный вариант кросс-метода для решения систем сеточных уравнений. Труды ВЦ СО РАН, серия: Вычислительная математика, выпуск 3, Новосибирск, 1995, с. 74-92.

9. Ларин ¡VI.Р. Многоуровневый итерационный метод для решения сеточных систем уравненнй. (в печати сб. трудов ВЦ СО РАН), 1995.

10. М. Larin, Quadri-storage scheme for sparse Cholesky factor and tree-based implementation of Cholesky method, Proceedings of the AMCA conference, Novosibirsk, June 1995, pp. 428-434.