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:   PetscObjectComposeFunction((PetscObject)pc, "PCSORSetSymmetric_C", NULL);
 17:   PetscObjectComposeFunction((PetscObject)pc, "PCSORSetOmega_C", NULL);
 18:   PetscObjectComposeFunction((PetscObject)pc, "PCSORSetIterations_C", NULL);
 19:   PetscObjectComposeFunction((PetscObject)pc, "PCSORGetSymmetric_C", NULL);
 20:   PetscObjectComposeFunction((PetscObject)pc, "PCSORGetOmega_C", NULL);
 21:   PetscObjectComposeFunction((PetscObject)pc, "PCSORGetIterations_C", NULL);
 22:   PetscFree(pc->data);
 23:   return 0;
 24: }

 26: static PetscErrorCode PCApply_SOR(PC pc, Vec x, Vec y)
 27: {
 28:   PC_SOR  *jac  = (PC_SOR *)pc->data;
 29:   PetscInt flag = jac->sym | SOR_ZERO_INITIAL_GUESS;

 31:   MatSOR(pc->pmat, x, jac->omega, (MatSORType)flag, jac->fshift, jac->its, jac->lits, y);
 32:   MatFactorGetError(pc->pmat, (MatFactorError *)&pc->failedreason);
 33:   return 0;
 34: }

 36: static PetscErrorCode PCApplyTranspose_SOR(PC pc, Vec x, Vec y)
 37: {
 38:   PC_SOR   *jac  = (PC_SOR *)pc->data;
 39:   PetscInt  flag = jac->sym | SOR_ZERO_INITIAL_GUESS;
 40:   PetscBool set, sym;

 42:   MatIsSymmetricKnown(pc->pmat, &set, &sym);
 44:   MatSOR(pc->pmat, x, jac->omega, (MatSORType)flag, jac->fshift, jac->its, jac->lits, y);
 45:   MatFactorGetError(pc->pmat, (MatFactorError *)&pc->failedreason);
 46:   return 0;
 47: }

 49: 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)
 50: {
 51:   PC_SOR    *jac   = (PC_SOR *)pc->data;
 52:   MatSORType stype = jac->sym;

 54:   PetscInfo(pc, "Warning, convergence criteria ignored, using %" PetscInt_FMT " iterations\n", its);
 55:   if (guesszero) stype = (MatSORType)(stype | SOR_ZERO_INITIAL_GUESS);
 56:   MatSOR(pc->pmat, b, jac->omega, stype, jac->fshift, its * jac->its, jac->lits, y);
 57:   MatFactorGetError(pc->pmat, (MatFactorError *)&pc->failedreason);
 58:   *outits = its;
 59:   *reason = PCRICHARDSON_CONVERGED_ITS;
 60:   return 0;
 61: }

 63: PetscErrorCode PCSetFromOptions_SOR(PC pc, PetscOptionItems *PetscOptionsObject)
 64: {
 65:   PC_SOR   *jac = (PC_SOR *)pc->data;
 66:   PetscBool flg;
 67:   PetscReal omega;

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

 91: PetscErrorCode PCView_SOR(PC pc, PetscViewer viewer)
 92: {
 93:   PC_SOR     *jac = (PC_SOR *)pc->data;
 94:   MatSORType  sym = jac->sym;
 95:   const char *sortype;
 96:   PetscBool   iascii;

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

116: static PetscErrorCode PCSORSetSymmetric_SOR(PC pc, MatSORType flag)
117: {
118:   PC_SOR *jac = (PC_SOR *)pc->data;

120:   jac->sym = flag;
121:   return 0;
122: }

124: static PetscErrorCode PCSORSetOmega_SOR(PC pc, PetscReal omega)
125: {
126:   PC_SOR *jac = (PC_SOR *)pc->data;

129:   jac->omega = omega;
130:   return 0;
131: }

133: static PetscErrorCode PCSORSetIterations_SOR(PC pc, PetscInt its, PetscInt lits)
134: {
135:   PC_SOR *jac = (PC_SOR *)pc->data;

137:   jac->its  = its;
138:   jac->lits = lits;
139:   return 0;
140: }

142: static PetscErrorCode PCSORGetSymmetric_SOR(PC pc, MatSORType *flag)
143: {
144:   PC_SOR *jac = (PC_SOR *)pc->data;

146:   *flag = jac->sym;
147:   return 0;
148: }

150: static PetscErrorCode PCSORGetOmega_SOR(PC pc, PetscReal *omega)
151: {
152:   PC_SOR *jac = (PC_SOR *)pc->data;

154:   *omega = jac->omega;
155:   return 0;
156: }

158: static PetscErrorCode PCSORGetIterations_SOR(PC pc, PetscInt *its, PetscInt *lits)
159: {
160:   PC_SOR *jac = (PC_SOR *)pc->data;

162:   if (its) *its = jac->its;
163:   if (lits) *lits = jac->lits;
164:   return 0;
165: }

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

171:    Logically Collective

173:    Input Parameter:
174: .  pc - the preconditioner context

176:    Output Parameter:
177: .  flag - one of the following
178: .vb
179:     SOR_FORWARD_SWEEP
180:     SOR_BACKWARD_SWEEP
181:     SOR_SYMMETRIC_SWEEP
182:     SOR_LOCAL_FORWARD_SWEEP
183:     SOR_LOCAL_BACKWARD_SWEEP
184:     SOR_LOCAL_SYMMETRIC_SWEEP
185: .ve

187:    Options Database Keys:
188: +  -pc_sor_symmetric - Activates symmetric version
189: .  -pc_sor_backward - Activates backward version
190: .  -pc_sor_local_forward - Activates local forward version
191: .  -pc_sor_local_symmetric - Activates local symmetric version
192: -  -pc_sor_local_backward - Activates local backward version

194:    Note:
195:    To use the Eisenstat trick with SSOR, employ the `PCEISENSTAT` preconditioner,
196:    which can be chosen with the option
197: .  -pc_type eisenstat - Activates Eisenstat trick

199:    Level: intermediate

201: .seealso: `PCSOR`, `PCEisenstatSetOmega()`, `PCSORSetIterations()`, `PCSORSetOmega()`, `PCSORSetSymmetric()`
202: @*/
203: PetscErrorCode PCSORGetSymmetric(PC pc, MatSORType *flag)
204: {
206:   PetscUseMethod(pc, "PCSORGetSymmetric_C", (PC, MatSORType *), (pc, flag));
207:   return 0;
208: }

210: /*@
211:    PCSORGetOmega - Gets the SOR relaxation coefficient, omega
212:    (where omega = 1.0 by default).

214:    Logically Collective

216:    Input Parameter:
217: .  pc - the preconditioner context

219:    Output Parameter:
220: .  omega - relaxation coefficient (0 < omega < 2).

222:    Options Database Key:
223: .  -pc_sor_omega <omega> - Sets omega

225:    Level: intermediate

227: .seealso: `PCSOR`, `PCSORSetSymmetric()`, `PCSORSetIterations()`, `PCEisenstatSetOmega()`, `PCSORSetOmega()`
228: @*/
229: PetscErrorCode PCSORGetOmega(PC pc, PetscReal *omega)
230: {
232:   PetscUseMethod(pc, "PCSORGetOmega_C", (PC, PetscReal *), (pc, omega));
233:   return 0;
234: }

236: /*@
237:    PCSORGetIterations - Gets the number of inner iterations to
238:    be used by the SOR preconditioner. The default is 1.

240:    Logically Collective

242:    Input Parameter:
243: .  pc - the preconditioner context

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

249:    Options Database Keys:
250: +  -pc_sor_its <its> - Sets number of iterations
251: -  -pc_sor_lits <lits> - Sets number of local iterations

253:    Level: intermediate

255:    Note:
256:     When run on one processor the number of smoothings is lits*its

258: .seealso: `PCSOR`, `PCSORSetOmega()`, `PCSORSetSymmetric()`, `PCSORSetIterations()`
259: @*/
260: PetscErrorCode PCSORGetIterations(PC pc, PetscInt *its, PetscInt *lits)
261: {
263:   PetscUseMethod(pc, "PCSORGetIterations_C", (PC, PetscInt *, PetscInt *), (pc, its, lits));
264:   return 0;
265: }

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

272:    Logically Collective

274:    Input Parameters:
275: +  pc - the preconditioner context
276: -  flag - one of the following
277: .vb
278:     SOR_FORWARD_SWEEP
279:     SOR_BACKWARD_SWEEP
280:     SOR_SYMMETRIC_SWEEP
281:     SOR_LOCAL_FORWARD_SWEEP
282:     SOR_LOCAL_BACKWARD_SWEEP
283:     SOR_LOCAL_SYMMETRIC_SWEEP
284: .ve

286:    Options Database Keys:
287: +  -pc_sor_symmetric - Activates symmetric version
288: .  -pc_sor_backward - Activates backward version
289: .  -pc_sor_local_forward - Activates local forward version
290: .  -pc_sor_local_symmetric - Activates local symmetric version
291: -  -pc_sor_local_backward - Activates local backward version

293:    Note:
294:    To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner,
295:    which can be chosen with the option
296: .  -pc_type eisenstat - Activates Eisenstat trick

298:    Level: intermediate

300: .seealso: `PCSOR`, `PCEisenstatSetOmega()`, `PCSORSetIterations()`, `PCSORSetOmega()`
301: @*/
302: PetscErrorCode PCSORSetSymmetric(PC pc, MatSORType flag)
303: {
306:   PetscTryMethod(pc, "PCSORSetSymmetric_C", (PC, MatSORType), (pc, flag));
307:   return 0;
308: }

310: /*@
311:    PCSORSetOmega - Sets the SOR relaxation coefficient, omega
312:    (where omega = 1.0 by default).

314:    Logically Collective

316:    Input Parameters:
317: +  pc - the preconditioner context
318: -  omega - relaxation coefficient (0 < omega < 2).

320:    Options Database Key:
321: .  -pc_sor_omega <omega> - Sets omega

323:    Level: intermediate

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

328: .seealso: `PCSOR`, `PCSORSetSymmetric()`, `PCSORSetIterations()`, `PCEisenstatSetOmega()`, `MatSetOption()`
329: @*/
330: PetscErrorCode PCSORSetOmega(PC pc, PetscReal omega)
331: {
334:   PetscTryMethod(pc, "PCSORSetOmega_C", (PC, PetscReal), (pc, omega));
335:   return 0;
336: }

338: /*@
339:    PCSORSetIterations - Sets the number of inner iterations to
340:    be used by the SOR preconditioner. The default is 1.

342:    Logically Collective

344:    Input Parameters:
345: +  pc - the preconditioner context
346: .  lits - number of local iterations, smoothings over just variables on processor
347: -  its - number of parallel iterations to use; each parallel iteration has lits local iterations

349:    Options Database Keys:
350: +  -pc_sor_its <its> - Sets number of iterations
351: -  -pc_sor_lits <lits> - Sets number of local iterations

353:    Level: intermediate

355:    Note:
356:     When run on one processor the number of smoothings is lits*its

358: .seealso: `PCSOR`, `PCSORSetOmega()`, `PCSORSetSymmetric()`
359: @*/
360: PetscErrorCode PCSORSetIterations(PC pc, PetscInt its, PetscInt lits)
361: {
364:   PetscTryMethod(pc, "PCSORSetIterations_C", (PC, PetscInt, PetscInt), (pc, its, lits));
365:   return 0;
366: }

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

371:    Options Database Keys:
372: +  -pc_sor_symmetric - Activates symmetric version
373: .  -pc_sor_backward - Activates backward version
374: .  -pc_sor_forward - Activates forward version
375: .  -pc_sor_local_forward - Activates local forward version
376: .  -pc_sor_local_symmetric - Activates local symmetric version  (default version)
377: .  -pc_sor_local_backward - Activates local backward version
378: .  -pc_sor_omega <omega> - Sets omega
379: .  -pc_sor_diagonal_shift <shift> - shift the diagonal entries; useful if the matrix has zeros on the diagonal
380: .  -pc_sor_its <its> - Sets number of iterations   (default 1)
381: -  -pc_sor_lits <lits> - Sets number of local iterations  (default 1)

383:    Level: beginner

385:    Notes:
386:    Only implemented for the `MATAIJ`  and `MATSEQBAIJ` matrix formats.

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

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

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

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

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

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

406: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCJACOBI`,
407:           `PCSORSetIterations()`, `PCSORSetSymmetric()`, `PCSORSetOmega()`, `PCEISENSTAT`, `MatSetOption()`
408: M*/

410: PETSC_EXTERN PetscErrorCode PCCreate_SOR(PC pc)
411: {
412:   PC_SOR *jac;

414:   PetscNew(&jac);

416:   pc->ops->apply           = PCApply_SOR;
417:   pc->ops->applytranspose  = PCApplyTranspose_SOR;
418:   pc->ops->applyrichardson = PCApplyRichardson_SOR;
419:   pc->ops->setfromoptions  = PCSetFromOptions_SOR;
420:   pc->ops->setup           = NULL;
421:   pc->ops->view            = PCView_SOR;
422:   pc->ops->destroy         = PCDestroy_SOR;
423:   pc->data                 = (void *)jac;
424:   jac->sym                 = SOR_LOCAL_SYMMETRIC_SWEEP;
425:   jac->omega               = 1.0;
426:   jac->fshift              = 0.0;
427:   jac->its                 = 1;
428:   jac->lits                = 1;

430:   PetscObjectComposeFunction((PetscObject)pc, "PCSORSetSymmetric_C", PCSORSetSymmetric_SOR);
431:   PetscObjectComposeFunction((PetscObject)pc, "PCSORSetOmega_C", PCSORSetOmega_SOR);
432:   PetscObjectComposeFunction((PetscObject)pc, "PCSORSetIterations_C", PCSORSetIterations_SOR);
433:   PetscObjectComposeFunction((PetscObject)pc, "PCSORGetSymmetric_C", PCSORGetSymmetric_SOR);
434:   PetscObjectComposeFunction((PetscObject)pc, "PCSORGetOmega_C", PCSORGetOmega_SOR);
435:   PetscObjectComposeFunction((PetscObject)pc, "PCSORGetIterations_C", PCSORGetIterations_SOR);
436:   return 0;
437: }