| | 1 | = Cg.c = |
| | 2 | {{{#!C |
| | 3 | /* |
| | 4 | Large Scale Computing |
| | 5 | |
| | 6 | solve Ax=b |
| | 7 | Conjugate Gradient Method |
| | 8 | General Sparse Matrix A |
| | 9 | |
| | 10 | solve_cg solve Ax=b with CG method |
| | 11 | this returns the number of iterations, or negative if failed. |
| | 12 | |
| | 13 | Note: A must be symmetric and positive-defined matrix |
| | 14 | */ |
| | 15 | #include<stdio.h> |
| | 16 | |
| | 17 | int solve_cg(int dim, /* dimension */ |
| | 18 | int nnz, /* # of non-zeros in the matrix */ |
| | 19 | double *nzval, /* array of nonzero values */ |
| | 20 | int *colidx, /* array of column indices of the nonzeros */ |
| | 21 | int *rowptr, /* the first column of nonzero in each row */ |
| | 22 | /* rowptr[j] stores the location of nzval[] and colidx[], */ |
| | 23 | /* which starts row j. This has dim+1 entry that is nnz */ |
| | 24 | double *x, /* solition */ |
| | 25 | double *b, /* right hand side */ |
| | 26 | double tol); /* tolerance */ |
| | 27 | }}} |
| | 28 | {{{#!C |
| | 29 | /* |
| | 30 | Large Scale Computing |
| | 31 | |
| | 32 | solve Ax=b |
| | 33 | Conjugate Gradient Method |
| | 34 | General Sparse Matrix A |
| | 35 | |
| | 36 | Note: A must be symmetric and positive-defined matrix |
| | 37 | */ |
| | 38 | #include<stdio.h> |
| | 39 | #include<stdlib.h> |
| | 40 | #include<math.h> |
| | 41 | |
| | 42 | int solve_cg(int dim, /* dimension */ |
| | 43 | int nnz, /* # of non-zeros in the matrix */ |
| | 44 | double *nzval, /* array of nonzero values */ |
| | 45 | int *colidx, /* array of column indices of the nonzeros */ |
| | 46 | int *rowptr, /* the first column of nonzero in each row */ |
| | 47 | double *x, /* solition */ |
| | 48 | double *b, /* right hand side */ |
| | 49 | double tol){ /* tolerance */ |
| | 50 | double *r; |
| | 51 | double *p; |
| | 52 | double *Ap; |
| | 53 | double rr,pAp,rr1,alpha,beta; |
| | 54 | int i,j,count; |
| | 55 | |
| | 56 | /* Allocate Working Spaces */ |
| | 57 | r = (double *)calloc(dim, sizeof(double)); |
| | 58 | p = (double *)calloc(dim, sizeof(double)); |
| | 59 | Ap = (double *)calloc(dim, sizeof(double)); |
| | 60 | if ((r == NULL) || (p == NULL) || (Ap == NULL)){ |
| | 61 | printf("memory allocation failed\n"); |
| | 62 | return -1; |
| | 63 | } |
| | 64 | |
| | 65 | /* compute r0 */ |
| | 66 | for (i = 0 ; i < dim ; i++){ |
| | 67 | r[i] = b[i]; |
| | 68 | for (j = rowptr[i] ; j < rowptr[i+1] ; j++){ |
| | 69 | r[i] -= nzval[j] * x[colidx[j]]; |
| | 70 | } |
| | 71 | } |
| | 72 | |
| | 73 | rr = 0.0; |
| | 74 | for (i = 0 ; i < dim ; i++){ |
| | 75 | p[i] = r[i]; /* p = r */ |
| | 76 | rr += r[i] * r[i]; /* rr = r.r */ |
| | 77 | } |
| | 78 | |
| | 79 | /* cg iteration */ |
| | 80 | count = 0; |
| | 81 | while(rr > tol){ |
| | 82 | // Ap = A*p |
| | 83 | for (i = 0 ; i < dim ; i++){ |
| | 84 | Ap[i] = 0.0; |
| | 85 | for (j = rowptr[i] ; j < rowptr[i+1] ; j++){ |
| | 86 | Ap[i] += nzval[j] * p[colidx[j]]; |
| | 87 | } |
| | 88 | } |
| | 89 | |
| | 90 | // alpha = r.r / p.Ap |
| | 91 | pAp = 0.0; |
| | 92 | for (i = 0 ; i < dim ; i++){ |
| | 93 | pAp += p[i] * Ap[i]; |
| | 94 | } |
| | 95 | alpha = rr / pAp; |
| | 96 | |
| | 97 | //Beta |
| | 98 | rr1 = 0.0; |
| | 99 | for (i = 0 ; i < dim ; i++){ |
| | 100 | x[i] += alpha * p[i]; |
| | 101 | r[i] -= alpha * Ap[i]; |
| | 102 | rr1 += r[i] * r[i]; |
| | 103 | } |
| | 104 | |
| | 105 | beta = rr1 / rr; |
| | 106 | for (i = 0 ; i < dim ; i++){ |
| | 107 | p[i] = r[i] + beta * p[i]; |
| | 108 | } |
| | 109 | |
| | 110 | rr = rr1; |
| | 111 | count++; |
| | 112 | } |
| | 113 | |
| | 114 | /* Deallocate Working Spaces */ |
| | 115 | free(r); |
| | 116 | free(p); |
| | 117 | free(Ap); |
| | 118 | |
| | 119 | return count; |
| | 120 | } |
| | 121 | }}} |
