Actual source code: jacobi.c

  1: #define PETSCKSP_DLL

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

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

 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 private/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:   PetscTruth useabs;             /* use the absolute values of the diagonal entries */
 65: } PC_Jacobi;

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

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

 84: PetscErrorCode  PCJacobiSetUseAbs_Jacobi(PC pc)
 85: {
 86:   PC_Jacobi *j;

 89:   j         = (PC_Jacobi*)pc->data;
 90:   j->useabs = PETSC_TRUE;
 91:   return(0);
 92: }

 95: /* -------------------------------------------------------------------------- */
 96: /*
 97:    PCSetUp_Jacobi - Prepares for the use of the Jacobi preconditioner
 98:                     by setting data structures and options.   

100:    Input Parameter:
101: .  pc - the preconditioner context

103:    Application Interface Routine: PCSetUp()

105:    Notes:
106:    The interface routine PCSetUp() is not usually called directly by
107:    the user, but instead is called by PCApply() if necessary.
108: */
111: static PetscErrorCode PCSetUp_Jacobi(PC pc)
112: {
113:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;
114:   Vec            diag,diagsqrt;
116:   PetscInt       n,i;
117:   PetscScalar    *x;
118:   PetscTruth     zeroflag = PETSC_FALSE;

121:   /*
122:        For most preconditioners the code would begin here something like

124:   if (pc->setupcalled == 0) { allocate space the first time this is ever called
125:     MatGetVecs(pc->mat,&jac->diag);
126:     PetscLogObjectParent(pc,jac->diag);
127:   }

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

138:   /*
139:     Here we set up the preconditioner; that is, we copy the diagonal values from
140:     the matrix and put them into a format to make them quick to apply as a preconditioner.
141:   */
142:   diag     = jac->diag;
143:   diagsqrt = jac->diagsqrt;

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

195:    Input Parameter:
196: .  pc - the preconditioner context
197: */
200: static PetscErrorCode PCSetUp_Jacobi_Symmetric(PC pc)
201: {
203:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

206:   MatGetVecs(pc->pmat,&jac->diagsqrt,0);
207:   PetscLogObjectParent(pc,jac->diagsqrt);
208:   PCSetUp_Jacobi(pc);
209:   return(0);
210: }
211: /* -------------------------------------------------------------------------- */
212: /*
213:    PCSetUp_Jacobi_NonSymmetric - Allocates the vector needed to store the
214:    inverse of the diagonal entries of the matrix.  This is used for left of
215:    right application of the Jacobi preconditioner.

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

228:   MatGetVecs(pc->pmat,&jac->diag,0);
229:   PetscLogObjectParent(pc,jac->diag);
230:   PCSetUp_Jacobi(pc);
231:   return(0);
232: }
233: /* -------------------------------------------------------------------------- */
234: /*
235:    PCApply_Jacobi - Applies the Jacobi preconditioner to a vector.

237:    Input Parameters:
238: .  pc - the preconditioner context
239: .  x - input vector

241:    Output Parameter:
242: .  y - output vector

244:    Application Interface Routine: PCApply()
245:  */
248: static PetscErrorCode PCApply_Jacobi(PC pc,Vec x,Vec y)
249: {
250:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

254:   if (!jac->diag) {
255:     PCSetUp_Jacobi_NonSymmetric(pc);
256:   }
257:   VecPointwiseMult(y,x,jac->diag);
258:   return(0);
259: }
260: /* -------------------------------------------------------------------------- */
261: /*
262:    PCApplySymmetricLeftOrRight_Jacobi - Applies the left or right part of a
263:    symmetric preconditioner to a vector.

265:    Input Parameters:
266: .  pc - the preconditioner context
267: .  x - input vector

269:    Output Parameter:
270: .  y - output vector

272:    Application Interface Routines: PCApplySymmetricLeft(), PCApplySymmetricRight()
273: */
276: static PetscErrorCode PCApplySymmetricLeftOrRight_Jacobi(PC pc,Vec x,Vec y)
277: {
279:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

282:   if (!jac->diagsqrt) {
283:     PCSetUp_Jacobi_Symmetric(pc);
284:   }
285:   VecPointwiseMult(y,x,jac->diagsqrt);
286:   return(0);
287: }
288: /* -------------------------------------------------------------------------- */
289: /*
290:    PCDestroy_Jacobi - Destroys the private context for the Jacobi preconditioner
291:    that was created with PCCreate_Jacobi().

293:    Input Parameter:
294: .  pc - the preconditioner context

296:    Application Interface Routine: PCDestroy()
297: */
300: static PetscErrorCode PCDestroy_Jacobi(PC pc)
301: {
302:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

306:   if (jac->diag)     {VecDestroy(jac->diag);}
307:   if (jac->diagsqrt) {VecDestroy(jac->diagsqrt);}

309:   /*
310:       Free the private data structure that was hanging off the PC
311:   */
312:   PetscFree(jac);
313:   return(0);
314: }

318: static PetscErrorCode PCSetFromOptions_Jacobi(PC pc)
319: {
320:   PC_Jacobi      *jac = (PC_Jacobi*)pc->data;

324:   PetscOptionsHead("Jacobi options");
325:     PetscOptionsTruth("-pc_jacobi_rowmax","Use row maximums for diagonal","PCJacobiSetUseRowMax",jac->userowmax,
326:                           &jac->userowmax,PETSC_NULL);
327:     PetscOptionsTruth("-pc_jacobi_abs","Use absolute values of diagaonal entries","PCJacobiSetUseAbs",jac->useabs,
328:                           &jac->useabs,PETSC_NULL);
329:   PetscOptionsTail();
330:   return(0);
331: }

333: /* -------------------------------------------------------------------------- */
334: /*
335:    PCCreate_Jacobi - Creates a Jacobi preconditioner context, PC_Jacobi, 
336:    and sets this as the private data within the generic preconditioning 
337:    context, PC, that was created within PCCreate().

339:    Input Parameter:
340: .  pc - the preconditioner context

342:    Application Interface Routine: PCCreate()
343: */

345: /*MC
346:      PCJACOBI - Jacobi (i.e. diagonal scaling preconditioning)

348:    Options Database Key:
349: +    -pc_jacobi_rowmax - use the maximum absolute value in each row as the scaling factor,
350:                         rather than the diagonal
351: -    -pc_jacobi_abs - use the absolute value of the diagaonl entry

353:    Level: beginner

355:   Concepts: Jacobi, diagonal scaling, preconditioners

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

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

362: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
363:            PCJacobiSetUseRowMax(), PCJacobiSetUseAbs()
364: M*/

369: PetscErrorCode  PCCreate_Jacobi(PC pc)
370: {
371:   PC_Jacobi      *jac;

375:   /*
376:      Creates the private data structure for this preconditioner and
377:      attach it to the PC object.
378:   */
379:   PetscNew(PC_Jacobi,&jac);
380:   pc->data  = (void*)jac;

382:   /*
383:      Logs the memory usage; this is not needed but allows PETSc to 
384:      monitor how much memory is being used for various purposes.
385:   */
386:   PetscLogObjectMemory(pc,sizeof(PC_Jacobi));

388:   /*
389:      Initialize the pointers to vectors to ZERO; these will be used to store
390:      diagonal entries of the matrix for fast preconditioner application.
391:   */
392:   jac->diag          = 0;
393:   jac->diagsqrt      = 0;
394:   jac->userowmax     = PETSC_FALSE;
395:   jac->useabs        = PETSC_FALSE;

397:   /*
398:       Set the pointers for the functions that are provided above.
399:       Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
400:       are called, they will automatically call these functions.  Note we
401:       choose not to provide a couple of these functions since they are
402:       not needed.
403:   */
404:   pc->ops->apply               = PCApply_Jacobi;
405:   pc->ops->applytranspose      = PCApply_Jacobi;
406:   pc->ops->setup               = PCSetUp_Jacobi;
407:   pc->ops->destroy             = PCDestroy_Jacobi;
408:   pc->ops->setfromoptions      = PCSetFromOptions_Jacobi;
409:   pc->ops->view                = 0;
410:   pc->ops->applyrichardson     = 0;
411:   pc->ops->applysymmetricleft  = PCApplySymmetricLeftOrRight_Jacobi;
412:   pc->ops->applysymmetricright = PCApplySymmetricLeftOrRight_Jacobi;
413:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCJacobiSetUseRowMax_C","PCJacobiSetUseRowMax_Jacobi",PCJacobiSetUseRowMax_Jacobi);
414:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCJacobiSetUseAbs_C","PCJacobiSetUseAbs_Jacobi",PCJacobiSetUseAbs_Jacobi);
415:   return(0);
416: }


422: /*@
423:    PCJacobiSetUseAbs - Causes the Jacobi preconditioner to use the 
424:       absolute value of the diagonal to for the preconditioner

426:    Collective on PC

428:    Input Parameters:
429: .  pc - the preconditioner context


432:    Options Database Key:
433: .  -pc_jacobi_abs

435:    Level: intermediate

437:    Concepts: Jacobi preconditioner

439: .seealso: PCJacobiaUseRowMax()

441: @*/
442: PetscErrorCode  PCJacobiSetUseAbs(PC pc)
443: {
444:   PetscErrorCode ierr,(*f)(PC);

448:   PetscObjectQueryFunction((PetscObject)pc,"PCJacobiSetUseAbs_C",(void (**)(void))&f);
449:   if (f) {
450:     (*f)(pc);
451:   }
452:   return(0);
453: }

457: /*@
458:    PCJacobiSetUseRowMax - Causes the Jacobi preconditioner to use the 
459:       maximum entry in each row as the diagonal preconditioner, instead of
460:       the diagonal entry

462:    Collective on PC

464:    Input Parameters:
465: .  pc - the preconditioner context


468:    Options Database Key:
469: .  -pc_jacobi_rowmax 

471:    Level: intermediate

473:    Concepts: Jacobi preconditioner

475: .seealso: PCJacobiaUseAbs()
476: @*/
477: PetscErrorCode  PCJacobiSetUseRowMax(PC pc)
478: {
479:   PetscErrorCode ierr,(*f)(PC);

483:   PetscObjectQueryFunction((PetscObject)pc,"PCJacobiSetUseRowMax_C",(void (**)(void))&f);
484:   if (f) {
485:     (*f)(pc);
486:   }
487:   return(0);
488: }