Метод решения

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

Мы будем решать систему (3.1) итерационным предобусловленным методом сопряжённых градиентов (PCG, preconditioned conjugate gradient). Метод применим к симметричным положительно определённым матрицам — именно такие возникают в задачах МКЭ с эллиптическими операторами и в неявных схемах теплопроводности. Распишем алгоритм по шагам.

Пусть x0x_0​=0{} = 0 — начальное приближение. Начальный остаток примет вид

r0\displaystyle r_0 =b\displaystyle {} = b​−A⋅x0.\displaystyle {} - A \cdot x_0.
(3.2)

Начальное направление поиска и вектор предобуславливания имеют вид

P⋅z0\displaystyle P \cdot z_0 =r0,\displaystyle {} = r_0,​p0\displaystyle p_0​=z0,\displaystyle {} = z_0,
(3.3)

где PP — матрица предобуславливателя: она приближает AA, но системы с ней решаются существенно дешевле. Далее на каждой итерации kk​=0,1,2,…{} = 0, 1, 2, \ldots вычисляем

αk\displaystyle \alpha_k =rkT⋅zkpkT⋅A⋅pk,\displaystyle {} = \frac{r_k^T \cdot z_k}{p_k^T \cdot A \cdot p_k},
(3.4)
xk+1\displaystyle x_{k+1} =xk\displaystyle {} = x_k​+αk⋅pk,\displaystyle {} + \alpha_k \cdot p_k,
(3.5)
rk+1\displaystyle r_{k+1} =rk\displaystyle {} = r_k​−αk⋅A⋅pk,\displaystyle {} - \alpha_k \cdot A \cdot p_k,
(3.6)
P⋅zk+1\displaystyle P \cdot z_{k+1} =rk+1,\displaystyle {} = r_{k+1},
(3.7)
βk\displaystyle \beta_k =rk+1T⋅zk+1rkT⋅zk,\displaystyle {} = \frac{r_{k+1}^T \cdot z_{k+1}}{r_k^T \cdot z_k},
(3.8)
pk+1\displaystyle p_{k+1} =zk+1\displaystyle {} = z_{k+1}​+βk⋅pk.\displaystyle {} + \beta_k \cdot p_k.
(3.9)

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

rk+1T⋅rk+1\displaystyle r_{k+1}^T \cdot r_{k+1} <ε.\displaystyle {} < \varepsilon.
(3.10)

Основная вычислительная стоимость каждой итерации — одно умножение разреженной матрицы на вектор A⋅pkA \cdot p_k. Именно поэтому перед началом итераций матрица конвертируется из COO в CSR: в этом формате умножение строки на вектор проходит только по ненулевым элементам.

Схема одной итерации предобусловленного метода сопряжённых градиентов
Рис. 3.2. Одна итерация PCG: шаг A⋅pkA \cdot p_k — самый дорогостоящий.

Далее поговорим о предобуславливателях.