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(pc,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(pc,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(pc,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 PC331: */
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,PC_Jacobi,&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 PC447: 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: PetscErrorCodePCJacobiSetUseAbs(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 PC480: 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: PetscErrorCodePCJacobiSetUseRowMax(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 PC512: 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: PetscErrorCodePCJacobiSetUseRowSum(PC pc)526: {
531: PetscTryMethod(pc,"PCJacobiSetUseRowSum_C",(PC),(pc));
532: return(0);
533: }