对称逐步超松弛预处理共轭梯度迭代
2015-07-23 19:39:06   来源：Fcode研讨团队   评论：0 点击：

在数值计算中，经常会遇到求解大型、稀疏、复系数的线性方程组的问题，通常采用直接法或者迭代法求解。常用的开源直接求解器有SuperLU，MUMPS、UMFPACK，这些求解器的具体算法和存储格式可从其手册获得，这里不再赘述。直接求解器需要大量的内存，当内存不够或者矩阵过于庞大时，迭代法是较好的选择。迭代求解器比较多，常用预处理共轭梯度迭代算法，这里介绍的是改进的对称逐步超松弛预处理共轭梯度迭代算法，算法参考林绍忠《用预处理共轭梯度法求解有限元方程组及程序设计》，其算法如下：
设线性方程组 Ax=b 的系数矩阵 A 为 n 阶对称正定矩阵，用对称逐步超松弛迭代(SSOR法)的分裂矩阵作为预处理矩阵M，
M=\left ( 2-\omega \right )^{-1}\left ( D/\omega +L \right )\left ( D/\omega \right )^{-1}\left ( D/\omega +L \right )^{T}
式中：D 为 A 的对角阵；L 为 A 的严格下三角矩阵；0<\omega<2 为松弛因子。
SSOR-PCG法的迭代格式为：

林绍忠对其进行了改写，节省了一定的计算量，其具体算法为：令W=D/\omega+L，V=(2-\omega)D/\omega，y=W^{-1}g，Z=W^{T}d，则：


于是，上述迭代公式可以改写为：

代码及部分演示代码（请自行补充输入输出部分）如下：

Subroutine ssor_pcg_csr(nd, nnz, a, b, x, irw, jcol, eps, itmax, w, it, reg_err)
  ! ***********************************************************************************
  ! **   本程序参考林绍忠《用预处理共轭梯度法求解有限元方程组及程序设计》改进的迭代格式
  ! **    采用SSOR预处理，预处理矩阵为:M=(2-W)^{-1}*(D/W+L)*(D/W)^{-1}*(D/W+L)^{T}
  ! **    在原来的程序基础上去掉ID映射的部分，而是一个小技巧来解决这个问题
  ! **                   aliouying
  ! **                   2013-7-31
  ! ***********************************************************************************
  Implicit None
  ! ***********************************************************************************
  ! **   输入和输出
  ! ***********************************************************************************
  Integer, Intent (In) :: nd, nnz
  Integer, Intent (In) :: irw(nnz), jcol(nd+1) ! CSR压缩存储索引矩阵pointB,pointE
  Complex (Kind=8), Intent (In) :: a(nnz), b(nd) ! 矩阵变量A*x=B
  Complex (Kind=8), Intent (Inout) :: x(nd) ! 输入初始估计值和输出最终结果
  Real (Kind=8), Intent (In) :: eps ! 收敛误差
  Integer, Intent (In) :: itmax ! 最大迭代次数
  Real (Kind=8), Intent (In) :: w ! 松弛因子
  Integer, Intent (Out) :: it ! 收敛总计迭代次数
  Real (Kind=8), Intent (Out) :: reg_err ! 迭代后的残差
  ! ***********************************************************************************
  ! **   程序中间变量
  ! ***********************************************************************************
  Complex (Kind=8) :: y(nd), z(nd), d(nd), v(nd), temp(nd)
  Integer :: i
  Complex (Kind=8) :: delta, tao, beta, delta1, delta0
  Complex (Kind=8), External :: m_dot_product
  ! ***********************************************************************************
  ! **   设置初始值
  ! ***********************************************************************************
  Forall (i=1:nd)
    v(i) = (2.0-w)*a(jcol(i))/w ! 这里V当成了一维向量，实际上它是一个对角矩阵，及A的对角阵计算来
    ! 若想节约内存，可以把用到V的地方用A(jcol())来代替，这里我主要是为
    ! 了程序跟理论公式的一致性
  End Forall
  ! ***********************************************************************************
  ! **   开始计算
  ! ***********************************************************************************
  Call matmulvec_csr(nd, nnz, a, x, y, irw, jcol)
  y = y - b
  Call forward_solve(nd, nnz, a, y, y, irw, jcol, w)
  Forall (i=1:nd) ! 这里可以直接用Z=V*Y代替，主要目的是为了方便替换V,而且forall的效率不错^_^
    z(i) = -v(i)*y(i)
  End Forall
  Call back_solve(nd, nnz, a, z, d, irw, jcol, w)
  it = 0
  delta = m_dot_product(y, v*y, nd)
  delta0 = delta
  beta = 1
  Do i = 1, itmax
    ! if( abs(delta) < EPS) then
    ! return
    ! endif
    tao = delta/m_dot_product(d, 2.0*z-v*d, nd)
    x = x + tao*d
    ! 多种收敛方式
    ! if( sqrt(sum(abs(tao*D/X)**2)) < EPS ) then
    ! return
    ! endif
    Call forward_solve(nd, nnz, a, z-v*d, temp, irw, jcol, w)
    y = y + tao*(temp+d)
    delta1 = m_dot_product(y, v*y, nd)
    beta = delta1/delta
    z = -v*y + beta*z
    Call back_solve(nd, nnz, a, z, d, irw, jcol, w)
    delta = delta1
    it = it + 1
    reg_err = abs(delta/delta0)
    ! write(*,*) 'PCG iteration:',i,'of res:',reg_err
    If (reg_err<eps) Then
      Write (*, *) 'PCG iteration:', i, 'of res:', reg_err
      Return
    End If
  End Do
  Return
