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 PC353: */
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(), PCPBJACOBI411: 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 PC471: 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: PetscErrorCodePCJacobiSetUseAbs(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 PC505: 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: PetscErrorCodePCJacobiGetUseAbs(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 PC539: 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: PetscErrorCodePCJacobiSetType(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 PC569: 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: PetscErrorCodePCJacobiGetType(PC pc,PCJacobiType *type)582: {
587: PetscUseMethod(pc,"PCJacobiGetType_C",(PC,PCJacobiType*),(pc,type));
588: return(0);
589: }