= 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 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 #include #include 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; } }}}