//=================================================================C
//  Test Program of CG for 3-Dimensional FDM  with symmetric       C
//    for symmetric solution on Y-Z-axis                           C
//-----------------------------------------------------------------C
//    Written by Yasunori Ushiro,   2007/06/04                     C
//        ( Tokyo Polytechnic University )                         C
//    後 保範（東京工芸大学）                                      C
//=================================================================C
#include <stdio.h>
#include <math.h>
#include <time.h>
// Global Define
#define NDX  51
#define NDY  51
#define NDZ  51
#define A(i,j,k,l)  A[l-1][k][j][i]
#define B(i,j,k)    B[k][j][i]
#define X(i,j,k)    X[k][j][i]
double A[4][NDZ][NDY][NDX], B[NDZ][NDY][NDX], X[NDZ][NDY][NDX] ;
double R[NDZ][NDY][NDX], P[NDZ][NDY][NDX], Q[NDZ][NDY][NDX] ;
FILE   *FT1 ;
// Function Definition
void SET3DS(int NX, int NY, int NZ) ;
void OUT3D(int NX, int NY, int NZ, int ID) ;
int CG3DS(int NX, int NY, int NZ, double EPS, int *ITER, double *ERR) ;
int CG3D2(int NX, int NY, int NZ, double EPS, int *ITER, double *ERR) ;
// Main Program
void main ()
{ int     k, NX, NY, NZ, ITER, IERR, T1, T2 ;
  double  EPS, ERR, CPU, F1, F2, FLOP, FLOPS ;
//  Open File
  FT1 = fopen("C-CG3DS.txt","w") ;
  EPS = 1.0e-8 ;
//  Size Input
  printf("Type In  NX,NY,NZ (<=300) \n") ;
  scanf("%d %d %d",&NX,&NY,&NZ) ;
  if(NX > NDX-1) NX = NDX-1 ;
  if(NY > NDY-1) NY = NDY-1 ;
  if(NZ > NDZ-1) NZ = NDZ-1 ;
  fprintf(FT1," Type  Loop   Error   Time(s)   MFLOPS \n") ;
  for (k=1; k<=2; k++) {
//  Set A,B,X
    printf("----- 3D FDM with CG -------------------------- \n") ;
    fprintf(FT1,"----- 3D FDM with CG -------------------------- \n") ;
    if(k == 1) {
      fprintf(FT1,"Symmetric solution CG for Y-axis \n") ; }
     else {
      fprintf(FT1,"Unsymmetric solution CG for Y-axis \n") ; }
    fprintf(FT1," CG Method for 3-Dimensional FDM,  NX=%d, NY=%d, NZ=%d \n",
          NX,NY,NZ) ;
    printf(" CG Method for 3-Dimensional FDM,  NX=%d, NY=%d, NZ=%d \n",
          NX,NY,NZ) ;
    SET3DS(NX, NY, NZ) ;
//  Solve Ax=b by CG
    ITER = NX*NY*NZ ;
    T1 = clock() ;
    F1 = 16 ;
    F2 = 23 ;
    if(k == 1) {
//   Symmetric solution
       IERR = CG3DS(NX, NY, NZ, EPS, &ITER, &ERR) ;  }
      else {
//   Unsymmetric solution
       IERR = CG3D2(NX, NY, NZ, EPS, &ITER, &ERR) ;  }
    T2 = clock() ;
//  Write
    CPU  = (double)(T2 - T1)/CLOCKS_PER_SEC ;
    FLOP = (NX-1)*(NY-1)*(NZ-1)*(F1+F2*ITER)*1.0e-6 ;
    FLOPS = 0.0 ;
    if(CPU > 0.0) FLOPS = FLOP/CPU ;
    fprintf(FT1,"%3d %5d %10.2e %8.1f %8.1f \n",
                  k,ITER,ERR,CPU,FLOPS) ;
    printf("%3d %5d %10.2e %8.1f %8.1f \n",k,ITER,ERR,CPU,FLOPS) ;
//  Output center-line
    OUT3D(NX, NY, NZ, 1) ;
   }
 }
