= Cg.c =
{{{#!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 */
}}}
{{{#!C
/*
  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;
}
}}}