Журнал вычислительной математики и математической физики, 2021, T. 61, № 11, стр. 1786-1813

2.1. Алгоритмы для симметричных СЛАУ

Зарождение и первый период развития, который можно отнести к 1950–1965 годам, итерационных методов в подпространствах Крылова, следует связать с именами авторов пионерских работ К. Ланцоша (1950, 1952), В.Е. Арнольди (1951), М.Р. Хестенеса и Е.Л. Штифеля (1951, 1952), Е.Х. Крейга (1954, 1955), Ю.В. Воробьёва (1958), Д.К. Фаддеева и В.Н. Фаддеевой (1958, 1960), В.М. Фридмана (1962), Дж. Голуба (1965) и Р. Варги (1962). Список последующих исследований насчитывается сотнями.

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

Пусть итерационное приближение ${{u}^{n}}$ и соответствующий вектор невязки ${{r}^{n}} = f - A{{u}^{n}}$ вычисляются по следующим формулам:

где ${{\alpha }_{k}}$ и ${{p}^{k}}$ суть итерационные параметры и направляющие векторы. В данном представлении мы опишем семейство алгоритмов сопряженных направлений, характеризуемых различными вариационными и проекционными свойствами, определяемыми соответствующими скалярными произведениями и нормами векторов: где показатели степени будем брать равными $\gamma = 0,\;1,\;2$. Предположим, что направляющие векторы являются ${{A}^{\gamma }}$-ортогональными, т.е. где ${{\delta }_{{n,k}}}$ есть символ Кронекера. Тогда невязки ${{r}^{n}}$ удовлетворяют соотношениям Отсюда следует, что при значении параметров функционалы ${{\Phi }_{\gamma }}$ принимают минимальные значения Здесь и далее индекс $\gamma $ у коэффициентов и векторов для краткости опускается, а матрица $A$ пока предполагается невырожденной.

Для выполнения условия (6) есть два пути. Первый из них состоит в проведении ${{A}^{\gamma }}$-ортогонализации Ланцоша, которую можно записать в следующем виде:

где величины $\mathop {\bar {\beta }}\nolimits_n $ выбираются по условию ${{\left\| p \right\|}_{\gamma }} = 1$. Другой путь (Хестенеса–Штифеля), имеющий наибольшее распространение в силу лучшей устойчивости, заключается в использовании двучленной рекурсии который совместно с формулой для невязки в (4) приводит к другим соотношениям для направляющих векторов: С помощью соотношений (4), (11) можно установить свойство ${{A}^{{\gamma - 1}}}$-ортогональности векторов невязок, а также более удобные (для $\gamma = 1,\;2)$ формулы итерационных коэффициентов: Для случая $\gamma = 0$, однако, вычисление ${{\alpha }_{k}}$ необходимо делать другим образом. Поскольку точное решение СЛАУ может быть представлено в виде разложения по базису то векторы ошибки итерационного приближения и соответствующей невязки записываются в следующей форме: Отсюда с помощью ${{A}^{\gamma }}$-ортогонализации векторов ${{p}^{k}}$ имеем Здесь использована ортогональность ${{p}^{n}}$ векторам ${{p}^{{n - 1}}}$ и ${{p}^{{n - 2}}}$, а также его связь с вектором $A{{p}^{{n - 1}}}$ из (12). Расчет коэффициента ${{\beta }_{{n - 1}}}$ при этом надо проводить по формуле (11).

Рассматриваемые итерационные процессы обладают оптимальными свойствами, обеспечивая минимизацию соответствующих функционалов ${{\Phi }_{\gamma }}({{r}^{n}})$ вида (7) в подпространствах Крылова

Данные алгоритмы при $\gamma = 0,\;1,\;2$ имеют названия методов минимальных итераций, или ошибок (см. [8]), а также сопряженных градиентов и сопряженных невязок соответственнo (см. [5], [9], [27]).

Описанные подходы допускают простое обобщение на СЛАУ, предобусловленные с помощью некоторых СПО матриц $B$. Для сохранения симметричности систем это сделать целесообразно путем двухстороннего предобусловливания с помощью формально вводимой матрицы ${{B}^{{1/2}}}$. В результате система (1) принимает вид

В результате применения формул сопряженных направлений к СЛАУ (17) после некоторых преобразований для $\gamma = 1,2$ получаем следующий итерационный процесс:

где новые векторы связаны со старыми соотношениями $\mathop {\hat {p}}\nolimits^n = {{B}^{{ - 1}}}{{p}^{n}}$, ${{\hat {r}}^{n}} = {{B}^{{ - 1}}}{{r}^{n}}$. Для метода минимальных итераций с $\gamma = 0$ вычисление параметров ${{\alpha }_{{n - 1}}},{{\beta }_{{n - 1}}}$ следует проводить по формулам (15), (11) с заменой величин ${{r}^{k}},{{p}^{k}}$ на ${{\bar {r}}^{k}} = {{B}^{{ - 1}}}{{r}^{k}}$, $\mathop {\bar {p}}\nolimits^k = {{B}^{{ - 1}}}{{p}^{k}}$ соответственно.

Замечание 1. На практике зачастую актуальным является решение серии СЛАУ с одинаковой матрицей, но разными последовательно определяемыми правыми частями. В этом случае при реализации второй и последующих систем можно сэкономить значительный объем вычислений, если после первого расчета запомнить найденные векторы ${{p}^{k}}$ и $A{{p}^{k}}$ и использовать их каким-то образом для расширения крыловского базиса при решении остальных СЛАУ. Этим вопросам посвящено много работ (см. [28] со списком литературы из 157 наименований), однако проблема существенного повышения “информативности” базиса итерационных процессов здесь пока требует дальнейших исследований.

Замечание 2. В книге [8] отмечался такой методический вопрос, что методы в подпространствах Крылова тесно связаны с проблемой моментов, играющей важную роль в теории операторов и многих прикладных задачах. Начала этой тематики заложены П.Л. Чебышевым, А.М. Марковым и Т. Стильтьесом.

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

Такое решение всегда существует и является единственным, а его формальное представление получается с помощью левой трансформации Гаусса исходной системы $({{A}^{{\text{т}}}}Au = {{A}^{{\text{т}}}}f)$ и имеет вид

где ${{A}^{ + }}$ означает псевдообратную, или обобщенно обратную, матрицы $A$.

Для решения таких СЛАУ в работах М.А. Саундерса и соавт. (см. обзор в [29]) разработаны методы SYMMLQ, MINRES, MINRES–QLP, основанные на $A$-ортогонализации Ланцоша, которая в силу рекурсий (10) записывается в матричном виде как

