Actual source code: jacobi.c
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>
53: const char *const PCJacobiTypes[] = {"DIAGONAL","ROWMAX","ROWSUM","PCJacobiType","PC_JACOBI_",NULL};
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: PetscBool fixdiag; /* fix zero diagonal terms */
67: } PC_Jacobi;
69: static PetscErrorCode PCJacobiSetType_Jacobi(PC pc,PCJacobiType type)
70: {
71: PC_Jacobi *j = (PC_Jacobi*)pc->data;
73: j->userowmax = PETSC_FALSE;
74: j->userowsum = PETSC_FALSE;
75: if (type == PC_JACOBI_ROWMAX) {
76: j->userowmax = PETSC_TRUE;
77: } else if (type == PC_JACOBI_ROWSUM) {
78: j->userowsum = PETSC_TRUE;
79: }
80: return 0;
81: }
83: static PetscErrorCode PCJacobiGetType_Jacobi(PC pc,PCJacobiType *type)
84: {
85: PC_Jacobi *j = (PC_Jacobi*)pc->data;
87: if (j->userowmax) {
88: *type = PC_JACOBI_ROWMAX;
89: } else if (j->userowsum) {
90: *type = PC_JACOBI_ROWSUM;
91: } else {
92: *type = PC_JACOBI_DIAGONAL;
93: }
94: return 0;
95: }
97: static PetscErrorCode PCJacobiSetUseAbs_Jacobi(PC pc,PetscBool flg)
98: {
99: PC_Jacobi *j = (PC_Jacobi*)pc->data;
101: j->useabs = flg;
102: return 0;
103: }
105: static PetscErrorCode PCJacobiGetUseAbs_Jacobi(PC pc,PetscBool *flg)
106: {
107: PC_Jacobi *j = (PC_Jacobi*)pc->data;
109: *flg = j->useabs;
110: return 0;
111: }
113: static PetscErrorCode PCJacobiSetFixDiagonal_Jacobi(PC pc,PetscBool flg)
114: {
115: PC_Jacobi *j = (PC_Jacobi*)pc->data;
117: j->fixdiag = flg;
118: return 0;
119: }
121: static PetscErrorCode PCJacobiGetFixDiagonal_Jacobi(PC pc,PetscBool *flg)
122: {
123: PC_Jacobi *j = (PC_Jacobi*)pc->data;
125: *flg = j->fixdiag;
126: return 0;
127: }
129: /* -------------------------------------------------------------------------- */
130: /*
131: PCSetUp_Jacobi - Prepares for the use of the Jacobi preconditioner
132: by setting data structures and options.
134: Input Parameter:
135: . pc - the preconditioner context
137: Application Interface Routine: PCSetUp()
139: Notes:
140: The interface routine PCSetUp() is not usually called directly by
141: the user, but instead is called by PCApply() if necessary.
142: */
143: static PetscErrorCode PCSetUp_Jacobi(PC pc)
144: {
145: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
146: Vec diag,diagsqrt;
147: PetscInt n,i;
148: PetscScalar *x;
149: PetscBool zeroflag = PETSC_FALSE;
151: /*
152: For most preconditioners the code would begin here something like
154: if (pc->setupcalled == 0) { allocate space the first time this is ever called
155: MatCreateVecs(pc->mat,&jac->diag);
156: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->diag);
157: }
159: But for this preconditioner we want to support use of both the matrix' diagonal
160: elements (for left or right preconditioning) and square root of diagonal elements
161: (for symmetric preconditioning). Hence we do not allocate space here, since we
162: don't know at this point which will be needed (diag and/or diagsqrt) until the user
163: applies the preconditioner, and we don't want to allocate BOTH unless we need
164: them both. Thus, the diag and diagsqrt are allocated in PCSetUp_Jacobi_NonSymmetric()
165: and PCSetUp_Jacobi_Symmetric(), respectively.
166: */
168: /*
169: Here we set up the preconditioner; that is, we copy the diagonal values from
170: the matrix and put them into a format to make them quick to apply as a preconditioner.
171: */
172: diag = jac->diag;
173: diagsqrt = jac->diagsqrt;
175: if (diag) {
176: PetscBool isspd;
178: if (jac->userowmax) {
179: MatGetRowMaxAbs(pc->pmat,diag,NULL);
180: } else if (jac->userowsum) {
181: MatGetRowSum(pc->pmat,diag);
182: } else {
183: MatGetDiagonal(pc->pmat,diag);
184: }
185: VecReciprocal(diag);
186: if (jac->useabs) {
187: VecAbs(diag);
188: }
189: MatGetOption(pc->pmat,MAT_SPD,&isspd);
190: if (jac->fixdiag && !isspd) {
191: VecGetLocalSize(diag,&n);
192: VecGetArray(diag,&x);
193: for (i=0; i<n; i++) {
194: if (x[i] == 0.0) {
195: x[i] = 1.0;
196: zeroflag = PETSC_TRUE;
197: }
198: }
199: VecRestoreArray(diag,&x);
200: }
201: }
202: if (diagsqrt) {
203: if (jac->userowmax) {
204: MatGetRowMaxAbs(pc->pmat,diagsqrt,NULL);
205: } else if (jac->userowsum) {
206: MatGetRowSum(pc->pmat,diagsqrt);
207: } else {
208: MatGetDiagonal(pc->pmat,diagsqrt);
209: }
210: VecGetLocalSize(diagsqrt,&n);
211: VecGetArray(diagsqrt,&x);
212: for (i=0; i<n; i++) {
213: if (x[i] != 0.0) x[i] = 1.0/PetscSqrtReal(PetscAbsScalar(x[i]));
214: else {
215: x[i] = 1.0;
216: zeroflag = PETSC_TRUE;
217: }
218: }
219: VecRestoreArray(diagsqrt,&x);
220: }
221: if (zeroflag) {
222: PetscInfo(pc,"Zero detected in diagonal of matrix, using 1 at those locations\n");
223: }
224: return 0;
225: }
226: /* -------------------------------------------------------------------------- */
227: /*
228: PCSetUp_Jacobi_Symmetric - Allocates the vector needed to store the
229: inverse of the square root of the diagonal entries of the matrix. This
230: is used for symmetric application of the Jacobi preconditioner.
232: Input Parameter:
233: . pc - the preconditioner context
234: */
235: static PetscErrorCode PCSetUp_Jacobi_Symmetric(PC pc)
236: {
237: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
239: MatCreateVecs(pc->pmat,&jac->diagsqrt,NULL);
240: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->diagsqrt);
241: PCSetUp_Jacobi(pc);
242: return 0;
243: }
244: /* -------------------------------------------------------------------------- */
245: /*
246: PCSetUp_Jacobi_NonSymmetric - Allocates the vector needed to store the
247: inverse of the diagonal entries of the matrix. This is used for left of
248: right application of the Jacobi preconditioner.
250: Input Parameter:
251: . pc - the preconditioner context
252: */
253: static PetscErrorCode PCSetUp_Jacobi_NonSymmetric(PC pc)
254: {
255: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
257: MatCreateVecs(pc->pmat,&jac->diag,NULL);
258: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->diag);
259: PCSetUp_Jacobi(pc);
260: return 0;
261: }
262: /* -------------------------------------------------------------------------- */
263: /*
264: PCApply_Jacobi - Applies the Jacobi preconditioner to a vector.
266: Input Parameters:
267: . pc - the preconditioner context
268: . x - input vector
270: Output Parameter:
271: . y - output vector
273: Application Interface Routine: PCApply()
274: */
275: static PetscErrorCode PCApply_Jacobi(PC pc,Vec x,Vec y)
276: {
277: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
279: if (!jac->diag) {
280: PCSetUp_Jacobi_NonSymmetric(pc);
281: }
282: VecPointwiseMult(y,x,jac->diag);
283: return 0;
284: }
285: /* -------------------------------------------------------------------------- */
286: /*
287: PCApplySymmetricLeftOrRight_Jacobi - Applies the left or right part of a
288: symmetric preconditioner to a vector.
290: Input Parameters:
291: . pc - the preconditioner context
292: . x - input vector
294: Output Parameter:
295: . y - output vector
297: Application Interface Routines: PCApplySymmetricLeft(), PCApplySymmetricRight()
298: */
299: static PetscErrorCode PCApplySymmetricLeftOrRight_Jacobi(PC pc,Vec x,Vec y)
300: {
301: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
303: if (!jac->diagsqrt) {
304: PCSetUp_Jacobi_Symmetric(pc);
305: }
306: VecPointwiseMult(y,x,jac->diagsqrt);
307: return 0;
308: }
310: /* -------------------------------------------------------------------------- */
311: static PetscErrorCode PCReset_Jacobi(PC pc)
312: {
313: 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: */
329: static PetscErrorCode PCDestroy_Jacobi(PC pc)
330: {
331: PCReset_Jacobi(pc);
332: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetType_C",NULL);
333: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiGetType_C",NULL);
334: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetUseAbs_C",NULL);
335: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiGetUseAbs_C",NULL);
336: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetFixDiagonal_C",NULL);
337: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiGetFixDiagonal_C",NULL);
339: /*
340: Free the private data structure that was hanging off the PC
341: */
342: PetscFree(pc->data);
343: return 0;
344: }
346: static PetscErrorCode PCSetFromOptions_Jacobi(PetscOptionItems *PetscOptionsObject,PC pc)
347: {
348: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
349: PetscBool flg;
350: PCJacobiType deflt,type;
352: PCJacobiGetType(pc,&deflt);
353: PetscOptionsHead(PetscOptionsObject,"Jacobi options");
354: PetscOptionsEnum("-pc_jacobi_type","How to construct diagonal matrix","PCJacobiSetType",PCJacobiTypes,(PetscEnum)deflt,(PetscEnum*)&type,&flg);
355: if (flg) {
356: PCJacobiSetType(pc,type);
357: }
358: PetscOptionsBool("-pc_jacobi_abs","Use absolute values of diagonal entries","PCJacobiSetUseAbs",jac->useabs,&jac->useabs,NULL);
359: PetscOptionsBool("-pc_jacobi_fixdiagonal","Fix null terms on diagonal","PCJacobiSetFixDiagonal",jac->fixdiag,&jac->fixdiag,NULL);
360: PetscOptionsTail();
361: return 0;
362: }
364: static PetscErrorCode PCView_Jacobi(PC pc, PetscViewer viewer)
365: {
366: PC_Jacobi *jac = (PC_Jacobi *) pc->data;
367: PetscBool iascii;
369: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
370: if (iascii) {
371: PCJacobiType type;
372: PetscBool useAbs,fixdiag;
373: PetscViewerFormat format;
375: PCJacobiGetType(pc, &type);
376: PCJacobiGetUseAbs(pc, &useAbs);
377: PCJacobiGetFixDiagonal(pc, &fixdiag);
378: PetscViewerASCIIPrintf(viewer, " type %s%s%s\n", PCJacobiTypes[type], useAbs ? ", using absolute value of entries" : "", !fixdiag ? ", not checking null diagonal entries" : "");
379: PetscViewerGetFormat(viewer, &format);
380: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
381: VecView(jac->diag, viewer);
382: }
383: }
384: return 0;
385: }
387: /* -------------------------------------------------------------------------- */
388: /*
389: PCCreate_Jacobi - Creates a Jacobi preconditioner context, PC_Jacobi,
390: and sets this as the private data within the generic preconditioning
391: context, PC, that was created within PCCreate().
393: Input Parameter:
394: . pc - the preconditioner context
396: Application Interface Routine: PCCreate()
397: */
399: /*MC
400: PCJACOBI - Jacobi (i.e. diagonal scaling preconditioning)
402: Options Database Key:
403: + -pc_jacobi_type <diagonal,rowmax,rowsum> - approach for forming the preconditioner
404: . -pc_jacobi_abs - use the absolute value of the diagonal entry
405: - -pc_jacobi_fixdiag - fix for zero diagonal terms by placing 1.0 in those locations
407: Level: beginner
409: Notes:
410: By using KSPSetPCSide(ksp,PC_SYMMETRIC) or -ksp_pc_side symmetric
411: can scale each side of the matrix by the square root of the diagonal entries.
413: Zero entries along the diagonal are replaced with the value 1.0
415: See PCPBJACOBI for a point-block Jacobi preconditioner
417: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
418: PCJacobiSetType(), PCJacobiSetUseAbs(), PCJacobiGetUseAbs(),
419: PCJacobiSetFixDiagonal(), PCJacobiGetFixDiagonal(), PCPBJACOBI
420: M*/
422: PETSC_EXTERN PetscErrorCode PCCreate_Jacobi(PC pc)
423: {
424: PC_Jacobi *jac;
426: /*
427: Creates the private data structure for this preconditioner and
428: attach it to the PC object.
429: */
430: PetscNewLog(pc,&jac);
431: pc->data = (void*)jac;
433: /*
434: Initialize the pointers to vectors to ZERO; these will be used to store
435: diagonal entries of the matrix for fast preconditioner application.
436: */
437: jac->diag = NULL;
438: jac->diagsqrt = NULL;
439: jac->userowmax = PETSC_FALSE;
440: jac->userowsum = PETSC_FALSE;
441: jac->useabs = PETSC_FALSE;
442: jac->fixdiag = PETSC_TRUE;
444: /*
445: Set the pointers for the functions that are provided above.
446: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
447: are called, they will automatically call these functions. Note we
448: choose not to provide a couple of these functions since they are
449: not needed.
450: */
451: pc->ops->apply = PCApply_Jacobi;
452: pc->ops->applytranspose = PCApply_Jacobi;
453: pc->ops->setup = PCSetUp_Jacobi;
454: pc->ops->reset = PCReset_Jacobi;
455: pc->ops->destroy = PCDestroy_Jacobi;
456: pc->ops->setfromoptions = PCSetFromOptions_Jacobi;
457: pc->ops->view = PCView_Jacobi;
458: pc->ops->applyrichardson = NULL;
459: pc->ops->applysymmetricleft = PCApplySymmetricLeftOrRight_Jacobi;
460: pc->ops->applysymmetricright = PCApplySymmetricLeftOrRight_Jacobi;
462: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetType_C",PCJacobiSetType_Jacobi);
463: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiGetType_C",PCJacobiGetType_Jacobi);
464: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetUseAbs_C",PCJacobiSetUseAbs_Jacobi);
465: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiGetUseAbs_C",PCJacobiGetUseAbs_Jacobi);
466: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetFixDiagonal_C",PCJacobiSetFixDiagonal_Jacobi);
467: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiGetFixDiagonal_C",PCJacobiGetFixDiagonal_Jacobi);
468: return 0;
469: }
471: /*@
472: PCJacobiSetUseAbs - Causes the Jacobi preconditioner to use the
473: absolute values of the diagonal divisors in the preconditioner
475: Logically Collective on PC
477: Input Parameters:
478: + pc - the preconditioner context
479: - flg - whether to use absolute values or not
481: Options Database Key:
482: . -pc_jacobi_abs <bool> - use absolute values
484: Notes:
485: This takes affect at the next construction of the preconditioner
487: Level: intermediate
489: .seealso: PCJacobiaSetType(), PCJacobiGetUseAbs()
491: @*/
492: PetscErrorCode PCJacobiSetUseAbs(PC pc,PetscBool flg)
493: {
495: PetscTryMethod(pc,"PCJacobiSetUseAbs_C",(PC,PetscBool),(pc,flg));
496: return 0;
497: }
499: /*@
500: PCJacobiGetUseAbs - Determines if the Jacobi preconditioner uses the
501: absolute values of the diagonal divisors in the preconditioner
503: Logically Collective on PC
505: Input Parameter:
506: . pc - the preconditioner context
508: Output Parameter:
509: . flg - whether to use absolute values or not
511: Level: intermediate
513: .seealso: PCJacobiaSetType(), PCJacobiSetUseAbs(), PCJacobiGetType()
515: @*/
516: PetscErrorCode PCJacobiGetUseAbs(PC pc,PetscBool *flg)
517: {
519: PetscUseMethod(pc,"PCJacobiGetUseAbs_C",(PC,PetscBool*),(pc,flg));
520: return 0;
521: }
523: /*@
524: PCJacobiSetFixDiagonal - Check for zero values on the diagonal and replace them with 1.0
526: Logically Collective on PC
528: Input Parameters:
529: + pc - the preconditioner context
530: - flg - the boolean flag
532: Options Database Key:
533: . -pc_jacobi_fixdiagonal <bool> - check for zero values on the diagonal
535: Notes:
536: This takes affect at the next construction of the preconditioner
538: Level: intermediate
540: .seealso: PCJacobiSetType(), PCJacobiGetFixDiagonal()
542: @*/
543: PetscErrorCode PCJacobiSetFixDiagonal(PC pc,PetscBool flg)
544: {
546: PetscTryMethod(pc,"PCJacobiSetFixDiagonal_C",(PC,PetscBool),(pc,flg));
547: return 0;
548: }
550: /*@
551: PCJacobiGetFixDiagonal - Determines if the Jacobi preconditioner checks for zero diagonal terms
553: Logically Collective on PC
555: Input Parameter:
556: . pc - the preconditioner context
558: Output Parameter:
559: . flg - the boolean flag
561: Options Database Key:
562: . -pc_jacobi_fixdiagonal <bool> - Fix 0 terms on diagonal by using 1
564: Level: intermediate
566: .seealso: PCJacobiSetType(), PCJacobiSetFixDiagonal()
568: @*/
569: PetscErrorCode PCJacobiGetFixDiagonal(PC pc,PetscBool *flg)
570: {
572: PetscUseMethod(pc,"PCJacobiGetFixDiagonal_C",(PC,PetscBool*),(pc,flg));
573: return 0;
574: }
576: /*@
577: PCJacobiSetType - Causes the Jacobi preconditioner to use either the diagonal, the maximum entry in each row,
578: of the sum of rows entries for the diagonal preconditioner
580: Logically Collective on PC
582: Input Parameters:
583: + pc - the preconditioner context
584: - type - PC_JACOBI_DIAGONAL, PC_JACOBI_ROWMAX, PC_JACOBI_ROWSUM
586: Options Database Key:
587: . -pc_jacobi_type <diagonal,rowmax,rowsum> - the type of diagonal matrix to use for Jacobi
589: Level: intermediate
591: .seealso: PCJacobiaUseAbs(), PCJacobiGetType()
592: @*/
593: PetscErrorCode PCJacobiSetType(PC pc,PCJacobiType type)
594: {
596: PetscTryMethod(pc,"PCJacobiSetType_C",(PC,PCJacobiType),(pc,type));
597: return 0;
598: }
600: /*@
601: PCJacobiGetType - Gets how the diagonal matrix is produced for the preconditioner
603: Not Collective on PC
605: Input Parameter:
606: . pc - the preconditioner context
608: Output Parameter:
609: . type - PC_JACOBI_DIAGONAL, PC_JACOBI_ROWMAX, PC_JACOBI_ROWSUM
611: Level: intermediate
613: .seealso: PCJacobiaUseAbs(), PCJacobiSetType()
614: @*/
615: PetscErrorCode PCJacobiGetType(PC pc,PCJacobiType *type)
616: {
618: PetscUseMethod(pc,"PCJacobiGetType_C",(PC,PCJacobiType*),(pc,type));
619: return 0;
620: }