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: *type = PCBJACOBI;
37: } else if (hasopsolve) {
38: *type = PCMAT;
39: } else {
40: *type = PCNONE;
41: }
42: } else {
43: if (hasopblock) {
44: *type = PCBJACOBI;
45: } else if (hasopsolve) {
46: *type = PCMAT;
47: } else {
48: *type = PCNONE;
49: }
50: }
51: } else *type = NULL;
52: PetscFunctionReturn(PETSC_SUCCESS);
53: }
55: /* do not log solves, setup and applications of preconditioners while constructing preconditioners; perhaps they should be logged separately from the regular solves */
56: PETSC_EXTERN PetscLogEvent KSP_Solve, KSP_SetUp;
58: static PetscErrorCode PCLogEventsDeactivatePush(void)
59: {
60: PetscFunctionBegin;
61: PetscCall(KSPInitializePackage());
62: PetscCall(PetscLogEventDeactivatePush(KSP_Solve));
63: PetscCall(PetscLogEventDeactivatePush(KSP_SetUp));
64: PetscCall(PetscLogEventDeactivatePush(PC_Apply));
65: PetscCall(PetscLogEventDeactivatePush(PC_SetUp));
66: PetscCall(PetscLogEventDeactivatePush(PC_SetUpOnBlocks));
67: PetscFunctionReturn(PETSC_SUCCESS);
68: }
70: static PetscErrorCode PCLogEventsDeactivatePop(void)
71: {
72: PetscFunctionBegin;
73: PetscCall(KSPInitializePackage());
74: PetscCall(PetscLogEventDeactivatePop(KSP_Solve));
75: PetscCall(PetscLogEventDeactivatePop(KSP_SetUp));
76: PetscCall(PetscLogEventDeactivatePop(PC_Apply));
77: PetscCall(PetscLogEventDeactivatePop(PC_SetUp));
78: PetscCall(PetscLogEventDeactivatePop(PC_SetUpOnBlocks));
79: PetscFunctionReturn(PETSC_SUCCESS);
80: }
82: /*@
83: 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
85: Collective
87: Input Parameter:
88: . pc - the `PC` preconditioner context
90: Level: developer
92: Notes:
93: Any options set, including those set with `KSPSetFromOptions()` remain.
95: 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`
97: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetUp()`
98: @*/
99: PetscErrorCode PCReset(PC pc)
100: {
101: PetscFunctionBegin;
103: PetscTryTypeMethod(pc, reset);
104: PetscCall(VecDestroy(&pc->diagonalscaleright));
105: PetscCall(VecDestroy(&pc->diagonalscaleleft));
106: PetscCall(MatDestroy(&pc->pmat));
107: PetscCall(MatDestroy(&pc->mat));
109: pc->setupcalled = PETSC_FALSE;
110: PetscFunctionReturn(PETSC_SUCCESS);
111: }
113: /*@
114: PCDestroy - Destroys `PC` context that was created with `PCCreate()`.
116: Collective
118: Input Parameter:
119: . pc - the `PC` preconditioner context
121: Level: developer
123: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetUp()`
124: @*/
125: PetscErrorCode PCDestroy(PC *pc)
126: {
127: PetscFunctionBegin;
128: if (!*pc) PetscFunctionReturn(PETSC_SUCCESS);
130: if (--((PetscObject)*pc)->refct > 0) {
131: *pc = NULL;
132: PetscFunctionReturn(PETSC_SUCCESS);
133: }
135: PetscCall(PCReset(*pc));
137: /* if memory was published with SAWs then destroy it */
138: PetscCall(PetscObjectSAWsViewOff((PetscObject)*pc));
139: PetscTryTypeMethod(*pc, destroy);
140: PetscCall(DMDestroy(&(*pc)->dm));
141: PetscCall(PetscHeaderDestroy(pc));
142: PetscFunctionReturn(PETSC_SUCCESS);
143: }
145: /*@
146: PCGetDiagonalScale - Indicates if the preconditioner applies an additional left and right
147: scaling as needed by certain time-stepping codes.
149: Logically Collective
151: Input Parameter:
152: . pc - the `PC` preconditioner context
154: Output Parameter:
155: . flag - `PETSC_TRUE` if it applies the scaling
157: Level: developer
159: Note:
160: If this returns `PETSC_TRUE` then the system solved via the Krylov method is, for left and right preconditioning,
162: $$
163: \begin{align*}
164: D M A D^{-1} y = D M b \\
165: D A M D^{-1} z = D b.
166: \end{align*}
167: $$
169: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCDiagonalScaleRight()`, `PCSetDiagonalScale()`
170: @*/
171: PetscErrorCode PCGetDiagonalScale(PC pc, PetscBool *flag)
172: {
173: PetscFunctionBegin;
175: PetscAssertPointer(flag, 2);
176: *flag = pc->diagonalscale;
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
180: /*@
181: PCSetDiagonalScale - Indicates the left scaling to use to apply an additional left and right
182: scaling as needed by certain time-stepping codes.
184: Logically Collective
186: Input Parameters:
187: + pc - the `PC` preconditioner context
188: - s - scaling vector
190: Level: intermediate
192: Notes:
193: The system solved via the Krylov method is, for left and right preconditioning,
194: $$
195: \begin{align*}
196: D M A D^{-1} y = D M b \\
197: D A M D^{-1} z = D b.
198: \end{align*}
199: $$
201: `PCDiagonalScaleLeft()` scales a vector by $D$. `PCDiagonalScaleRight()` scales a vector by $D^{-1}$.
203: .seealso: [](ch_ksp), `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCDiagonalScaleRight()`, `PCGetDiagonalScale()`
204: @*/
205: PetscErrorCode PCSetDiagonalScale(PC pc, Vec s)
206: {
207: PetscFunctionBegin;
210: pc->diagonalscale = PETSC_TRUE;
212: PetscCall(PetscObjectReference((PetscObject)s));
213: PetscCall(VecDestroy(&pc->diagonalscaleleft));
215: pc->diagonalscaleleft = s;
217: PetscCall(VecDuplicate(s, &pc->diagonalscaleright));
218: PetscCall(VecCopy(s, pc->diagonalscaleright));
219: PetscCall(VecReciprocal(pc->diagonalscaleright));
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
223: /*@
224: PCDiagonalScaleLeft - Scales a vector by the left scaling as needed by certain time-stepping codes.
226: Logically Collective
228: Input Parameters:
229: + pc - the `PC` preconditioner context
230: . in - input vector
231: - out - scaled vector (maybe the same as in)
233: Level: intermediate
235: Notes:
236: The system solved via the Krylov method is, for left and right preconditioning,
238: $$
239: \begin{align*}
240: D M A D^{-1} y = D M b \\
241: D A M D^{-1} z = D b.
242: \end{align*}
243: $$
245: `PCDiagonalScaleLeft()` scales a vector by $D$. `PCDiagonalScaleRight()` scales a vector by $D^{-1}$.
247: If diagonal scaling is turned off and `in` is not `out` then `in` is copied to `out`
249: .seealso: [](ch_ksp), `PCCreate()`, `PCSetUp()`, `PCSetDiagonalScale()`, `PCDiagonalScaleRight()`, `MatDiagonalScale()`
250: @*/
251: PetscErrorCode PCDiagonalScaleLeft(PC pc, Vec in, Vec out)
252: {
253: PetscFunctionBegin;
257: if (pc->diagonalscale) PetscCall(VecPointwiseMult(out, pc->diagonalscaleleft, in));
258: else if (in != out) PetscCall(VecCopy(in, out));
259: PetscFunctionReturn(PETSC_SUCCESS);
260: }
262: /*@
263: PCDiagonalScaleRight - Scales a vector by the right scaling as needed by certain time-stepping codes.
265: Logically Collective
267: Input Parameters:
268: + pc - the `PC` preconditioner context
269: . in - input vector
270: - out - scaled vector (maybe the same as in)
272: Level: intermediate
274: Notes:
275: The system solved via the Krylov method is, for left and right preconditioning,
277: $$
278: \begin{align*}
279: D M A D^{-1} y = D M b \\
280: D A M D^{-1} z = D b.
281: \end{align*}
282: $$
284: `PCDiagonalScaleLeft()` scales a vector by $D$. `PCDiagonalScaleRight()` scales a vector by $D^{-1}$.
286: If diagonal scaling is turned off and `in` is not `out` then `in` is copied to `out`
288: .seealso: [](ch_ksp), `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCSetDiagonalScale()`, `MatDiagonalScale()`
289: @*/
290: PetscErrorCode PCDiagonalScaleRight(PC pc, Vec in, Vec out)
291: {
292: PetscFunctionBegin;
296: if (pc->diagonalscale) {
297: PetscCall(VecPointwiseMult(out, pc->diagonalscaleright, in));
298: } else if (in != out) {
299: PetscCall(VecCopy(in, out));
300: }
301: PetscFunctionReturn(PETSC_SUCCESS);
302: }
304: /*@
305: PCSetUseAmat - Sets a flag to indicate that when the preconditioner needs to apply (part of) the
306: operator during the preconditioning process it applies the Amat provided to `TSSetRHSJacobian()`,
307: `TSSetIJacobian()`, `SNESSetJacobian()`, `KSPSetOperators()` or `PCSetOperators()` not the Pmat.
309: Logically Collective
311: Input Parameters:
312: + pc - the `PC` preconditioner context
313: - flg - `PETSC_TRUE` to use the Amat, `PETSC_FALSE` to use the Pmat (default is false)
315: Options Database Key:
316: . -pc_use_amat (true|false) - use the amat argument to `KSPSetOperators()` or `PCSetOperators()` to apply the operator
318: Level: intermediate
320: Note:
321: For the common case in which the linear system matrix and the matrix used to construct the
322: preconditioner are identical, this routine has no affect.
324: .seealso: [](ch_ksp), `PC`, `PCGetUseAmat()`, `PCBJACOBI`, `PCMG`, `PCFIELDSPLIT`, `PCCOMPOSITE`,
325: `KSPSetOperators()`, `PCSetOperators()`
326: @*/
327: PetscErrorCode PCSetUseAmat(PC pc, PetscBool flg)
328: {
329: PetscFunctionBegin;
331: pc->useAmat = flg;
332: PetscFunctionReturn(PETSC_SUCCESS);
333: }
335: /*@
336: PCSetErrorIfFailure - Causes `PC` to generate an error if a floating point exception, for example a zero pivot, is detected.
338: Logically Collective
340: Input Parameters:
341: + pc - iterative context obtained from `PCCreate()`
342: - flg - `PETSC_TRUE` indicates you want the error generated
344: Level: advanced
346: Notes:
347: Normally PETSc continues if a linear solver fails due to a failed setup of a preconditioner, you can call `KSPGetConvergedReason()` after a `KSPSolve()`
348: 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
349: detected.
351: This is propagated into `KSP`s used by this `PC`, which then propagate it into `PC`s used by those `KSP`s
353: .seealso: [](ch_ksp), `PC`, `KSPSetErrorIfNotConverged()`, `PCGetInitialGuessNonzero()`, `PCSetInitialGuessKnoll()`, `PCGetInitialGuessKnoll()`
354: @*/
355: PetscErrorCode PCSetErrorIfFailure(PC pc, PetscBool flg)
356: {
357: PetscFunctionBegin;
360: pc->erroriffailure = flg;
361: PetscFunctionReturn(PETSC_SUCCESS);
362: }
364: /*@
365: PCGetUseAmat - Gets a flag to indicate that when the preconditioner needs to apply (part of) the
366: operator during the preconditioning process it applies the Amat provided to `TSSetRHSJacobian()`,
367: `TSSetIJacobian()`, `SNESSetJacobian()`, `KSPSetOperators()` or `PCSetOperators()` not the Pmat.
369: Logically Collective
371: Input Parameter:
372: . pc - the `PC` preconditioner context
374: Output Parameter:
375: . flg - `PETSC_TRUE` to use the Amat, `PETSC_FALSE` to use the Pmat (default is false)
377: Level: intermediate
379: Note:
380: For the common case in which the linear system matrix and the matrix used to construct the
381: preconditioner are identical, this routine is does nothing.
383: .seealso: [](ch_ksp), `PC`, `PCSetUseAmat()`, `PCBJACOBI`, `PCMG`, `PCFIELDSPLIT`, `PCCOMPOSITE`
384: @*/
385: PetscErrorCode PCGetUseAmat(PC pc, PetscBool *flg)
386: {
387: PetscFunctionBegin;
389: *flg = pc->useAmat;
390: PetscFunctionReturn(PETSC_SUCCESS);
391: }
393: /*@
394: PCSetKSPNestLevel - sets the amount of nesting the `KSP` that contains this `PC` has
396: Collective
398: Input Parameters:
399: + pc - the `PC`
400: - level - the nest level
402: Level: developer
404: .seealso: [](ch_ksp), `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSP`, `KSPGMRES`, `KSPType`, `KSPGetNestLevel()`, `PCGetKSPNestLevel()`, `KSPSetNestLevel()`
405: @*/
406: PetscErrorCode PCSetKSPNestLevel(PC pc, PetscInt level)
407: {
408: PetscFunctionBegin;
411: pc->kspnestlevel = level;
412: PetscFunctionReturn(PETSC_SUCCESS);
413: }
415: /*@
416: PCGetKSPNestLevel - gets the amount of nesting the `KSP` that contains this `PC` has
418: Not Collective
420: Input Parameter:
421: . pc - the `PC`
423: Output Parameter:
424: . level - the nest level
426: Level: developer
428: .seealso: [](ch_ksp), `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSP`, `KSPGMRES`, `KSPType`, `KSPSetNestLevel()`, `PCSetKSPNestLevel()`, `KSPGetNestLevel()`
429: @*/
430: PetscErrorCode PCGetKSPNestLevel(PC pc, PetscInt *level)
431: {
432: PetscFunctionBegin;
434: PetscAssertPointer(level, 2);
435: *level = pc->kspnestlevel;
436: PetscFunctionReturn(PETSC_SUCCESS);
437: }
439: /*@
440: PCCreate - Creates a preconditioner context, `PC`
442: Collective
444: Input Parameter:
445: . comm - MPI communicator
447: Output Parameter:
448: . newpc - location to put the `PC` preconditioner context
450: Level: developer
452: Notes:
453: 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`.
455: Use `PCSetType()` or `PCSetFromOptions()` with the option `-pc_type pctype` to set the `PCType` for this `PC`
457: 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`
458: in parallel. For dense matrices it is always `PCNONE`.
460: .seealso: [](ch_ksp), `PC`, `PCType`, `PCSetType`, `PCSetUp()`, `PCApply()`, `PCDestroy()`, `KSP`, `KSPGetPC()`
461: @*/
462: PetscErrorCode PCCreate(MPI_Comm comm, PC *newpc)
463: {
464: PC pc;
466: PetscFunctionBegin;
467: PetscAssertPointer(newpc, 2);
468: PetscCall(PCInitializePackage());
470: PetscCall(PetscHeaderCreate(pc, PC_CLASSID, "PC", "Preconditioner", "PC", comm, PCDestroy, PCView));
471: pc->mat = NULL;
472: pc->pmat = NULL;
473: pc->setupcalled = PETSC_FALSE;
474: pc->setfromoptionscalled = 0;
475: pc->data = NULL;
476: pc->diagonalscale = PETSC_FALSE;
477: pc->diagonalscaleleft = NULL;
478: pc->diagonalscaleright = NULL;
480: pc->modifysubmatrices = NULL;
481: pc->modifysubmatricesP = NULL;
483: *newpc = pc;
484: PetscFunctionReturn(PETSC_SUCCESS);
485: }
487: /*@
488: PCApply - Applies the preconditioner to a vector.
490: Collective
492: Input Parameters:
493: + pc - the `PC` preconditioner context
494: - x - input vector
496: Output Parameter:
497: . y - output vector
499: Level: developer
501: .seealso: [](ch_ksp), `PC`, `PCApplyTranspose()`, `PCApplyBAorAB()`
502: @*/
503: PetscErrorCode PCApply(PC pc, Vec x, Vec y)
504: {
505: PetscInt m, n, mv, nv;
507: PetscFunctionBegin;
511: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
512: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
513: /* use pmat to check vector sizes since for KSPLSQR the pmat may be of a different size than mat */
514: PetscCall(MatGetLocalSize(pc->pmat, &m, &n));
515: PetscCall(VecGetLocalSize(x, &mv));
516: PetscCall(VecGetLocalSize(y, &nv));
517: /* check pmat * y = x is feasible */
518: 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);
519: 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);
520: PetscCall(VecSetErrorIfLocked(y, 3));
522: PetscCall(PCSetUp(pc));
523: PetscCall(VecLockReadPush(x));
524: PetscCall(PetscLogEventBegin(PC_Apply, pc, x, y, 0));
525: PetscUseTypeMethod(pc, apply, x, y);
526: PetscCall(PetscLogEventEnd(PC_Apply, pc, x, y, 0));
527: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
528: PetscCall(VecLockReadPop(x));
529: PetscFunctionReturn(PETSC_SUCCESS);
530: }
532: static PetscErrorCode PCMatApplyTranspose_Private(PC pc, Mat X, Mat Y, PetscBool transpose)
533: {
534: Mat A;
535: Vec cy, cx;
536: PetscInt m1, M1, m2, M2, n1, N1, n2, N2, m3, M3, n3, N3;
537: PetscBool match;
539: PetscFunctionBegin;
543: PetscCheckSameComm(pc, 1, X, 2);
544: PetscCheckSameComm(pc, 1, Y, 3);
545: PetscCheck(Y != X, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "Y and X must be different matrices");
546: PetscCall(PCGetOperators(pc, NULL, &A));
547: PetscCall(MatGetLocalSize(A, &m3, &n3));
548: PetscCall(MatGetLocalSize(X, &m2, &n2));
549: PetscCall(MatGetLocalSize(Y, &m1, &n1));
550: PetscCall(MatGetSize(A, &M3, &N3));
551: PetscCall(MatGetSize(X, &M2, &N2));
552: PetscCall(MatGetSize(Y, &M1, &N1));
553: 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);
554: 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);
555: 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);
556: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)Y, &match, MATSEQDENSE, MATMPIDENSE, ""));
557: PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided block of output vectors not stored in a dense Mat");
558: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)X, &match, MATSEQDENSE, MATMPIDENSE, ""));
559: PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided block of input vectors not stored in a dense Mat");
560: PetscCall(PCSetUp(pc));
561: if (!transpose && pc->ops->matapply) {
562: PetscCall(PetscLogEventBegin(PC_MatApply, pc, X, Y, 0));
563: PetscUseTypeMethod(pc, matapply, X, Y);
564: PetscCall(PetscLogEventEnd(PC_MatApply, pc, X, Y, 0));
565: } else if (transpose && pc->ops->matapplytranspose) {
566: PetscCall(PetscLogEventBegin(PC_MatApply, pc, X, Y, 0));
567: PetscUseTypeMethod(pc, matapplytranspose, X, Y);
568: PetscCall(PetscLogEventEnd(PC_MatApply, pc, X, Y, 0));
569: } else {
570: PetscCall(PetscInfo(pc, "PC type %s applying column by column\n", ((PetscObject)pc)->type_name));
571: for (n1 = 0; n1 < N1; ++n1) {
572: PetscCall(MatDenseGetColumnVecRead(X, n1, &cx));
573: PetscCall(MatDenseGetColumnVecWrite(Y, n1, &cy));
574: if (!transpose) PetscCall(PCApply(pc, cx, cy));
575: else PetscCall(PCApplyTranspose(pc, cx, cy));
576: PetscCall(MatDenseRestoreColumnVecWrite(Y, n1, &cy));
577: PetscCall(MatDenseRestoreColumnVecRead(X, n1, &cx));
578: }
579: }
580: PetscFunctionReturn(PETSC_SUCCESS);
581: }
583: /*@
584: PCMatApply - Applies the preconditioner to multiple vectors stored as a `MATDENSE`. Like `PCApply()`, `Y` and `X` must be different matrices.
586: Collective
588: Input Parameters:
589: + pc - the `PC` preconditioner context
590: - X - block of input vectors
592: Output Parameter:
593: . Y - block of output vectors
595: Level: developer
597: .seealso: [](ch_ksp), `PC`, `PCApply()`, `KSPMatSolve()`
598: @*/
599: PetscErrorCode PCMatApply(PC pc, Mat X, Mat Y)
600: {
601: PetscFunctionBegin;
602: PetscCall(PCMatApplyTranspose_Private(pc, X, Y, PETSC_FALSE));
603: PetscFunctionReturn(PETSC_SUCCESS);
604: }
606: /*@
607: PCMatApplyTranspose - Applies the transpose of preconditioner to multiple vectors stored as a `MATDENSE`. Like `PCApplyTranspose()`, `Y` and `X` must be different matrices.
609: Collective
611: Input Parameters:
612: + pc - the `PC` preconditioner context
613: - X - block of input vectors
615: Output Parameter:
616: . Y - block of output vectors
618: Level: developer
620: .seealso: [](ch_ksp), `PC`, `PCApplyTranspose()`, `KSPMatSolveTranspose()`
621: @*/
622: PetscErrorCode PCMatApplyTranspose(PC pc, Mat X, Mat Y)
623: {
624: PetscFunctionBegin;
625: PetscCall(PCMatApplyTranspose_Private(pc, X, Y, PETSC_TRUE));
626: PetscFunctionReturn(PETSC_SUCCESS);
627: }
629: /*@
630: PCApplySymmetricLeft - Applies the left part of a symmetric preconditioner to a vector.
632: Collective
634: Input Parameters:
635: + pc - the `PC` preconditioner context
636: - x - input vector
638: Output Parameter:
639: . y - output vector
641: Level: developer
643: Note:
644: Currently, this routine is implemented only for `PCICC` and `PCJACOBI` preconditioners.
646: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplySymmetricRight()`
647: @*/
648: PetscErrorCode PCApplySymmetricLeft(PC pc, Vec x, Vec y)
649: {
650: PetscFunctionBegin;
654: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
655: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
656: PetscCall(PCSetUp(pc));
657: PetscCall(VecLockReadPush(x));
658: PetscCall(PetscLogEventBegin(PC_ApplySymmetricLeft, pc, x, y, 0));
659: PetscUseTypeMethod(pc, applysymmetricleft, x, y);
660: PetscCall(PetscLogEventEnd(PC_ApplySymmetricLeft, pc, x, y, 0));
661: PetscCall(VecLockReadPop(x));
662: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
663: PetscFunctionReturn(PETSC_SUCCESS);
664: }
666: /*@
667: PCApplySymmetricRight - Applies the right part of a symmetric preconditioner to a vector.
669: Collective
671: Input Parameters:
672: + pc - the `PC` preconditioner context
673: - x - input vector
675: Output Parameter:
676: . y - output vector
678: Level: developer
680: Note:
681: Currently, this routine is implemented only for `PCICC` and `PCJACOBI` preconditioners.
683: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplySymmetricLeft()`
684: @*/
685: PetscErrorCode PCApplySymmetricRight(PC pc, Vec x, Vec y)
686: {
687: PetscFunctionBegin;
691: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
692: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
693: PetscCall(PCSetUp(pc));
694: PetscCall(VecLockReadPush(x));
695: PetscCall(PetscLogEventBegin(PC_ApplySymmetricRight, pc, x, y, 0));
696: PetscUseTypeMethod(pc, applysymmetricright, x, y);
697: PetscCall(PetscLogEventEnd(PC_ApplySymmetricRight, pc, x, y, 0));
698: PetscCall(VecLockReadPop(x));
699: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
700: PetscFunctionReturn(PETSC_SUCCESS);
701: }
703: /*@
704: PCApplyTranspose - Applies the transpose of preconditioner to a vector.
706: Collective
708: Input Parameters:
709: + pc - the `PC` preconditioner context
710: - x - input vector
712: Output Parameter:
713: . y - output vector
715: Level: developer
717: Note:
718: For complex numbers this applies the non-Hermitian transpose.
720: Developer Note:
721: We need to implement a `PCApplyHermitianTranspose()`
723: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplyBAorAB()`, `PCApplyBAorABTranspose()`, `PCApplyTransposeExists()`
724: @*/
725: PetscErrorCode PCApplyTranspose(PC pc, Vec x, Vec y)
726: {
727: PetscFunctionBegin;
731: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
732: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
733: PetscCall(PCSetUp(pc));
734: PetscCall(VecLockReadPush(x));
735: PetscCall(PetscLogEventBegin(PC_Apply, pc, x, y, 0));
736: PetscUseTypeMethod(pc, applytranspose, x, y);
737: PetscCall(PetscLogEventEnd(PC_Apply, pc, x, y, 0));
738: PetscCall(VecLockReadPop(x));
739: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
740: PetscFunctionReturn(PETSC_SUCCESS);
741: }
743: /*@
744: PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation
746: Collective
748: Input Parameter:
749: . pc - the `PC` preconditioner context
751: Output Parameter:
752: . flg - `PETSC_TRUE` if a transpose operation is defined
754: Level: developer
756: .seealso: [](ch_ksp), `PC`, `PCApplyTranspose()`
757: @*/
758: PetscErrorCode PCApplyTransposeExists(PC pc, PetscBool *flg)
759: {
760: PetscFunctionBegin;
762: PetscAssertPointer(flg, 2);
763: if (pc->ops->applytranspose) *flg = PETSC_TRUE;
764: else *flg = PETSC_FALSE;
765: PetscFunctionReturn(PETSC_SUCCESS);
766: }
768: /*@
769: PCApplyBAorAB - Applies the preconditioner and operator to a vector. $y = B*A*x $ or $ y = A*B*x$.
771: Collective
773: Input Parameters:
774: + pc - the `PC` preconditioner context
775: . side - indicates the preconditioner side, one of `PC_LEFT`, `PC_RIGHT`, or `PC_SYMMETRIC`
776: . x - input vector
777: - work - work vector
779: Output Parameter:
780: . y - output vector
782: Level: developer
784: Note:
785: 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.
786: The specific `KSPSolve()` method must also be written to handle the post-solve "correction" for the diagonal scaling.
788: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplyTranspose()`, `PCApplyBAorABTranspose()`
789: @*/
790: PetscErrorCode PCApplyBAorAB(PC pc, PCSide side, Vec x, Vec y, Vec work)
791: {
792: PetscFunctionBegin;
798: PetscCheckSameComm(pc, 1, x, 3);
799: PetscCheckSameComm(pc, 1, y, 4);
800: PetscCheckSameComm(pc, 1, work, 5);
801: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
802: PetscCheck(side == PC_LEFT || side == PC_SYMMETRIC || side == PC_RIGHT, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Side must be right, left, or symmetric");
803: PetscCheck(!pc->diagonalscale || side != PC_SYMMETRIC, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot include diagonal scaling with symmetric preconditioner application");
804: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 3, PETSC_TRUE));
806: PetscCall(PCSetUp(pc));
807: if (pc->diagonalscale) {
808: if (pc->ops->applyBA) {
809: Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
810: PetscCall(VecDuplicate(x, &work2));
811: PetscCall(PCDiagonalScaleRight(pc, x, work2));
812: PetscUseTypeMethod(pc, applyBA, side, work2, y, work);
813: PetscCall(PCDiagonalScaleLeft(pc, y, y));
814: PetscCall(VecDestroy(&work2));
815: } else if (side == PC_RIGHT) {
816: PetscCall(PCDiagonalScaleRight(pc, x, y));
817: PetscCall(PCApply(pc, y, work));
818: PetscCall(MatMult(pc->mat, work, y));
819: PetscCall(PCDiagonalScaleLeft(pc, y, y));
820: } else if (side == PC_LEFT) {
821: PetscCall(PCDiagonalScaleRight(pc, x, y));
822: PetscCall(MatMult(pc->mat, y, work));
823: PetscCall(PCApply(pc, work, y));
824: PetscCall(PCDiagonalScaleLeft(pc, y, y));
825: } else PetscCheck(side != PC_SYMMETRIC, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot provide diagonal scaling with symmetric application of preconditioner");
826: } else {
827: if (pc->ops->applyBA) {
828: PetscUseTypeMethod(pc, applyBA, side, x, y, work);
829: } else if (side == PC_RIGHT) {
830: PetscCall(PCApply(pc, x, work));
831: PetscCall(MatMult(pc->mat, work, y));
832: } else if (side == PC_LEFT) {
833: PetscCall(MatMult(pc->mat, x, work));
834: PetscCall(PCApply(pc, work, y));
835: } else if (side == PC_SYMMETRIC) {
836: /* There's an extra copy here; maybe should provide 2 work vectors instead? */
837: PetscCall(PCApplySymmetricRight(pc, x, work));
838: PetscCall(MatMult(pc->mat, work, y));
839: PetscCall(VecCopy(y, work));
840: PetscCall(PCApplySymmetricLeft(pc, work, y));
841: }
842: }
843: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 4, PETSC_FALSE));
844: PetscFunctionReturn(PETSC_SUCCESS);
845: }
847: /*@
848: PCApplyBAorABTranspose - Applies the transpose of the preconditioner
849: and operator to a vector. That is, applies $B^T * A^T$ with left preconditioning,
850: NOT $(B*A)^T = A^T*B^T$.
852: Collective
854: Input Parameters:
855: + pc - the `PC` preconditioner context
856: . side - indicates the preconditioner side, one of `PC_LEFT`, `PC_RIGHT`, or `PC_SYMMETRIC`
857: . x - input vector
858: - work - work vector
860: Output Parameter:
861: . y - output vector
863: Level: developer
865: Note:
866: 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
867: defined by $B^T$. This is why this has the funny form that it computes $B^T * A^T$
869: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplyTranspose()`, `PCApplyBAorAB()`
870: @*/
871: PetscErrorCode PCApplyBAorABTranspose(PC pc, PCSide side, Vec x, Vec y, Vec work)
872: {
873: PetscFunctionBegin;
878: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
879: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 3, PETSC_TRUE));
880: if (pc->ops->applyBAtranspose) {
881: PetscUseTypeMethod(pc, applyBAtranspose, side, x, y, work);
882: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 4, PETSC_FALSE));
883: PetscFunctionReturn(PETSC_SUCCESS);
884: }
885: PetscCheck(side == PC_LEFT || side == PC_RIGHT, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Side must be right or left");
887: PetscCall(PCSetUp(pc));
888: if (side == PC_RIGHT) {
889: PetscCall(PCApplyTranspose(pc, x, work));
890: PetscCall(MatMultTranspose(pc->mat, work, y));
891: } else if (side == PC_LEFT) {
892: PetscCall(MatMultTranspose(pc->mat, x, work));
893: PetscCall(PCApplyTranspose(pc, work, y));
894: }
895: /* add support for PC_SYMMETRIC */
896: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 4, PETSC_FALSE));
897: PetscFunctionReturn(PETSC_SUCCESS);
898: }
900: /*@
901: PCApplyRichardsonExists - Determines whether a particular preconditioner has a
902: built-in fast application of Richardson's method.
904: Not Collective
906: Input Parameter:
907: . pc - the preconditioner
909: Output Parameter:
910: . exists - `PETSC_TRUE` or `PETSC_FALSE`
912: Level: developer
914: .seealso: [](ch_ksp), `PC`, `KSPRICHARDSON`, `PCApplyRichardson()`
915: @*/
916: PetscErrorCode PCApplyRichardsonExists(PC pc, PetscBool *exists)
917: {
918: PetscFunctionBegin;
920: PetscAssertPointer(exists, 2);
921: if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
922: else *exists = PETSC_FALSE;
923: PetscFunctionReturn(PETSC_SUCCESS);
924: }
926: /*@
927: PCApplyRichardson - Applies several steps of Richardson iteration with
928: the particular preconditioner. This routine is usually used by the
929: Krylov solvers and not the application code directly.
931: Collective
933: Input Parameters:
934: + pc - the `PC` preconditioner context
935: . b - the right-hand side
936: . w - one work vector
937: . rtol - relative decrease in residual norm convergence criteria
938: . abstol - absolute residual norm convergence criteria
939: . dtol - divergence residual norm increase criteria
940: . its - the number of iterations to apply.
941: - guesszero - if the input x contains nonzero initial guess
943: Output Parameters:
944: + outits - number of iterations actually used (for SOR this always equals its)
945: . reason - the reason the apply terminated
946: - y - the solution (also contains initial guess if guesszero is `PETSC_FALSE`
948: Level: developer
950: Notes:
951: Most preconditioners do not support this function. Use the command
952: `PCApplyRichardsonExists()` to determine if one does.
954: Except for the `PCMG` this routine ignores the convergence tolerances
955: and always runs for the number of iterations
957: .seealso: [](ch_ksp), `PC`, `PCApplyRichardsonExists()`
958: @*/
959: PetscErrorCode PCApplyRichardson(PC pc, Vec b, Vec y, Vec w, PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt its, PetscBool guesszero, PetscInt *outits, PCRichardsonConvergedReason *reason)
960: {
961: PetscFunctionBegin;
966: PetscCheck(b != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "b and y must be different vectors");
967: PetscCall(PCSetUp(pc));
968: PetscUseTypeMethod(pc, applyrichardson, b, y, w, rtol, abstol, dtol, its, guesszero, outits, reason);
969: PetscFunctionReturn(PETSC_SUCCESS);
970: }
972: /*@
973: PCSetFailedReason - Sets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail
975: Logically Collective
977: Input Parameters:
978: + pc - the `PC` preconditioner context
979: - reason - the reason it failed
981: Level: advanced
983: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCFailedReason`
984: @*/
985: PetscErrorCode PCSetFailedReason(PC pc, PCFailedReason reason)
986: {
987: PetscFunctionBegin;
989: pc->failedreason = reason;
990: PetscFunctionReturn(PETSC_SUCCESS);
991: }
993: /*@
994: PCGetFailedReason - Gets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail
996: Not Collective
998: Input Parameter:
999: . pc - the `PC` preconditioner context
1001: Output Parameter:
1002: . reason - the reason it failed
1004: Level: advanced
1006: Note:
1007: After a call to `KSPCheckDot()` or `KSPCheckNorm()` inside a `KSPSolve()` or a call to `PCReduceFailedReason()`
1008: this is the maximum reason over all MPI processes in the `PC` communicator and hence logically collective.
1009: Otherwise it returns the local value.
1011: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCSetFailedReason()`, `PCFailedReason`
1012: @*/
1013: PetscErrorCode PCGetFailedReason(PC pc, PCFailedReason *reason)
1014: {
1015: PetscFunctionBegin;
1017: *reason = pc->failedreason;
1018: PetscFunctionReturn(PETSC_SUCCESS);
1019: }
1021: /*@
1022: PCReduceFailedReason - Reduce the failed reason among the MPI processes that share the `PC`
1024: Collective
1026: Input Parameter:
1027: . pc - the `PC` preconditioner context
1029: Level: advanced
1031: Note:
1032: Different MPI processes may have different reasons or no reason, see `PCGetFailedReason()`. This routine
1033: makes them have a common value (failure if any MPI process had a failure).
1035: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCGetFailedReason()`, `PCSetFailedReason()`, `PCFailedReason`
1036: @*/
1037: PetscErrorCode PCReduceFailedReason(PC pc)
1038: {
1039: PetscInt buf;
1041: PetscFunctionBegin;
1043: buf = (PetscInt)pc->failedreason;
1044: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &buf, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)pc)));
1045: pc->failedreason = (PCFailedReason)buf;
1046: PetscFunctionReturn(PETSC_SUCCESS);
1047: }
1049: /*
1050: a setupcall of 0 indicates never setup,
1051: 1 indicates has been previously setup
1052: -1 indicates a PCSetUp() was attempted and failed
1053: */
1054: /*@
1055: PCSetUp - Prepares for the use of a preconditioner. Performs all the one-time operations needed before the preconditioner
1056: can be used with `PCApply()`
1058: Collective
1060: Input Parameter:
1061: . pc - the `PC` preconditioner context
1063: Level: developer
1065: Notes:
1066: For example, for `PCLU` this will compute the factorization.
1068: This is called automatically by `KSPSetUp()` or `PCApply()` so rarely needs to be called directly.
1070: For nested preconditioners, such as `PCFIELDSPLIT` or `PCBJACOBI` this may not finish the construction of the preconditioner
1071: on the inner levels, the routine `PCSetUpOnBlocks()` may compute more of the preconditioner in those situations.
1073: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `KSPSetUp()`, `PCSetUpOnBlocks()`
1074: @*/
1075: PetscErrorCode PCSetUp(PC pc)
1076: {
1077: const char *def;
1078: PetscObjectState matstate, matnonzerostate;
1080: PetscFunctionBegin;
1082: PetscCheck(pc->mat, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Matrix must be set first");
1084: if (pc->setupcalled && pc->reusepreconditioner) {
1085: PetscCall(PetscInfo(pc, "Leaving PC with identical preconditioner since reuse preconditioner is set\n"));
1086: PetscFunctionReturn(PETSC_SUCCESS);
1087: }
1089: PetscCall(PetscObjectStateGet((PetscObject)pc->pmat, &matstate));
1090: PetscCall(MatGetNonzeroState(pc->pmat, &matnonzerostate));
1091: if (!pc->setupcalled) {
1092: //PetscCall(PetscInfo(pc, "Setting up PC for first time\n"));
1093: pc->flag = DIFFERENT_NONZERO_PATTERN;
1094: } else if (matstate == pc->matstate) PetscFunctionReturn(PETSC_SUCCESS);
1095: else {
1096: if (matnonzerostate != pc->matnonzerostate) {
1097: PetscCall(PetscInfo(pc, "Setting up PC with different nonzero pattern\n"));
1098: pc->flag = DIFFERENT_NONZERO_PATTERN;
1099: } else {
1100: //PetscCall(PetscInfo(pc, "Setting up PC with same nonzero pattern\n"));
1101: pc->flag = SAME_NONZERO_PATTERN;
1102: }
1103: }
1104: pc->matstate = matstate;
1105: pc->matnonzerostate = matnonzerostate;
1107: if (!((PetscObject)pc)->type_name) {
1108: PetscCall(PCGetDefaultType_Private(pc, &def));
1109: PetscCall(PCSetType(pc, def));
1110: }
1112: PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
1113: PetscCall(MatSetErrorIfFailure(pc->mat, pc->erroriffailure));
1114: PetscCall(PetscLogEventBegin(PC_SetUp, pc, 0, 0, 0));
1115: if (pc->ops->setup) {
1116: PetscCall(PCLogEventsDeactivatePush());
1117: PetscUseTypeMethod(pc, setup);
1118: PetscCall(PCLogEventsDeactivatePop());
1119: }
1120: PetscCall(PetscLogEventEnd(PC_SetUp, pc, 0, 0, 0));
1121: if (pc->postsetup) PetscCall((*pc->postsetup)(pc));
1122: if (!pc->setupcalled) pc->setupcalled = PETSC_TRUE;
1123: PetscFunctionReturn(PETSC_SUCCESS);
1124: }
1126: /*@
1127: PCSetUpOnBlocks - Sets up the preconditioner for each block in
1128: the block Jacobi, overlapping Schwarz, and fieldsplit methods.
1130: Collective
1132: Input Parameter:
1133: . pc - the `PC` preconditioner context
1135: Level: developer
1137: Notes:
1138: For nested preconditioners such as `PCBJACOBI`, `PCSetUp()` is not called on each sub-`KSP` when `PCSetUp()` is
1139: called on the outer `PC`, this routine ensures it is called.
1141: It calls `PCSetUp()` if not yet called.
1143: .seealso: [](ch_ksp), `PC`, `PCSetUp()`, `PCCreate()`, `PCApply()`, `PCDestroy()`
1144: @*/
1145: PetscErrorCode PCSetUpOnBlocks(PC pc)
1146: {
1147: PetscFunctionBegin;
1149: if (!pc->setupcalled) PetscCall(PCSetUp(pc)); /* "if" to prevent -info extra prints */
1150: if (!pc->ops->setuponblocks) PetscFunctionReturn(PETSC_SUCCESS);
1151: PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
1152: PetscCall(PetscLogEventBegin(PC_SetUpOnBlocks, pc, 0, 0, 0));
1153: PetscCall(PCLogEventsDeactivatePush());
1154: PetscUseTypeMethod(pc, setuponblocks);
1155: PetscCall(PCLogEventsDeactivatePop());
1156: PetscCall(PetscLogEventEnd(PC_SetUpOnBlocks, pc, 0, 0, 0));
1157: PetscFunctionReturn(PETSC_SUCCESS);
1158: }
1160: /*@C
1161: PCSetModifySubMatrices - Sets a user-defined routine for modifying the
1162: submatrices that arise within certain subdomain-based preconditioners such as `PCASM`
1164: Logically Collective
1166: Input Parameters:
1167: + pc - the `PC` preconditioner context
1168: . func - routine for modifying the submatrices, see `PCModifySubMatricesFn`
1169: - ctx - optional user-defined context (may be `NULL`)
1171: Level: advanced
1173: Notes:
1174: The basic submatrices are extracted from the matrix used to construct the preconditioner as
1175: usual; the user can then alter these (for example, to set different boundary
1176: conditions for each submatrix) before they are used for the local solves.
1178: `PCSetModifySubMatrices()` MUST be called before `KSPSetUp()` and
1179: `KSPSolve()`.
1181: A routine set by `PCSetModifySubMatrices()` is currently called within
1182: `PCBJACOBI`, `PCASM`, `PCGASM`, and `PCHPDDM`.
1183: All other preconditioners ignore this routine.
1185: .seealso: [](ch_ksp), `PC`, `PCModifySubMatricesFn`, `PCBJACOBI`, `PCASM`, `PCModifySubMatrices()`
1186: @*/
1187: PetscErrorCode PCSetModifySubMatrices(PC pc, PCModifySubMatricesFn *func, PetscCtx ctx)
1188: {
1189: PetscFunctionBegin;
1191: pc->modifysubmatrices = func;
1192: pc->modifysubmatricesP = ctx;
1193: PetscFunctionReturn(PETSC_SUCCESS);
1194: }
1196: /*@C
1197: PCModifySubMatrices - Calls an optional user-defined routine within
1198: certain preconditioners if one has been set with `PCSetModifySubMatrices()`.
1200: Collective
1202: Input Parameters:
1203: + pc - the `PC` preconditioner context
1204: . nsub - the number of local submatrices
1205: . row - an array of index sets that contain the global row numbers
1206: that comprise each local submatrix
1207: . col - an array of index sets that contain the global column numbers
1208: that comprise each local submatrix
1209: . submat - array of local submatrices
1210: - ctx - optional user-defined context for private data for the
1211: user-defined routine (may be `NULL`)
1213: Output Parameter:
1214: . submat - array of local submatrices (the entries of which may
1215: have been modified)
1217: Level: developer
1219: Note:
1220: The user should NOT generally call this routine, as it will
1221: automatically be called within certain preconditioners.
1223: .seealso: [](ch_ksp), `PC`, `PCModifySubMatricesFn`, `PCSetModifySubMatrices()`
1224: @*/
1225: PetscErrorCode PCModifySubMatrices(PC pc, PetscInt nsub, const IS row[], const IS col[], Mat submat[], PetscCtx ctx)
1226: {
1227: PetscFunctionBegin;
1229: if (!pc->modifysubmatrices) PetscFunctionReturn(PETSC_SUCCESS);
1230: PetscCall(PetscLogEventBegin(PC_ModifySubMatrices, pc, 0, 0, 0));
1231: PetscCall((*pc->modifysubmatrices)(pc, nsub, row, col, submat, ctx));
1232: PetscCall(PetscLogEventEnd(PC_ModifySubMatrices, pc, 0, 0, 0));
1233: PetscFunctionReturn(PETSC_SUCCESS);
1234: }
1236: /*@
1237: PCSetOperators - Sets the matrix associated with the linear system and
1238: a (possibly) different one from which the preconditioner will be constructed.
1240: Logically Collective
1242: Input Parameters:
1243: + pc - the `PC` preconditioner context
1244: . Amat - the matrix that defines the linear system
1245: - Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
1247: Level: advanced
1249: Notes:
1250: Using this routine directly is rarely needed, the preferred, and equivalent, usage is `KSPSetOperators()`.
1252: Passing a `NULL` for `Amat` or `Pmat` removes the matrix that is currently used.
1254: If you wish to replace either `Amat` or `Pmat` but leave the other one untouched then
1255: first call `KSPGetOperators()` to get the one you wish to keep, call `PetscObjectReference()`
1256: on it and then pass it back in in your call to `KSPSetOperators()`.
1258: More Notes about Repeated Solution of Linear Systems:
1259: PETSc does NOT reset the matrix entries of either `Amat` or `Pmat`
1260: to zero after a linear solve; the user is completely responsible for
1261: matrix assembly. See the routine `MatZeroEntries()` if desiring to
1262: zero all elements of a matrix.
1264: .seealso: [](ch_ksp), `PC`, `PCGetOperators()`, `MatZeroEntries()`
1265: @*/
1266: PetscErrorCode PCSetOperators(PC pc, Mat Amat, Mat Pmat)
1267: {
1268: PetscInt m1, n1, m2, n2;
1270: PetscFunctionBegin;
1274: if (Amat) PetscCheckSameComm(pc, 1, Amat, 2);
1275: if (Pmat) PetscCheckSameComm(pc, 1, Pmat, 3);
1276: if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
1277: PetscCall(MatGetLocalSize(Amat, &m1, &n1));
1278: PetscCall(MatGetLocalSize(pc->mat, &m2, &n2));
1279: 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);
1280: PetscCall(MatGetLocalSize(Pmat, &m1, &n1));
1281: PetscCall(MatGetLocalSize(pc->pmat, &m2, &n2));
1282: 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);
1283: }
1285: if (Pmat != pc->pmat) {
1286: /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
1287: pc->matnonzerostate = -1;
1288: pc->matstate = -1;
1289: }
1291: /* reference first in case the matrices are the same */
1292: if (Amat) PetscCall(PetscObjectReference((PetscObject)Amat));
1293: PetscCall(MatDestroy(&pc->mat));
1294: if (Pmat) PetscCall(PetscObjectReference((PetscObject)Pmat));
1295: PetscCall(MatDestroy(&pc->pmat));
1296: pc->mat = Amat;
1297: pc->pmat = Pmat;
1298: PetscFunctionReturn(PETSC_SUCCESS);
1299: }
1301: /*@
1302: PCSetReusePreconditioner - reuse the current preconditioner even if the operator in the preconditioner `PC` has changed.
1304: Logically Collective
1306: Input Parameters:
1307: + pc - the `PC` preconditioner context
1308: - flag - `PETSC_TRUE` do not compute a new preconditioner, `PETSC_FALSE` do compute a new preconditioner
1310: Level: intermediate
1312: Note:
1313: Normally if a matrix inside a `PC` changes the `PC` automatically updates itself using information from the changed matrix. This option
1314: prevents this.
1316: .seealso: [](ch_ksp), `PC`, `PCGetOperators()`, `MatZeroEntries()`, `PCGetReusePreconditioner()`, `KSPSetReusePreconditioner()`
1317: @*/
1318: PetscErrorCode PCSetReusePreconditioner(PC pc, PetscBool flag)
1319: {
1320: PetscFunctionBegin;
1323: pc->reusepreconditioner = flag;
1324: PetscTryMethod(pc, "PCSetReusePreconditioner_C", (PC, PetscBool), (pc, flag));
1325: PetscFunctionReturn(PETSC_SUCCESS);
1326: }
1328: /*@
1329: PCGetReusePreconditioner - Determines if the `PC` reuses the current preconditioner even if the operator in the preconditioner has changed.
1331: Not Collective
1333: Input Parameter:
1334: . pc - the `PC` preconditioner context
1336: Output Parameter:
1337: . flag - `PETSC_TRUE` do not compute a new preconditioner, `PETSC_FALSE` do compute a new preconditioner
1339: Level: intermediate
1341: .seealso: [](ch_ksp), `PC`, `PCGetOperators()`, `MatZeroEntries()`, `PCSetReusePreconditioner()`
1342: @*/
1343: PetscErrorCode PCGetReusePreconditioner(PC pc, PetscBool *flag)
1344: {
1345: PetscFunctionBegin;
1347: PetscAssertPointer(flag, 2);
1348: *flag = pc->reusepreconditioner;
1349: PetscFunctionReturn(PETSC_SUCCESS);
1350: }
1352: /*@
1353: PCGetOperators - Gets the matrix associated with the linear system and
1354: possibly a different one which is used to construct the preconditioner.
1356: Not Collective, though parallel `Mat`s are returned if `pc` is parallel
1358: Input Parameter:
1359: . pc - the `PC` preconditioner context
1361: Output Parameters:
1362: + Amat - the matrix defining the linear system
1363: - Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.
1365: Level: intermediate
1367: Note:
1368: Does not increase the reference count of the matrices, so you should not destroy them
1370: Alternative usage: If the operators have NOT been set with `KSPSetOperators()`/`PCSetOperators()` then the operators
1371: are created in `PC` and returned to the user. In this case, if both operators
1372: mat and pmat are requested, two DIFFERENT operators will be returned. If
1373: only one is requested both operators in the PC will be the same (i.e. as
1374: if one had called `KSPSetOperators()`/`PCSetOperators()` with the same argument for both Mats).
1375: The user must set the sizes of the returned matrices and their type etc just
1376: as if the user created them with `MatCreate()`. For example,
1378: .vb
1379: KSP/PCGetOperators(ksp/pc,&Amat,NULL); is equivalent to
1380: set size, type, etc of Amat
1382: MatCreate(comm,&mat);
1383: KSP/PCSetOperators(ksp/pc,Amat,Amat);
1384: PetscObjectDereference((PetscObject)mat);
1385: set size, type, etc of Amat
1386: .ve
1388: and
1390: .vb
1391: KSP/PCGetOperators(ksp/pc,&Amat,&Pmat); is equivalent to
1392: set size, type, etc of Amat and Pmat
1394: MatCreate(comm,&Amat);
1395: MatCreate(comm,&Pmat);
1396: KSP/PCSetOperators(ksp/pc,Amat,Pmat);
1397: PetscObjectDereference((PetscObject)Amat);
1398: PetscObjectDereference((PetscObject)Pmat);
1399: set size, type, etc of Amat and Pmat
1400: .ve
1402: The rationale for this support is so that when creating a `TS`, `SNES`, or `KSP` the hierarchy
1403: of underlying objects (i.e. `SNES`, `KSP`, `PC`, `Mat`) and their lifespans can be completely
1404: managed by the top most level object (i.e. the `TS`, `SNES`, or `KSP`). Another way to look
1405: at this is when you create a `SNES` you do not NEED to create a `KSP` and attach it to
1406: the `SNES` object (the `SNES` object manages it for you). Similarly when you create a KSP
1407: you do not need to attach a `PC` to it (the `KSP` object manages the `PC` object for you).
1408: Thus, why should YOU have to create the `Mat` and attach it to the `SNES`/`KSP`/`PC`, when
1409: it can be created for you?
1411: .seealso: [](ch_ksp), `PC`, `PCSetOperators()`, `KSPGetOperators()`, `KSPSetOperators()`, `PCGetOperatorsSet()`
1412: @*/
1413: PetscErrorCode PCGetOperators(PC pc, Mat *Amat, Mat *Pmat)
1414: {
1415: PetscFunctionBegin;
1417: if (Amat) {
1418: if (!pc->mat) {
1419: if (pc->pmat && !Pmat) { /* Pmat has been set, but user did not request it, so use for Amat */
1420: pc->mat = pc->pmat;
1421: PetscCall(PetscObjectReference((PetscObject)pc->mat));
1422: } else { /* both Amat and Pmat are empty */
1423: PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &pc->mat));
1424: if (!Pmat) { /* user did NOT request Pmat, so make same as Amat */
1425: pc->pmat = pc->mat;
1426: PetscCall(PetscObjectReference((PetscObject)pc->pmat));
1427: }
1428: }
1429: }
1430: *Amat = pc->mat;
1431: }
1432: if (Pmat) {
1433: if (!pc->pmat) {
1434: if (pc->mat && !Amat) { /* Amat has been set but was not requested, so use for pmat */
1435: pc->pmat = pc->mat;
1436: PetscCall(PetscObjectReference((PetscObject)pc->pmat));
1437: } else {
1438: PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &pc->pmat));
1439: if (!Amat) { /* user did NOT request Amat, so make same as Pmat */
1440: pc->mat = pc->pmat;
1441: PetscCall(PetscObjectReference((PetscObject)pc->mat));
1442: }
1443: }
1444: }
1445: *Pmat = pc->pmat;
1446: }
1447: PetscFunctionReturn(PETSC_SUCCESS);
1448: }
1450: /*@
1451: PCGetOperatorsSet - Determines if the matrix associated with the linear system and
1452: possibly a different one associated with the preconditioner have been set in the `PC`.
1454: Not Collective, though the results on all processes should be the same
1456: Input Parameter:
1457: . pc - the `PC` preconditioner context
1459: Output Parameters:
1460: + mat - the matrix associated with the linear system was set
1461: - pmat - matrix associated with the preconditioner was set, usually the same
1463: Level: intermediate
1465: .seealso: [](ch_ksp), `PC`, `PCSetOperators()`, `KSPGetOperators()`, `KSPSetOperators()`, `PCGetOperators()`
1466: @*/
1467: PetscErrorCode PCGetOperatorsSet(PC pc, PetscBool *mat, PetscBool *pmat)
1468: {
1469: PetscFunctionBegin;
1471: if (mat) *mat = (pc->mat) ? PETSC_TRUE : PETSC_FALSE;
1472: if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1473: PetscFunctionReturn(PETSC_SUCCESS);
1474: }
1476: /*@
1477: PCFactorGetMatrix - Gets the factored matrix from the
1478: preconditioner context. This routine is valid only for the `PCLU`,
1479: `PCILU`, `PCCHOLESKY`, and `PCICC` methods.
1481: Not Collective though `mat` is parallel if `pc` is parallel
1483: Input Parameter:
1484: . pc - the `PC` preconditioner context
1486: Output Parameters:
1487: . mat - the factored matrix
1489: Level: advanced
1491: Note:
1492: Does not increase the reference count for `mat` so DO NOT destroy it
1494: .seealso: [](ch_ksp), `PC`, `PCLU`, `PCILU`, `PCCHOLESKY`, `PCICC`
1495: @*/
1496: PetscErrorCode PCFactorGetMatrix(PC pc, Mat *mat)
1497: {
1498: PetscFunctionBegin;
1500: PetscAssertPointer(mat, 2);
1501: PetscCall(PCFactorSetUpMatSolverType(pc));
1502: PetscUseTypeMethod(pc, getfactoredmatrix, mat);
1503: PetscFunctionReturn(PETSC_SUCCESS);
1504: }
1506: /*@
1507: PCSetOptionsPrefix - Sets the prefix used for searching for all
1508: `PC` options in the database.
1510: Logically Collective
1512: Input Parameters:
1513: + pc - the `PC` preconditioner context
1514: - prefix - the prefix string to prepend to all `PC` option requests
1516: Note:
1517: A hyphen (-) must NOT be given at the beginning of the prefix name.
1518: The first character of all runtime options is AUTOMATICALLY the
1519: hyphen.
1521: Level: advanced
1523: .seealso: [](ch_ksp), `PC`, `PCSetFromOptions`, `PCAppendOptionsPrefix()`, `PCGetOptionsPrefix()`
1524: @*/
1525: PetscErrorCode PCSetOptionsPrefix(PC pc, const char prefix[])
1526: {
1527: PetscFunctionBegin;
1529: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc, prefix));
1530: PetscFunctionReturn(PETSC_SUCCESS);
1531: }
1533: /*@
1534: PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1535: `PC` options in the database.
1537: Logically Collective
1539: Input Parameters:
1540: + pc - the `PC` preconditioner context
1541: - prefix - the prefix string to prepend to all `PC` option requests
1543: Note:
1544: A hyphen (-) must NOT be given at the beginning of the prefix name.
1545: The first character of all runtime options is AUTOMATICALLY the
1546: hyphen.
1548: Level: advanced
1550: .seealso: [](ch_ksp), `PC`, `PCSetFromOptions`, `PCSetOptionsPrefix()`, `PCGetOptionsPrefix()`
1551: @*/
1552: PetscErrorCode PCAppendOptionsPrefix(PC pc, const char prefix[])
1553: {
1554: PetscFunctionBegin;
1556: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)pc, prefix));
1557: PetscFunctionReturn(PETSC_SUCCESS);
1558: }
1560: /*@
1561: PCGetOptionsPrefix - Gets the prefix used for searching for all
1562: PC options in the database.
1564: Not Collective
1566: Input Parameter:
1567: . pc - the `PC` preconditioner context
1569: Output Parameter:
1570: . prefix - pointer to the prefix string used, is returned
1572: Level: advanced
1574: .seealso: [](ch_ksp), `PC`, `PCSetFromOptions`, `PCSetOptionsPrefix()`, `PCAppendOptionsPrefix()`
1575: @*/
1576: PetscErrorCode PCGetOptionsPrefix(PC pc, const char *prefix[])
1577: {
1578: PetscFunctionBegin;
1580: PetscAssertPointer(prefix, 2);
1581: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, prefix));
1582: PetscFunctionReturn(PETSC_SUCCESS);
1583: }
1585: /*
1586: Indicates the right-hand side will be changed by KSPSolve(), this occurs for a few
1587: preconditioners including BDDC and Eisentat that transform the equations before applying
1588: the Krylov methods
1589: */
1590: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC pc, PetscBool *change)
1591: {
1592: PetscFunctionBegin;
1594: PetscAssertPointer(change, 2);
1595: *change = PETSC_FALSE;
1596: PetscTryMethod(pc, "PCPreSolveChangeRHS_C", (PC, PetscBool *), (pc, change));
1597: PetscFunctionReturn(PETSC_SUCCESS);
1598: }
1600: /*@
1601: PCPreSolve - Optional pre-solve phase, intended for any preconditioner-specific actions that must be performed before
1602: the iterative solve itself. Used in conjunction with `PCPostSolve()`
1604: Collective
1606: Input Parameters:
1607: + pc - the `PC` preconditioner context
1608: - ksp - the Krylov subspace context
1610: Level: developer
1612: Notes:
1613: `KSPSolve()` calls this directly, so is rarely called by the user.
1615: Certain preconditioners, such as the `PCType` of `PCEISENSTAT`, change the formulation of the linear system to be solved iteratively.
1616: This function performs that transformation. `PCPostSolve()` then transforms the system back to its original form after the solve.
1617: `PCPostSolve()` also transforms the resulting solution of the transformed system to the solution of the original problem.
1619: `KSPSetPostSolve()` provides an alternative way to provide such transformations.
1621: .seealso: [](ch_ksp), `PC`, `PCPostSolve()`, `KSP`, `PCSetPostSetUp()`, `KSPSetPreSolve()`, `KSPSetPostSolve()`
1622: @*/
1623: PetscErrorCode PCPreSolve(PC pc, KSP ksp)
1624: {
1625: Vec x, rhs;
1627: PetscFunctionBegin;
1630: pc->presolvedone++;
1631: PetscCheck(pc->presolvedone <= 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot embed PCPreSolve() more than twice");
1632: PetscCall(KSPGetSolution(ksp, &x));
1633: PetscCall(KSPGetRhs(ksp, &rhs));
1634: PetscTryTypeMethod(pc, presolve, ksp, rhs, x);
1635: PetscFunctionReturn(PETSC_SUCCESS);
1636: }
1638: /*@C
1639: PCSetPostSetUp - Sets function called at the end of `PCSetUp()` to adjust the computed preconditioner
1641: Logically Collective
1643: Input Parameters:
1644: + pc - the preconditioner object
1645: - postsetup - the function to call after `PCSetUp()`
1647: Calling sequence of `postsetup`:
1648: . pc - the `PC` context
1650: Level: developer
1652: .seealso: [](ch_ksp), `PC`, `PCSetUp()`
1653: @*/
1654: PetscErrorCode PCSetPostSetUp(PC pc, PetscErrorCode (*postsetup)(PC pc))
1655: {
1656: PetscFunctionBegin;
1658: pc->postsetup = postsetup;
1659: PetscFunctionReturn(PETSC_SUCCESS);
1660: }
1662: /*@
1663: PCPostSolve - Optional post-solve phase, intended for any
1664: preconditioner-specific actions that must be performed after
1665: the iterative solve itself.
1667: Collective
1669: Input Parameters:
1670: + pc - the `PC` preconditioner context
1671: - ksp - the `KSP` Krylov subspace context
1673: Example Usage:
1674: .vb
1675: PCPreSolve(pc,ksp);
1676: KSPSolve(ksp,b,x);
1677: PCPostSolve(pc,ksp);
1678: .ve
1680: Level: developer
1682: Note:
1683: `KSPSolve()` calls this routine directly, so it is rarely called by the user.
1685: .seealso: [](ch_ksp), `PC`, `KSPSetPostSolve()`, `KSPSetPreSolve()`, `PCPreSolve()`, `KSPSolve()`
1686: @*/
1687: PetscErrorCode PCPostSolve(PC pc, KSP ksp)
1688: {
1689: Vec x, rhs;
1691: PetscFunctionBegin;
1694: pc->presolvedone--;
1695: PetscCall(KSPGetSolution(ksp, &x));
1696: PetscCall(KSPGetRhs(ksp, &rhs));
1697: PetscTryTypeMethod(pc, postsolve, ksp, rhs, x);
1698: PetscFunctionReturn(PETSC_SUCCESS);
1699: }
1701: /*@
1702: PCLoad - Loads a `PC` that has been stored in binary with `PCView()`.
1704: Collective
1706: Input Parameters:
1707: + newdm - the newly loaded `PC`, this needs to have been created with `PCCreate()` or
1708: some related function before a call to `PCLoad()`.
1709: - viewer - binary file viewer `PETSCVIEWERBINARY`, obtained from `PetscViewerBinaryOpen()`
1711: Level: intermediate
1713: Note:
1714: The type is determined by the data in the file, any `PCType` set into the `PC` before this call is ignored.
1716: .seealso: [](ch_ksp), `PC`, `PetscViewerBinaryOpen()`, `PCView()`, `MatLoad()`, `VecLoad()`, `PETSCVIEWERBINARY`
1717: @*/
1718: PetscErrorCode PCLoad(PC newdm, PetscViewer viewer)
1719: {
1720: PetscBool isbinary;
1721: PetscInt classid;
1722: char type[256];
1724: PetscFunctionBegin;
1727: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1728: PetscCheck(isbinary, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1730: PetscCall(PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT));
1731: PetscCheck(classid == PC_FILE_CLASSID, PetscObjectComm((PetscObject)newdm), PETSC_ERR_ARG_WRONG, "Not PC next in file");
1732: PetscCall(PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR));
1733: PetscCall(PCSetType(newdm, type));
1734: PetscTryTypeMethod(newdm, load, viewer);
1735: PetscFunctionReturn(PETSC_SUCCESS);
1736: }
1738: #include <petscdraw.h>
1739: #if defined(PETSC_HAVE_SAWS)
1740: #include <petscviewersaws.h>
1741: #endif
1743: /*@
1744: PCViewFromOptions - View (print or provide information about) the `PC`, based on options in the options database
1746: Collective
1748: Input Parameters:
1749: + A - the `PC` context
1750: . obj - Optional object that provides the options prefix
1751: - name - command line option name
1753: Level: developer
1755: .seealso: [](ch_ksp), `PC`, `PCView`, `PetscObjectViewFromOptions()`, `PCCreate()`
1756: @*/
1757: PetscErrorCode PCViewFromOptions(PC A, PetscObject obj, const char name[])
1758: {
1759: PetscFunctionBegin;
1761: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
1762: PetscFunctionReturn(PETSC_SUCCESS);
1763: }
1765: /*@
1766: PCView - Prints information about the `PC`
1768: Collective
1770: Input Parameters:
1771: + pc - the `PC` preconditioner context
1772: - viewer - optional `PetscViewer` visualization context
1774: Level: intermediate
1776: Notes:
1777: The available visualization contexts include
1778: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
1779: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
1780: output where only the first processor opens
1781: the file. All other processors send their
1782: data to the first processor to print.
1784: The user can open an alternative visualization contexts with
1785: `PetscViewerASCIIOpen()` (output to a specified file).
1787: .seealso: [](ch_ksp), `PC`, `PetscViewer`, `PetscViewerType`, `KSPView()`, `PetscViewerASCIIOpen()`
1788: @*/
1789: PetscErrorCode PCView(PC pc, PetscViewer viewer)
1790: {
1791: PCType cstr;
1792: PetscViewerFormat format;
1793: PetscBool isascii, isstring, isbinary, isdraw, pop = PETSC_FALSE;
1794: #if defined(PETSC_HAVE_SAWS)
1795: PetscBool issaws;
1796: #endif
1798: PetscFunctionBegin;
1800: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc), &viewer));
1802: PetscCheckSameComm(pc, 1, viewer, 2);
1804: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1805: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
1806: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1807: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1808: #if defined(PETSC_HAVE_SAWS)
1809: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws));
1810: #endif
1812: if (isascii) {
1813: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)pc, viewer));
1814: if (!pc->setupcalled) PetscCall(PetscViewerASCIIPrintf(viewer, " PC has not been set up so information may be incomplete\n"));
1815: PetscCall(PetscViewerASCIIPushTab(viewer));
1816: PetscTryTypeMethod(pc, view, viewer);
1817: PetscCall(PetscViewerASCIIPopTab(viewer));
1818: if (pc->mat) {
1819: PetscCall(PetscViewerGetFormat(viewer, &format));
1820: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
1821: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
1822: pop = PETSC_TRUE;
1823: }
1824: if (pc->pmat == pc->mat) {
1825: PetscCall(PetscViewerASCIIPrintf(viewer, " linear system matrix, which is also used to construct the preconditioner:\n"));
1826: PetscCall(PetscViewerASCIIPushTab(viewer));
1827: PetscCall(MatView(pc->mat, viewer));
1828: PetscCall(PetscViewerASCIIPopTab(viewer));
1829: } else {
1830: if (pc->pmat) {
1831: PetscCall(PetscViewerASCIIPrintf(viewer, " linear system matrix, followed by the matrix used to construct the preconditioner:\n"));
1832: } else {
1833: PetscCall(PetscViewerASCIIPrintf(viewer, " linear system matrix:\n"));
1834: }
1835: PetscCall(PetscViewerASCIIPushTab(viewer));
1836: PetscCall(MatView(pc->mat, viewer));
1837: if (pc->pmat) PetscCall(MatView(pc->pmat, viewer));
1838: PetscCall(PetscViewerASCIIPopTab(viewer));
1839: }
1840: if (pop) PetscCall(PetscViewerPopFormat(viewer));
1841: }
1842: } else if (isstring) {
1843: PetscCall(PCGetType(pc, &cstr));
1844: PetscCall(PetscViewerStringSPrintf(viewer, " PCType: %-7.7s", cstr));
1845: PetscTryTypeMethod(pc, view, viewer);
1846: if (pc->mat) PetscCall(MatView(pc->mat, viewer));
1847: if (pc->pmat && pc->pmat != pc->mat) PetscCall(MatView(pc->pmat, viewer));
1848: } else if (isbinary) {
1849: PetscInt classid = PC_FILE_CLASSID;
1850: MPI_Comm comm;
1851: PetscMPIInt rank;
1852: char type[256];
1854: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
1855: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1856: if (rank == 0) {
1857: PetscCall(PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT));
1858: PetscCall(PetscStrncpy(type, ((PetscObject)pc)->type_name, 256));
1859: PetscCall(PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR));
1860: }
1861: PetscTryTypeMethod(pc, view, viewer);
1862: } else if (isdraw) {
1863: PetscDraw draw;
1864: char str[25];
1865: PetscReal x, y, bottom, h;
1866: PetscInt n;
1868: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1869: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
1870: if (pc->mat) {
1871: PetscCall(MatGetSize(pc->mat, &n, NULL));
1872: PetscCall(PetscSNPrintf(str, 25, "PC: %s (%" PetscInt_FMT ")", ((PetscObject)pc)->type_name, n));
1873: } else {
1874: PetscCall(PetscSNPrintf(str, 25, "PC: %s", ((PetscObject)pc)->type_name));
1875: }
1876: PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h));
1877: bottom = y - h;
1878: PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
1879: PetscTryTypeMethod(pc, view, viewer);
1880: PetscCall(PetscDrawPopCurrentPoint(draw));
1881: #if defined(PETSC_HAVE_SAWS)
1882: } else if (issaws) {
1883: PetscMPIInt rank;
1885: PetscCall(PetscObjectName((PetscObject)pc));
1886: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1887: if (!((PetscObject)pc)->amsmem && rank == 0) PetscCall(PetscObjectViewSAWs((PetscObject)pc, viewer));
1888: if (pc->mat) PetscCall(MatView(pc->mat, viewer));
1889: if (pc->pmat && pc->pmat != pc->mat) PetscCall(MatView(pc->pmat, viewer));
1890: #endif
1891: }
1892: PetscFunctionReturn(PETSC_SUCCESS);
1893: }
1895: /*@C
1896: PCRegister - Adds a method (`PCType`) to the PETSc preconditioner package.
1898: Not collective. No Fortran Support
1900: Input Parameters:
1901: + sname - name of a new user-defined solver
1902: - function - routine to create the method context which will be stored in a `PC` when `PCSetType()` is called
1904: Example Usage:
1905: .vb
1906: PCRegister("my_solver", MySolverCreate);
1907: .ve
1909: Then, your solver can be chosen with the procedural interface via
1910: .vb
1911: PCSetType(pc, "my_solver")
1912: .ve
1913: or at runtime via the option
1914: .vb
1915: -pc_type my_solver
1916: .ve
1918: Level: advanced
1920: Note:
1921: A simpler alternative to using `PCRegister()` for an application specific preconditioner is to use a `PC` of `PCType` `PCSHELL` and
1922: provide your customizations with `PCShellSetContext()` and `PCShellSetApply()`
1924: `PCRegister()` may be called multiple times to add several user-defined preconditioners.
1926: .seealso: [](ch_ksp), `PC`, `PCType`, `PCRegisterAll()`, `PCSetType()`, `PCShellSetContext()`, `PCShellSetApply()`, `PCSHELL`
1927: @*/
1928: PetscErrorCode PCRegister(const char sname[], PetscErrorCode (*function)(PC))
1929: {
1930: PetscFunctionBegin;
1931: PetscCall(PCInitializePackage());
1932: PetscCall(PetscFunctionListAdd(&PCList, sname, function));
1933: PetscFunctionReturn(PETSC_SUCCESS);
1934: }
1936: static PetscErrorCode MatMult_PC(Mat A, Vec X, Vec Y)
1937: {
1938: PC pc;
1940: PetscFunctionBegin;
1941: PetscCall(MatShellGetContext(A, &pc));
1942: PetscCall(PCApply(pc, X, Y));
1943: PetscFunctionReturn(PETSC_SUCCESS);
1944: }
1946: /*@
1947: PCComputeOperator - Computes the explicit preconditioned operator as a matrix `Mat`.
1949: Collective
1951: Input Parameters:
1952: + pc - the `PC` preconditioner object
1953: - mattype - the `MatType` to be used for the operator
1955: Output Parameter:
1956: . mat - the explicit preconditioned operator
1958: Level: advanced
1960: Note:
1961: This computation is done by applying the operators to columns of the identity matrix.
1962: This routine is costly in general, and is recommended for use only with relatively small systems.
1963: Currently, this routine uses a dense matrix format when `mattype` == `NULL`
1965: Developer Note:
1966: This should be called `PCCreateExplicitOperator()`
1968: .seealso: [](ch_ksp), `PC`, `KSPComputeOperator()`, `MatType`
1969: @*/
1970: PetscErrorCode PCComputeOperator(PC pc, MatType mattype, Mat *mat)
1971: {
1972: PetscInt N, M, m, n;
1973: Mat A, Apc;
1975: PetscFunctionBegin;
1977: PetscAssertPointer(mat, 3);
1978: PetscCall(PCGetOperators(pc, &A, NULL));
1979: PetscCall(MatGetLocalSize(A, &m, &n));
1980: PetscCall(MatGetSize(A, &M, &N));
1981: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)pc), m, n, M, N, pc, &Apc));
1982: PetscCall(MatShellSetOperation(Apc, MATOP_MULT, (PetscErrorCodeFn *)MatMult_PC));
1983: PetscCall(MatComputeOperator(Apc, mattype, mat));
1984: PetscCall(MatDestroy(&Apc));
1985: PetscFunctionReturn(PETSC_SUCCESS);
1986: }
1988: /*@
1989: PCSetCoordinates - sets the coordinates of all the nodes (degrees of freedom in the vector) on the local process
1991: Collective
1993: Input Parameters:
1994: + pc - the `PC` preconditioner context
1995: . dim - the dimension of the coordinates 1, 2, or 3
1996: . nloc - the blocked size of the coordinates array
1997: - coords - the coordinates array
1999: Level: intermediate
2001: Notes:
2002: `coords` is an array of the dim coordinates for the nodes on
2003: the local processor, of size `dim`*`nloc`.
2004: If there are 108 equations (dofs) on a processor
2005: for a 3d displacement finite element discretization of elasticity (so
2006: that there are nloc = 36 = 108/3 nodes) then the array must have 108
2007: double precision values (ie, 3 * 36). These x y z coordinates
2008: should be ordered for nodes 0 to N-1 like so: [ 0.x, 0.y, 0.z, 1.x,
2009: ... , N-1.z ].
2011: The information provided here can be used by some preconditioners, such as `PCGAMG`, to produce a better preconditioner.
2012: See also `MatSetNearNullSpace()`.
2014: .seealso: [](ch_ksp), `PC`, `MatSetNearNullSpace()`
2015: @*/
2016: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[])
2017: {
2018: PetscFunctionBegin;
2021: PetscTryMethod(pc, "PCSetCoordinates_C", (PC, PetscInt, PetscInt, PetscReal[]), (pc, dim, nloc, coords));
2022: PetscFunctionReturn(PETSC_SUCCESS);
2023: }
2025: /*@
2026: PCGetInterpolations - Gets interpolation matrices for all levels (except level 0)
2028: Logically Collective
2030: Input Parameter:
2031: . pc - the precondition context
2033: Output Parameters:
2034: + num_levels - the number of levels
2035: - interpolations - the interpolation matrices (size of `num_levels`-1)
2037: Level: advanced
2039: Developer Note:
2040: Why is this here instead of in `PCMG` etc?
2042: .seealso: [](ch_ksp), `PC`, `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetInterpolation()`, `PCGetCoarseOperators()`
2043: @*/
2044: PetscErrorCode PCGetInterpolations(PC pc, PetscInt *num_levels, Mat *interpolations[])
2045: {
2046: PetscFunctionBegin;
2048: PetscAssertPointer(num_levels, 2);
2049: PetscAssertPointer(interpolations, 3);
2050: PetscUseMethod(pc, "PCGetInterpolations_C", (PC, PetscInt *, Mat *[]), (pc, num_levels, interpolations));
2051: PetscFunctionReturn(PETSC_SUCCESS);
2052: }
2054: /*@
2055: PCGetCoarseOperators - Gets coarse operator matrices for all levels (except the finest level)
2057: Logically Collective
2059: Input Parameter:
2060: . pc - the precondition context
2062: Output Parameters:
2063: + num_levels - the number of levels
2064: - coarseOperators - the coarse operator matrices (size of `num_levels`-1)
2066: Level: advanced
2068: Developer Note:
2069: Why is this here instead of in `PCMG` etc?
2071: .seealso: [](ch_ksp), `PC`, `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetRScale()`, `PCMGGetInterpolation()`, `PCGetInterpolations()`
2072: @*/
2073: PetscErrorCode PCGetCoarseOperators(PC pc, PetscInt *num_levels, Mat *coarseOperators[])
2074: {
2075: PetscFunctionBegin;
2077: PetscAssertPointer(num_levels, 2);
2078: PetscAssertPointer(coarseOperators, 3);
2079: PetscUseMethod(pc, "PCGetCoarseOperators_C", (PC, PetscInt *, Mat *[]), (pc, num_levels, coarseOperators));
2080: PetscFunctionReturn(PETSC_SUCCESS);
2081: }