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