Version 1 (modified by cmaggio, 7 years ago) (diff)

migrated content from ccs wiki

# Cg.c

```/*
Large Scale Computing

solve Ax=b
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
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;
}
```