Bcg.c
/*
Large Scale Computing
solve Ax=b
Bi-Conjugate Gradient Method
General Sparse Matrix A
solve_cg solve Ax=b with BCG method
this returns the number of iterations, or negative if failed.
*/
#include<stdio.h>
int solve_bcg(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
Bi-Conjugate Gradient Method
General Sparse Matrix A
*/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
int solve_bcg(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 *rt;
double *pt;
double *ATpt;
double *nzvalT; /* AT array of nonzero values */
int *colidxT; /* AT array of column indices of the nonzeros */
int *rowptrT; /* AT the first column of nonzero in each row */
double rr,ptAp,rr1,alpha,beta;
int i,j,count;
int col, relpos;
int *marker;
/* Allocate Working Spaces */
r = (double *)calloc(dim, sizeof(double));
rt = (double *)calloc(dim, sizeof(double));
p = (double *)calloc(dim, sizeof(double));
pt = (double *)calloc(dim, sizeof(double));
Ap = (double *)calloc(dim, sizeof(double));
ATpt = (double *)calloc(dim, sizeof(double));
nzvalT = (double *)calloc(nnz, sizeof(double));
colidxT = (int *)calloc(nnz, sizeof(int));
rowptrT = (int *)calloc(dim+1, sizeof(int));
if ((r == NULL) || (rt == NULL) ||
(p == NULL) || (pt == NULL) ||
(Ap == NULL) || (ATpt == NULL) ||
(nzvalT == NULL) || (colidxT == NULL) ||
(rowptrT == NULL)){
printf("memory allocation failed\n");
return -1;
}
/*
Transpose A
*/
/* Allocate marker */
marker = (int *)calloc(dim, sizeof(int));
for (i = 0 ; i < dim ; i++) marker[i] = 0;
/* Get counts of each column of A, and set up column pointers */
for (i = 0 ; i < dim ; i++){
for (j = rowptr[i] ; j < rowptr[i+1]; j++) ++marker[colidx[j]];
}
rowptrT[0] = 0;
for (j = 0 ; j < dim ; j++) {
rowptrT[j+1] = rowptrT[j] + marker[j];
marker[j] = rowptrT[j];
}
/* Transfer the matrix into the compressed column storage. */
for (i = 0 ; i < dim ; i++) {
for (j = rowptr[i] ; j < rowptr[i+1] ; j++) {
col = colidx[j];
relpos = marker[col];
colidxT[relpos] = i;
nzvalT[relpos] = nzval[j];
++marker[col];
}
}
free(marker);
/*
Bi-Conjugate Gradient
*/
/* 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]; /* p0 = r0 */
rt[i] = r[i]; /* rt0 = r0 */
pt[i] = r[i]; /* pt0 = r0 */
rr += r[i] * r[i]; /* rr = r.r */
}
/* bcg iteration */
count = 0;
while(rr > tol){
for (i = 0 ; i < dim ; i++){
/* Ap = A*p */
Ap[i] = 0.0;
for (j = rowptr[i] ; j < rowptr[i+1] ; j++){
Ap[i] += nzval[j] * p[colidx[j]];
}
/* ATpt = A^T*pt */
ATpt[i] = 0.0;
for (j = rowptrT[i] ; j < rowptrT[i+1] ; j++){
ATpt[i] += nzvalT[j] * pt[colidxT[j]];
}
}
// alpha = r.r / pt.Ap
ptAp = 0.0;
for (i = 0 ; i < dim ; i++){
ptAp += pt[i] * Ap[i];
}
alpha = rr / ptAp;
//Beta
rr1 = 0.0;
for (i = 0 ; i < dim ; i++){
x[i] += alpha * p[i];
r[i] -= alpha * Ap[i];
rt[i] -= alpha * ATpt[i];
rr1 += r[i] * rt[i];
}
beta = rr1 / rr;
for (i = 0 ; i < dim ; i++){
p[i] = r[i] + beta * p[i];
pt[i] = rt[i] + beta * pt[i];
}
rr = rr1;
count++;
}
/* Deallocate Working Spaces */
free(r);
free(p);
free(Ap);
free(rt);
free(pt);
free(ATpt);
free(nzvalT);
free(colidxT);
free(rowptrT);
return count;
}