Actual source code: cp.c
2: #include <petsc/private/pcimpl.h>
3: #include <../src/mat/impls/aij/seq/aij.h>
5: /*
6: Private context (data structure) for the CP preconditioner.
7: */
8: typedef struct {
9: PetscInt n,m;
10: Vec work;
11: PetscScalar *d; /* sum of squares of each column */
12: PetscScalar *a; /* non-zeros by column */
13: PetscInt *i,*j; /* offsets of nonzeros by column, non-zero indices by column */
14: } PC_CP;
16: static PetscErrorCode PCSetUp_CP(PC pc)
17: {
18: PC_CP *cp = (PC_CP*)pc->data;
19: PetscInt i,j,*colcnt;
20: PetscBool flg;
21: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)pc->pmat->data;
23: PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQAIJ,&flg);
26: MatGetLocalSize(pc->pmat,&cp->m,&cp->n);
29: if (!cp->work) MatCreateVecs(pc->pmat,&cp->work,NULL);
30: if (!cp->d) PetscMalloc1(cp->n,&cp->d);
31: if (cp->a && pc->flag != SAME_NONZERO_PATTERN) {
32: PetscFree3(cp->a,cp->i,cp->j);
33: cp->a = NULL;
34: }
36: /* convert to column format */
37: if (!cp->a) {
38: PetscMalloc3(aij->nz,&cp->a,cp->n+1,&cp->i,aij->nz,&cp->j);
39: }
40: PetscCalloc1(cp->n,&colcnt);
42: for (i=0; i<aij->nz; i++) colcnt[aij->j[i]]++;
43: cp->i[0] = 0;
44: for (i=0; i<cp->n; i++) cp->i[i+1] = cp->i[i] + colcnt[i];
45: PetscArrayzero(colcnt,cp->n);
46: for (i=0; i<cp->m; i++) { /* over rows */
47: for (j=aij->i[i]; j<aij->i[i+1]; j++) { /* over columns in row */
48: cp->j[cp->i[aij->j[j]]+colcnt[aij->j[j]]] = i;
49: cp->a[cp->i[aij->j[j]]+colcnt[aij->j[j]]++] = aij->a[j];
50: }
51: }
52: PetscFree(colcnt);
54: /* compute sum of squares of each column d[] */
55: for (i=0; i<cp->n; i++) { /* over columns */
56: cp->d[i] = 0.;
57: for (j=cp->i[i]; j<cp->i[i+1]; j++) cp->d[i] += cp->a[j]*cp->a[j]; /* over rows in column */
58: cp->d[i] = 1.0/cp->d[i];
59: }
60: return 0;
61: }
62: /* -------------------------------------------------------------------------- */
63: static PetscErrorCode PCApply_CP(PC pc,Vec bb,Vec xx)
64: {
65: PC_CP *cp = (PC_CP*)pc->data;
66: PetscScalar *b,*x,xt;
67: PetscInt i,j;
69: VecCopy(bb,cp->work);
70: VecGetArray(cp->work,&b);
71: VecGetArray(xx,&x);
73: for (i=0; i<cp->n; i++) { /* over columns */
74: xt = 0.;
75: for (j=cp->i[i]; j<cp->i[i+1]; j++) xt += cp->a[j]*b[cp->j[j]]; /* over rows in column */
76: xt *= cp->d[i];
77: x[i] = xt;
78: for (j=cp->i[i]; j<cp->i[i+1]; j++) b[cp->j[j]] -= xt*cp->a[j]; /* over rows in column updating b*/
79: }
80: for (i=cp->n-1; i>-1; i--) { /* over columns */
81: xt = 0.;
82: for (j=cp->i[i]; j<cp->i[i+1]; j++) xt += cp->a[j]*b[cp->j[j]]; /* over rows in column */
83: xt *= cp->d[i];
84: x[i] = xt;
85: for (j=cp->i[i]; j<cp->i[i+1]; j++) b[cp->j[j]] -= xt*cp->a[j]; /* over rows in column updating b*/
86: }
88: VecRestoreArray(cp->work,&b);
89: VecRestoreArray(xx,&x);
90: return 0;
91: }
92: /* -------------------------------------------------------------------------- */
93: static PetscErrorCode PCReset_CP(PC pc)
94: {
95: PC_CP *cp = (PC_CP*)pc->data;
97: PetscFree(cp->d);
98: VecDestroy(&cp->work);
99: PetscFree3(cp->a,cp->i,cp->j);
100: return 0;
101: }
103: static PetscErrorCode PCDestroy_CP(PC pc)
104: {
105: PC_CP *cp = (PC_CP*)pc->data;
107: PCReset_CP(pc);
108: PetscFree(cp->d);
109: PetscFree3(cp->a,cp->i,cp->j);
110: PetscFree(pc->data);
111: return 0;
112: }
114: static PetscErrorCode PCSetFromOptions_CP(PetscOptionItems *PetscOptionsObject,PC pc)
115: {
116: return 0;
117: }
119: /*MC
120: PCCP - a "column-projection" preconditioner
122: This is a terrible preconditioner and is not recommended, ever!
124: Loops over the entries of x computing dx_i (e_i is the unit vector in the ith direction) to
125: $
126: $ min || b - A(x + dx_i e_i ||_2
127: $ dx_i
128: $
129: $ That is, it changes a single entry of x to minimize the new residual norm.
130: $ Let A_i represent the ith column of A, then the minimization can be written as
131: $
132: $ min || r - (dx_i) A e_i ||_2
133: $ dx_i
134: $ or min || r - (dx_i) A_i ||_2
135: $ dx_i
136: $
137: $ take the derivative with respect to dx_i to obtain
138: $ dx_i = (A_i^T A_i)^(-1) A_i^T r
139: $
140: $ This algorithm can be thought of as Gauss-Seidel on the normal equations
142: Notes:
143: This proceedure can also be done with block columns or any groups of columns
144: but this is not coded.
146: These "projections" can be done simultaneously for all columns (similar to Jacobi)
147: or sequentially (similar to Gauss-Seidel/SOR). This is only coded for SOR type.
149: This is related to, but not the same as "row projection" methods.
151: This is currently coded only for SeqAIJ matrices in sequential (SOR) form.
153: Level: intermediate
155: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCJACOBI, PCSOR
157: M*/
159: PETSC_EXTERN PetscErrorCode PCCreate_CP(PC pc)
160: {
161: PC_CP *cp;
163: PetscNewLog(pc,&cp);
164: pc->data = (void*)cp;
166: pc->ops->apply = PCApply_CP;
167: pc->ops->applytranspose = PCApply_CP;
168: pc->ops->setup = PCSetUp_CP;
169: pc->ops->reset = PCReset_CP;
170: pc->ops->destroy = PCDestroy_CP;
171: pc->ops->setfromoptions = PCSetFromOptions_CP;
172: pc->ops->view = NULL;
173: pc->ops->applyrichardson = NULL;
174: return 0;
175: }