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) PetscCall(MatMultDiagonalBlock(pc->pmat, x, y));
50: else PetscCall(VecPointwiseMult(y, x, eis->diag));
51: } else PetscCall(VecCopy(x, y));
52: PetscFunctionReturn(PETSC_SUCCESS);
53: }
55: static PetscErrorCode PCApplyTranspose_Eisenstat(PC pc, Vec x, Vec y)
56: {
57: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
58: PetscBool hasop, set, sym;
60: PetscFunctionBegin;
61: PetscCall(MatIsSymmetricKnown(eis->A, &set, &sym));
62: PetscCheck(set && sym, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Can only apply transpose of Eisenstat if matrix is symmetric");
63: if (eis->usediag) {
64: PetscCall(MatHasOperation(pc->pmat, MATOP_MULT_DIAGONAL_BLOCK, &hasop));
65: if (hasop) PetscCall(MatMultDiagonalBlock(pc->pmat, x, y));
66: else PetscCall(VecPointwiseMult(y, x, eis->diag));
67: } else PetscCall(VecCopy(x, y));
68: PetscFunctionReturn(PETSC_SUCCESS);
69: }
71: static PetscErrorCode PCPreSolve_Eisenstat(PC pc, KSP ksp, Vec b, Vec x)
72: {
73: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
74: PetscBool nonzero;
76: PetscFunctionBegin;
77: if (pc->presolvedone < 2) {
78: PetscCheck(pc->mat == pc->pmat, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot have different mat and pmat");
79: /* swap shell matrix and true matrix */
80: eis->A = pc->mat;
81: pc->mat = eis->shell;
82: }
84: if (!eis->b[pc->presolvedone - 1]) PetscCall(VecDuplicate(b, &eis->b[pc->presolvedone - 1]));
86: /* if nonzero initial guess, modify x */
87: PetscCall(KSPGetInitialGuessNonzero(ksp, &nonzero));
88: if (nonzero) {
89: PetscCall(VecCopy(x, eis->b[pc->presolvedone - 1]));
90: PetscCall(MatSOR(eis->A, eis->b[pc->presolvedone - 1], eis->omega, SOR_APPLY_UPPER, 0.0, 1, 1, x));
91: PetscCall(MatFactorGetError(eis->A, (MatFactorError *)&pc->failedreason));
92: }
94: /* save true b, other option is to swap pointers */
95: PetscCall(VecCopy(b, eis->b[pc->presolvedone - 1]));
97: /* modify b by (L + D/omega)^{-1} */
98: 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));
99: PetscCall(MatFactorGetError(eis->A, (MatFactorError *)&pc->failedreason));
100: PetscFunctionReturn(PETSC_SUCCESS);
101: }
103: static PetscErrorCode PCPostSolve_Eisenstat(PC pc, KSP ksp, Vec b, Vec x)
104: {
105: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
107: PetscFunctionBegin;
108: /* get back true b */
109: PetscCall(VecCopy(eis->b[pc->presolvedone], b));
111: /* modify x by (U + D/omega)^{-1} */
112: PetscCall(VecCopy(x, eis->b[pc->presolvedone]));
113: PetscCall(MatSOR(eis->A, eis->b[pc->presolvedone], eis->omega, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_BACKWARD_SWEEP), 0.0, 1, 1, x));
114: PetscCall(MatFactorGetError(eis->A, (MatFactorError *)&pc->failedreason));
115: if (!pc->presolvedone) pc->mat = eis->A;
116: PetscFunctionReturn(PETSC_SUCCESS);
117: }
119: static PetscErrorCode PCReset_Eisenstat(PC pc)
120: {
121: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
123: PetscFunctionBegin;
124: PetscCall(VecDestroy(&eis->b[0]));
125: PetscCall(VecDestroy(&eis->b[1]));
126: PetscCall(MatDestroy(&eis->shell));
127: PetscCall(VecDestroy(&eis->diag));
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: static PetscErrorCode PCDestroy_Eisenstat(PC pc)
132: {
133: PetscFunctionBegin;
134: PetscCall(PCReset_Eisenstat(pc));
135: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetOmega_C", NULL));
136: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetNoDiagonalScaling_C", NULL));
137: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetOmega_C", NULL));
138: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetNoDiagonalScaling_C", NULL));
139: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", NULL));
140: PetscCall(PetscFree(pc->data));
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: static PetscErrorCode PCSetFromOptions_Eisenstat(PC pc, PetscOptionItems PetscOptionsObject)
145: {
146: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
147: PetscBool set, flg;
148: PetscReal omega;
150: PetscFunctionBegin;
151: PetscOptionsHeadBegin(PetscOptionsObject, "Eisenstat SSOR options");
152: PetscCall(PetscOptionsReal("-pc_eisenstat_omega", "Relaxation factor 0 < omega < 2", "PCEisenstatSetOmega", eis->omega, &omega, &flg));
153: if (flg) PetscCall(PCEisenstatSetOmega(pc, omega));
154: PetscCall(PetscOptionsBool("-pc_eisenstat_no_diagonal_scaling", "Do not use standard diagonal scaling", "PCEisenstatSetNoDiagonalScaling", eis->usediag ? PETSC_FALSE : PETSC_TRUE, &flg, &set));
155: if (set) PetscCall(PCEisenstatSetNoDiagonalScaling(pc, flg));
156: PetscOptionsHeadEnd();
157: PetscFunctionReturn(PETSC_SUCCESS);
158: }
160: static PetscErrorCode PCView_Eisenstat(PC pc, PetscViewer viewer)
161: {
162: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
163: PetscBool isascii;
165: PetscFunctionBegin;
166: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
167: if (isascii) {
168: PetscCall(PetscViewerASCIIPrintf(viewer, " omega = %g\n", (double)eis->omega));
169: if (eis->usediag) {
170: PetscCall(PetscViewerASCIIPrintf(viewer, " Using diagonal scaling (default)\n"));
171: } else {
172: PetscCall(PetscViewerASCIIPrintf(viewer, " Not using diagonal scaling\n"));
173: }
174: }
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: static PetscErrorCode PCSetUp_Eisenstat(PC pc)
179: {
180: PetscInt M, N, m, n;
181: PetscBool set, sym;
182: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
184: PetscFunctionBegin;
185: if (!pc->setupcalled) {
186: PetscCall(MatGetSize(pc->mat, &M, &N));
187: PetscCall(MatGetLocalSize(pc->mat, &m, &n));
188: PetscCall(MatIsSymmetricKnown(pc->mat, &set, &sym));
189: PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &eis->shell));
190: PetscCall(MatSetSizes(eis->shell, m, n, M, N));
191: PetscCall(MatSetType(eis->shell, MATSHELL));
192: PetscCall(MatSetUp(eis->shell));
193: PetscCall(MatShellSetContext(eis->shell, pc));
194: PetscCall(MatShellSetOperation(eis->shell, MATOP_MULT, (PetscErrorCodeFn *)PCMult_Eisenstat));
195: if (set && sym) PetscCall(MatShellSetOperation(eis->shell, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)PCMult_Eisenstat));
196: PetscCall(MatShellSetOperation(eis->shell, MATOP_NORM, (PetscErrorCodeFn *)PCNorm_Eisenstat));
197: }
198: if (!eis->usediag) PetscFunctionReturn(PETSC_SUCCESS);
199: if (!pc->setupcalled) PetscCall(MatCreateVecs(pc->pmat, &eis->diag, NULL));
200: PetscCall(MatGetDiagonal(pc->pmat, eis->diag));
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: static PetscErrorCode PCEisenstatSetOmega_Eisenstat(PC pc, PetscReal omega)
205: {
206: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
208: PetscFunctionBegin;
209: PetscCheck(omega > 0.0 && omega < 2.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Relaxation out of range");
210: eis->omega = omega;
211: PetscFunctionReturn(PETSC_SUCCESS);
212: }
214: static PetscErrorCode PCEisenstatSetNoDiagonalScaling_Eisenstat(PC pc, PetscBool flg)
215: {
216: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
218: PetscFunctionBegin;
219: eis->usediag = flg;
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
223: static PetscErrorCode PCEisenstatGetOmega_Eisenstat(PC pc, PetscReal *omega)
224: {
225: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
227: PetscFunctionBegin;
228: *omega = eis->omega;
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
232: static PetscErrorCode PCEisenstatGetNoDiagonalScaling_Eisenstat(PC pc, PetscBool *flg)
233: {
234: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
236: PetscFunctionBegin;
237: *flg = eis->usediag;
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
241: /*@
242: PCEisenstatSetOmega - Sets the SSOR relaxation coefficient, omega,
243: to use with Eisenstat's trick (where omega = 1.0 by default)
245: Logically Collective
247: Input Parameters:
248: + pc - the preconditioner context
249: - omega - relaxation coefficient (0 < omega < 2)
251: Options Database Key:
252: . -pc_eisenstat_omega omega - Sets omega
254: Level: intermediate
256: Notes:
257: The Eisenstat trick implementation of SSOR requires about 50% of the
258: usual amount of floating point operations used for SSOR + Krylov method;
259: however, the preconditioned problem must be solved with both left
260: and right preconditioning.
262: To use SSOR without the Eisenstat trick, employ the `PCSOR` preconditioner,
263: which can be chosen with the database options `-pc_type sor -pc_sor_symmetric`
265: .seealso: [](ch_ksp), `PCSORSetOmega()`, `PCEISENSTAT`
266: @*/
267: PetscErrorCode PCEisenstatSetOmega(PC pc, PetscReal omega)
268: {
269: PetscFunctionBegin;
272: PetscTryMethod(pc, "PCEisenstatSetOmega_C", (PC, PetscReal), (pc, omega));
273: PetscFunctionReturn(PETSC_SUCCESS);
274: }
276: /*@
277: PCEisenstatSetNoDiagonalScaling - Causes the Eisenstat preconditioner, `PCEISENSTAT`
278: not to do additional diagonal preconditioning. For matrices with a constant
279: along the diagonal, this may save a small amount of work.
281: Logically Collective
283: Input Parameters:
284: + pc - the preconditioner context
285: - flg - `PETSC_TRUE` turns off diagonal scaling inside the algorithm
287: Options Database Key:
288: . -pc_eisenstat_no_diagonal_scaling - Activates `PCEisenstatSetNoDiagonalScaling()`
290: Level: intermediate
292: .seealso: [](ch_ksp), `PCEisenstatSetOmega()`, `PCEISENSTAT`
293: @*/
294: PetscErrorCode PCEisenstatSetNoDiagonalScaling(PC pc, PetscBool flg)
295: {
296: PetscFunctionBegin;
298: PetscTryMethod(pc, "PCEisenstatSetNoDiagonalScaling_C", (PC, PetscBool), (pc, flg));
299: PetscFunctionReturn(PETSC_SUCCESS);
300: }
302: /*@
303: PCEisenstatGetOmega - Gets the SSOR relaxation coefficient, omega,
304: to use with Eisenstat's trick (where omega = 1.0 by default).
306: Logically Collective
308: Input Parameter:
309: . pc - the preconditioner context
311: Output Parameter:
312: . omega - relaxation coefficient (0 < omega < 2)
314: Options Database Key:
315: . -pc_eisenstat_omega omega - Sets omega
317: Notes:
318: The Eisenstat trick implementation of SSOR requires about 50% of the
319: usual amount of floating point operations used for SSOR + Krylov method;
320: however, the preconditioned problem must be solved with both left
321: and right preconditioning.
323: To use SSOR without the Eisenstat trick, employ the PCSOR preconditioner,
324: which can be chosen with the database options `-pc_type sor -pc_sor_symmetric`
326: Level: intermediate
328: .seealso: [](ch_ksp), `PCEISENSTAT`, `PCSORGetOmega()`, `PCEisenstatSetOmega()`
329: @*/
330: PetscErrorCode PCEisenstatGetOmega(PC pc, PetscReal *omega)
331: {
332: PetscFunctionBegin;
334: PetscUseMethod(pc, "PCEisenstatGetOmega_C", (PC, PetscReal *), (pc, omega));
335: PetscFunctionReturn(PETSC_SUCCESS);
336: }
338: /*@
339: PCEisenstatGetNoDiagonalScaling - Tells if the Eisenstat preconditioner
340: not to do additional diagonal preconditioning. For matrices with a constant
341: along the diagonal, this may save a small amount of work.
343: Logically Collective
345: Input Parameter:
346: . pc - the preconditioner context
348: Output Parameter:
349: . flg - `PETSC_TRUE` means there is no diagonal scaling applied
351: Options Database Key:
352: . -pc_eisenstat_no_diagonal_scaling - Activates `PCEisenstatSetNoDiagonalScaling()`
354: Level: intermediate
356: .seealso: `PCEISENSTAT`, `PCEisenstatGetOmega()`
357: @*/
358: PetscErrorCode PCEisenstatGetNoDiagonalScaling(PC pc, PetscBool *flg)
359: {
360: PetscFunctionBegin;
362: PetscUseMethod(pc, "PCEisenstatGetNoDiagonalScaling_C", (PC, PetscBool *), (pc, flg));
363: PetscFunctionReturn(PETSC_SUCCESS);
364: }
366: static PetscErrorCode PCPreSolveChangeRHS_Eisenstat(PC pc, PetscBool *change)
367: {
368: PetscFunctionBegin;
369: *change = PETSC_TRUE;
370: PetscFunctionReturn(PETSC_SUCCESS);
371: }
373: /*MC
374: PCEISENSTAT - An implementation of SSOR (symmetric successive over relaxation, symmetric Gauss-Seidel)
375: preconditioning that incorporates Eisenstat's trick to reduce the amount of computation needed.
377: Options Database Keys:
378: + -pc_eisenstat_omega omega - Sets omega
379: - -pc_eisenstat_no_diagonal_scaling - Activates `PCEisenstatSetNoDiagonalScaling()`
381: Level: beginner
383: Notes:
384: Only implemented for the `MATAIJ` matrix format.
386: Not a true parallel SOR, in parallel this implementation corresponds to block Jacobi with SOR on each block.
388: Developer Note:
389: Since this algorithm runs the Krylov method on a transformed linear system the implementation provides `PCPreSolve()` and `PCPostSolve()`
390: routines that `KSP` uses to set up the transformed linear system.
392: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCEisenstatGetOmega()`,
393: `PCEisenstatSetNoDiagonalScaling()`, `PCEisenstatSetOmega()`, `PCSOR`
394: M*/
396: PETSC_EXTERN PetscErrorCode PCCreate_Eisenstat(PC pc)
397: {
398: PC_Eisenstat *eis;
400: PetscFunctionBegin;
401: PetscCall(PetscNew(&eis));
403: pc->ops->apply = PCApply_Eisenstat;
404: pc->ops->applytranspose = PCApplyTranspose_Eisenstat;
405: pc->ops->presolve = PCPreSolve_Eisenstat;
406: pc->ops->postsolve = PCPostSolve_Eisenstat;
407: pc->ops->applyrichardson = NULL;
408: pc->ops->setfromoptions = PCSetFromOptions_Eisenstat;
409: pc->ops->destroy = PCDestroy_Eisenstat;
410: pc->ops->reset = PCReset_Eisenstat;
411: pc->ops->view = PCView_Eisenstat;
412: pc->ops->setup = PCSetUp_Eisenstat;
414: pc->data = eis;
415: eis->omega = 1.0;
416: eis->b[0] = NULL;
417: eis->b[1] = NULL;
418: eis->diag = NULL;
419: eis->usediag = PETSC_TRUE;
421: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetOmega_C", PCEisenstatSetOmega_Eisenstat));
422: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetNoDiagonalScaling_C", PCEisenstatSetNoDiagonalScaling_Eisenstat));
423: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetOmega_C", PCEisenstatGetOmega_Eisenstat));
424: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetNoDiagonalScaling_C", PCEisenstatGetNoDiagonalScaling_Eisenstat));
425: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", PCPreSolveChangeRHS_Eisenstat));
426: PetscFunctionReturn(PETSC_SUCCESS);
427: }