Cholesky 预处理共轭梯度法¶

约 174 个字 60 行代码预计阅读时间 1 分钟

Cholesky预处理： 预处理矩阵 $\mathbf{M} = LL^T$，其中 $L$ 为 $A$ 的下三角Cholesky分解

每步的PCG主要分量公式如下：

预处理向量 $z$ 的计算
（两步三角系统求解，逐分量写为）：
先解 $L y = r$： $$ y_1 = \frac{r_1}{L_{11}} $$ $$ y_i = \frac{1}{L_{ii}}\left(r_i - \sum_{j=1}^{i-1} L_{ij} y_j\right),\quad i=2,\ldots,n $$
再解 $L^T z = y$： $$ z_n = \frac{y_n}{L_{nn}} $$ $$ z_i = \frac{1}{L_{ii}}\left(y_i - \sum_{j=i+1}^{n} L_{ji} z_j\right),\quad i=n-1,\ldots,1 $$
其余 $\alpha_k$、$x^{(k+1)}$、$r^{(k+1)}$、$\beta_k$、$p^{(k+1)}$ 的表达方式 完全和Jacobi预处理一致，只是预处理向量$z$的获得方式不同。

subroutine pcg_c(a, b, x, x0, n)
    ! Cholesky预处理共轭梯度法（适用于对称正定A）
    integer, intent(in) :: n
    real*8, intent(in) :: a(n,n), b(n), x0(n)
    real*8, intent(out) :: x(n)

    real*8 :: x1(n), x2(n), r0(n), r1(n)
    real*8 :: z0(n), z1(n), y(n)
    real*8 :: p0(n), p1(n)
    real*8 :: L(n, n)
    integer :: i, k, iter_max
    real*8 :: eps, tmp1, tmp2, tmp3, afa, beta, dr_s

    iter_max = 1000
    eps = 1.0d-7

    write(*,501)

    ! --- 1. Cholesky分解 A = LL^T ---
    call cholesky(a, L, n)

    ! --- 2. 初始化 ---
    x1 = x0
    r0 = b - matmul(a, x1)
    call lowtri(L, r0, y, n)
    call uptri(transpose(L), y, z0, n)
    p0 = z0

    do k = 1, iter_max
        tmp1 = dot_product(r0, z0)
        tmp2 = dot_product(p0, matmul(a, p0))
        afa = tmp1 / tmp2

        x2 = x1 + afa * p0

        dr_s = sqrt(dot_product(r0, r0))
        if (dr_s < eps) exit

        r1 = r0 - afa * matmul(a, p0)
        call lowtri(L, r1, y, n)
        call uptri(transpose(L), y, z1, n)

        tmp3 = dot_product(r1, z1)
        beta = tmp3 / tmp1

        p1 = z1 + beta * p0

        write(*,503) k, (x2(i), i=1, n)

        ! 全部更新
        r0 = r1
        z0 = z1
        p0 = p1
        x1 = x2
    end do
    x = x2

501 format(//,18x,"Cholesky-PCG Iterative",//)
503 format(i5, 4(3x, f12.8))
end subroutine pcg_c