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) PetscMalloc3(aij->nz, &cp->a, cp->n + 1, &cp->i, aij->nz, &cp->j);
38: PetscCalloc1(cp->n, &colcnt);
40: for (i = 0; i < aij->nz; i++) colcnt[aij->j[i]]++;
41: cp->i[0] = 0;
42: for (i = 0; i < cp->n; i++) cp->i[i + 1] = cp->i[i] + colcnt[i];
43: PetscArrayzero(colcnt, cp->n);
44: for (i = 0; i < cp->m; i++) { /* over rows */
45: for (j = aij->i[i]; j < aij->i[i + 1]; j++) { /* over columns in row */
46: cp->j[cp->i[aij->j[j]] + colcnt[aij->j[j]]] = i;
47: cp->a[cp->i[aij->j[j]] + colcnt[aij->j[j]]++] = aij->a[j];
48: }
49: }
50: PetscFree(colcnt);
52: /* compute sum of squares of each column d[] */
53: for (i = 0; i < cp->n; i++) { /* over columns */
54: cp->d[i] = 0.;
55: for (j = cp->i[i]; j < cp->i[i + 1]; j++) cp->d[i] += cp->a[j] * cp->a[j]; /* over rows in column */
56: cp->d[i] = 1.0 / cp->d[i];
57: }
58: return 0;
59: }
61: static PetscErrorCode PCApply_CP(PC pc, Vec bb, Vec xx)
62: {
63: PC_CP *cp = (PC_CP *)pc->data;
64: PetscScalar *b, *x, xt;
65: PetscInt i, j;
67: VecCopy(bb, cp->work);
68: VecGetArray(cp->work, &b);
69: VecGetArray(xx, &x);
71: for (i = 0; i < cp->n; i++) { /* over columns */
72: xt = 0.;
73: for (j = cp->i[i]; j < cp->i[i + 1]; j++) xt += cp->a[j] * b[cp->j[j]]; /* over rows in column */
74: xt *= cp->d[i];
75: x[i] = xt;
76: 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*/
77: }
78: for (i = cp->n - 1; i > -1; i--) { /* over columns */
79: xt = 0.;
80: for (j = cp->i[i]; j < cp->i[i + 1]; j++) xt += cp->a[j] * b[cp->j[j]]; /* over rows in column */
81: xt *= cp->d[i];
82: x[i] = xt;
83: 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*/
84: }
86: VecRestoreArray(cp->work, &b);
87: VecRestoreArray(xx, &x);
88: return 0;
89: }
91: static PetscErrorCode PCReset_CP(PC pc)
92: {
93: PC_CP *cp = (PC_CP *)pc->data;
95: PetscFree(cp->d);
96: VecDestroy(&cp->work);
97: PetscFree3(cp->a, cp->i, cp->j);
98: return 0;
99: }
101: static PetscErrorCode PCDestroy_CP(PC pc)
102: {
103: PC_CP *cp = (PC_CP *)pc->data;
105: PCReset_CP(pc);
106: PetscFree(cp->d);
107: PetscFree3(cp->a, cp->i, cp->j);
108: PetscFree(pc->data);
109: return 0;
110: }
112: static PetscErrorCode PCSetFromOptions_CP(PC pc, PetscOptionItems *PetscOptionsObject)
113: {
114: return 0;
115: }
117: /*MC
118: PCCP - a "column-projection" preconditioner
120: This is a terrible preconditioner and is not recommended, ever!
122: Loops over the entries of x computing dx_i (e_i is the unit vector in the ith direction) to
123: .vb
125: min || b - A(x + dx_i e_i ||_2
126: dx_i
128: That is, it changes a single entry of x to minimize the new residual norm.
129: Let A_i represent the ith column of A, then the minimization can be written as
131: min || r - (dx_i) A e_i ||_2
132: dx_i
133: or min || r - (dx_i) A_i ||_2
134: dx_i
136: take the derivative with respect to dx_i to obtain
137: dx_i = (A_i^T A_i)^(-1) A_i^T r
139: This algorithm can be thought of as Gauss-Seidel on the normal equations
140: .ve
142: Notes:
143: This procedure 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 `MATSEQAIJ` matrices
153: Level: intermediate
155: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCJACOBI`, `PCSOR`
156: M*/
158: PETSC_EXTERN PetscErrorCode PCCreate_CP(PC pc)
159: {
160: PC_CP *cp;
162: PetscNew(&cp);
163: pc->data = (void *)cp;
165: pc->ops->apply = PCApply_CP;
166: pc->ops->applytranspose = PCApply_CP;
167: pc->ops->setup = PCSetUp_CP;
168: pc->ops->reset = PCReset_CP;
169: pc->ops->destroy = PCDestroy_CP;
170: pc->ops->setfromoptions = PCSetFromOptions_CP;
171: pc->ops->view = NULL;
172: pc->ops->applyrichardson = NULL;
173: return 0;
174: }