#include <stdio.h>
#include <math.h>
// Global Define
#define NDX  301
#define NDY  301
#define A(i,j,l)  A[l-1][j][i]
#define B(i,j)    B[j][i]
#define X(i,j)    X[j][i]
#define R(i,j)    R[j][i]
#define P(i,j)    P[j][i]
#define Q(i,j)    Q[j][i]
extern double A[3][NDY][NDX], B[NDY][NDX], X[NDY][NDX] ;
extern double R[NDY][NDX], P[NDY][NDX], Q[NDY][NDX] ;
extern FILE *FT1 ;
//=================================================================C
int CG2DS(int NX, int NY, double EPS, int *ITER, double *ERR)
//=================================================================C
//  Solve Ax=b by CG NO.2 with 2 dimensional FDM                   C
//    Alpha=(R,R)/(P,AP),   Beta=new(R,R)/old(R,R)                 C
//-----------------------------------------------------------------C
//    NX           I*4, In,  Grid Numbers on X-axis                C
//    NY           I*4, In,  Grid Numbers on Y-axis                C
//    EPS          R*8, In,  if ||r||/||b|| <= EPS --> return      C
//    ITER         I*4, I/O, Number of Iteration                   C
//    ERR          R*8, Out, ERR=||r||/||b||                       C
//    return       I*4, Out, IERR=0,  Normal Return                C
//                               =1,  No Convergent                C
//-----------------------------------------------------------------C
//    Written by Yasunori Ushiro,   2007/06/04                     C
//        ( Tokyo Polytechnic University )                         C
//    後 保範（東京工芸大学）                                      C
//=================================================================C
{ int     i, j, k, IERR ;
  double  BN, C0, C1, Alpha, Beta ;
//  P=R=B-A*X
  IERR = 0 ;
  BN   = 0.0 ;
  C0   = 0.0 ;
  for (j=1; j<=NY-1; j++) {
    for (i=1; i<=NX-1; i++) {
       BN     = BN + B(i,j)*B(i,j) ;
       R(i,j) = B(i,j) - A(i,j,3)*X(i,j)
              - ( A(i,j,1)*X(i,j-1) + A(i,j+1,1)*X(i,j+1) )
              - ( A(i,j,2)*X(i-1,j) + A(i+1,j,2)*X(i+1,j) ) ;
       C0     = C0 + R(i,j)*R(i,j) ;
       P(i,j) = R(i,j) ;  }
     }
//  Main Loop
  for (k=1; k<=*ITER; k++) {
//   Q=A*P, C=(P,Q), Alpha=C0/(P,Q)
    Alpha = 0.0 ;
    for (j=1; j<=NY-1; j++) {
      for (i=1; i<=NX-1; i++) {
        Q(i,j) = A(i,j,3)*P(i,j)
               + ( A(i,j,1)*P(i,j-1) + A(i,j+1,1)*P(i,j+1) )
               + ( A(i,j,2)*P(i-1,j) + A(i+1,j,2)*P(i+1,j) ) ;
        Alpha  = Alpha + P(i,j)*Q(i,j) ;  }
      }
    Alpha = C0/Alpha ;
//   X=X+Alpha*P,  R=R-Alpha*Q
    C1   = 0.0 ;
    for (j=1; j<=NY-1; j++) {
      for (i=1; i<=NX-1; i++) {
        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)*R(i,j) ;  }
      }
//   if(ERR <= EPS) Convergent,  Beta=C1/C0
    *ERR = sqrt(C1/BN) ;
    if(*ERR <= EPS) goto M100 ;
    Beta = C1/C0 ;
    C0   = C1 ;
//   P=R+Beta*P
    for (j=1; j<=NY-1; j++) {
      for (i=1; i<=NX-1; i++) {
        P(i,j) = R(i,j) + Beta*P(i,j) ;  }
      }
   }
   IERR = 1 ;
 M100: *ITER = k ;
 return (IERR) ; }