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: }