//=================================================================C
void SET3DS(int NX, int NY, int NZ)
//=================================================================C
//  Set A and B by 3 Dimensional FDM  with symmetric               C
//    -div(K.grad(X)) = f,   K=1.0, f=10,0 and 1x1 square          C
//-----------------------------------------------------------------C
//    NX           I*4, In,  Grid Numbers on X-axis                C
//    NY           I*4, In,  Grid Numbers on Y-axis                C
//    NZ           I*4, In,  Grid Numbers on Y-axis                C
//-----------------------------------------------------------------C
//    Written by Yasunori Ushiro,   2007/06/04                     C
//        ( Tokyo Polytechnic University )                         C
//=================================================================C
{ int    i, j, k ;
  double DX, DY, DZ, HX, HY, HZ, F ;
//  Clear A,B,X
  for (k=0; k<=NZ; k++) {
   for (j=0; j<=NY; j++) {
    for (i=0; i<=NX; i++) {
      A(i,j,k,1) = 0.0 ;
      A(i,j,k,2) = 0.0 ;
      A(i,j,k,3) = 0.0 ;
      A(i,j,k,4) = 0.0 ;
      B(i,j,k)   = 0.0 ;
      X(i,j,k)   = 0.0 ; }
   } }
//  Initial Data
  DX = NX ;
  DY = NY ;
  DZ = NZ ;
  HX = 1.0/DX ;
  HY = 1.0/DY ;
  HZ = 1.0/DZ ;
  F  = 100.0 ;
//  Set A,B
  for (k=1; k<=NZ-1; k++) {
   for (j=1; j<=NY-1; j++) {
    for (i=1; i<=NX-1; i++) {
      A(i,j,k,1) = -HX*HY*DZ ;
      if(k == 1) A(i,j,k,1) = 0.0 ;
      A(i,j,k,2) = -HZ*HX*DY ;
      if(j == 1) A(i,j,k,2) = 0.0 ;
      A(i,j,k,3) = -HY*HZ*DX ;
      if(i == 1) A(i,j,k,3) = 0.0 ;
      A(i,j,k,4) = 2.0*(HY*HZ*DX+HZ*HX*DY+HX*HY*DZ) ;
      B(i,j,k)   = HX*HY*HZ*F ;  }
   } }
 }
//=================================================================C
void OUT3D(int NX, int NY, int NZ, int ID)
//-----------------------------------------------------------------C
//    Written by Yasunori Ushiro,   2007/06/04                     C
//        ( Tokyo Polytechnic University )                         C
//=================================================================C
{ int    i, j, k, NO, *a, *b, *c, *d ;
  double DX, XV ;
//  Output Y
  fprintf(FT1,"Output on YZ-center line \n","  X-axis  Solution(X) \n") ;
  DX = NX ;
  for (i=0; i<=NX; i++) {
    XV = i/DX ;
    fprintf(FT1,"%10.6f %10.6f \n",XV,X(i,NY/2,NZ/2)) ; }
//  Output Un-symmetric Solution
  if(ID == 1) {
    fprintf(FT1," ** Un-symmetrix Solution ** \n") ; 
    fprintf(FT1,"  i  j  k     X(i,j,k)        X(i,NY-j,k)   "
                "   X(i,j,NZ-k)     X(i,NY-j,NZ-k)\n") ;
    NO = 0 ;
    for (k=1; k<=NZ/2-1; k++) {
     for (j=1; j<=NY/2-1; j++) {
      for (i=1; i<=NX/2-1; i++) {
       if(X(i,j,k) != X(i,NY-j,k) || X(i,j,k) != X(i,j,NZ-k)
          || X(i,j,k) != X(i,NY-j,NZ-k) ) {
         a = &X(i,j,k) ;      b = &X(i,NY-j,k) ;
         c = &X(i,j,NZ-k) ;   d = &X(i,NY-j,NZ-k) ;
//  Little Endium  (pentium)
         fprintf(FT1,"%3d%3d%3d  %08X%08X %08X%08X %08X%08X %08X%08X\n",
                 i,j,k,*(a+1),*a,*(b+1),*b,*(c+1),*c,*(d+1),*d ) ;
//  Big Endium  (SR8000)
//       fprintf(FT1,"%3d%3d%3d  %08X%08X %08X%08X %08X%08X %08X%08X\n",i,j,k,
//               i,j,k,*a,*(a+1),*b,*(b+1),*c,*(c+1),*d,*(d+1) ) ;
         NO = NO + 1 ;  }
	}  }  }
     printf("Un-symmetric solution numbers= %d \n",NO) ;
     fprintf(FT1,"Un-symmetric solution numbers= %d \n",NO) ;
   }
  }