= Ex19.c =
{{{#!C
/*
  Large Scale Computing
  2D Heat/Mass Transfer
  ex19.c : Numerical Analysis
  Solve for
  d2c/dx2 + d2c/dy2 + 1 = 0
  With the boundary conditions of c = 0 
  along lines of x=1,-1, and y=1, and -1
 
  Conjugate Gradient Method
 */
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
 
#define NUM 100  /* number of points for x and y*/
#define TOL 1.0e-12  /* convergence tolerance */
 
int main(void ){
  int i,j,n;
  double x,y,dx;
  double cnc[NUM * NUM]; /* define 1D array of length NUM * NUM */
  double r[NUM * NUM];
  double p[NUM * NUM];
  double Ap[NUM * NUM];
  double rr,pAp,rr1,alpha,beta;
  int k;
  FILE *fp;
 
  /*assuming dx = dy : domain is 2x2 */
  dx = 2.0 / (double)(NUM - 1);
 
  /* Initialize */
  for (n = 0 ; n < NUM * NUM ; n++){
    cnc[n] = 0.0;
    r[n] = 0.0;
  }
 
  /* computing for C with Conjugate Gradient Method */
  /* r0 = b - Ax0 */
  for (n = 0 ; n < NUM * NUM ; n++){
    i = n % NUM; /* i is residual */
    j = n / NUM; /* j is division */
    if (i == 0) continue; /*skip left end*/
    if (i == NUM-1) continue; /*skip right end*/
    if (j == 0) continue; /*skip bottom end*/
    if (j == NUM-1) continue; /*skip top end*/
    r[n] = -dx * dx;
    r[n] -= -4.0 * cnc[n];
    r[n] -= cnc[n + 1];
    r[n] -= cnc[n + NUM];
    r[n] -= cnc[n - 1];
    r[n] -= cnc[n - NUM];
  }
 
  rr = 0.0;
  for (n = 0 ; n < NUM * NUM ; n++){
    p[n] = r[n];    /* p = r */
    rr += r[n] * r[n];  /* rr = r.r */
  }
 
  k = 0;
  while(rr > TOL){    
    // Ap = A*p
    for (n = 0 ; n < NUM * NUM ; n++){
      Ap[n] = 0.0;
      i = n % NUM; /* i is residual */
      j = n / NUM; /* j is division */
      if (i == 0) continue; /*skip left end*/
      if (i == NUM-1) continue; /*skip right end*/
      if (j == 0) continue; /*skip bottom end*/
      if (j == NUM-1) continue; /*skip top end*/
 
      /* computing A*p */
      Ap[n] = -4.0 * p[n];
      Ap[n] += p[n + 1];
      Ap[n] += p[n + NUM];
      Ap[n] += p[n - 1];
      Ap[n] += p[n - NUM];
    }    
 
    // alpha = r.r / p.Ap
    pAp = 0.0;
    for (n = 0 ; n < NUM * NUM ; n++){
      pAp += p[n] * Ap[n];
    }
    alpha = rr / pAp;
 
    //Beta
    rr1 = 0.0;
    for (n = 0 ; n < NUM * NUM ; n++){
      cnc[n] += alpha * p[n];
      r[n] -= alpha * Ap[n];
      rr1 += r[n] * r[n];
    }  
 
    beta = rr1 / rr;
 
    for (n = 0 ; n < NUM * NUM ; n++){
      p[n] = r[n] + beta * p[n];
    }
 
    rr = rr1;
    k++;
  }
 
  printf("# of Iteration=%d\n",k);
 
  /* Output Result */
  fp = fopen("res.dat","w");
  if (fp == NULL){
    printf("File could not create\n");
    exit(-1);
  }
 
  for (i = 0 ;i < NUM ; i++){
    for (j = 0 ; j < NUM ; j++){
      x = -1.0 + 2.0 * (double)i / (double)(NUM - 1);
      y = -1.0 + 2.0 * (double)j / (double)(NUM - 1);
      fprintf(fp,"%e %e %e\n",x,y,cnc[i + NUM * j]);
    }
  }
  fclose(fp);
 
  printf("Done\n");
}
}}}