Ex21.c

/*
  Large Scale Computing
  1D Heat/Mass Transfer
  ex21.c : use cg.c
 */
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include"cg.h"
 
int main(int argc, char **argv){
  double a,b,d;
  double *phi;
  double *rhs;
  double dx;
  double tol = 1.0e-12;
  int i,num,itc;
  double *nzval; /* non-zero elements of A */
  int *colind; /* column indices of nonzeros */
  int *rowptr; /* row info */
  int nnz;
  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));
  rhs = (double *)calloc(num, sizeof(double));
  nzval = (double *)calloc(3 * num, sizeof(double));
  colind = (int *)calloc(3 * num, sizeof(int));
  rowptr = (int *)calloc(num + 1, sizeof(int));
  if ((phi == NULL) || (rhs == NULL) ||
      (nzval == NULL) || (colind == NULL) || (rowptr == NULL)){
    printf("Memory Allocation Failed\n");
    exit(-1);    
  }  
 
  /** Setup **/
  dx = 1.0 / (double)(num - 1);
 
  /* Boundary Condition*/
  phi[0] = a;
  phi[num-1] = b;
 
  /* Setup Matrix A and RHS */
  nnz = 0;
  for (i = 1 ; i < num - 1 ; i++){
 
    /* set general RHS */
    rhs[i] = -dx * dx / d * (dx * (double)i); 
 
    /* set first nonzero column of row index(i,j) */
    rowptr[i - 1] = nnz;
 
    /* west */
    if (i == 1){
      rhs[i] -= phi[0]; 
    }else{
      nzval[nnz] = 1.0;
      colind[nnz] = i-2;
      nnz++;
    } 
 
    /* diagonal Element */
    nzval[nnz] = -2.0;
    colind[nnz] = i - 1;
    nnz++;
 
    /* east */
    if (i == num - 2){
      rhs[i] -= phi[num-1]; 
    }else{
      nzval[nnz] = 1.0;
      colind[nnz] = i;
      nnz++;
    } 
  }
 
  /* last element of rowptr si nnz */
  rowptr[num-2] = nnz;
 
  /* solve with CG */
  itc = solve_cg(num-2, nnz, nzval, colind, rowptr, &phi[1], &rhs[1], tol);
  if (itc < 0){
    printf("calculation failed\n");
    exit(-1);
  }
 
  printf("Number of Iteration=%d\n",itc);
 
  /* 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(rhs);
  free(nzval);
  free(colind);
  free(rowptr);
 
  printf("Done\n");
}

cg.c