Actual source code: jacobi.c

      2: /*  --------------------------------------------------------------------

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

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

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

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

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

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

     44:     -------------------------------------------------------------------- */

     46: /*
     47:    Include files needed for the Jacobi preconditioner:
     48:      pcimpl.h - private include file intended for use by all preconditioners
     49: */

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

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

     55: /*
     56:    Private context (data structure) for the Jacobi preconditioner.
     57: */
     58: typedef struct {
     59:   Vec diag;                      /* vector containing the reciprocals of the diagonal elements of the preconditioner matrix */
     60:   Vec diagsqrt;                  /* vector containing the reciprocals of the square roots of
     61:                                     the diagonal elements of the preconditioner matrix (used
     62:                                     only for symmetric preconditioner application) */
     63:   PetscBool userowmax;           /* set with PCJacobiSetType() */
     64:   PetscBool userowsum;
     65:   PetscBool useabs;              /* use the absolute values of the diagonal entries */
     66: } PC_Jacobi;

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

     73:   j->userowmax = PETSC_FALSE;
     74:   j->userowsum = PETSC_FALSE;
     75:   if (type == PC_JACOBI_ROWMAX) {
     76:     j->userowmax = PETSC_TRUE;
     77:   } else if (type == PC_JACOBI_ROWSUM) {
     78:     j->userowsum = PETSC_TRUE;
     79:   }
     80:   return(0);
     81: }

     83: static PetscErrorCode  PCJacobiGetType_Jacobi(PC pc,PCJacobiType *type)
     84: {
     85:   PC_Jacobi *j = (PC_Jacobi*)pc->data;

     88:   if (j->userowmax) {
     89:     *type = PC_JACOBI_ROWMAX;
     90:   } else if (j->userowsum) {
     91:     *type = PC_JACOBI_ROWSUM;
     92:   } else {
     93:     *type = PC_JACOBI_DIAGONAL;
     94:   }
     95:   return(0);
     96: }

     98: static PetscErrorCode  PCJacobiSetUseAbs_Jacobi(PC pc,PetscBool flg)
     99: {
    100:   PC_Jacobi *j = (PC_Jacobi*)pc->data;

    103:   j->useabs = flg;
    104:   return(0);
    105: }

    107: static PetscErrorCode  PCJacobiGetUseAbs_Jacobi(PC pc,PetscBool *flg)
    108: {
    109:   PC_Jacobi *j = (PC_Jacobi*)pc->data;

    112:   *flg = j->useabs;
    113:   return(0);
    114: }

    116: /* -------------------------------------------------------------------------- */
    117: /*
    118:    PCSetUp_Jacobi - Prepares for the use of the Jacobi preconditioner
    119:                     by setting data structures and options.

    121:    Input Parameter:
    122: .  pc - the preconditioner context

    124:    Application Interface Routine: PCSetUp()

    126:    Notes:
    127:    The interface routine PCSetUp() is not usually called directly by
    128:    the user, but instead is called by PCApply() if necessary.
    129: */
    130: static PetscErrorCode PCSetUp_Jacobi(PC pc)
    131: {
    132:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;
    133:   Vec            diag,diagsqrt;
    135:   PetscInt       n,i;
    136:   PetscScalar    *x;
    137:   PetscBool      zeroflag = PETSC_FALSE;

    140:   /*
    141:        For most preconditioners the code would begin here something like

    143:   if (pc->setupcalled == 0) { allocate space the first time this is ever called
    144:     MatCreateVecs(pc->mat,&jac->diag);
    145:     PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->diag);
    146:   }

    148:     But for this preconditioner we want to support use of both the matrix' diagonal
    149:     elements (for left or right preconditioning) and square root of diagonal elements
    150:     (for symmetric preconditioning).  Hence we do not allocate space here, since we
    151:     don't know at this point which will be needed (diag and/or diagsqrt) until the user
    152:     applies the preconditioner, and we don't want to allocate BOTH unless we need
    153:     them both.  Thus, the diag and diagsqrt are allocated in PCSetUp_Jacobi_NonSymmetric()
    154:     and PCSetUp_Jacobi_Symmetric(), respectively.
    155:   */

    157:   /*
    158:     Here we set up the preconditioner; that is, we copy the diagonal values from
    159:     the matrix and put them into a format to make them quick to apply as a preconditioner.
    160:   */
    161:   diag     = jac->diag;
    162:   diagsqrt = jac->diagsqrt;

    164:   if (diag) {
    165:     if (jac->userowmax) {
    166:       MatGetRowMaxAbs(pc->pmat,diag,NULL);
    167:     } else if (jac->userowsum) {
    168:       MatGetRowSum(pc->pmat,diag);
    169:     } else {
    170:       MatGetDiagonal(pc->pmat,diag);
    171:     }
    172:     VecReciprocal(diag);
    173:     VecGetLocalSize(diag,&n);
    174:     if (jac->useabs) {
    175:       VecAbs(diag);
    176:     }
    177:     VecGetArray(diag,&x);
    178:     for (i=0; i<n; i++) {
    179:       if (x[i] == 0.0) {
    180:         x[i]     = 1.0;
    181:         zeroflag = PETSC_TRUE;
    182:       }
    183:     }
    184:     VecRestoreArray(diag,&x);
    185:   }
    186:   if (diagsqrt) {
    187:     if (jac->userowmax) {
    188:       MatGetRowMaxAbs(pc->pmat,diagsqrt,NULL);
    189:     } else if (jac->userowsum) {
    190:       MatGetRowSum(pc->pmat,diagsqrt);
    191:     } else {
    192:       MatGetDiagonal(pc->pmat,diagsqrt);
    193:     }
    194:     VecGetLocalSize(diagsqrt,&n);
    195:     VecGetArray(diagsqrt,&x);
    196:     for (i=0; i<n; i++) {
    197:       if (x[i] != 0.0) x[i] = 1.0/PetscSqrtReal(PetscAbsScalar(x[i]));
    198:       else {
    199:         x[i]     = 1.0;
    200:         zeroflag = PETSC_TRUE;
    201:       }
    202:     }
    203:     VecRestoreArray(diagsqrt,&x);
    204:   }
    205:   if (zeroflag) {
    206:     PetscInfo(pc,"Zero detected in diagonal of matrix, using 1 at those locations\n");
    207:   }
    208:   return(0);
    209: }
    210: /* -------------------------------------------------------------------------- */
    211: /*
    212:    PCSetUp_Jacobi_Symmetric - Allocates the vector needed to store the
    213:    inverse of the square root of the diagonal entries of the matrix.  This
    214:    is used for symmetric application of the Jacobi preconditioner.

    216:    Input Parameter:
    217: .  pc - the preconditioner context
    218: */
    219: static PetscErrorCode PCSetUp_Jacobi_Symmetric(PC pc)
    220: {
    222:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

    225:   MatCreateVecs(pc->pmat,&jac->diagsqrt,NULL);
    226:   PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->diagsqrt);
    227:   PCSetUp_Jacobi(pc);
    228:   return(0);
    229: }
    230: /* -------------------------------------------------------------------------- */
    231: /*
    232:    PCSetUp_Jacobi_NonSymmetric - Allocates the vector needed to store the
    233:    inverse of the diagonal entries of the matrix.  This is used for left of
    234:    right application of the Jacobi preconditioner.

    236:    Input Parameter:
    237: .  pc - the preconditioner context
    238: */
    239: static PetscErrorCode PCSetUp_Jacobi_NonSymmetric(PC pc)
    240: {
    242:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

    245:   MatCreateVecs(pc->pmat,&jac->diag,NULL);
    246:   PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->diag);
    247:   PCSetUp_Jacobi(pc);
    248:   return(0);
    249: }
    250: /* -------------------------------------------------------------------------- */
    251: /*
    252:    PCApply_Jacobi - Applies the Jacobi preconditioner to a vector.

    254:    Input Parameters:
    255: .  pc - the preconditioner context
    256: .  x - input vector

    258:    Output Parameter:
    259: .  y - output vector

    261:    Application Interface Routine: PCApply()
    262:  */
    263: static PetscErrorCode PCApply_Jacobi(PC pc,Vec x,Vec y)
    264: {
    265:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

    269:   if (!jac->diag) {
    270:     PCSetUp_Jacobi_NonSymmetric(pc);
    271:   }
    272:   VecPointwiseMult(y,x,jac->diag);
    273:   return(0);
    274: }
    275: /* -------------------------------------------------------------------------- */
    276: /*
    277:    PCApplySymmetricLeftOrRight_Jacobi - Applies the left or right part of a
    278:    symmetric preconditioner to a vector.

    280:    Input Parameters:
    281: .  pc - the preconditioner context
    282: .  x - input vector

    284:    Output Parameter:
    285: .  y - output vector

    287:    Application Interface Routines: PCApplySymmetricLeft(), PCApplySymmetricRight()
    288: */
    289: static PetscErrorCode PCApplySymmetricLeftOrRight_Jacobi(PC pc,Vec x,Vec y)
    290: {
    292:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

    295:   if (!jac->diagsqrt) {
    296:     PCSetUp_Jacobi_Symmetric(pc);
    297:   }
    298:   VecPointwiseMult(y,x,jac->diagsqrt);
    299:   return(0);
    300: }
    301: /* -------------------------------------------------------------------------- */
    302: static PetscErrorCode PCReset_Jacobi(PC pc)
    303: {
    304:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

    308:   VecDestroy(&jac->diag);
    309:   VecDestroy(&jac->diagsqrt);
    310:   return(0);
    311: }

    313: /*
    314:    PCDestroy_Jacobi - Destroys the private context for the Jacobi preconditioner
    315:    that was created with PCCreate_Jacobi().

    317:    Input Parameter:
    318: .  pc - the preconditioner context

    320:    Application Interface Routine: PCDestroy()
    321: */
    322: static PetscErrorCode PCDestroy_Jacobi(PC pc)
    323: {

    327:   PCReset_Jacobi(pc);

    329:   /*
    330:       Free the private data structure that was hanging off the PC
    331:   */
    332:   PetscFree(pc->data);
    333:   return(0);
    334: }

    336: static PetscErrorCode PCSetFromOptions_Jacobi(PetscOptionItems *PetscOptionsObject,PC pc)
    337: {
    338:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;
    340:   PetscBool      flg;
    341:   PCJacobiType   deflt,type;

    344:   PCJacobiGetType(pc,&deflt);
    345:   PetscOptionsHead(PetscOptionsObject,"Jacobi options");
    346:   PetscOptionsEnum("-pc_jacobi_type","How to construct diagonal matrix","PCJacobiSetType",PCJacobiTypes,(PetscEnum)deflt,(PetscEnum*)&type,&flg);
    347:   if (flg) {
    348:     PCJacobiSetType(pc,type);
    349:   }
    350:   PetscOptionsBool("-pc_jacobi_abs","Use absolute values of diagaonal entries","PCJacobiSetUseAbs",jac->useabs,&jac->useabs,NULL);
    351:   PetscOptionsTail();
    352:   return(0);
    353: }

    355: static PetscErrorCode PCView_Jacobi(PC pc, PetscViewer viewer)
    356: {
    357:   PC_Jacobi     *jac = (PC_Jacobi *) pc->data;
    358:   PetscBool      iascii;

    362:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
    363:   if (iascii) {
    364:     PCJacobiType      type;
    365:     PetscBool         useAbs;
    366:     PetscViewerFormat format;

    368:     PCJacobiGetType(pc, &type);
    369:     PCJacobiGetUseAbs(pc, &useAbs);
    370:     PetscViewerASCIIPrintf(viewer, "  type %s%s\n", PCJacobiTypes[type], useAbs ? ", using absolute value of entries" : "");
    371:     PetscViewerGetFormat(viewer, &format);
    372:     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
    373:       VecView(jac->diag, viewer);
    374:     }
    375:   }
    376:   return(0);
    377: }

    379: /* -------------------------------------------------------------------------- */
    380: /*
    381:    PCCreate_Jacobi - Creates a Jacobi preconditioner context, PC_Jacobi,
    382:    and sets this as the private data within the generic preconditioning
    383:    context, PC, that was created within PCCreate().

    385:    Input Parameter:
    386: .  pc - the preconditioner context

    388:    Application Interface Routine: PCCreate()
    389: */

    391: /*MC
    392:      PCJACOBI - Jacobi (i.e. diagonal scaling preconditioning)

    394:    Options Database Key:
    395: +    -pc_jacobi_type <diagonal,rowmax,rowsum> - approach for forming the preconditioner
    396: -    -pc_jacobi_abs - use the absolute value of the diagonal entry

    398:    Level: beginner

    400:   Notes:
    401:     By using KSPSetPCSide(ksp,PC_SYMMETRIC) or -ksp_pc_side symmetric
    402:          can scale each side of the matrix by the square root of the diagonal entries.

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

    406:          See PCPBJACOBI for a point-block Jacobi preconditioner

    408: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
    409:            PCJacobiSetType(), PCJacobiSetUseAbs(), PCJacobiGetUseAbs(), PCPBJACOBI
    410: M*/

    412: PETSC_EXTERN PetscErrorCode PCCreate_Jacobi(PC pc)
    413: {
    414:   PC_Jacobi      *jac;

    418:   /*
    419:      Creates the private data structure for this preconditioner and
    420:      attach it to the PC object.
    421:   */
    422:   PetscNewLog(pc,&jac);
    423:   pc->data = (void*)jac;

    425:   /*
    426:      Initialize the pointers to vectors to ZERO; these will be used to store
    427:      diagonal entries of the matrix for fast preconditioner application.
    428:   */
    429:   jac->diag      = NULL;
    430:   jac->diagsqrt  = NULL;
    431:   jac->userowmax = PETSC_FALSE;
    432:   jac->userowsum = PETSC_FALSE;
    433:   jac->useabs    = PETSC_FALSE;

    435:   /*
    436:       Set the pointers for the functions that are provided above.
    437:       Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
    438:       are called, they will automatically call these functions.  Note we
    439:       choose not to provide a couple of these functions since they are
    440:       not needed.
    441:   */
    442:   pc->ops->apply               = PCApply_Jacobi;
    443:   pc->ops->applytranspose      = PCApply_Jacobi;
    444:   pc->ops->setup               = PCSetUp_Jacobi;
    445:   pc->ops->reset               = PCReset_Jacobi;
    446:   pc->ops->destroy             = PCDestroy_Jacobi;
    447:   pc->ops->setfromoptions      = PCSetFromOptions_Jacobi;
    448:   pc->ops->view                = PCView_Jacobi;
    449:   pc->ops->applyrichardson     = NULL;
    450:   pc->ops->applysymmetricleft  = PCApplySymmetricLeftOrRight_Jacobi;
    451:   pc->ops->applysymmetricright = PCApplySymmetricLeftOrRight_Jacobi;

    453:   PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetType_C",PCJacobiSetType_Jacobi);
    454:   PetscObjectComposeFunction((PetscObject)pc,"PCJacobiGetType_C",PCJacobiGetType_Jacobi);
    455:   PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetUseAbs_C",PCJacobiSetUseAbs_Jacobi);
    456:   PetscObjectComposeFunction((PetscObject)pc,"PCJacobiGetUseAbs_C",PCJacobiGetUseAbs_Jacobi);
    457:   return(0);
    458: }

    460: /*@
    461:    PCJacobiSetUseAbs - Causes the Jacobi preconditioner to use the
    462:       absolute values of the digonal divisors in the preconditioner

    464:    Logically Collective on PC

    466:    Input Parameters:
    467: +  pc - the preconditioner context
    468: -  flg - whether to use absolute values or not

    470:    Options Database Key:
    471: .  -pc_jacobi_abs

    473:    Notes:
    474:     This takes affect at the next construction of the preconditioner

    476:    Level: intermediate

    478: .seealso: PCJacobiaSetType(), PCJacobiGetUseAbs()

    480: @*/
    481: PetscErrorCode  PCJacobiSetUseAbs(PC pc,PetscBool flg)
    482: {

    487:   PetscTryMethod(pc,"PCJacobiSetUseAbs_C",(PC,PetscBool),(pc,flg));
    488:   return(0);
    489: }

    491: /*@
    492:    PCJacobiGetUseAbs - Determines if the Jacobi preconditioner uses the
    493:       absolute values of the digonal divisors in the preconditioner

    495:    Logically Collective on PC

    497:    Input Parameter:
    498: .  pc - the preconditioner context

    500:    Output Parameter:
    501: .  flg - whether to use absolute values or not

    503:    Options Database Key:
    504: .  -pc_jacobi_abs

    506:    Level: intermediate

    508: .seealso: PCJacobiaSetType(), PCJacobiSetUseAbs(), PCJacobiGetType()

    510: @*/
    511: PetscErrorCode  PCJacobiGetUseAbs(PC pc,PetscBool *flg)
    512: {

    517:   PetscUseMethod(pc,"PCJacobiGetUseAbs_C",(PC,PetscBool*),(pc,flg));
    518:   return(0);
    519: }

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

    525:    Logically Collective on PC

    527:    Input Parameters:
    528: +  pc - the preconditioner context
    529: -  type - PC_JACOBI_DIAGONAL, PC_JACOBI_ROWMAX, PC_JACOBI_ROWSUM

    531:    Options Database Key:
    532: .  -pc_jacobi_type <diagonal,rowmax,rowsum>

    534:    Level: intermediate

    536: .seealso: PCJacobiaUseAbs(), PCJacobiGetType()
    537: @*/
    538: PetscErrorCode  PCJacobiSetType(PC pc,PCJacobiType type)
    539: {

    544:   PetscTryMethod(pc,"PCJacobiSetType_C",(PC,PCJacobiType),(pc,type));
    545:   return(0);
    546: }

    548: /*@
    549:    PCJacobiGetType - Gets how the diagonal matrix is produced for the preconditioner

    551:    Not Collective on PC

    553:    Input Parameter:
    554: .  pc - the preconditioner context

    556:    Output Parameter:
    557: .  type - PC_JACOBI_DIAGONAL, PC_JACOBI_ROWMAX, PC_JACOBI_ROWSUM

    559:    Level: intermediate

    561: .seealso: PCJacobiaUseAbs(), PCJacobiSetType()
    562: @*/
    563: PetscErrorCode  PCJacobiGetType(PC pc,PCJacobiType *type)
    564: {

    569:   PetscUseMethod(pc,"PCJacobiGetType_C",(PC,PCJacobiType*),(pc,type));
    570:   return(0);
    571: }