首页 > 算法 > 正文

阅读排行

对称逐步超松弛预处理共轭梯度迭代
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+LV=(2-\omega)D/\omegay=W^{-1}gZ=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

相关热词搜索:共轭梯度 超松弛 迭代

上一篇:盛金公式解一元三次方程
下一篇:简单数据分割与“打点法”排序

分享到: 收藏