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:   if (guesszero) stype = (MatSORType)(stype | SOR_ZERO_INITIAL_GUESS);
 59:   PetscCall(MatSOR(pc->pmat, b, jac->omega, stype, jac->fshift, its * jac->its, jac->lits, y));
 60:   PetscCall(MatFactorGetError(pc->pmat, (MatFactorError *)&pc->failedreason));
 61:   *outits = its;
 62:   *reason = PCRICHARDSON_CONVERGED_ITS;
 63:   PetscFunctionReturn(PETSC_SUCCESS);
 64: }

 66: static PetscErrorCode PCSetFromOptions_SOR(PC pc, PetscOptionItems *PetscOptionsObject)
 67: {
 68:   PC_SOR   *jac = (PC_SOR *)pc->data;
 69:   PetscBool flg;
 70:   PetscReal omega;

 72:   PetscFunctionBegin;
 73:   PetscOptionsHeadBegin(PetscOptionsObject, "(S)SOR options");
 74:   PetscCall(PetscOptionsReal("-pc_sor_omega", "relaxation factor (0 < omega < 2)", "PCSORSetOmega", jac->omega, &omega, &flg));
 75:   if (flg) PetscCall(PCSORSetOmega(pc, omega));
 76:   PetscCall(PetscOptionsReal("-pc_sor_diagonal_shift", "Add to the diagonal entries", "", jac->fshift, &jac->fshift, NULL));
 77:   PetscCall(PetscOptionsInt("-pc_sor_its", "number of inner SOR iterations", "PCSORSetIterations", jac->its, &jac->its, NULL));
 78:   PetscCall(PetscOptionsInt("-pc_sor_lits", "number of local inner SOR iterations", "PCSORSetIterations", jac->lits, &jac->lits, NULL));
 79:   PetscCall(PetscOptionsBoolGroupBegin("-pc_sor_symmetric", "SSOR, not SOR", "PCSORSetSymmetric", &flg));
 80:   if (flg) PetscCall(PCSORSetSymmetric(pc, SOR_SYMMETRIC_SWEEP));
 81:   PetscCall(PetscOptionsBoolGroup("-pc_sor_backward", "use backward sweep instead of forward", "PCSORSetSymmetric", &flg));
 82:   if (flg) PetscCall(PCSORSetSymmetric(pc, SOR_BACKWARD_SWEEP));
 83:   PetscCall(PetscOptionsBoolGroup("-pc_sor_forward", "use forward sweep", "PCSORSetSymmetric", &flg));
 84:   if (flg) PetscCall(PCSORSetSymmetric(pc, SOR_FORWARD_SWEEP));
 85:   PetscCall(PetscOptionsBoolGroup("-pc_sor_local_symmetric", "use SSOR separately on each processor", "PCSORSetSymmetric", &flg));
 86:   if (flg) PetscCall(PCSORSetSymmetric(pc, SOR_LOCAL_SYMMETRIC_SWEEP));
 87:   PetscCall(PetscOptionsBoolGroup("-pc_sor_local_backward", "use backward sweep locally", "PCSORSetSymmetric", &flg));
 88:   if (flg) PetscCall(PCSORSetSymmetric(pc, SOR_LOCAL_BACKWARD_SWEEP));
 89:   PetscCall(PetscOptionsBoolGroupEnd("-pc_sor_local_forward", "use forward sweep locally", "PCSORSetSymmetric", &flg));
 90:   if (flg) PetscCall(PCSORSetSymmetric(pc, SOR_LOCAL_FORWARD_SWEEP));
 91:   PetscOptionsHeadEnd();
 92:   PetscFunctionReturn(PETSC_SUCCESS);
 93: }

 95: static PetscErrorCode PCView_SOR(PC pc, PetscViewer viewer)
 96: {
 97:   PC_SOR     *jac = (PC_SOR *)pc->data;
 98:   MatSORType  sym = jac->sym;
 99:   const char *sortype;
100:   PetscBool   iascii;

102:   PetscFunctionBegin;
103:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
104:   if (iascii) {
105:     if (sym & SOR_ZERO_INITIAL_GUESS) PetscCall(PetscViewerASCIIPrintf(viewer, "  zero initial guess\n"));
106:     if (sym == SOR_APPLY_UPPER) sortype = "apply_upper";
107:     else if (sym == SOR_APPLY_LOWER) sortype = "apply_lower";
108:     else if (sym & SOR_EISENSTAT) sortype = "Eisenstat";
109:     else if ((sym & SOR_SYMMETRIC_SWEEP) == SOR_SYMMETRIC_SWEEP) sortype = "symmetric";
110:     else if (sym & SOR_BACKWARD_SWEEP) sortype = "backward";
111:     else if (sym & SOR_FORWARD_SWEEP) sortype = "forward";
112:     else if ((sym & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) sortype = "local_symmetric";
113:     else if (sym & SOR_LOCAL_FORWARD_SWEEP) sortype = "local_forward";
114:     else if (sym & SOR_LOCAL_BACKWARD_SWEEP) sortype = "local_backward";
115:     else sortype = "unknown";
116:     PetscCall(PetscViewerASCIIPrintf(viewer, "  type = %s, iterations = %" PetscInt_FMT ", local iterations = %" PetscInt_FMT ", omega = %g\n", sortype, jac->its, jac->lits, (double)jac->omega));
117:   }
118:   PetscFunctionReturn(PETSC_SUCCESS);
119: }

121: static PetscErrorCode PCSORSetSymmetric_SOR(PC pc, MatSORType flag)
122: {
123:   PC_SOR *jac = (PC_SOR *)pc->data;

125:   PetscFunctionBegin;
126:   jac->sym = flag;
127:   PetscFunctionReturn(PETSC_SUCCESS);
128: }

130: static PetscErrorCode PCSORSetOmega_SOR(PC pc, PetscReal omega)
131: {
132:   PC_SOR *jac = (PC_SOR *)pc->data;

134:   PetscFunctionBegin;
135:   PetscCheck(omega > 0.0 && omega < 2.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Relaxation out of range");
136:   jac->omega = omega;
137:   PetscFunctionReturn(PETSC_SUCCESS);
138: }

140: static PetscErrorCode PCSORSetIterations_SOR(PC pc, PetscInt its, PetscInt lits)
141: {
142:   PC_SOR *jac = (PC_SOR *)pc->data;

144:   PetscFunctionBegin;
145:   jac->its  = its;
146:   jac->lits = lits;
147:   PetscFunctionReturn(PETSC_SUCCESS);
148: }

150: static PetscErrorCode PCSORGetSymmetric_SOR(PC pc, MatSORType *flag)
151: {
152:   PC_SOR *jac = (PC_SOR *)pc->data;

154:   PetscFunctionBegin;
155:   *flag = jac->sym;
156:   PetscFunctionReturn(PETSC_SUCCESS);
157: }

159: static PetscErrorCode PCSORGetOmega_SOR(PC pc, PetscReal *omega)
160: {
161:   PC_SOR *jac = (PC_SOR *)pc->data;

163:   PetscFunctionBegin;
164:   *omega = jac->omega;
165:   PetscFunctionReturn(PETSC_SUCCESS);
166: }

168: static PetscErrorCode PCSORGetIterations_SOR(PC pc, PetscInt *its, PetscInt *lits)
169: {
170:   PC_SOR *jac = (PC_SOR *)pc->data;

172:   PetscFunctionBegin;
173:   if (its) *its = jac->its;
174:   if (lits) *lits = jac->lits;
175:   PetscFunctionReturn(PETSC_SUCCESS);
176: }

178: /*@
179:   PCSORGetSymmetric - Gets the form the SOR preconditioner is using;   backward, or forward relaxation.  The local variants perform SOR on
180:   each processor.  By default forward relaxation is used.

182:   Logically Collective

184:   Input Parameter:
185: . pc - the preconditioner context

187:   Output Parameter:
188: . flag - one of the following
189: .vb
190:     SOR_FORWARD_SWEEP
191:     SOR_BACKWARD_SWEEP
192:     SOR_SYMMETRIC_SWEEP
193:     SOR_LOCAL_FORWARD_SWEEP
194:     SOR_LOCAL_BACKWARD_SWEEP
195:     SOR_LOCAL_SYMMETRIC_SWEEP
196: .ve

198:   Options Database Keys:
199: + -pc_sor_symmetric       - Activates symmetric version
200: . -pc_sor_backward        - Activates backward version
201: . -pc_sor_local_forward   - Activates local forward version
202: . -pc_sor_local_symmetric - Activates local symmetric version
203: - -pc_sor_local_backward  - Activates local backward version

205:   Note:
206:   To use the Eisenstat trick with SSOR, employ the `PCEISENSTAT` preconditioner,
207:   which can be chosen with the option
208: .  -pc_type eisenstat - Activates Eisenstat trick

210:   Level: intermediate

212: .seealso: [](ch_ksp), `PCSOR`, `PCEisenstatSetOmega()`, `PCSORSetIterations()`, `PCSORSetOmega()`, `PCSORSetSymmetric()`
213: @*/
214: PetscErrorCode PCSORGetSymmetric(PC pc, MatSORType *flag)
215: {
216:   PetscFunctionBegin;
218:   PetscUseMethod(pc, "PCSORGetSymmetric_C", (PC, MatSORType *), (pc, flag));
219:   PetscFunctionReturn(PETSC_SUCCESS);
220: }

222: /*@
223:   PCSORGetOmega - Gets the SOR relaxation coefficient, omega
224:   (where omega = 1.0 by default).

226:   Logically Collective

228:   Input Parameter:
229: . pc - the preconditioner context

231:   Output Parameter:
232: . omega - relaxation coefficient (0 < omega < 2).

234:   Options Database Key:
235: . -pc_sor_omega <omega> - Sets omega

237:   Level: intermediate

239: .seealso: [](ch_ksp), `PCSOR`, `PCSORSetSymmetric()`, `PCSORSetIterations()`, `PCEisenstatSetOmega()`, `PCSORSetOmega()`
240: @*/
241: PetscErrorCode PCSORGetOmega(PC pc, PetscReal *omega)
242: {
243:   PetscFunctionBegin;
245:   PetscUseMethod(pc, "PCSORGetOmega_C", (PC, PetscReal *), (pc, omega));
246:   PetscFunctionReturn(PETSC_SUCCESS);
247: }

249: /*@
250:   PCSORGetIterations - Gets the number of inner iterations to
251:   be used by the SOR preconditioner. The default is 1.

253:   Logically Collective

255:   Input Parameter:
256: . pc - the preconditioner context

258:   Output Parameters:
259: + lits - number of local iterations, smoothings over just variables on processor
260: - its  - number of parallel iterations to use; each parallel iteration has lits local iterations

262:   Options Database Keys:
263: + -pc_sor_its <its>   - Sets number of iterations
264: - -pc_sor_lits <lits> - Sets number of local iterations

266:   Level: intermediate

268:   Note:
269:   When run on one processor the number of smoothings is lits*its

271: .seealso: [](ch_ksp), `PCSOR`, `PCSORSetOmega()`, `PCSORSetSymmetric()`, `PCSORSetIterations()`
272: @*/
273: PetscErrorCode PCSORGetIterations(PC pc, PetscInt *its, PetscInt *lits)
274: {
275:   PetscFunctionBegin;
277:   PetscUseMethod(pc, "PCSORGetIterations_C", (PC, PetscInt *, PetscInt *), (pc, its, lits));
278:   PetscFunctionReturn(PETSC_SUCCESS);
279: }

281: /*@
282:   PCSORSetSymmetric - Sets the SOR preconditioner to use symmetric (SSOR),
283:   backward, or forward relaxation.  The local variants perform SOR on
284:   each processor.  By default forward relaxation is used.

286:   Logically Collective

288:   Input Parameters:
289: + pc   - the preconditioner context
290: - flag - one of the following
291: .vb
292:     SOR_FORWARD_SWEEP
293:     SOR_BACKWARD_SWEEP
294:     SOR_SYMMETRIC_SWEEP
295:     SOR_LOCAL_FORWARD_SWEEP
296:     SOR_LOCAL_BACKWARD_SWEEP
297:     SOR_LOCAL_SYMMETRIC_SWEEP
298: .ve

300:   Options Database Keys:
301: + -pc_sor_symmetric       - Activates symmetric version
302: . -pc_sor_backward        - Activates backward version
303: . -pc_sor_local_forward   - Activates local forward version
304: . -pc_sor_local_symmetric - Activates local symmetric version
305: - -pc_sor_local_backward  - Activates local backward version

307:   Note:
308:   To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner,
309:   which can be chosen with the option
310: .  -pc_type eisenstat - Activates Eisenstat trick

312:   Level: intermediate

314: .seealso: [](ch_ksp), `PCSOR`, `PCEisenstatSetOmega()`, `PCSORSetIterations()`, `PCSORSetOmega()`
315: @*/
316: PetscErrorCode PCSORSetSymmetric(PC pc, MatSORType flag)
317: {
318:   PetscFunctionBegin;
321:   PetscTryMethod(pc, "PCSORSetSymmetric_C", (PC, MatSORType), (pc, flag));
322:   PetscFunctionReturn(PETSC_SUCCESS);
323: }

325: /*@
326:   PCSORSetOmega - Sets the SOR relaxation coefficient, omega
327:   (where omega = 1.0 by default).

329:   Logically Collective

331:   Input Parameters:
332: + pc    - the preconditioner context
333: - omega - relaxation coefficient (0 < omega < 2).

335:   Options Database Key:
336: . -pc_sor_omega <omega> - Sets omega

338:   Level: intermediate

340:   Note:
341:   If omega != 1, you will need to set the `MAT_USE_INODE`S option to `PETSC_FALSE` on the matrix.

343: .seealso: [](ch_ksp), `PCSOR`, `PCSORSetSymmetric()`, `PCSORSetIterations()`, `PCEisenstatSetOmega()`, `MatSetOption()`
344: @*/
345: PetscErrorCode PCSORSetOmega(PC pc, PetscReal omega)
346: {
347:   PetscFunctionBegin;
350:   PetscTryMethod(pc, "PCSORSetOmega_C", (PC, PetscReal), (pc, omega));
351:   PetscFunctionReturn(PETSC_SUCCESS);
352: }

354: /*@
355:   PCSORSetIterations - Sets the number of inner iterations to
356:   be used by the SOR preconditioner. The default is 1.

358:   Logically Collective

360:   Input Parameters:
361: + pc   - the preconditioner context
362: . lits - number of local iterations, smoothings over just variables on processor
363: - its  - number of parallel iterations to use; each parallel iteration has lits local iterations

365:   Options Database Keys:
366: + -pc_sor_its <its>   - Sets number of iterations
367: - -pc_sor_lits <lits> - Sets number of local iterations

369:   Level: intermediate

371:   Note:
372:   When run on one processor the number of smoothings is lits*its

374: .seealso: [](ch_ksp), `PCSOR`, `PCSORSetOmega()`, `PCSORSetSymmetric()`
375: @*/
376: PetscErrorCode PCSORSetIterations(PC pc, PetscInt its, PetscInt lits)
377: {
378:   PetscFunctionBegin;
381:   PetscTryMethod(pc, "PCSORSetIterations_C", (PC, PetscInt, PetscInt), (pc, its, lits));
382:   PetscFunctionReturn(PETSC_SUCCESS);
383: }

385: /*MC
386:      PCSOR - (S)SOR (successive over relaxation, Gauss-Seidel) preconditioning

388:    Options Database Keys:
389: +  -pc_sor_symmetric - Activates symmetric version
390: .  -pc_sor_backward - Activates backward version
391: .  -pc_sor_forward - Activates forward version
392: .  -pc_sor_local_forward - Activates local forward version
393: .  -pc_sor_local_symmetric - Activates local symmetric version  (default version)
394: .  -pc_sor_local_backward - Activates local backward version
395: .  -pc_sor_omega <omega> - Sets omega
396: .  -pc_sor_diagonal_shift <shift> - shift the diagonal entries; useful if the matrix has zeros on the diagonal
397: .  -pc_sor_its <its> - Sets number of iterations   (default 1)
398: -  -pc_sor_lits <lits> - Sets number of local iterations  (default 1)

400:    Level: beginner

402:    Notes:
403:    Only implemented for the `MATAIJ`  and `MATSEQBAIJ` matrix formats.

405:    Not a true parallel SOR, in parallel this implementation corresponds to block
406:    Jacobi with SOR on each block.

408:           For `MATAIJ` matrices if a diagonal entry is zero (and the diagonal shift is zero) then by default the inverse of that
409:           zero will be used and hence the `KSPSolve()` will terminate with `KSP_DIVERGED_NANORINF`. If the option
410:           `KSPSetErrorIfNotConverged()` or -ksp_error_if_not_converged the code will terminate as soon as it detects the
411:           zero pivot.

413:           For `MATSEQBAIJ` matrices this implements point-block SOR, but the omega, its, lits options are not supported.

415:           For `MATSEQBAIJ` the diagonal blocks are inverted using dense LU with partial pivoting. If a zero pivot is detected
416:           the computation is stopped with an error

418:           If used with `KSPRICHARDSON` and no monitors the convergence test is skipped to improve speed, thus it always iterates
419:           the maximum number of iterations you've selected for `KSP`. It is usually used in this mode as a smoother for multigrid.

421:           If omega != 1, you will need to set the `MAT_USE_INODES` option to `PETSC_FALSE` on the matrix.

423: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCJACOBI`,
424:           `PCSORSetIterations()`, `PCSORSetSymmetric()`, `PCSORSetOmega()`, `PCEISENSTAT`, `MatSetOption()`
425: M*/

427: PETSC_EXTERN PetscErrorCode PCCreate_SOR(PC pc)
428: {
429:   PC_SOR *jac;

431:   PetscFunctionBegin;
432:   PetscCall(PetscNew(&jac));

434:   pc->ops->apply           = PCApply_SOR;
435:   pc->ops->applytranspose  = PCApplyTranspose_SOR;
436:   pc->ops->applyrichardson = PCApplyRichardson_SOR;
437:   pc->ops->setfromoptions  = PCSetFromOptions_SOR;
438:   pc->ops->setup           = NULL;
439:   pc->ops->view            = PCView_SOR;
440:   pc->ops->destroy         = PCDestroy_SOR;
441:   pc->data                 = (void *)jac;
442:   jac->sym                 = SOR_LOCAL_SYMMETRIC_SWEEP;
443:   jac->omega               = 1.0;
444:   jac->fshift              = 0.0;
445:   jac->its                 = 1;
446:   jac->lits                = 1;

448:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORSetSymmetric_C", PCSORSetSymmetric_SOR));
449:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORSetOmega_C", PCSORSetOmega_SOR));
450:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORSetIterations_C", PCSORSetIterations_SOR));
451:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORGetSymmetric_C", PCSORGetSymmetric_SOR));
452:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORGetOmega_C", PCSORGetOmega_SOR));
453:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORGetIterations_C", PCSORGetIterations_SOR));
454:   PetscFunctionReturn(PETSC_SUCCESS);
455: }