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> /*I "petscpc.h" I*/
8
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;
15
PCMult_Eisenstat(Mat mat,Vec b,Vec x)16 static PetscErrorCode PCMult_Eisenstat(Mat mat, Vec b, Vec x)
17 {
18 PC pc;
19 PC_Eisenstat *eis;
20
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 }
28
PCNorm_Eisenstat(Mat mat,NormType type,PetscReal * nrm)29 static PetscErrorCode PCNorm_Eisenstat(Mat mat, NormType type, PetscReal *nrm)
30 {
31 PC pc;
32 PC_Eisenstat *eis;
33
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 }
40
PCApply_Eisenstat(PC pc,Vec x,Vec y)41 static PetscErrorCode PCApply_Eisenstat(PC pc, Vec x, Vec y)
42 {
43 PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
44 PetscBool hasop;
45
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 }
57
PCApplyTranspose_Eisenstat(PC pc,Vec x,Vec y)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;
62
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 }
76
PCPreSolve_Eisenstat(PC pc,KSP ksp,Vec b,Vec x)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;
81
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 }
89
90 if (!eis->b[pc->presolvedone - 1]) PetscCall(VecDuplicate(b, &eis->b[pc->presolvedone - 1]));
91
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 }
99
100 /* save true b, other option is to swap pointers */
101 PetscCall(VecCopy(b, eis->b[pc->presolvedone - 1]));
102
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 }
108
PCPostSolve_Eisenstat(PC pc,KSP ksp,Vec b,Vec x)109 static PetscErrorCode PCPostSolve_Eisenstat(PC pc, KSP ksp, Vec b, Vec x)
110 {
111 PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
112
113 PetscFunctionBegin;
114 /* get back true b */
115 PetscCall(VecCopy(eis->b[pc->presolvedone], b));
116
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 }
124
PCReset_Eisenstat(PC pc)125 static PetscErrorCode PCReset_Eisenstat(PC pc)
126 {
127 PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
128
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 }
136
PCDestroy_Eisenstat(PC pc)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 }
149
PCSetFromOptions_Eisenstat(PC pc,PetscOptionItems PetscOptionsObject)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;
155
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 }
165
PCView_Eisenstat(PC pc,PetscViewer viewer)166 static PetscErrorCode PCView_Eisenstat(PC pc, PetscViewer viewer)
167 {
168 PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
169 PetscBool isascii;
170
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 }
183
PCSetUp_Eisenstat(PC pc)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;
189
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 }
209
PCEisenstatSetOmega_Eisenstat(PC pc,PetscReal omega)210 static PetscErrorCode PCEisenstatSetOmega_Eisenstat(PC pc, PetscReal omega)
211 {
212 PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
213
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 }
219
PCEisenstatSetNoDiagonalScaling_Eisenstat(PC pc,PetscBool flg)220 static PetscErrorCode PCEisenstatSetNoDiagonalScaling_Eisenstat(PC pc, PetscBool flg)
221 {
222 PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
223
224 PetscFunctionBegin;
225 eis->usediag = flg;
226 PetscFunctionReturn(PETSC_SUCCESS);
227 }
228
PCEisenstatGetOmega_Eisenstat(PC pc,PetscReal * omega)229 static PetscErrorCode PCEisenstatGetOmega_Eisenstat(PC pc, PetscReal *omega)
230 {
231 PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
232
233 PetscFunctionBegin;
234 *omega = eis->omega;
235 PetscFunctionReturn(PETSC_SUCCESS);
236 }
237
PCEisenstatGetNoDiagonalScaling_Eisenstat(PC pc,PetscBool * flg)238 static PetscErrorCode PCEisenstatGetNoDiagonalScaling_Eisenstat(PC pc, PetscBool *flg)
239 {
240 PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
241
242 PetscFunctionBegin;
243 *flg = eis->usediag;
244 PetscFunctionReturn(PETSC_SUCCESS);
245 }
246
247 /*@
248 PCEisenstatSetOmega - Sets the SSOR relaxation coefficient, omega,
249 to use with Eisenstat's trick (where omega = 1.0 by default)
250
251 Logically Collective
252
253 Input Parameters:
254 + pc - the preconditioner context
255 - omega - relaxation coefficient (0 < omega < 2)
256
257 Options Database Key:
258 . -pc_eisenstat_omega <omega> - Sets omega
259
260 Level: intermediate
261
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.
267
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`
270
271 .seealso: [](ch_ksp), `PCSORSetOmega()`, `PCEISENSTAT`
272 @*/
PCEisenstatSetOmega(PC pc,PetscReal omega)273 PetscErrorCode PCEisenstatSetOmega(PC pc, PetscReal omega)
274 {
275 PetscFunctionBegin;
276 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
277 PetscValidLogicalCollectiveReal(pc, omega, 2);
278 PetscTryMethod(pc, "PCEisenstatSetOmega_C", (PC, PetscReal), (pc, omega));
279 PetscFunctionReturn(PETSC_SUCCESS);
280 }
281
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.
286
287 Logically Collective
288
289 Input Parameters:
290 + pc - the preconditioner context
291 - flg - `PETSC_TRUE` turns off diagonal scaling inside the algorithm
292
293 Options Database Key:
294 . -pc_eisenstat_no_diagonal_scaling - Activates `PCEisenstatSetNoDiagonalScaling()`
295
296 Level: intermediate
297
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.
301
302 .seealso: [](ch_ksp), `PCEisenstatSetOmega()`, `PCEISENSTAT`
303 @*/
PCEisenstatSetNoDiagonalScaling(PC pc,PetscBool flg)304 PetscErrorCode PCEisenstatSetNoDiagonalScaling(PC pc, PetscBool flg)
305 {
306 PetscFunctionBegin;
307 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
308 PetscTryMethod(pc, "PCEisenstatSetNoDiagonalScaling_C", (PC, PetscBool), (pc, flg));
309 PetscFunctionReturn(PETSC_SUCCESS);
310 }
311
312 /*@
313 PCEisenstatGetOmega - Gets the SSOR relaxation coefficient, omega,
314 to use with Eisenstat's trick (where omega = 1.0 by default).
315
316 Logically Collective
317
318 Input Parameter:
319 . pc - the preconditioner context
320
321 Output Parameter:
322 . omega - relaxation coefficient (0 < omega < 2)
323
324 Options Database Key:
325 . -pc_eisenstat_omega <omega> - Sets omega
326
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.
332
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`
335
336 Level: intermediate
337
338 .seealso: [](ch_ksp), `PCEISENSTAT`, `PCSORGetOmega()`, `PCEisenstatSetOmega()`
339 @*/
PCEisenstatGetOmega(PC pc,PetscReal * omega)340 PetscErrorCode PCEisenstatGetOmega(PC pc, PetscReal *omega)
341 {
342 PetscFunctionBegin;
343 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
344 PetscUseMethod(pc, "PCEisenstatGetOmega_C", (PC, PetscReal *), (pc, omega));
345 PetscFunctionReturn(PETSC_SUCCESS);
346 }
347
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.
352
353 Logically Collective
354
355 Input Parameter:
356 . pc - the preconditioner context
357
358 Output Parameter:
359 . flg - `PETSC_TRUE` means there is no diagonal scaling applied
360
361 Options Database Key:
362 . -pc_eisenstat_no_diagonal_scaling - Activates `PCEisenstatSetNoDiagonalScaling()`
363
364 Level: intermediate
365
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.
369
370 .seealso: , `PCEISENSTAT`, `PCEisenstatGetOmega()`
371 @*/
PCEisenstatGetNoDiagonalScaling(PC pc,PetscBool * flg)372 PetscErrorCode PCEisenstatGetNoDiagonalScaling(PC pc, PetscBool *flg)
373 {
374 PetscFunctionBegin;
375 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
376 PetscUseMethod(pc, "PCEisenstatGetNoDiagonalScaling_C", (PC, PetscBool *), (pc, flg));
377 PetscFunctionReturn(PETSC_SUCCESS);
378 }
379
PCPreSolveChangeRHS_Eisenstat(PC pc,PetscBool * change)380 static PetscErrorCode PCPreSolveChangeRHS_Eisenstat(PC pc, PetscBool *change)
381 {
382 PetscFunctionBegin;
383 *change = PETSC_TRUE;
384 PetscFunctionReturn(PETSC_SUCCESS);
385 }
386
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.
390
391 Options Database Keys:
392 + -pc_eisenstat_omega <omega> - Sets omega
393 - -pc_eisenstat_no_diagonal_scaling - Activates `PCEisenstatSetNoDiagonalScaling()`
394
395 Level: beginner
396
397 Notes:
398 Only implemented for the `MATAIJ` matrix format.
399
400 Not a true parallel SOR, in parallel this implementation corresponds to block Jacobi with SOR on each block.
401
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.
405
406 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCEisenstatGetOmega()`,
407 `PCEisenstatSetNoDiagonalScaling()`, `PCEisenstatSetOmega()`, `PCSOR`
408 M*/
409
PCCreate_Eisenstat(PC pc)410 PETSC_EXTERN PetscErrorCode PCCreate_Eisenstat(PC pc)
411 {
412 PC_Eisenstat *eis;
413
414 PetscFunctionBegin;
415 PetscCall(PetscNew(&eis));
416
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;
427
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;
434
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 }
442