有限元方法中的数值技术：预处理共轭梯度法 PCG (2)
Cholesky预处理： 预处理矩阵 M=LLT\mathbf{M} = LL^TM=LLT，其中 LLL 为 AAA 的下三角Cholesky分解
每步的PCG主要分量公式如下：
-
预处理向量 zzz 的计算
（两步三角系统求解，逐分量写为）：-
先解 Ly=rL y = rLy=r：
y1=r1L11 y_1 = \frac{r_1}{L_{11}} y1=L11r1
yi=1Lii(ri−∑j=1i−1Lijyj),i=2,…,n y_i = \frac{1}{L_{ii}}\left(r_i - \sum_{j=1}^{i-1} L_{ij} y_j\right),\quad i=2,\ldots,n yi=Lii1(ri−j=1∑i−1Lijyj),i=2,…,n -
再解 LTz=yL^T z = yLTz=y：
zn=ynLnn z_n = \frac{y_n}{L_{nn}} zn=Lnnyn
zi=1Lii(yi−∑j=i+1nLjizj),i=n−1,…,1 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 zi=Lii1(yi−j=i+1∑nLjizj),i=n−1,…,1
-
-
其余 αk\alpha_kαk、x(k+1)x^{(k+1)}x(k+1)、r(k+1)r^{(k+1)}r(k+1)、βk\beta_kβk、p(k+1)p^{(k+1)}p(k+1) 的表达方式完全和Jacobi预处理一致，只是预处理向量zzz的获得方式不同。
subroutine pcg_c(a, b, x, x0, n)! Cholesky预处理共轭梯度法（适用于对称正定A）integer, intent(in) :: nreal*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_maxreal*8 :: eps, tmp1, tmp2, tmp3, afa, beta, dr_siter_max = 1000eps = 1.0d-7write(*,501)! --- 1. Cholesky分解 A = LL^T ---call cholesky(a, L, n)! --- 2. 初始化 ---x1 = x0r0 = b - matmul(a, x1)call lowtri(L, r0, y, n)call uptri(transpose(L), y, z0, n)p0 = z0do k = 1, iter_maxtmp1 = dot_product(r0, z0)tmp2 = dot_product(p0, matmul(a, p0))afa = tmp1 / tmp2x2 = x1 + afa * p0dr_s = sqrt(dot_product(r0, r0))if (dr_s < eps) exitr1 = 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 / tmp1p1 = z1 + beta * p0write(*,503) k, (x2(i), i=1, n)! 全部更新r0 = r1z0 = z1p0 = p1x1 = x2end dox = x2501 format(//,18x,"Cholesky-PCG Iterative",//)
503 format(i5, 4(3x, f12.8))
end subroutine pcg_c