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: MatShellGetContext(mat, &pc);
23: eis = (PC_Eisenstat *)pc->data;
24: MatSOR(eis->A, b, eis->omega, SOR_EISENSTAT, 0.0, 1, 1, x);
25: MatFactorGetError(eis->A, (MatFactorError *)&pc->failedreason);
26: return 0;
27: }
29: static PetscErrorCode PCNorm_Eisenstat(Mat mat, NormType type, PetscReal *nrm)
30: {
31: PC pc;
32: PC_Eisenstat *eis;
34: MatShellGetContext(mat, &pc);
35: eis = (PC_Eisenstat *)pc->data;
36: MatNorm(eis->A, type, nrm);
37: return 0;
38: }
40: static PetscErrorCode PCApply_Eisenstat(PC pc, Vec x, Vec y)
41: {
42: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
43: PetscBool hasop;
45: if (eis->usediag) {
46: MatHasOperation(pc->pmat, MATOP_MULT_DIAGONAL_BLOCK, &hasop);
47: if (hasop) {
48: MatMultDiagonalBlock(pc->pmat, x, y);
49: } else {
50: VecPointwiseMult(y, x, eis->diag);
51: }
52: } else VecCopy(x, y);
53: return 0;
54: }
56: static PetscErrorCode PCApplyTranspose_Eisenstat(PC pc, Vec x, Vec y)
57: {
58: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
59: PetscBool hasop, set, sym;
61: MatIsSymmetricKnown(eis->A, &set, &sym);
63: if (eis->usediag) {
64: MatHasOperation(pc->pmat, MATOP_MULT_DIAGONAL_BLOCK, &hasop);
65: if (hasop) {
66: MatMultDiagonalBlock(pc->pmat, x, y);
67: } else {
68: VecPointwiseMult(y, x, eis->diag);
69: }
70: } else VecCopy(x, y);
71: return 0;
72: }
74: static PetscErrorCode PCPreSolve_Eisenstat(PC pc, KSP ksp, Vec b, Vec x)
75: {
76: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
77: PetscBool nonzero;
79: if (pc->presolvedone < 2) {
81: /* swap shell matrix and true matrix */
82: eis->A = pc->mat;
83: pc->mat = eis->shell;
84: }
86: if (!eis->b[pc->presolvedone - 1]) { VecDuplicate(b, &eis->b[pc->presolvedone - 1]); }
88: /* if nonzero initial guess, modify x */
89: KSPGetInitialGuessNonzero(ksp, &nonzero);
90: if (nonzero) {
91: VecCopy(x, eis->b[pc->presolvedone - 1]);
92: MatSOR(eis->A, eis->b[pc->presolvedone - 1], eis->omega, SOR_APPLY_UPPER, 0.0, 1, 1, x);
93: MatFactorGetError(eis->A, (MatFactorError *)&pc->failedreason);
94: }
96: /* save true b, other option is to swap pointers */
97: VecCopy(b, eis->b[pc->presolvedone - 1]);
99: /* modify b by (L + D/omega)^{-1} */
100: MatSOR(eis->A, eis->b[pc->presolvedone - 1], eis->omega, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_FORWARD_SWEEP), 0.0, 1, 1, b);
101: MatFactorGetError(eis->A, (MatFactorError *)&pc->failedreason);
102: return 0;
103: }
105: static PetscErrorCode PCPostSolve_Eisenstat(PC pc, KSP ksp, Vec b, Vec x)
106: {
107: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
109: /* get back true b */
110: VecCopy(eis->b[pc->presolvedone], b);
112: /* modify x by (U + D/omega)^{-1} */
113: VecCopy(x, eis->b[pc->presolvedone]);
114: MatSOR(eis->A, eis->b[pc->presolvedone], eis->omega, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_BACKWARD_SWEEP), 0.0, 1, 1, x);
115: MatFactorGetError(eis->A, (MatFactorError *)&pc->failedreason);
116: if (!pc->presolvedone) pc->mat = eis->A;
117: return 0;
118: }
120: static PetscErrorCode PCReset_Eisenstat(PC pc)
121: {
122: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
124: VecDestroy(&eis->b[0]);
125: VecDestroy(&eis->b[1]);
126: MatDestroy(&eis->shell);
127: VecDestroy(&eis->diag);
128: return 0;
129: }
131: static PetscErrorCode PCDestroy_Eisenstat(PC pc)
132: {
133: PCReset_Eisenstat(pc);
134: PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetOmega_C", NULL);
135: PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetNoDiagonalScaling_C", NULL);
136: PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetOmega_C", NULL);
137: PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetNoDiagonalScaling_C", NULL);
138: PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", NULL);
139: PetscFree(pc->data);
140: return 0;
141: }
143: static PetscErrorCode PCSetFromOptions_Eisenstat(PC pc, PetscOptionItems *PetscOptionsObject)
144: {
145: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
146: PetscBool set, flg;
148: PetscOptionsHeadBegin(PetscOptionsObject, "Eisenstat SSOR options");
149: PetscOptionsReal("-pc_eisenstat_omega", "Relaxation factor 0 < omega < 2", "PCEisenstatSetOmega", eis->omega, &eis->omega, NULL);
150: PetscOptionsBool("-pc_eisenstat_no_diagonal_scaling", "Do not use standard diagonal scaling", "PCEisenstatSetNoDiagonalScaling", eis->usediag ? PETSC_FALSE : PETSC_TRUE, &flg, &set);
151: if (set) PCEisenstatSetNoDiagonalScaling(pc, flg);
152: PetscOptionsHeadEnd();
153: return 0;
154: }
156: static PetscErrorCode PCView_Eisenstat(PC pc, PetscViewer viewer)
157: {
158: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
159: PetscBool iascii;
161: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
162: if (iascii) {
163: PetscViewerASCIIPrintf(viewer, " omega = %g\n", (double)eis->omega);
164: if (eis->usediag) {
165: PetscViewerASCIIPrintf(viewer, " Using diagonal scaling (default)\n");
166: } else {
167: PetscViewerASCIIPrintf(viewer, " Not using diagonal scaling\n");
168: }
169: }
170: return 0;
171: }
173: static PetscErrorCode PCSetUp_Eisenstat(PC pc)
174: {
175: PetscInt M, N, m, n;
176: PetscBool set, sym;
177: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
179: if (!pc->setupcalled) {
180: MatGetSize(pc->mat, &M, &N);
181: MatGetLocalSize(pc->mat, &m, &n);
182: MatIsSymmetricKnown(pc->mat, &set, &sym);
183: MatCreate(PetscObjectComm((PetscObject)pc), &eis->shell);
184: MatSetSizes(eis->shell, m, n, M, N);
185: MatSetType(eis->shell, MATSHELL);
186: MatSetUp(eis->shell);
187: MatShellSetContext(eis->shell, pc);
188: MatShellSetOperation(eis->shell, MATOP_MULT, (void (*)(void))PCMult_Eisenstat);
189: if (set && sym) MatShellSetOperation(eis->shell, MATOP_MULT_TRANSPOSE, (void (*)(void))PCMult_Eisenstat);
190: MatShellSetOperation(eis->shell, MATOP_NORM, (void (*)(void))PCNorm_Eisenstat);
191: }
192: if (!eis->usediag) return 0;
193: if (!pc->setupcalled) { MatCreateVecs(pc->pmat, &eis->diag, NULL); }
194: MatGetDiagonal(pc->pmat, eis->diag);
195: return 0;
196: }
198: /* --------------------------------------------------------------------*/
200: static PetscErrorCode PCEisenstatSetOmega_Eisenstat(PC pc, PetscReal omega)
201: {
202: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
205: eis->omega = omega;
206: return 0;
207: }
209: static PetscErrorCode PCEisenstatSetNoDiagonalScaling_Eisenstat(PC pc, PetscBool flg)
210: {
211: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
213: eis->usediag = flg;
214: return 0;
215: }
217: static PetscErrorCode PCEisenstatGetOmega_Eisenstat(PC pc, PetscReal *omega)
218: {
219: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
221: *omega = eis->omega;
222: return 0;
223: }
225: static PetscErrorCode PCEisenstatGetNoDiagonalScaling_Eisenstat(PC pc, PetscBool *flg)
226: {
227: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
229: *flg = eis->usediag;
230: return 0;
231: }
233: /*@
234: PCEisenstatSetOmega - Sets the SSOR relaxation coefficient, omega,
235: to use with Eisenstat's trick (where omega = 1.0 by default)
237: Logically Collective
239: Input Parameters:
240: + pc - the preconditioner context
241: - omega - relaxation coefficient (0 < omega < 2)
243: Options Database Key:
244: . -pc_eisenstat_omega <omega> - Sets omega
246: Notes:
247: The Eisenstat trick implementation of SSOR requires about 50% of the
248: usual amount of floating point operations used for SSOR + Krylov method;
249: however, the preconditioned problem must be solved with both left
250: and right preconditioning.
252: To use SSOR without the Eisenstat trick, employ the `PCSOR` preconditioner,
253: which can be chosen with the database options
254: $ -pc_type sor -pc_sor_symmetric
256: Level: intermediate
258: .seealso: `PCSORSetOmega()`, `PCEISENSTAT`
259: @*/
260: PetscErrorCode PCEisenstatSetOmega(PC pc, PetscReal omega)
261: {
264: PetscTryMethod(pc, "PCEisenstatSetOmega_C", (PC, PetscReal), (pc, omega));
265: return 0;
266: }
268: /*@
269: PCEisenstatSetNoDiagonalScaling - Causes the Eisenstat preconditioner, `PCEISENSTAT`
270: not to do additional diagonal preconditioning. For matrices with a constant
271: along the diagonal, this may save a small amount of work.
273: Logically Collective
275: Input Parameters:
276: + pc - the preconditioner context
277: - flg - `PETSC_TRUE` turns off diagonal scaling inside the algorithm
279: Options Database Key:
280: . -pc_eisenstat_no_diagonal_scaling - Activates `PCEisenstatSetNoDiagonalScaling()`
282: Level: intermediate
284: Note:
285: If you use the `KSPSetDiagonalScaling()` or -ksp_diagonal_scale option then you will
286: likely want to use this routine since it will save you some unneeded flops.
288: .seealso: `PCEisenstatSetOmega()`, `PCEISENSTAT`
289: @*/
290: PetscErrorCode PCEisenstatSetNoDiagonalScaling(PC pc, PetscBool flg)
291: {
293: PetscTryMethod(pc, "PCEisenstatSetNoDiagonalScaling_C", (PC, PetscBool), (pc, flg));
294: return 0;
295: }
297: /*@
298: PCEisenstatGetOmega - Gets the SSOR relaxation coefficient, omega,
299: to use with Eisenstat's trick (where omega = 1.0 by default).
301: Logically Collective
303: Input Parameter:
304: . pc - the preconditioner context
306: Output Parameter:
307: . omega - relaxation coefficient (0 < omega < 2)
309: Options Database Key:
310: . -pc_eisenstat_omega <omega> - Sets omega
312: Notes:
313: The Eisenstat trick implementation of SSOR requires about 50% of the
314: usual amount of floating point operations used for SSOR + Krylov method;
315: however, the preconditioned problem must be solved with both left
316: and right preconditioning.
318: To use SSOR without the Eisenstat trick, employ the PCSOR preconditioner,
319: which can be chosen with the database options
320: $ -pc_type sor -pc_sor_symmetric
322: Level: intermediate
324: .seealso: `PCEISENSTAT`, `PCSORGetOmega()`, `PCEisenstatSetOmega()`
325: @*/
326: PetscErrorCode PCEisenstatGetOmega(PC pc, PetscReal *omega)
327: {
329: PetscUseMethod(pc, "PCEisenstatGetOmega_C", (PC, PetscReal *), (pc, omega));
330: return 0;
331: }
333: /*@
334: PCEisenstatGetNoDiagonalScaling - Tells if the Eisenstat preconditioner
335: not to do additional diagonal preconditioning. For matrices with a constant
336: along the diagonal, this may save a small amount of work.
338: Logically Collective
340: Input Parameter:
341: . pc - the preconditioner context
343: Output Parameter:
344: . flg - `PETSC_TRUE` means there is no diagonal scaling applied
346: Options Database Key:
347: . -pc_eisenstat_no_diagonal_scaling - Activates `PCEisenstatSetNoDiagonalScaling()`
349: Level: intermediate
351: Note:
352: If you use the KSPSetDiagonalScaling() or -ksp_diagonal_scale option then you will
353: likely want to use this routine since it will save you some unneeded flops.
355: .seealso: , `PCEISENSTAT`, `PCEisenstatGetOmega()`
356: @*/
357: PetscErrorCode PCEisenstatGetNoDiagonalScaling(PC pc, PetscBool *flg)
358: {
360: PetscUseMethod(pc, "PCEisenstatGetNoDiagonalScaling_C", (PC, PetscBool *), (pc, flg));
361: return 0;
362: }
364: static PetscErrorCode PCPreSolveChangeRHS_Eisenstat(PC pc, PetscBool *change)
365: {
366: *change = PETSC_TRUE;
367: return 0;
368: }
370: /*MC
371: PCEISENSTAT - An implementation of SSOR (symmetric successive over relaxation, symmetric Gauss-Seidel)
372: preconditioning that incorporates Eisenstat's trick to reduce the amount of computation needed.
374: Options Database Keys:
375: + -pc_eisenstat_omega <omega> - Sets omega
376: - -pc_eisenstat_no_diagonal_scaling - Activates `PCEisenstatSetNoDiagonalScaling()`
378: Level: beginner
380: Notes:
381: Only implemented for the `MATAIJ` matrix format.
383: Not a true parallel SOR, in parallel this implementation corresponds to block Jacobi with SOR on each block.
385: Developer Note:
386: Since this algorithm runs the Krylov method on a transformed linear system the implementation provides `PCPreSolve()` and `PCPostSolve()`
387: routines that `KSP` uses to set up the transformed linear system.
389: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCEisenstatGetOmega()`,
390: `PCEisenstatSetNoDiagonalScaling()`, `PCEisenstatSetOmega()`, `PCSOR`
391: M*/
393: PETSC_EXTERN PetscErrorCode PCCreate_Eisenstat(PC pc)
394: {
395: PC_Eisenstat *eis;
397: PetscNew(&eis);
399: pc->ops->apply = PCApply_Eisenstat;
400: pc->ops->applytranspose = PCApplyTranspose_Eisenstat;
401: pc->ops->presolve = PCPreSolve_Eisenstat;
402: pc->ops->postsolve = PCPostSolve_Eisenstat;
403: pc->ops->applyrichardson = NULL;
404: pc->ops->setfromoptions = PCSetFromOptions_Eisenstat;
405: pc->ops->destroy = PCDestroy_Eisenstat;
406: pc->ops->reset = PCReset_Eisenstat;
407: pc->ops->view = PCView_Eisenstat;
408: pc->ops->setup = PCSetUp_Eisenstat;
410: pc->data = eis;
411: eis->omega = 1.0;
412: eis->b[0] = NULL;
413: eis->b[1] = NULL;
414: eis->diag = NULL;
415: eis->usediag = PETSC_TRUE;
417: PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetOmega_C", PCEisenstatSetOmega_Eisenstat);
418: PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetNoDiagonalScaling_C", PCEisenstatSetNoDiagonalScaling_Eisenstat);
419: PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetOmega_C", PCEisenstatGetOmega_Eisenstat);
420: PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetNoDiagonalScaling_C", PCEisenstatGetNoDiagonalScaling_Eisenstat);
421: PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", PCPreSolveChangeRHS_Eisenstat);
422: return 0;
423: }