где ${{P}_{k}} = [{{p}^{0}},{{p}^{1}},\; \ldots ,\;{{p}^{{k - 1}}}]$. При этом приближенное решение ищется в форме ${{u}^{k}} = {{P}_{k}}{{q}^{k}}$ и сводится к определению вектора ${{q}^{k}} \in {{\mathcal{R}}^{k}}$ по условию минимума невязки Данная задача реализуется с помощью LQ- или QR-преобразований трехдиагональной матрицы ${{T}_{k}}$ с ортогональной матрицей $Q$ на основе устойчивых операций отражения Хаусхолдера. Приведенные авторами результаты представительных численных экспериментов демонстрируют хорошую эффективность алгоритмов, а их программные реализации доступны в Интернете.

Скорость сходимости итераций в методах крыловского типа для положительно полуопределенного СЛАУ характеризуется оценкой числа итераций $n(\varepsilon )$, необходимой для подавления начальной ошибки в ${{\varepsilon }^{{ - 1}}}$ раз:

где $\kappa $ в общем случае есть эффективное число обусловленности, которое для вырожденных матриц выражается через ${{\lambda }_{{min}}}$ – минимальное ненулевое собственное число.

Для симметричных СЛАУ проблема устойчивого численного решения в основном была решена в XX веке, в том числе для знаконеопределенных, вырожденных и несовместных систем. Что касается классического метода сопряженных градиентов, то здесь особо следует отметить оригинальную теорию сходимости И.Е. Капорина (см. [30]), основанную на введенном им $K$-числе обусловленности, выражаемого через след и определитель матрицы.

Более конкретно рассмотрим предобусловленный алгоритм сопряженных градиентов в следующем виде:

При этом для выполнения условия

в [31] доказывается достаточность проведения числа итераций где $K$ – число обусловленности предобусловленной матрицы ${{B}^{{ - 1}}}A$, определяемое как

В [32] для ряда предобусловливателей проведена оптимизация итерационных алгоритмов в плане минимизации $K$ – числа обусловленности.

2.2. Крыловские алгоритмы для несимметричных СЛАУ

Изложение конкретных подходов мы продолжим с достаточно широкого класса мультипредобусловленных алгоритмов полусопряженных направлений (SCD – Semi-Conjugate Direction, см. [31]). В общей блочной форме такие итерационные методы в подпространствах Крылова записываются следующим образом:

Здесь $p_{1}^{n},\; \ldots ,\;p_{{{{M}_{n}}}}^{n}$ – направляющие векторы, составляющие на $n$-й итерации матрицу ${{P}_{n}}$, а ${{\bar {\alpha }}_{n}}$ – вектор итерационных параметров. Относительно векторов $p_{k}^{n}$ в соотношениях (23) пока предполагается выполнение только условий ортогональности Однако, если при этом коэффициенты ${{\bar {\alpha }}_{n}} = {\text{\{ }}{{\alpha }_{{n,l}}}{\text{\} }}$ определить по формулам то из (23) получаем следующие выражения для функционалов невязки: которые достигают своих минимумов в блочных подпространствах Крылова при $\gamma = 1$, а в случае симметричности матрицы $A$ и для $\gamma = 0$. Здесь следует отметить, что для несимметричных матриц при $\gamma = 0$ применение алгоритмов фактически ограничивается СЛАУ с матрицами, имеющими положительно определенную симметричную часть ${{A}_{s}} = 0.5(A + {{A}^{{\text{т}}}})$, так как для кососимметричной матрицы, например, $(Au,u) = 0$ при $u = {\text{\{ }}1{\text{\} }}$.

Свойства ортогональности направляющих векторов (24) можно обеспечить, если их определить с помощью “мультипредобусловленных” рекуррентных соотношений, в которых каждому вектору $p_{l}^{{n + 1}}$ соответствует “своя” предобусловливающая матрица ${{B}_{{n + 1,l}}}$:

Если матрица $A$ является симметричной, то данные рекурсии превращаются из длинных в двучленные, и мы приходим к методам сопряженных градиентов или сопряженных невязок ($\gamma = 0,\;1$ соответственно, в мультипредобусловленных или в классических вариантах). Для несимметричных СЛАУ эти алгоритмы называются методами полусопряженных градиентов и полусопряженных невязок (SCG или SCR – Semi-Conjugate Gradient или Semi-Graduate Residual). В блочных (мультипредобусловленных) методах сопряженных направлений для решения симметричных СЛАУ формулы (23)–(25) не меняются, а вычисление направляющих матриц ${{P}_{n}}$ выполняется по следующим двучленным рекуррентным соотношениям: Данные формулы, как и (28), содержат “встроенные” в них предобусловливающие матрицы. Если положить ${{B}_{{n + 1,l}}} = I$ и ${{M}_{n}} = 1$ для всех $n$, то получим крыловский процесс “в чистом виде”, без предобусловливания. Отметим, что формулы (28) реализуют алгоритм ортогонализации Грама–Шмидта, для повышения устойчивости которого целесообразно перейти к модифицированному методу ортогонализации (MGS, см. [5], [33]).

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

где $t$ – номер рестарта (далее вычисления до $n = {{n}_{{t + 1}}}$ осуществляются по обычным рекурсиям вида (23)). Оба данных подхода приводят к существенному замедлению итерационного процесса.

Для устранения такого стагнирующего эффекта предлагается добавить второй уровень итераций с применением методов наименьших квадратов (МНК) (см. [34]). Пусть нам известны рестартовые приближения ${{u}^{{{{n}_{0}}}}},{{u}^{{{{n}_{1}}}}},\; \ldots ,\;{{u}^{{{{n}_{t}}}}}$, ${{n}_{0}} = 0$. Тогда для коррекции итерационного вектора ${{u}^{{{{n}_{t}}}}}$, который является начальным для очередного рестартового периода метода (23)–(29), будем использовать следующую линейную комбинацию:

вектор коэффициентов $\bar {b}$ которой определяется по условию минимума нормы $\left\| {{{r}^{{{{u}_{t}}}}}} \right\|$ невязки из обобщенного решения переопределенной алгебраической системы

Решение этой СЛАУ можно получить, например, с помощью $QR$- или $SVD$-разложений матрицы ${{W}_{t}}$. Нормальное решение с минимальной нормой $\left\| {\bar {b}} \right\|$ определяется после применения к (31) левой трансформации Гаусса

