Actual source code: jacobi.c

  1: /*$Id: jacobi.c,v 1.75 2001/08/07 03:03:32 balay Exp $*/

  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/sles/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;

 66: EXTERN_C_BEGIN
 67: #undef __FUNCT__  
 69: int PCJacobiSetUseRowMax_Jacobi(PC pc)
 70: {
 71:   PC_Jacobi *j;

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

 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: */
 94: #undef __FUNCT__  
 96: static int PCSetUp_Jacobi(PC pc)
 97: {
 98:   PC_Jacobi     *jac = (PC_Jacobi*)pc->data;
 99:   Vec           diag,diagsqrt;
100:   int           ierr,n,i,zeroflag=0;
101:   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:     VecDuplicate(pc->vec,&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 locationsn");
165:   }

167:   return(0);
168: }
169: /* -------------------------------------------------------------------------- */
170: /*
171:    PCSetUp_Jacobi_Symmetric - Allocates the vector needed to store the
172:    inverse of the square root of the diagonal entries of the matrix.  This
173:    is used for symmetric application of the Jacobi preconditioner.

175:    Input Parameter:
176: .  pc - the preconditioner context
177: */
178: #undef __FUNCT__  
180: static int PCSetUp_Jacobi_Symmetric(PC pc)
181: {
182:   int        ierr;
183:   PC_Jacobi  *jac = (PC_Jacobi*)pc->data;


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

198:    Input Parameter:
199: .  pc - the preconditioner context
200: */
201: #undef __FUNCT__  
203: static int PCSetUp_Jacobi_NonSymmetric(PC pc)
204: {
205:   int        ierr;
206:   PC_Jacobi  *jac = (PC_Jacobi*)pc->data;


210:   VecDuplicate(pc->vec,&jac->diag);
211:   PetscLogObjectParent(pc,jac->diag);
212:   PCSetUp_Jacobi(pc);
213:   return(0);
214: }
215: /* -------------------------------------------------------------------------- */
216: /*
217:    PCApply_Jacobi - Applies the Jacobi preconditioner to a vector.

219:    Input Parameters:
220: .  pc - the preconditioner context
221: .  x - input vector

223:    Output Parameter:
224: .  y - output vector

226:    Application Interface Routine: PCApply()
227:  */
228: #undef __FUNCT__  
230: static int PCApply_Jacobi(PC pc,Vec x,Vec y)
231: {
232:   PC_Jacobi *jac = (PC_Jacobi*)pc->data;
233:   int       ierr;

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

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

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

254:    Application Interface Routines: PCApplySymmetricLeft(), PCApplySymmetricRight()
255: */
256: #undef __FUNCT__  
258: static int PCApplySymmetricLeftOrRight_Jacobi(PC pc,Vec x,Vec y)
259: {
260:   int       ierr;
261:   PC_Jacobi *jac = (PC_Jacobi*)pc->data;

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

275:    Input Parameter:
276: .  pc - the preconditioner context

278:    Application Interface Routine: PCDestroy()
279: */
280: #undef __FUNCT__  
282: static int PCDestroy_Jacobi(PC pc)
283: {
284:   PC_Jacobi *jac = (PC_Jacobi*)pc->data;
285:   int       ierr;

288:   if (jac->diag)     {VecDestroy(jac->diag);}
289:   if (jac->diagsqrt) {VecDestroy(jac->diagsqrt);}

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

298: #undef __FUNCT__  
300: static int PCSetFromOptions_Jacobi(PC pc)
301: {
302:   PC_Jacobi  *jac = (PC_Jacobi*)pc->data;
303:   int        ierr;

306:   PetscOptionsHead("Jacobi options");
307:     PetscOptionsLogical("-pc_jacobi_rowmax","Use row maximums for diagonal","PCJacobiSetUseRowMax",jac->userowmax,
308:                           &jac->userowmax,PETSC_NULL);
309:   PetscOptionsTail();
310:   return(0);
311: }

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

319:    Input Parameter:
320: .  pc - the preconditioner context

322:    Application Interface Routine: PCCreate()
323: */

325: /*MC
326:      PCJacobi - Jacobi (i.e. diagonal scaling preconditioning)

328:    Options Database Key:
329: .    -pc_jacobi_rowmax - use the maximum absolute value in each row as the scaling factor,
330:                         rather than the diagonal

332:    Level: beginner

334:   Concepts: Jacobi, diagonal scaling, preconditioners

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

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

341: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
342:            PCJacobiSetUseRowMax(), 
343: M*/

345: EXTERN_C_BEGIN
346: #undef __FUNCT__  
348: int PCCreate_Jacobi(PC pc)
349: {
350:   PC_Jacobi *jac;
351:   int       ierr;


355:   /*
356:      Creates the private data structure for this preconditioner and
357:      attach it to the PC object.
358:   */
359:   ierr      = PetscNew(PC_Jacobi,&jac);
360:   pc->data  = (void*)jac;

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

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

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

398: #undef __FUNCT__  
400: /*@
401:    PCJacobiSetUseRowMax - Causes the Jacobi preconditioner to use the 
402:       maximum entry in each row as the diagonal preconditioner, instead of
403:       the diagonal entry

405:    Collective on PC

407:    Input Parameters:
408: .  pc - the preconditioner context


411:    Options Database Key:
412: .  -pc_jacobi_rowmax 

414:    Level: intermediate

416:    Concepts: Jacobi preconditioner

418: @*/
419: int PCJacobiSetUseRowMax(PC pc)
420: {
421:   int ierr,(*f)(PC);

425:   PetscObjectQueryFunction((PetscObject)pc,"PCJacobiSetRowMax_C",(void (**)(void))&f);
426:   if (f) {
427:     (*f)(pc);
428:   }
429:   return(0);
430: }