Actual source code: jacobi.c

petsc-3.8.2 2017-11-09
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>

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

 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:     VecGetArray(diag,&x);
175:     if (jac->useabs) {
176:       for (i=0; i<n; i++) x[i] = PetscAbsScalar(x[i]);
177:     }
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,0);
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,0);
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: /* -------------------------------------------------------------------------- */
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_type <diagonal,rowmax,rowsum> - approach for forming the preconditioner
372: -    -pc_jacobi_abs - use the absolute value of the diagonal entry

374:    Level: beginner

376:   Concepts: Jacobi, diagonal scaling, preconditioners

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

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

383:          See PCPBJACOBI for a point-block Jacobi preconditioner

385: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
386:            PCJacobiSetType(), PCJacobiSetUseAbs(), PCJacobiGetUseAbs(), PCPBJACOBI
387: M*/

389: PETSC_EXTERN PetscErrorCode PCCreate_Jacobi(PC pc)
390: {
391:   PC_Jacobi      *jac;

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

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

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

430:   PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetType_C",PCJacobiSetType_Jacobi);
431:   PetscObjectComposeFunction((PetscObject)pc,"PCJacobiGetType_C",PCJacobiGetType_Jacobi);
432:   PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetUseAbs_C",PCJacobiSetUseAbs_Jacobi);
433:   PetscObjectComposeFunction((PetscObject)pc,"PCJacobiGetUseAbs_C",PCJacobiGetUseAbs_Jacobi);
434:   return(0);
435: }

437: /*@
438:    PCJacobiSetUseAbs - Causes the Jacobi preconditioner to use the
439:       absolute values of the digonal divisors in the preconditioner

441:    Logically Collective on PC

443:    Input Parameters:
444: +  pc - the preconditioner context
445: -  flg - whether to use absolute values or not

447:    Options Database Key:
448: .  -pc_jacobi_abs

450:    Notes: This takes affect at the next construction of the preconditioner

452:    Level: intermediate

454:    Concepts: Jacobi preconditioner

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

458: @*/
459: PetscErrorCode  PCJacobiSetUseAbs(PC pc,PetscBool flg)
460: {

465:   PetscTryMethod(pc,"PCJacobiSetUseAbs_C",(PC,PetscBool),(pc,flg));
466:   return(0);
467: }

469: /*@
470:    PCJacobiGetUseAbs - Determines if the Jacobi preconditioner uses the
471:       absolute values of the digonal divisors in the preconditioner

473:    Logically Collective on PC

475:    Input Parameter:
476: .  pc - the preconditioner context

478:    Output Parameter:
479: .  flg - whether to use absolute values or not

481:    Options Database Key:
482: .  -pc_jacobi_abs

484:    Level: intermediate

486:    Concepts: Jacobi preconditioner

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

490: @*/
491: PetscErrorCode  PCJacobiGetUseAbs(PC pc,PetscBool *flg)
492: {

497:   PetscUseMethod(pc,"PCJacobiGetUseAbs_C",(PC,PetscBool*),(pc,flg));
498:   return(0);
499: }

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

505:    Logically Collective on PC

507:    Input Parameters:
508: +  pc - the preconditioner context
509: -  type - PC_JACOBI_DIAGONAL, PC_JACOBI_ROWMAX, PC_JACOBI_ROWSUM

511:    Options Database Key:
512: .  -pc_jacobi_type <diagonal,rowmax,rowsum>

514:    Level: intermediate

516:    Concepts: Jacobi preconditioner

518: .seealso: PCJacobiaUseAbs(), PCJacobiGetType()
519: @*/
520: PetscErrorCode  PCJacobiSetType(PC pc,PCJacobiType type)
521: {

526:   PetscTryMethod(pc,"PCJacobiSetType_C",(PC,PCJacobiType),(pc,type));
527:   return(0);
528: }

530: /*@
531:    PCJacobiGetType - Gets how the diagonal matrix is produced for the preconditioner

533:    Not Collective on PC

535:    Input Parameter:
536: .  pc - the preconditioner context

538:    Output Parameter:
539: .  type - PC_JACOBI_DIAGONAL, PC_JACOBI_ROWMAX, PC_JACOBI_ROWSUM

541:    Level: intermediate

543:    Concepts: Jacobi preconditioner

545: .seealso: PCJacobiaUseAbs(), PCJacobiSetType()
546: @*/
547: PetscErrorCode  PCJacobiGetType(PC pc,PCJacobiType *type)
548: {

553:   PetscUseMethod(pc,"PCJacobiGetType_C",(PC,PCJacobiType*),(pc,type));
554:   return(0);
555: }