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 isascii;
171: PetscFunctionBegin;
172: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
173: if (isascii) {
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, (PetscErrorCodeFn *)PCMult_Eisenstat));
201: if (set && sym) PetscCall(MatShellSetOperation(eis->shell, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)PCMult_Eisenstat));
202: PetscCall(MatShellSetOperation(eis->shell, MATOP_NORM, (PetscErrorCodeFn *)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: static PetscErrorCode PCEisenstatSetOmega_Eisenstat(PC pc, PetscReal omega)
211: {
212: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
214: PetscFunctionBegin;
215: PetscCheck(omega > 0.0 && omega < 2.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Relaxation out of range");
216: eis->omega = omega;
217: PetscFunctionReturn(PETSC_SUCCESS);
218: }
220: static PetscErrorCode PCEisenstatSetNoDiagonalScaling_Eisenstat(PC pc, PetscBool flg)
221: {
222: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
224: PetscFunctionBegin;
225: eis->usediag = flg;
226: PetscFunctionReturn(PETSC_SUCCESS);
227: }
229: static PetscErrorCode PCEisenstatGetOmega_Eisenstat(PC pc, PetscReal *omega)
230: {
231: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
233: PetscFunctionBegin;
234: *omega = eis->omega;
235: PetscFunctionReturn(PETSC_SUCCESS);
236: }
238: static PetscErrorCode PCEisenstatGetNoDiagonalScaling_Eisenstat(PC pc, PetscBool *flg)
239: {
240: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
242: PetscFunctionBegin;
243: *flg = eis->usediag;
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }
247: /*@
248: PCEisenstatSetOmega - Sets the SSOR relaxation coefficient, omega,
249: to use with Eisenstat's trick (where omega = 1.0 by default)
251: Logically Collective
253: Input Parameters:
254: + pc - the preconditioner context
255: - omega - relaxation coefficient (0 < omega < 2)
257: Options Database Key:
258: . -pc_eisenstat_omega <omega> - Sets omega
260: Level: intermediate
262: Notes:
263: The Eisenstat trick implementation of SSOR requires about 50% of the
264: usual amount of floating point operations used for SSOR + Krylov method;
265: however, the preconditioned problem must be solved with both left
266: and right preconditioning.
268: To use SSOR without the Eisenstat trick, employ the `PCSOR` preconditioner,
269: which can be chosen with the database options `-pc_type sor -pc_sor_symmetric`
271: .seealso: [](ch_ksp), `PCSORSetOmega()`, `PCEISENSTAT`
272: @*/
273: PetscErrorCode PCEisenstatSetOmega(PC pc, PetscReal omega)
274: {
275: PetscFunctionBegin;
278: PetscTryMethod(pc, "PCEisenstatSetOmega_C", (PC, PetscReal), (pc, omega));
279: PetscFunctionReturn(PETSC_SUCCESS);
280: }
282: /*@
283: PCEisenstatSetNoDiagonalScaling - Causes the Eisenstat preconditioner, `PCEISENSTAT`
284: not to do additional diagonal preconditioning. For matrices with a constant
285: along the diagonal, this may save a small amount of work.
287: Logically Collective
289: Input Parameters:
290: + pc - the preconditioner context
291: - flg - `PETSC_TRUE` turns off diagonal scaling inside the algorithm
293: Options Database Key:
294: . -pc_eisenstat_no_diagonal_scaling - Activates `PCEisenstatSetNoDiagonalScaling()`
296: Level: intermediate
298: Note:
299: If you use the `KSPSetDiagonalScaling()` or -ksp_diagonal_scale option then you will
300: likely want to use this routine since it will save you some unneeded flops.
302: .seealso: [](ch_ksp), `PCEisenstatSetOmega()`, `PCEISENSTAT`
303: @*/
304: PetscErrorCode PCEisenstatSetNoDiagonalScaling(PC pc, PetscBool flg)
305: {
306: PetscFunctionBegin;
308: PetscTryMethod(pc, "PCEisenstatSetNoDiagonalScaling_C", (PC, PetscBool), (pc, flg));
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
312: /*@
313: PCEisenstatGetOmega - Gets the SSOR relaxation coefficient, omega,
314: to use with Eisenstat's trick (where omega = 1.0 by default).
316: Logically Collective
318: Input Parameter:
319: . pc - the preconditioner context
321: Output Parameter:
322: . omega - relaxation coefficient (0 < omega < 2)
324: Options Database Key:
325: . -pc_eisenstat_omega <omega> - Sets omega
327: Notes:
328: The Eisenstat trick implementation of SSOR requires about 50% of the
329: usual amount of floating point operations used for SSOR + Krylov method;
330: however, the preconditioned problem must be solved with both left
331: and right preconditioning.
333: To use SSOR without the Eisenstat trick, employ the PCSOR preconditioner,
334: which can be chosen with the database options `-pc_type sor -pc_sor_symmetric`
336: Level: intermediate
338: .seealso: [](ch_ksp), `PCEISENSTAT`, `PCSORGetOmega()`, `PCEisenstatSetOmega()`
339: @*/
340: PetscErrorCode PCEisenstatGetOmega(PC pc, PetscReal *omega)
341: {
342: PetscFunctionBegin;
344: PetscUseMethod(pc, "PCEisenstatGetOmega_C", (PC, PetscReal *), (pc, omega));
345: PetscFunctionReturn(PETSC_SUCCESS);
346: }
348: /*@
349: PCEisenstatGetNoDiagonalScaling - Tells if the Eisenstat preconditioner
350: not to do additional diagonal preconditioning. For matrices with a constant
351: along the diagonal, this may save a small amount of work.
353: Logically Collective
355: Input Parameter:
356: . pc - the preconditioner context
358: Output Parameter:
359: . flg - `PETSC_TRUE` means there is no diagonal scaling applied
361: Options Database Key:
362: . -pc_eisenstat_no_diagonal_scaling - Activates `PCEisenstatSetNoDiagonalScaling()`
364: Level: intermediate
366: Note:
367: If you use the KSPSetDiagonalScaling() or -ksp_diagonal_scale option then you will
368: likely want to use this routine since it will save you some unneeded flops.
370: .seealso: , `PCEISENSTAT`, `PCEisenstatGetOmega()`
371: @*/
372: PetscErrorCode PCEisenstatGetNoDiagonalScaling(PC pc, PetscBool *flg)
373: {
374: PetscFunctionBegin;
376: PetscUseMethod(pc, "PCEisenstatGetNoDiagonalScaling_C", (PC, PetscBool *), (pc, flg));
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
380: static PetscErrorCode PCPreSolveChangeRHS_Eisenstat(PC pc, PetscBool *change)
381: {
382: PetscFunctionBegin;
383: *change = PETSC_TRUE;
384: PetscFunctionReturn(PETSC_SUCCESS);
385: }
387: /*MC
388: PCEISENSTAT - An implementation of SSOR (symmetric successive over relaxation, symmetric Gauss-Seidel)
389: preconditioning that incorporates Eisenstat's trick to reduce the amount of computation needed.
391: Options Database Keys:
392: + -pc_eisenstat_omega <omega> - Sets omega
393: - -pc_eisenstat_no_diagonal_scaling - Activates `PCEisenstatSetNoDiagonalScaling()`
395: Level: beginner
397: Notes:
398: Only implemented for the `MATAIJ` matrix format.
400: Not a true parallel SOR, in parallel this implementation corresponds to block Jacobi with SOR on each block.
402: Developer Note:
403: Since this algorithm runs the Krylov method on a transformed linear system the implementation provides `PCPreSolve()` and `PCPostSolve()`
404: routines that `KSP` uses to set up the transformed linear system.
406: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCEisenstatGetOmega()`,
407: `PCEisenstatSetNoDiagonalScaling()`, `PCEisenstatSetOmega()`, `PCSOR`
408: M*/
410: PETSC_EXTERN PetscErrorCode PCCreate_Eisenstat(PC pc)
411: {
412: PC_Eisenstat *eis;
414: PetscFunctionBegin;
415: PetscCall(PetscNew(&eis));
417: pc->ops->apply = PCApply_Eisenstat;
418: pc->ops->applytranspose = PCApplyTranspose_Eisenstat;
419: pc->ops->presolve = PCPreSolve_Eisenstat;
420: pc->ops->postsolve = PCPostSolve_Eisenstat;
421: pc->ops->applyrichardson = NULL;
422: pc->ops->setfromoptions = PCSetFromOptions_Eisenstat;
423: pc->ops->destroy = PCDestroy_Eisenstat;
424: pc->ops->reset = PCReset_Eisenstat;
425: pc->ops->view = PCView_Eisenstat;
426: pc->ops->setup = PCSetUp_Eisenstat;
428: pc->data = eis;
429: eis->omega = 1.0;
430: eis->b[0] = NULL;
431: eis->b[1] = NULL;
432: eis->diag = NULL;
433: eis->usediag = PETSC_TRUE;
435: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetOmega_C", PCEisenstatSetOmega_Eisenstat));
436: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetNoDiagonalScaling_C", PCEisenstatSetNoDiagonalScaling_Eisenstat));
437: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetOmega_C", PCEisenstatGetOmega_Eisenstat));
438: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetNoDiagonalScaling_C", PCEisenstatGetNoDiagonalScaling_Eisenstat));
439: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", PCPreSolveChangeRHS_Eisenstat));
440: PetscFunctionReturn(PETSC_SUCCESS);
441: }