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

migrated content from ccs wiki

# Bcg.c

```/*
Large Scale Computing

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

/*
*/
/* 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;
}
```
Note: See TracWiki for help on using the wiki.