= Bcg.c = {{{#!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 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 */ }}} {{{#!C /* Large Scale Computing solve Ax=b Bi-Conjugate Gradient Method General Sparse Matrix A */ #include #include #include 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; } }}}