Более облегченный формат СЛАУ, в смысле уменьшения ее числа обусловленности, следует после умножения системы (31) слева на матрицу $V_{t}^{{\text{т}}}$: Если матрица ${{V}_{t}}$ имеет полный ранг, то матрицы $A$ и ${{C}_{t}}$ будут невырожденными одновременно. При этом для вектора коррекции из (30) имеем ${{{v}}^{{{{n}_{t}}}}} = {{B}_{t}}{{r}^{{{{n}_{t}}}}} \equiv {{V}_{t}}{{(V_{t}^{{\text{т}}}A{{V}_{t}})}^{{ - 1}}}V_{t}^{{\text{т}}}{{r}^{{{{n}_{t}}}}}$, где матрица ${{B}_{t}} = {{V}_{t}}{{\hat {A}}^{{ - 1}}}V_{t}^{{\text{т}}}$, $\hat {A} = V_{t}^{{\text{т}}}AV$, является малоранговой аппроксимацией матрицы ${{A}^{{ - 1}}}$. В рассматриваемом подходе все рестартовые векторы запоминаются в корректированном виде, а соответствующие невязки вычисляются по формуле ${{r}^{{{{n}_{t}}}}} = f - A{{u}^{{{{n}_{t}}}}}$. Если для какой-то рассматриваемой матрицы обратной не существует, то используется обобщенная обратная матрица. Многочисленные эксперименты с использованием МНК для ускорения крыловских процессов с рестартами показывают его высокую эффективность.

Отметим еще следующую возможность повышения производительности методов SCD с рестартами: при проведении итераций первого рестартового периода запомнить все направляющие векторы ${{p}^{n}}$, а также векторы $A{{p}^{n}}$, и при вычислениях последующих рестартовых периодов новые векторы ${{p}^{n}}$ и $A{{p}^{k}}$ не рассматривать, а использовать старые.

Описанный класс методов полусопряженных направлений с динамическим мультипредобусловливанием по скорости сходимости итераций эквивалентен другим известным алгоритмам решения несимметричных СЛАУ в подпространствах Крылова, среди которых наиболее популярен обобщенный метод минимальных невязок GMRES, основанный на ортогонализации Арнольди и существующий в различных вариантах. Данный алгоритм, предложенный в 1986 г. Й. Саадом и М. Шульцем (см. [5]), получил широкую и заслуженную популярность. Среди многочисленных его исследований отметим результаты Л.А. Книжнермана (см. [35]) по оценкам погрешности ортогонализации Арнольди.

Другой принцип построения крыловских итерационных процессов для решения несимметричных СЛАУ основан на построении последовательностей биортогональных векторов. В этом случае направляющие векторы вычисляются из коротких (двучленных) рекурсий, но на каждой итерации требуется двукратное вычисление векторно-матричного произведения. На идеях биортогонализации были построены алгоритмы BiCG, CGS, BiCGStab в разных модификациях, затем появились их аналоги с $A$-биортогонализацией направляющих векторов (см. обзор в [36]). Примыкающее к ним семейство алгоритмов IDR (Induced Dimension Reduction) (см. [37] и цитируемую там литературу) базируется на подпространствах Сонневельда. Среди других подходов следует рекомендовать методы квазиминимальных невязок QMR и наименьших квадратов (LSQR, LSMR) (см. [38]). В последние годы возродились также алгоритмы, основанные на бидиагональных преобразованиях Крейга–Голуба–Кахана (изначально предложенных в 1955 и в 1965 г.), содержательное исследование, обобщение и сравнительный анализ которых даны в [25], [39]. И наконец, отметим еще “гибкие” методы сопряженных градиентов (FCG) (см. [40]), являющиеся обобщением классических методов CG (в том числе предобусловленных) на несимметричный случай.

3.1. Алгоритмы декомпозиции областей

При решении сеточных СЛАУ, особенно узлового типа, каждая компонента искомого вектора ассоциируется с одним узлом сетки, а структура и портрет матрицы (конфигурация ее ненулевых элементов) изоморфны сеточному графу не направленному для симметричных матриц и направленному – для несимметричных). В этом случае матричная блочная структура вида (2) геометрически наглядно представляется как декомпозиция сеточной расчетной области $\Omega $ на сеточные подобласти ${{\Omega }_{q}}$, в каждой из которых формируется своя подзадача, алгебраически соответствующая диагональному блоку ${{A}_{{q,q}}}$. Исторический приоритет в методах декомпозиции принадлежит теории многомерных краевых задач еще с XIX века, когда решение дифференциальных уравнений в сложных расчетных областях с помощью альтернирующего метода Шварца сводилось к исследованию последовательности вспомогательных задач в пересекающихся подобластях с использованием условий Дирихле на вводимых внутренних границах. Позднее такие подходы были обобщены на непрерывном уровне в работах В.В. Смелова, П.Л. Лионса, А.М. Мацокина, О. Видлунда, Ф. Натафа и соавт., Ю.В. Василевского, М.А. Ольшанского и многих других авторов (см. обзоры литературы в [5], [7], [17]–[19], [52]–[57]), а также перенесены на дискретный, или сеточный, уровень. В эпоху многопроцессорных суперкомпьютеров методы декомпозиции стали главным инструментом для распараллеливания вычислений при решении многомерных начально-краевых задач со сложными реальными данными, требующими высокого разрешения и оперативности расчетов.

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

В целом реализация МДО осуществляется как двухуровневый итерационный процесс в подпространствах Крылова: верхний уровень – это аддитивный блочный метод Шварца–Якоби вида (34), который для простоты в стационарном варианте можно записать как

где ${{B}_{1}}$ есть соответствующий предобусловливатель, который фактически будет использоваться в крыловском процессе. На каждой $n$-й итерации (37) осуществляется синхронное, или параллельное, решение вспомогательных подзадач в подобластях ${{\Omega }_{q}}$, которое формально заключается в обращении блоков предобусловливателя ${{B}_{{q,q}}}$. Это может быть сделано или прямыми, или итерационными методами, среди которых для больших алгебраических систем естественно выбрать предобусловленные алгоритмы в соответствующих подпространствах Крылова. Отметим, что в последнем случае мы приходим к задаче многократного решения СЛАУ с одинаковыми матрицами, но разными правыми частями, что можно сделать экономично путем переиспользования базисов крыловских подпространств, как это указывалось выше в замечании 1.

С точки зрения распараллеливания вычислений на МВС, предпочтительнее разбивать сеточную расчетную область на большее число подобластей ${{\Omega }_{q}}$, $q = 1,\;2,\; \ldots ,\;P$, чтобы все соответствующие подсистемы решать синхронно на $P$ процессорах. Однако при этом будет очевидным образом расти количество внешних итераций, которые можно оценить величиной $O({{P}^{{1/d}}})$, где $d$ есть геометрическая размерность области (в предположении, что декомпозиция формирует $d$‑мерную “макросеть” из подобластей). Число внешних итераций можно сократить за счет увеличения размеров пересечения контактирующих подобластей (over-lapping), но тогда реализация каждого шага будет становиться более трудоемкой (на практике пересечение сеточных подобластей выбирается размером в несколько шагов $h$).

