Actual source code: jacobi.c

      1: #define PETSCKSP_DLL

      3: /*  --------------------------------------------------------------------

      5:      This file implements a Jacobi preconditioner for matrices that use
      6:      the Mat interface (various matrix formats).  Actually, the only
      7:      matrix operation used here is MatGetDiagonal(), which extracts
      8:      diagonal elements of the preconditioning matrix.

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

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

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

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

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

     45:     -------------------------------------------------------------------- */

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

     52:  #include src/ksp/pc/pcimpl.h

     54: /*
     55:    Private context (data structure) for the Jacobi preconditioner.
     56: */
     57: typedef struct {
     58:   Vec        diag;               /* vector containing the reciprocals of the diagonal elements
     59:                                     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:   PetscTruth userowmax;
     64: } PC_Jacobi;

     69: PetscErrorCode PETSCKSP_DLLEXPORT PCJacobiSetUseRowMax_Jacobi(PC pc)
     70: {
     71:   PC_Jacobi *j;

     74:   j            = (PC_Jacobi*)pc->data;
     75:   j->userowmax = PETSC_TRUE;
     76:   return(0);
     77: }

     80: /* -------------------------------------------------------------------------- */
     81: /*
     82:    PCSetUp_Jacobi - Prepares for the use of the Jacobi preconditioner
     83:                     by setting data structures and options.

     85:    Input Parameter:
     86: .  pc - the preconditioner context

     88:    Application Interface Routine: PCSetUp()

     90:    Notes:
     91:    The interface routine PCSetUp() is not usually called directly by
     92:    the user, but instead is called by PCApply() if necessary.
     93: */
     96: static PetscErrorCode PCSetUp_Jacobi(PC pc)
     97: {
     98:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;
     99:   Vec            diag,diagsqrt;
    101:   PetscInt       n,i,zeroflag=0;
    102:   PetscScalar   *x;

    105:   /*
    106:        For most preconditioners the code would begin here something like

    108:   if (pc->setupcalled == 0) { allocate space the first time this is ever called
    109:     MatGetVecs(pc->mat,&jac->diag);
    110:     PetscLogObjectParent(pc,jac->diag);
    111:   }

    113:     But for this preconditioner we want to support use of both the matrix' diagonal
    114:     elements (for left or right preconditioning) and square root of diagonal elements
    115:     (for symmetric preconditioning).  Hence we do not allocate space here, since we
    116:     don't know at this point which will be needed (diag and/or diagsqrt) until the user
    117:     applies the preconditioner, and we don't want to allocate BOTH unless we need
    118:     them both.  Thus, the diag and diagsqrt are allocated in PCSetUp_Jacobi_NonSymmetric()
    119:     and PCSetUp_Jacobi_Symmetric(), respectively.
    120:   */

    122:   /*
    123:     Here we set up the preconditioner; that is, we copy the diagonal values from
    124:     the matrix and put them into a format to make them quick to apply as a preconditioner.
    125:   */
    126:   diag     = jac->diag;
    127:   diagsqrt = jac->diagsqrt;

    129:   if (diag) {
    130:     if (jac->userowmax) {
    131:       MatGetRowMax(pc->pmat,diag);
    132:     } else {
    133:       MatGetDiagonal(pc->pmat,diag);
    134:     }
    135:     VecReciprocal(diag);
    136:     VecGetLocalSize(diag,&n);
    137:     VecGetArray(diag,&x);
    138:     for (i=0; i<n; i++) {
    139:       if (x[i] == 0.0) {
    140:         x[i]     = 1.0;
    141:         zeroflag = 1;
    142:       }
    143:     }
    144:     VecRestoreArray(diag,&x);
    145:   }
    146:   if (diagsqrt) {
    147:     if (jac->userowmax) {
    148:       MatGetRowMax(pc->pmat,diagsqrt);
    149:     } else {
    150:       MatGetDiagonal(pc->pmat,diagsqrt);
    151:     }
    152:     VecGetLocalSize(diagsqrt,&n);
    153:     VecGetArray(diagsqrt,&x);
    154:     for (i=0; i<n; i++) {
    155:       if (x[i] != 0.0) x[i] = 1.0/sqrt(PetscAbsScalar(x[i]));
    156:       else {
    157:         x[i]     = 1.0;
    158:         zeroflag = 1;
    159:       }
    160:     }
    161:     VecRestoreArray(diagsqrt,&x);
    162:   }
    163:   if (zeroflag) {
    164:     PetscLogInfo((pc,"PCSetUp_Jacobi:Zero detected in diagonal of matrix, using 1 at those locations\n"));
    165:   }
    166:   return(0);
    167: }
    168: /* -------------------------------------------------------------------------- */
    169: /*
    170:    PCSetUp_Jacobi_Symmetric - Allocates the vector needed to store the
    171:    inverse of the square root of the diagonal entries of the matrix.  This
    172:    is used for symmetric application of the Jacobi preconditioner.

    174:    Input Parameter:
    175: .  pc - the preconditioner context
    176: */
    179: static PetscErrorCode PCSetUp_Jacobi_Symmetric(PC pc)
    180: {
    182:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

    185:   MatGetVecs(pc->pmat,&jac->diagsqrt,0);
    186:   PetscLogObjectParent(pc,jac->diagsqrt);
    187:   PCSetUp_Jacobi(pc);
    188:   return(0);
    189: }
    190: /* -------------------------------------------------------------------------- */
    191: /*
    192:    PCSetUp_Jacobi_NonSymmetric - Allocates the vector needed to store the
    193:    inverse of the diagonal entries of the matrix.  This is used for left of
    194:    right application of the Jacobi preconditioner.

    196:    Input Parameter:
    197: .  pc - the preconditioner context
    198: */
    201: static PetscErrorCode PCSetUp_Jacobi_NonSymmetric(PC pc)
    202: {
    204:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

    207:   MatGetVecs(pc->pmat,&jac->diag,0);
    208:   PetscLogObjectParent(pc,jac->diag);
    209:   PCSetUp_Jacobi(pc);
    210:   return(0);
    211: }
    212: /* -------------------------------------------------------------------------- */
    213: /*
    214:    PCApply_Jacobi - Applies the Jacobi preconditioner to a vector.

    216:    Input Parameters:
    217: .  pc - the preconditioner context
    218: .  x - input vector

    220:    Output Parameter:
    221: .  y - output vector

    223:    Application Interface Routine: PCApply()
    224:  */
    227: static PetscErrorCode PCApply_Jacobi(PC pc,Vec x,Vec y)
    228: {
    229:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

    233:   if (!jac->diag) {
    234:     PCSetUp_Jacobi_NonSymmetric(pc);
    235:   }
    236:   VecPointwiseMult(y,x,jac->diag);
    237:   return(0);
    238: }
    239: /* -------------------------------------------------------------------------- */
    240: /*
    241:    PCApplySymmetricLeftOrRight_Jacobi - Applies the left or right part of a
    242:    symmetric preconditioner to a vector.

    244:    Input Parameters:
    245: .  pc - the preconditioner context
    246: .  x - input vector

    248:    Output Parameter:
    249: .  y - output vector

    251:    Application Interface Routines: PCApplySymmetricLeft(), PCApplySymmetricRight()
    252: */
    255: static PetscErrorCode PCApplySymmetricLeftOrRight_Jacobi(PC pc,Vec x,Vec y)
    256: {
    258:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

    261:   if (!jac->diagsqrt) {
    262:     PCSetUp_Jacobi_Symmetric(pc);
    263:   }
    264:   VecPointwiseMult(y,x,jac->diagsqrt);
    265:   return(0);
    266: }
    267: /* -------------------------------------------------------------------------- */
    268: /*
    269:    PCDestroy_Jacobi - Destroys the private context for the Jacobi preconditioner
    270:    that was created with PCCreate_Jacobi().

    272:    Input Parameter:
    273: .  pc - the preconditioner context

    275:    Application Interface Routine: PCDestroy()
    276: */
    279: static PetscErrorCode PCDestroy_Jacobi(PC pc)
    280: {
    281:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

    285:   if (jac->diag)     {VecDestroy(jac->diag);}
    286:   if (jac->diagsqrt) {VecDestroy(jac->diagsqrt);}

    288:   /*
    289:       Free the private data structure that was hanging off the PC
    290:   */
    291:   PetscFree(jac);
    292:   return(0);
    293: }

    297: static PetscErrorCode PCSetFromOptions_Jacobi(PC pc)
    298: {
    299:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

    303:   PetscOptionsHead("Jacobi options");
    304:     PetscOptionsTruth("-pc_jacobi_rowmax","Use row maximums for diagonal","PCJacobiSetUseRowMax",jac->userowmax,
    305:                           &jac->userowmax,PETSC_NULL);
    306:   PetscOptionsTail();
    307:   return(0);
    308: }

    310: /* -------------------------------------------------------------------------- */
    311: /*
    312:    PCCreate_Jacobi - Creates a Jacobi preconditioner context, PC_Jacobi,
    313:    and sets this as the private data within the generic preconditioning
    314:    context, PC, that was created within PCCreate().

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

    319:    Application Interface Routine: PCCreate()
    320: */

    322: /*MC
    323:      PCJacobi - Jacobi (i.e. diagonal scaling preconditioning)

    325:    Options Database Key:
    326: .    -pc_jacobi_rowmax - use the maximum absolute value in each row as the scaling factor,
    327:                         rather than the diagonal

    329:    Level: beginner

    331:   Concepts: Jacobi, diagonal scaling, preconditioners

    333:   Notes: By using KSPSetPreconditionerSide(ksp,PC_SYMMETRIC) or -ksp_symmetric_pc you
    334:          can scale each side of the matrix by the squareroot of the diagonal entries.

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

    338: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
    339:            PCJacobiSetUseRowMax(),
    340: M*/

    345: PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_Jacobi(PC pc)
    346: {
    347:   PC_Jacobi      *jac;

    351:   /*
    352:      Creates the private data structure for this preconditioner and
    353:      attach it to the PC object.
    354:   */
    355:   PetscNew(PC_Jacobi,&jac);
    356:   pc->data  = (void*)jac;

    358:   /*
    359:      Logs the memory usage; this is not needed but allows PETSc to
    360:      monitor how much memory is being used for various purposes.
    361:   */
    362:   PetscLogObjectMemory(pc,sizeof(PC_Jacobi));

    364:   /*
    365:      Initialize the pointers to vectors to ZERO; these will be used to store
    366:      diagonal entries of the matrix for fast preconditioner application.
    367:   */
    368:   jac->diag          = 0;
    369:   jac->diagsqrt      = 0;
    370:   jac->userowmax     = PETSC_FALSE;

    372:   /*
    373:       Set the pointers for the functions that are provided above.
    374:       Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
    375:       are called, they will automatically call these functions.  Note we
    376:       choose not to provide a couple of these functions since they are
    377:       not needed.
    378:   */
    379:   pc->ops->apply               = PCApply_Jacobi;
    380:   pc->ops->applytranspose      = PCApply_Jacobi;
    381:   pc->ops->setup               = PCSetUp_Jacobi;
    382:   pc->ops->destroy             = PCDestroy_Jacobi;
    383:   pc->ops->setfromoptions      = PCSetFromOptions_Jacobi;
    384:   pc->ops->view                = 0;
    385:   pc->ops->applyrichardson     = 0;
    386:   pc->ops->applysymmetricleft  = PCApplySymmetricLeftOrRight_Jacobi;
    387:   pc->ops->applysymmetricright = PCApplySymmetricLeftOrRight_Jacobi;
    388:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCJacobiSetUseRowMax_C","PCJacobiSetUseRowMax_Jacobi",
    389:                     PCJacobiSetUseRowMax_Jacobi);
    390:   return(0);
    391: }

    396: /*@
    397:    PCJacobiSetUseRowMax - Causes the Jacobi preconditioner to use the
    398:       maximum entry in each row as the diagonal preconditioner, instead of
    399:       the diagonal entry

    401:    Collective on PC

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

    407:    Options Database Key:
    408: .  -pc_jacobi_rowmax

    410:    Level: intermediate

    412:    Concepts: Jacobi preconditioner

    414: @*/
    415: PetscErrorCode PETSCKSP_DLLEXPORT PCJacobiSetUseRowMax(PC pc)
    416: {
    417:   PetscErrorCode ierr,(*f)(PC);

    421:   PetscObjectQueryFunction((PetscObject)pc,"PCJacobiSetUseRowMax_C",(void (**)(void))&f);
    422:   if (f) {
    423:     (*f)(pc);
    424:   }
    425:   return(0);
    426: }