Actual source code: jacobi.c

      1: /*
      2:      This file implements a Jacobi preconditioner in PETSc as part of PC.
      3:      You can use this as a starting point for implementing your own
      4:      preconditioner that is not provided with PETSc. (You might also consider
      5:      just using PCSHELL)

      7:      The following basic routines are required for each preconditioner.
      8:           PCCreate_XXX()          - Creates a preconditioner context
      9:           PCSetFromOptions_XXX()  - Sets runtime options
     10:           PCApply_XXX()           - Applies the preconditioner
     11:           PCDestroy_XXX()         - Destroys the preconditioner context
     12:      where the suffix "_XXX" denotes a particular implementation, in
     13:      this case we use _Jacobi (e.g., PCCreate_Jacobi, PCApply_Jacobi).
     14:      These routines are actually called via the common user interface
     15:      routines PCCreate(), PCSetFromOptions(), PCApply(), and PCDestroy(),
     16:      so the application code interface remains identical for all
     17:      preconditioners.

     19:      Another key routine is:
     20:           PCSetUp_XXX()           - Prepares for the use of a preconditioner
     21:      by setting data structures and options.   The interface routine PCSetUp()
     22:      is not usually called directly by the user, but instead is called by
     23:      PCApply() if necessary.

     25:      Additional basic routines are:
     26:           PCView_XXX()            - Prints details of runtime options that
     27:                                     have actually been used.
     28:      These are called by application codes via the interface routines
     29:      PCView().

     31:      The various types of solvers (preconditioners, Krylov subspace methods,
     32:      nonlinear solvers, timesteppers) are all organized similarly, so the
     33:      above description applies to these categories also.  One exception is
     34:      that the analogues of PCApply() for these components are KSPSolve(),
     35:      SNESSolve(), and TSSolve().

     37:      Additional optional functionality unique to preconditioners is left and
     38:      right symmetric preconditioner application via PCApplySymmetricLeft()
     39:      and PCApplySymmetricRight().  The Jacobi implementation is
     40:      PCApplySymmetricLeftOrRight_Jacobi().
     41: */

     43: /*
     44:    Include files needed for the Jacobi preconditioner:
     45:      pcimpl.h - private include file intended for use by all preconditioners
     46: */

     48: #include <petsc/private/pcimpl.h>

     50: const char *const PCJacobiTypes[] = {"DIAGONAL", "ROWL1", "ROWMAX", "ROWSUM", "PCJacobiType", "PC_JACOBI_", NULL};

     52: /*
     53:    Private context (data structure) for the Jacobi preconditioner.
     54: */
     55: typedef struct {
     56:   Vec          diag;     /* vector containing the reciprocals of the diagonal elements of the matrix used to construct the preconditioner */
     57:   Vec          diagsqrt; /* vector containing the reciprocals of the square roots of
     58:                                     the diagonal elements of the matrix used to compute the preconditioner (used
     59:                                     only for symmetric preconditioner application) */
     60:   PCJacobiType type;
     61:   PetscBool    useabs;  /* use the absolute values of the diagonal entries */
     62:   PetscBool    fixdiag; /* fix zero diagonal terms */
     63:   PetscReal    scale;   /* for scaling rowl1 off-diagonals */
     64: } PC_Jacobi;

     66: static PetscErrorCode PCReset_Jacobi(PC);

     68: static PetscErrorCode PCJacobiSetType_Jacobi(PC pc, PCJacobiType type)
     69: {
     70:   PC_Jacobi   *j = (PC_Jacobi *)pc->data;
     71:   PCJacobiType old_type;

     73:   PetscFunctionBegin;
     74:   PetscCall(PCJacobiGetType(pc, &old_type));
     75:   if (old_type == type) PetscFunctionReturn(PETSC_SUCCESS);
     76:   PetscCall(PCReset_Jacobi(pc));
     77:   j->type = type;
     78:   PetscFunctionReturn(PETSC_SUCCESS);
     79: }

     81: static PetscErrorCode PCJacobiGetUseAbs_Jacobi(PC pc, PetscBool *flg)
     82: {
     83:   PC_Jacobi *j = (PC_Jacobi *)pc->data;

     85:   PetscFunctionBegin;
     86:   *flg = j->useabs;
     87:   PetscFunctionReturn(PETSC_SUCCESS);
     88: }

     90: static PetscErrorCode PCJacobiSetUseAbs_Jacobi(PC pc, PetscBool flg)
     91: {
     92:   PC_Jacobi *j = (PC_Jacobi *)pc->data;

     94:   PetscFunctionBegin;
     95:   j->useabs = flg;
     96:   PetscFunctionReturn(PETSC_SUCCESS);
     97: }

     99: static PetscErrorCode PCJacobiGetType_Jacobi(PC pc, PCJacobiType *type)
    100: {
    101:   PC_Jacobi *j = (PC_Jacobi *)pc->data;

    103:   PetscFunctionBegin;
    104:   *type = j->type;
    105:   PetscFunctionReturn(PETSC_SUCCESS);
    106: }

    108: static PetscErrorCode PCJacobiSetRowl1Scale_Jacobi(PC pc, PetscReal flg)
    109: {
    110:   PC_Jacobi *j = (PC_Jacobi *)pc->data;

    112:   PetscFunctionBegin;
    113:   j->scale = flg;
    114:   PetscFunctionReturn(PETSC_SUCCESS);
    115: }

    117: static PetscErrorCode PCJacobiGetRowl1Scale_Jacobi(PC pc, PetscReal *flg)
    118: {
    119:   PC_Jacobi *j = (PC_Jacobi *)pc->data;

    121:   PetscFunctionBegin;
    122:   *flg = j->scale;
    123:   PetscFunctionReturn(PETSC_SUCCESS);
    124: }

    126: static PetscErrorCode PCJacobiSetFixDiagonal_Jacobi(PC pc, PetscBool flg)
    127: {
    128:   PC_Jacobi *j = (PC_Jacobi *)pc->data;

    130:   PetscFunctionBegin;
    131:   j->fixdiag = flg;
    132:   PetscFunctionReturn(PETSC_SUCCESS);
    133: }

    135: static PetscErrorCode PCJacobiGetFixDiagonal_Jacobi(PC pc, PetscBool *flg)
    136: {
    137:   PC_Jacobi *j = (PC_Jacobi *)pc->data;

    139:   PetscFunctionBegin;
    140:   *flg = j->fixdiag;
    141:   PetscFunctionReturn(PETSC_SUCCESS);
    142: }

    144: static PetscErrorCode PCJacobiGetDiagonal_Jacobi(PC pc, Vec diag, Vec diagsqrt)
    145: {
    146:   PC_Jacobi *j    = (PC_Jacobi *)pc->data;
    147:   MPI_Comm   comm = PetscObjectComm((PetscObject)pc);

    149:   PetscFunctionBegin;
    150:   PetscCheck(j->diag || j->diagsqrt, comm, PETSC_ERR_ARG_WRONGSTATE, "Jacobi diagonal has not been created yet. Use PCApply to force creation");
    151:   PetscCheck(!diag || (diag && j->diag), comm, PETSC_ERR_ARG_WRONGSTATE, "Jacobi diagonal not available. Check if PC is non-symmetric");
    152:   PetscCheck(!diagsqrt || (diagsqrt && j->diagsqrt), comm, PETSC_ERR_ARG_WRONGSTATE, "Jacobi diagonal squareroot not available. Check if PC is symmetric");

    154:   if (diag) PetscCall(VecCopy(j->diag, diag));
    155:   if (diagsqrt) PetscCall(VecCopy(j->diagsqrt, diagsqrt));
    156:   PetscFunctionReturn(PETSC_SUCCESS);
    157: }

    159: /*
    160:    PCSetUp_Jacobi - Prepares for the use of the Jacobi preconditioner
    161:                     by setting data structures and options.

    163:    Input Parameter:
    164: .  pc - the preconditioner context

    166:    Application Interface Routine: PCSetUp()

    168:    Note:
    169:    The interface routine PCSetUp() is not usually called directly by
    170:    the user, but instead is called by PCApply() if necessary.
    171: */
    172: static PetscErrorCode PCSetUp_Jacobi(PC pc)
    173: {
    174:   PC_Jacobi *jac = (PC_Jacobi *)pc->data;
    175:   Vec        diag, diagsqrt;
    176:   PetscInt   n, i;
    177:   PetscBool  zeroflag = PETSC_FALSE, negflag = PETSC_FALSE;

    179:   PetscFunctionBegin;
    180:   /*
    181:        For most preconditioners the code would begin here something like

    183:   if (!pc->setupcalled) { allocate space the first time this is ever called
    184:     PetscCall(MatCreateVecs(pc->mat,&jac->diag));
    185:   }

    187:     But for this preconditioner we want to support use of both the matrix' diagonal
    188:     elements (for left or right preconditioning) and square root of diagonal elements
    189:     (for symmetric preconditioning).  Hence we do not allocate space here, since we
    190:     don't know at this point which will be needed (diag and/or diagsqrt) until the user
    191:     applies the preconditioner, and we don't want to allocate BOTH unless we need
    192:     them both.  Thus, the diag and diagsqrt are allocated in PCSetUp_Jacobi_NonSymmetric()
    193:     and PCSetUp_Jacobi_Symmetric(), respectively.
    194:   */

    196:   /*
    197:     Here we set up the preconditioner; that is, we copy the diagonal values from
    198:     the matrix and put them into a format to make them quick to apply as a preconditioner.
    199:   */
    200:   diag     = jac->diag;
    201:   diagsqrt = jac->diagsqrt;

    203:   if (diag) {
    204:     PetscBool isset, isspd;

    206:     PetscCall(VecLockReadPop(diag));
    207:     switch (jac->type) {
    208:     case PC_JACOBI_DIAGONAL:
    209:       PetscCall(MatGetDiagonal(pc->pmat, diag));
    210:       break;
    211:     case PC_JACOBI_ROWMAX:
    212:       PetscCall(MatGetRowMaxAbs(pc->pmat, diag, NULL));
    213:       break;
    214:     case PC_JACOBI_ROWL1:
    215:       PetscCall(MatGetRowSumAbs(pc->pmat, diag));
    216:       // fix negative rows (eg, negative definite) -- this could be done for all, not needed for userowmax
    217:       PetscCall(MatIsSPDKnown(pc->pmat, &isset, &isspd));
    218:       if (jac->fixdiag && (!isset || !isspd)) {
    219:         PetscScalar       *x2;
    220:         const PetscScalar *x;
    221:         Vec                true_diag;
    222:         PetscCall(VecDuplicate(diag, &true_diag));
    223:         PetscCall(MatGetDiagonal(pc->pmat, true_diag));
    224:         PetscCall(VecGetLocalSize(diag, &n));
    225:         PetscCall(VecGetArrayWrite(diag, &x2));
    226:         PetscCall(VecGetArrayRead(true_diag, &x)); // to make more general -todo
    227:         for (i = 0; i < n; i++) {
    228:           if (PetscRealPart(x[i]) < 0.0) {
    229:             x2[i]   = -x2[i]; // flip sign to keep DA > 0
    230:             negflag = PETSC_TRUE;
    231:           }
    232:         }
    233:         PetscCall(VecRestoreArrayRead(true_diag, &x));
    234:         PetscCall(VecRestoreArrayWrite(diag, &x2));
    235:         PetscCheck(!jac->useabs || !negflag, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Jacobi use_abs and l1 not compatible with negative diagonal");
    236:         PetscCall(VecDestroy(&true_diag));
    237:       }
    238:       if (jac->scale != 1.0) {
    239:         Vec true_diag;
    240:         PetscCall(VecDuplicate(diag, &true_diag));
    241:         PetscCall(MatGetDiagonal(pc->pmat, true_diag));
    242:         PetscCall(VecAXPY(diag, -1, true_diag)); // subtract off diag
    243:         PetscCall(VecScale(diag, jac->scale));   // scale off-diag
    244:         PetscCall(VecAXPY(diag, 1, true_diag));  // add diag back in
    245:         PetscCall(VecDestroy(&true_diag));
    246:       }
    247:       break;
    248:     case PC_JACOBI_ROWSUM:
    249:       PetscCall(MatGetRowSum(pc->pmat, diag));
    250:       break;
    251:     }
    252:     PetscCall(VecReciprocal(diag));
    253:     if (jac->useabs) PetscCall(VecAbs(diag));
    254:     PetscCall(MatIsSPDKnown(pc->pmat, &isset, &isspd));
    255:     if (jac->fixdiag && (!isset || !isspd)) {
    256:       PetscScalar *x;
    257:       PetscCall(VecGetLocalSize(diag, &n));
    258:       PetscCall(VecGetArray(diag, &x));
    259:       for (i = 0; i < n; i++) {
    260:         if (x[i] == 0.0) {
    261:           x[i]     = 1.0;
    262:           zeroflag = PETSC_TRUE;
    263:         }
    264:       }
    265:       PetscCall(VecRestoreArray(diag, &x));
    266:     }
    267:     PetscCall(VecLockReadPush(diag));
    268:   }
    269:   if (diagsqrt) {
    270:     PetscScalar *x;

    272:     PetscCall(VecLockReadPop(diagsqrt));
    273:     switch (jac->type) {
    274:     case PC_JACOBI_DIAGONAL:
    275:       PetscCall(MatGetDiagonal(pc->pmat, diagsqrt));
    276:       break;
    277:     case PC_JACOBI_ROWMAX:
    278:       PetscCall(MatGetRowMaxAbs(pc->pmat, diagsqrt, NULL));
    279:       break;
    280:     case PC_JACOBI_ROWL1:
    281:       PetscCall(MatGetRowSumAbs(pc->pmat, diagsqrt));
    282:       break;
    283:     case PC_JACOBI_ROWSUM:
    284:       PetscCall(MatGetRowSum(pc->pmat, diagsqrt));
    285:       break;
    286:     }
    287:     PetscCall(VecGetLocalSize(diagsqrt, &n));
    288:     PetscCall(VecGetArray(diagsqrt, &x));
    289:     for (i = 0; i < n; i++) {
    290:       if (PetscRealPart(x[i]) < 0.0) x[i] = 1.0 / PetscSqrtReal(PetscAbsScalar(-x[i]));
    291:       else if (PetscRealPart(x[i]) > 0.0) x[i] = 1.0 / PetscSqrtReal(PetscAbsScalar(x[i]));
    292:       else {
    293:         x[i]     = 1.0;
    294:         zeroflag = PETSC_TRUE;
    295:       }
    296:     }
    297:     PetscCall(VecRestoreArray(diagsqrt, &x));
    298:     PetscCall(VecLockReadPush(diagsqrt));
    299:   }
    300:   if (zeroflag) PetscCall(PetscInfo(pc, "Zero detected in diagonal of matrix, using 1 at those locations\n"));
    301:   PetscFunctionReturn(PETSC_SUCCESS);
    302: }

    304: /*
    305:    PCSetUp_Jacobi_Symmetric - Allocates the vector needed to store the
    306:    inverse of the square root of the diagonal entries of the matrix.  This
    307:    is used for symmetric application of the Jacobi preconditioner.

    309:    Input Parameter:
    310: .  pc - the preconditioner context
    311: */
    312: static PetscErrorCode PCSetUp_Jacobi_Symmetric(PC pc)
    313: {
    314:   PC_Jacobi *jac = (PC_Jacobi *)pc->data;

    316:   PetscFunctionBegin;
    317:   PetscCall(MatCreateVecs(pc->pmat, &jac->diagsqrt, NULL));
    318:   PetscCall(VecLockReadPush(jac->diagsqrt));
    319:   PetscCall(PCSetUp_Jacobi(pc));
    320:   PetscFunctionReturn(PETSC_SUCCESS);
    321: }

    323: /*
    324:    PCSetUp_Jacobi_NonSymmetric - Allocates the vector needed to store the
    325:    inverse of the diagonal entries of the matrix.  This is used for left of
    326:    right application of the Jacobi preconditioner.

    328:    Input Parameter:
    329: .  pc - the preconditioner context
    330: */
    331: static PetscErrorCode PCSetUp_Jacobi_NonSymmetric(PC pc)
    332: {
    333:   PC_Jacobi *jac = (PC_Jacobi *)pc->data;

    335:   PetscFunctionBegin;
    336:   PetscCall(MatCreateVecs(pc->pmat, &jac->diag, NULL));
    337:   PetscCall(VecLockReadPush(jac->diag));
    338:   PetscCall(PCSetUp_Jacobi(pc));
    339:   PetscFunctionReturn(PETSC_SUCCESS);
    340: }

    342: /*
    343:    PCApply_Jacobi - Applies the Jacobi preconditioner to a vector.

    345:    Input Parameters:
    346: .  pc - the preconditioner context
    347: .  x - input vector

    349:    Output Parameter:
    350: .  y - output vector

    352:    Application Interface Routine: PCApply()
    353:  */
    354: static PetscErrorCode PCApply_Jacobi(PC pc, Vec x, Vec y)
    355: {
    356:   PC_Jacobi *jac = (PC_Jacobi *)pc->data;

    358:   PetscFunctionBegin;
    359:   if (!jac->diag) PetscCall(PCSetUp_Jacobi_NonSymmetric(pc));
    360:   PetscCall(VecPointwiseMult(y, x, jac->diag));
    361:   PetscFunctionReturn(PETSC_SUCCESS);
    362: }

    364: /*
    365:    PCApplySymmetricLeftOrRight_Jacobi - Applies the left or right part of a
    366:    symmetric preconditioner to a vector.

    368:    Input Parameters:
    369: .  pc - the preconditioner context
    370: .  x - input vector

    372:    Output Parameter:
    373: .  y - output vector

    375:    Application Interface Routines: PCApplySymmetricLeft(), PCApplySymmetricRight()
    376: */
    377: static PetscErrorCode PCApplySymmetricLeftOrRight_Jacobi(PC pc, Vec x, Vec y)
    378: {
    379:   PC_Jacobi *jac = (PC_Jacobi *)pc->data;

    381:   PetscFunctionBegin;
    382:   if (!jac->diagsqrt) PetscCall(PCSetUp_Jacobi_Symmetric(pc));
    383:   PetscCall(VecPointwiseMult(y, x, jac->diagsqrt));
    384:   PetscFunctionReturn(PETSC_SUCCESS);
    385: }

    387: static PetscErrorCode PCReset_Jacobi(PC pc)
    388: {
    389:   PC_Jacobi *jac = (PC_Jacobi *)pc->data;

    391:   PetscFunctionBegin;
    392:   if (jac->diag) PetscCall(VecLockReadPop(jac->diag));
    393:   if (jac->diagsqrt) PetscCall(VecLockReadPop(jac->diagsqrt));
    394:   PetscCall(VecDestroy(&jac->diag));
    395:   PetscCall(VecDestroy(&jac->diagsqrt));
    396:   PetscFunctionReturn(PETSC_SUCCESS);
    397: }

    399: /*
    400:    PCDestroy_Jacobi - Destroys the private context for the Jacobi preconditioner
    401:    that was created with PCCreate_Jacobi().

    403:    Input Parameter:
    404: .  pc - the preconditioner context

    406:    Application Interface Routine: PCDestroy()
    407: */
    408: static PetscErrorCode PCDestroy_Jacobi(PC pc)
    409: {
    410:   PetscFunctionBegin;
    411:   PetscCall(PCReset_Jacobi(pc));
    412:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiSetType_C", NULL));
    413:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetType_C", NULL));
    414:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiSetUseAbs_C", NULL));
    415:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetUseAbs_C", NULL));
    416:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiSetRowl1Scale_C", NULL));
    417:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetRowl1Scale_C", NULL));
    418:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiSetFixDiagonal_C", NULL));
    419:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetFixDiagonal_C", NULL));
    420:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetDiagonal_C", NULL));

    422:   /*
    423:       Free the private data structure that was hanging off the PC
    424:   */
    425:   PetscCall(PetscFree(pc->data));
    426:   PetscFunctionReturn(PETSC_SUCCESS);
    427: }

    429: static PetscErrorCode PCSetFromOptions_Jacobi(PC pc, PetscOptionItems PetscOptionsObject)
    430: {
    431:   PC_Jacobi   *jac = (PC_Jacobi *)pc->data;
    432:   PetscBool    flg;
    433:   PCJacobiType deflt, type;

    435:   PetscFunctionBegin;
    436:   PetscCall(PCJacobiGetType(pc, &deflt));
    437:   PetscOptionsHeadBegin(PetscOptionsObject, "Jacobi options");
    438:   PetscCall(PetscOptionsEnum("-pc_jacobi_type", "How to construct diagonal matrix", "PCJacobiSetType", PCJacobiTypes, (PetscEnum)deflt, (PetscEnum *)&type, &flg));
    439:   if (flg) PetscCall(PCJacobiSetType(pc, type));
    440:   PetscCall(PetscOptionsBool("-pc_jacobi_abs", "Use absolute values of diagonal entries", "PCJacobiSetUseAbs", jac->useabs, &jac->useabs, NULL));
    441:   PetscCall(PetscOptionsBool("-pc_jacobi_fixdiagonal", "Fix null terms on diagonal", "PCJacobiSetFixDiagonal", jac->fixdiag, &jac->fixdiag, NULL));
    442:   PetscCall(PetscOptionsRangeReal("-pc_jacobi_rowl1_scale", "scaling of off-diagonal elements for rowl1", "PCJacobiSetRowl1Scale", jac->scale, &jac->scale, NULL, 0.0, 1.0));
    443:   PetscOptionsHeadEnd();
    444:   PetscFunctionReturn(PETSC_SUCCESS);
    445: }

    447: static PetscErrorCode PCView_Jacobi(PC pc, PetscViewer viewer)
    448: {
    449:   PC_Jacobi *jac = (PC_Jacobi *)pc->data;
    450:   PetscBool  isascii;

    452:   PetscFunctionBegin;
    453:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    454:   if (isascii) {
    455:     PCJacobiType      type;
    456:     PetscBool         useAbs, fixdiag;
    457:     PetscViewerFormat format;
    458:     PetscReal         scale;

    460:     PetscCall(PCJacobiGetType(pc, &type));
    461:     PetscCall(PCJacobiGetUseAbs(pc, &useAbs));
    462:     PetscCall(PCJacobiGetFixDiagonal(pc, &fixdiag));
    463:     PetscCall(PCJacobiGetRowl1Scale(pc, &scale));
    464:     if (type == PC_JACOBI_ROWL1)
    465:       PetscCall(PetscViewerASCIIPrintf(viewer, "  type %s%s%s (l1-norm off-diagonal scaling %e)\n", PCJacobiTypes[type], useAbs ? ", using absolute value of entries" : "", !fixdiag ? ", not checking null diagonal entries" : "", (double)scale));
    466:     else PetscCall(PetscViewerASCIIPrintf(viewer, "  type %s%s%s\n", PCJacobiTypes[type], useAbs ? ", using absolute value of entries" : "", !fixdiag ? ", not checking null diagonal entries" : ""));
    467:     PetscCall(PetscViewerGetFormat(viewer, &format));
    468:     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL && jac->diag) PetscCall(VecView(jac->diag, viewer));
    469:   }
    470:   PetscFunctionReturn(PETSC_SUCCESS);
    471: }

    473: /*
    474:    PCCreate_Jacobi - Creates a Jacobi preconditioner context, PC_Jacobi,
    475:    and sets this as the private data within the generic preconditioning
    476:    context, PC, that was created within PCCreate().

    478:    Input Parameter:
    479: .  pc - the preconditioner context

    481:    Application Interface Routine: PCCreate()
    482: */

    484: /*MC
    485:    PCJACOBI - Jacobi (i.e. diagonal scaling preconditioning)

    487:    Options Database Keys:
    488: +  -pc_jacobi_type (diagonal|rowl1|rowmax|rowsum) - approach for forming the preconditioner
    489: .  -pc_jacobi_abs (true|false)                    - use the absolute value of the diagonal entry
    490: .  -pc_jacobi_rowl1_scale (true|false)            - scaling of off-diagonal terms
    491: -  -pc_jacobi_fixdiag (true|false)                - fix for zero diagonal terms by placing 1.0 in those locations

    493:   Level: beginner

    495:   Notes:
    496:   By using `KSPSetPCSide`(ksp,`PC_SYMMETRIC`) or -ksp_pc_side symmetric
    497:   can scale each side of the matrix by the square root of the diagonal entries.

    499:   Zero entries along the diagonal are replaced with the value 1.0

    501:   See `PCPBJACOBI` for fixed-size point block, `PCVPBJACOBI` for variable-sized point block, and `PCBJACOBI` for large size block preconditioners

    503: .seealso:  `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
    504:            `PCJacobiSetType()`, `PCJacobiSetUseAbs()`, `PCJacobiGetUseAbs()`, `PCASM`,
    505:            `PCJacobiSetFixDiagonal()`, `PCJacobiGetFixDiagonal()`,
    506:            `PCPBJACOBI`, `PCBJACOBI`, `PCVPBJACOBI`
    507: M*/

    509: PETSC_EXTERN PetscErrorCode PCCreate_Jacobi(PC pc)
    510: {
    511:   PC_Jacobi *jac;

    513:   PetscFunctionBegin;
    514:   /*
    515:      Creates the private data structure for this preconditioner and
    516:      attach it to the PC object.
    517:   */
    518:   PetscCall(PetscNew(&jac));
    519:   pc->data = (void *)jac;

    521:   /*
    522:      Initialize the pointers to vectors to ZERO; these will be used to store
    523:      diagonal entries of the matrix for fast preconditioner application.
    524:   */
    525:   jac->diag     = NULL;
    526:   jac->diagsqrt = NULL;
    527:   jac->type     = PC_JACOBI_DIAGONAL;
    528:   jac->useabs   = PETSC_FALSE;
    529:   jac->fixdiag  = PETSC_TRUE;
    530:   jac->scale    = 1.0;

    532:   /*
    533:       Set the pointers for the functions that are provided above.
    534:       Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
    535:       are called, they will automatically call these functions.  Note we
    536:       choose not to provide a couple of these functions since they are
    537:       not needed.
    538:   */
    539:   pc->ops->apply               = PCApply_Jacobi;
    540:   pc->ops->applytranspose      = PCApply_Jacobi;
    541:   pc->ops->setup               = PCSetUp_Jacobi;
    542:   pc->ops->reset               = PCReset_Jacobi;
    543:   pc->ops->destroy             = PCDestroy_Jacobi;
    544:   pc->ops->setfromoptions      = PCSetFromOptions_Jacobi;
    545:   pc->ops->view                = PCView_Jacobi;
    546:   pc->ops->applyrichardson     = NULL;
    547:   pc->ops->applysymmetricleft  = PCApplySymmetricLeftOrRight_Jacobi;
    548:   pc->ops->applysymmetricright = PCApplySymmetricLeftOrRight_Jacobi;

    550:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiSetType_C", PCJacobiSetType_Jacobi));
    551:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetType_C", PCJacobiGetType_Jacobi));
    552:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiSetRowl1Scale_C", PCJacobiSetRowl1Scale_Jacobi));
    553:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetRowl1Scale_C", PCJacobiGetRowl1Scale_Jacobi));
    554:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiSetUseAbs_C", PCJacobiSetUseAbs_Jacobi));
    555:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetUseAbs_C", PCJacobiGetUseAbs_Jacobi));
    556:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiSetFixDiagonal_C", PCJacobiSetFixDiagonal_Jacobi));
    557:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetFixDiagonal_C", PCJacobiGetFixDiagonal_Jacobi));
    558:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetDiagonal_C", PCJacobiGetDiagonal_Jacobi));
    559:   PetscFunctionReturn(PETSC_SUCCESS);
    560: }

    562: /*@
    563:   PCJacobiSetUseAbs - Causes the Jacobi preconditioner `PCJACOBI` to use the
    564:   absolute values of the diagonal divisors in the preconditioner

    566:   Logically Collective

    568:   Input Parameters:
    569: + pc  - the preconditioner context
    570: - flg - whether to use absolute values or not

    572:   Options Database Key:
    573: . -pc_jacobi_abs (true|false) - use absolute values

    575:   Note:
    576:   This takes affect at the next construction of the preconditioner

    578:   Level: intermediate

    580: .seealso: [](ch_ksp), `PCJACOBI`, `PCJacobiSetType()`, `PCJacobiGetUseAbs()`
    581: @*/
    582: PetscErrorCode PCJacobiSetUseAbs(PC pc, PetscBool flg)
    583: {
    584:   PetscFunctionBegin;
    586:   PetscTryMethod(pc, "PCJacobiSetUseAbs_C", (PC, PetscBool), (pc, flg));
    587:   PetscFunctionReturn(PETSC_SUCCESS);
    588: }

    590: /*@
    591:   PCJacobiGetUseAbs - Determines if the Jacobi preconditioner `PCJACOBI` uses the
    592:   absolute values of the diagonal divisors in the preconditioner

    594:   Logically Collective

    596:   Input Parameter:
    597: . pc - the preconditioner context

    599:   Output Parameter:
    600: . flg - whether to use absolute values or not

    602:   Level: intermediate

    604: .seealso: [](ch_ksp), `PCJACOBI`, `PCJacobiSetType()`, `PCJacobiSetUseAbs()`, `PCJacobiGetType()`
    605: @*/
    606: PetscErrorCode PCJacobiGetUseAbs(PC pc, PetscBool *flg)
    607: {
    608:   PetscFunctionBegin;
    610:   PetscUseMethod(pc, "PCJacobiGetUseAbs_C", (PC, PetscBool *), (pc, flg));
    611:   PetscFunctionReturn(PETSC_SUCCESS);
    612: }

    614: /*@
    615:   PCJacobiSetRowl1Scale - Set scaling of off-diagonal of operator when computing l1 row norms, eg,
    616:   Remark 6.1 in "Multigrid Smoothers for Ultraparallel Computing", Baker et al, with 0.5 scaling

    618:   Logically Collective

    620:   Input Parameters:
    621: + pc    - the preconditioner context
    622: - scale - scaling

    624:   Options Database Key:
    625: . -pc_jacobi_rowl1_scale scale - scaling of off diagonal values

    627:   Level: intermediate

    629: .seealso: [](ch_ksp), `PCJACOBI`, `PCJacobiSetType()`, `PCJacobiGetRowl1Scale()`
    630: @*/
    631: PetscErrorCode PCJacobiSetRowl1Scale(PC pc, PetscReal scale)
    632: {
    633:   PetscFunctionBegin;
    635:   PetscTryMethod(pc, "PCJacobiSetRowl1Scale_C", (PC, PetscReal), (pc, scale));
    636:   PetscFunctionReturn(PETSC_SUCCESS);
    637: }

    639: /*@
    640:   PCJacobiGetRowl1Scale - Get scaling of off-diagonal elements summed into l1-norm diagonal

    642:   Logically Collective

    644:   Input Parameter:
    645: . pc - the preconditioner context

    647:   Output Parameter:
    648: . scale - scaling

    650:   Level: intermediate

    652: .seealso: [](ch_ksp), `PCJACOBI`, `PCJacobiSetType()`, `PCJacobiSetRowl1Scale()`, `PCJacobiGetType()`
    653: @*/
    654: PetscErrorCode PCJacobiGetRowl1Scale(PC pc, PetscReal *scale)
    655: {
    656:   PetscFunctionBegin;
    658:   PetscUseMethod(pc, "PCJacobiGetRowl1Scale_C", (PC, PetscReal *), (pc, scale));
    659:   PetscFunctionReturn(PETSC_SUCCESS);
    660: }

    662: /*@
    663:   PCJacobiSetFixDiagonal - Check for zero values on the diagonal and replace them with 1.0

    665:   Logically Collective

    667:   Input Parameters:
    668: + pc  - the preconditioner context
    669: - flg - the boolean flag

    671:   Options Database Key:
    672: . -pc_jacobi_fixdiagonal (true|false) - check for zero values on the diagonal

    674:   Note:
    675:   This takes affect at the next construction of the preconditioner

    677:   Level: intermediate

    679: .seealso: [](ch_ksp), `PCJACOBI`, `PCJacobiSetType()`, `PCJacobiGetFixDiagonal()`, `PCJacobiSetUseAbs()`
    680: @*/
    681: PetscErrorCode PCJacobiSetFixDiagonal(PC pc, PetscBool flg)
    682: {
    683:   PetscFunctionBegin;
    685:   PetscTryMethod(pc, "PCJacobiSetFixDiagonal_C", (PC, PetscBool), (pc, flg));
    686:   PetscFunctionReturn(PETSC_SUCCESS);
    687: }

    689: /*@
    690:   PCJacobiGetFixDiagonal - Determines if the Jacobi preconditioner `PCJACOBI` checks for zero diagonal terms

    692:   Logically Collective

    694:   Input Parameter:
    695: . pc - the preconditioner context

    697:   Output Parameter:
    698: . flg - the boolean flag

    700:   Options Database Key:
    701: . -pc_jacobi_fixdiagonal (true|false) - Fix 0 terms on diagonal by using 1

    703:   Level: intermediate

    705: .seealso: [](ch_ksp), `PCJACOBI`, `PCJacobiSetType()`, `PCJacobiSetFixDiagonal()`
    706: @*/
    707: PetscErrorCode PCJacobiGetFixDiagonal(PC pc, PetscBool *flg)
    708: {
    709:   PetscFunctionBegin;
    711:   PetscUseMethod(pc, "PCJacobiGetFixDiagonal_C", (PC, PetscBool *), (pc, flg));
    712:   PetscFunctionReturn(PETSC_SUCCESS);
    713: }

    715: /*@
    716:   PCJacobiGetDiagonal - Returns copy of the diagonal and/or diagonal squareroot `Vec`

    718:   Logically Collective

    720:   Input Parameter:
    721: . pc - the preconditioner context

    723:   Output Parameters:
    724: + diagonal      - Copy of `Vec` of the inverted diagonal
    725: - diagonal_sqrt - Copy of `Vec` of the inverted square root diagonal

    727:   Level: developer

    729: .seealso: [](ch_ksp), `PCJACOBI`, `PCJacobiSetType()`
    730: @*/
    731: PetscErrorCode PCJacobiGetDiagonal(PC pc, Vec diagonal, Vec diagonal_sqrt)
    732: {
    733:   PetscFunctionBegin;
    735:   PetscUseMethod(pc, "PCJacobiGetDiagonal_C", (PC, Vec, Vec), (pc, diagonal, diagonal_sqrt));
    736:   PetscFunctionReturn(PETSC_SUCCESS);
    737: }

    739: /*@
    740:   PCJacobiSetType - Causes the Jacobi preconditioner to use either the diagonal, the maximum entry in each row,
    741:   of the sum of rows entries for the diagonal preconditioner

    743:   Logically Collective

    745:   Input Parameters:
    746: + pc   - the preconditioner context
    747: - type - `PC_JACOBI_DIAGONAL`, `PC_JACOBI_ROWL1`, `PC_JACOBI_ROWMAX`, `PC_JACOBI_ROWSUM`

    749:   Options Database Key:
    750: . -pc_jacobi_type [diagonal, rowl1, rowmax, rowsum] - the type of diagonal matrix to use for Jacobi

    752:   Level: intermediate

    754:   Developer Notes:
    755:   Why is there a separate function for using the absolute value?

    757: .seealso: [](ch_ksp), `PCJACOBI`, `PCJacobiSetUseAbs()`, `PCJacobiGetType()`
    758: @*/
    759: PetscErrorCode PCJacobiSetType(PC pc, PCJacobiType type)
    760: {
    761:   PetscFunctionBegin;
    763:   PetscTryMethod(pc, "PCJacobiSetType_C", (PC, PCJacobiType), (pc, type));
    764:   PetscFunctionReturn(PETSC_SUCCESS);
    765: }

    767: /*@
    768:   PCJacobiGetType - Gets how the diagonal matrix is produced for the preconditioner

    770:   Not Collective

    772:   Input Parameter:
    773: . pc - the preconditioner context

    775:   Output Parameter:
    776: . type - `PC_JACOBI_DIAGONAL`, `PC_JACOBI_ROWL1`, `PC_JACOBI_ROWMAX`, `PC_JACOBI_ROWSUM`

    778:   Level: intermediate

    780: .seealso: [](ch_ksp), `PCJACOBI`, `PCJacobiSetUseAbs()`, `PCJacobiSetType()`
    781: @*/
    782: PetscErrorCode PCJacobiGetType(PC pc, PCJacobiType *type)
    783: {
    784:   PetscFunctionBegin;
    786:   PetscUseMethod(pc, "PCJacobiGetType_C", (PC, PCJacobiType *), (pc, type));
    787:   PetscFunctionReturn(PETSC_SUCCESS);
    788: }