C=================================================================C
      SUBROUTINE CG2D2(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=(R,R)/(P,AP),   Beta=new(R,R)/old(R,R)                 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,   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,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
      C0   = 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)
          C0     = C0 + R(i,j)**2
          P(i,j) = R(i,j)
        end do
      end do
C  Main Loop
      do k=1,ITER
C   Q=A*P, Alpha=C0/(P,Q)
        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)
            Alpha  = Alpha + P(i,j)*Q(i,j)
          end do
        end do
        Alpha = C0/Alpha
C   X=X+Alpha*P,  R=R-Alpha*Q
        C1  = 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)
            C1     = C1 + R(i,j)**2
          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 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