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
58: of the preconditioner matrix */
59: Vec diagsqrt; /* vector containing the reciprocals of the square roots of
60: the diagonal elements of the preconditioner matrix (used
61: only for symmetric preconditioner application) */
62: PetscBool userowmax;
63: PetscBool userowsum;
64: PetscBool useabs; /* use the absolute values of the diagonal entries */
65: } PC_Jacobi;
67: EXTERN_C_BEGIN
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: }
79: EXTERN_C_END
81: EXTERN_C_BEGIN
84: PetscErrorCode PCJacobiSetUseRowSum_Jacobi(PC pc) 85: {
86: PC_Jacobi *j;
89: j = (PC_Jacobi*)pc->data;
90: j->userowsum = PETSC_TRUE;
91: return(0);
92: }
93: EXTERN_C_END
95: EXTERN_C_BEGIN
98: PetscErrorCode PCJacobiSetUseAbs_Jacobi(PC pc) 99: {
100: PC_Jacobi *j;
103: j = (PC_Jacobi*)pc->data;
104: j->useabs = PETSC_TRUE;
105: return(0);
106: }
107: EXTERN_C_END
109: /* -------------------------------------------------------------------------- */
110: /*
111: PCSetUp_Jacobi - Prepares for the use of the Jacobi preconditioner
112: by setting data structures and options.
114: Input Parameter:
115: . pc - the preconditioner context
117: Application Interface Routine: PCSetUp()
119: Notes:
120: The interface routine PCSetUp() is not usually called directly by
121: the user, but instead is called by PCApply() if necessary.
122: */
125: static PetscErrorCode PCSetUp_Jacobi(PC pc)126: {
127: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
128: Vec diag,diagsqrt;
130: PetscInt n,i;
131: PetscScalar *x;
132: PetscBool zeroflag = PETSC_FALSE;
135: /*
136: For most preconditioners the code would begin here something like
138: if (pc->setupcalled == 0) { allocate space the first time this is ever called
139: MatGetVecs(pc->mat,&jac->diag);
140: PetscLogObjectParent(pc,jac->diag);
141: }
143: But for this preconditioner we want to support use of both the matrix' diagonal
144: elements (for left or right preconditioning) and square root of diagonal elements
145: (for symmetric preconditioning). Hence we do not allocate space here, since we
146: don't know at this point which will be needed (diag and/or diagsqrt) until the user
147: applies the preconditioner, and we don't want to allocate BOTH unless we need
148: them both. Thus, the diag and diagsqrt are allocated in PCSetUp_Jacobi_NonSymmetric()
149: and PCSetUp_Jacobi_Symmetric(), respectively.
150: */
152: /*
153: Here we set up the preconditioner; that is, we copy the diagonal values from
154: the matrix and put them into a format to make them quick to apply as a preconditioner.
155: */
156: diag = jac->diag;
157: diagsqrt = jac->diagsqrt;
159: if (diag) {
160: if (jac->userowmax) {
161: MatGetRowMaxAbs(pc->pmat,diag,PETSC_NULL);
162: } else if (jac->userowsum) {
163: MatGetRowSum(pc->pmat,diag);
164: } else {
165: MatGetDiagonal(pc->pmat,diag);
166: }
167: VecReciprocal(diag);
168: VecGetLocalSize(diag,&n);
169: VecGetArray(diag,&x);
170: if (jac->useabs) {
171: for (i=0; i<n; i++) {
172: x[i] = PetscAbsScalar(x[i]);
173: }
174: }
175: for (i=0; i<n; i++) {
176: if (x[i] == 0.0) {
177: x[i] = 1.0;
178: zeroflag = PETSC_TRUE;
179: }
180: }
181: VecRestoreArray(diag,&x);
182: }
183: if (diagsqrt) {
184: if (jac->userowmax) {
185: MatGetRowMaxAbs(pc->pmat,diagsqrt,PETSC_NULL);
186: } else if (jac->userowsum) {
187: MatGetRowSum(pc->pmat,diagsqrt);
188: } else {
189: MatGetDiagonal(pc->pmat,diagsqrt);
190: }
191: VecGetLocalSize(diagsqrt,&n);
192: VecGetArray(diagsqrt,&x);
193: for (i=0; i<n; i++) {
194: if (x[i] != 0.0) x[i] = 1.0/PetscSqrtReal(PetscAbsScalar(x[i]));
195: else {
196: x[i] = 1.0;
197: zeroflag = PETSC_TRUE;
198: }
199: }
200: VecRestoreArray(diagsqrt,&x);
201: }
202: if (zeroflag) {
203: PetscInfo(pc,"Zero detected in diagonal of matrix, using 1 at those locations\n");
204: }
205: return(0);
206: }
207: /* -------------------------------------------------------------------------- */
208: /*
209: PCSetUp_Jacobi_Symmetric - Allocates the vector needed to store the
210: inverse of the square root of the diagonal entries of the matrix. This
211: is used for symmetric application of the Jacobi preconditioner.
213: Input Parameter:
214: . pc - the preconditioner context
215: */
218: static PetscErrorCode PCSetUp_Jacobi_Symmetric(PC pc)219: {
221: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
224: MatGetVecs(pc->pmat,&jac->diagsqrt,0);
225: PetscLogObjectParent(pc,jac->diagsqrt);
226: PCSetUp_Jacobi(pc);
227: return(0);
228: }
229: /* -------------------------------------------------------------------------- */
230: /*
231: PCSetUp_Jacobi_NonSymmetric - Allocates the vector needed to store the
232: inverse of the diagonal entries of the matrix. This is used for left of
233: right application of the Jacobi preconditioner.
235: Input Parameter:
236: . pc - the preconditioner context
237: */
240: static PetscErrorCode PCSetUp_Jacobi_NonSymmetric(PC pc)241: {
243: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
246: MatGetVecs(pc->pmat,&jac->diag,0);
247: PetscLogObjectParent(pc,jac->diag);
248: PCSetUp_Jacobi(pc);
249: return(0);
250: }
251: /* -------------------------------------------------------------------------- */
252: /*
253: PCApply_Jacobi - Applies the Jacobi preconditioner to a vector.
255: Input Parameters:
256: . pc - the preconditioner context
257: . x - input vector
259: Output Parameter:
260: . y - output vector
262: Application Interface Routine: PCApply()
263: */
266: static PetscErrorCode PCApply_Jacobi(PC pc,Vec x,Vec y)267: {
268: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
272: if (!jac->diag) {
273: PCSetUp_Jacobi_NonSymmetric(pc);
274: }
275: VecPointwiseMult(y,x,jac->diag);
276: return(0);
277: }
278: /* -------------------------------------------------------------------------- */
279: /*
280: PCApplySymmetricLeftOrRight_Jacobi - Applies the left or right part of a
281: symmetric preconditioner to a vector.
283: Input Parameters:
284: . pc - the preconditioner context
285: . x - input vector
287: Output Parameter:
288: . y - output vector
290: Application Interface Routines: PCApplySymmetricLeft(), PCApplySymmetricRight()
291: */
294: static PetscErrorCode PCApplySymmetricLeftOrRight_Jacobi(PC pc,Vec x,Vec y)295: {
297: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
300: if (!jac->diagsqrt) {
301: PCSetUp_Jacobi_Symmetric(pc);
302: }
303: VecPointwiseMult(y,x,jac->diagsqrt);
304: return(0);
305: }
306: /* -------------------------------------------------------------------------- */
309: static PetscErrorCode PCReset_Jacobi(PC pc)310: {
311: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
315: VecDestroy(&jac->diag);
316: VecDestroy(&jac->diagsqrt);
317: return(0);
318: }
320: /*
321: PCDestroy_Jacobi - Destroys the private context for the Jacobi preconditioner
322: that was created with PCCreate_Jacobi().
324: Input Parameter:
325: . pc - the preconditioner context
327: Application Interface Routine: PCDestroy()
328: */
331: static PetscErrorCode PCDestroy_Jacobi(PC pc)332: {
336: PCReset_Jacobi(pc);
338: /*
339: Free the private data structure that was hanging off the PC340: */
341: PetscFree(pc->data);
342: return(0);
343: }
347: static PetscErrorCode PCSetFromOptions_Jacobi(PC pc)348: {
349: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
353: PetscOptionsHead("Jacobi options");
354: PetscOptionsBool("-pc_jacobi_rowmax","Use row maximums for diagonal","PCJacobiSetUseRowMax",jac->userowmax,
355: &jac->userowmax,PETSC_NULL);
356: PetscOptionsBool("-pc_jacobi_rowsum","Use row sums for diagonal","PCJacobiSetUseRowSum",jac->userowsum,
357: &jac->userowsum,PETSC_NULL);
358: PetscOptionsBool("-pc_jacobi_abs","Use absolute values of diagaonal entries","PCJacobiSetUseAbs",jac->useabs,
359: &jac->useabs,PETSC_NULL);
360: PetscOptionsTail();
361: return(0);
362: }
364: /* -------------------------------------------------------------------------- */
365: /*
366: PCCreate_Jacobi - Creates a Jacobi preconditioner context, PC_Jacobi,
367: and sets this as the private data within the generic preconditioning
368: context, PC, that was created within PCCreate().
370: Input Parameter:
371: . pc - the preconditioner context
373: Application Interface Routine: PCCreate()
374: */
376: /*MC
377: PCJACOBI - Jacobi (i.e. diagonal scaling preconditioning)
379: Options Database Key:
380: + -pc_jacobi_rowmax - use the maximum absolute value in each row as the scaling factor,
381: rather than the diagonal
382: . -pc_jacobi_rowsum - use the row sums as the scaling factor,
383: rather than the diagonal
384: - -pc_jacobi_abs - use the absolute value of the diagaonl entry
386: Level: beginner
388: Concepts: Jacobi, diagonal scaling, preconditioners
390: Notes: By using KSPSetPCSide(ksp,PC_SYMMETRIC) or -ksp_pc_side symmetric
391: can scale each side of the matrix by the squareroot of the diagonal entries.
393: Zero entries along the diagonal are replaced with the value 1.0
395: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
396: PCJacobiSetUseRowMax(), PCJacobiSetUseRowSum(), PCJacobiSetUseAbs()
397: M*/
399: EXTERN_C_BEGIN
402: PetscErrorCode PCCreate_Jacobi(PC pc)403: {
404: PC_Jacobi *jac;
408: /*
409: Creates the private data structure for this preconditioner and
410: attach it to the PC object.
411: */
412: PetscNewLog(pc,PC_Jacobi,&jac);
413: pc->data = (void*)jac;
415: /*
416: Initialize the pointers to vectors to ZERO; these will be used to store
417: diagonal entries of the matrix for fast preconditioner application.
418: */
419: jac->diag = 0;
420: jac->diagsqrt = 0;
421: jac->userowmax = PETSC_FALSE;
422: jac->userowsum = PETSC_FALSE;
423: jac->useabs = PETSC_FALSE;
425: /*
426: Set the pointers for the functions that are provided above.
427: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
428: are called, they will automatically call these functions. Note we
429: choose not to provide a couple of these functions since they are
430: not needed.
431: */
432: pc->ops->apply = PCApply_Jacobi;
433: pc->ops->applytranspose = PCApply_Jacobi;
434: pc->ops->setup = PCSetUp_Jacobi;
435: pc->ops->reset = PCReset_Jacobi;
436: pc->ops->destroy = PCDestroy_Jacobi;
437: pc->ops->setfromoptions = PCSetFromOptions_Jacobi;
438: pc->ops->view = 0;
439: pc->ops->applyrichardson = 0;
440: pc->ops->applysymmetricleft = PCApplySymmetricLeftOrRight_Jacobi;
441: pc->ops->applysymmetricright = PCApplySymmetricLeftOrRight_Jacobi;
442: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCJacobiSetUseRowMax_C","PCJacobiSetUseRowMax_Jacobi",PCJacobiSetUseRowMax_Jacobi);
443: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCJacobiSetUseRowSum_C","PCJacobiSetUseRowSum_Jacobi",PCJacobiSetUseRowSum_Jacobi);
444: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCJacobiSetUseAbs_C","PCJacobiSetUseAbs_Jacobi",PCJacobiSetUseAbs_Jacobi);
445: return(0);
446: }
447: EXTERN_C_END
452: /*@
453: PCJacobiSetUseAbs - Causes the Jacobi preconditioner to use the
454: absolute value of the diagonal to for the preconditioner
456: Logically Collective on PC458: Input Parameters:
459: . pc - the preconditioner context
462: Options Database Key:
463: . -pc_jacobi_abs
465: Level: intermediate
467: Concepts: Jacobi preconditioner
469: .seealso: PCJacobiaUseRowMax(), PCJacobiaUseRowSum()
471: @*/
472: PetscErrorCodePCJacobiSetUseAbs(PC pc)473: {
478: PetscTryMethod(pc,"PCJacobiSetUseAbs_C",(PC),(pc));
479: return(0);
480: }
484: /*@
485: PCJacobiSetUseRowMax - Causes the Jacobi preconditioner to use the
486: maximum entry in each row as the diagonal preconditioner, instead of
487: the diagonal entry
489: Logically Collective on PC491: Input Parameters:
492: . pc - the preconditioner context
495: Options Database Key:
496: . -pc_jacobi_rowmax
498: Level: intermediate
500: Concepts: Jacobi preconditioner
502: .seealso: PCJacobiaUseAbs()
503: @*/
504: PetscErrorCodePCJacobiSetUseRowMax(PC pc)505: {
510: PetscTryMethod(pc,"PCJacobiSetUseRowMax_C",(PC),(pc));
511: return(0);
512: }
516: /*@
517: PCJacobiSetUseRowSum - Causes the Jacobi preconditioner to use the
518: sum of each row as the diagonal preconditioner, instead of
519: the diagonal entry
521: Logical Collective on PC523: Input Parameters:
524: . pc - the preconditioner context
527: Options Database Key:
528: . -pc_jacobi_rowsum
530: Level: intermediate
532: Concepts: Jacobi preconditioner
534: .seealso: PCJacobiaUseAbs(), PCJacobiaUseRowSum()
535: @*/
536: PetscErrorCodePCJacobiSetUseRowSum(PC pc)537: {
542: PetscTryMethod(pc,"PCJacobiSetUseRowSum_C",(PC),(pc));
543: return(0);
544: }