Actual source code: precon.c
1: /*
2: The PC (preconditioner) interface routines, callable by users.
3: */
4: #include <petsc/private/pcimpl.h>
5: #include <petscdm.h>
7: /* Logging support */
8: PetscClassId PC_CLASSID;
9: PetscLogEvent PC_SetUp, PC_SetUpOnBlocks, PC_Apply, PC_MatApply, PC_ApplyCoarse, PC_ApplySymmetricLeft;
10: PetscLogEvent PC_ApplySymmetricRight, PC_ModifySubMatrices, PC_ApplyOnBlocks, PC_ApplyTransposeOnBlocks;
11: PetscInt PetscMGLevelId;
12: PetscLogStage PCMPIStage;
14: PETSC_INTERN PetscErrorCode PCGetDefaultType_Private(PC pc, const char *type[])
15: {
16: PetscMPIInt size;
17: PetscBool hasopblock, hasopsolve, flg1, flg2, set, flg3, isnormal;
19: PetscFunctionBegin;
20: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
21: if (pc->pmat) {
22: PetscCall(MatHasOperation(pc->pmat, MATOP_GET_DIAGONAL_BLOCK, &hasopblock));
23: PetscCall(MatHasOperation(pc->pmat, MATOP_SOLVE, &hasopsolve));
24: if (size == 1) {
25: PetscCall(MatGetFactorAvailable(pc->pmat, "petsc", MAT_FACTOR_ICC, &flg1));
26: PetscCall(MatGetFactorAvailable(pc->pmat, "petsc", MAT_FACTOR_ILU, &flg2));
27: PetscCall(MatIsSymmetricKnown(pc->pmat, &set, &flg3));
28: PetscCall(PetscObjectTypeCompareAny((PetscObject)pc->pmat, &isnormal, MATNORMAL, MATNORMALHERMITIAN, NULL));
29: if (flg1 && (!flg2 || (set && flg3))) {
30: *type = PCICC;
31: } else if (flg2) {
32: *type = PCILU;
33: } else if (isnormal) {
34: *type = PCNONE;
35: } else if (hasopblock) { /* likely is a parallel matrix run on one processor */
36: if (pc->kspnestlevel > 0) {
37: Mat D;
39: PetscCall(MatGetDiagonalBlock(pc->pmat, &D));
40: PetscCall(PetscObjectTypeCompare((PetscObject)D, ((PetscObject)pc->pmat)->type_name, &flg1)); /* make sure there is no recursive call to PCGetDefaultType_Private() */
41: } else flg1 = PETSC_FALSE;
42: if (!flg1) *type = PCBJACOBI;
43: else *type = PCNONE;
44: } else if (hasopsolve) {
45: *type = PCMAT;
46: } else {
47: *type = PCNONE;
48: }
49: } else {
50: if (hasopblock) {
51: *type = PCBJACOBI;
52: } else if (hasopsolve) {
53: *type = PCMAT;
54: } else {
55: *type = PCNONE;
56: }
57: }
58: } else *type = NULL;
59: PetscFunctionReturn(PETSC_SUCCESS);
60: }
62: /* do not log solves, setup and applications of preconditioners while constructing preconditioners; perhaps they should be logged separately from the regular solves */
63: PETSC_EXTERN PetscLogEvent KSP_Solve, KSP_SetUp;
65: static PetscErrorCode PCLogEventsDeactivatePush(void)
66: {
67: PetscFunctionBegin;
68: PetscCall(KSPInitializePackage());
69: PetscCall(PetscLogEventDeactivatePush(KSP_Solve));
70: PetscCall(PetscLogEventDeactivatePush(KSP_SetUp));
71: PetscCall(PetscLogEventDeactivatePush(PC_Apply));
72: PetscCall(PetscLogEventDeactivatePush(PC_SetUp));
73: PetscCall(PetscLogEventDeactivatePush(PC_SetUpOnBlocks));
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
77: static PetscErrorCode PCLogEventsDeactivatePop(void)
78: {
79: PetscFunctionBegin;
80: PetscCall(KSPInitializePackage());
81: PetscCall(PetscLogEventDeactivatePop(KSP_Solve));
82: PetscCall(PetscLogEventDeactivatePop(KSP_SetUp));
83: PetscCall(PetscLogEventDeactivatePop(PC_Apply));
84: PetscCall(PetscLogEventDeactivatePop(PC_SetUp));
85: PetscCall(PetscLogEventDeactivatePop(PC_SetUpOnBlocks));
86: PetscFunctionReturn(PETSC_SUCCESS);
87: }
89: /*@
90: PCReset - Resets a `PC` context to the state it was in before `PCSetUp()` was called, and removes any allocated `Vec` and `Mat` from its data structure
92: Collective
94: Input Parameter:
95: . pc - the `PC` preconditioner context
97: Level: developer
99: Notes:
100: Any options set, including those set with `KSPSetFromOptions()` remain.
102: This allows a `PC` to be reused for a different sized linear system but using the same options that have been previously set in `pc`
104: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetUp()`
105: @*/
106: PetscErrorCode PCReset(PC pc)
107: {
108: PetscFunctionBegin;
110: PetscTryTypeMethod(pc, reset);
111: PetscCall(VecDestroy(&pc->diagonalscaleright));
112: PetscCall(VecDestroy(&pc->diagonalscaleleft));
113: PetscCall(MatDestroy(&pc->pmat));
114: PetscCall(MatDestroy(&pc->mat));
116: pc->setupcalled = PETSC_FALSE;
117: PetscFunctionReturn(PETSC_SUCCESS);
118: }
120: /*@
121: PCDestroy - Destroys `PC` context that was created with `PCCreate()`.
123: Collective
125: Input Parameter:
126: . pc - the `PC` preconditioner context
128: Level: developer
130: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetUp()`
131: @*/
132: PetscErrorCode PCDestroy(PC *pc)
133: {
134: PetscFunctionBegin;
135: if (!*pc) PetscFunctionReturn(PETSC_SUCCESS);
137: if (--((PetscObject)*pc)->refct > 0) {
138: *pc = NULL;
139: PetscFunctionReturn(PETSC_SUCCESS);
140: }
142: PetscCall(PCReset(*pc));
144: /* if memory was published with SAWs then destroy it */
145: PetscCall(PetscObjectSAWsViewOff((PetscObject)*pc));
146: PetscTryTypeMethod(*pc, destroy);
147: PetscCall(DMDestroy(&(*pc)->dm));
148: PetscCall(PetscHeaderDestroy(pc));
149: PetscFunctionReturn(PETSC_SUCCESS);
150: }
152: /*@
153: PCGetDiagonalScale - Indicates if the preconditioner applies an additional left and right
154: scaling as needed by certain time-stepping codes.
156: Logically Collective
158: Input Parameter:
159: . pc - the `PC` preconditioner context
161: Output Parameter:
162: . flag - `PETSC_TRUE` if it applies the scaling
164: Level: developer
166: Note:
167: If this returns `PETSC_TRUE` then the system solved via the Krylov method is, for left and right preconditioning,
169: $$
170: \begin{align*}
171: D M A D^{-1} y = D M b \\
172: D A M D^{-1} z = D b.
173: \end{align*}
174: $$
176: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCDiagonalScaleRight()`, `PCSetDiagonalScale()`
177: @*/
178: PetscErrorCode PCGetDiagonalScale(PC pc, PetscBool *flag)
179: {
180: PetscFunctionBegin;
182: PetscAssertPointer(flag, 2);
183: *flag = pc->diagonalscale;
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: /*@
188: PCSetDiagonalScale - Indicates the left scaling to use to apply an additional left and right
189: scaling as needed by certain time-stepping codes.
191: Logically Collective
193: Input Parameters:
194: + pc - the `PC` preconditioner context
195: - s - scaling vector
197: Level: intermediate
199: Notes:
200: The system solved via the Krylov method is, for left and right preconditioning,
201: $$
202: \begin{align*}
203: D M A D^{-1} y = D M b \\
204: D A M D^{-1} z = D b.
205: \end{align*}
206: $$
208: `PCDiagonalScaleLeft()` scales a vector by $D$. `PCDiagonalScaleRight()` scales a vector by $D^{-1}$.
210: .seealso: [](ch_ksp), `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCDiagonalScaleRight()`, `PCGetDiagonalScale()`
211: @*/
212: PetscErrorCode PCSetDiagonalScale(PC pc, Vec s)
213: {
214: PetscFunctionBegin;
217: pc->diagonalscale = PETSC_TRUE;
219: PetscCall(PetscObjectReference((PetscObject)s));
220: PetscCall(VecDestroy(&pc->diagonalscaleleft));
222: pc->diagonalscaleleft = s;
224: PetscCall(VecDuplicate(s, &pc->diagonalscaleright));
225: PetscCall(VecCopy(s, pc->diagonalscaleright));
226: PetscCall(VecReciprocal(pc->diagonalscaleright));
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
230: /*@
231: PCDiagonalScaleLeft - Scales a vector by the left scaling as needed by certain time-stepping codes.
233: Logically Collective
235: Input Parameters:
236: + pc - the `PC` preconditioner context
237: . in - input vector
238: - out - scaled vector (maybe the same as in)
240: Level: intermediate
242: Notes:
243: The system solved via the Krylov method is, for left and right preconditioning,
245: $$
246: \begin{align*}
247: D M A D^{-1} y = D M b \\
248: D A M D^{-1} z = D b.
249: \end{align*}
250: $$
252: `PCDiagonalScaleLeft()` scales a vector by $D$. `PCDiagonalScaleRight()` scales a vector by $D^{-1}$.
254: If diagonal scaling is turned off and `in` is not `out` then `in` is copied to `out`
256: .seealso: [](ch_ksp), `PCCreate()`, `PCSetUp()`, `PCSetDiagonalScale()`, `PCDiagonalScaleRight()`, `MatDiagonalScale()`
257: @*/
258: PetscErrorCode PCDiagonalScaleLeft(PC pc, Vec in, Vec out)
259: {
260: PetscFunctionBegin;
264: if (pc->diagonalscale) PetscCall(VecPointwiseMult(out, pc->diagonalscaleleft, in));
265: else if (in != out) PetscCall(VecCopy(in, out));
266: PetscFunctionReturn(PETSC_SUCCESS);
267: }
269: /*@
270: PCDiagonalScaleRight - Scales a vector by the right scaling as needed by certain time-stepping codes.
272: Logically Collective
274: Input Parameters:
275: + pc - the `PC` preconditioner context
276: . in - input vector
277: - out - scaled vector (maybe the same as in)
279: Level: intermediate
281: Notes:
282: The system solved via the Krylov method is, for left and right preconditioning,
284: $$
285: \begin{align*}
286: D M A D^{-1} y = D M b \\
287: D A M D^{-1} z = D b.
288: \end{align*}
289: $$
291: `PCDiagonalScaleLeft()` scales a vector by $D$. `PCDiagonalScaleRight()` scales a vector by $D^{-1}$.
293: If diagonal scaling is turned off and `in` is not `out` then `in` is copied to `out`
295: .seealso: [](ch_ksp), `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCSetDiagonalScale()`, `MatDiagonalScale()`
296: @*/
297: PetscErrorCode PCDiagonalScaleRight(PC pc, Vec in, Vec out)
298: {
299: PetscFunctionBegin;
303: if (pc->diagonalscale) {
304: PetscCall(VecPointwiseMult(out, pc->diagonalscaleright, in));
305: } else if (in != out) {
306: PetscCall(VecCopy(in, out));
307: }
308: PetscFunctionReturn(PETSC_SUCCESS);
309: }
311: /*@
312: PCSetUseAmat - Sets a flag to indicate that when the preconditioner needs to apply (part of) the
313: operator during the preconditioning process it applies the `Amat` provided to `TSSetRHSJacobian()`,
314: `TSSetIJacobian()`, `SNESSetJacobian()`, `KSPSetOperators()` or `PCSetOperators()` not the `Pmat`.
316: Logically Collective
318: Input Parameters:
319: + pc - the `PC` preconditioner context
320: - flg - `PETSC_TRUE` to use the `Amat`, `PETSC_FALSE` to use the `Pmat` (default is `PETSC_FALSE`)
322: Options Database Key:
323: . -pc_use_amat (true|false) - use the `Amat` argument to `KSPSetOperators()` or `PCSetOperators()` to apply the operator
325: Level: intermediate
327: Note:
328: For the common case in which the linear system matrix and the matrix used to construct the
329: preconditioner are identical, this routine has no effect.
331: .seealso: [](ch_ksp), `PC`, `PCGetUseAmat()`, `PCBJACOBI`, `PCMG`, `PCFIELDSPLIT`, `PCCOMPOSITE`,
332: `KSPSetOperators()`, `PCSetOperators()`
333: @*/
334: PetscErrorCode PCSetUseAmat(PC pc, PetscBool flg)
335: {
336: PetscFunctionBegin;
338: pc->useAmat = flg;
339: PetscFunctionReturn(PETSC_SUCCESS);
340: }
342: /*@
343: PCSetErrorIfFailure - Causes `PC` to generate an error if a floating point exception, for example a zero pivot, is detected.
345: Logically Collective
347: Input Parameters:
348: + pc - iterative context obtained from `PCCreate()`
349: - flg - `PETSC_TRUE` indicates you want the error generated
351: Level: advanced
353: Notes:
354: Normally PETSc continues if a linear solver fails due to a failed setup of a preconditioner, you can call `KSPGetConvergedReason()` after a `KSPSolve()`
355: 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
356: detected.
358: This is propagated into `KSP`s used by this `PC`, which then propagate it into `PC`s used by those `KSP`s
360: .seealso: [](ch_ksp), `PC`, `KSPSetErrorIfNotConverged()`, `PCGetInitialGuessNonzero()`, `PCSetInitialGuessKnoll()`, `PCGetInitialGuessKnoll()`
361: @*/
362: PetscErrorCode PCSetErrorIfFailure(PC pc, PetscBool flg)
363: {
364: PetscFunctionBegin;
367: pc->erroriffailure = flg;
368: PetscFunctionReturn(PETSC_SUCCESS);
369: }
371: /*@
372: PCGetUseAmat - Gets the flag that indicates that when the preconditioner needs to apply (part of) the
373: operator during the preconditioning process it applies the `Amat` provided to `TSSetRHSJacobian()`,
374: `TSSetIJacobian()`, `SNESSetJacobian()`, `KSPSetOperators()` or `PCSetOperators()` not the `Pmat`.
376: Logically Collective
378: Input Parameter:
379: . pc - the `PC` preconditioner context
381: Output Parameter:
382: . flg - `PETSC_TRUE` to use the `Amat`, `PETSC_FALSE` to use the `Pmat`
384: Level: intermediate
386: .seealso: [](ch_ksp), `PC`, `PCSetUseAmat()`, `PCBJACOBI`, `PCMG`, `PCFIELDSPLIT`, `PCCOMPOSITE`
387: @*/
388: PetscErrorCode PCGetUseAmat(PC pc, PetscBool *flg)
389: {
390: PetscFunctionBegin;
392: *flg = pc->useAmat;
393: PetscFunctionReturn(PETSC_SUCCESS);
394: }
396: /*@
397: PCSetKSPNestLevel - sets the amount of nesting the `KSP` that contains this `PC` has
399: Collective
401: Input Parameters:
402: + pc - the `PC`
403: - level - the nest level
405: Level: developer
407: .seealso: [](ch_ksp), `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSP`, `KSPGMRES`, `KSPType`, `KSPGetNestLevel()`, `PCGetKSPNestLevel()`, `KSPSetNestLevel()`
408: @*/
409: PetscErrorCode PCSetKSPNestLevel(PC pc, PetscInt level)
410: {
411: PetscFunctionBegin;
414: pc->kspnestlevel = level;
415: PetscFunctionReturn(PETSC_SUCCESS);
416: }
418: /*@
419: PCGetKSPNestLevel - gets the amount of nesting the `KSP` that contains this `PC` has
421: Not Collective
423: Input Parameter:
424: . pc - the `PC`
426: Output Parameter:
427: . level - the nest level
429: Level: developer
431: .seealso: [](ch_ksp), `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSP`, `KSPGMRES`, `KSPType`, `KSPSetNestLevel()`, `PCSetKSPNestLevel()`, `KSPGetNestLevel()`
432: @*/
433: PetscErrorCode PCGetKSPNestLevel(PC pc, PetscInt *level)
434: {
435: PetscFunctionBegin;
437: PetscAssertPointer(level, 2);
438: *level = pc->kspnestlevel;
439: PetscFunctionReturn(PETSC_SUCCESS);
440: }
442: /*@
443: PCCreate - Creates a preconditioner context, `PC`
445: Collective
447: Input Parameter:
448: . comm - MPI communicator
450: Output Parameter:
451: . newpc - location to put the `PC` preconditioner context
453: Level: developer
455: Notes:
456: This is rarely called directly by users since `KSP` manages the `PC` objects it uses. Use `KSPGetPC()` to access the `PC` used by a `KSP`.
458: Use `PCSetType()` or `PCSetFromOptions()` with the option `-pc_type pctype` to set the `PCType` for this `PC`
460: The default preconditioner type `PCType` for sparse matrices is `PCILU` or `PCICC` with 0 fill on one process and block Jacobi (`PCBJACOBI`) with `PCILU` or `PCICC`
461: in parallel. For dense matrices it is always `PCNONE`.
463: .seealso: [](ch_ksp), `PC`, `PCType`, `PCSetType`, `PCSetUp()`, `PCApply()`, `PCDestroy()`, `KSP`, `KSPGetPC()`
464: @*/
465: PetscErrorCode PCCreate(MPI_Comm comm, PC *newpc)
466: {
467: PC pc;
469: PetscFunctionBegin;
470: PetscAssertPointer(newpc, 2);
471: PetscCall(PCInitializePackage());
473: PetscCall(PetscHeaderCreate(pc, PC_CLASSID, "PC", "Preconditioner", "PC", comm, PCDestroy, PCView));
474: PetscCall(PCParametersInitialize(pc));
475: *newpc = pc;
476: PetscFunctionReturn(PETSC_SUCCESS);
477: }
479: /*@
480: PCParametersInitialize - Sets the base defaults for parameters in `pc`, updating a parameter's current value when it matches its previously recorded default.
482: Logically collective
484: Input Parameter:
485: . pc - the `PC` object
487: Level: developer
489: Notes:
491: The base defaults are the non-type-specific values established when the `PC` is created. A `PCType` constructor may subsequently replace them with type-specific defaults.
493: Developer Notes:
495: `PCCreate()` calls this routine to establish the base defaults. `PCSetType()` calls it before constructing a new `PCType`, so the recorded defaults associated with the previous type are replaced before the new type installs its own defaults.
497: Default tracking is based on value equality, not on whether a setter was called. Consequently, an explicitly assigned value that equals the recorded default may be updated when the type changes.
499: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCDestroy()`, `PetscObjectParameterSetDefault()`
500: @*/
501: PetscErrorCode PCParametersInitialize(PC pc)
502: {
503: PetscObjectParameterSetDefault(pc, useAmat, PETSC_FALSE);
504: return PETSC_SUCCESS;
505: }
507: /*@
508: PCApply - Applies the preconditioner to a vector.
510: Collective
512: Input Parameters:
513: + pc - the `PC` preconditioner context
514: - x - input vector
516: Output Parameter:
517: . y - output vector
519: Level: developer
521: .seealso: [](ch_ksp), `PC`, `PCApplyTranspose()`, `PCApplyBAorAB()`
522: @*/
523: PetscErrorCode PCApply(PC pc, Vec x, Vec y)
524: {
525: PetscInt m, n, mv, nv;
527: PetscFunctionBegin;
531: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
532: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
533: /* use pmat to check vector sizes since for KSPLSQR the pmat may be of a different size than mat */
534: PetscCall(MatGetLocalSize(pc->pmat, &m, &n));
535: PetscCall(VecGetLocalSize(x, &mv));
536: PetscCall(VecGetLocalSize(y, &nv));
537: /* check pmat * y = x is feasible */
538: PetscCheck(mv == m, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Preconditioner number of local rows %" PetscInt_FMT " does not equal input vector size %" PetscInt_FMT, m, mv);
539: PetscCheck(nv == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Preconditioner number of local columns %" PetscInt_FMT " does not equal output vector size %" PetscInt_FMT, n, nv);
540: PetscCall(VecSetErrorIfLocked(y, 3));
542: PetscCall(PCSetUp(pc));
543: PetscCall(VecLockReadPush(x));
544: PetscCall(PetscLogEventBegin(PC_Apply, pc, x, y, 0));
545: PetscUseTypeMethod(pc, apply, x, y);
546: PetscCall(PetscLogEventEnd(PC_Apply, pc, x, y, 0));
547: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
548: PetscCall(VecLockReadPop(x));
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }
552: static PetscErrorCode PCMatApplyTranspose_Private(PC pc, Mat X, Mat Y, PetscBool transpose)
553: {
554: Mat A;
555: Vec cy, cx;
556: PetscInt m1, M1, m2, M2, n1, N1, n2, N2, m3, M3, n3, N3;
557: PetscBool match;
559: PetscFunctionBegin;
563: PetscCheckSameComm(pc, 1, X, 2);
564: PetscCheckSameComm(pc, 1, Y, 3);
565: PetscCheck(Y != X, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "Y and X must be different matrices");
566: PetscCall(PCGetOperators(pc, NULL, &A));
567: PetscCall(MatGetLocalSize(A, &m3, &n3));
568: PetscCall(MatGetLocalSize(X, &m2, &n2));
569: PetscCall(MatGetLocalSize(Y, &m1, &n1));
570: PetscCall(MatGetSize(A, &M3, &N3));
571: PetscCall(MatGetSize(X, &M2, &N2));
572: PetscCall(MatGetSize(Y, &M1, &N1));
573: PetscCheck(n1 == n2 && N1 == N2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible number of columns between block of input vectors (n,N) = (%" PetscInt_FMT ",%" PetscInt_FMT ") and block of output vectors (n,N) = (%" PetscInt_FMT ",%" PetscInt_FMT ")", n2, N2, n1, N1);
574: PetscCheck(m2 == m3 && M2 == M3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible layout between block of input vectors (m,M) = (%" PetscInt_FMT ",%" PetscInt_FMT ") and Pmat (m,M)x(n,N) = (%" PetscInt_FMT ",%" PetscInt_FMT ")x(%" PetscInt_FMT ",%" PetscInt_FMT ")", m2, M2, m3, M3, n3, N3);
575: PetscCheck(m1 == n3 && M1 == N3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible layout between block of output vectors (m,M) = (%" PetscInt_FMT ",%" PetscInt_FMT ") and Pmat (m,M)x(n,N) = (%" PetscInt_FMT ",%" PetscInt_FMT ")x(%" PetscInt_FMT ",%" PetscInt_FMT ")", m1, M1, m3, M3, n3, N3);
576: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)Y, &match, MATSEQDENSE, MATMPIDENSE, ""));
577: PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided block of output vectors not stored in a dense Mat");
578: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)X, &match, MATSEQDENSE, MATMPIDENSE, ""));
579: PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided block of input vectors not stored in a dense Mat");
580: PetscCall(PCSetUp(pc));
581: if (!transpose && pc->ops->matapply) {
582: PetscCall(PetscLogEventBegin(PC_MatApply, pc, X, Y, 0));
583: PetscUseTypeMethod(pc, matapply, X, Y);
584: PetscCall(PetscLogEventEnd(PC_MatApply, pc, X, Y, 0));
585: } else if (transpose && pc->ops->matapplytranspose) {
586: PetscCall(PetscLogEventBegin(PC_MatApply, pc, X, Y, 0));
587: PetscUseTypeMethod(pc, matapplytranspose, X, Y);
588: PetscCall(PetscLogEventEnd(PC_MatApply, pc, X, Y, 0));
589: } else {
590: PetscCall(PetscInfo(pc, "PC type %s applying column by column\n", ((PetscObject)pc)->type_name));
591: for (n1 = 0; n1 < N1; ++n1) {
592: PetscCall(MatDenseGetColumnVecRead(X, n1, &cx));
593: PetscCall(MatDenseGetColumnVecWrite(Y, n1, &cy));
594: if (!transpose) PetscCall(PCApply(pc, cx, cy));
595: else PetscCall(PCApplyTranspose(pc, cx, cy));
596: PetscCall(MatDenseRestoreColumnVecWrite(Y, n1, &cy));
597: PetscCall(MatDenseRestoreColumnVecRead(X, n1, &cx));
598: }
599: }
600: PetscFunctionReturn(PETSC_SUCCESS);
601: }
603: /*@
604: PCMatApply - Applies the preconditioner to multiple vectors stored as a `MATDENSE`. Like `PCApply()`, `Y` and `X` must be different matrices.
606: Collective
608: Input Parameters:
609: + pc - the `PC` preconditioner context
610: - X - block of input vectors
612: Output Parameter:
613: . Y - block of output vectors
615: Level: developer
617: .seealso: [](ch_ksp), `PC`, `PCApply()`, `KSPMatSolve()`
618: @*/
619: PetscErrorCode PCMatApply(PC pc, Mat X, Mat Y)
620: {
621: PetscFunctionBegin;
622: PetscCall(PCMatApplyTranspose_Private(pc, X, Y, PETSC_FALSE));
623: PetscFunctionReturn(PETSC_SUCCESS);
624: }
626: /*@
627: PCMatApplyTranspose - Applies the transpose of preconditioner to multiple vectors stored as a `MATDENSE`. Like `PCApplyTranspose()`, `Y` and `X` must be different matrices.
629: Collective
631: Input Parameters:
632: + pc - the `PC` preconditioner context
633: - X - block of input vectors
635: Output Parameter:
636: . Y - block of output vectors
638: Level: developer
640: .seealso: [](ch_ksp), `PC`, `PCApplyTranspose()`, `KSPMatSolveTranspose()`
641: @*/
642: PetscErrorCode PCMatApplyTranspose(PC pc, Mat X, Mat Y)
643: {
644: PetscFunctionBegin;
645: PetscCall(PCMatApplyTranspose_Private(pc, X, Y, PETSC_TRUE));
646: PetscFunctionReturn(PETSC_SUCCESS);
647: }
649: /*@
650: PCApplySymmetricLeft - Applies the left part of a symmetric preconditioner to a vector.
652: Collective
654: Input Parameters:
655: + pc - the `PC` preconditioner context
656: - x - input vector
658: Output Parameter:
659: . y - output vector
661: Level: developer
663: Note:
664: Currently, this routine is implemented only for `PCICC` and `PCJACOBI` preconditioners.
666: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplySymmetricRight()`
667: @*/
668: PetscErrorCode PCApplySymmetricLeft(PC pc, Vec x, Vec y)
669: {
670: PetscFunctionBegin;
674: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
675: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
676: PetscCall(PCSetUp(pc));
677: PetscCall(VecLockReadPush(x));
678: PetscCall(PetscLogEventBegin(PC_ApplySymmetricLeft, pc, x, y, 0));
679: PetscUseTypeMethod(pc, applysymmetricleft, x, y);
680: PetscCall(PetscLogEventEnd(PC_ApplySymmetricLeft, pc, x, y, 0));
681: PetscCall(VecLockReadPop(x));
682: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
683: PetscFunctionReturn(PETSC_SUCCESS);
684: }
686: /*@
687: PCApplySymmetricRight - Applies the right part of a symmetric preconditioner to a vector.
689: Collective
691: Input Parameters:
692: + pc - the `PC` preconditioner context
693: - x - input vector
695: Output Parameter:
696: . y - output vector
698: Level: developer
700: Note:
701: Currently, this routine is implemented only for `PCICC` and `PCJACOBI` preconditioners.
703: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplySymmetricLeft()`
704: @*/
705: PetscErrorCode PCApplySymmetricRight(PC pc, Vec x, Vec y)
706: {
707: PetscFunctionBegin;
711: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
712: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
713: PetscCall(PCSetUp(pc));
714: PetscCall(VecLockReadPush(x));
715: PetscCall(PetscLogEventBegin(PC_ApplySymmetricRight, pc, x, y, 0));
716: PetscUseTypeMethod(pc, applysymmetricright, x, y);
717: PetscCall(PetscLogEventEnd(PC_ApplySymmetricRight, pc, x, y, 0));
718: PetscCall(VecLockReadPop(x));
719: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
720: PetscFunctionReturn(PETSC_SUCCESS);
721: }
723: /*@
724: PCApplyTranspose - Applies the transpose of preconditioner to a vector.
726: Collective
728: Input Parameters:
729: + pc - the `PC` preconditioner context
730: - x - input vector
732: Output Parameter:
733: . y - output vector
735: Level: developer
737: Note:
738: For complex numbers this applies the non-Hermitian transpose.
740: Developer Note:
741: We need to implement a `PCApplyHermitianTranspose()`
743: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplyBAorAB()`, `PCApplyBAorABTranspose()`, `PCApplyTransposeExists()`
744: @*/
745: PetscErrorCode PCApplyTranspose(PC pc, Vec x, Vec y)
746: {
747: PetscFunctionBegin;
751: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
752: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
753: PetscCall(PCSetUp(pc));
754: PetscCall(VecLockReadPush(x));
755: PetscCall(PetscLogEventBegin(PC_Apply, pc, x, y, 0));
756: PetscUseTypeMethod(pc, applytranspose, x, y);
757: PetscCall(PetscLogEventEnd(PC_Apply, pc, x, y, 0));
758: PetscCall(VecLockReadPop(x));
759: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
760: PetscFunctionReturn(PETSC_SUCCESS);
761: }
763: /*@
764: PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation
766: Collective
768: Input Parameter:
769: . pc - the `PC` preconditioner context
771: Output Parameter:
772: . flg - `PETSC_TRUE` if a transpose operation is defined
774: Level: developer
776: .seealso: [](ch_ksp), `PC`, `PCApplyTranspose()`
777: @*/
778: PetscErrorCode PCApplyTransposeExists(PC pc, PetscBool *flg)
779: {
780: PetscFunctionBegin;
782: PetscAssertPointer(flg, 2);
783: if (pc->ops->applytranspose) *flg = PETSC_TRUE;
784: else *flg = PETSC_FALSE;
785: PetscFunctionReturn(PETSC_SUCCESS);
786: }
788: /*@
789: PCApplyBAorAB - Applies the preconditioner and operator to a vector. $y = B*A*x $ or $ y = A*B*x$.
791: Collective
793: Input Parameters:
794: + pc - the `PC` preconditioner context
795: . side - indicates the preconditioner side, one of `PC_LEFT`, `PC_RIGHT`, or `PC_SYMMETRIC`
796: . x - input vector
797: - work - work vector
799: Output Parameter:
800: . y - output vector
802: Level: developer
804: Note:
805: 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.
806: The specific `KSPSolve()` method must also be written to handle the post-solve "correction" for the diagonal scaling.
808: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplyTranspose()`, `PCApplyBAorABTranspose()`
809: @*/
810: PetscErrorCode PCApplyBAorAB(PC pc, PCSide side, Vec x, Vec y, Vec work)
811: {
812: PetscFunctionBegin;
818: PetscCheckSameComm(pc, 1, x, 3);
819: PetscCheckSameComm(pc, 1, y, 4);
820: PetscCheckSameComm(pc, 1, work, 5);
821: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
822: PetscCheck(side == PC_LEFT || side == PC_SYMMETRIC || side == PC_RIGHT, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Side must be right, left, or symmetric");
823: PetscCheck(!pc->diagonalscale || side != PC_SYMMETRIC, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot include diagonal scaling with symmetric preconditioner application");
824: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 3, PETSC_TRUE));
826: PetscCall(PCSetUp(pc));
827: if (pc->diagonalscale) {
828: if (pc->ops->applyBA) {
829: Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
830: PetscCall(VecDuplicate(x, &work2));
831: PetscCall(PCDiagonalScaleRight(pc, x, work2));
832: PetscUseTypeMethod(pc, applyBA, side, work2, y, work);
833: PetscCall(PCDiagonalScaleLeft(pc, y, y));
834: PetscCall(VecDestroy(&work2));
835: } else if (side == PC_RIGHT) {
836: PetscCall(PCDiagonalScaleRight(pc, x, y));
837: PetscCall(PCApply(pc, y, work));
838: PetscCall(MatMult(pc->mat, work, y));
839: PetscCall(PCDiagonalScaleLeft(pc, y, y));
840: } else if (side == PC_LEFT) {
841: PetscCall(PCDiagonalScaleRight(pc, x, y));
842: PetscCall(MatMult(pc->mat, y, work));
843: PetscCall(PCApply(pc, work, y));
844: PetscCall(PCDiagonalScaleLeft(pc, y, y));
845: } else PetscCheck(side != PC_SYMMETRIC, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot provide diagonal scaling with symmetric application of preconditioner");
846: } else {
847: if (pc->ops->applyBA) {
848: PetscUseTypeMethod(pc, applyBA, side, x, y, work);
849: } else if (side == PC_RIGHT) {
850: PetscCall(PCApply(pc, x, work));
851: PetscCall(MatMult(pc->mat, work, y));
852: } else if (side == PC_LEFT) {
853: PetscCall(MatMult(pc->mat, x, work));
854: PetscCall(PCApply(pc, work, y));
855: } else if (side == PC_SYMMETRIC) {
856: /* There's an extra copy here; maybe should provide 2 work vectors instead? */
857: PetscCall(PCApplySymmetricRight(pc, x, work));
858: PetscCall(MatMult(pc->mat, work, y));
859: PetscCall(VecCopy(y, work));
860: PetscCall(PCApplySymmetricLeft(pc, work, y));
861: }
862: }
863: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 4, PETSC_FALSE));
864: PetscFunctionReturn(PETSC_SUCCESS);
865: }
867: /*@
868: PCApplyBAorABTranspose - Applies the transpose of the preconditioner
869: and operator to a vector. That is, applies $B^T * A^T$ with left preconditioning,
870: NOT $(B*A)^T = A^T*B^T$.
872: Collective
874: Input Parameters:
875: + pc - the `PC` preconditioner context
876: . side - indicates the preconditioner side, one of `PC_LEFT`, `PC_RIGHT`, or `PC_SYMMETRIC`
877: . x - input vector
878: - work - work vector
880: Output Parameter:
881: . y - output vector
883: Level: developer
885: Note:
886: This routine is used internally so that the same Krylov code can be used to solve $A x = b$ and $A^T x = b$, with a preconditioner
887: defined by $B^T$. This is why this has the funny form that it computes $B^T * A^T$
889: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplyTranspose()`, `PCApplyBAorAB()`
890: @*/
891: PetscErrorCode PCApplyBAorABTranspose(PC pc, PCSide side, Vec x, Vec y, Vec work)
892: {
893: PetscFunctionBegin;
898: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
899: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 3, PETSC_TRUE));
900: if (pc->ops->applyBAtranspose) {
901: PetscUseTypeMethod(pc, applyBAtranspose, side, x, y, work);
902: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 4, PETSC_FALSE));
903: PetscFunctionReturn(PETSC_SUCCESS);
904: }
905: PetscCheck(side == PC_LEFT || side == PC_RIGHT, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Side must be right or left");
907: PetscCall(PCSetUp(pc));
908: if (side == PC_RIGHT) {
909: PetscCall(PCApplyTranspose(pc, x, work));
910: PetscCall(MatMultTranspose(pc->mat, work, y));
911: } else if (side == PC_LEFT) {
912: PetscCall(MatMultTranspose(pc->mat, x, work));
913: PetscCall(PCApplyTranspose(pc, work, y));
914: }
915: /* add support for PC_SYMMETRIC */
916: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 4, PETSC_FALSE));
917: PetscFunctionReturn(PETSC_SUCCESS);
918: }
920: /*@
921: PCApplyRichardsonExists - Determines whether a particular preconditioner has a
922: built-in fast application of Richardson's method.
924: Not Collective
926: Input Parameter:
927: . pc - the preconditioner
929: Output Parameter:
930: . exists - `PETSC_TRUE` or `PETSC_FALSE`
932: Level: developer
934: .seealso: [](ch_ksp), `PC`, `KSPRICHARDSON`, `PCApplyRichardson()`
935: @*/
936: PetscErrorCode PCApplyRichardsonExists(PC pc, PetscBool *exists)
937: {
938: PetscFunctionBegin;
940: PetscAssertPointer(exists, 2);
941: if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
942: else *exists = PETSC_FALSE;
943: PetscFunctionReturn(PETSC_SUCCESS);
944: }
946: /*@
947: PCApplyRichardson - Applies several steps of Richardson iteration with
948: the particular preconditioner. This routine is usually used by the
949: Krylov solvers and not the application code directly.
951: Collective
953: Input Parameters:
954: + pc - the `PC` preconditioner context
955: . b - the right-hand side
956: . w - one work vector
957: . rtol - relative decrease in residual norm convergence criteria
958: . abstol - absolute residual norm convergence criteria
959: . dtol - divergence residual norm increase criteria
960: . its - the number of iterations to apply.
961: - guesszero - if the input x contains nonzero initial guess
963: Output Parameters:
964: + outits - number of iterations actually used (for SOR this always equals its)
965: . reason - the reason the apply terminated
966: - y - the solution (also contains initial guess if guesszero is `PETSC_FALSE`
968: Level: developer
970: Notes:
971: Most preconditioners do not support this function. Use the command
972: `PCApplyRichardsonExists()` to determine if one does.
974: Except for the `PCMG` this routine ignores the convergence tolerances
975: and always runs for the number of iterations
977: .seealso: [](ch_ksp), `PC`, `PCApplyRichardsonExists()`
978: @*/
979: PetscErrorCode PCApplyRichardson(PC pc, Vec b, Vec y, Vec w, PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt its, PetscBool guesszero, PetscInt *outits, PCRichardsonConvergedReason *reason)
980: {
981: PetscFunctionBegin;
986: PetscCheck(b != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "b and y must be different vectors");
987: PetscCall(PCSetUp(pc));
988: PetscUseTypeMethod(pc, applyrichardson, b, y, w, rtol, abstol, dtol, its, guesszero, outits, reason);
989: PetscFunctionReturn(PETSC_SUCCESS);
990: }
992: /*@
993: PCSetFailedReason - Sets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail
995: Logically Collective
997: Input Parameters:
998: + pc - the `PC` preconditioner context
999: - reason - the reason it failed
1001: Level: advanced
1003: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCFailedReason`
1004: @*/
1005: PetscErrorCode PCSetFailedReason(PC pc, PCFailedReason reason)
1006: {
1007: PetscFunctionBegin;
1009: pc->failedreason = reason;
1010: PetscFunctionReturn(PETSC_SUCCESS);
1011: }
1013: /*@
1014: PCGetFailedReason - Gets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail
1016: Not Collective
1018: Input Parameter:
1019: . pc - the `PC` preconditioner context
1021: Output Parameter:
1022: . reason - the reason it failed
1024: Level: advanced
1026: Note:
1027: After a call to `KSPCheckDot()` or `KSPCheckNorm()` inside a `KSPSolve()` or a call to `PCReduceFailedReason()`
1028: this is the maximum reason over all MPI processes in the `PC` communicator and hence logically collective.
1029: Otherwise it returns the local value.
1031: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCSetFailedReason()`, `PCFailedReason`
1032: @*/
1033: PetscErrorCode PCGetFailedReason(PC pc, PCFailedReason *reason)
1034: {
1035: PetscFunctionBegin;
1037: *reason = pc->failedreason;
1038: PetscFunctionReturn(PETSC_SUCCESS);
1039: }
1041: /*@
1042: PCReduceFailedReason - Reduce the failed reason among the MPI processes that share the `PC`
1044: Collective
1046: Input Parameter:
1047: . pc - the `PC` preconditioner context
1049: Level: advanced
1051: Note:
1052: Different MPI processes may have different reasons or no reason, see `PCGetFailedReason()`. This routine
1053: makes them have a common value (failure if any MPI process had a failure).
1055: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCGetFailedReason()`, `PCSetFailedReason()`, `PCFailedReason`
1056: @*/
1057: PetscErrorCode PCReduceFailedReason(PC pc)
1058: {
1059: PetscInt buf;
1061: PetscFunctionBegin;
1063: buf = (PetscInt)pc->failedreason;
1064: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &buf, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)pc)));
1065: pc->failedreason = (PCFailedReason)buf;
1066: PetscFunctionReturn(PETSC_SUCCESS);
1067: }
1069: /*
1070: a setupcall of 0 indicates never setup,
1071: 1 indicates has been previously setup
1072: -1 indicates a PCSetUp() was attempted and failed
1073: */
1074: /*@
1075: PCSetUp - Prepares for the use of a preconditioner. Performs all the one-time operations needed before the preconditioner
1076: can be used with `PCApply()`
1078: Collective
1080: Input Parameter:
1081: . pc - the `PC` preconditioner context
1083: Level: developer
1085: Notes:
1086: For example, for `PCLU` this will compute the factorization.
1088: This is called automatically by `KSPSetUp()` or `PCApply()` so rarely needs to be called directly.
1090: For nested preconditioners, such as `PCFIELDSPLIT` or `PCBJACOBI` this may not finish the construction of the preconditioner
1091: on the inner levels, the routine `PCSetUpOnBlocks()` may compute more of the preconditioner in those situations.
1093: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `KSPSetUp()`, `PCSetUpOnBlocks()`
1094: @*/
1095: PetscErrorCode PCSetUp(PC pc)
1096: {
1097: const char *def;
1098: PetscObjectState matstate, matnonzerostate;
1100: PetscFunctionBegin;
1102: PetscCheck(pc->mat, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Missing A matrix");
1103: PetscCheck(pc->pmat, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Missing P matrix");
1105: if (pc->setupcalled && pc->reusepreconditioner) {
1106: PetscCall(PetscInfo(pc, "Leaving PC with identical preconditioner since reuse preconditioner is set\n"));
1107: PetscFunctionReturn(PETSC_SUCCESS);
1108: }
1110: PetscCall(PetscObjectStateGet((PetscObject)pc->pmat, &matstate));
1111: PetscCall(MatGetNonzeroState(pc->pmat, &matnonzerostate));
1112: if (!pc->setupcalled) {
1113: //PetscCall(PetscInfo(pc, "Setting up PC for first time\n"));
1114: pc->flag = DIFFERENT_NONZERO_PATTERN;
1115: } else if (matstate == pc->matstate) PetscFunctionReturn(PETSC_SUCCESS);
1116: else {
1117: if (matnonzerostate != pc->matnonzerostate) {
1118: PetscCall(PetscInfo(pc, "Setting up PC with different nonzero pattern\n"));
1119: pc->flag = DIFFERENT_NONZERO_PATTERN;
1120: } else {
1121: //PetscCall(PetscInfo(pc, "Setting up PC with same nonzero pattern\n"));
1122: pc->flag = SAME_NONZERO_PATTERN;
1123: }
1124: }
1125: pc->matstate = matstate;
1126: pc->matnonzerostate = matnonzerostate;
1128: if (!((PetscObject)pc)->type_name) {
1129: PetscCall(PCGetDefaultType_Private(pc, &def));
1130: PetscCall(PCSetType(pc, def));
1131: }
1133: PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
1134: PetscCall(MatSetErrorIfFailure(pc->mat, pc->erroriffailure));
1135: PetscCall(PetscLogEventBegin(PC_SetUp, pc, 0, 0, 0));
1136: if (pc->ops->setup) {
1137: PetscCall(PCLogEventsDeactivatePush());
1138: PetscUseTypeMethod(pc, setup);
1139: PetscCall(PCLogEventsDeactivatePop());
1140: }
1141: PetscCall(PetscLogEventEnd(PC_SetUp, pc, 0, 0, 0));
1142: if (pc->postsetup) PetscCall((*pc->postsetup)(pc));
1143: if (!pc->setupcalled) pc->setupcalled = PETSC_TRUE;
1144: PetscFunctionReturn(PETSC_SUCCESS);
1145: }
1147: /*@
1148: PCSetUpOnBlocks - Sets up the preconditioner for each block in
1149: the block Jacobi, overlapping Schwarz, and fieldsplit methods.
1151: Collective
1153: Input Parameter:
1154: . pc - the `PC` preconditioner context
1156: Level: developer
1158: Notes:
1159: For nested preconditioners such as `PCBJACOBI`, `PCSetUp()` is not called on each sub-`KSP` when `PCSetUp()` is
1160: called on the outer `PC`, this routine ensures it is called.
1162: It calls `PCSetUp()` if not yet called.
1164: .seealso: [](ch_ksp), `PC`, `PCSetUp()`, `PCCreate()`, `PCApply()`, `PCDestroy()`
1165: @*/
1166: PetscErrorCode PCSetUpOnBlocks(PC pc)
1167: {
1168: PetscFunctionBegin;
1170: if (!pc->setupcalled) PetscCall(PCSetUp(pc)); /* "if" to prevent -info extra prints */
1171: if (!pc->ops->setuponblocks) PetscFunctionReturn(PETSC_SUCCESS);
1172: PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
1173: PetscCall(PetscLogEventBegin(PC_SetUpOnBlocks, pc, 0, 0, 0));
1174: PetscCall(PCLogEventsDeactivatePush());
1175: PetscUseTypeMethod(pc, setuponblocks);
1176: PetscCall(PCLogEventsDeactivatePop());
1177: PetscCall(PetscLogEventEnd(PC_SetUpOnBlocks, pc, 0, 0, 0));
1178: PetscFunctionReturn(PETSC_SUCCESS);
1179: }
1181: /*@C
1182: PCSetModifySubMatrices - Sets a user-defined routine for modifying the
1183: submatrices that arise within certain subdomain-based preconditioners such as `PCASM`
1185: Logically Collective
1187: Input Parameters:
1188: + pc - the `PC` preconditioner context
1189: . func - routine for modifying the submatrices, see `PCModifySubMatricesFn`
1190: - ctx - optional user-defined context (may be `NULL`)
1192: Level: advanced
1194: Notes:
1195: The basic submatrices are extracted from the matrix used to construct the preconditioner as
1196: usual; the user can then alter these (for example, to set different boundary
1197: conditions for each submatrix) before they are used for the local solves.
1199: `PCSetModifySubMatrices()` MUST be called before `KSPSetUp()` and
1200: `KSPSolve()`.
1202: A routine set by `PCSetModifySubMatrices()` is currently called within
1203: `PCBJACOBI`, `PCASM`, `PCGASM`, and `PCHPDDM`.
1204: All other preconditioners ignore this routine.
1206: .seealso: [](ch_ksp), `PC`, `PCModifySubMatricesFn`, `PCBJACOBI`, `PCASM`, `PCModifySubMatrices()`
1207: @*/
1208: PetscErrorCode PCSetModifySubMatrices(PC pc, PCModifySubMatricesFn *func, PetscCtx ctx)
1209: {
1210: PetscFunctionBegin;
1212: pc->modifysubmatrices = func;
1213: pc->modifysubmatricesP = ctx;
1214: PetscFunctionReturn(PETSC_SUCCESS);
1215: }
1217: /*@C
1218: PCModifySubMatrices - Calls an optional user-defined routine within
1219: certain preconditioners if one has been set with `PCSetModifySubMatrices()`.
1221: Collective
1223: Input Parameters:
1224: + pc - the `PC` preconditioner context
1225: . nsub - the number of local submatrices
1226: . row - an array of index sets that contain the global row numbers
1227: that comprise each local submatrix
1228: . col - an array of index sets that contain the global column numbers
1229: that comprise each local submatrix
1230: . submat - array of local submatrices
1231: - ctx - optional user-defined context for private data for the
1232: user-defined routine (may be `NULL`)
1234: Output Parameter:
1235: . submat - array of local submatrices (the entries of which may
1236: have been modified)
1238: Level: developer
1240: Note:
1241: The user should NOT generally call this routine, as it will
1242: automatically be called within certain preconditioners.
1244: .seealso: [](ch_ksp), `PC`, `PCModifySubMatricesFn`, `PCSetModifySubMatrices()`
1245: @*/
1246: PetscErrorCode PCModifySubMatrices(PC pc, PetscInt nsub, const IS row[], const IS col[], Mat submat[], PetscCtx ctx)
1247: {
1248: PetscFunctionBegin;
1250: if (!pc->modifysubmatrices) PetscFunctionReturn(PETSC_SUCCESS);
1251: PetscCall(PetscLogEventBegin(PC_ModifySubMatrices, pc, 0, 0, 0));
1252: PetscCall((*pc->modifysubmatrices)(pc, nsub, row, col, submat, ctx));
1253: PetscCall(PetscLogEventEnd(PC_ModifySubMatrices, pc, 0, 0, 0));
1254: PetscFunctionReturn(PETSC_SUCCESS);
1255: }
1257: /*@
1258: PCSetOperators - Sets the matrix associated with the linear system and
1259: a (possibly) different one from which the preconditioner will be constructed.
1261: Logically Collective
1263: Input Parameters:
1264: + pc - the `PC` preconditioner context
1265: . Amat - the matrix that defines the linear system
1266: - Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
1268: Level: advanced
1270: Notes:
1271: Using this routine directly is rarely needed, the preferred, and equivalent, usage is `KSPSetOperators()`.
1273: Passing a `NULL` for `Amat` or `Pmat` removes the matrix that is currently used.
1275: If you wish to replace either `Amat` or `Pmat` but leave the other one untouched then
1276: first call `KSPGetOperators()` to get the one you wish to keep, call `PetscObjectReference()`
1277: on it and then pass it back in in your call to `KSPSetOperators()`.
1279: More Notes about Repeated Solution of Linear Systems:
1280: PETSc does NOT reset the matrix entries of either `Amat` or `Pmat`
1281: to zero after a linear solve; the user is completely responsible for
1282: matrix assembly. See the routine `MatZeroEntries()` if desiring to
1283: zero all elements of a matrix.
1285: .seealso: [](ch_ksp), `PC`, `PCGetOperators()`, `MatZeroEntries()`
1286: @*/
1287: PetscErrorCode PCSetOperators(PC pc, Mat Amat, Mat Pmat)
1288: {
1289: PetscInt m1, n1, m2, n2;
1291: PetscFunctionBegin;
1295: if (Amat) PetscCheckSameComm(pc, 1, Amat, 2);
1296: if (Pmat) PetscCheckSameComm(pc, 1, Pmat, 3);
1297: if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
1298: PetscCall(MatGetLocalSize(Amat, &m1, &n1));
1299: PetscCall(MatGetLocalSize(pc->mat, &m2, &n2));
1300: PetscCheck(m1 == m2 && n1 == n2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Cannot change local size of Amat after use old sizes %" PetscInt_FMT " %" PetscInt_FMT " new sizes %" PetscInt_FMT " %" PetscInt_FMT, m2, n2, m1, n1);
1301: PetscCall(MatGetLocalSize(Pmat, &m1, &n1));
1302: PetscCall(MatGetLocalSize(pc->pmat, &m2, &n2));
1303: PetscCheck(m1 == m2 && n1 == n2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Cannot change local size of Pmat after use old sizes %" PetscInt_FMT " %" PetscInt_FMT " new sizes %" PetscInt_FMT " %" PetscInt_FMT, m2, n2, m1, n1);
1304: }
1306: if (Pmat != pc->pmat) {
1307: /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
1308: pc->matnonzerostate = -1;
1309: pc->matstate = -1;
1310: }
1312: /* reference first in case the matrices are the same */
1313: PetscCall(PetscObjectReference((PetscObject)Amat));
1314: PetscCall(MatDestroy(&pc->mat));
1315: PetscCall(PetscObjectReference((PetscObject)Pmat));
1316: PetscCall(MatDestroy(&pc->pmat));
1317: pc->mat = Amat;
1318: pc->pmat = Pmat;
1319: PetscFunctionReturn(PETSC_SUCCESS);
1320: }
1322: /*@
1323: PCSetReusePreconditioner - reuse the current preconditioner even if the operator in the preconditioner `PC` has changed.
1325: Logically Collective
1327: Input Parameters:
1328: + pc - the `PC` preconditioner context
1329: - flag - `PETSC_TRUE` do not compute a new preconditioner, `PETSC_FALSE` do compute a new preconditioner
1331: Level: intermediate
1333: Note:
1334: Normally if a matrix inside a `PC` changes the `PC` automatically updates itself using information from the changed matrix. This option
1335: prevents this.
1337: .seealso: [](ch_ksp), `PC`, `PCGetOperators()`, `MatZeroEntries()`, `PCGetReusePreconditioner()`, `KSPSetReusePreconditioner()`
1338: @*/
1339: PetscErrorCode PCSetReusePreconditioner(PC pc, PetscBool flag)
1340: {
1341: PetscFunctionBegin;
1344: pc->reusepreconditioner = flag;
1345: PetscTryMethod(pc, "PCSetReusePreconditioner_C", (PC, PetscBool), (pc, flag));
1346: PetscFunctionReturn(PETSC_SUCCESS);
1347: }
1349: /*@
1350: PCGetReusePreconditioner - Determines if the `PC` reuses the current preconditioner even if the operator in the preconditioner has changed.
1352: Not Collective
1354: Input Parameter:
1355: . pc - the `PC` preconditioner context
1357: Output Parameter:
1358: . flag - `PETSC_TRUE` do not compute a new preconditioner, `PETSC_FALSE` do compute a new preconditioner
1360: Level: intermediate
1362: .seealso: [](ch_ksp), `PC`, `PCGetOperators()`, `MatZeroEntries()`, `PCSetReusePreconditioner()`
1363: @*/
1364: PetscErrorCode PCGetReusePreconditioner(PC pc, PetscBool *flag)
1365: {
1366: PetscFunctionBegin;
1368: PetscAssertPointer(flag, 2);
1369: *flag = pc->reusepreconditioner;
1370: PetscFunctionReturn(PETSC_SUCCESS);
1371: }
1373: /*@
1374: PCGetOperators - Gets the matrix associated with the linear system and
1375: possibly a different one which is used to construct the preconditioner.
1377: Not Collective, though parallel `Mat`s are returned if `pc` is parallel
1379: Input Parameter:
1380: . pc - the `PC` preconditioner context
1382: Output Parameters:
1383: + Amat - the matrix defining the linear system
1384: - Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.
1386: Level: intermediate
1388: Note:
1389: Does not increase the reference count of the matrices, so you should not destroy them
1391: Alternative usage: If the operators have NOT been set with `KSPSetOperators()`/`PCSetOperators()` then the operators
1392: are created in `PC` and returned to the user. In this case, if both operators
1393: mat and pmat are requested, two DIFFERENT operators will be returned. If
1394: only one is requested both operators in the PC will be the same (i.e. as
1395: if one had called `KSPSetOperators()`/`PCSetOperators()` with the same argument for both Mats).
1396: The user must set the sizes of the returned matrices and their type etc just
1397: as if the user created them with `MatCreate()`. For example,
1399: .vb
1400: KSP/PCGetOperators(ksp/pc,&Amat,NULL); is equivalent to
1401: set size, type, etc of Amat
1403: MatCreate(comm,&mat);
1404: KSP/PCSetOperators(ksp/pc,Amat,Amat);
1405: PetscObjectDereference((PetscObject)mat);
1406: set size, type, etc of Amat
1407: .ve
1409: and
1411: .vb
1412: KSP/PCGetOperators(ksp/pc,&Amat,&Pmat); is equivalent to
1413: set size, type, etc of Amat and Pmat
1415: MatCreate(comm,&Amat);
1416: MatCreate(comm,&Pmat);
1417: KSP/PCSetOperators(ksp/pc,Amat,Pmat);
1418: PetscObjectDereference((PetscObject)Amat);
1419: PetscObjectDereference((PetscObject)Pmat);
1420: set size, type, etc of Amat and Pmat
1421: .ve
1423: The rationale for this support is so that when creating a `TS`, `SNES`, or `KSP` the hierarchy
1424: of underlying objects (i.e. `SNES`, `KSP`, `PC`, `Mat`) and their lifespans can be completely
1425: managed by the top most level object (i.e. the `TS`, `SNES`, or `KSP`). Another way to look
1426: at this is when you create a `SNES` you do not NEED to create a `KSP` and attach it to
1427: the `SNES` object (the `SNES` object manages it for you). Similarly when you create a KSP
1428: you do not need to attach a `PC` to it (the `KSP` object manages the `PC` object for you).
1429: Thus, why should YOU have to create the `Mat` and attach it to the `SNES`/`KSP`/`PC`, when
1430: it can be created for you?
1432: .seealso: [](ch_ksp), `PC`, `PCSetOperators()`, `KSPGetOperators()`, `KSPSetOperators()`, `PCGetOperatorsSet()`
1433: @*/
1434: PetscErrorCode PCGetOperators(PC pc, Mat *Amat, Mat *Pmat)
1435: {
1436: PetscFunctionBegin;
1438: if (Amat) {
1439: if (!pc->mat) {
1440: if (pc->pmat && !Pmat) { /* Pmat has been set, but user did not request it, so use for Amat */
1441: pc->mat = pc->pmat;
1442: PetscCall(PetscObjectReference((PetscObject)pc->mat));
1443: } else { /* both Amat and Pmat are empty */
1444: PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &pc->mat));
1445: if (!Pmat) { /* user did NOT request Pmat, so make same as Amat */
1446: pc->pmat = pc->mat;
1447: PetscCall(PetscObjectReference((PetscObject)pc->pmat));
1448: }
1449: }
1450: }
1451: *Amat = pc->mat;
1452: }
1453: if (Pmat) {
1454: if (!pc->pmat) {
1455: if (pc->mat && !Amat) { /* Amat has been set but was not requested, so use for pmat */
1456: pc->pmat = pc->mat;
1457: PetscCall(PetscObjectReference((PetscObject)pc->pmat));
1458: } else {
1459: PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &pc->pmat));
1460: if (!Amat) { /* user did NOT request Amat, so make same as Pmat */
1461: pc->mat = pc->pmat;
1462: PetscCall(PetscObjectReference((PetscObject)pc->mat));
1463: }
1464: }
1465: }
1466: *Pmat = pc->pmat;
1467: }
1468: PetscFunctionReturn(PETSC_SUCCESS);
1469: }
1471: /*@
1472: PCGetOperatorsSet - Determines if the matrix associated with the linear system and
1473: possibly a different one associated with the preconditioner have been set in the `PC`.
1475: Not Collective, though the results on all processes should be the same
1477: Input Parameter:
1478: . pc - the `PC` preconditioner context
1480: Output Parameters:
1481: + mat - the matrix associated with the linear system was set
1482: - pmat - matrix associated with the preconditioner was set, usually the same
1484: Level: intermediate
1486: .seealso: [](ch_ksp), `PC`, `PCSetOperators()`, `KSPGetOperators()`, `KSPSetOperators()`, `PCGetOperators()`
1487: @*/
1488: PetscErrorCode PCGetOperatorsSet(PC pc, PetscBool *mat, PetscBool *pmat)
1489: {
1490: PetscFunctionBegin;
1492: if (mat) *mat = (pc->mat) ? PETSC_TRUE : PETSC_FALSE;
1493: if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1494: PetscFunctionReturn(PETSC_SUCCESS);
1495: }
1497: /*@
1498: PCFactorGetMatrix - Gets the factored matrix from the
1499: preconditioner context. This routine is valid only for the `PCLU`,
1500: `PCILU`, `PCCHOLESKY`, and `PCICC` methods.
1502: Not Collective though `mat` is parallel if `pc` is parallel
1504: Input Parameter:
1505: . pc - the `PC` preconditioner context
1507: Output Parameters:
1508: . mat - the factored matrix
1510: Level: advanced
1512: Note:
1513: Does not increase the reference count for `mat` so DO NOT destroy it
1515: .seealso: [](ch_ksp), `PC`, `PCLU`, `PCILU`, `PCCHOLESKY`, `PCICC`
1516: @*/
1517: PetscErrorCode PCFactorGetMatrix(PC pc, Mat *mat)
1518: {
1519: PetscFunctionBegin;
1521: PetscAssertPointer(mat, 2);
1522: PetscCall(PCFactorSetUpMatSolverType(pc));
1523: PetscUseTypeMethod(pc, getfactoredmatrix, mat);
1524: PetscFunctionReturn(PETSC_SUCCESS);
1525: }
1527: /*@
1528: PCSetOptionsPrefix - Sets the prefix used for searching for all
1529: `PC` options in the database.
1531: Logically Collective
1533: Input Parameters:
1534: + pc - the `PC` preconditioner context
1535: - prefix - the prefix string to prepend to all `PC` option requests
1537: Level: advanced
1539: Note:
1540: A hyphen (-) must NOT be given at the beginning of the prefix name.
1541: The first character of all runtime options is AUTOMATICALLY the
1542: hyphen.
1544: .seealso: [](ch_ksp), `PC`, `PCSetFromOptions()`, `PCAppendOptionsPrefix()`, `PCGetOptionsPrefix()`
1545: @*/
1546: PetscErrorCode PCSetOptionsPrefix(PC pc, const char prefix[])
1547: {
1548: PetscFunctionBegin;
1550: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc, prefix));
1551: PetscFunctionReturn(PETSC_SUCCESS);
1552: }
1554: /*@
1555: PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1556: `PC` options in the database.
1558: Logically Collective
1560: Input Parameters:
1561: + pc - the `PC` preconditioner context
1562: - prefix - the prefix string to prepend to all `PC` option requests
1564: Level: advanced
1566: Note:
1567: A hyphen (-) must NOT be given at the beginning of the prefix name.
1568: The first character of all runtime options is AUTOMATICALLY the
1569: hyphen.
1571: .seealso: [](ch_ksp), `PC`, `PCSetFromOptions()`, `PCSetOptionsPrefix()`, `PCGetOptionsPrefix()`
1572: @*/
1573: PetscErrorCode PCAppendOptionsPrefix(PC pc, const char prefix[])
1574: {
1575: PetscFunctionBegin;
1577: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)pc, prefix));
1578: PetscFunctionReturn(PETSC_SUCCESS);
1579: }
1581: /*@
1582: PCGetOptionsPrefix - Gets the prefix used for searching for all
1583: `PC` options in the database.
1585: Not Collective
1587: Input Parameter:
1588: . pc - the `PC` preconditioner context
1590: Output Parameter:
1591: . prefix - pointer to the prefix string used, is returned
1593: Level: advanced
1595: .seealso: [](ch_ksp), `PC`, `PCSetFromOptions()`, `PCSetOptionsPrefix()`, `PCAppendOptionsPrefix()`
1596: @*/
1597: PetscErrorCode PCGetOptionsPrefix(PC pc, const char *prefix[])
1598: {
1599: PetscFunctionBegin;
1601: PetscAssertPointer(prefix, 2);
1602: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, prefix));
1603: PetscFunctionReturn(PETSC_SUCCESS);
1604: }
1606: /*
1607: Indicates the right-hand side will be changed by KSPSolve(), this occurs for a few
1608: preconditioners including BDDC and Eisentat that transform the equations before applying
1609: the Krylov methods
1610: */
1611: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC pc, PetscBool *change)
1612: {
1613: PetscFunctionBegin;
1615: PetscAssertPointer(change, 2);
1616: *change = PETSC_FALSE;
1617: PetscTryMethod(pc, "PCPreSolveChangeRHS_C", (PC, PetscBool *), (pc, change));
1618: PetscFunctionReturn(PETSC_SUCCESS);
1619: }
1621: /*@
1622: PCPreSolve - Optional pre-solve phase, intended for any preconditioner-specific actions that must be performed before
1623: the iterative solve itself. Used in conjunction with `PCPostSolve()`
1625: Collective
1627: Input Parameters:
1628: + pc - the `PC` preconditioner context
1629: - ksp - the Krylov subspace context
1631: Level: developer
1633: Notes:
1634: `KSPSolve()` calls this directly, so is rarely called by the user.
1636: Certain preconditioners, such as the `PCType` of `PCEISENSTAT`, change the formulation of the linear system to be solved iteratively.
1637: This function performs that transformation. `PCPostSolve()` then transforms the system back to its original form after the solve.
1638: `PCPostSolve()` also transforms the resulting solution of the transformed system to the solution of the original problem.
1640: `KSPSetPostSolve()` provides an alternative way to provide such transformations.
1642: .seealso: [](ch_ksp), `PC`, `PCPostSolve()`, `KSP`, `PCSetPostSetUp()`, `KSPSetPreSolve()`, `KSPSetPostSolve()`
1643: @*/
1644: PetscErrorCode PCPreSolve(PC pc, KSP ksp)
1645: {
1646: Vec x, rhs;
1648: PetscFunctionBegin;
1651: pc->presolvedone++;
1652: PetscCheck(pc->presolvedone <= 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot embed PCPreSolve() more than twice");
1653: PetscCall(KSPGetSolution(ksp, &x));
1654: PetscCall(KSPGetRhs(ksp, &rhs));
1655: PetscTryTypeMethod(pc, presolve, ksp, rhs, x);
1656: PetscFunctionReturn(PETSC_SUCCESS);
1657: }
1659: /*@C
1660: PCSetPostSetUp - Sets function called at the end of `PCSetUp()` to adjust the computed preconditioner
1662: Logically Collective
1664: Input Parameters:
1665: + pc - the preconditioner object
1666: - postsetup - the function to call after `PCSetUp()`
1668: Calling sequence of `postsetup`:
1669: . pc - the `PC` context
1671: Level: developer
1673: .seealso: [](ch_ksp), `PC`, `PCSetUp()`
1674: @*/
1675: PetscErrorCode PCSetPostSetUp(PC pc, PetscErrorCode (*postsetup)(PC pc))
1676: {
1677: PetscFunctionBegin;
1679: pc->postsetup = postsetup;
1680: PetscFunctionReturn(PETSC_SUCCESS);
1681: }
1683: /*@
1684: PCPostSolve - Optional post-solve phase, intended for any
1685: preconditioner-specific actions that must be performed after
1686: the iterative solve itself.
1688: Collective
1690: Input Parameters:
1691: + pc - the `PC` preconditioner context
1692: - ksp - the `KSP` Krylov subspace context
1694: Example Usage:
1695: .vb
1696: PCPreSolve(pc,ksp);
1697: KSPSolve(ksp,b,x);
1698: PCPostSolve(pc,ksp);
1699: .ve
1701: Level: developer
1703: Note:
1704: `KSPSolve()` calls this routine directly, so it is rarely called by the user.
1706: .seealso: [](ch_ksp), `PC`, `KSPSetPostSolve()`, `KSPSetPreSolve()`, `PCPreSolve()`, `KSPSolve()`
1707: @*/
1708: PetscErrorCode PCPostSolve(PC pc, KSP ksp)
1709: {
1710: Vec x, rhs;
1712: PetscFunctionBegin;
1715: pc->presolvedone--;
1716: PetscCall(KSPGetSolution(ksp, &x));
1717: PetscCall(KSPGetRhs(ksp, &rhs));
1718: PetscTryTypeMethod(pc, postsolve, ksp, rhs, x);
1719: PetscFunctionReturn(PETSC_SUCCESS);
1720: }
1722: /*@
1723: PCLoad - Loads a `PC` that has been stored in binary with `PCView()`.
1725: Collective
1727: Input Parameters:
1728: + newdm - the newly loaded `PC`, this needs to have been created with `PCCreate()` or
1729: some related function before a call to `PCLoad()`.
1730: - viewer - binary file viewer `PETSCVIEWERBINARY`, obtained from `PetscViewerBinaryOpen()`
1732: Level: intermediate
1734: Note:
1735: The type is determined by the data in the file, any `PCType` set into the `PC` before this call is ignored.
1737: .seealso: [](ch_ksp), `PC`, `PetscViewerBinaryOpen()`, `PCView()`, `MatLoad()`, `VecLoad()`, `PETSCVIEWERBINARY`
1738: @*/
1739: PetscErrorCode PCLoad(PC newdm, PetscViewer viewer)
1740: {
1741: PetscBool isbinary;
1742: PetscInt classid;
1743: char type[256];
1745: PetscFunctionBegin;
1748: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1749: PetscCheck(isbinary, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1751: PetscCall(PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT));
1752: PetscCheck(classid == PC_FILE_CLASSID, PetscObjectComm((PetscObject)newdm), PETSC_ERR_ARG_WRONG, "Not PC next in file");
1753: PetscCall(PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR));
1754: PetscCall(PCSetType(newdm, type));
1755: PetscTryTypeMethod(newdm, load, viewer);
1756: PetscFunctionReturn(PETSC_SUCCESS);
1757: }
1759: #include <petscdraw.h>
1760: #if PetscDefined(HAVE_SAWS)
1761: #include <petscviewersaws.h>
1762: #endif
1764: /*@
1765: PCViewFromOptions - View (print or provide information about) the `PC`, based on options in the options database
1767: Collective
1769: Input Parameters:
1770: + A - the `PC` context
1771: . obj - Optional object that provides the options prefix
1772: - name - command line option name
1774: Options Database Key:
1775: . -name [viewertype][:...] - option name and values. See `PetscObjectViewFromOptions()` for the possible arguments
1777: Level: developer
1779: .seealso: [](ch_ksp), `PC`, `PCView`, `PetscObjectViewFromOptions()`, `PCCreate()`
1780: @*/
1781: PetscErrorCode PCViewFromOptions(PC A, PetscObject obj, const char name[])
1782: {
1783: PetscFunctionBegin;
1785: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
1786: PetscFunctionReturn(PETSC_SUCCESS);
1787: }
1789: /*@
1790: PCView - Prints information about the `PC`
1792: Collective
1794: Input Parameters:
1795: + pc - the `PC` preconditioner context
1796: - viewer - optional `PetscViewer` visualization context
1798: Level: intermediate
1800: Notes:
1801: The available visualization contexts include
1802: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
1803: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
1804: output where only the first processor opens
1805: the file. All other processors send their
1806: data to the first processor to print.
1808: The user can open an alternative visualization contexts with
1809: `PetscViewerASCIIOpen()` (output to a specified file).
1811: .seealso: [](ch_ksp), `PC`, `PetscViewer`, `PetscViewerType`, `KSPView()`, `PetscViewerASCIIOpen()`
1812: @*/
1813: PetscErrorCode PCView(PC pc, PetscViewer viewer)
1814: {
1815: PCType cstr;
1816: PetscViewerFormat format;
1817: PetscBool isascii, isstring, isbinary, isdraw, pop = PETSC_FALSE;
1818: #if PetscDefined(HAVE_SAWS)
1819: PetscBool issaws;
1820: #endif
1822: PetscFunctionBegin;
1824: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc), &viewer));
1826: PetscCheckSameComm(pc, 1, viewer, 2);
1828: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1829: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
1830: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1831: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1832: #if PetscDefined(HAVE_SAWS)
1833: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws));
1834: #endif
1836: if (isascii) {
1837: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)pc, viewer));
1838: if (!pc->setupcalled) PetscCall(PetscViewerASCIIPrintf(viewer, " PC has not been set up so information may be incomplete\n"));
1839: PetscCall(PetscViewerASCIIPushTab(viewer));
1840: PetscTryTypeMethod(pc, view, viewer);
1841: PetscCall(PetscViewerASCIIPopTab(viewer));
1842: if (pc->mat) {
1843: PetscCall(PetscViewerGetFormat(viewer, &format));
1844: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
1845: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
1846: pop = PETSC_TRUE;
1847: }
1848: if (pc->pmat == pc->mat) {
1849: PetscCall(PetscViewerASCIIPrintf(viewer, " linear system matrix, which is also used to construct the preconditioner:\n"));
1850: PetscCall(PetscViewerASCIIPushTab(viewer));
1851: PetscCall(MatView(pc->mat, viewer));
1852: PetscCall(PetscViewerASCIIPopTab(viewer));
1853: } else {
1854: if (pc->pmat) {
1855: PetscCall(PetscViewerASCIIPrintf(viewer, " linear system matrix, followed by the matrix used to construct the preconditioner:\n"));
1856: } else {
1857: PetscCall(PetscViewerASCIIPrintf(viewer, " linear system matrix:\n"));
1858: }
1859: PetscCall(PetscViewerASCIIPushTab(viewer));
1860: PetscCall(MatView(pc->mat, viewer));
1861: if (pc->pmat) PetscCall(MatView(pc->pmat, viewer));
1862: PetscCall(PetscViewerASCIIPopTab(viewer));
1863: }
1864: if (pop) PetscCall(PetscViewerPopFormat(viewer));
1865: }
1866: } else if (isstring) {
1867: PetscCall(PCGetType(pc, &cstr));
1868: PetscCall(PetscViewerStringSPrintf(viewer, " PCType: %-7.7s", cstr));
1869: PetscTryTypeMethod(pc, view, viewer);
1870: if (pc->mat) PetscCall(MatView(pc->mat, viewer));
1871: if (pc->pmat && pc->pmat != pc->mat) PetscCall(MatView(pc->pmat, viewer));
1872: } else if (isbinary) {
1873: PetscInt classid = PC_FILE_CLASSID;
1874: MPI_Comm comm;
1875: PetscMPIInt rank;
1876: char type[256];
1878: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
1879: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1880: if (rank == 0) {
1881: PetscCall(PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT));
1882: PetscCall(PetscStrncpy(type, ((PetscObject)pc)->type_name, 256));
1883: PetscCall(PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR));
1884: }
1885: PetscTryTypeMethod(pc, view, viewer);
1886: } else if (isdraw) {
1887: PetscDraw draw;
1888: char str[25];
1889: PetscReal x, y, bottom, h;
1890: PetscInt n;
1892: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1893: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
1894: if (pc->mat) {
1895: PetscCall(MatGetSize(pc->mat, &n, NULL));
1896: PetscCall(PetscSNPrintf(str, 25, "PC: %s (%" PetscInt_FMT ")", ((PetscObject)pc)->type_name, n));
1897: } else {
1898: PetscCall(PetscSNPrintf(str, 25, "PC: %s", ((PetscObject)pc)->type_name));
1899: }
1900: PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h));
1901: bottom = y - h;
1902: PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
1903: PetscTryTypeMethod(pc, view, viewer);
1904: PetscCall(PetscDrawPopCurrentPoint(draw));
1905: #if PetscDefined(HAVE_SAWS)
1906: } else if (issaws) {
1907: PetscMPIInt rank;
1909: PetscCall(PetscObjectName((PetscObject)pc));
1910: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1911: if (!((PetscObject)pc)->amsmem && rank == 0) PetscCall(PetscObjectViewSAWs((PetscObject)pc, viewer));
1912: if (pc->mat) PetscCall(MatView(pc->mat, viewer));
1913: if (pc->pmat && pc->pmat != pc->mat) PetscCall(MatView(pc->pmat, viewer));
1914: #endif
1915: }
1916: PetscFunctionReturn(PETSC_SUCCESS);
1917: }
1919: /*@C
1920: PCRegister - Adds a method (`PCType`) to the PETSc preconditioner package.
1922: Not collective. No Fortran Support
1924: Input Parameters:
1925: + sname - name of a new user-defined solver
1926: - function - routine to create the method context which will be stored in a `PC` when `PCSetType()` is called
1928: Example Usage:
1929: .vb
1930: PCRegister("my_solver", MySolverCreate);
1931: .ve
1933: Then, your solver can be chosen with the procedural interface via
1934: .vb
1935: PCSetType(pc, "my_solver")
1936: .ve
1937: or at runtime via the option
1938: .vb
1939: -pc_type my_solver
1940: .ve
1942: Level: advanced
1944: Note:
1945: A simpler alternative to using `PCRegister()` for an application specific preconditioner is to use a `PC` of `PCType` `PCSHELL` and
1946: provide your customizations with `PCShellSetContext()` and `PCShellSetApply()`
1948: `PCRegister()` may be called multiple times to add several user-defined preconditioners.
1950: .seealso: [](ch_ksp), `PC`, `PCType`, `PCRegisterAll()`, `PCSetType()`, `PCShellSetContext()`, `PCShellSetApply()`, `PCSHELL`
1951: @*/
1952: PetscErrorCode PCRegister(const char sname[], PetscErrorCode (*function)(PC))
1953: {
1954: PetscFunctionBegin;
1955: PetscCall(PCInitializePackage());
1956: PetscCall(PetscFunctionListAdd(&PCList, sname, function));
1957: PetscFunctionReturn(PETSC_SUCCESS);
1958: }
1960: static PetscErrorCode MatMult_PC(Mat A, Vec X, Vec Y)
1961: {
1962: PC pc;
1964: PetscFunctionBegin;
1965: PetscCall(MatShellGetContext(A, &pc));
1966: PetscCall(PCApply(pc, X, Y));
1967: PetscFunctionReturn(PETSC_SUCCESS);
1968: }
1970: /*@
1971: PCComputeOperator - Computes the explicit preconditioned operator as a matrix `Mat`.
1973: Collective
1975: Input Parameters:
1976: + pc - the `PC` preconditioner object
1977: - mattype - the `MatType` to be used for the operator
1979: Output Parameter:
1980: . mat - the explicit preconditioned operator
1982: Level: advanced
1984: Note:
1985: This computation is done by applying the operators to columns of the identity matrix.
1986: This routine is costly in general, and is recommended for use only with relatively small systems.
1987: Currently, this routine uses a dense matrix format when `mattype` == `NULL`
1989: Developer Note:
1990: This should be called `PCCreateExplicitOperator()`
1992: .seealso: [](ch_ksp), `PC`, `KSPComputeOperator()`, `MatType`
1993: @*/
1994: PetscErrorCode PCComputeOperator(PC pc, MatType mattype, Mat *mat)
1995: {
1996: PetscInt N, M, m, n;
1997: Mat A, Apc;
1999: PetscFunctionBegin;
2001: PetscAssertPointer(mat, 3);
2002: PetscCall(PCGetOperators(pc, &A, NULL));
2003: PetscCall(MatGetLocalSize(A, &m, &n));
2004: PetscCall(MatGetSize(A, &M, &N));
2005: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)pc), m, n, M, N, pc, &Apc));
2006: PetscCall(MatShellSetOperation(Apc, MATOP_MULT, (PetscErrorCodeFn *)MatMult_PC));
2007: PetscCall(MatComputeOperator(Apc, mattype, mat));
2008: PetscCall(MatDestroy(&Apc));
2009: PetscFunctionReturn(PETSC_SUCCESS);
2010: }
2012: /*@
2013: PCSetCoordinates - sets the coordinates of all the nodes (degrees of freedom in the vector) on the local process
2015: Collective
2017: Input Parameters:
2018: + pc - the `PC` preconditioner context
2019: . dim - the dimension of the coordinates 1, 2, or 3
2020: . nloc - the blocked size of the coordinates array
2021: - coords - the coordinates array
2023: Level: intermediate
2025: Notes:
2026: `coords` is an array of the dim coordinates for the nodes on
2027: the local processor, of size `dim`*`nloc`.
2028: If there are 108 equations (dofs) on a processor
2029: for a 3d displacement finite element discretization of elasticity (so
2030: that there are nloc = 36 = 108/3 nodes) then the array must have 108
2031: double precision values (ie, 3 * 36). These x y z coordinates
2032: should be ordered for nodes 0 to N-1 like so: [ 0.x, 0.y, 0.z, 1.x,
2033: ... , N-1.z ].
2035: The information provided here can be used by some preconditioners, such as `PCGAMG`, to produce a better preconditioner.
2036: See also `MatSetNearNullSpace()`.
2038: .seealso: [](ch_ksp), `PC`, `MatSetNearNullSpace()`
2039: @*/
2040: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[])
2041: {
2042: PetscFunctionBegin;
2045: PetscTryMethod(pc, "PCSetCoordinates_C", (PC, PetscInt, PetscInt, PetscReal[]), (pc, dim, nloc, coords));
2046: PetscFunctionReturn(PETSC_SUCCESS);
2047: }
2049: /*@
2050: PCGetInterpolations - Gets interpolation matrices for all levels (except level 0)
2052: Logically Collective
2054: Input Parameter:
2055: . pc - the precondition context
2057: Output Parameters:
2058: + num_levels - the number of levels
2059: - interpolations - the interpolation matrices (size of `num_levels`-1)
2061: Level: advanced
2063: Developer Note:
2064: Why is this here instead of in `PCMG` etc?
2066: .seealso: [](ch_ksp), `PC`, `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetInterpolation()`, `PCGetCoarseOperators()`
2067: @*/
2068: PetscErrorCode PCGetInterpolations(PC pc, PetscInt *num_levels, Mat *interpolations[])
2069: {
2070: PetscFunctionBegin;
2072: PetscAssertPointer(num_levels, 2);
2073: PetscAssertPointer(interpolations, 3);
2074: PetscUseMethod(pc, "PCGetInterpolations_C", (PC, PetscInt *, Mat *[]), (pc, num_levels, interpolations));
2075: PetscFunctionReturn(PETSC_SUCCESS);
2076: }
2078: /*@
2079: PCGetCoarseOperators - Gets coarse operator matrices for all levels (except the finest level)
2081: Logically Collective
2083: Input Parameter:
2084: . pc - the precondition context
2086: Output Parameters:
2087: + num_levels - the number of levels
2088: - coarseOperators - the coarse operator matrices (size of `num_levels`-1)
2090: Level: advanced
2092: Developer Note:
2093: Why is this here instead of in `PCMG` etc?
2095: .seealso: [](ch_ksp), `PC`, `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetRScale()`, `PCMGGetInterpolation()`, `PCGetInterpolations()`
2096: @*/
2097: PetscErrorCode PCGetCoarseOperators(PC pc, PetscInt *num_levels, Mat *coarseOperators[])
2098: {
2099: PetscFunctionBegin;
2101: PetscAssertPointer(num_levels, 2);
2102: PetscAssertPointer(coarseOperators, 3);
2103: PetscUseMethod(pc, "PCGetCoarseOperators_C", (PC, PetscInt *, Mat *[]), (pc, num_levels, coarseOperators));
2104: PetscFunctionReturn(PETSC_SUCCESS);
2105: }