Для ускорения итераций крыловского типа в МДО существуют различные подходы, которые могут быть интерпретированы как применение некоторых дополнительных предобусловливателей, в том числе сглаживающего типа, в силу чего такие методы иногда называются гибридными (см. [58], [59]). Рассмотрим в качестве примера метод грубосеточной коррекции, основанный на принципе дефляции для выбора начального приближения ${{u}^{0}}$.

Пусть у нас имеется какой-то информационный базис из векторов ${{v}_{1}},\; \ldots ,\;{{v}_{m}}$, составляющих матрицу $V = ({{{v}}_{1}}\; \ldots \;{{{v}}_{m}}) \in {{\mathcal{R}}^{{N,m}}}$, $m < N$. Выбирая какой-то произвольный вектор ${{u}^{{ - 1}}}$, будем искать начальное приближение и соответствующий вектор невязки в виде

где вектор неизвестных коэффициентов $c = {{({{c}_{1}},\; \ldots ,\;{{c}_{m}})}^{{\text{т}}}}$ определяется по условию минимума невязки ${{r}^{0}}$ из переопределенного уравнения Умножая обе части (38) слева на ${{V}^{{\text{т}}}}{{A}^{\gamma }}$ с параметром $\gamma = 0$ или $\gamma = 1$, получаем семейство обобщенных решений где случай $\gamma = 1$ соответствует нормальному решению, т.е. минимизирующему евклидовые нормы невязки ${{r}^{0}} = {{r}^{{ - 1}}} - AVc$ и самого решения $c$. Отсюда для откорректированной невязки получаем выражения которые соответствуют минимуму функционала Нетрудно видеть, что при этом векторы ${{{v}}_{s}}$ играют роль направляющих векторов в подпространствах Крылова, т.е. они обеспечивают вектор невязки ${{r}^{0}}$ ортогональными и вариационными свойствами, аналогичными тем, которые имеются в методах сопряженных градиентов и сопряженных невязок при $\gamma = 0$ и $\gamma = 1$ соответственно.

Следуя этим соображениям, далее для методов сопряженных направлений определим через ${{r}^{0}}$ начальный направляющий вектор из соответствующих ортогональных условий:

Введенная здесь матрица ${{B}_{2}}$ играет роль второго предобусловлителя, в дополнение к ${{B}_{1}}$ из (37), и может вместе с ним использоваться на последующих итерациях по формулам мульти-предобусловленных сопряженных направлений (см. [17]). Отметим, что в (40) степени матрицы $A$ подобраны таким образом, чтобы обеспечить симметричность ${{B}_{2}}$.

Вопрос выбора базисных векторов ${{{v}}_{1}},\; \ldots ,\;{{{v}}_{m}}$, эффективно дополняющих подпространство Крылова, актуален во многих ускоряющих итерационных процедурах. Так, алгоритм грубосеточной коррекции, активно используемый в МДО, основан в простейшем случае на “геометрической” декомпозиции сеточной области ${{\Omega }_{q}}$, а его компоненты равны 1 или 0 в зависимости от того, принадлежит ли соответствующий узел $q$-й подобласти или нет. Дальнейшее развитие этого подхода заключается в переходе от кусочно-постоянных базисных функций ${{{v}}_{q}}$ к кусочно-полиномиальным более высоких порядков (см. [53]).

3.2. Многосеточные методы

Изначально в пионерских работах Р.П. Федоренко и Н.С. Бахвалова многосеточные подходы базировались на спектрально-аппроксимационных принципах, с раздельным подавлением компонент ошибки, соответствующих низким и высоким частотам. Дальнейшее их развитие осуществлялось в геометрическом плане, в алгебраическом (наиболее популярное название – AMG, от Algebraic Multi-Grid), а также в комбинаторном (см. [5], [56], [60]–[67]).

Данные подходы занимают исключительное место в вычислительной алгебре, поскольку являются асимптотически оптимальными по порядку, т.е. требуют объема компьютерных ресурсов, пропорционального размерности СЛАУ. Наиболее общая алгебраическая методология базируется на универсальном принципе предобусловливания матриц. А реализация технологий с различным количеством вложенных сеток трактуется как рекурсивное применение двухсеточного алгоритма. Для простоты рассмотрим последовательность вложенных сеток ${{\Omega }_{m}} \subset \; \ldots \; \subset {{\Omega }_{2}} \subset {{\Omega }_{1}}$, т.е. узлы более редкой сетки являются узлами “ближайшей” густой сетки ${{\Omega }_{k}}$. Будем предполагать, что сетка ${{\Omega }_{k}}$ получается из ${{\Omega }_{{k + 1}}}$ путем “сгущения вдвое”, так что соответствующие количества узлов соотносятся как ${{N}_{k}} \approx {{2}^{d}}{{N}_{{k + 1}}}$, где $d$ есть размерность области. Для решения исходной СЛАУ $Au = f$ на $\Omega $, которую можно переобозначить как ${{A}_{1}}{{u}_{1}} = {{f}_{1}}$, строится неявно предобусловливатель ${{B}_{1}}$ с помощью приближенного решения системы со специально конструируемой матрицей ${{A}_{2}}$ для сетки ${{\Omega }_{2}}$ и т.д., процесс “зацикливается”.

В современной трактовке достаточно общего вида, многосеточный метод AMG может быть представлен следующими этапами итерационного решения алгебраической системы ${{A}_{k}}{{u}_{k}} = {{f}_{k}}$ на ${{\Omega }_{k}}$.

1. На “густой” сетке ${{\Omega }_{k}}$ по заданному приближению для вектора $u_{k}^{n}$ вычисляется невязка

2. Для вектора $r_{k}^{n}$ производится препроцессинг (предварительное сглаживание), как правило, путем проведения нескольких итераций с помощью какого-либо простого алгоритма. Более конкретно, данная процедура реализуется за два шага. На первом вычисляется направляющий вектор

где ${{\hat {S}}_{k}}$ – оператор, или матрица, этого этапа (рresmoothing). Фактически ${{\hat {S}}_{k}}$ есть некоторая аппроксимация матрицы ${{A}_{k}}$ (формальное представление результатов нескольких циклов применяемого итерационного метода сглаживания). На втором шаге определяется очередное (сглаженное) приближенное решение и соответствующая невязка:

3. По определенному на густой сетке ${{\Omega }_{k}}$ вектору $\tilde {r}_{k}^{n}$ формируется вектор невязки $\tilde {r}_{{k + 1}}^{n}$ на редкой сетке ${{\Omega }_{{k + 1}}}$:

