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: PetscFunctionBegin;
24: PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATSEQAIJ, &flg));
25: PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Currently only handles SeqAIJ matrices");
27: PetscCall(MatGetLocalSize(pc->pmat, &cp->m, &cp->n));
28: PetscCheck(cp->m == cp->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Currently only for square matrices");
30: if (!cp->work) PetscCall(MatCreateVecs(pc->pmat, &cp->work, NULL));
31: if (!cp->d) PetscCall(PetscMalloc1(cp->n, &cp->d));
32: if (cp->a && pc->flag != SAME_NONZERO_PATTERN) {
33: PetscCall(PetscFree3(cp->a, cp->i, cp->j));
34: cp->a = NULL;
35: }
37: /* convert to column format */
38: if (!cp->a) PetscCall(PetscMalloc3(aij->nz, &cp->a, cp->n + 1, &cp->i, aij->nz, &cp->j));
39: PetscCall(PetscCalloc1(cp->n, &colcnt));
41: for (i = 0; i < aij->nz; i++) colcnt[aij->j[i]]++;
42: cp->i[0] = 0;
43: for (i = 0; i < cp->n; i++) cp->i[i + 1] = cp->i[i] + colcnt[i];
44: PetscCall(PetscArrayzero(colcnt, cp->n));
45: for (i = 0; i < cp->m; i++) { /* over rows */
46: for (j = aij->i[i]; j < aij->i[i + 1]; j++) { /* over columns in row */
47: cp->j[cp->i[aij->j[j]] + colcnt[aij->j[j]]] = i;
48: cp->a[cp->i[aij->j[j]] + colcnt[aij->j[j]]++] = aij->a[j];
49: }
50: }
51: PetscCall(PetscFree(colcnt));
53: /* compute sum of squares of each column d[] */
54: for (i = 0; i < cp->n; i++) { /* over columns */
55: cp->d[i] = 0.;
56: for (j = cp->i[i]; j < cp->i[i + 1]; j++) cp->d[i] += cp->a[j] * cp->a[j]; /* over rows in column */
57: cp->d[i] = 1.0 / cp->d[i];
58: }
59: PetscFunctionReturn(PETSC_SUCCESS);
60: }
62: static PetscErrorCode PCApply_CP(PC pc, Vec bb, Vec xx)
63: {
64: PC_CP *cp = (PC_CP *)pc->data;
65: PetscScalar *b, *x, xt;
66: PetscInt i, j;
68: PetscFunctionBegin;
69: PetscCall(VecCopy(bb, cp->work));
70: PetscCall(VecGetArray(cp->work, &b));
71: PetscCall(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: PetscCall(VecRestoreArray(cp->work, &b));
89: PetscCall(VecRestoreArray(xx, &x));
90: PetscFunctionReturn(PETSC_SUCCESS);
91: }
93: static PetscErrorCode PCReset_CP(PC pc)
94: {
95: PC_CP *cp = (PC_CP *)pc->data;
97: PetscFunctionBegin;
98: PetscCall(PetscFree(cp->d));
99: PetscCall(VecDestroy(&cp->work));
100: PetscCall(PetscFree3(cp->a, cp->i, cp->j));
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
104: static PetscErrorCode PCDestroy_CP(PC pc)
105: {
106: PC_CP *cp = (PC_CP *)pc->data;
108: PetscFunctionBegin;
109: PetscCall(PCReset_CP(pc));
110: PetscCall(PetscFree(cp->d));
111: PetscCall(PetscFree3(cp->a, cp->i, cp->j));
112: PetscCall(PetscFree(pc->data));
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: static PetscErrorCode PCSetFromOptions_CP(PC pc, PetscOptionItems *PetscOptionsObject)
117: {
118: PetscFunctionBegin;
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: /*MC
123: PCCP - a "column-projection" preconditioner
125: This is a terrible preconditioner and is not recommended, ever!
127: Loops over the entries of x computing dx_i (e_i is the unit vector in the ith direction) to
128: .vb
130: min || b - A(x + dx_i e_i ||_2
131: dx_i
133: That is, it changes a single entry of x to minimize the new residual norm.
134: Let A_i represent the ith column of A, then the minimization can be written as
136: min || r - (dx_i) A e_i ||_2
137: dx_i
138: or min || r - (dx_i) A_i ||_2
139: dx_i
141: take the derivative with respect to dx_i to obtain
142: dx_i = (A_i^T A_i)^(-1) A_i^T r
144: This algorithm can be thought of as Gauss-Seidel on the normal equations
145: .ve
147: Notes:
148: This procedure can also be done with block columns or any groups of columns
149: but this is not coded.
151: These "projections" can be done simultaneously for all columns (similar to Jacobi)
152: or sequentially (similar to Gauss-Seidel/SOR). This is only coded for SOR type.
154: This is related to, but not the same as "row projection" methods.
156: This is currently coded only for `MATSEQAIJ` matrices
158: Level: intermediate
160: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCJACOBI`, `PCSOR`
161: M*/
163: PETSC_EXTERN PetscErrorCode PCCreate_CP(PC pc)
164: {
165: PC_CP *cp;
167: PetscFunctionBegin;
168: PetscCall(PetscNew(&cp));
169: pc->data = (void *)cp;
171: pc->ops->apply = PCApply_CP;
172: pc->ops->applytranspose = PCApply_CP;
173: pc->ops->setup = PCSetUp_CP;
174: pc->ops->reset = PCReset_CP;
175: pc->ops->destroy = PCDestroy_CP;
176: pc->ops->setfromoptions = PCSetFromOptions_CP;
177: pc->ops->view = NULL;
178: pc->ops->applyrichardson = NULL;
179: PetscFunctionReturn(PETSC_SUCCESS);
180: }