Actual source code: jacobi.c

petsc-3.7.3 2016-08-01
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: 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;

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

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

 87: static PetscErrorCode  PCJacobiGetType_Jacobi(PC pc,PCJacobiType *type)
 88: {
 89:   PC_Jacobi *j = (PC_Jacobi*)pc->data;

 92:   if (j->userowmax) {
 93:     *type = PC_JACOBI_ROWMAX;
 94:   } else if (j->userowsum) {
 95:     *type = PC_JACOBI_ROWSUM;
 96:   } else {
 97:     *type = PC_JACOBI_DIAGONAL;
 98:   }
 99:   return(0);
100: }

104: static PetscErrorCode  PCJacobiSetUseAbs_Jacobi(PC pc,PetscBool flg)
105: {
106:   PC_Jacobi *j = (PC_Jacobi*)pc->data;

109:   j->useabs = flg;
110:   return(0);
111: }

115: static PetscErrorCode  PCJacobiGetUseAbs_Jacobi(PC pc,PetscBool *flg)
116: {
117:   PC_Jacobi *j = (PC_Jacobi*)pc->data;

120:   *flg = j->useabs;
121:   return(0);
122: }

124: /* -------------------------------------------------------------------------- */
125: /*
126:    PCSetUp_Jacobi - Prepares for the use of the Jacobi preconditioner
127:                     by setting data structures and options.

129:    Input Parameter:
130: .  pc - the preconditioner context

132:    Application Interface Routine: PCSetUp()

134:    Notes:
135:    The interface routine PCSetUp() is not usually called directly by
136:    the user, but instead is called by PCApply() if necessary.
137: */
140: static PetscErrorCode PCSetUp_Jacobi(PC pc)
141: {
142:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;
143:   Vec            diag,diagsqrt;
145:   PetscInt       n,i;
146:   PetscScalar    *x;
147:   PetscBool      zeroflag = PETSC_FALSE;

150:   /*
151:        For most preconditioners the code would begin here something like

153:   if (pc->setupcalled == 0) { allocate space the first time this is ever called
154:     MatCreateVecs(pc->mat,&jac->diag);
155:     PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->diag);
156:   }

158:     But for this preconditioner we want to support use of both the matrix' diagonal
159:     elements (for left or right preconditioning) and square root of diagonal elements
160:     (for symmetric preconditioning).  Hence we do not allocate space here, since we
161:     don't know at this point which will be needed (diag and/or diagsqrt) until the user
162:     applies the preconditioner, and we don't want to allocate BOTH unless we need
163:     them both.  Thus, the diag and diagsqrt are allocated in PCSetUp_Jacobi_NonSymmetric()
164:     and PCSetUp_Jacobi_Symmetric(), respectively.
165:   */

167:   /*
168:     Here we set up the preconditioner; that is, we copy the diagonal values from
169:     the matrix and put them into a format to make them quick to apply as a preconditioner.
170:   */
171:   diag     = jac->diag;
172:   diagsqrt = jac->diagsqrt;

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

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

237:   MatCreateVecs(pc->pmat,&jac->diagsqrt,0);
238:   PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->diagsqrt);
239:   PCSetUp_Jacobi(pc);
240:   return(0);
241: }
242: /* -------------------------------------------------------------------------- */
243: /*
244:    PCSetUp_Jacobi_NonSymmetric - Allocates the vector needed to store the
245:    inverse of the diagonal entries of the matrix.  This is used for left of
246:    right application of the Jacobi preconditioner.

248:    Input Parameter:
249: .  pc - the preconditioner context
250: */
253: static PetscErrorCode PCSetUp_Jacobi_NonSymmetric(PC pc)
254: {
256:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

259:   MatCreateVecs(pc->pmat,&jac->diag,0);
260:   PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->diag);
261:   PCSetUp_Jacobi(pc);
262:   return(0);
263: }
264: /* -------------------------------------------------------------------------- */
265: /*
266:    PCApply_Jacobi - Applies the Jacobi preconditioner to a vector.

268:    Input Parameters:
269: .  pc - the preconditioner context
270: .  x - input vector

272:    Output Parameter:
273: .  y - output vector

275:    Application Interface Routine: PCApply()
276:  */
279: static PetscErrorCode PCApply_Jacobi(PC pc,Vec x,Vec y)
280: {
281:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

285:   if (!jac->diag) {
286:     PCSetUp_Jacobi_NonSymmetric(pc);
287:   }
288:   VecPointwiseMult(y,x,jac->diag);
289:   return(0);
290: }
291: /* -------------------------------------------------------------------------- */
292: /*
293:    PCApplySymmetricLeftOrRight_Jacobi - Applies the left or right part of a
294:    symmetric preconditioner to a vector.

296:    Input Parameters:
297: .  pc - the preconditioner context
298: .  x - input vector

300:    Output Parameter:
301: .  y - output vector

303:    Application Interface Routines: PCApplySymmetricLeft(), PCApplySymmetricRight()
304: */
307: static PetscErrorCode PCApplySymmetricLeftOrRight_Jacobi(PC pc,Vec x,Vec y)
308: {
310:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

313:   if (!jac->diagsqrt) {
314:     PCSetUp_Jacobi_Symmetric(pc);
315:   }
316:   VecPointwiseMult(y,x,jac->diagsqrt);
317:   return(0);
318: }
319: /* -------------------------------------------------------------------------- */
322: static PetscErrorCode PCReset_Jacobi(PC pc)
323: {
324:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

328:   VecDestroy(&jac->diag);
329:   VecDestroy(&jac->diagsqrt);
330:   return(0);
331: }