где ${{R}_{k}}$ есть некоторый оператор сужения, или ограничения (этап restriction). В простейшем случае данная операция есть проектирование вектора с сетки ${{\Omega }_{k}}$ на ${{\Omega }_{{k + 1}}}$.

4. На редкой сетке вычисляется направляющий вектор $\hat {p}_{{k + 1}}^{n}$ из решения СЛАУ

где ${{A}_{{k + 1}}}$ – матрица СЛАУ для сетки ${{\Omega }_{{k + 1}}}$, в некотором смысле наследующая аппроксимационные и алгебраические свойства оператора ${{A}_{k}}$ на сетке ${{\Omega }_{k}}$. Существуют различные приемы для построения матрицы ${{A}_{k}}$, которая иногда называется матрицей грубосеточной коррекции.

5. Найденный из решения системы (45) вектор $\hat {p}_{{k + 1}}^{n}$ продолжается с редкой сетки $\Omega $ на густую сетку ${{\Omega }_{k}}$ (этап prolongation):

В некотором смысле согласованным определением вводимых операторов является такое, которое удовлетворяет соотношению $A_{k}^{{ - 1}} \approx {{P}_{k}}A_{{k + 1}}^{{ - 1}}{{R}_{k}}$. Если, например, матрица ${{A}_{k}}$ симметрична, то для наследования этого свойства на редкой сетке ${{\Omega }_{{k + 1}}}$ целесообразно выбирать ${{P}_{k}} = R_{k}^{{\text{т}}}$. В частности, при этом определяется так называемая Галёркинская аппроксимация ${{A}_{{k + 1}}} = P_{k}^{{\text{т}}}{{A}_{k}}{{P}_{k}}$.

6. Для вектора определяются соответствующиe векторы приближенного решения и невязки на густой сетке (resudial update)

7. Проводится “пост-процессинг”, т.е. повторное сглаживание, для вновь полученных направляющего вектора и итерационного приближения решения на густой сетке ${{\Omega }_{k}}$, с одновременным вычислением нового направляющего вектора $\hat {p}_{k}^{n}$ из решения вспомогательной СЛАУ с матрицей (оператор повторного сглаживания, этап postsmoothing):

аналогично (42).

8. Итоговый направляющий вектор определяется как результат первой итерации AMG:

где ${{B}_{k}}$ есть предобусловливающая матрица данной двухсеточной стадии метода, явный вид которой можно выразить через операторы сглаживания, ограничения, грубосеточной коррекции и продолжения, однако ее конкретное представление не суть важнo, так как общая схема вычислений и ее реализация описываются формулами (41)–(48). Частных вариантов двухсеточного метода может быть достаточно много, и все они (при организации стационарного итерационного процесса) имеют вид представляющий собой двухуровневый итерационный процесс, поскольку операция умножения на $B_{k}^{{ - 1}}$ включает решение СЛАУ с матрицей ${{A}_{{k + 1}}}$, что при больших значениях ${{N}_{{k + 1}}}$ прямым алгоритмом делать нeцелесообразно. Понятно, что на основе алгоритма (49) можно построить предобусловленный метод в подпространствах Крылова.

Но идея многосеточного подхода состоит в использовании рекурсивного принципа: для приближенного решения уравнений на редкой сетке ${{\Omega }_{{k + 1}}}$ применяются приведенные выше приемы с использованием еще более редкой сетки ${{\Omega }_{{k + 2}}}$ и т.д. При этом на каждой сетке ${{\Omega }_{k}}$ весь цикл вычислений по формулам (41)–(48) можно повторять заданное число $\gamma $ раз, что может быть записано как обобщение соотношения (49):

На практике число таких повторов выбирается равным $\gamma = 1$ или $\gamma = 2$. В первом случае получаемая вычислительная схема называется $V$-циклом, а во втором – $W$-циклом. В [64] и [65] предложено развитие такого подхода: $K$-цикл, в котором при $\gamma = 2$ используется не стационарный итерационный процесс, а крыловского типа (более конкретно, “гибкий” метод сопряженных градиентов FCG из [40]). В последней формуле надо иметь в виду, что предобусловливатель ${{B}_{k}}$ выражается рекурсивно для $k = 1,\;2,\; \ldots ,\;m - 1$, а на последней $m$-й сетке СЛАУ вида (45) с матрицей ${{A}_{{m + 1}}}$ решается непосредственно (прямым или итерационным методом).

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

3.3. Алгоритмы неполной факторизации

Наиболее классическими приемами предобусловливания СЛАУ являются многочисленные явные и неявные методы приближенной факторизации матриц, обширный анализ которых без аспектов распараллеливания приведен в [10]. Современные подходы к этим технологиям, включая перспективные аппроксимации матриц, обратных к треугольным множителям из неполных матричных разложений, приведены в [68]–[70] и цитируемой там литературе.

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

Здесь $L$ и $U$ – нижняя и верхняя треугольные части исходной матрицы $A = D + L + U$, а $D$ – некоторая блочно-диагональная или диагональная в частном случае матрица. Из общего требования $B \approx A$ блочно-диагональная матрица $G$ строится с помощью соотношения где черта над матрицей обозначает ее какую-то ленточную аппроксимацию.

На основе формул (51), (52) конструируются различные методы симметричной последовательной верхней релаксации (при этом $G = 2{{\omega }^{{ - 1}}}D)$, где $\omega $ есть итерационный параметр, а также явной и неявной неполной факторизации, содержательный анализ которых представлен в книгах [2], [5], [10], [68], а также в специальных обзорных статьях [51], [69], [70]. Широкое семейство алгоритмов основано на неполном LU-разложении исходной матрицы, которое в симметричном случае сводится к приближенному разложению Холеского. Качество аппроксимации матрицы и скорость сходимости повышаются за счет “уплотнения” матричных множителей, но при этом “дорожает” каждая итерация.

Схожий подход, но ориентированный на распараллеливание, заключается в конструировании разреженной аппроксимации обратной матрицы, устраняющей необходимость решения треугольных систем на каждой итерации. Описанные методы исследованы как в поточечном, так и в блочных вариантах. Оригинальные результаты в этих направлениях получены А.Ю. Ереминым, И.Е. Капориным, Л.Ю. Колотилиной, И.Н. Коньшиным и другими авторами в работах [31], [32], [71], [72]. Другое перспективное направление основано на факторизации, в том числе приближенной, иерархических матриц в HSS-форматах с использованием малоранговых блочных структур в LU-разложении (см. [73], [74]).

Наконец, отметим еще перспективное семейство алгоритмов построения предобусловливателей, базирующееся на сетевом программировании или на методах преобразования графов (см. [75]–[77]). Существенным является и такой момент, что если традиционно неполное треугольное разложение матриц реализуется для невырожденных матриц, то сейчас имеются его модификации и для положительно-полуопределенных случаев.

