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: Note:
293: If you use the `KSPSetDiagonalScaling()` or -ksp_diagonal_scale option then you will
294: likely want to use this routine since it will save you some unneeded flops.
296: .seealso: [](ch_ksp), `PCEisenstatSetOmega()`, `PCEISENSTAT`
297: @*/
298: PetscErrorCode PCEisenstatSetNoDiagonalScaling(PC pc, PetscBool flg)
299: {
300: PetscFunctionBegin;
302: PetscTryMethod(pc, "PCEisenstatSetNoDiagonalScaling_C", (PC, PetscBool), (pc, flg));
303: PetscFunctionReturn(PETSC_SUCCESS);
304: }
306: /*@
307: PCEisenstatGetOmega - Gets the SSOR relaxation coefficient, omega,
308: to use with Eisenstat's trick (where omega = 1.0 by default).
310: Logically Collective
312: Input Parameter:
313: . pc - the preconditioner context
315: Output Parameter:
316: . omega - relaxation coefficient (0 < omega < 2)
318: Options Database Key:
319: . -pc_eisenstat_omega omega - Sets omega
321: Notes:
322: The Eisenstat trick implementation of SSOR requires about 50% of the
323: usual amount of floating point operations used for SSOR + Krylov method;
324: however, the preconditioned problem must be solved with both left
325: and right preconditioning.
327: To use SSOR without the Eisenstat trick, employ the PCSOR preconditioner,
328: which can be chosen with the database options `-pc_type sor -pc_sor_symmetric`
330: Level: intermediate
332: .seealso: [](ch_ksp), `PCEISENSTAT`, `PCSORGetOmega()`, `PCEisenstatSetOmega()`
333: @*/
334: PetscErrorCode PCEisenstatGetOmega(PC pc, PetscReal *omega)
335: {
336: PetscFunctionBegin;
338: PetscUseMethod(pc, "PCEisenstatGetOmega_C", (PC, PetscReal *), (pc, omega));
339: PetscFunctionReturn(PETSC_SUCCESS);
340: }
342: /*@
343: PCEisenstatGetNoDiagonalScaling - Tells if the Eisenstat preconditioner
344: not to do additional diagonal preconditioning. For matrices with a constant
345: along the diagonal, this may save a small amount of work.
347: Logically Collective
349: Input Parameter:
350: . pc - the preconditioner context
352: Output Parameter:
353: . flg - `PETSC_TRUE` means there is no diagonal scaling applied
355: Options Database Key:
356: . -pc_eisenstat_no_diagonal_scaling - Activates `PCEisenstatSetNoDiagonalScaling()`
358: Level: intermediate
360: Note:
361: If you use the KSPSetDiagonalScaling() or -ksp_diagonal_scale option then you will
362: likely want to use this routine since it will save you some unneeded flops.
364: .seealso: `PCEISENSTAT`, `PCEisenstatGetOmega()`
365: @*/
366: PetscErrorCode PCEisenstatGetNoDiagonalScaling(PC pc, PetscBool *flg)
367: {
368: PetscFunctionBegin;
370: PetscUseMethod(pc, "PCEisenstatGetNoDiagonalScaling_C", (PC, PetscBool *), (pc, flg));
371: PetscFunctionReturn(PETSC_SUCCESS);
372: }
374: static PetscErrorCode PCPreSolveChangeRHS_Eisenstat(PC pc, PetscBool *change)
375: {
376: PetscFunctionBegin;
377: *change = PETSC_TRUE;
378: PetscFunctionReturn(PETSC_SUCCESS);
379: }
381: /*MC
382: PCEISENSTAT - An implementation of SSOR (symmetric successive over relaxation, symmetric Gauss-Seidel)
383: preconditioning that incorporates Eisenstat's trick to reduce the amount of computation needed.
385: Options Database Keys:
386: + -pc_eisenstat_omega omega - Sets omega
387: - -pc_eisenstat_no_diagonal_scaling - Activates `PCEisenstatSetNoDiagonalScaling()`
389: Level: beginner
391: Notes:
392: Only implemented for the `MATAIJ` matrix format.
394: Not a true parallel SOR, in parallel this implementation corresponds to block Jacobi with SOR on each block.
396: Developer Note:
397: Since this algorithm runs the Krylov method on a transformed linear system the implementation provides `PCPreSolve()` and `PCPostSolve()`
398: routines that `KSP` uses to set up the transformed linear system.
400: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCEisenstatGetOmega()`,
401: `PCEisenstatSetNoDiagonalScaling()`, `PCEisenstatSetOmega()`, `PCSOR`
402: M*/
404: PETSC_EXTERN PetscErrorCode PCCreate_Eisenstat(PC pc)
405: {
406: PC_Eisenstat *eis;
408: PetscFunctionBegin;
409: PetscCall(PetscNew(&eis));
411: pc->ops->apply = PCApply_Eisenstat;
412: pc->ops->applytranspose = PCApplyTranspose_Eisenstat;
413: pc->ops->presolve = PCPreSolve_Eisenstat;
414: pc->ops->postsolve = PCPostSolve_Eisenstat;
415: pc->ops->applyrichardson = NULL;
416: pc->ops->setfromoptions = PCSetFromOptions_Eisenstat;
417: pc->ops->destroy = PCDestroy_Eisenstat;
418: pc->ops->reset = PCReset_Eisenstat;
419: pc->ops->view = PCView_Eisenstat;
420: pc->ops->setup = PCSetUp_Eisenstat;
422: pc->data = eis;
423: eis->omega = 1.0;
424: eis->b[0] = NULL;
425: eis->b[1] = NULL;
426: eis->diag = NULL;
427: eis->usediag = PETSC_TRUE;
429: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetOmega_C", PCEisenstatSetOmega_Eisenstat));
430: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetNoDiagonalScaling_C", PCEisenstatSetNoDiagonalScaling_Eisenstat));
431: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetOmega_C", PCEisenstatGetOmega_Eisenstat));
432: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetNoDiagonalScaling_C", PCEisenstatGetNoDiagonalScaling_Eisenstat));
433: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", PCPreSolveChangeRHS_Eisenstat));
434: PetscFunctionReturn(PETSC_SUCCESS);
435: }