Cg.c
/*
Large Scale Computing
solve Ax=b
Conjugate Gradient Method
General Sparse Matrix A
solve_cg solve Ax=b with CG method
this returns the number of iterations, or negative if failed.
Note: A must be symmetric and positive-defined matrix
*/
#include<stdio.h>
int solve_cg(int dim, /* dimension */
int nnz, /* # of non-zeros in the matrix */
double *nzval, /* array of nonzero values */
int *colidx, /* array of column indices of the nonzeros */
int *rowptr, /* the first column of nonzero in each row */
/* rowptr[j] stores the location of nzval[] and colidx[], */
/* which starts row j. This has dim+1 entry that is nnz */
double *x, /* solition */
double *b, /* right hand side */
double tol); /* tolerance */
/*
Large Scale Computing
solve Ax=b
Conjugate Gradient Method
General Sparse Matrix A
Note: A must be symmetric and positive-defined matrix
*/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
int solve_cg(int dim, /* dimension */
int nnz, /* # of non-zeros in the matrix */
double *nzval, /* array of nonzero values */
int *colidx, /* array of column indices of the nonzeros */
int *rowptr, /* the first column of nonzero in each row */
double *x, /* solition */
double *b, /* right hand side */
double tol){ /* tolerance */
double *r;
double *p;
double *Ap;
double rr,pAp,rr1,alpha,beta;
int i,j,count;
/* Allocate Working Spaces */
r = (double *)calloc(dim, sizeof(double));
p = (double *)calloc(dim, sizeof(double));
Ap = (double *)calloc(dim, sizeof(double));
if ((r == NULL) || (p == NULL) || (Ap == NULL)){
printf("memory allocation failed\n");
return -1;
}
/* compute r0 */
for (i = 0 ; i < dim ; i++){
r[i] = b[i];
for (j = rowptr[i] ; j < rowptr[i+1] ; j++){
r[i] -= nzval[j] * x[colidx[j]];
}
}
rr = 0.0;
for (i = 0 ; i < dim ; i++){
p[i] = r[i]; /* p = r */
rr += r[i] * r[i]; /* rr = r.r */
}
/* cg iteration */
count = 0;
while(rr > tol){
// Ap = A*p
for (i = 0 ; i < dim ; i++){
Ap[i] = 0.0;
for (j = rowptr[i] ; j < rowptr[i+1] ; j++){
Ap[i] += nzval[j] * p[colidx[j]];
}
}
// alpha = r.r / p.Ap
pAp = 0.0;
for (i = 0 ; i < dim ; i++){
pAp += p[i] * Ap[i];
}
alpha = rr / pAp;
//Beta
rr1 = 0.0;
for (i = 0 ; i < dim ; i++){
x[i] += alpha * p[i];
r[i] -= alpha * Ap[i];
rr1 += r[i] * r[i];
}
beta = rr1 / rr;
for (i = 0 ; i < dim ; i++){
p[i] = r[i] + beta * p[i];
}
rr = rr1;
count++;
}
/* Deallocate Working Spaces */
free(r);
free(p);
free(Ap);
return count;
}