End Subroutine ssor_pcg_csr
! ================this is split line================================================!
Subroutine matmulvec_csr(nd, nnz, mat, vec, vec_out, irw, jcol)
  Implicit None
  ! ***********************************************************************************
  ! **   输入和输出
  ! ***********************************************************************************
  Integer, Intent (In) :: nd, nnz
  Integer, Intent (In) :: irw(nnz), jcol(nd+1) ! CSR压缩存储索引矩阵pointB,pointE
  Complex (Kind=8), Intent (In) :: mat(nnz), vec(nd) ! 矩阵变量MAT和向量VEC
  Complex (Kind=8) :: vec_out(nd) ! 输出最终结果
  ! ***********************************************************************************
  ! **   程序中间变量
  ! ***********************************************************************************
  Integer :: i, j
  vec_out = (0.0D0, 0.0D0)
  Do i = 1, nd
    vec_out(i) = vec_out(i) + mat(jcol(i))*vec(i)
    Do j = jcol(i) + 1, jcol(i+1) - 1
      vec_out(i) = vec_out(i) + mat(j)*vec(irw(j))
      vec_out(irw(j)) = vec_out(irw(j)) + mat(j)*vec(i)
    End Do
  End Do
  Return
End Subroutine matmulvec_csr
! ================this is split line================================================!
Subroutine forward_solve(nd, nnz, mat, vec, vec_out, irw, jcol, w)
  Implicit None
  ! ***********************************************************************************
  ! **   输入和输出
  ! ***********************************************************************************
  Integer, Intent (In) :: nd, nnz
  Integer, Intent (In) :: irw(nnz), jcol(nd+1) ! CSR压缩存储索引矩阵pointB,pointE
  Complex (Kind=8), Intent (In) :: mat(nnz), vec(nd) ! 矩阵变量MAT和向量VEC
  Real (Kind=8), Intent (In) :: w
  Complex (Kind=8), Intent (Out) :: vec_out(nd) ! 输出最终结果
  ! ***********************************************************************************
  ! **   程序中间变量
  ! ***********************************************************************************
  Integer :: i, j
  vec_out = vec
  Do i = 1, nd
    vec_out(i) = vec_out(i)/mat(jcol(i))*w
    Do j = jcol(i) + 1, jcol(i+1) - 1
      vec_out(irw(j)) = vec_out(irw(j)) - mat(j)*vec_out(i)
    End Do
  End Do
  Return
