Actual source code: eisen.c
1: /*
2: Defines a Eisenstat trick SSOR preconditioner. This uses about
3: %50 of the usual amount of floating point ops used for SSOR + Krylov
4: method. But it requires actually solving the preconditioned problem
5: with both left and right preconditioning.
6: */
7: #include <petsc/private/pcimpl.h>
9: typedef struct {
10: Mat shell, A;
11: Vec b[2], diag; /* temporary storage for true right-hand side */
12: PetscReal omega;
13: PetscBool usediag; /* indicates preconditioner should include diagonal scaling*/
14: } PC_Eisenstat;
16: static PetscErrorCode PCMult_Eisenstat(Mat mat, Vec b, Vec x)
17: {
18: PC pc;
19: PC_Eisenstat *eis;
21: PetscFunctionBegin;
22: PetscCall(MatShellGetContext(mat, &pc));
23: eis = (PC_Eisenstat *)pc->data;
24: PetscCall(MatSOR(eis->A, b, eis->omega, SOR_EISENSTAT, 0.0, 1, 1, x));
25: PetscCall(MatFactorGetError(eis->A, (MatFactorError *)&pc->failedreason));
26: PetscFunctionReturn(PETSC_SUCCESS);
27: }
29: static PetscErrorCode PCNorm_Eisenstat(Mat mat, NormType type, PetscReal *nrm)
30: {
31: PC pc;
32: PC_Eisenstat *eis;
34: PetscFunctionBegin;
35: PetscCall(MatShellGetContext(mat, &pc));
36: eis = (PC_Eisenstat *)pc->data;
37: PetscCall(MatNorm(eis->A, type, nrm));
38: PetscFunctionReturn(PETSC_SUCCESS);
39: }
41: static PetscErrorCode PCApply_Eisenstat(PC pc, Vec x, Vec y)
42: {
43: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
44: PetscBool hasop;
46: PetscFunctionBegin;
47: if (eis->usediag) {
48: PetscCall(MatHasOperation(pc->pmat, MATOP_MULT_DIAGONAL_BLOCK, &hasop));
49: if (hasop) {
50: PetscCall(MatMultDiagonalBlock(pc->pmat, x, y));
51: } else {
52: PetscCall(VecPointwiseMult(y, x, eis->diag));
53: }
54: } else PetscCall(VecCopy(x, y));
55: PetscFunctionReturn(PETSC_SUCCESS);
56: }
58: static PetscErrorCode PCApplyTranspose_Eisenstat(PC pc, Vec x, Vec y)
59: {
60: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
61: PetscBool hasop, set, sym;
63: PetscFunctionBegin;
64: PetscCall(MatIsSymmetricKnown(eis->A, &set, &sym));
65: PetscCheck(set && sym, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Can only apply transpose of Eisenstat if matrix is symmetric");
66: if (eis->usediag) {
67: PetscCall(MatHasOperation(pc->pmat, MATOP_MULT_DIAGONAL_BLOCK, &hasop));
68: if (hasop) {
69: PetscCall(MatMultDiagonalBlock(pc->pmat, x, y));
70: } else {
71: PetscCall(VecPointwiseMult(y, x, eis->diag));
72: }
73: } else PetscCall(VecCopy(x, y));
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
77: static PetscErrorCode PCPreSolve_Eisenstat(PC pc, KSP ksp, Vec b, Vec x)
78: {
79: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
80: PetscBool nonzero;
82: PetscFunctionBegin;
83: if (pc->presolvedone < 2) {
84: PetscCheck(pc->mat == pc->pmat, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot have different mat and pmat");
85: /* swap shell matrix and true matrix */
86: eis->A = pc->mat;
87: pc->mat = eis->shell;
88: }
90: if (!eis->b[pc->presolvedone - 1]) { PetscCall(VecDuplicate(b, &eis->b[pc->presolvedone - 1])); }
92: /* if nonzero initial guess, modify x */
93: PetscCall(KSPGetInitialGuessNonzero(ksp, &nonzero));
94: if (nonzero) {
95: PetscCall(VecCopy(x, eis->b[pc->presolvedone - 1]));
96: PetscCall(MatSOR(eis->A, eis->b[pc->presolvedone - 1], eis->omega, SOR_APPLY_UPPER, 0.0, 1, 1, x));
97: PetscCall(MatFactorGetError(eis->A, (MatFactorError *)&pc->failedreason));
98: }
100: /* save true b, other option is to swap pointers */
101: PetscCall(VecCopy(b, eis->b[pc->presolvedone - 1]));
103: /* modify b by (L + D/omega)^{-1} */
104: PetscCall(MatSOR(eis->A, eis->b[pc->presolvedone - 1], eis->omega, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_FORWARD_SWEEP), 0.0, 1, 1, b));
105: PetscCall(MatFactorGetError(eis->A, (MatFactorError *)&pc->failedreason));
106: PetscFunctionReturn(PETSC_SUCCESS);
107: }
109: static PetscErrorCode PCPostSolve_Eisenstat(PC pc, KSP ksp, Vec b, Vec x)
110: {
111: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
113: PetscFunctionBegin;
114: /* get back true b */
115: PetscCall(VecCopy(eis->b[pc->presolvedone], b));
117: /* modify x by (U + D/omega)^{-1} */
118: PetscCall(VecCopy(x, eis->b[pc->presolvedone]));
119: PetscCall(MatSOR(eis->A, eis->b[pc->presolvedone], eis->omega, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_BACKWARD_SWEEP), 0.0, 1, 1, x));
120: PetscCall(MatFactorGetError(eis->A, (MatFactorError *)&pc->failedreason));
121: if (!pc->presolvedone) pc->mat = eis->A;
122: PetscFunctionReturn(PETSC_SUCCESS);
123: }
125: static PetscErrorCode PCReset_Eisenstat(PC pc)
126: {
127: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
129: PetscFunctionBegin;
130: PetscCall(VecDestroy(&eis->b[0]));
131: PetscCall(VecDestroy(&eis->b[1]));
132: PetscCall(MatDestroy(&eis->shell));
133: PetscCall(VecDestroy(&eis->diag));
134: PetscFunctionReturn(PETSC_SUCCESS);
135: }
137: static PetscErrorCode PCDestroy_Eisenstat(PC pc)
138: {
139: PetscFunctionBegin;
140: PetscCall(PCReset_Eisenstat(pc));
141: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetOmega_C", NULL));
142: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetNoDiagonalScaling_C", NULL));
143: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetOmega_C", NULL));
144: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetNoDiagonalScaling_C", NULL));
145: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", NULL));
146: PetscCall(PetscFree(pc->data));
147: PetscFunctionReturn(PETSC_SUCCESS);
148: }
150: static PetscErrorCode PCSetFromOptions_Eisenstat(PC pc, PetscOptionItems *PetscOptionsObject)
151: {
152: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
153: PetscBool set, flg;
154: PetscReal omega;
156: PetscFunctionBegin;
157: PetscOptionsHeadBegin(PetscOptionsObject, "Eisenstat SSOR options");
158: PetscCall(PetscOptionsReal("-pc_eisenstat_omega", "Relaxation factor 0 < omega < 2", "PCEisenstatSetOmega", eis->omega, &omega, &flg));
159: if (flg) PetscCall(PCEisenstatSetOmega(pc, omega));
160: PetscCall(PetscOptionsBool("-pc_eisenstat_no_diagonal_scaling", "Do not use standard diagonal scaling", "PCEisenstatSetNoDiagonalScaling", eis->usediag ? PETSC_FALSE : PETSC_TRUE, &flg, &set));
161: if (set) PetscCall(PCEisenstatSetNoDiagonalScaling(pc, flg));
162: PetscOptionsHeadEnd();
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: static PetscErrorCode PCView_Eisenstat(PC pc, PetscViewer viewer)
167: {
168: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
169: PetscBool iascii;
171: PetscFunctionBegin;
172: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
173: if (iascii) {
174: PetscCall(PetscViewerASCIIPrintf(viewer, " omega = %g\n", (double)eis->omega));
175: if (eis->usediag) {
176: PetscCall(PetscViewerASCIIPrintf(viewer, " Using diagonal scaling (default)\n"));
177: } else {
178: PetscCall(PetscViewerASCIIPrintf(viewer, " Not using diagonal scaling\n"));
179: }
180: }
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
184: static PetscErrorCode PCSetUp_Eisenstat(PC pc)
185: {
186: PetscInt M, N, m, n;
187: PetscBool set, sym;
188: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
190: PetscFunctionBegin;
191: if (!pc->setupcalled) {
192: PetscCall(MatGetSize(pc->mat, &M, &N));
193: PetscCall(MatGetLocalSize(pc->mat, &m, &n));
194: PetscCall(MatIsSymmetricKnown(pc->mat, &set, &sym));
195: PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &eis->shell));
196: PetscCall(MatSetSizes(eis->shell, m, n, M, N));
197: PetscCall(MatSetType(eis->shell, MATSHELL));
198: PetscCall(MatSetUp(eis->shell));
199: PetscCall(MatShellSetContext(eis->shell, pc));
200: PetscCall(MatShellSetOperation(eis->shell, MATOP_MULT, (void (*)(void))PCMult_Eisenstat));
201: if (set && sym) PetscCall(MatShellSetOperation(eis->shell, MATOP_MULT_TRANSPOSE, (void (*)(void))PCMult_Eisenstat));
202: PetscCall(MatShellSetOperation(eis->shell, MATOP_NORM, (void (*)(void))PCNorm_Eisenstat));
203: }
204: if (!eis->usediag) PetscFunctionReturn(PETSC_SUCCESS);
205: if (!pc->setupcalled) { PetscCall(MatCreateVecs(pc->pmat, &eis->diag, NULL)); }
206: PetscCall(MatGetDiagonal(pc->pmat, eis->diag));
207: PetscFunctionReturn(PETSC_SUCCESS);
208: }
210: /* --------------------------------------------------------------------*/
212: static PetscErrorCode PCEisenstatSetOmega_Eisenstat(PC pc, PetscReal omega)
213: {
214: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
216: PetscFunctionBegin;
217: PetscCheck(omega > 0.0 && omega < 2.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Relaxation out of range");
218: eis->omega = omega;
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: static PetscErrorCode PCEisenstatSetNoDiagonalScaling_Eisenstat(PC pc, PetscBool flg)
223: {
224: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
226: PetscFunctionBegin;
227: eis->usediag = flg;
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: static PetscErrorCode PCEisenstatGetOmega_Eisenstat(PC pc, PetscReal *omega)
232: {
233: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
235: PetscFunctionBegin;
236: *omega = eis->omega;
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: static PetscErrorCode PCEisenstatGetNoDiagonalScaling_Eisenstat(PC pc, PetscBool *flg)
241: {
242: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
244: PetscFunctionBegin;
245: *flg = eis->usediag;
246: PetscFunctionReturn(PETSC_SUCCESS);
247: }
249: /*@
250: PCEisenstatSetOmega - Sets the SSOR relaxation coefficient, omega,
251: to use with Eisenstat's trick (where omega = 1.0 by default)
253: Logically Collective
255: Input Parameters:
256: + pc - the preconditioner context
257: - omega - relaxation coefficient (0 < omega < 2)
259: Options Database Key:
260: . -pc_eisenstat_omega <omega> - Sets omega
262: Level: intermediate
264: Notes:
265: The Eisenstat trick implementation of SSOR requires about 50% of the
266: usual amount of floating point operations used for SSOR + Krylov method;
267: however, the preconditioned problem must be solved with both left
268: and right preconditioning.
270: To use SSOR without the Eisenstat trick, employ the `PCSOR` preconditioner,
271: which can be chosen with the database options `-pc_type sor -pc_sor_symmetric`
273: .seealso: [](ch_ksp), `PCSORSetOmega()`, `PCEISENSTAT`
274: @*/
275: PetscErrorCode PCEisenstatSetOmega(PC pc, PetscReal omega)
276: {
277: PetscFunctionBegin;
280: PetscTryMethod(pc, "PCEisenstatSetOmega_C", (PC, PetscReal), (pc, omega));
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }
284: /*@
285: PCEisenstatSetNoDiagonalScaling - Causes the Eisenstat preconditioner, `PCEISENSTAT`
286: not to do additional diagonal preconditioning. For matrices with a constant
287: along the diagonal, this may save a small amount of work.
289: Logically Collective
291: Input Parameters:
292: + pc - the preconditioner context
293: - flg - `PETSC_TRUE` turns off diagonal scaling inside the algorithm
295: Options Database Key:
296: . -pc_eisenstat_no_diagonal_scaling - Activates `PCEisenstatSetNoDiagonalScaling()`
298: Level: intermediate
300: Note:
301: If you use the `KSPSetDiagonalScaling()` or -ksp_diagonal_scale option then you will
302: likely want to use this routine since it will save you some unneeded flops.
304: .seealso: [](ch_ksp), `PCEisenstatSetOmega()`, `PCEISENSTAT`
305: @*/
306: PetscErrorCode PCEisenstatSetNoDiagonalScaling(PC pc, PetscBool flg)
307: {
308: PetscFunctionBegin;
310: PetscTryMethod(pc, "PCEisenstatSetNoDiagonalScaling_C", (PC, PetscBool), (pc, flg));
311: PetscFunctionReturn(PETSC_SUCCESS);
312: }
314: /*@
315: PCEisenstatGetOmega - Gets the SSOR relaxation coefficient, omega,
316: to use with Eisenstat's trick (where omega = 1.0 by default).
318: Logically Collective
320: Input Parameter:
321: . pc - the preconditioner context
323: Output Parameter:
324: . omega - relaxation coefficient (0 < omega < 2)
326: Options Database Key:
327: . -pc_eisenstat_omega <omega> - Sets omega
329: Notes:
330: The Eisenstat trick implementation of SSOR requires about 50% of the
331: usual amount of floating point operations used for SSOR + Krylov method;
332: however, the preconditioned problem must be solved with both left
333: and right preconditioning.
335: To use SSOR without the Eisenstat trick, employ the PCSOR preconditioner,
336: which can be chosen with the database options `-pc_type sor -pc_sor_symmetric`
338: Level: intermediate
340: .seealso: [](ch_ksp), `PCEISENSTAT`, `PCSORGetOmega()`, `PCEisenstatSetOmega()`
341: @*/
342: PetscErrorCode PCEisenstatGetOmega(PC pc, PetscReal *omega)
343: {
344: PetscFunctionBegin;
346: PetscUseMethod(pc, "PCEisenstatGetOmega_C", (PC, PetscReal *), (pc, omega));
347: PetscFunctionReturn(PETSC_SUCCESS);
348: }
350: /*@
351: PCEisenstatGetNoDiagonalScaling - Tells if the Eisenstat preconditioner
352: not to do additional diagonal preconditioning. For matrices with a constant
353: along the diagonal, this may save a small amount of work.
355: Logically Collective
357: Input Parameter:
358: . pc - the preconditioner context
360: Output Parameter:
361: . flg - `PETSC_TRUE` means there is no diagonal scaling applied
363: Options Database Key:
364: . -pc_eisenstat_no_diagonal_scaling - Activates `PCEisenstatSetNoDiagonalScaling()`
366: Level: intermediate
368: Note:
369: If you use the KSPSetDiagonalScaling() or -ksp_diagonal_scale option then you will
370: likely want to use this routine since it will save you some unneeded flops.
372: .seealso: , `PCEISENSTAT`, `PCEisenstatGetOmega()`
373: @*/
374: PetscErrorCode PCEisenstatGetNoDiagonalScaling(PC pc, PetscBool *flg)
375: {
376: PetscFunctionBegin;
378: PetscUseMethod(pc, "PCEisenstatGetNoDiagonalScaling_C", (PC, PetscBool *), (pc, flg));
379: PetscFunctionReturn(PETSC_SUCCESS);
380: }
382: static PetscErrorCode PCPreSolveChangeRHS_Eisenstat(PC pc, PetscBool *change)
383: {
384: PetscFunctionBegin;
385: *change = PETSC_TRUE;
386: PetscFunctionReturn(PETSC_SUCCESS);
387: }
389: /*MC
390: PCEISENSTAT - An implementation of SSOR (symmetric successive over relaxation, symmetric Gauss-Seidel)
391: preconditioning that incorporates Eisenstat's trick to reduce the amount of computation needed.
393: Options Database Keys:
394: + -pc_eisenstat_omega <omega> - Sets omega
395: - -pc_eisenstat_no_diagonal_scaling - Activates `PCEisenstatSetNoDiagonalScaling()`
397: Level: beginner
399: Notes:
400: Only implemented for the `MATAIJ` matrix format.
402: Not a true parallel SOR, in parallel this implementation corresponds to block Jacobi with SOR on each block.
404: Developer Note:
405: Since this algorithm runs the Krylov method on a transformed linear system the implementation provides `PCPreSolve()` and `PCPostSolve()`
406: routines that `KSP` uses to set up the transformed linear system.
408: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCEisenstatGetOmega()`,
409: `PCEisenstatSetNoDiagonalScaling()`, `PCEisenstatSetOmega()`, `PCSOR`
410: M*/
412: PETSC_EXTERN PetscErrorCode PCCreate_Eisenstat(PC pc)
413: {
414: PC_Eisenstat *eis;
416: PetscFunctionBegin;
417: PetscCall(PetscNew(&eis));
419: pc->ops->apply = PCApply_Eisenstat;
420: pc->ops->applytranspose = PCApplyTranspose_Eisenstat;
421: pc->ops->presolve = PCPreSolve_Eisenstat;
422: pc->ops->postsolve = PCPostSolve_Eisenstat;
423: pc->ops->applyrichardson = NULL;
424: pc->ops->setfromoptions = PCSetFromOptions_Eisenstat;
425: pc->ops->destroy = PCDestroy_Eisenstat;
426: pc->ops->reset = PCReset_Eisenstat;
427: pc->ops->view = PCView_Eisenstat;
428: pc->ops->setup = PCSetUp_Eisenstat;
430: pc->data = eis;
431: eis->omega = 1.0;
432: eis->b[0] = NULL;
433: eis->b[1] = NULL;
434: eis->diag = NULL;
435: eis->usediag = PETSC_TRUE;
437: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetOmega_C", PCEisenstatSetOmega_Eisenstat));
438: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetNoDiagonalScaling_C", PCEisenstatSetNoDiagonalScaling_Eisenstat));
439: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetOmega_C", PCEisenstatGetOmega_Eisenstat));
440: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetNoDiagonalScaling_C", PCEisenstatGetNoDiagonalScaling_Eisenstat));
441: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", PCPreSolveChangeRHS_Eisenstat));
442: PetscFunctionReturn(PETSC_SUCCESS);
443: }