2次元FDM用CG法NO.1(FORTRAN)プログラム

    2003/08/24 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------

1. 概要

 2次元差分法で離散化した疎行列を係数とする連立一次方程式Ax=bに対してCG法で
反復解xを求める。
α=(p,r)/(p,Ap),β=-(r,Ap)/(p,Ap)を使用した、直交性に比較的強いバージョン

2. プログラム

C=================================================================C
      SUBROUTINE CG2D1(A,B,NX,NY,X,R,P,Q,EPS,ITER,ERR,IERR)
C=================================================================C
C  Solve Ax=b by CG NO.1 with 2 dimensional FDM                   C
C    Alpha=(P,R)/(P,AP),   Beta=-(R,AP)/(P,AP)                    C
C-----------------------------------------------------------------C
C    A(0:NX,0:NY,3)  R*8, In, A Coefficient Matrix                C
C                             with symmetric                      C
C    B(0:NX,0:NY) R*8, In,  A Right-hand Vector(b)                C
C    NX           I*4, In,  Grid Numbers on X-axis                C
C    NY           I*4, In,  Grid Numbers on Y-axis                C
C    X(0:NX,0:NY) R*8, I/O, Initial and Solution vector           C
C    R(0:NX,0:NY) R*8, Out, Residual vector                       C
C    P(0:NX,0:NY) R*8, Wk,  Work Vector (P)                       C
C    Q(0:NX,0:NY) R*8, Wk,  Work Vector (Q)                       C
C    EPS          R*8, In,  if ||r||/||b|| <= EPS --> return      C
C    ITER         I*4, I/O, Number of Iteration                   C
C    ERR          R*8, Out, ERR=||r||/||b||                       C
C    IERR         I*4, Out, IERR=0,  Normal Return                C
C                               =1,  No Convergent                C
C-----------------------------------------------------------------C
C    Written by Yasunori Ushiro,   2003/08/17                     C
C        ( Hitachi Ltd. and Waseda University )                   C
C=================================================================C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(0:NX,0:NY,3), B(0:NX,0:NY)
      DIMENSION X(0:NX,0:NY), R(0:NX,0:NY)
      DIMENSION P(0:NX,0:NY), Q(0:NX,0:NY)
C  P=R=B-A*X
      IERR = 0
      BN   = 0.0
      do j=1,NY-1
        do i=1,NX-1
          BN     = BN + B(i,j)**2
          R(i,j) = B(i,j) - A(i,j,1)*X(i,j-1) - A(i,j,2)*X(i-1,j)
     1                    - A(i,j,3)*X(i,j) - A(i+1,j,2)*X(i+1,j)
     2                    - A(i,j+1,1)*X(i,j+1)
          P(i,j) = R(i,j)
        end do
      end do
C  Main Loop
      do k=1,ITER
C   Q=A*P, C=(P,Q), Alpha=(P,R)/C
        C     = 0.0
        Alpha = 0.0
        do j=1,NY-1
         do i=1,NX-1
            Q(i,j) = A(i,j,1)*P(i,j-1) + A(i,j,2)*P(i-1,j)
     1             + A(i,j,3)*P(i,j) + A(i+1,j,2)*P(i+1,j)
     2             + A(i,j+1,1)*P(i,j+1)
            C      = C + P(i,j)*Q(i,j)
            Alpha  = Alpha + P(i,j)*R(i,j)
          end do
        end do
        Alpha = Alpha/C
C   X=X+Alpha*P,  R=R-Alpha*Q
        RN   = 0.0
        Beta = 0.0
        do j=1,NY-1
          do i=1,NX-1
            X(i,j) = X(i,j) + Alpha*P(i,j)
            R(i,j) = R(i,j) - Alpha*Q(i,j)
            RN     = RN + R(i,j)**2
            Beta   = Beta + R(i,j)*Q(i,j)
          end do
        end do
C   if(ERR <= EPS) Convergent,  Beta=-(R,Q)/(P,Q)
        ERR = SQRT(RN/BN)
        if(ERR.le.EPS) go to 100
        Beta = -Beta/C
C   P=R+Beta*P
        do j=1,NY-1
          do i=1,NX-1
            P(i,j) = R(i,j) + Beta*P(i,j)
          end do
        end do
      end do
      IERR = 1
C
  100 continue
      ITER = k
C
      RETURN
      END