Один из общих подходов к ускорению итераций заключается в принципе компенсации, или согласования строчных сумм, который заключается в подборе матрицы $G$ в (51), (52) по выполнению условий

на некотором заданном наборе из $m$ пробных векторов (см. [49], [78] и цитируемую там литературу). В некоторых работах такой прием также называется фильтрацией.

Для удовлетворения условия (53) матрица $G$ преобразуется к виду

где $\theta \in \left[ {0,\theta \in \left[ {0,1} \right]} \right],$ – некоторый компенсирующий параметр, а $S$ – блочно-диагональная матрица, формируемая по условию

Среди методов неполной факторизации для решения сеточных СЛАУ наиболее быстрыми являются неявные, которые базируются на алгоритмах прогонки решения вспомогательных трехдиагональных систем. Такие алгоритмы эффективно распараллеливаются на МВС с помощью многопотоковых реализаций рекурсивных прогоночных процедур, основанных или на методах редукции исходных СЛАУ, или на методах окаймления. Отметим, что в [57] предложен двухпотоковый блочный вариант попеременно-треугольной факторизации исходной матрицы, когда каждый из множителей не является нижней или верхней треугольной матрицей, а состоит из блочных строк различной ориентации: одни являются нижними треугольными, а остальные – верхними треугольными (такое разложение названо авторами “twisted decomposition”, а в русскоязычной литературе данный подход традиционно определяется как алгоритм встречных прогонок, см. обзор в [27]).

Если матрица $A$ симметричная, т.е. $D = {{D}^{{\text{т}}}}$ и $L = {{U}^{{\text{т}}}}$, то естественно выбирать матрицы $G$ и $B$ также симметричными, а для исходной СЛАУ – использовать двустороннее предобусловливание с сохранением симметрии итоговой системы:

Предобусловленная матрица $\bar {A}$ после несложных преобразований может быть приведена к виду

которое допускает реализацию векторно-матричного умножения в виде во многих случаях дающем существенную экономию вычислений. Отметим, что матричные представления (56)–(58) справедливы и без требований ${{L}_{B}} = U_{B}^{{\text{т}}}$, ${{L}_{G}} = U_{G}^{{\text{т}}}$, т.е. данные формулы приемлемы и для несимметричных СЛАУ.

Непосредственное обращение треугольных матриц, которое присутствует в описанных выше алгоритмах, плохо распараллеливается. Существуют специальные технологические приемы, позволяющие в определенной степени этот недостаток устранить (см., например, [79] и обзор литературы в [80]). Мы рассмотрим другой подход, использующий вместо $L$, $U$ переменно-треугольные матрицы $\hat {L}$, $\hat {U}$, имеющиe на части строк ненулевые элементы только слева от главной диагонали, а на остальных – только справа. Построение таких преобусловливателей мы проиллюстрируем на примере блочно-трехдиагональной матрицы

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

Пoскольку в данном случае $A = D + \hat {L} + \hat {U}$, формулы неполной факторизации (51)–(58) для определения предобусловливателя $B$ остаются в силе. Надо только заменить $L$ и $U$ на $\hat {L}$ и $\hat {U}$ соответственно. При этом трехдиагональные блоки матрицы $\bar {G} = \{ {{\bar {G}}_{i}}\} $ находятся по формулам, аналогичным (28), (29), но только с помощью встречных рекурсий (блочных прогонок, см. [25]), которые мы выпишем для произвольного нечетного ${{N}_{x}} = 2m + 1$:

где ${{S}_{i}}$ суть диагональные матрицы для $l = 1$ и трехдиагональные при $l = 2$.

При использовании матрицы $\bar {B}$ в качестве предобусловливателя для какого-либо итерационного процесса на каждом шаге требуется решать вспомогательную систему вида

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

В рекурсиях (61) и (62) этапы вычислений с увеличением индекса $i$ и с его уменьшением можно назвать прямыми и обратными прогонками соответственно. Очевидно, что их реализации могут выполняться параллельно на двух потоках. Данный ограниченный параллелизм основан на односменной триангуляции вида (60) для матриц $L$, $U$. Данный подход можно обобщить на $P$ – сменную триангуляцию, т.е. определить матрицы $\hat {L}$, $\hat {U}$ с условием $\hat {L} + \hat {U} = L + U$ как состоящие из $P$ блочных строк, каждая из которых последовательно меняет тип триангуляции. Формулы блочных прогонок вида (61), (62) принимают более сложный вид, но могут распараллеливаться на $P$ процессорах.

Матрица в (59) соответствует, например, пятиточечным уравнениям Лапласа на прямоугольной сетке с числом узлов $N = {{N}_{x}} \cdot {{N}_{y}}$. В этом случае ${{D}_{i}}$ – трехдиагональные, а ${{L}_{i}}$ и ${{U}_{i}}$ – диагональные блоки. При переходе к трехмерным задачам на параллелепипедоидальной сетке с количеством узлов $N = {{N}_{x}} \cdot {{N}_{y}} \cdot {{N}_{z}}$ матрица $A$ может быть также представлена в блочно-трехдиагональном виде (59), но каждый диагональный блок ${{D}_{k}}$ будет теперь иметь размерность $N = {{N}_{x}} \cdot {{N}_{y}}$ и представлять собой матрицу, структурно совпадающую с (59).

Развитие рассмотренных выше подходов для решения СЛАУ на трехмерных сетках получило название методов вложенных факторизаций (nested factorization), предложенных в 1983 г. Дж. Апплеярдом и И. Чешире, а впоследствии исследованных многими авторами (см. [57], [82], [83] и цитируемую там литературу). Пусть матрица $A$ представлена в виде

где $D$ – диагональная, а ${{L}_{k}}$ и ${{U}_{k}},k = 1,2,3$, – нижние и верхние треугольные матрицы. Определим предобусловливающую матрицу $B$ путем рекурсивной факторизации вида (51): в результате чего получаем В зависимости от способа определения матрицы $M$ могут быть сформированы различные предобусловливатели $B$. Мы остановимся на простом случае, когда матрица (63) соответствует семиточечной аппроксимации задачи Дирихле для уравнения Пуассона в параллелепипеде на параллелепипедной сетке. Тогда при естественной упорядоченности узлов матрицы $M$, $T$ и $Р$ можно сформировать диагональной, трехдиагональной и пятидиагональной соответственно, а предобусловливатель можно определить по формулам Здесь ${{\theta }_{2}}$ и ${{\theta }_{3}}$ – итерационные (релаксирующие) параметры, а ${{S}_{2}}$ и ${{S}_{3}}$ – диагональные матрицы, определяемые из условий согласования строчных сумм, аналогично рассмотренному выше: Отметим, что вместо (67) можно использовать согласование нe строчных, а столбцовых сумм (см. [83] и приведенную там литературу):

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

