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

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

Мы будем решать систему (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,p0\displaystyle = r_0, \qquad 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 — самый дорогостоящий.

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