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

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

1. 概要

 3次元差分法で離散化した疎行列を係数とする連立一次方程式Ax=bに対してCG法で
反復解xを求める。
α=(r,r)/(p,Ap),β=new(r,r)/old(r,r)を使用した、演算量の少ないバージョン

2. プログラム

C=================================================================C
      SUBROUTINE CG3D2(A,B,NX,NY,NZ,X,R,P,Q,EPS,ITER,ERR,IERR)
C=================================================================C
C  Solve Ax=b by CG NO.2 with 3 dimensional FDM                   C
C    Alpha=(R,R)/(P,AP),   Beta=new(R,R)/old(R,R)                 C
C-----------------------------------------------------------------C
C    A(0:NX,0:NY,0:NZ,4)  R*8, In, A Coefficient Matrix           C
C                                    with symmetric               C
C    B(0:NX,0:NY,0:NZ) 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    NZ                I*4, In,  Grid Numbers on Z-axis           C
C    X(0:NX,0:NY,0:NZ) R*8, I/O, Initial and Solution vector      C
C    R(0:NX,0:NY,0:NZ) R*8, Out, Residual vector                  C
C    P(0:NX,0:NY,0;NZ) R*8, Wk,  Work Vector (P)                  C
C    Q(0:NX,0:NY,0:NZ) 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,0:NZ,4), B(0:NX,0:NY,0:NZ)
      DIMENSION X(0:NX,0:NY,0:NZ), R(0:NX,0:NY,0:NZ)
      DIMENSION P(0:NX,0:NY,0:NZ), Q(0:NX,0:NY,0:NZ)
C  P=R=B-A*X
      IERR = 0
      BN   = 0.0
      C0   = 0.0 
      do k=1,NZ-1
       do j=1,NY-1
        do i=1,NX-1
          BN       = BN + B(i,j,k)**2
          R(i,j,k) = B(i,j,k) - A(i,j,k,1)*X(i,j,k-1)
     1             - A(i,j,k,2)*X(i,j-1,k) - A(i,j,k,3)*X(i-1,j,k)
     2             - A(i,j,k,4)*X(i,j,k) - A(i+1,j,k,3)*X(i+1,j,k)
     3             - A(i,j+1,k,2)*X(i,j+1,k) - A(i,j,k+1,1)*X(i,j,k+1)
          C0       = C0 + R(i,j,k)**2
          P(i,j,k) = R(i,j,k)
        end do
       end do
      end do
C  Main Loop
      do kk=1,ITER
C   Q=A*P, Alpha=C0/(P,Q)
        Alpha = 0.0
        do k=1,NZ-1
         do j=1,NY-1
          do i=1,NX-1
            Q(i,j,k) = A(i,j,k,1)*P(i,j,k-1) + A(i,j,k,2)*P(i,j-1,k)
     1               + A(i,j,k,3)*P(i-1,j,k) + A(i,j,k,4)*P(i,j,k)
     2               + A(i+1,j,k,3)*P(i+1,j,k) + A(i,j+1,k,2)*P(i,j+1,k)
     3               + A(i,j,k+1,1)*P(i,j,k+1)
            Alpha    = Alpha + P(i,j,k)*Q(i,j,k)
          end do
         end do
        end do
        Alpha = C0/Alpha
C   X=X+Alpha*P,  R=R-Alpha*Q
        C1   = 0.0
        do k=1,NZ-1
         do j=1,NY-1
          do i=1,NX-1
            X(i,j,k) = X(i,j,k) + Alpha*P(i,j,k)
            R(i,j,k) = R(i,j,k) - Alpha*Q(i,j,k)
            C1       = C1 + R(i,j,k)**2 
          end do
         end do
        end do
C  if(ERR <= EPS) Convergent,  Beta=C1/C0
        ERR = SQRT(C1/BN)
        if(ERR.le.EPS) go to 100
        Beta = C1/C0
        C0   = C1 
C   P=R+Beta*P
        do k=1,NZ-1
         do j=1,NY-1
          do i=1,NX-1
            P(i,j,k) = R(i,j,k) + Beta*P(i,j,k)
          end do
         end do
        end do
      end do
      IERR = 1
C
  100 continue
      ITER = kk
C
      RETURN
      END