В данном случае ${{D}_{3}}$ есть блочно-диагональная матрица с пятидиагональными блоками ${{D}_{{3,i}}}$ размерности ${{N}_{y}}{{N}_{z}}$, каждый из которых соответствует плоской задаче в сечении $x = {\text{const}}$ и имеет структуру матрицы $A$ в (59). Тогда матрицу $P$ в (64) определяем формулой что соответствует определению ${{L}_{3}}{{P}^{{ - 1}}}{{U}_{3}} = 0$ в (64). Отметим, что ${{L}_{3}}$ и ${{U}_{3}}$ можно определить как переменно-треугольные, и тогда реализацию алгоритма можно распараллеливать с помощью встречных блочных прогонок.

3.4. Методы решения седловых СЛАУ

Рассмотрим СЛАУ с матрицей седлового типа

где $A \in {{\mathcal{R}}^{{N,N}}}$, $D \in {{\mathcal{R}}^{{{{N}_{1}},{{N}_{1}}}}}$, $C \in {{\mathcal{R}}^{{{{N}_{2}},{{N}_{1}}}}}$, $N = {{N}_{1}} + {{N}_{2}}$ и $f \in {{\mathcal{R}}^{{{{N}_{1}}}}}$ и $p,g \in {{\mathcal{R}}^{{{{N}_{2}}}}}$. Матрица $D$ предполагается симметричной и положительно-полуопределенной. При этом необходимым и достаточным условием невырожденности $A$ является $ker(D) \cap ker(C) = {\text{\{ }}0{\text{\} }}$, а достаточным условием является положительная определенность матрицы $D$ на $ker(C)$ и $ker({{C}^{{\text{т}}}}) = {\text{\{ }}0{\text{\} }}$, что мы будем предполагать выполненным в дальнейшем.

Задачи с седловой точкой представляют собой широко распространенную форму математических постановок в проблемах моделирования, включая начально-краевые смешанные формулировки для дифференциальных классических или обобщенных уравнений, оптимизационные задачи и вычислительную алгебру (см. [84]–[87] и цитируемую там литературу). Мы делаем акцент на методах решения седловых СЛАУ с большими разреженными матрицами, возникающими при сеточных аппроксимациях многомерных краевых задач, которые появляются во многих актуальных приложениях: электромагнетизме, гидро-газодинамике, упруго-пластичности, многофазной фильтрации в пористых средах, в задачах оптимизации и т.д.

В силу особенностей блочной структуры седловых алгебраических систем, методам их решения посвящено значительное количество работ: обзоры Дж. Голуба и соавт. [84], Ф. Брецци [85], П. Василевского [72] и И. Нотея [86], монография Ю.В. Быченкова и Е.В. Чижонкова [87], цикл статей Ч. Грейфа и соавт. (в том числе для несимметричных СЛАУ седлового типа, см. [68], [85]) и М. Ариоли с соавт. [25], [75] (см. также обзоры в [90]–[92]).

Отметим, что без потери общности мы можем рассматривать седловую СЛАУ в виде

Действительно, если мы возьмем какое-либо частное решение системы $C\hat {u} = g$, то вектор $u = \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u} + \hat {u}$, являющийся решением СЛАУ (70), будет удовлетворять системе Заметим также, что любые решения СЛАУ (71) удовлетворяют одновременно системе где ${v},\tilde {f} \in {{\mathcal{R}}^{N}}$, а ${{K}_{0}} \in {{\mathcal{R}}^{{{{N}_{1}},{{N}_{1}}}}}$ – произвольная невырожденная матрица. Поскольку последняя система формально является регуляризацией, или обобщением, СЛАУ (71), мы в дальнейшем останавливаемся на алгоритмах решения именно уравнения (72).

С использованием дополнения Шура

матрица системы (72) факторизуется в виде

Если матрицы $\bar {D}$, $S$ заменить их приближениями (предобусловливателями) ${{B}_{d}}$ и ${{B}_{s}}$, то получаем предобусловливатель матрицы $B$ в виде неполной блочной факторизации

Несколько более грубая аппроксимация, при использовании в (73) только одного (левого или правого) множителя приводит к неполному блочно-треугольному предобусловливанию с матрицей

или к неполному предобусловливателю Узавы

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

Здесь $Q_{d}^{{(1)}}$, $Q_{d}^{{(2)}}$ – некоторые аппроксимации матрицы, обратной или обобщенно обратной к предобусловливателю ${{B}_{d}}$. В частности, если $Q_{d}^{{(1)}} = B_{d}^{{ - 1}}$, $Q_{d}^{{(2)}} = 0$ или $Q_{d}^{{(1)}} = 0$, $Q_{d}^{{(2)}} = B_{d}^{{ - 1}}$, то из (76) мы получаем или алгоритм Узавы с предобусловливателем ${{B}_{3}}$ из (75) (при этом третья стадия опускается, т.е. $u_{d}^{{n + 1}} = \hat {u}_{d}^{{n + 1}}$), или неполное блочно-треугольное предобусловливание с матрицей ${{B}_{2}}$ из (74) (при этом первая стадия в (76) опускается и $u_{d}^{{n + 1}} = u_{d}^{n}$).

Если матрицa ${{Q}_{d}} = Q_{d}^{{(1)}} + Q_{d}^{{(2)}} - Q_{d}^{{(2)}}\bar {D}Q_{d}^{{(1)}}$ является невырожденной, то итерационному процессу (76) соответствует предобусловливатель

В этом случае для $Q_{d}^{{(1)}} = Q_{d}^{{(2)}} = B_{d}^{{ - 1}}$ из (77) следует так называемый симметризованный метод Узавы с предобусловливающей матрицей

Одним из перспективных блочно-диагональных предобусловливателей для решения СЛАУ (72) является следующий:

где ${{K}_{1}}$ и ${{K}_{2}}$ – некоторые симметричные невырожденные матрицы.

Cобственные числа и векторы предобусловленной матрицы $\bar {A} = B_{6}^{{ - 1}}\tilde {A}$ из “возмущенной” СЛАУ (72) определяются спектральной задачей

имеющей в компонентном виде следующее представление:

Анализ данной спектральной задачи позволяет получать интересные результаты для различных частных случаев. В [81], [82] показано, в частности, что уникальные спектральные свойства матрицы $\bar {A}$ получаются при сильной вырожденности блока $\tilde {D}$, актуальной для алгоритмов решения систем уравнений Максвелла. Например, справедливы следующие результаты.

