C=================================================================C
      SUBROUTINE CG3DS(A,B,NX,NY,NZ,X,R,P,Q,EPS,ITER,ERR,IERR)
C=================================================================C
C  Solve Ax=b by CG NO.1 with 3 dimensional FDM                   C
C    Alpha=(P,R)/(P,AP),   Beta=-(R,AP)/(P,AP)                    C
C     for symmetric solution on x, y and z  axis                  C
C     with pair computation                                       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,   2007/06/04                     C
C        ( Tokyo Polytechnic University )                         C
C     後 保範（東京工芸大学）                                     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
      do k=1,NZ-1
       do j=1,NY-1
        do i=1,NX-1
          BN       = BN + B(i,j,k)**2
C          R(i,j,k) = B(i,j,k) - A(i,j,k,1)*X(i,j,k-1)
C     1             - A(i,j,k,2)*X(i,j-1,k) - A(i,j,k,3)*X(i-1,j,k)
C     2             - A(i,j,k,4)*X(i,j,k) - A(i+1,j,k,3)*X(i+1,j,k)
C     3             - A(i,j+1,k,2)*X(i,j+1,k) - A(i,j,k+1,1)*X(i,j,k+1)
C   Symmetric computation
          R(i,j,k) = B(i,j,k) - A(i,j,k,4)*X(i,j,k)
     1             - (A(i,j,k,3)*X(i-1,j,k) + A(i+1,j,k,3)*X(i+1,j,k))
	2             - (A(i,j,k,2)*X(i,j-1,k) + A(i,j+1,k,2)*X(i,j+1,k))
     3             - (A(i,j,k,1)*X(i,j,k-1) + A(i,j,k+1,1)*X(i,j,k+1))
          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, C=(P,Q), Alpha=(P,R)/C
        C     = 0.0
        Alpha = 0.0
        do k=1,NZ-1
         do j=1,NY-1
          do i=1,NX-1
C            Q(i,j,k) = A(i,j,k,1)*P(i,j,k-1) + A(i,j,k,2)*P(i,j-1,k)
C     1               + A(i,j,k,3)*P(i-1,j,k) + A(i,j,k,4)*P(i,j,k)
C     2               + A(i+1,j,k,3)*P(i+1,j,k) + A(i,j+1,k,2)*P(i,j+1,k)
C     3               + A(i,j,k+1,1)*P(i,j,k+1)
C   Symmetric computation
            Q(i,j,k) = A(i,j,k,4)*P(i,j,k)
     1               + (A(i,j,k,3)*P(i-1,j,k) + A(i+1,j,k,3)*P(i+1,j,k))
     2               + (A(i,j,k,2)*P(i,j-1,k) + A(i,j+1,k,2)*P(i,j+1,k))
     3               + (A(i,j,k,1)*P(i,j,k-1) + A(i,j,k+1,1)*P(i,j,k+1))
            C        = C + P(i,j,k)*Q(i,j,k)
            Alpha    = Alpha + P(i,j,k)*R(i,j,k)
          end do
         end do
        end do
        Alpha = Alpha/C
C   X=X+Alpha*P,  R=R-Alpha*Q
        RN   = 0.0
        Beta = 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)
            RN       = RN + R(i,j,k)**2
            Beta     = Beta + R(i,j,k)*Q(i,j,k)
          end do
         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 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