Actual source code: sor.c
1: /*
2: Defines a (S)SOR preconditioner for any Mat implementation
3: */
4: #include <petsc/private/pcimpl.h>
6: typedef struct {
7: PetscInt its; /* inner iterations, number of sweeps */
8: PetscInt lits; /* local inner iterations, number of sweeps applied by the local matrix mat->A */
9: MatSORType sym; /* forward, reverse, symmetric etc. */
10: PetscReal omega;
11: PetscReal fshift;
12: } PC_SOR;
14: static PetscErrorCode PCDestroy_SOR(PC pc)
15: {
16: PetscFunctionBegin;
17: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORSetSymmetric_C", NULL));
18: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORSetOmega_C", NULL));
19: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORSetIterations_C", NULL));
20: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORGetSymmetric_C", NULL));
21: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORGetOmega_C", NULL));
22: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORGetIterations_C", NULL));
23: PetscCall(PetscFree(pc->data));
24: PetscFunctionReturn(PETSC_SUCCESS);
25: }
27: static PetscErrorCode PCApply_SOR(PC pc, Vec x, Vec y)
28: {
29: PC_SOR *jac = (PC_SOR *)pc->data;
30: PetscInt flag = jac->sym | SOR_ZERO_INITIAL_GUESS;
32: PetscFunctionBegin;
33: PetscCall(MatSOR(pc->pmat, x, jac->omega, (MatSORType)flag, jac->fshift, jac->its, jac->lits, y));
34: PetscCall(MatFactorGetError(pc->pmat, (MatFactorError *)&pc->failedreason));
35: PetscFunctionReturn(PETSC_SUCCESS);
36: }
38: static PetscErrorCode PCApplyTranspose_SOR(PC pc, Vec x, Vec y)
39: {
40: PC_SOR *jac = (PC_SOR *)pc->data;
41: PetscInt flag = jac->sym | SOR_ZERO_INITIAL_GUESS;
42: PetscBool set, sym;
44: PetscFunctionBegin;
45: PetscCall(MatIsSymmetricKnown(pc->pmat, &set, &sym));
46: PetscCheck(set && sym && (jac->sym == SOR_SYMMETRIC_SWEEP || jac->sym == SOR_LOCAL_SYMMETRIC_SWEEP), PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Can only apply transpose of SOR if matrix is symmetric and sweep is symmetric");
47: PetscCall(MatSOR(pc->pmat, x, jac->omega, (MatSORType)flag, jac->fshift, jac->its, jac->lits, y));
48: PetscCall(MatFactorGetError(pc->pmat, (MatFactorError *)&pc->failedreason));
49: PetscFunctionReturn(PETSC_SUCCESS);
50: }
52: static PetscErrorCode PCApplyRichardson_SOR(PC pc, Vec b, Vec y, Vec w, PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt its, PetscBool guesszero, PetscInt *outits, PCRichardsonConvergedReason *reason)
53: {
54: PC_SOR *jac = (PC_SOR *)pc->data;
55: MatSORType stype = jac->sym;
57: PetscFunctionBegin;
58: PetscCall(PetscInfo(pc, "Warning, convergence criteria ignored, using %" PetscInt_FMT " iterations\n", its));
59: if (guesszero) stype = (MatSORType)(stype | SOR_ZERO_INITIAL_GUESS);
60: PetscCall(MatSOR(pc->pmat, b, jac->omega, stype, jac->fshift, its * jac->its, jac->lits, y));
61: PetscCall(MatFactorGetError(pc->pmat, (MatFactorError *)&pc->failedreason));
62: *outits = its;
63: *reason = PCRICHARDSON_CONVERGED_ITS;
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
67: static PetscErrorCode PCSetFromOptions_SOR(PC pc, PetscOptionItems *PetscOptionsObject)
68: {
69: PC_SOR *jac = (PC_SOR *)pc->data;
70: PetscBool flg;
71: PetscReal omega;
73: PetscFunctionBegin;
74: PetscOptionsHeadBegin(PetscOptionsObject, "(S)SOR options");
75: PetscCall(PetscOptionsReal("-pc_sor_omega", "relaxation factor (0 < omega < 2)", "PCSORSetOmega", jac->omega, &omega, &flg));
76: if (flg) PetscCall(PCSORSetOmega(pc, omega));
77: PetscCall(PetscOptionsReal("-pc_sor_diagonal_shift", "Add to the diagonal entries", "", jac->fshift, &jac->fshift, NULL));
78: PetscCall(PetscOptionsInt("-pc_sor_its", "number of inner SOR iterations", "PCSORSetIterations", jac->its, &jac->its, NULL));
79: PetscCall(PetscOptionsInt("-pc_sor_lits", "number of local inner SOR iterations", "PCSORSetIterations", jac->lits, &jac->lits, NULL));
80: PetscCall(PetscOptionsBoolGroupBegin("-pc_sor_symmetric", "SSOR, not SOR", "PCSORSetSymmetric", &flg));
81: if (flg) PetscCall(PCSORSetSymmetric(pc, SOR_SYMMETRIC_SWEEP));
82: PetscCall(PetscOptionsBoolGroup("-pc_sor_backward", "use backward sweep instead of forward", "PCSORSetSymmetric", &flg));
83: if (flg) PetscCall(PCSORSetSymmetric(pc, SOR_BACKWARD_SWEEP));
84: PetscCall(PetscOptionsBoolGroup("-pc_sor_forward", "use forward sweep", "PCSORSetSymmetric", &flg));
85: if (flg) PetscCall(PCSORSetSymmetric(pc, SOR_FORWARD_SWEEP));
86: PetscCall(PetscOptionsBoolGroup("-pc_sor_local_symmetric", "use SSOR separately on each processor", "PCSORSetSymmetric", &flg));
87: if (flg) PetscCall(PCSORSetSymmetric(pc, SOR_LOCAL_SYMMETRIC_SWEEP));
88: PetscCall(PetscOptionsBoolGroup("-pc_sor_local_backward", "use backward sweep locally", "PCSORSetSymmetric", &flg));
89: if (flg) PetscCall(PCSORSetSymmetric(pc, SOR_LOCAL_BACKWARD_SWEEP));
90: PetscCall(PetscOptionsBoolGroupEnd("-pc_sor_local_forward", "use forward sweep locally", "PCSORSetSymmetric", &flg));
91: if (flg) PetscCall(PCSORSetSymmetric(pc, SOR_LOCAL_FORWARD_SWEEP));
92: PetscOptionsHeadEnd();
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
96: static PetscErrorCode PCView_SOR(PC pc, PetscViewer viewer)
97: {
98: PC_SOR *jac = (PC_SOR *)pc->data;
99: MatSORType sym = jac->sym;
100: const char *sortype;
101: PetscBool iascii;
103: PetscFunctionBegin;
104: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
105: if (iascii) {
106: if (sym & SOR_ZERO_INITIAL_GUESS) PetscCall(PetscViewerASCIIPrintf(viewer, " zero initial guess\n"));
107: if (sym == SOR_APPLY_UPPER) sortype = "apply_upper";
108: else if (sym == SOR_APPLY_LOWER) sortype = "apply_lower";
109: else if (sym & SOR_EISENSTAT) sortype = "Eisenstat";
110: else if ((sym & SOR_SYMMETRIC_SWEEP) == SOR_SYMMETRIC_SWEEP) sortype = "symmetric";
111: else if (sym & SOR_BACKWARD_SWEEP) sortype = "backward";
112: else if (sym & SOR_FORWARD_SWEEP) sortype = "forward";
113: else if ((sym & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) sortype = "local_symmetric";
114: else if (sym & SOR_LOCAL_FORWARD_SWEEP) sortype = "local_forward";
115: else if (sym & SOR_LOCAL_BACKWARD_SWEEP) sortype = "local_backward";
116: else sortype = "unknown";
117: PetscCall(PetscViewerASCIIPrintf(viewer, " type = %s, iterations = %" PetscInt_FMT ", local iterations = %" PetscInt_FMT ", omega = %g\n", sortype, jac->its, jac->lits, (double)jac->omega));
118: }
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: static PetscErrorCode PCSORSetSymmetric_SOR(PC pc, MatSORType flag)
123: {
124: PC_SOR *jac = (PC_SOR *)pc->data;
126: PetscFunctionBegin;
127: jac->sym = flag;
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: static PetscErrorCode PCSORSetOmega_SOR(PC pc, PetscReal omega)
132: {
133: PC_SOR *jac = (PC_SOR *)pc->data;
135: PetscFunctionBegin;
136: PetscCheck(omega > 0.0 && omega < 2.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Relaxation out of range");
137: jac->omega = omega;
138: PetscFunctionReturn(PETSC_SUCCESS);
139: }
141: static PetscErrorCode PCSORSetIterations_SOR(PC pc, PetscInt its, PetscInt lits)
142: {
143: PC_SOR *jac = (PC_SOR *)pc->data;
145: PetscFunctionBegin;
146: jac->its = its;
147: jac->lits = lits;
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: static PetscErrorCode PCSORGetSymmetric_SOR(PC pc, MatSORType *flag)
152: {
153: PC_SOR *jac = (PC_SOR *)pc->data;
155: PetscFunctionBegin;
156: *flag = jac->sym;
157: PetscFunctionReturn(PETSC_SUCCESS);
158: }
160: static PetscErrorCode PCSORGetOmega_SOR(PC pc, PetscReal *omega)
161: {
162: PC_SOR *jac = (PC_SOR *)pc->data;
164: PetscFunctionBegin;
165: *omega = jac->omega;
166: PetscFunctionReturn(PETSC_SUCCESS);
167: }
169: static PetscErrorCode PCSORGetIterations_SOR(PC pc, PetscInt *its, PetscInt *lits)
170: {
171: PC_SOR *jac = (PC_SOR *)pc->data;
173: PetscFunctionBegin;
174: if (its) *its = jac->its;
175: if (lits) *lits = jac->lits;
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: /*@
180: PCSORGetSymmetric - Gets the form the SOR preconditioner is using; backward, or forward relaxation. The local variants perform SOR on
181: each processor. By default forward relaxation is used.
183: Logically Collective
185: Input Parameter:
186: . pc - the preconditioner context
188: Output Parameter:
189: . flag - one of the following
190: .vb
191: SOR_FORWARD_SWEEP
192: SOR_BACKWARD_SWEEP
193: SOR_SYMMETRIC_SWEEP
194: SOR_LOCAL_FORWARD_SWEEP
195: SOR_LOCAL_BACKWARD_SWEEP
196: SOR_LOCAL_SYMMETRIC_SWEEP
197: .ve
199: Options Database Keys:
200: + -pc_sor_symmetric - Activates symmetric version
201: . -pc_sor_backward - Activates backward version
202: . -pc_sor_local_forward - Activates local forward version
203: . -pc_sor_local_symmetric - Activates local symmetric version
204: - -pc_sor_local_backward - Activates local backward version
206: Note:
207: To use the Eisenstat trick with SSOR, employ the `PCEISENSTAT` preconditioner,
208: which can be chosen with the option
209: . -pc_type eisenstat - Activates Eisenstat trick
211: Level: intermediate
213: .seealso: [](ch_ksp), `PCSOR`, `PCEisenstatSetOmega()`, `PCSORSetIterations()`, `PCSORSetOmega()`, `PCSORSetSymmetric()`
214: @*/
215: PetscErrorCode PCSORGetSymmetric(PC pc, MatSORType *flag)
216: {
217: PetscFunctionBegin;
219: PetscUseMethod(pc, "PCSORGetSymmetric_C", (PC, MatSORType *), (pc, flag));
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
223: /*@
224: PCSORGetOmega - Gets the SOR relaxation coefficient, omega
225: (where omega = 1.0 by default).
227: Logically Collective
229: Input Parameter:
230: . pc - the preconditioner context
232: Output Parameter:
233: . omega - relaxation coefficient (0 < omega < 2).
235: Options Database Key:
236: . -pc_sor_omega <omega> - Sets omega
238: Level: intermediate
240: .seealso: [](ch_ksp), `PCSOR`, `PCSORSetSymmetric()`, `PCSORSetIterations()`, `PCEisenstatSetOmega()`, `PCSORSetOmega()`
241: @*/
242: PetscErrorCode PCSORGetOmega(PC pc, PetscReal *omega)
243: {
244: PetscFunctionBegin;
246: PetscUseMethod(pc, "PCSORGetOmega_C", (PC, PetscReal *), (pc, omega));
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
250: /*@
251: PCSORGetIterations - Gets the number of inner iterations to
252: be used by the SOR preconditioner. The default is 1.
254: Logically Collective
256: Input Parameter:
257: . pc - the preconditioner context
259: Output Parameters:
260: + lits - number of local iterations, smoothings over just variables on processor
261: - its - number of parallel iterations to use; each parallel iteration has lits local iterations
263: Options Database Keys:
264: + -pc_sor_its <its> - Sets number of iterations
265: - -pc_sor_lits <lits> - Sets number of local iterations
267: Level: intermediate
269: Note:
270: When run on one processor the number of smoothings is lits*its
272: .seealso: [](ch_ksp), `PCSOR`, `PCSORSetOmega()`, `PCSORSetSymmetric()`, `PCSORSetIterations()`
273: @*/
274: PetscErrorCode PCSORGetIterations(PC pc, PetscInt *its, PetscInt *lits)
275: {
276: PetscFunctionBegin;
278: PetscUseMethod(pc, "PCSORGetIterations_C", (PC, PetscInt *, PetscInt *), (pc, its, lits));
279: PetscFunctionReturn(PETSC_SUCCESS);
280: }
282: /*@
283: PCSORSetSymmetric - Sets the SOR preconditioner to use symmetric (SSOR),
284: backward, or forward relaxation. The local variants perform SOR on
285: each processor. By default forward relaxation is used.
287: Logically Collective
289: Input Parameters:
290: + pc - the preconditioner context
291: - flag - one of the following
292: .vb
293: SOR_FORWARD_SWEEP
294: SOR_BACKWARD_SWEEP
295: SOR_SYMMETRIC_SWEEP
296: SOR_LOCAL_FORWARD_SWEEP
297: SOR_LOCAL_BACKWARD_SWEEP
298: SOR_LOCAL_SYMMETRIC_SWEEP
299: .ve
301: Options Database Keys:
302: + -pc_sor_symmetric - Activates symmetric version
303: . -pc_sor_backward - Activates backward version
304: . -pc_sor_local_forward - Activates local forward version
305: . -pc_sor_local_symmetric - Activates local symmetric version
306: - -pc_sor_local_backward - Activates local backward version
308: Note:
309: To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner,
310: which can be chosen with the option
311: . -pc_type eisenstat - Activates Eisenstat trick
313: Level: intermediate
315: .seealso: [](ch_ksp), `PCSOR`, `PCEisenstatSetOmega()`, `PCSORSetIterations()`, `PCSORSetOmega()`
316: @*/
317: PetscErrorCode PCSORSetSymmetric(PC pc, MatSORType flag)
318: {
319: PetscFunctionBegin;
322: PetscTryMethod(pc, "PCSORSetSymmetric_C", (PC, MatSORType), (pc, flag));
323: PetscFunctionReturn(PETSC_SUCCESS);
324: }
326: /*@
327: PCSORSetOmega - Sets the SOR relaxation coefficient, omega
328: (where omega = 1.0 by default).
330: Logically Collective
332: Input Parameters:
333: + pc - the preconditioner context
334: - omega - relaxation coefficient (0 < omega < 2).
336: Options Database Key:
337: . -pc_sor_omega <omega> - Sets omega
339: Level: intermediate
341: Note:
342: If omega != 1, you will need to set the `MAT_USE_INODE`S option to `PETSC_FALSE` on the matrix.
344: .seealso: [](ch_ksp), `PCSOR`, `PCSORSetSymmetric()`, `PCSORSetIterations()`, `PCEisenstatSetOmega()`, `MatSetOption()`
345: @*/
346: PetscErrorCode PCSORSetOmega(PC pc, PetscReal omega)
347: {
348: PetscFunctionBegin;
351: PetscTryMethod(pc, "PCSORSetOmega_C", (PC, PetscReal), (pc, omega));
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: /*@
356: PCSORSetIterations - Sets the number of inner iterations to
357: be used by the SOR preconditioner. The default is 1.
359: Logically Collective
361: Input Parameters:
362: + pc - the preconditioner context
363: . lits - number of local iterations, smoothings over just variables on processor
364: - its - number of parallel iterations to use; each parallel iteration has lits local iterations
366: Options Database Keys:
367: + -pc_sor_its <its> - Sets number of iterations
368: - -pc_sor_lits <lits> - Sets number of local iterations
370: Level: intermediate
372: Note:
373: When run on one processor the number of smoothings is lits*its
375: .seealso: [](ch_ksp), `PCSOR`, `PCSORSetOmega()`, `PCSORSetSymmetric()`
376: @*/
377: PetscErrorCode PCSORSetIterations(PC pc, PetscInt its, PetscInt lits)
378: {
379: PetscFunctionBegin;
382: PetscTryMethod(pc, "PCSORSetIterations_C", (PC, PetscInt, PetscInt), (pc, its, lits));
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: /*MC
387: PCSOR - (S)SOR (successive over relaxation, Gauss-Seidel) preconditioning
389: Options Database Keys:
390: + -pc_sor_symmetric - Activates symmetric version
391: . -pc_sor_backward - Activates backward version
392: . -pc_sor_forward - Activates forward version
393: . -pc_sor_local_forward - Activates local forward version
394: . -pc_sor_local_symmetric - Activates local symmetric version (default version)
395: . -pc_sor_local_backward - Activates local backward version
396: . -pc_sor_omega <omega> - Sets omega
397: . -pc_sor_diagonal_shift <shift> - shift the diagonal entries; useful if the matrix has zeros on the diagonal
398: . -pc_sor_its <its> - Sets number of iterations (default 1)
399: - -pc_sor_lits <lits> - Sets number of local iterations (default 1)
401: Level: beginner
403: Notes:
404: Only implemented for the `MATAIJ` and `MATSEQBAIJ` matrix formats.
406: Not a true parallel SOR, in parallel this implementation corresponds to block
407: Jacobi with SOR on each block.
409: For `MATAIJ` matrices if a diagonal entry is zero (and the diagonal shift is zero) then by default the inverse of that
410: zero will be used and hence the `KSPSolve()` will terminate with `KSP_DIVERGED_NANORINF`. If the option
411: `KSPSetErrorIfNotConverged()` or -ksp_error_if_not_converged the code will terminate as soon as it detects the
412: zero pivot.
414: For `MATSEQBAIJ` matrices this implements point-block SOR, but the omega, its, lits options are not supported.
416: For `MATSEQBAIJ` the diagonal blocks are inverted using dense LU with partial pivoting. If a zero pivot is detected
417: the computation is stopped with an error
419: If used with `KSPRICHARDSON` and no monitors the convergence test is skipped to improve speed, thus it always iterates
420: the maximum number of iterations you've selected for `KSP`. It is usually used in this mode as a smoother for multigrid.
422: If omega != 1, you will need to set the `MAT_USE_INODES` option to `PETSC_FALSE` on the matrix.
424: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCJACOBI`,
425: `PCSORSetIterations()`, `PCSORSetSymmetric()`, `PCSORSetOmega()`, `PCEISENSTAT`, `MatSetOption()`
426: M*/
428: PETSC_EXTERN PetscErrorCode PCCreate_SOR(PC pc)
429: {
430: PC_SOR *jac;
432: PetscFunctionBegin;
433: PetscCall(PetscNew(&jac));
435: pc->ops->apply = PCApply_SOR;
436: pc->ops->applytranspose = PCApplyTranspose_SOR;
437: pc->ops->applyrichardson = PCApplyRichardson_SOR;
438: pc->ops->setfromoptions = PCSetFromOptions_SOR;
439: pc->ops->setup = NULL;
440: pc->ops->view = PCView_SOR;
441: pc->ops->destroy = PCDestroy_SOR;
442: pc->data = (void *)jac;
443: jac->sym = SOR_LOCAL_SYMMETRIC_SWEEP;
444: jac->omega = 1.0;
445: jac->fshift = 0.0;
446: jac->its = 1;
447: jac->lits = 1;
449: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORSetSymmetric_C", PCSORSetSymmetric_SOR));
450: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORSetOmega_C", PCSORSetOmega_SOR));
451: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORSetIterations_C", PCSORSetIterations_SOR));
452: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORGetSymmetric_C", PCSORGetSymmetric_SOR));
453: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORGetOmega_C", PCSORGetOmega_SOR));
454: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORGetIterations_C", PCSORGetIterations_SOR));
455: PetscFunctionReturn(PETSC_SUCCESS);
456: }