Утверждение 1. Пусть ${{B}_{6}}$ есть СПО-матрица, а ${\text{\{ }}{{z}_{i}},\;i = 1,\;2,\; \ldots ,\;{{N}_{1}} - {{N}_{2}}{\text{\} }}$ – базис ядра матрицы $C$. Тогда предобусловленная матрица $B_{6}^{{ - 1}}\tilde {A}$ имеет ${{N}_{1}} - {{N}_{2}}$ линейно независимых векторов $(z_{i}^{{\text{т}}}0) \in {{\mathcal{R}}^{N}}$, которым соответствуют ${{N}_{1}} - {{N}_{2}}$ кратных собственных значений $\lambda = 1$.

Утверждение 2. Пусть ${{K}_{1}} = {{K}_{2}} = K$, а $\tilde {D}$ – симметричная положительно-полуопределенная матрица с размерностью ядра, равной $r$. Тогда $\lambda = 1$ есть собственное значение матрицы $B_{6}^{{ - 1}}\tilde {A}$ с кратностью ${{N}_{1}}$, а $\lambda = - 1$ – собственное значение кратности $r$. Остальные ${{N}_{2}} - r$ собственных значений принадлежат интервалу $\lambda \in ( - 1,\;0)$ и удовлетворяют соотношению

где $\mu $ суть ${{N}_{2}} - r$ положительных обобщенных собственных значенийПусть ${\text{\{ }}{{z}_{i}},\;i = 1,\;2,\; \ldots ,\;{{N}_{1}} - {{N}_{2}}{\text{\} }}$ есть базис ядра $C$, ${\text{\{ }}{{x}_{i}},\;i = 1,\;2,\; \ldots ,\;{{N}_{2}} - 2{\text{\} }}$ – базис ядра $\tilde {D}$, а ${\text{\{ }}{{y}_{i}},\;i = 1,\;2,\; \ldots ,\;{{N}_{2}} - r{\text{\} }}$ – совокупность линейно независимых векторов, дополняющих $ker(\tilde {D}) \cup ker(C)$ в базиcе ${{\mathcal{R}}^{{{{N}_{1}}}}}$. Тогда ${{N}_{1}} - {{N}_{2}}$ векторов $({{z}_{i}},0)$, $r$ векторов (${{x}_{i}}$, ${{K}^{{ - 1}}}C{{x}_{i}}$) и ${{N}_{2}}$ – $r$ векторов (${{y}_{i}}$, ${{K}^{{ - 1}}}{{y}_{i}}$) суть линейно независимые собственные векторы, соответствующие собственным значениям $\lambda = 1$, а $r$ векторов (${{x}_{i}}$, ${{K}^{{ - 1}}}{{x}_{i}}$) – собственные векторы, соответствующие $\lambda = - 1$.

В целом наличие в предобусловливателе ${{B}_{6}}$ различных матриц ${{K}_{0}}$, ${{K}_{1}}$, ${{K}_{2}}$ предоставляет широкие возможности для конструирования эффективных алгоритмов в конкретных случаях.

Мы рассмотрим далее семейство итерационных методов для решения седловых симметричных СЛАУ с матрицей $\tilde {A}$ из (72), основанное на эффективном подходе $G - K$-бидиагонализации Голуба–Кахана, которая изначально была предложена для сингулярного разложения прямоугольных матриц, но затем в работах М. Саундерса, М. Ариоли, Ч. Грейфа и других авторов успешно применялась для решения алгебраических систем, в том числе с учетом блочной седловой структуры.

Без ограничения общности исследуемую СЛАУ запишем в виде

Легко проверить, что если в (72) функцию $u$ заменить на $u + {{\tilde {D}}^{{ - 1}}}f$, то эта система примет форму (80) с правой частью $g = - C{{\tilde {D}}^{{ - 1}}}f$. Предполагается, что в (80) $\tilde {D}$ и $K$ являются СПО-матрицами, а также выполняется неравенство ${{N}_{1}} \geqslant {{N}_{2}}$.

Метод $G - K$-бидиагонализации базируется на построении $\tilde {D}$-ортогональных векторов ${{v}_{k}}$ и $K$-ортогональных векторов ${{q}_{k}}$, которые удовлетворяют условиям

где $V = [{{{v}}_{1}},\; \ldots ,\;{{{v}}_{{{{N}_{1}}}}}] \in {{\mathcal{R}}^{{{{N}_{1}},{{N}_{2}}}}}$, $Q = [{{q}_{1}},\; \ldots ,\;{{q}_{{{{N}_{1}}}}}]$, а $B \in {{\mathcal{R}}^{{{{N}_{2}},{{N}_{2}}}}}$ есть двухдиагональная матрица Вводя новые неизвестные функции и умножая систему (80) на блочно-диагональную матрицу block-diag(${{V}^{{\text{т}}}},{{Q}^{{\text{т}}}}$), получаем

Из (83) cледует, что вектор $u$ зависит от ${{N}_{2}}$ столбцов матрицы $V$, поскольку ${{z}_{2}} = 0$. Таким образом, СЛАУ (83) редуцируется к виду

Отсюда, полагая находим начальный вектор ${{{v}}_{1}}$: Дальнейшие векторы ${{{v}}_{n}}$, ${{q}_{i}}$ и элементы ${{\alpha }_{i}}$, ${{\beta }_{i}}$ матрицы $B$ вычисляются из рекуррентных соотношений, $n = 1,\;2,\; \ldots $: Последовательные приближения ${{u}^{n}}$ в соответствии с (82) определяются по первым $n$ столбцам матрицы $V$, т.е. где $z_{1}^{j}$ – компоненты вектора ${{z}_{1}}$ из (84), вычисляемые по формулам Опуская детали вывода формул (см. подробнее [87]), мы приведем итоговые рекуррентные соотношения для итерационного решения: где ${{d}_{n}}$ есть $n$-й столбец $D = Q{{B}^{{ - 1}}}$. Отметим, что описанный алгоритм обладает следующими свойствами оптимизации: на каждом $n$-м шаге ошибка итерационного приближения достигает своего минимума: Как отмечается в [71], данный метод обладает высокой скоростью сходимости при решении СЛАУ, возникающих при конечно-элементных аппроксимациях непрерывных многомерных краевых задач седлового типа, т.е. в смешанных постановках. При этом существенным фактором является то, что на каждой итерации требуется решать возникающие алгебраические подсистемы с матрицами $\tilde {D}$ и $K$, которые в некотором смысле играют роль предобусловливателей. Их приближенное обращение приводит фактически к двухуровневым итерационным процессам в определенных подпространствах.