End Subroutine forward_solve
! ================this is split line================================================!
Subroutine back_solve(nd, nnz, mat, vec, vec_out, irw, jcol, w)
  Implicit None
  ! ***********************************************************************************
  ! **   输入和输出
  ! ***********************************************************************************
  Integer, Intent (In) :: nd, nnz
  Integer, Intent (In) :: irw(nnz), jcol(nd+1) ! CSR压缩存储索引矩阵pointB,pointE
  Complex (Kind=8), Intent (In) :: mat(nnz), vec(nd) ! 矩阵变量MAT和向量VEC
  Real (Kind=8), Intent (In) :: w
  Complex (Kind=8), Intent (Out) :: vec_out(nd) ! 输出最终结果
  ! ***********************************************************************************
  ! **   程序中间变量
  ! ***********************************************************************************
  Integer :: i, j
  Complex (Kind=8) :: sum
  Do i = nd, 1, -1
    sum = (0.0D0, 0.0D0)
    Do j = jcol(i+1) - 1, jcol(i) + 1, -1
      sum = sum + mat(j)*vec_out(irw(j))
    End Do
    vec_out(i) = (vec(i)-sum)/mat(jcol(i))*w
  End Do
  Return
End Subroutine back_solve
! ================this is split line================================================!
Function m_dot_product(a_in, b_in, n)
  ! ***********************************************************************************
  ! **   向量的内积
  ! ***********************************************************************************
  Implicit None
  Integer :: n
  Complex (Kind=8) :: a_in(n), b_in(n)
  Complex (Kind=8) :: m_dot_product
  Integer :: i
  m_dot_product = 0.0
  Do i = 1, n
    m_dot_product = m_dot_product + a_in(i)*b_in(i)
  End Do
  Return
End Function m_dot_product
! ================this is split line================================================!
Subroutine id_change(nd, nnz, irw, jcol, icol, jrow, id_ch)
  ! ***********************************************************************************
  ! **   矩阵A影射,由下三角CSC压缩存储影射到下三角CSR压缩存储
  ! ***********************************************************************************
  Implicit None
  ! ***********************************************************************************
  ! **   输入和输出
  ! ***********************************************************************************
  Integer, Intent (In) :: nd, nnz
  Integer, Intent (In) :: irw(nnz), jcol(nd+1) ! CSR压缩存储索引矩阵pointB,pointE
  Integer, Intent (Out) :: icol(nnz), jrow(nd) ! CSC压缩存储索引矩阵pointB,pointE
  Integer, Intent (Out) :: id_ch(nnz)
  ! ***********************************************************************************
  ! **   程序中间变量
  ! ***********************************************************************************
  Integer :: i, j, k, now_col_nnz, n
  icol(1) = 1
  jrow(1) = 1
  now_col_nnz = 1
  id_ch(1) = 1
  n = 1
  Do i = 2, nd
    now_col_nnz = 0
    Do j = 1, i - 1
      Do k = jcol(j) + 1, jcol(j+1) - 1
        If (irw(k)==i) Then
          now_col_nnz = now_col_nnz + 1
          n = n + 1
          icol(n) = j
          id_ch(n) = k
        End If
      End Do
    End Do
    now_col_nnz = now_col_nnz + 1
    n = n + 1
    icol(n) = i
    id_ch(n) = jcol(i)
    jrow(i) = jrow(i-1) + now_col_nnz
  End Do
  Return
End Subroutine id_change

! 下面是测试程序

Program www_fcode_cn
  Implicit None
  Complex (Kind=8) , allocatable :: x(:) , a(:) , rhs(:)
  integer , allocatable :: jcol(:) , irw(:)
  Real (Kind=8) :: w, eps, reg_err
  Integer :: it , itmax , n , nnz
  itmax = 1000
  eps = 1.0D-20
  w = 1
  ! 获得 n 和 nnz
  Allocate (x(n),a(nnz),rhs(n),irw(nnz),jcol(n+1))
  ! 获得 a rhs
  x = 0.0D0
  ! 调用程序
  Call ssor_pcg_csr(n, nnz, a, rhs, x, irw ,jcol, eps, itmax, w, it, reg_err)
  ! 输出解
  ! write(*,*) X,IT,REG_ERR(IT)
  Stop
End Program www_fcode_cn