Actual source code: precon.c
2: /*
3: The PC (preconditioner) interface routines, callable by users.
4: */
5: #include <petsc/private/pcimpl.h>
6: #include <petscdm.h>
8: /* Logging support */
9: PetscClassId PC_CLASSID;
10: PetscLogEvent PC_SetUp, PC_SetUpOnBlocks, PC_Apply, PC_MatApply, PC_ApplyCoarse, PC_ApplyMultiple, PC_ApplySymmetricLeft;
11: PetscLogEvent PC_ApplySymmetricRight, PC_ModifySubMatrices, PC_ApplyOnBlocks, PC_ApplyTransposeOnBlocks;
12: PetscInt PetscMGLevelId;
13: PetscLogStage PCMPIStage;
15: PetscErrorCode PCGetDefaultType_Private(PC pc, const char *type[])
16: {
17: PetscMPIInt size;
18: PetscBool hasop, flg1, flg2, set, flg3, isnormal;
20: MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
21: if (pc->pmat) {
22: MatHasOperation(pc->pmat, MATOP_GET_DIAGONAL_BLOCK, &hasop);
23: if (size == 1) {
24: MatGetFactorAvailable(pc->pmat, "petsc", MAT_FACTOR_ICC, &flg1);
25: MatGetFactorAvailable(pc->pmat, "petsc", MAT_FACTOR_ILU, &flg2);
26: MatIsSymmetricKnown(pc->pmat, &set, &flg3);
27: PetscObjectTypeCompareAny((PetscObject)pc->pmat, &isnormal, MATNORMAL, MATNORMALHERMITIAN, NULL);
28: if (flg1 && (!flg2 || (set && flg3))) {
29: *type = PCICC;
30: } else if (flg2) {
31: *type = PCILU;
32: } else if (isnormal) {
33: *type = PCNONE;
34: } else if (hasop) { /* likely is a parallel matrix run on one processor */
35: *type = PCBJACOBI;
36: } else {
37: *type = PCNONE;
38: }
39: } else {
40: if (hasop) {
41: *type = PCBJACOBI;
42: } else {
43: *type = PCNONE;
44: }
45: }
46: } else {
47: if (size == 1) {
48: *type = PCILU;
49: } else {
50: *type = PCBJACOBI;
51: }
52: }
53: return 0;
54: }
56: /*@
57: PCReset - Resets a PC context to the pcsetupcalled = 0 state and removes any allocated `Vec`s and `Mat`s
59: Collective
61: Input Parameter:
62: . pc - the preconditioner context
64: Level: developer
66: Note:
67: This allows a `PC` to be reused for a different sized linear system but using the same options that have been previously set in the PC
69: .seealso: `PC`, `PCCreate()`, `PCSetUp()`
70: @*/
71: PetscErrorCode PCReset(PC pc)
72: {
74: PetscTryTypeMethod(pc, reset);
75: VecDestroy(&pc->diagonalscaleright);
76: VecDestroy(&pc->diagonalscaleleft);
77: MatDestroy(&pc->pmat);
78: MatDestroy(&pc->mat);
80: pc->setupcalled = 0;
81: return 0;
82: }
84: /*@C
85: PCDestroy - Destroys `PC` context that was created with `PCCreate()`.
87: Collective
89: Input Parameter:
90: . pc - the preconditioner context
92: Level: developer
94: .seealso: `PC`, `PCCreate()`, `PCSetUp()`
95: @*/
96: PetscErrorCode PCDestroy(PC *pc)
97: {
98: if (!*pc) return 0;
100: if (--((PetscObject)(*pc))->refct > 0) {
101: *pc = NULL;
102: return 0;
103: }
105: PCReset(*pc);
107: /* if memory was published with SAWs then destroy it */
108: PetscObjectSAWsViewOff((PetscObject)*pc);
109: PetscTryTypeMethod((*pc), destroy);
110: DMDestroy(&(*pc)->dm);
111: PetscHeaderDestroy(pc);
112: return 0;
113: }
115: /*@C
116: PCGetDiagonalScale - Indicates if the preconditioner applies an additional left and right
117: scaling as needed by certain time-stepping codes.
119: Logically Collective
121: Input Parameter:
122: . pc - the preconditioner context
124: Output Parameter:
125: . flag - `PETSC_TRUE` if it applies the scaling
127: Level: developer
129: Note:
130: If this returns `PETSC_TRUE` then the system solved via the Krylov method is
131: .vb
132: D M A D^{-1} y = D M b for left preconditioning or
133: D A M D^{-1} z = D b for right preconditioning
134: .ve
136: .seealso: `PC`, `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCDiagonalScaleRight()`, `PCSetDiagonalScale()`
137: @*/
138: PetscErrorCode PCGetDiagonalScale(PC pc, PetscBool *flag)
139: {
142: *flag = pc->diagonalscale;
143: return 0;
144: }
146: /*@
147: PCSetDiagonalScale - Indicates the left scaling to use to apply an additional left and right
148: scaling as needed by certain time-stepping codes.
150: Logically Collective
152: Input Parameters:
153: + pc - the preconditioner context
154: - s - scaling vector
156: Level: intermediate
158: Notes:
159: The system solved via the Krylov method is
160: .vb
161: D M A D^{-1} y = D M b for left preconditioning or
162: D A M D^{-1} z = D b for right preconditioning
163: .ve
165: `PCDiagonalScaleLeft()` scales a vector by D. `PCDiagonalScaleRight()` scales a vector by D^{-1}.
167: .seealso: `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCDiagonalScaleRight()`, `PCGetDiagonalScale()`
168: @*/
169: PetscErrorCode PCSetDiagonalScale(PC pc, Vec s)
170: {
173: pc->diagonalscale = PETSC_TRUE;
175: PetscObjectReference((PetscObject)s);
176: VecDestroy(&pc->diagonalscaleleft);
178: pc->diagonalscaleleft = s;
180: VecDuplicate(s, &pc->diagonalscaleright);
181: VecCopy(s, pc->diagonalscaleright);
182: VecReciprocal(pc->diagonalscaleright);
183: return 0;
184: }
186: /*@
187: PCDiagonalScaleLeft - Scales a vector by the left scaling as needed by certain time-stepping codes.
189: Logically Collective
191: Input Parameters:
192: + pc - the preconditioner context
193: . in - input vector
194: - out - scaled vector (maybe the same as in)
196: Level: intermediate
198: Notes:
199: The system solved via the Krylov method is
200: .vb
201: D M A D^{-1} y = D M b for left preconditioning or
202: D A M D^{-1} z = D b for right preconditioning
203: .ve
205: `PCDiagonalScaleLeft()` scales a vector by D. `PCDiagonalScaleRight()` scales a vector by D^{-1}.
207: If diagonal scaling is turned off and in is not out then in is copied to out
209: .seealso: `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleSet()`, `PCDiagonalScaleRight()`, `PCDiagonalScale()`
210: @*/
211: PetscErrorCode PCDiagonalScaleLeft(PC pc, Vec in, Vec out)
212: {
216: if (pc->diagonalscale) {
217: VecPointwiseMult(out, pc->diagonalscaleleft, in);
218: } else if (in != out) {
219: VecCopy(in, out);
220: }
221: return 0;
222: }
224: /*@
225: PCDiagonalScaleRight - Scales a vector by the right scaling as needed by certain time-stepping codes.
227: Logically Collective
229: Input Parameters:
230: + pc - the preconditioner context
231: . in - input vector
232: - out - scaled vector (maybe the same as in)
234: Level: intermediate
236: Notes:
237: The system solved via the Krylov method is
238: .vb
239: D M A D^{-1} y = D M b for left preconditioning or
240: D A M D^{-1} z = D b for right preconditioning
241: .ve
243: `PCDiagonalScaleLeft()` scales a vector by D. `PCDiagonalScaleRight()` scales a vector by D^{-1}.
245: If diagonal scaling is turned off and in is not out then in is copied to out
247: .seealso: `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCDiagonalScaleSet()`, `PCDiagonalScale()`
248: @*/
249: PetscErrorCode PCDiagonalScaleRight(PC pc, Vec in, Vec out)
250: {
254: if (pc->diagonalscale) {
255: VecPointwiseMult(out, pc->diagonalscaleright, in);
256: } else if (in != out) {
257: VecCopy(in, out);
258: }
259: return 0;
260: }
262: /*@
263: PCSetUseAmat - Sets a flag to indicate that when the preconditioner needs to apply (part of) the
264: operator during the preconditioning process it applies the Amat provided to `TSSetRHSJacobian()`,
265: `TSSetIJacobian()`, `SNESSetJacobian()`, `KSPSetOperators()` or `PCSetOperators()` not the Pmat.
267: Logically Collective
269: Input Parameters:
270: + pc - the preconditioner context
271: - flg - `PETSC_TRUE` to use the Amat, `PETSC_FALSE` to use the Pmat (default is false)
273: Options Database Key:
274: . -pc_use_amat <true,false> - use the amat to apply the operator
276: Note:
277: For the common case in which the linear system matrix and the matrix used to construct the
278: preconditioner are identical, this routine is does nothing.
280: Level: intermediate
282: .seealso: `PC`, `PCGetUseAmat()`, `PCBJACOBI`, `PGMG`, `PCFIELDSPLIT`, `PCCOMPOSITE`
283: @*/
284: PetscErrorCode PCSetUseAmat(PC pc, PetscBool flg)
285: {
287: pc->useAmat = flg;
288: return 0;
289: }
291: /*@
292: PCSetErrorIfFailure - Causes `PC` to generate an error if a FPE, for example a zero pivot, is detected.
294: Logically Collective
296: Input Parameters:
297: + pc - iterative context obtained from PCCreate()
298: - flg - `PETSC_TRUE` indicates you want the error generated
300: Level: advanced
302: Notes:
303: Normally PETSc continues if a linear solver fails due to a failed setup of a preconditioner, you can call `KSPGetConvergedReason()` after a `KSPSolve()`
304: to determine if it has converged or failed. Or use -ksp_error_if_not_converged to cause the program to terminate as soon as lack of convergence is
305: detected.
307: This is propagated into KSPs used by this PC, which then propagate it into PCs used by those KSPs
309: .seealso: `PC`, `KSPSetErrorIfNotConverged()`, `PCGetInitialGuessNonzero()`, `PCSetInitialGuessKnoll()`, `PCGetInitialGuessKnoll()`
310: @*/
311: PetscErrorCode PCSetErrorIfFailure(PC pc, PetscBool flg)
312: {
315: pc->erroriffailure = flg;
316: return 0;
317: }
319: /*@
320: PCGetUseAmat - Gets a flag to indicate that when the preconditioner needs to apply (part of) the
321: operator during the preconditioning process it applies the Amat provided to `TSSetRHSJacobian()`,
322: `TSSetIJacobian()`, `SNESSetJacobian()`, `KSPSetOperators()` or `PCSetOperators()` not the Pmat.
324: Logically Collective
326: Input Parameter:
327: . pc - the preconditioner context
329: Output Parameter:
330: . flg - `PETSC_TRUE` to use the Amat, `PETSC_FALSE` to use the Pmat (default is false)
332: Note:
333: For the common case in which the linear system matrix and the matrix used to construct the
334: preconditioner are identical, this routine is does nothing.
336: Level: intermediate
338: .seealso: `PC`, `PCSetUseAmat()`, `PCBJACOBI`, `PGMG`, `PCFIELDSPLIT`, `PCCOMPOSITE`
339: @*/
340: PetscErrorCode PCGetUseAmat(PC pc, PetscBool *flg)
341: {
343: *flg = pc->useAmat;
344: return 0;
345: }
347: /*@
348: PCCreate - Creates a preconditioner context, `PC`
350: Collective
352: Input Parameter:
353: . comm - MPI communicator
355: Output Parameter:
356: . pc - location to put the preconditioner context
358: Note:
359: The default preconditioner for sparse matrices is `PCILU` or `PCICC` with 0 fill on one process and block Jacobi (`PCBJACOBI`) with `PCILU` or `PCICC`
360: in parallel. For dense matrices it is always `PCNONE`.
362: Level: developer
364: .seealso: `PC`, `PCSetUp()`, `PCApply()`, `PCDestroy()`
365: @*/
366: PetscErrorCode PCCreate(MPI_Comm comm, PC *newpc)
367: {
368: PC pc;
371: *newpc = NULL;
372: PCInitializePackage();
374: PetscHeaderCreate(pc, PC_CLASSID, "PC", "Preconditioner", "PC", comm, PCDestroy, PCView);
376: pc->mat = NULL;
377: pc->pmat = NULL;
378: pc->setupcalled = 0;
379: pc->setfromoptionscalled = 0;
380: pc->data = NULL;
381: pc->diagonalscale = PETSC_FALSE;
382: pc->diagonalscaleleft = NULL;
383: pc->diagonalscaleright = NULL;
385: pc->modifysubmatrices = NULL;
386: pc->modifysubmatricesP = NULL;
388: *newpc = pc;
389: return 0;
390: }
392: /*@
393: PCApply - Applies the preconditioner to a vector.
395: Collective
397: Input Parameters:
398: + pc - the preconditioner context
399: - x - input vector
401: Output Parameter:
402: . y - output vector
404: Level: developer
406: .seealso: `PC`, `PCApplyTranspose()`, `PCApplyBAorAB()`
407: @*/
408: PetscErrorCode PCApply(PC pc, Vec x, Vec y)
409: {
410: PetscInt m, n, mv, nv;
416: if (pc->erroriffailure) VecValidValues_Internal(x, 2, PETSC_TRUE);
417: /* use pmat to check vector sizes since for KSPLSQR the pmat may be of a different size than mat */
418: MatGetLocalSize(pc->pmat, &m, &n);
419: VecGetLocalSize(x, &mv);
420: VecGetLocalSize(y, &nv);
421: /* check pmat * y = x is feasible */
424: VecSetErrorIfLocked(y, 3);
426: PCSetUp(pc);
427: VecLockReadPush(x);
428: PetscLogEventBegin(PC_Apply, pc, x, y, 0);
429: PetscUseTypeMethod(pc, apply, x, y);
430: PetscLogEventEnd(PC_Apply, pc, x, y, 0);
431: if (pc->erroriffailure) VecValidValues_Internal(y, 3, PETSC_FALSE);
432: VecLockReadPop(x);
433: return 0;
434: }
436: /*@
437: PCMatApply - Applies the preconditioner to multiple vectors stored as a `MATDENSE`. Like `PCApply()`, Y and X must be different matrices.
439: Collective
441: Input Parameters:
442: + pc - the preconditioner context
443: - X - block of input vectors
445: Output Parameter:
446: . Y - block of output vectors
448: Level: developer
450: .seealso: `PC`, `PCApply()`, `KSPMatSolve()`
451: @*/
452: PetscErrorCode PCMatApply(PC pc, Mat X, Mat Y)
453: {
454: Mat A;
455: Vec cy, cx;
456: PetscInt m1, M1, m2, M2, n1, N1, n2, N2, m3, M3, n3, N3;
457: PetscBool match;
465: PCGetOperators(pc, NULL, &A);
466: MatGetLocalSize(A, &m3, &n3);
467: MatGetLocalSize(X, &m2, &n2);
468: MatGetLocalSize(Y, &m1, &n1);
469: MatGetSize(A, &M3, &N3);
470: MatGetSize(X, &M2, &N2);
471: MatGetSize(Y, &M1, &N1);
475: PetscObjectBaseTypeCompareAny((PetscObject)Y, &match, MATSEQDENSE, MATMPIDENSE, "");
477: PetscObjectBaseTypeCompareAny((PetscObject)X, &match, MATSEQDENSE, MATMPIDENSE, "");
479: PCSetUp(pc);
480: if (pc->ops->matapply) {
481: PetscLogEventBegin(PC_MatApply, pc, X, Y, 0);
482: PetscUseTypeMethod(pc, matapply, X, Y);
483: PetscLogEventEnd(PC_MatApply, pc, X, Y, 0);
484: } else {
485: PetscInfo(pc, "PC type %s applying column by column\n", ((PetscObject)pc)->type_name);
486: for (n1 = 0; n1 < N1; ++n1) {
487: MatDenseGetColumnVecRead(X, n1, &cx);
488: MatDenseGetColumnVecWrite(Y, n1, &cy);
489: PCApply(pc, cx, cy);
490: MatDenseRestoreColumnVecWrite(Y, n1, &cy);
491: MatDenseRestoreColumnVecRead(X, n1, &cx);
492: }
493: }
494: return 0;
495: }
497: /*@
498: PCApplySymmetricLeft - Applies the left part of a symmetric preconditioner to a vector.
500: Collective
502: Input Parameters:
503: + pc - the preconditioner context
504: - x - input vector
506: Output Parameter:
507: . y - output vector
509: Note:
510: Currently, this routine is implemented only for `PCICC` and `PCJACOBI` preconditioners.
512: Level: developer
514: .seealso: `PC`, `PCApply()`, `PCApplySymmetricRight()`
515: @*/
516: PetscErrorCode PCApplySymmetricLeft(PC pc, Vec x, Vec y)
517: {
522: if (pc->erroriffailure) VecValidValues_Internal(x, 2, PETSC_TRUE);
523: PCSetUp(pc);
524: VecLockReadPush(x);
525: PetscLogEventBegin(PC_ApplySymmetricLeft, pc, x, y, 0);
526: PetscUseTypeMethod(pc, applysymmetricleft, x, y);
527: PetscLogEventEnd(PC_ApplySymmetricLeft, pc, x, y, 0);
528: VecLockReadPop(x);
529: if (pc->erroriffailure) VecValidValues_Internal(y, 3, PETSC_FALSE);
530: return 0;
531: }
533: /*@
534: PCApplySymmetricRight - Applies the right part of a symmetric preconditioner to a vector.
536: Collective
538: Input Parameters:
539: + pc - the preconditioner context
540: - x - input vector
542: Output Parameter:
543: . y - output vector
545: Level: developer
547: Note:
548: Currently, this routine is implemented only for `PCICC` and `PCJACOBI` preconditioners.
550: .seealso: `PC`, `PCApply()`, `PCApplySymmetricLeft()`
551: @*/
552: PetscErrorCode PCApplySymmetricRight(PC pc, Vec x, Vec y)
553: {
558: if (pc->erroriffailure) VecValidValues_Internal(x, 2, PETSC_TRUE);
559: PCSetUp(pc);
560: VecLockReadPush(x);
561: PetscLogEventBegin(PC_ApplySymmetricRight, pc, x, y, 0);
562: PetscUseTypeMethod(pc, applysymmetricright, x, y);
563: PetscLogEventEnd(PC_ApplySymmetricRight, pc, x, y, 0);
564: VecLockReadPop(x);
565: if (pc->erroriffailure) VecValidValues_Internal(y, 3, PETSC_FALSE);
566: return 0;
567: }
569: /*@
570: PCApplyTranspose - Applies the transpose of preconditioner to a vector.
572: Collective
574: Input Parameters:
575: + pc - the preconditioner context
576: - x - input vector
578: Output Parameter:
579: . y - output vector
581: Note:
582: For complex numbers this applies the non-Hermitian transpose.
584: Developer Note:
585: We need to implement a `PCApplyHermitianTranspose()`
587: Level: developer
589: .seealso: `PC`, `PCApply()`, `PCApplyBAorAB()`, `PCApplyBAorABTranspose()`, `PCApplyTransposeExists()`
590: @*/
591: PetscErrorCode PCApplyTranspose(PC pc, Vec x, Vec y)
592: {
597: if (pc->erroriffailure) VecValidValues_Internal(x, 2, PETSC_TRUE);
598: PCSetUp(pc);
599: VecLockReadPush(x);
600: PetscLogEventBegin(PC_Apply, pc, x, y, 0);
601: PetscUseTypeMethod(pc, applytranspose, x, y);
602: PetscLogEventEnd(PC_Apply, pc, x, y, 0);
603: VecLockReadPop(x);
604: if (pc->erroriffailure) VecValidValues_Internal(y, 3, PETSC_FALSE);
605: return 0;
606: }
608: /*@
609: PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation
611: Collective
613: Input Parameter:
614: . pc - the preconditioner context
616: Output Parameter:
617: . flg - `PETSC_TRUE` if a transpose operation is defined
619: Level: developer
621: .seealso: `PC`, `PCApplyTranspose()`
622: @*/
623: PetscErrorCode PCApplyTransposeExists(PC pc, PetscBool *flg)
624: {
627: if (pc->ops->applytranspose) *flg = PETSC_TRUE;
628: else *flg = PETSC_FALSE;
629: return 0;
630: }
632: /*@
633: PCApplyBAorAB - Applies the preconditioner and operator to a vector. y = B*A*x or y = A*B*x.
635: Collective
637: Input Parameters:
638: + pc - the preconditioner context
639: . side - indicates the preconditioner side, one of `PC_LEFT`, `PC_RIGHT`, or `PC_SYMMETRIC`
640: . x - input vector
641: - work - work vector
643: Output Parameter:
644: . y - output vector
646: Level: developer
648: Note:
649: If the `PC` has had `PCSetDiagonalScale()` set then D M A D^{-1} for left preconditioning or D A M D^{-1} is actually applied. Note that the
650: specific `KSPSolve()` method must also be written to handle the post-solve "correction" for the diagonal scaling.
652: .seealso: `PC`, `PCApply()`, `PCApplyTranspose()`, `PCApplyBAorABTranspose()`
653: @*/
654: PetscErrorCode PCApplyBAorAB(PC pc, PCSide side, Vec x, Vec y, Vec work)
655: {
667: if (pc->erroriffailure) VecValidValues_Internal(x, 3, PETSC_TRUE);
669: PCSetUp(pc);
670: if (pc->diagonalscale) {
671: if (pc->ops->applyBA) {
672: Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
673: VecDuplicate(x, &work2);
674: PCDiagonalScaleRight(pc, x, work2);
675: PetscUseTypeMethod(pc, applyBA, side, work2, y, work);
676: PCDiagonalScaleLeft(pc, y, y);
677: VecDestroy(&work2);
678: } else if (side == PC_RIGHT) {
679: PCDiagonalScaleRight(pc, x, y);
680: PCApply(pc, y, work);
681: MatMult(pc->mat, work, y);
682: PCDiagonalScaleLeft(pc, y, y);
683: } else if (side == PC_LEFT) {
684: PCDiagonalScaleRight(pc, x, y);
685: MatMult(pc->mat, y, work);
686: PCApply(pc, work, y);
687: PCDiagonalScaleLeft(pc, y, y);
689: } else {
690: if (pc->ops->applyBA) {
691: PetscUseTypeMethod(pc, applyBA, side, x, y, work);
692: } else if (side == PC_RIGHT) {
693: PCApply(pc, x, work);
694: MatMult(pc->mat, work, y);
695: } else if (side == PC_LEFT) {
696: MatMult(pc->mat, x, work);
697: PCApply(pc, work, y);
698: } else if (side == PC_SYMMETRIC) {
699: /* There's an extra copy here; maybe should provide 2 work vectors instead? */
700: PCApplySymmetricRight(pc, x, work);
701: MatMult(pc->mat, work, y);
702: VecCopy(y, work);
703: PCApplySymmetricLeft(pc, work, y);
704: }
705: }
706: if (pc->erroriffailure) VecValidValues_Internal(y, 4, PETSC_FALSE);
707: return 0;
708: }
710: /*@
711: PCApplyBAorABTranspose - Applies the transpose of the preconditioner
712: and operator to a vector. That is, applies tr(B) * tr(A) with left preconditioning,
713: NOT tr(B*A) = tr(A)*tr(B).
715: Collective
717: Input Parameters:
718: + pc - the preconditioner context
719: . side - indicates the preconditioner side, one of `PC_LEFT`, `PC_RIGHT`, or `PC_SYMMETRIC`
720: . x - input vector
721: - work - work vector
723: Output Parameter:
724: . y - output vector
726: Note:
727: This routine is used internally so that the same Krylov code can be used to solve A x = b and A' x = b, with a preconditioner
728: defined by B'. This is why this has the funny form that it computes tr(B) * tr(A)
730: Level: developer
732: .seealso: `PC`, `PCApply()`, `PCApplyTranspose()`, `PCApplyBAorAB()`
733: @*/
734: PetscErrorCode PCApplyBAorABTranspose(PC pc, PCSide side, Vec x, Vec y, Vec work)
735: {
741: if (pc->erroriffailure) VecValidValues_Internal(x, 3, PETSC_TRUE);
742: if (pc->ops->applyBAtranspose) {
743: PetscUseTypeMethod(pc, applyBAtranspose, side, x, y, work);
744: if (pc->erroriffailure) VecValidValues_Internal(y, 4, PETSC_FALSE);
745: return 0;
746: }
749: PCSetUp(pc);
750: if (side == PC_RIGHT) {
751: PCApplyTranspose(pc, x, work);
752: MatMultTranspose(pc->mat, work, y);
753: } else if (side == PC_LEFT) {
754: MatMultTranspose(pc->mat, x, work);
755: PCApplyTranspose(pc, work, y);
756: }
757: /* add support for PC_SYMMETRIC */
758: if (pc->erroriffailure) VecValidValues_Internal(y, 4, PETSC_FALSE);
759: return 0;
760: }
762: /*@
763: PCApplyRichardsonExists - Determines whether a particular preconditioner has a
764: built-in fast application of Richardson's method.
766: Not Collective
768: Input Parameter:
769: . pc - the preconditioner
771: Output Parameter:
772: . exists - `PETSC_TRUE` or `PETSC_FALSE`
774: Level: developer
776: .seealso: `PC`, `PCRICHARDSON`, `PCApplyRichardson()`
777: @*/
778: PetscErrorCode PCApplyRichardsonExists(PC pc, PetscBool *exists)
779: {
782: if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
783: else *exists = PETSC_FALSE;
784: return 0;
785: }
787: /*@
788: PCApplyRichardson - Applies several steps of Richardson iteration with
789: the particular preconditioner. This routine is usually used by the
790: Krylov solvers and not the application code directly.
792: Collective
794: Input Parameters:
795: + pc - the preconditioner context
796: . b - the right hand side
797: . w - one work vector
798: . rtol - relative decrease in residual norm convergence criteria
799: . abstol - absolute residual norm convergence criteria
800: . dtol - divergence residual norm increase criteria
801: . its - the number of iterations to apply.
802: - guesszero - if the input x contains nonzero initial guess
804: Output Parameters:
805: + outits - number of iterations actually used (for SOR this always equals its)
806: . reason - the reason the apply terminated
807: - y - the solution (also contains initial guess if guesszero is `PETSC_FALSE`
809: Notes:
810: Most preconditioners do not support this function. Use the command
811: `PCApplyRichardsonExists()` to determine if one does.
813: Except for the `PCMG` this routine ignores the convergence tolerances
814: and always runs for the number of iterations
816: Level: developer
818: .seealso: `PC`, `PCApplyRichardsonExists()`
819: @*/
820: PetscErrorCode PCApplyRichardson(PC pc, Vec b, Vec y, Vec w, PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt its, PetscBool guesszero, PetscInt *outits, PCRichardsonConvergedReason *reason)
821: {
827: PCSetUp(pc);
828: PetscUseTypeMethod(pc, applyrichardson, b, y, w, rtol, abstol, dtol, its, guesszero, outits, reason);
829: return 0;
830: }
832: /*@
833: PCSetFailedReason - Sets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail
835: Logically Collective
837: Input Parameters:
838: + pc - the preconditioner context
839: - reason - the reason it failedx
841: Level: advanced
843: .seealso: `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCFailedReason`
844: @*/
845: PetscErrorCode PCSetFailedReason(PC pc, PCFailedReason reason)
846: {
847: pc->failedreason = reason;
848: return 0;
849: }
851: /*@
852: PCGetFailedReason - Gets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail
854: Logically Collective
856: Input Parameter:
857: . pc - the preconditioner context
859: Output Parameter:
860: . reason - the reason it failed
862: Level: advanced
864: Note:
865: This is the maximum over reason over all ranks in the PC communicator. It is only valid after
866: a call `KSPCheckDot()` or `KSPCheckNorm()` inside a `KSPSolve()`. It is not valid immediately after a `PCSetUp()`
867: or `PCApply()`, then use `PCGetFailedReasonRank()`
869: .seealso: PC`, ``PCCreate()`, `PCApply()`, `PCDestroy()`, `PCGetFailedReasonRank()`, `PCSetFailedReason()`
870: @*/
871: PetscErrorCode PCGetFailedReason(PC pc, PCFailedReason *reason)
872: {
873: if (pc->setupcalled < 0) *reason = (PCFailedReason)pc->setupcalled;
874: else *reason = pc->failedreason;
875: return 0;
876: }
878: /*@
879: PCGetFailedReasonRank - Gets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail on this MPI rank
881: Not Collective
883: Input Parameter:
884: . pc - the preconditioner context
886: Output Parameter:
887: . reason - the reason it failed
889: Note:
890: Different ranks may have different reasons or no reason, see `PCGetFailedReason()`
892: Level: advanced
894: .seealso: `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCGetFailedReason()`, `PCSetFailedReason()`
895: @*/
896: PetscErrorCode PCGetFailedReasonRank(PC pc, PCFailedReason *reason)
897: {
898: if (pc->setupcalled < 0) *reason = (PCFailedReason)pc->setupcalled;
899: else *reason = pc->failedreason;
900: return 0;
901: }
903: /* Next line needed to deactivate KSP_Solve logging */
904: #include <petsc/private/kspimpl.h>
906: /*
907: a setupcall of 0 indicates never setup,
908: 1 indicates has been previously setup
909: -1 indicates a PCSetUp() was attempted and failed
910: */
911: /*@
912: PCSetUp - Prepares for the use of a preconditioner.
914: Collective
916: Input Parameter:
917: . pc - the preconditioner context
919: Level: developer
921: .seealso: `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`
922: @*/
923: PetscErrorCode PCSetUp(PC pc)
924: {
925: const char *def;
926: PetscObjectState matstate, matnonzerostate;
931: if (pc->setupcalled && pc->reusepreconditioner) {
932: PetscInfo(pc, "Leaving PC with identical preconditioner since reuse preconditioner is set\n");
933: return 0;
934: }
936: PetscObjectStateGet((PetscObject)pc->pmat, &matstate);
937: MatGetNonzeroState(pc->pmat, &matnonzerostate);
938: if (!pc->setupcalled) {
939: PetscInfo(pc, "Setting up PC for first time\n");
940: pc->flag = DIFFERENT_NONZERO_PATTERN;
941: } else if (matstate == pc->matstate) {
942: PetscInfo(pc, "Leaving PC with identical preconditioner since operator is unchanged\n");
943: return 0;
944: } else {
945: if (matnonzerostate > pc->matnonzerostate) {
946: PetscInfo(pc, "Setting up PC with different nonzero pattern\n");
947: pc->flag = DIFFERENT_NONZERO_PATTERN;
948: } else {
949: PetscInfo(pc, "Setting up PC with same nonzero pattern\n");
950: pc->flag = SAME_NONZERO_PATTERN;
951: }
952: }
953: pc->matstate = matstate;
954: pc->matnonzerostate = matnonzerostate;
956: if (!((PetscObject)pc)->type_name) {
957: PCGetDefaultType_Private(pc, &def);
958: PCSetType(pc, def);
959: }
961: MatSetErrorIfFailure(pc->pmat, pc->erroriffailure);
962: MatSetErrorIfFailure(pc->mat, pc->erroriffailure);
963: PetscLogEventBegin(PC_SetUp, pc, 0, 0, 0);
964: if (pc->ops->setup) {
965: /* do not log solves and applications of preconditioners while constructing preconditioners; perhaps they should be logged separately from the regular solves */
966: KSPInitializePackage();
967: PetscLogEventDeactivatePush(KSP_Solve);
968: PetscLogEventDeactivatePush(PC_Apply);
969: PetscUseTypeMethod(pc, setup);
970: PetscLogEventDeactivatePop(KSP_Solve);
971: PetscLogEventDeactivatePop(PC_Apply);
972: }
973: PetscLogEventEnd(PC_SetUp, pc, 0, 0, 0);
974: if (!pc->setupcalled) pc->setupcalled = 1;
975: return 0;
976: }
978: /*@
979: PCSetUpOnBlocks - Sets up the preconditioner for each block in
980: the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
981: methods.
983: Collective
985: Input Parameter:
986: . pc - the preconditioner context
988: Level: developer
990: Note:
991: For nested preconditioners such as `PCBJACOBI` `PCSetUp()` is not called on each sub-`KSP` when `PCSetUp()` is
992: called on the outer `PC`, this routine ensures it is called.
994: .seealso: `PC`, `PCSetUp()`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCSetUp()`
995: @*/
996: PetscErrorCode PCSetUpOnBlocks(PC pc)
997: {
999: if (!pc->ops->setuponblocks) return 0;
1000: PetscLogEventBegin(PC_SetUpOnBlocks, pc, 0, 0, 0);
1001: PetscUseTypeMethod(pc, setuponblocks);
1002: PetscLogEventEnd(PC_SetUpOnBlocks, pc, 0, 0, 0);
1003: return 0;
1004: }
1006: /*@C
1007: PCSetModifySubMatrices - Sets a user-defined routine for modifying the
1008: submatrices that arise within certain subdomain-based preconditioners.
1009: The basic submatrices are extracted from the preconditioner matrix as
1010: usual; the user can then alter these (for example, to set different boundary
1011: conditions for each submatrix) before they are used for the local solves.
1013: Logically Collective
1015: Input Parameters:
1016: + pc - the preconditioner context
1017: . func - routine for modifying the submatrices
1018: - ctx - optional user-defined context (may be null)
1020: Calling sequence of func:
1021: $ func (PC pc,PetscInt nsub,IS *row,IS *col,Mat *submat,void *ctx);
1023: + row - an array of index sets that contain the global row numbers
1024: that comprise each local submatrix
1025: . col - an array of index sets that contain the global column numbers
1026: that comprise each local submatrix
1027: . submat - array of local submatrices
1028: - ctx - optional user-defined context for private data for the
1029: user-defined func routine (may be null)
1031: Notes:
1032: `PCSetModifySubMatrices()` MUST be called before `KSPSetUp()` and
1033: `KSPSolve()`.
1035: A routine set by `PCSetModifySubMatrices()` is currently called within
1036: the block Jacobi (`PCBJACOBI`) and additive Schwarz (`PCASM`)
1037: preconditioners. All other preconditioners ignore this routine.
1039: Level: advanced
1041: .seealso: `PC`, `PCBJACOBI`, `PCASM`, `PCModifySubMatrices()`
1042: @*/
1043: PetscErrorCode PCSetModifySubMatrices(PC pc, PetscErrorCode (*func)(PC, PetscInt, const IS[], const IS[], Mat[], void *), void *ctx)
1044: {
1046: pc->modifysubmatrices = func;
1047: pc->modifysubmatricesP = ctx;
1048: return 0;
1049: }
1051: /*@C
1052: PCModifySubMatrices - Calls an optional user-defined routine within
1053: certain preconditioners if one has been set with `PCSetModifySubMatrices()`.
1055: Collective
1057: Input Parameters:
1058: + pc - the preconditioner context
1059: . nsub - the number of local submatrices
1060: . row - an array of index sets that contain the global row numbers
1061: that comprise each local submatrix
1062: . col - an array of index sets that contain the global column numbers
1063: that comprise each local submatrix
1064: . submat - array of local submatrices
1065: - ctx - optional user-defined context for private data for the
1066: user-defined routine (may be null)
1068: Output Parameter:
1069: . submat - array of local submatrices (the entries of which may
1070: have been modified)
1072: Notes:
1073: The user should NOT generally call this routine, as it will
1074: automatically be called within certain preconditioners (currently
1075: block Jacobi, additive Schwarz) if set.
1077: The basic submatrices are extracted from the preconditioner matrix
1078: as usual; the user can then alter these (for example, to set different
1079: boundary conditions for each submatrix) before they are used for the
1080: local solves.
1082: Level: developer
1084: .seealso: `PC`, `PCSetModifySubMatrices()`
1085: @*/
1086: PetscErrorCode PCModifySubMatrices(PC pc, PetscInt nsub, const IS row[], const IS col[], Mat submat[], void *ctx)
1087: {
1089: if (!pc->modifysubmatrices) return 0;
1090: PetscLogEventBegin(PC_ModifySubMatrices, pc, 0, 0, 0);
1091: (*pc->modifysubmatrices)(pc, nsub, row, col, submat, ctx);
1092: PetscLogEventEnd(PC_ModifySubMatrices, pc, 0, 0, 0);
1093: return 0;
1094: }
1096: /*@
1097: PCSetOperators - Sets the matrix associated with the linear system and
1098: a (possibly) different one associated with the preconditioner.
1100: Logically Collective
1102: Input Parameters:
1103: + pc - the preconditioner context
1104: . Amat - the matrix that defines the linear system
1105: - Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
1107: Notes:
1108: Passing a NULL for Amat or Pmat removes the matrix that is currently used.
1110: If you wish to replace either Amat or Pmat but leave the other one untouched then
1111: first call `KSPGetOperators()` to get the one you wish to keep, call `PetscObjectReference()`
1112: on it and then pass it back in in your call to `KSPSetOperators()`.
1114: More Notes about Repeated Solution of Linear Systems:
1115: PETSc does NOT reset the matrix entries of either Amat or Pmat
1116: to zero after a linear solve; the user is completely responsible for
1117: matrix assembly. See the routine `MatZeroEntries()` if desiring to
1118: zero all elements of a matrix.
1120: Level: intermediate
1122: .seealso: `PC`, `PCGetOperators()`, `MatZeroEntries()`
1123: @*/
1124: PetscErrorCode PCSetOperators(PC pc, Mat Amat, Mat Pmat)
1125: {
1126: PetscInt m1, n1, m2, n2;
1133: if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
1134: MatGetLocalSize(Amat, &m1, &n1);
1135: MatGetLocalSize(pc->mat, &m2, &n2);
1137: MatGetLocalSize(Pmat, &m1, &n1);
1138: MatGetLocalSize(pc->pmat, &m2, &n2);
1140: }
1142: if (Pmat != pc->pmat) {
1143: /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
1144: pc->matnonzerostate = -1;
1145: pc->matstate = -1;
1146: }
1148: /* reference first in case the matrices are the same */
1149: if (Amat) PetscObjectReference((PetscObject)Amat);
1150: MatDestroy(&pc->mat);
1151: if (Pmat) PetscObjectReference((PetscObject)Pmat);
1152: MatDestroy(&pc->pmat);
1153: pc->mat = Amat;
1154: pc->pmat = Pmat;
1155: return 0;
1156: }
1158: /*@
1159: PCSetReusePreconditioner - reuse the current preconditioner even if the operator in the preconditioner has changed.
1161: Logically Collective
1163: Input Parameters:
1164: + pc - the preconditioner context
1165: - flag - `PETSC_TRUE` do not compute a new preconditioner, `PETSC_FALSE` do compute a new preconditioner
1167: Level: intermediate
1169: Note:
1170: Normally if a matrix inside a `PC` changes the `PC` automatically updates itself using information from the changed matrix. This option
1171: prevents this.
1173: .seealso: `PC`, `PCGetOperators()`, `MatZeroEntries()`, `PCGetReusePreconditioner()`, `KSPSetReusePreconditioner()`
1174: @*/
1175: PetscErrorCode PCSetReusePreconditioner(PC pc, PetscBool flag)
1176: {
1179: pc->reusepreconditioner = flag;
1180: return 0;
1181: }
1183: /*@
1184: PCGetReusePreconditioner - Determines if the `PC` reuses the current preconditioner even if the operator in the preconditioner has changed.
1186: Not Collective
1188: Input Parameter:
1189: . pc - the preconditioner context
1191: Output Parameter:
1192: . flag - `PETSC_TRUE` do not compute a new preconditioner, `PETSC_FALSE` do compute a new preconditioner
1194: Level: intermediate
1196: .seealso: `PC`, `PCGetOperators()`, `MatZeroEntries()`, `PCSetReusePreconditioner()`
1197: @*/
1198: PetscErrorCode PCGetReusePreconditioner(PC pc, PetscBool *flag)
1199: {
1202: *flag = pc->reusepreconditioner;
1203: return 0;
1204: }
1206: /*@
1207: PCGetOperators - Gets the matrix associated with the linear system and
1208: possibly a different one associated with the preconditioner.
1210: Not collective, though parallel mats are returned if pc is parallel
1212: Input Parameter:
1213: . pc - the preconditioner context
1215: Output Parameters:
1216: + Amat - the matrix defining the linear system
1217: - Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.
1219: Level: intermediate
1221: Note:
1222: Does not increase the reference count of the matrices, so you should not destroy them
1224: Alternative usage: If the operators have NOT been set with `KSPSetOperators()`/`PCSetOperators()` then the operators
1225: are created in `PC` and returned to the user. In this case, if both operators
1226: mat and pmat are requested, two DIFFERENT operators will be returned. If
1227: only one is requested both operators in the PC will be the same (i.e. as
1228: if one had called `KSPSetOperators()`/`PCSetOperators()` with the same argument for both Mats).
1229: The user must set the sizes of the returned matrices and their type etc just
1230: as if the user created them with `MatCreate()`. For example,
1232: .vb
1233: KSP/PCGetOperators(ksp/pc,&Amat,NULL); is equivalent to
1234: set size, type, etc of Amat
1236: MatCreate(comm,&mat);
1237: KSP/PCSetOperators(ksp/pc,Amat,Amat);
1238: PetscObjectDereference((PetscObject)mat);
1239: set size, type, etc of Amat
1240: .ve
1242: and
1244: .vb
1245: KSP/PCGetOperators(ksp/pc,&Amat,&Pmat); is equivalent to
1246: set size, type, etc of Amat and Pmat
1248: MatCreate(comm,&Amat);
1249: MatCreate(comm,&Pmat);
1250: KSP/PCSetOperators(ksp/pc,Amat,Pmat);
1251: PetscObjectDereference((PetscObject)Amat);
1252: PetscObjectDereference((PetscObject)Pmat);
1253: set size, type, etc of Amat and Pmat
1254: .ve
1256: The rationale for this support is so that when creating a `TS`, `SNES`, or `KSP` the hierarchy
1257: of underlying objects (i.e. `SNES`, `KSP`, `PC`, `Mat`) and their lifespans can be completely
1258: managed by the top most level object (i.e. the `TS`, `SNES`, or `KSP`). Another way to look
1259: at this is when you create a `SNES` you do not NEED to create a `KSP` and attach it to
1260: the `SNES` object (the `SNES` object manages it for you). Similarly when you create a KSP
1261: you do not need to attach a `PC` to it (the `KSP` object manages the `PC` object for you).
1262: Thus, why should YOU have to create the `Mat` and attach it to the `SNES`/`KSP`/`PC`, when
1263: it can be created for you?
1265: .seealso: `PC`, `PCSetOperators()`, `KSPGetOperators()`, `KSPSetOperators()`, `PCGetOperatorsSet()`
1266: @*/
1267: PetscErrorCode PCGetOperators(PC pc, Mat *Amat, Mat *Pmat)
1268: {
1270: if (Amat) {
1271: if (!pc->mat) {
1272: if (pc->pmat && !Pmat) { /* Apmat has been set, but user did not request it, so use for Amat */
1273: pc->mat = pc->pmat;
1274: PetscObjectReference((PetscObject)pc->mat);
1275: } else { /* both Amat and Pmat are empty */
1276: MatCreate(PetscObjectComm((PetscObject)pc), &pc->mat);
1277: if (!Pmat) { /* user did NOT request Pmat, so make same as Amat */
1278: pc->pmat = pc->mat;
1279: PetscObjectReference((PetscObject)pc->pmat);
1280: }
1281: }
1282: }
1283: *Amat = pc->mat;
1284: }
1285: if (Pmat) {
1286: if (!pc->pmat) {
1287: if (pc->mat && !Amat) { /* Amat has been set but was not requested, so use for pmat */
1288: pc->pmat = pc->mat;
1289: PetscObjectReference((PetscObject)pc->pmat);
1290: } else {
1291: MatCreate(PetscObjectComm((PetscObject)pc), &pc->pmat);
1292: if (!Amat) { /* user did NOT request Amat, so make same as Pmat */
1293: pc->mat = pc->pmat;
1294: PetscObjectReference((PetscObject)pc->mat);
1295: }
1296: }
1297: }
1298: *Pmat = pc->pmat;
1299: }
1300: return 0;
1301: }
1303: /*@C
1304: PCGetOperatorsSet - Determines if the matrix associated with the linear system and
1305: possibly a different one associated with the preconditioner have been set in the `PC`.
1307: Not collective, though the results on all processes should be the same
1309: Input Parameter:
1310: . pc - the preconditioner context
1312: Output Parameters:
1313: + mat - the matrix associated with the linear system was set
1314: - pmat - matrix associated with the preconditioner was set, usually the same
1316: Level: intermediate
1318: .seealso: `PC`, `PCSetOperators()`, `KSPGetOperators()`, `KSPSetOperators()`, `PCGetOperators()`
1319: @*/
1320: PetscErrorCode PCGetOperatorsSet(PC pc, PetscBool *mat, PetscBool *pmat)
1321: {
1323: if (mat) *mat = (pc->mat) ? PETSC_TRUE : PETSC_FALSE;
1324: if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1325: return 0;
1326: }
1328: /*@
1329: PCFactorGetMatrix - Gets the factored matrix from the
1330: preconditioner context. This routine is valid only for the `PCLU`,
1331: `PCILU`, `PCCHOLESKY`, and `PCICC` methods.
1333: Not Collective though mat is parallel if pc is parallel
1335: Input Parameter:
1336: . pc - the preconditioner context
1338: Output parameters:
1339: . mat - the factored matrix
1341: Level: advanced
1343: Note:
1344: Does not increase the reference count for the matrix so DO NOT destroy it
1346: .seealso: `PC`, `PCLU`, `PCILU`, `PCCHOLESKY`, `PCICC`
1347: @*/
1348: PetscErrorCode PCFactorGetMatrix(PC pc, Mat *mat)
1349: {
1352: PetscUseTypeMethod(pc, getfactoredmatrix, mat);
1353: return 0;
1354: }
1356: /*@C
1357: PCSetOptionsPrefix - Sets the prefix used for searching for all
1358: `PC` options in the database.
1360: Logically Collective
1362: Input Parameters:
1363: + pc - the preconditioner context
1364: - prefix - the prefix string to prepend to all `PC` option requests
1366: Note:
1367: A hyphen (-) must NOT be given at the beginning of the prefix name.
1368: The first character of all runtime options is AUTOMATICALLY the
1369: hyphen.
1371: Level: advanced
1373: .seealso: `PC`, `PCSetFromOptions`, `PCAppendOptionsPrefix()`, `PCGetOptionsPrefix()`
1374: @*/
1375: PetscErrorCode PCSetOptionsPrefix(PC pc, const char prefix[])
1376: {
1378: PetscObjectSetOptionsPrefix((PetscObject)pc, prefix);
1379: return 0;
1380: }
1382: /*@C
1383: PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1384: `PC` options in the database.
1386: Logically Collective
1388: Input Parameters:
1389: + pc - the preconditioner context
1390: - prefix - the prefix string to prepend to all `PC` option requests
1392: Note:
1393: A hyphen (-) must NOT be given at the beginning of the prefix name.
1394: The first character of all runtime options is AUTOMATICALLY the
1395: hyphen.
1397: Level: advanced
1399: .seealso: `PC`, `PCSetFromOptions`, `PCSetOptionsPrefix()`, `PCGetOptionsPrefix()`
1400: @*/
1401: PetscErrorCode PCAppendOptionsPrefix(PC pc, const char prefix[])
1402: {
1404: PetscObjectAppendOptionsPrefix((PetscObject)pc, prefix);
1405: return 0;
1406: }
1408: /*@C
1409: PCGetOptionsPrefix - Gets the prefix used for searching for all
1410: PC options in the database.
1412: Not Collective
1414: Input Parameter:
1415: . pc - the preconditioner context
1417: Output Parameter:
1418: . prefix - pointer to the prefix string used, is returned
1420: Fortran Note:
1421: The user should pass in a string 'prefix' of
1422: sufficient length to hold the prefix.
1424: Level: advanced
1426: .seealso: `PC`, `PCSetFromOptions`, `PCSetOptionsPrefix()`, `PCAppendOptionsPrefix()`
1427: @*/
1428: PetscErrorCode PCGetOptionsPrefix(PC pc, const char *prefix[])
1429: {
1432: PetscObjectGetOptionsPrefix((PetscObject)pc, prefix);
1433: return 0;
1434: }
1436: /*
1437: Indicates the right hand side will be changed by KSPSolve(), this occurs for a few
1438: preconditioners including BDDC and Eisentat that transform the equations before applying
1439: the Krylov methods
1440: */
1441: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC pc, PetscBool *change)
1442: {
1445: *change = PETSC_FALSE;
1446: PetscTryMethod(pc, "PCPreSolveChangeRHS_C", (PC, PetscBool *), (pc, change));
1447: return 0;
1448: }
1450: /*@
1451: PCPreSolve - Optional pre-solve phase, intended for any
1452: preconditioner-specific actions that must be performed before
1453: the iterative solve itself.
1455: Collective
1457: Input Parameters:
1458: + pc - the preconditioner context
1459: - ksp - the Krylov subspace context
1461: Level: developer
1463: Sample of Usage:
1464: .vb
1465: PCPreSolve(pc,ksp);
1466: KSPSolve(ksp,b,x);
1467: PCPostSolve(pc,ksp);
1468: .ve
1470: Notes:
1471: The pre-solve phase is distinct from the `PCSetUp()` phase.
1473: `KSPSolve()` calls this directly, so is rarely called by the user.
1475: .seealso: `PC`, `PCPostSolve()`
1476: @*/
1477: PetscErrorCode PCPreSolve(PC pc, KSP ksp)
1478: {
1479: Vec x, rhs;
1483: pc->presolvedone++;
1485: KSPGetSolution(ksp, &x);
1486: KSPGetRhs(ksp, &rhs);
1488: if (pc->ops->presolve) PetscUseTypeMethod(pc, presolve, ksp, rhs, x);
1489: else if (pc->presolve) (pc->presolve)(pc, ksp);
1490: return 0;
1491: }
1493: /*@C
1494: PCSetPreSolve - Sets function used by `PCPreSolve()` which is intended for any
1495: preconditioner-specific actions that must be performed before
1496: the iterative solve itself.
1498: Logically Collective
1500: Input Parameters:
1501: + pc - the preconditioner object
1502: - presolve - the function to call before the solve
1504: Calling sequence of presolve:
1505: $ func(PC pc,KSP ksp)
1507: + pc - the PC context
1508: - ksp - the KSP context
1510: Level: developer
1512: .seealso: `PC`, `PCSetUp()`, `PCPreSolve()`
1513: @*/
1514: PetscErrorCode PCSetPreSolve(PC pc, PetscErrorCode (*presolve)(PC, KSP))
1515: {
1517: pc->presolve = presolve;
1518: return 0;
1519: }
1521: /*@
1522: PCPostSolve - Optional post-solve phase, intended for any
1523: preconditioner-specific actions that must be performed after
1524: the iterative solve itself.
1526: Collective
1528: Input Parameters:
1529: + pc - the preconditioner context
1530: - ksp - the Krylov subspace context
1532: Sample of Usage:
1533: .vb
1534: PCPreSolve(pc,ksp);
1535: KSPSolve(ksp,b,x);
1536: PCPostSolve(pc,ksp);
1537: .ve
1539: Note:
1540: `KSPSolve()` calls this routine directly, so it is rarely called by the user.
1542: Level: developer
1544: .seealso: `PC`, `PCSetPostSolve()`, `PCSetPresolve()`, `PCPreSolve()`, `KSPSolve()`
1545: @*/
1546: PetscErrorCode PCPostSolve(PC pc, KSP ksp)
1547: {
1548: Vec x, rhs;
1552: pc->presolvedone--;
1553: KSPGetSolution(ksp, &x);
1554: KSPGetRhs(ksp, &rhs);
1555: PetscTryTypeMethod(pc, postsolve, ksp, rhs, x);
1556: return 0;
1557: }
1559: /*@C
1560: PCLoad - Loads a `PC` that has been stored in binary with `PCView()`.
1562: Collective
1564: Input Parameters:
1565: + newdm - the newly loaded `PC`, this needs to have been created with `PCCreate()` or
1566: some related function before a call to `PCLoad()`.
1567: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()`
1569: Level: intermediate
1571: Note:
1572: The type is determined by the data in the file, any type set into the PC before this call is ignored.
1574: .seealso: `PC`, `PetscViewerBinaryOpen()`, `PCView()`, `MatLoad()`, `VecLoad()`
1575: @*/
1576: PetscErrorCode PCLoad(PC newdm, PetscViewer viewer)
1577: {
1578: PetscBool isbinary;
1579: PetscInt classid;
1580: char type[256];
1584: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary);
1587: PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT);
1589: PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR);
1590: PCSetType(newdm, type);
1591: PetscTryTypeMethod(newdm, load, viewer);
1592: return 0;
1593: }
1595: #include <petscdraw.h>
1596: #if defined(PETSC_HAVE_SAWS)
1597: #include <petscviewersaws.h>
1598: #endif
1600: /*@C
1601: PCViewFromOptions - View from the `PC` based on options in the database
1603: Collective
1605: Input Parameters:
1606: + A - the PC context
1607: . obj - Optional object that provides the options prefix
1608: - name - command line option
1610: Level: intermediate
1612: .seealso: `PC`, `PCView`, `PetscObjectViewFromOptions()`, `PCCreate()`
1613: @*/
1614: PetscErrorCode PCViewFromOptions(PC A, PetscObject obj, const char name[])
1615: {
1617: PetscObjectViewFromOptions((PetscObject)A, obj, name);
1618: return 0;
1619: }
1621: /*@C
1622: PCView - Prints information about the `PC`
1624: Collective
1626: Input Parameters:
1627: + PC - the `PC` context
1628: - viewer - optional visualization context
1630: Notes:
1631: The available visualization contexts include
1632: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
1633: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
1634: output where only the first processor opens
1635: the file. All other processors send their
1636: data to the first processor to print.
1638: The user can open an alternative visualization contexts with
1639: `PetscViewerASCIIOpen()` (output to a specified file).
1641: Level: developer
1643: .seealso: `PC`, `PetscViewer`, `KSPView()`, `PetscViewerASCIIOpen()`
1644: @*/
1645: PetscErrorCode PCView(PC pc, PetscViewer viewer)
1646: {
1647: PCType cstr;
1648: PetscBool iascii, isstring, isbinary, isdraw;
1649: #if defined(PETSC_HAVE_SAWS)
1650: PetscBool issaws;
1651: #endif
1654: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc), &viewer);
1658: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1659: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring);
1660: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary);
1661: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
1662: #if defined(PETSC_HAVE_SAWS)
1663: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws);
1664: #endif
1666: if (iascii) {
1667: PetscObjectPrintClassNamePrefixType((PetscObject)pc, viewer);
1668: if (!pc->setupcalled) PetscViewerASCIIPrintf(viewer, " PC has not been set up so information may be incomplete\n");
1669: PetscViewerASCIIPushTab(viewer);
1670: PetscTryTypeMethod(pc, view, viewer);
1671: PetscViewerASCIIPopTab(viewer);
1672: if (pc->mat) {
1673: PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);
1674: if (pc->pmat == pc->mat) {
1675: PetscViewerASCIIPrintf(viewer, " linear system matrix = precond matrix:\n");
1676: PetscViewerASCIIPushTab(viewer);
1677: MatView(pc->mat, viewer);
1678: PetscViewerASCIIPopTab(viewer);
1679: } else {
1680: if (pc->pmat) {
1681: PetscViewerASCIIPrintf(viewer, " linear system matrix followed by preconditioner matrix:\n");
1682: } else {
1683: PetscViewerASCIIPrintf(viewer, " linear system matrix:\n");
1684: }
1685: PetscViewerASCIIPushTab(viewer);
1686: MatView(pc->mat, viewer);
1687: if (pc->pmat) MatView(pc->pmat, viewer);
1688: PetscViewerASCIIPopTab(viewer);
1689: }
1690: PetscViewerPopFormat(viewer);
1691: }
1692: } else if (isstring) {
1693: PCGetType(pc, &cstr);
1694: PetscViewerStringSPrintf(viewer, " PCType: %-7.7s", cstr);
1695: PetscTryTypeMethod(pc, view, viewer);
1696: if (pc->mat) MatView(pc->mat, viewer);
1697: if (pc->pmat && pc->pmat != pc->mat) MatView(pc->pmat, viewer);
1698: } else if (isbinary) {
1699: PetscInt classid = PC_FILE_CLASSID;
1700: MPI_Comm comm;
1701: PetscMPIInt rank;
1702: char type[256];
1704: PetscObjectGetComm((PetscObject)pc, &comm);
1705: MPI_Comm_rank(comm, &rank);
1706: if (rank == 0) {
1707: PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT);
1708: PetscStrncpy(type, ((PetscObject)pc)->type_name, 256);
1709: PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR);
1710: }
1711: PetscTryTypeMethod(pc, view, viewer);
1712: } else if (isdraw) {
1713: PetscDraw draw;
1714: char str[25];
1715: PetscReal x, y, bottom, h;
1716: PetscInt n;
1718: PetscViewerDrawGetDraw(viewer, 0, &draw);
1719: PetscDrawGetCurrentPoint(draw, &x, &y);
1720: if (pc->mat) {
1721: MatGetSize(pc->mat, &n, NULL);
1722: PetscSNPrintf(str, 25, "PC: %s (%" PetscInt_FMT ")", ((PetscObject)pc)->type_name, n);
1723: } else {
1724: PetscSNPrintf(str, 25, "PC: %s", ((PetscObject)pc)->type_name);
1725: }
1726: PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h);
1727: bottom = y - h;
1728: PetscDrawPushCurrentPoint(draw, x, bottom);
1729: PetscTryTypeMethod(pc, view, viewer);
1730: PetscDrawPopCurrentPoint(draw);
1731: #if defined(PETSC_HAVE_SAWS)
1732: } else if (issaws) {
1733: PetscMPIInt rank;
1735: PetscObjectName((PetscObject)pc);
1736: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
1737: if (!((PetscObject)pc)->amsmem && rank == 0) PetscObjectViewSAWs((PetscObject)pc, viewer);
1738: if (pc->mat) MatView(pc->mat, viewer);
1739: if (pc->pmat && pc->pmat != pc->mat) MatView(pc->pmat, viewer);
1740: #endif
1741: }
1742: return 0;
1743: }
1745: /*@C
1746: PCRegister - Adds a method to the preconditioner package.
1748: Not collective
1750: Input Parameters:
1751: + name_solver - name of a new user-defined solver
1752: - routine_create - routine to create method context
1754: Note:
1755: `PCRegister()` may be called multiple times to add several user-defined preconditioners.
1757: Sample usage:
1758: .vb
1759: PCRegister("my_solver", MySolverCreate);
1760: .ve
1762: Then, your solver can be chosen with the procedural interface via
1763: $ PCSetType(pc,"my_solver")
1764: or at runtime via the option
1765: $ -pc_type my_solver
1767: Level: advanced
1769: .seealso: `PC`, `PCType`, `PCRegisterAll()`
1770: @*/
1771: PetscErrorCode PCRegister(const char sname[], PetscErrorCode (*function)(PC))
1772: {
1773: PCInitializePackage();
1774: PetscFunctionListAdd(&PCList, sname, function);
1775: return 0;
1776: }
1778: static PetscErrorCode MatMult_PC(Mat A, Vec X, Vec Y)
1779: {
1780: PC pc;
1782: MatShellGetContext(A, &pc);
1783: PCApply(pc, X, Y);
1784: return 0;
1785: }
1787: /*@
1788: PCComputeOperator - Computes the explicit preconditioned operator.
1790: Collective
1792: Input Parameters:
1793: + pc - the preconditioner object
1794: - mattype - the matrix type to be used for the operator
1796: Output Parameter:
1797: . mat - the explicit preconditioned operator
1799: Note:
1800: This computation is done by applying the operators to columns of the identity matrix.
1801: This routine is costly in general, and is recommended for use only with relatively small systems.
1802: Currently, this routine uses a dense matrix format when mattype == NULL
1804: Level: advanced
1806: .seealso: `PC`, `KSPComputeOperator()`, `MatType`
1808: @*/
1809: PetscErrorCode PCComputeOperator(PC pc, MatType mattype, Mat *mat)
1810: {
1811: PetscInt N, M, m, n;
1812: Mat A, Apc;
1816: PCGetOperators(pc, &A, NULL);
1817: MatGetLocalSize(A, &m, &n);
1818: MatGetSize(A, &M, &N);
1819: MatCreateShell(PetscObjectComm((PetscObject)pc), m, n, M, N, pc, &Apc);
1820: MatShellSetOperation(Apc, MATOP_MULT, (void (*)(void))MatMult_PC);
1821: MatComputeOperator(Apc, mattype, mat);
1822: MatDestroy(&Apc);
1823: return 0;
1824: }
1826: /*@
1827: PCSetCoordinates - sets the coordinates of all the nodes on the local process
1829: Collective
1831: Input Parameters:
1832: + pc - the solver context
1833: . dim - the dimension of the coordinates 1, 2, or 3
1834: . nloc - the blocked size of the coordinates array
1835: - coords - the coordinates array
1837: Level: intermediate
1839: Note:
1840: coords is an array of the dim coordinates for the nodes on
1841: the local processor, of size dim*nloc.
1842: If there are 108 equation on a processor
1843: for a displacement finite element discretization of elasticity (so
1844: that there are nloc = 36 = 108/3 nodes) then the array must have 108
1845: double precision values (ie, 3 * 36). These x y z coordinates
1846: should be ordered for nodes 0 to N-1 like so: [ 0.x, 0.y, 0.z, 1.x,
1847: ... , N-1.z ].
1849: .seealso: `PC`, `MatSetNearNullSpace()`
1850: @*/
1851: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[])
1852: {
1855: PetscTryMethod(pc, "PCSetCoordinates_C", (PC, PetscInt, PetscInt, PetscReal *), (pc, dim, nloc, coords));
1856: return 0;
1857: }
1859: /*@
1860: PCGetInterpolations - Gets interpolation matrices for all levels (except level 0)
1862: Logically Collective
1864: Input Parameter:
1865: . pc - the precondition context
1867: Output Parameters:
1868: + num_levels - the number of levels
1869: - interpolations - the interpolation matrices (size of num_levels-1)
1871: Level: advanced
1873: Developer Note:
1874: Why is this here instead of in `PCMG` etc?
1876: .seealso: `PC`, `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetInterpolation()`, `PCGetCoarseOperators()`
1877: @*/
1878: PetscErrorCode PCGetInterpolations(PC pc, PetscInt *num_levels, Mat *interpolations[])
1879: {
1883: PetscUseMethod(pc, "PCGetInterpolations_C", (PC, PetscInt *, Mat *[]), (pc, num_levels, interpolations));
1884: return 0;
1885: }
1887: /*@
1888: PCGetCoarseOperators - Gets coarse operator matrices for all levels (except the finest level)
1890: Logically Collective
1892: Input Parameter:
1893: . pc - the precondition context
1895: Output Parameters:
1896: + num_levels - the number of levels
1897: - coarseOperators - the coarse operator matrices (size of num_levels-1)
1899: Level: advanced
1901: Developer Note:
1902: Why is this here instead of in `PCMG` etc?
1904: .seealso: `PC`, `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetRScale()`, `PCMGGetInterpolation()`, `PCGetInterpolations()`
1905: @*/
1906: PetscErrorCode PCGetCoarseOperators(PC pc, PetscInt *num_levels, Mat *coarseOperators[])
1907: {
1911: PetscUseMethod(pc, "PCGetCoarseOperators_C", (PC, PetscInt *, Mat *[]), (pc, num_levels, coarseOperators));
1912: return 0;
1913: }