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   isascii;

    102:   PetscFunctionBegin;
    103:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    104:   if (isascii) {
    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: }