Actual source code: jacobi.c

petsc-dev 2014-04-23
Report Typos and Errors
  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>   /*I "petscpc.h" I*/

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

 68: static PetscErrorCode  PCJacobiSetUseRowMax_Jacobi(PC pc)
 69: {
 70:   PC_Jacobi *j;

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

 80: static PetscErrorCode  PCJacobiSetUseRowSum_Jacobi(PC pc)
 81: {
 82:   PC_Jacobi *j;

 85:   j            = (PC_Jacobi*)pc->data;
 86:   j->userowsum = PETSC_TRUE;
 87:   return(0);
 88: }

 92: static PetscErrorCode  PCJacobiSetUseAbs_Jacobi(PC pc)
 93: {
 94:   PC_Jacobi *j;

 97:   j         = (PC_Jacobi*)pc->data;
 98:   j->useabs = PETSC_TRUE;
 99:   return(0);
100: }

102: /* -------------------------------------------------------------------------- */
103: /*
104:    PCSetUp_Jacobi - Prepares for the use of the Jacobi preconditioner
105:                     by setting data structures and options.

107:    Input Parameter:
108: .  pc - the preconditioner context

110:    Application Interface Routine: PCSetUp()

112:    Notes:
113:    The interface routine PCSetUp() is not usually called directly by
114:    the user, but instead is called by PCApply() if necessary.
115: */
118: static PetscErrorCode PCSetUp_Jacobi(PC pc)
119: {
120:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;
121:   Vec            diag,diagsqrt;
123:   PetscInt       n,i;
124:   PetscScalar    *x;
125:   PetscBool      zeroflag = PETSC_FALSE;

128:   /*
129:        For most preconditioners the code would begin here something like

131:   if (pc->setupcalled == 0) { allocate space the first time this is ever called
132:     MatGetVecs(pc->mat,&jac->diag);
133:     PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->diag);
134:   }

136:     But for this preconditioner we want to support use of both the matrix' diagonal
137:     elements (for left or right preconditioning) and square root of diagonal elements
138:     (for symmetric preconditioning).  Hence we do not allocate space here, since we
139:     don't know at this point which will be needed (diag and/or diagsqrt) until the user
140:     applies the preconditioner, and we don't want to allocate BOTH unless we need
141:     them both.  Thus, the diag and diagsqrt are allocated in PCSetUp_Jacobi_NonSymmetric()
142:     and PCSetUp_Jacobi_Symmetric(), respectively.
143:   */

145:   /*
146:     Here we set up the preconditioner; that is, we copy the diagonal values from
147:     the matrix and put them into a format to make them quick to apply as a preconditioner.
148:   */
149:   diag     = jac->diag;
150:   diagsqrt = jac->diagsqrt;

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

204:    Input Parameter:
205: .  pc - the preconditioner context
206: */
209: static PetscErrorCode PCSetUp_Jacobi_Symmetric(PC pc)
210: {
212:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

215:   MatGetVecs(pc->pmat,&jac->diagsqrt,0);
216:   PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->diagsqrt);
217:   PCSetUp_Jacobi(pc);
218:   return(0);
219: }
220: /* -------------------------------------------------------------------------- */
221: /*
222:    PCSetUp_Jacobi_NonSymmetric - Allocates the vector needed to store the
223:    inverse of the diagonal entries of the matrix.  This is used for left of
224:    right application of the Jacobi preconditioner.

226:    Input Parameter:
227: .  pc - the preconditioner context
228: */
231: static PetscErrorCode PCSetUp_Jacobi_NonSymmetric(PC pc)
232: {
234:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

237:   MatGetVecs(pc->pmat,&jac->diag,0);
238:   PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->diag);
239:   PCSetUp_Jacobi(pc);
240:   return(0);
241: }
242: /* -------------------------------------------------------------------------- */
243: /*
244:    PCApply_Jacobi - Applies the Jacobi preconditioner to a vector.

246:    Input Parameters:
247: .  pc - the preconditioner context
248: .  x - input vector

250:    Output Parameter:
251: .  y - output vector

253:    Application Interface Routine: PCApply()
254:  */
257: static PetscErrorCode PCApply_Jacobi(PC pc,Vec x,Vec y)
258: {
259:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

263:   if (!jac->diag) {
264:     PCSetUp_Jacobi_NonSymmetric(pc);
265:   }
266:   VecPointwiseMult(y,x,jac->diag);
267:   return(0);
268: }
269: /* -------------------------------------------------------------------------- */
270: /*
271:    PCApplySymmetricLeftOrRight_Jacobi - Applies the left or right part of a
272:    symmetric preconditioner to a vector.

274:    Input Parameters:
275: .  pc - the preconditioner context
276: .  x - input vector

278:    Output Parameter:
279: .  y - output vector

281:    Application Interface Routines: PCApplySymmetricLeft(), PCApplySymmetricRight()
282: */
285: static PetscErrorCode PCApplySymmetricLeftOrRight_Jacobi(PC pc,Vec x,Vec y)
286: {
288:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

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

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

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

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

318:    Application Interface Routine: PCDestroy()
319: */
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: }

338: static PetscErrorCode PCSetFromOptions_Jacobi(PC pc)
339: {
340:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

344:   PetscOptionsHead("Jacobi options");
345:   PetscOptionsBool("-pc_jacobi_rowmax","Use row maximums for diagonal","PCJacobiSetUseRowMax",jac->userowmax,
346:                           &jac->userowmax,NULL);
347:   PetscOptionsBool("-pc_jacobi_rowsum","Use row sums for diagonal","PCJacobiSetUseRowSum",jac->userowsum,
348:                           &jac->userowsum,NULL);
349:   PetscOptionsBool("-pc_jacobi_abs","Use absolute values of diagaonal entries","PCJacobiSetUseAbs",jac->useabs,
350:                           &jac->useabs,NULL);
351:   PetscOptionsTail();
352:   return(0);
353: }

