Actual source code: eisen.c
2: /*
3: Defines a Eisenstat trick SSOR preconditioner. This uses about
4: %50 of the usual amount of floating point ops used for SSOR + Krylov
5: method. But it requires actually solving the preconditioned problem
6: with both left and right preconditioning.
7: */
8: #include <petsc/private/pcimpl.h>
10: typedef struct {
11: Mat shell, A;
12: Vec b[2], diag; /* temporary storage for true right hand side */
13: PetscReal omega;
14: PetscBool usediag; /* indicates preconditioner should include diagonal scaling*/
15: } PC_Eisenstat;
17: static PetscErrorCode PCMult_Eisenstat(Mat mat, Vec b, Vec x)
18: {
19: PC pc;
20: PC_Eisenstat *eis;
22: PetscFunctionBegin;
23: PetscCall(MatShellGetContext(mat, &pc));
24: eis = (PC_Eisenstat *)pc->data;
25: PetscCall(MatSOR(eis->A, b, eis->omega, SOR_EISENSTAT, 0.0, 1, 1, x));
26: PetscCall(MatFactorGetError(eis->A, (MatFactorError *)&pc->failedreason));
27: PetscFunctionReturn(PETSC_SUCCESS);
28: }
30: static PetscErrorCode PCNorm_Eisenstat(Mat mat, NormType type, PetscReal *nrm)
31: {
32: PC pc;
33: PC_Eisenstat *eis;
35: PetscFunctionBegin;
36: PetscCall(MatShellGetContext(mat, &pc));
37: eis = (PC_Eisenstat *)pc->data;
38: PetscCall(MatNorm(eis->A, type, nrm));
39: PetscFunctionReturn(PETSC_SUCCESS);
40: }
42: static PetscErrorCode PCApply_Eisenstat(PC pc, Vec x, Vec y)
43: {
44: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
45: PetscBool hasop;
47: PetscFunctionBegin;
48: if (eis->usediag) {
49: PetscCall(MatHasOperation(pc->pmat, MATOP_MULT_DIAGONAL_BLOCK, &hasop));
50: if (hasop) {
51: PetscCall(MatMultDiagonalBlock(pc->pmat, x, y));
52: } else {
53: PetscCall(VecPointwiseMult(y, x, eis->diag));
54: }
55: } else PetscCall(VecCopy(x, y));
56: PetscFunctionReturn(PETSC_SUCCESS);
57: }
59: static PetscErrorCode PCApplyTranspose_Eisenstat(PC pc, Vec x, Vec y)
60: {
61: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
62: PetscBool hasop, set, sym;
64: PetscFunctionBegin;
65: PetscCall(MatIsSymmetricKnown(eis->A, &set, &sym));
66: PetscCheck(set && sym, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Can only apply transpose of Eisenstat if matrix is symmetric");
67: if (eis->usediag) {
68: PetscCall(MatHasOperation(pc->pmat, MATOP_MULT_DIAGONAL_BLOCK, &hasop));
69: if (hasop) {
70: PetscCall(MatMultDiagonalBlock(pc->pmat, x, y));
71: } else {
72: PetscCall(VecPointwiseMult(y, x, eis->diag));
73: }
74: } else PetscCall(VecCopy(x, y));
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: static PetscErrorCode PCPreSolve_Eisenstat(PC pc, KSP ksp, Vec b, Vec x)
79: {
80: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
81: PetscBool nonzero;
83: PetscFunctionBegin;
84: if (pc->presolvedone < 2) {
85: PetscCheck(pc->mat == pc->pmat, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot have different mat and pmat");
86: /* swap shell matrix and true matrix */
87: eis->A = pc->mat;
88: pc->mat = eis->shell;
89: }
91: if (!eis->b[pc->presolvedone - 1]) { PetscCall(VecDuplicate(b, &eis->b[pc->presolvedone - 1])); }
93: /* if nonzero initial guess, modify x */
94: PetscCall(KSPGetInitialGuessNonzero(ksp, &nonzero));
95: if (nonzero) {
96: PetscCall(VecCopy(x, eis->b[pc->presolvedone - 1]));
97: PetscCall(MatSOR(eis->A, eis->b[pc->presolvedone - 1], eis->omega, SOR_APPLY_UPPER, 0.0, 1, 1, x));
98: PetscCall(MatFactorGetError(eis->A, (MatFactorError *)&pc->failedreason));
99: }
101: /* save true b, other option is to swap pointers */
102: PetscCall(VecCopy(b, eis->b[pc->presolvedone - 1]));
104: /* modify b by (L + D/omega)^{-1} */
105: 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));
106: PetscCall(MatFactorGetError(eis->A, (MatFactorError *)&pc->failedreason));
107: PetscFunctionReturn(PETSC_SUCCESS);
108: }
110: static PetscErrorCode PCPostSolve_Eisenstat(PC pc, KSP ksp, Vec b, Vec x)
111: {
112: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
114: PetscFunctionBegin;
115: /* get back true b */
116: PetscCall(VecCopy(eis->b[pc->presolvedone], b));
118: /* modify x by (U + D/omega)^{-1} */
119: PetscCall(VecCopy(x, eis->b[pc->presolvedone]));
120: PetscCall(MatSOR(eis->A, eis->b[pc->presolvedone], eis->omega, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_BACKWARD_SWEEP), 0.0, 1, 1, x));
121: PetscCall(MatFactorGetError(eis->A, (MatFactorError *)&pc->failedreason));
122: if (!pc->presolvedone) pc->mat = eis->A;
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: static PetscErrorCode PCReset_Eisenstat(PC pc)
127: {
128: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
130: PetscFunctionBegin;
131: PetscCall(VecDestroy(&eis->b[0]));
132: PetscCall(VecDestroy(&eis->b[1]));
133: PetscCall(MatDestroy(&eis->shell));
134: PetscCall(VecDestroy(&eis->diag));
135: PetscFunctionReturn(PETSC_SUCCESS);
136: }
138: static PetscErrorCode PCDestroy_Eisenstat(PC pc)
139: {
140: PetscFunctionBegin;
141: PetscCall(PCReset_Eisenstat(pc));
142: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetOmega_C", NULL));
143: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetNoDiagonalScaling_C", NULL));
144: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetOmega_C", NULL));
145: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetNoDiagonalScaling_C", NULL));
146: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", NULL));
147: PetscCall(PetscFree(pc->data));
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: static PetscErrorCode PCSetFromOptions_Eisenstat(PC pc, PetscOptionItems *PetscOptionsObject)
152: {
153: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
154: PetscBool set, flg;
155: PetscReal omega;
157: PetscFunctionBegin;
158: PetscOptionsHeadBegin(PetscOptionsObject, "Eisenstat SSOR options");
159: PetscCall(PetscOptionsReal("-pc_eisenstat_omega", "Relaxation factor 0 < omega < 2", "PCEisenstatSetOmega", eis->omega, &omega, &flg));
160: if (flg) PetscCall(PCEisenstatSetOmega(pc, omega));
161: PetscCall(PetscOptionsBool("-pc_eisenstat_no_diagonal_scaling", "Do not use standard diagonal scaling", "PCEisenstatSetNoDiagonalScaling", eis->usediag ? PETSC_FALSE : PETSC_TRUE, &flg, &set));
162: if (set) PetscCall(PCEisenstatSetNoDiagonalScaling(pc, flg));
163: PetscOptionsHeadEnd();
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }
167: static PetscErrorCode PCView_Eisenstat(PC pc, PetscViewer viewer)
168: {
169: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
170: PetscBool iascii;
172: PetscFunctionBegin;
173: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
174: if (iascii) {
175: PetscCall(PetscViewerASCIIPrintf(viewer, " omega = %g\n", (double)eis->omega));
176: if (eis->usediag) {
177: PetscCall(PetscViewerASCIIPrintf(viewer, " Using diagonal scaling (default)\n"));
178: } else {
179: PetscCall(PetscViewerASCIIPrintf(viewer, " Not using diagonal scaling\n"));
180: }
181: }
182: PetscFunctionReturn(PETSC_SUCCESS);
183: }
185: static PetscErrorCode PCSetUp_Eisenstat(PC pc)
186: {
187: PetscInt M, N, m, n;
188: PetscBool set, sym;
189: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
191: PetscFunctionBegin;
192: if (!pc->setupcalled) {
193: PetscCall(MatGetSize(pc->mat, &M, &N));
194: PetscCall(MatGetLocalSize(pc->mat, &m, &n));
195: PetscCall(MatIsSymmetricKnown(pc->mat, &set, &sym));
196: PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &eis->shell));
197: PetscCall(MatSetSizes(eis->shell, m, n, M, N));
198: PetscCall(MatSetType(eis->shell, MATSHELL));
199: PetscCall(MatSetUp(eis->shell));
200: PetscCall(MatShellSetContext(eis->shell, pc));
201: PetscCall(MatShellSetOperation(eis->shell, MATOP_MULT, (void (*)(void))PCMult_Eisenstat));
202: if (set && sym) PetscCall(MatShellSetOperation(eis->shell, MATOP_MULT_TRANSPOSE, (void (*)(void))PCMult_Eisenstat));
203: PetscCall(MatShellSetOperation(eis->shell, MATOP_NORM, (void (*)(void))PCNorm_Eisenstat));
204: }
205: if (!eis->usediag) PetscFunctionReturn(PETSC_SUCCESS);
206: if (!pc->setupcalled) { PetscCall(MatCreateVecs(pc->pmat, &eis->diag, NULL)); }
207: PetscCall(MatGetDiagonal(pc->pmat, eis->diag));
208: PetscFunctionReturn(PETSC_SUCCESS);
209: }
211: /* --------------------------------------------------------------------*/
213: static PetscErrorCode PCEisenstatSetOmega_Eisenstat(PC pc, PetscReal omega)
214: {
215: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
217: PetscFunctionBegin;
218: PetscCheck(omega > 0.0 && omega < 2.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Relaxation out of range");
219: eis->omega = omega;
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
223: static PetscErrorCode PCEisenstatSetNoDiagonalScaling_Eisenstat(PC pc, PetscBool flg)
224: {
225: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
227: PetscFunctionBegin;
228: eis->usediag = flg;
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
232: static PetscErrorCode PCEisenstatGetOmega_Eisenstat(PC pc, PetscReal *omega)
233: {
234: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
236: PetscFunctionBegin;
237: *omega = eis->omega;
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
241: static PetscErrorCode PCEisenstatGetNoDiagonalScaling_Eisenstat(PC pc, PetscBool *flg)
242: {
243: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
245: PetscFunctionBegin;
246: *flg = eis->usediag;
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
250: /*@
251: PCEisenstatSetOmega - Sets the SSOR relaxation coefficient, omega,
252: to use with Eisenstat's trick (where omega = 1.0 by default)
254: Logically Collective
256: Input Parameters:
257: + pc - the preconditioner context
258: - omega - relaxation coefficient (0 < omega < 2)
260: Options Database Key:
261: . -pc_eisenstat_omega <omega> - Sets omega
263: Level: intermediate
265: Notes:
266: The Eisenstat trick implementation of SSOR requires about 50% of the
267: usual amount of floating point operations used for SSOR + Krylov method;
268: however, the preconditioned problem must be solved with both left
269: and right preconditioning.
271: To use SSOR without the Eisenstat trick, employ the `PCSOR` preconditioner,
272: which can be chosen with the database options `-pc_type sor -pc_sor_symmetric`
274: .seealso: `PCSORSetOmega()`, `PCEISENSTAT`
275: @*/
276: PetscErrorCode PCEisenstatSetOmega(PC pc, PetscReal omega)
277: {
278: PetscFunctionBegin;
281: PetscTryMethod(pc, "PCEisenstatSetOmega_C", (PC, PetscReal), (pc, omega));
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }
285: /*@
286: PCEisenstatSetNoDiagonalScaling - Causes the Eisenstat preconditioner, `PCEISENSTAT`
287: not to do additional diagonal preconditioning. For matrices with a constant
288: along the diagonal, this may save a small amount of work.
290: Logically Collective
292: Input Parameters:
293: + pc - the preconditioner context
294: - flg - `PETSC_TRUE` turns off diagonal scaling inside the algorithm
296: Options Database Key:
297: . -pc_eisenstat_no_diagonal_scaling - Activates `PCEisenstatSetNoDiagonalScaling()`
299: Level: intermediate
301: Note:
302: If you use the `KSPSetDiagonalScaling()` or -ksp_diagonal_scale option then you will
303: likely want to use this routine since it will save you some unneeded flops.
305: .seealso: `PCEisenstatSetOmega()`, `PCEISENSTAT`
306: @*/
307: PetscErrorCode PCEisenstatSetNoDiagonalScaling(PC pc, PetscBool flg)
308: {
309: PetscFunctionBegin;
311: PetscTryMethod(pc, "PCEisenstatSetNoDiagonalScaling_C", (PC, PetscBool), (pc, flg));
312: PetscFunctionReturn(PETSC_SUCCESS);
313: }
315: /*@
316: PCEisenstatGetOmega - Gets the SSOR relaxation coefficient, omega,
317: to use with Eisenstat's trick (where omega = 1.0 by default).
319: Logically Collective
321: Input Parameter:
322: . pc - the preconditioner context
324: Output Parameter:
325: . omega - relaxation coefficient (0 < omega < 2)
327: Options Database Key:
328: . -pc_eisenstat_omega <omega> - Sets omega
330: Notes:
331: The Eisenstat trick implementation of SSOR requires about 50% of the
332: usual amount of floating point operations used for SSOR + Krylov method;
333: however, the preconditioned problem must be solved with both left
334: and right preconditioning.
336: To use SSOR without the Eisenstat trick, employ the PCSOR preconditioner,
337: which can be chosen with the database options `-pc_type sor -pc_sor_symmetric`
339: Level: intermediate
341: .seealso: `PCEISENSTAT`, `PCSORGetOmega()`, `PCEisenstatSetOmega()`
342: @*/
343: PetscErrorCode PCEisenstatGetOmega(PC pc, PetscReal *omega)
344: {
345: PetscFunctionBegin;
347: PetscUseMethod(pc, "PCEisenstatGetOmega_C", (PC, PetscReal *), (pc, omega));
348: PetscFunctionReturn(PETSC_SUCCESS);
349: }
351: /*@
352: PCEisenstatGetNoDiagonalScaling - Tells if the Eisenstat preconditioner
353: not to do additional diagonal preconditioning. For matrices with a constant
354: along the diagonal, this may save a small amount of work.
356: Logically Collective
358: Input Parameter:
359: . pc - the preconditioner context
361: Output Parameter:
362: . flg - `PETSC_TRUE` means there is no diagonal scaling applied
364: Options Database Key:
365: . -pc_eisenstat_no_diagonal_scaling - Activates `PCEisenstatSetNoDiagonalScaling()`
367: Level: intermediate
369: Note:
370: If you use the KSPSetDiagonalScaling() or -ksp_diagonal_scale option then you will
371: likely want to use this routine since it will save you some unneeded flops.
373: .seealso: , `PCEISENSTAT`, `PCEisenstatGetOmega()`
374: @*/
375: PetscErrorCode PCEisenstatGetNoDiagonalScaling(PC pc, PetscBool *flg)
376: {
377: PetscFunctionBegin;
379: PetscUseMethod(pc, "PCEisenstatGetNoDiagonalScaling_C", (PC, PetscBool *), (pc, flg));
380: PetscFunctionReturn(PETSC_SUCCESS);
381: }
383: static PetscErrorCode PCPreSolveChangeRHS_Eisenstat(PC pc, PetscBool *change)
384: {
385: PetscFunctionBegin;
386: *change = PETSC_TRUE;
387: PetscFunctionReturn(PETSC_SUCCESS);
388: }
390: /*MC
391: PCEISENSTAT - An implementation of SSOR (symmetric successive over relaxation, symmetric Gauss-Seidel)
392: preconditioning that incorporates Eisenstat's trick to reduce the amount of computation needed.
394: Options Database Keys:
395: + -pc_eisenstat_omega <omega> - Sets omega
396: - -pc_eisenstat_no_diagonal_scaling - Activates `PCEisenstatSetNoDiagonalScaling()`
398: Level: beginner
400: Notes:
401: Only implemented for the `MATAIJ` matrix format.
403: Not a true parallel SOR, in parallel this implementation corresponds to block Jacobi with SOR on each block.
405: Developer Note:
406: Since this algorithm runs the Krylov method on a transformed linear system the implementation provides `PCPreSolve()` and `PCPostSolve()`
407: routines that `KSP` uses to set up the transformed linear system.
409: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCEisenstatGetOmega()`,
410: `PCEisenstatSetNoDiagonalScaling()`, `PCEisenstatSetOmega()`, `PCSOR`
411: M*/
413: PETSC_EXTERN PetscErrorCode PCCreate_Eisenstat(PC pc)
414: {
415: PC_Eisenstat *eis;
417: PetscFunctionBegin;
418: PetscCall(PetscNew(&eis));
420: pc->ops->apply = PCApply_Eisenstat;
421: pc->ops->applytranspose = PCApplyTranspose_Eisenstat;
422: pc->ops->presolve = PCPreSolve_Eisenstat;
423: pc->ops->postsolve = PCPostSolve_Eisenstat;
424: pc->ops->applyrichardson = NULL;
425: pc->ops->setfromoptions = PCSetFromOptions_Eisenstat;
426: pc->ops->destroy = PCDestroy_Eisenstat;
427: pc->ops->reset = PCReset_Eisenstat;
428: pc->ops->view = PCView_Eisenstat;
429: pc->ops->setup = PCSetUp_Eisenstat;
431: pc->data = eis;
432: eis->omega = 1.0;
433: eis->b[0] = NULL;
434: eis->b[1] = NULL;
435: eis->diag = NULL;
436: eis->usediag = PETSC_TRUE;
438: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetOmega_C", PCEisenstatSetOmega_Eisenstat));
439: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetNoDiagonalScaling_C", PCEisenstatSetNoDiagonalScaling_Eisenstat));
440: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetOmega_C", PCEisenstatGetOmega_Eisenstat));
441: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetNoDiagonalScaling_C", PCEisenstatGetNoDiagonalScaling_Eisenstat));
442: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", PCPreSolveChangeRHS_Eisenstat));
443: PetscFunctionReturn(PETSC_SUCCESS);
444: }