333: /*
334:    PCDestroy_Jacobi - Destroys the private context for the Jacobi preconditioner
335:    that was created with PCCreate_Jacobi().

337:    Input Parameter:
338: .  pc - the preconditioner context

340:    Application Interface Routine: PCDestroy()
341: */
344: static PetscErrorCode PCDestroy_Jacobi(PC pc)
345: {

349:   PCReset_Jacobi(pc);

351:   /*
352:       Free the private data structure that was hanging off the PC
353:   */
354:   PetscFree(pc->data);
355:   return(0);
356: }

360: static PetscErrorCode PCSetFromOptions_Jacobi(PetscOptionItems *PetscOptionsObject,PC pc)
361: {
362:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;
364:   PetscBool      flg;
365:   PCJacobiType   deflt,type;

368:   PCJacobiGetType(pc,&deflt);
369:   PetscOptionsHead(PetscOptionsObject,"Jacobi options");
370:   PetscOptionsEnum("-pc_jacobi_type","How to construct diagonal matrix","PCJacobiSetType",PCJacobiTypes,(PetscEnum)deflt,(PetscEnum*)&type,&flg);
371:   if (flg) {
372:     PCJacobiSetType(pc,type);
373:   }
374:   PetscOptionsBool("-pc_jacobi_abs","Use absolute values of diagaonal entries","PCJacobiSetUseAbs",jac->useabs,&jac->useabs,NULL);
375:   PetscOptionsTail();
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:   Concepts: Jacobi, diagonal scaling, preconditioners

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

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

407:          See PCPBJACOBI for a point-block Jacobi preconditioner

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

415: PETSC_EXTERN PetscErrorCode PCCreate_Jacobi(PC pc)
416: {
417:   PC_Jacobi      *jac;

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

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

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

456:   PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetType_C",PCJacobiSetType_Jacobi);
457:   PetscObjectComposeFunction((PetscObject)pc,"PCJacobiGetType_C",PCJacobiGetType_Jacobi);
458:   PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetUseAbs_C",PCJacobiSetUseAbs_Jacobi);
459:   PetscObjectComposeFunction((PetscObject)pc,"PCJacobiGetUseAbs_C",PCJacobiGetUseAbs_Jacobi);
460:   return(0);
461: }

465: /*@
466:    PCJacobiSetUseAbs - Causes the Jacobi preconditioner to use the
467:       absolute values of the digonal divisors in the preconditioner

469:    Logically Collective on PC

471:    Input Parameters:
472: +  pc - the preconditioner context
473: -  flg - whether to use absolute values or not

475:    Options Database Key:
476: .  -pc_jacobi_abs

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

480:    Level: intermediate

482:    Concepts: Jacobi preconditioner

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

486: @*/
487: PetscErrorCode  PCJacobiSetUseAbs(PC pc,PetscBool flg)
488: {

493:   PetscTryMethod(pc,"PCJacobiSetUseAbs_C",(PC,PetscBool),(pc,flg));
494:   return(0);
495: }

499: /*@
500:    PCJacobiGetUseAbs - Determines if the Jacobi preconditioner uses the
501:       absolute values of the digonal divisors in the preconditioner

503:    Logically Collective on PC

505:    Input Parameter:
506: .  pc - the preconditioner context

508:    Output Parameter:
509: .  flg - whether to use absolute values or not

511:    Options Database Key:
512: .  -pc_jacobi_abs

514:    Level: intermediate

516:    Concepts: Jacobi preconditioner

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

520: @*/
521: PetscErrorCode  PCJacobiGetUseAbs(PC pc,PetscBool *flg)
522: {

527:   PetscUseMethod(pc,"PCJacobiGetUseAbs_C",(PC,PetscBool*),(pc,flg));
528:   return(0);
529: }

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

537:    Logically Collective on PC

539:    Input Parameters:
540: +  pc - the preconditioner context
541: -  type - PC_JACOBI_DIAGONAL, PC_JACOBI_ROWMAX, PC_JACOBI_ROWSUM

543:    Options Database Key:
544: .  -pc_jacobi_type <diagonal,rowmax,rowsum>

546:    Level: intermediate

548:    Concepts: Jacobi preconditioner

550: .seealso: PCJacobiaUseAbs(), PCJacobiGetType()
551: @*/
552: PetscErrorCode  PCJacobiSetType(PC pc,PCJacobiType type)
553: {

558:   PetscTryMethod(pc,"PCJacobiSetType_C",(PC,PCJacobiType),(pc,type));
559:   return(0);
560: }

564: /*@
565:    PCJacobiGetType - Gets how the diagonal matrix is produced for the preconditioner

567:    Not Collective on PC

569:    Input Parameter:
570: .  pc - the preconditioner context

572:    Output Parameter:
573: .  type - PC_JACOBI_DIAGONAL, PC_JACOBI_ROWMAX, PC_JACOBI_ROWSUM

575:    Level: intermediate

577:    Concepts: Jacobi preconditioner

579: .seealso: PCJacobiaUseAbs(), PCJacobiSetType()
580: @*/
581: PetscErrorCode  PCJacobiGetType(PC pc,PCJacobiType *type)
582: {

587:   PetscUseMethod(pc,"PCJacobiGetType_C",(PC,PCJacobiType*),(pc,type));
588:   return(0);
589: }