Actual source code: cp.c
1: #include <petsc/private/pcimpl.h>
2: #include <../src/mat/impls/aij/seq/aij.h>
4: /*
5: Private context (data structure) for the CP preconditioner.
6: */
7: typedef struct {
8: PetscInt n, m;
9: Vec work;
10: PetscScalar *d; /* sum of squares of each column */
11: PetscScalar *a; /* non-zeros by column */
12: PetscInt *i, *j; /* offsets of nonzeros by column, non-zero indices by column */
13: } PC_CP;
15: static PetscErrorCode PCSetUp_CP(PC pc)
16: {
17: PC_CP *cp = (PC_CP *)pc->data;
18: PetscInt i, j, *colcnt;
19: PetscBool flg;
20: Mat_SeqAIJ *aij = (Mat_SeqAIJ *)pc->pmat->data;
22: PetscFunctionBegin;
23: PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATSEQAIJ, &flg));
24: PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Currently only handles SeqAIJ matrices");
26: PetscCall(MatGetLocalSize(pc->pmat, &cp->m, &cp->n));
27: PetscCheck(cp->m == cp->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Currently only for square matrices");
29: if (!cp->work) PetscCall(MatCreateVecs(pc->pmat, &cp->work, NULL));
30: if (!cp->d) PetscCall(PetscMalloc1(cp->n, &cp->d));
31: if (cp->a && pc->flag != SAME_NONZERO_PATTERN) {
32: PetscCall(PetscFree3(cp->a, cp->i, cp->j));
33: cp->a = NULL;
34: }
36: /* convert to column format */
37: if (!cp->a) PetscCall(PetscMalloc3(aij->nz, &cp->a, cp->n + 1, &cp->i, aij->nz, &cp->j));
38: PetscCall(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: PetscCall(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: PetscCall(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: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
68: PetscCall(VecCopy(bb, cp->work));
69: PetscCall(VecGetArray(cp->work, &b));
70: PetscCall(VecGetArray(xx, &x));
72: for (i = 0; i < cp->n; i++) { /* over columns */
73: xt = 0.;
74: for (j = cp->i[i]; j < cp->i[i + 1]; j++) xt += cp->a[j] * b[cp->j[j]]; /* over rows in column */
75: xt *= cp->d[i];
76: x[i] = xt;
77: 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*/
78: }
79: for (i = cp->n - 1; i > -1; i--) { /* over columns */
80: xt = 0.;
81: for (j = cp->i[i]; j < cp->i[i + 1]; j++) xt += cp->a[j] * b[cp->j[j]]; /* over rows in column */
82: xt *= cp->d[i];
83: x[i] = xt;
84: 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*/
85: }
87: PetscCall(VecRestoreArray(cp->work, &b));
88: PetscCall(VecRestoreArray(xx, &x));
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
92: static PetscErrorCode PCReset_CP(PC pc)
93: {
94: PC_CP *cp = (PC_CP *)pc->data;
96: PetscFunctionBegin;
97: PetscCall(PetscFree(cp->d));
98: PetscCall(VecDestroy(&cp->work));
99: PetscCall(PetscFree3(cp->a, cp->i, cp->j));
100: PetscFunctionReturn(PETSC_SUCCESS);
101: }
103: static PetscErrorCode PCDestroy_CP(PC pc)
104: {
105: PC_CP *cp = (PC_CP *)pc->data;
107: PetscFunctionBegin;
108: PetscCall(PCReset_CP(pc));
109: PetscCall(PetscFree(cp->d));
110: PetscCall(PetscFree3(cp->a, cp->i, cp->j));
111: PetscCall(PetscFree(pc->data));
112: PetscFunctionReturn(PETSC_SUCCESS);
113: }
115: static PetscErrorCode PCSetFromOptions_CP(PC pc, PetscOptionItems *PetscOptionsObject)
116: {
117: PetscFunctionBegin;
118: PetscFunctionReturn(PETSC_SUCCESS);
119: }
121: /*MC
122: PCCP - a "column-projection" preconditioner
124: This is a terrible preconditioner and is not recommended, ever!
126: Loops over the entries of x computing dx_i (e_i is the unit vector in the ith direction) to
127: .vb
129: min || b - A(x + dx_i e_i ||_2
130: dx_i
132: That is, it changes a single entry of x to minimize the new residual norm.
133: Let A_i represent the ith column of A, then the minimization can be written as
135: min || r - (dx_i) A e_i ||_2
136: dx_i
137: or min || r - (dx_i) A_i ||_2
138: dx_i
140: take the derivative with respect to dx_i to obtain
141: dx_i = (A_i^T A_i)^(-1) A_i^T r
143: This algorithm can be thought of as Gauss-Seidel on the normal equations
144: .ve
146: Notes:
147: This procedure can also be done with block columns or any groups of columns
148: but this is not coded.
150: These "projections" can be done simultaneously for all columns (similar to Jacobi)
151: or sequentially (similar to Gauss-Seidel/SOR). This is only coded for SOR type.
153: This is related to, but not the same as "row projection" methods.
155: This is currently coded only for `MATSEQAIJ` matrices
157: Level: intermediate
159: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PCJACOBI`, `PCSOR`
160: M*/
162: PETSC_EXTERN PetscErrorCode PCCreate_CP(PC pc)
163: {
164: PC_CP *cp;
166: PetscFunctionBegin;
167: PetscCall(PetscNew(&cp));
168: pc->data = (void *)cp;
170: pc->ops->apply = PCApply_CP;
171: pc->ops->applytranspose = PCApply_CP;
172: pc->ops->setup = PCSetUp_CP;
173: pc->ops->reset = PCReset_CP;
174: pc->ops->destroy = PCDestroy_CP;
175: pc->ops->setfromoptions = PCSetFromOptions_CP;
176: pc->ops->view = NULL;
177: pc->ops->applyrichardson = NULL;
178: PetscFunctionReturn(PETSC_SUCCESS);
179: }