= Ex13.c = 
a parallel version of ex3.c
{{{#!C
/*
  Large Scale Computing
  Heat/Mass Transfer
  ex13.c 
 */
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<omp.h>
 
int main(int argc, char **argv){
  double a,b,d;
  double *phi,*phi0;
  double *rhs;
  double dx;
  double err,merr,share_err;
  double tol = 1.0e-8;
  int i,num,itc;
  int nthreads;
  FILE *fp;
 
  if (argc < 5){
    printf("Usage:%s [NUM] [A] [B] [D]\n",argv[0]);
    exit(-1);
  }
 
  /* set parameters */
  num = atoi(argv[1]);
  a = atof(argv[2]);
  b = atof(argv[3]);
  d = atof(argv[4]);
 
  printf("num=%d A=%e B=%e D=%e\n",num,a,b,d);
 
  /* Memory Allocation and Check*/
  phi = (double *)calloc(num, sizeof(double));
  if (phi == (double *)NULL){
    printf("Memory Allocation Failed\n");
    exit(-1);    
  }
 
  phi0 = (double *)calloc(num, sizeof(double));
  if (phi0 == (double *)NULL){
    printf("Memory Allocation Failed\n");
    exit(-1);    
  }
 
  rhs = (double *)calloc(num, sizeof(double));
  if (rhs == (double *)NULL){
    printf("Memory Allocation Failed\n");
    exit(-1);    
  }  
 
  /** Setup **/
  dx = 1.0 / (double)(num - 1);
  for (i = 0 ; i < num ; i++){
    rhs[i] = -dx * dx / d * (dx * (double)i); 
    phi[i] = phi0[i] = 0.0;
  }
 
  /* Boundary Condition*/
  phi[0] = phi0[0] = a;
  phi[num-1] = phi0[num-1] = b;
 
  /* Solve with Jacobi Method*/
  itc = 0;
 
  while(1){
    share_err = 0.0;
 
#pragma omp parallel private(merr,err)
    {
      nthreads = omp_get_num_threads();  
 
      /* update phi */
#pragma omp for
      for (i = 1 ; i < (num - 1) ; i++){
	phi[i] = -0.5 * (rhs[i] - phi0[i + 1] - phi0[i - 1]);
      }
 
      /* Check Convergence */
      merr = 0.0;
#pragma omp for 
      for (i = 0 ; i < num ; i++){
	err = fabs(phi[i] - phi0[i]);
	if (merr < err) merr = err;
	phi0[i] = phi[i];
      }
 
#pragma omp critical
      {
	/* find max in whole threads */
	if (share_err < merr) share_err = merr;
      }
    }
    itc++;
 
    if (share_err < tol) break;
    if ((itc % 10000) == 0) printf("itc=%d err=%e\n",itc,share_err);
  }
  printf("Number of Iteration=%d\n",itc);
  printf("Number of Threads=%d\n",nthreads);
 
 
  /* Output Result */
  fp = fopen("res.dat","w");
  if (fp == NULL){
    printf("File could not create\n");
    exit(-1);
  }
 
  for (i = 0 ; i < num ; i++){
    fprintf(fp,"%e %e\n",dx * (double)i,phi[i]);
  }
  fclose(fp);
 
  /* END : Deallocate Memory */
  free(phi);
  free(phi0);
  free(rhs);
  printf("Done\n");
}
}}}