355: /* -------------------------------------------------------------------------- */
356: /*
357:    PCCreate_Jacobi - Creates a Jacobi preconditioner context, PC_Jacobi,
358:    and sets this as the private data within the generic preconditioning
359:    context, PC, that was created within PCCreate().

361:    Input Parameter:
362: .  pc - the preconditioner context

364:    Application Interface Routine: PCCreate()
365: */

367: /*MC
368:      PCJACOBI - Jacobi (i.e. diagonal scaling preconditioning)

370:    Options Database Key:
371: +    -pc_jacobi_rowmax - use the maximum absolute value in each row as the scaling factor,
372:                         rather than the diagonal
373: .    -pc_jacobi_rowsum - use the row sums as the scaling factor,
374:                         rather than the diagonal
375: -    -pc_jacobi_abs - use the absolute value of the diagaonl entry

377:    Level: beginner

379:   Concepts: Jacobi, diagonal scaling, preconditioners

381:   Notes: By using KSPSetPCSide(ksp,PC_SYMMETRIC) or -ksp_pc_side symmetric
382:          can scale each side of the matrix by the squareroot of the diagonal entries.

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

386: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
387:            PCJacobiSetUseRowMax(), PCJacobiSetUseRowSum(), PCJacobiSetUseAbs()
388: M*/

392: PETSC_EXTERN PetscErrorCode PCCreate_Jacobi(PC pc)
393: {
394:   PC_Jacobi      *jac;

398:   /*
399:      Creates the private data structure for this preconditioner and
400:      attach it to the PC object.
401:   */
402:   PetscNewLog(pc,&jac);
403:   pc->data = (void*)jac;

405:   /*
406:      Initialize the pointers to vectors to ZERO; these will be used to store
407:      diagonal entries of the matrix for fast preconditioner application.
408:   */
409:   jac->diag      = 0;
410:   jac->diagsqrt  = 0;
411:   jac->userowmax = PETSC_FALSE;
412:   jac->userowsum = PETSC_FALSE;
413:   jac->useabs    = PETSC_FALSE;

415:   /*
416:       Set the pointers for the functions that are provided above.
417:       Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
418:       are called, they will automatically call these functions.  Note we
419:       choose not to provide a couple of these functions since they are
420:       not needed.
421:   */
422:   pc->ops->apply               = PCApply_Jacobi;
423:   pc->ops->applytranspose      = PCApply_Jacobi;
424:   pc->ops->setup               = PCSetUp_Jacobi;
425:   pc->ops->reset               = PCReset_Jacobi;
426:   pc->ops->destroy             = PCDestroy_Jacobi;
427:   pc->ops->setfromoptions      = PCSetFromOptions_Jacobi;
428:   pc->ops->view                = 0;
429:   pc->ops->applyrichardson     = 0;
430:   pc->ops->applysymmetricleft  = PCApplySymmetricLeftOrRight_Jacobi;
431:   pc->ops->applysymmetricright = PCApplySymmetricLeftOrRight_Jacobi;

433:   PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetUseRowMax_C",PCJacobiSetUseRowMax_Jacobi);
434:   PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetUseRowSum_C",PCJacobiSetUseRowSum_Jacobi);
435:   PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetUseAbs_C",PCJacobiSetUseAbs_Jacobi);
436:   return(0);
437: }

441: /*@
442:    PCJacobiSetUseAbs - Causes the Jacobi preconditioner to use the
443:       absolute value of the diagonal to for the preconditioner

445:    Logically Collective on PC

447:    Input Parameters:
448: .  pc - the preconditioner context


451:    Options Database Key:
452: .  -pc_jacobi_abs

454:    Level: intermediate

456:    Concepts: Jacobi preconditioner

458: .seealso: PCJacobiaUseRowMax(), PCJacobiaUseRowSum()

460: @*/
461: PetscErrorCode  PCJacobiSetUseAbs(PC pc)
462: {

467:   PetscTryMethod(pc,"PCJacobiSetUseAbs_C",(PC),(pc));
468:   return(0);
469: }

473: /*@
474:    PCJacobiSetUseRowMax - Causes the Jacobi preconditioner to use the
475:       maximum entry in each row as the diagonal preconditioner, instead of
476:       the diagonal entry

478:    Logically Collective on PC

480:    Input Parameters:
481: .  pc - the preconditioner context


484:    Options Database Key:
485: .  -pc_jacobi_rowmax

487:    Level: intermediate

489:    Concepts: Jacobi preconditioner

491: .seealso: PCJacobiaUseAbs()
492: @*/
493: PetscErrorCode  PCJacobiSetUseRowMax(PC pc)
494: {

499:   PetscTryMethod(pc,"PCJacobiSetUseRowMax_C",(PC),(pc));
500:   return(0);
501: }

505: /*@
506:    PCJacobiSetUseRowSum - Causes the Jacobi preconditioner to use the
507:       sum of each row as the diagonal preconditioner, instead of
508:       the diagonal entry

510:    Logical Collective on PC

512:    Input Parameters:
513: .  pc - the preconditioner context


516:    Options Database Key:
517: .  -pc_jacobi_rowsum

519:    Level: intermediate

521:    Concepts: Jacobi preconditioner

523: .seealso: PCJacobiaUseAbs(), PCJacobiaUseRowSum()
524: @*/
525: PetscErrorCode  PCJacobiSetUseRowSum(PC pc)
526: {

531:   PetscTryMethod(pc,"PCJacobiSetUseRowSum_C",(PC),(pc));
532:   return(0);
533: }