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 = 0;
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) {
258: PetscCall(VecPointwiseMult(out, pc->diagonalscaleleft, in));
259: } else if (in != out) {
260: PetscCall(VecCopy(in, out));
261: }
262: PetscFunctionReturn(PETSC_SUCCESS);
263: }
265: /*@
266: PCDiagonalScaleRight - Scales a vector by the right scaling as needed by certain time-stepping codes.
268: Logically Collective
270: Input Parameters:
271: + pc - the `PC` preconditioner context
272: . in - input vector
273: - out - scaled vector (maybe the same as in)
275: Level: intermediate
277: Notes:
278: The system solved via the Krylov method is, for left and right preconditioning,
280: $$
281: \begin{align*}
282: D M A D^{-1} y = D M b \\
283: D A M D^{-1} z = D b.
284: \end{align*}
285: $$
287: `PCDiagonalScaleLeft()` scales a vector by $D$. `PCDiagonalScaleRight()` scales a vector by $D^{-1}$.
289: If diagonal scaling is turned off and `in` is not `out` then `in` is copied to `out`
291: .seealso: [](ch_ksp), `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCSetDiagonalScale()`, `MatDiagonalScale()`
292: @*/
293: PetscErrorCode PCDiagonalScaleRight(PC pc, Vec in, Vec out)
294: {
295: PetscFunctionBegin;
299: if (pc->diagonalscale) {
300: PetscCall(VecPointwiseMult(out, pc->diagonalscaleright, in));
301: } else if (in != out) {
302: PetscCall(VecCopy(in, out));
303: }
304: PetscFunctionReturn(PETSC_SUCCESS);
305: }
307: /*@
308: PCSetUseAmat - Sets a flag to indicate that when the preconditioner needs to apply (part of) the
309: operator during the preconditioning process it applies the Amat provided to `TSSetRHSJacobian()`,
310: `TSSetIJacobian()`, `SNESSetJacobian()`, `KSPSetOperators()` or `PCSetOperators()` not the Pmat.
312: Logically Collective
314: Input Parameters:
315: + pc - the `PC` preconditioner context
316: - flg - `PETSC_TRUE` to use the Amat, `PETSC_FALSE` to use the Pmat (default is false)
318: Options Database Key:
319: . -pc_use_amat <true,false> - use the amat argument to `KSPSetOperators()` or `PCSetOperators()` to apply the operator
321: Level: intermediate
323: Note:
324: For the common case in which the linear system matrix and the matrix used to construct the
325: preconditioner are identical, this routine has no affect.
327: .seealso: [](ch_ksp), `PC`, `PCGetUseAmat()`, `PCBJACOBI`, `PCMG`, `PCFIELDSPLIT`, `PCCOMPOSITE`,
328: `KSPSetOperators()`, `PCSetOperators()`
329: @*/
330: PetscErrorCode PCSetUseAmat(PC pc, PetscBool flg)
331: {
332: PetscFunctionBegin;
334: pc->useAmat = flg;
335: PetscFunctionReturn(PETSC_SUCCESS);
336: }
338: /*@
339: PCSetErrorIfFailure - Causes `PC` to generate an error if a floating point exception, for example a zero pivot, is detected.
341: Logically Collective
343: Input Parameters:
344: + pc - iterative context obtained from `PCCreate()`
345: - flg - `PETSC_TRUE` indicates you want the error generated
347: Level: advanced
349: Notes:
350: Normally PETSc continues if a linear solver fails due to a failed setup of a preconditioner, you can call `KSPGetConvergedReason()` after a `KSPSolve()`
351: 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
352: detected.
354: This is propagated into `KSP`s used by this `PC`, which then propagate it into `PC`s used by those `KSP`s
356: .seealso: [](ch_ksp), `PC`, `KSPSetErrorIfNotConverged()`, `PCGetInitialGuessNonzero()`, `PCSetInitialGuessKnoll()`, `PCGetInitialGuessKnoll()`
357: @*/
358: PetscErrorCode PCSetErrorIfFailure(PC pc, PetscBool flg)
359: {
360: PetscFunctionBegin;
363: pc->erroriffailure = flg;
364: PetscFunctionReturn(PETSC_SUCCESS);
365: }
367: /*@
368: PCGetUseAmat - Gets a flag to indicate that when the preconditioner needs to apply (part of) the
369: operator during the preconditioning process it applies the Amat provided to `TSSetRHSJacobian()`,
370: `TSSetIJacobian()`, `SNESSetJacobian()`, `KSPSetOperators()` or `PCSetOperators()` not the Pmat.
372: Logically Collective
374: Input Parameter:
375: . pc - the `PC` preconditioner context
377: Output Parameter:
378: . flg - `PETSC_TRUE` to use the Amat, `PETSC_FALSE` to use the Pmat (default is false)
380: Level: intermediate
382: Note:
383: For the common case in which the linear system matrix and the matrix used to construct the
384: preconditioner are identical, this routine is does nothing.
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: pc->mat = NULL;
475: pc->pmat = NULL;
476: pc->setupcalled = 0;
477: pc->setfromoptionscalled = 0;
478: pc->data = NULL;
479: pc->diagonalscale = PETSC_FALSE;
480: pc->diagonalscaleleft = NULL;
481: pc->diagonalscaleright = NULL;
483: pc->modifysubmatrices = NULL;
484: pc->modifysubmatricesP = NULL;
486: *newpc = pc;
487: PetscFunctionReturn(PETSC_SUCCESS);
488: }
490: /*@
491: PCApply - Applies the preconditioner to a vector.
493: Collective
495: Input Parameters:
496: + pc - the `PC` preconditioner context
497: - x - input vector
499: Output Parameter:
500: . y - output vector
502: Level: developer
504: .seealso: [](ch_ksp), `PC`, `PCApplyTranspose()`, `PCApplyBAorAB()`
505: @*/
506: PetscErrorCode PCApply(PC pc, Vec x, Vec y)
507: {
508: PetscInt m, n, mv, nv;
510: PetscFunctionBegin;
514: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
515: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
516: /* use pmat to check vector sizes since for KSPLSQR the pmat may be of a different size than mat */
517: PetscCall(MatGetLocalSize(pc->pmat, &m, &n));
518: PetscCall(VecGetLocalSize(x, &mv));
519: PetscCall(VecGetLocalSize(y, &nv));
520: /* check pmat * y = x is feasible */
521: 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);
522: 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);
523: PetscCall(VecSetErrorIfLocked(y, 3));
525: PetscCall(PCSetUp(pc));
526: PetscCall(VecLockReadPush(x));
527: PetscCall(PetscLogEventBegin(PC_Apply, pc, x, y, 0));
528: PetscUseTypeMethod(pc, apply, x, y);
529: PetscCall(PetscLogEventEnd(PC_Apply, pc, x, y, 0));
530: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
531: PetscCall(VecLockReadPop(x));
532: PetscFunctionReturn(PETSC_SUCCESS);
533: }
535: static PetscErrorCode PCMatApplyTranspose_Private(PC pc, Mat X, Mat Y, PetscBool transpose)
536: {
537: Mat A;
538: Vec cy, cx;
539: PetscInt m1, M1, m2, M2, n1, N1, n2, N2, m3, M3, n3, N3;
540: PetscBool match;
542: PetscFunctionBegin;
546: PetscCheckSameComm(pc, 1, X, 2);
547: PetscCheckSameComm(pc, 1, Y, 3);
548: PetscCheck(Y != X, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "Y and X must be different matrices");
549: PetscCall(PCGetOperators(pc, NULL, &A));
550: PetscCall(MatGetLocalSize(A, &m3, &n3));
551: PetscCall(MatGetLocalSize(X, &m2, &n2));
552: PetscCall(MatGetLocalSize(Y, &m1, &n1));
553: PetscCall(MatGetSize(A, &M3, &N3));
554: PetscCall(MatGetSize(X, &M2, &N2));
555: PetscCall(MatGetSize(Y, &M1, &N1));
556: 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);
557: 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);
558: 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);
559: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)Y, &match, MATSEQDENSE, MATMPIDENSE, ""));
560: PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided block of output vectors not stored in a dense Mat");
561: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)X, &match, MATSEQDENSE, MATMPIDENSE, ""));
562: PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided block of input vectors not stored in a dense Mat");
563: PetscCall(PCSetUp(pc));
564: if (!transpose && pc->ops->matapply) {
565: PetscCall(PetscLogEventBegin(PC_MatApply, pc, X, Y, 0));
566: PetscUseTypeMethod(pc, matapply, X, Y);
567: PetscCall(PetscLogEventEnd(PC_MatApply, pc, X, Y, 0));
568: } else if (transpose && pc->ops->matapplytranspose) {
569: PetscCall(PetscLogEventBegin(PC_MatApply, pc, X, Y, 0));
570: PetscUseTypeMethod(pc, matapplytranspose, X, Y);
571: PetscCall(PetscLogEventEnd(PC_MatApply, pc, X, Y, 0));
572: } else {
573: PetscCall(PetscInfo(pc, "PC type %s applying column by column\n", ((PetscObject)pc)->type_name));
574: for (n1 = 0; n1 < N1; ++n1) {
575: PetscCall(MatDenseGetColumnVecRead(X, n1, &cx));
576: PetscCall(MatDenseGetColumnVecWrite(Y, n1, &cy));
577: if (!transpose) PetscCall(PCApply(pc, cx, cy));
578: else PetscCall(PCApplyTranspose(pc, cx, cy));
579: PetscCall(MatDenseRestoreColumnVecWrite(Y, n1, &cy));
580: PetscCall(MatDenseRestoreColumnVecRead(X, n1, &cx));
581: }
582: }
583: PetscFunctionReturn(PETSC_SUCCESS);
584: }
586: /*@
587: PCMatApply - Applies the preconditioner to multiple vectors stored as a `MATDENSE`. Like `PCApply()`, `Y` and `X` must be different matrices.
589: Collective
591: Input Parameters:
592: + pc - the `PC` preconditioner context
593: - X - block of input vectors
595: Output Parameter:
596: . Y - block of output vectors
598: Level: developer
600: .seealso: [](ch_ksp), `PC`, `PCApply()`, `KSPMatSolve()`
601: @*/
602: PetscErrorCode PCMatApply(PC pc, Mat X, Mat Y)
603: {
604: PetscFunctionBegin;
605: PetscCall(PCMatApplyTranspose_Private(pc, X, Y, PETSC_FALSE));
606: PetscFunctionReturn(PETSC_SUCCESS);
607: }
609: /*@
610: PCMatApplyTranspose - Applies the transpose of preconditioner to multiple vectors stored as a `MATDENSE`. Like `PCApplyTranspose()`, `Y` and `X` must be different matrices.
612: Collective
614: Input Parameters:
615: + pc - the `PC` preconditioner context
616: - X - block of input vectors
618: Output Parameter:
619: . Y - block of output vectors
621: Level: developer
623: .seealso: [](ch_ksp), `PC`, `PCApplyTranspose()`, `KSPMatSolveTranspose()`
624: @*/
625: PetscErrorCode PCMatApplyTranspose(PC pc, Mat X, Mat Y)
626: {
627: PetscFunctionBegin;
628: PetscCall(PCMatApplyTranspose_Private(pc, X, Y, PETSC_TRUE));
629: PetscFunctionReturn(PETSC_SUCCESS);
630: }
632: /*@
633: PCApplySymmetricLeft - Applies the left part of a symmetric preconditioner to a vector.
635: Collective
637: Input Parameters:
638: + pc - the `PC` preconditioner context
639: - x - input vector
641: Output Parameter:
642: . y - output vector
644: Level: developer
646: Note:
647: Currently, this routine is implemented only for `PCICC` and `PCJACOBI` preconditioners.
649: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplySymmetricRight()`
650: @*/
651: PetscErrorCode PCApplySymmetricLeft(PC pc, Vec x, Vec y)
652: {
653: PetscFunctionBegin;
657: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
658: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
659: PetscCall(PCSetUp(pc));
660: PetscCall(VecLockReadPush(x));
661: PetscCall(PetscLogEventBegin(PC_ApplySymmetricLeft, pc, x, y, 0));
662: PetscUseTypeMethod(pc, applysymmetricleft, x, y);
663: PetscCall(PetscLogEventEnd(PC_ApplySymmetricLeft, pc, x, y, 0));
664: PetscCall(VecLockReadPop(x));
665: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
666: PetscFunctionReturn(PETSC_SUCCESS);
667: }
669: /*@
670: PCApplySymmetricRight - Applies the right part of a symmetric preconditioner to a vector.
672: Collective
674: Input Parameters:
675: + pc - the `PC` preconditioner context
676: - x - input vector
678: Output Parameter:
679: . y - output vector
681: Level: developer
683: Note:
684: Currently, this routine is implemented only for `PCICC` and `PCJACOBI` preconditioners.
686: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplySymmetricLeft()`
687: @*/
688: PetscErrorCode PCApplySymmetricRight(PC pc, Vec x, Vec y)
689: {
690: PetscFunctionBegin;
694: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
695: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
696: PetscCall(PCSetUp(pc));
697: PetscCall(VecLockReadPush(x));
698: PetscCall(PetscLogEventBegin(PC_ApplySymmetricRight, pc, x, y, 0));
699: PetscUseTypeMethod(pc, applysymmetricright, x, y);
700: PetscCall(PetscLogEventEnd(PC_ApplySymmetricRight, pc, x, y, 0));
701: PetscCall(VecLockReadPop(x));
702: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
703: PetscFunctionReturn(PETSC_SUCCESS);
704: }
706: /*@
707: PCApplyTranspose - Applies the transpose of preconditioner to a vector.
709: Collective
711: Input Parameters:
712: + pc - the `PC` preconditioner context
713: - x - input vector
715: Output Parameter:
716: . y - output vector
718: Level: developer
720: Note:
721: For complex numbers this applies the non-Hermitian transpose.
723: Developer Note:
724: We need to implement a `PCApplyHermitianTranspose()`
726: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplyBAorAB()`, `PCApplyBAorABTranspose()`, `PCApplyTransposeExists()`
727: @*/
728: PetscErrorCode PCApplyTranspose(PC pc, Vec x, Vec y)
729: {
730: PetscFunctionBegin;
734: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
735: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
736: PetscCall(PCSetUp(pc));
737: PetscCall(VecLockReadPush(x));
738: PetscCall(PetscLogEventBegin(PC_Apply, pc, x, y, 0));
739: PetscUseTypeMethod(pc, applytranspose, x, y);
740: PetscCall(PetscLogEventEnd(PC_Apply, pc, x, y, 0));
741: PetscCall(VecLockReadPop(x));
742: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
743: PetscFunctionReturn(PETSC_SUCCESS);
744: }
746: /*@
747: PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation
749: Collective
751: Input Parameter:
752: . pc - the `PC` preconditioner context
754: Output Parameter:
755: . flg - `PETSC_TRUE` if a transpose operation is defined
757: Level: developer
759: .seealso: [](ch_ksp), `PC`, `PCApplyTranspose()`
760: @*/
761: PetscErrorCode PCApplyTransposeExists(PC pc, PetscBool *flg)
762: {
763: PetscFunctionBegin;
765: PetscAssertPointer(flg, 2);
766: if (pc->ops->applytranspose) *flg = PETSC_TRUE;
767: else *flg = PETSC_FALSE;
768: PetscFunctionReturn(PETSC_SUCCESS);
769: }
771: /*@
772: PCApplyBAorAB - Applies the preconditioner and operator to a vector. $y = B*A*x $ or $ y = A*B*x$.
774: Collective
776: Input Parameters:
777: + pc - the `PC` preconditioner context
778: . side - indicates the preconditioner side, one of `PC_LEFT`, `PC_RIGHT`, or `PC_SYMMETRIC`
779: . x - input vector
780: - work - work vector
782: Output Parameter:
783: . y - output vector
785: Level: developer
787: Note:
788: 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.
789: The specific `KSPSolve()` method must also be written to handle the post-solve "correction" for the diagonal scaling.
791: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplyTranspose()`, `PCApplyBAorABTranspose()`
792: @*/
793: PetscErrorCode PCApplyBAorAB(PC pc, PCSide side, Vec x, Vec y, Vec work)
794: {
795: PetscFunctionBegin;
801: PetscCheckSameComm(pc, 1, x, 3);
802: PetscCheckSameComm(pc, 1, y, 4);
803: PetscCheckSameComm(pc, 1, work, 5);
804: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
805: PetscCheck(side == PC_LEFT || side == PC_SYMMETRIC || side == PC_RIGHT, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Side must be right, left, or symmetric");
806: PetscCheck(!pc->diagonalscale || side != PC_SYMMETRIC, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot include diagonal scaling with symmetric preconditioner application");
807: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 3, PETSC_TRUE));
809: PetscCall(PCSetUp(pc));
810: if (pc->diagonalscale) {
811: if (pc->ops->applyBA) {
812: Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
813: PetscCall(VecDuplicate(x, &work2));
814: PetscCall(PCDiagonalScaleRight(pc, x, work2));
815: PetscUseTypeMethod(pc, applyBA, side, work2, y, work);
816: PetscCall(PCDiagonalScaleLeft(pc, y, y));
817: PetscCall(VecDestroy(&work2));
818: } else if (side == PC_RIGHT) {
819: PetscCall(PCDiagonalScaleRight(pc, x, y));
820: PetscCall(PCApply(pc, y, work));
821: PetscCall(MatMult(pc->mat, work, y));
822: PetscCall(PCDiagonalScaleLeft(pc, y, y));
823: } else if (side == PC_LEFT) {
824: PetscCall(PCDiagonalScaleRight(pc, x, y));
825: PetscCall(MatMult(pc->mat, y, work));
826: PetscCall(PCApply(pc, work, y));
827: PetscCall(PCDiagonalScaleLeft(pc, y, y));
828: } else PetscCheck(side != PC_SYMMETRIC, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot provide diagonal scaling with symmetric application of preconditioner");
829: } else {
830: if (pc->ops->applyBA) {
831: PetscUseTypeMethod(pc, applyBA, side, x, y, work);
832: } else if (side == PC_RIGHT) {
833: PetscCall(PCApply(pc, x, work));
834: PetscCall(MatMult(pc->mat, work, y));
835: } else if (side == PC_LEFT) {
836: PetscCall(MatMult(pc->mat, x, work));
837: PetscCall(PCApply(pc, work, y));
838: } else if (side == PC_SYMMETRIC) {
839: /* There's an extra copy here; maybe should provide 2 work vectors instead? */
840: PetscCall(PCApplySymmetricRight(pc, x, work));
841: PetscCall(MatMult(pc->mat, work, y));
842: PetscCall(VecCopy(y, work));
843: PetscCall(PCApplySymmetricLeft(pc, work, y));
844: }
845: }
846: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 4, PETSC_FALSE));
847: PetscFunctionReturn(PETSC_SUCCESS);
848: }
850: /*@
851: PCApplyBAorABTranspose - Applies the transpose of the preconditioner
852: and operator to a vector. That is, applies $B^T * A^T$ with left preconditioning,
853: NOT $(B*A)^T = A^T*B^T$.
855: Collective
857: Input Parameters:
858: + pc - the `PC` preconditioner context
859: . side - indicates the preconditioner side, one of `PC_LEFT`, `PC_RIGHT`, or `PC_SYMMETRIC`
860: . x - input vector
861: - work - work vector
863: Output Parameter:
864: . y - output vector
866: Level: developer
868: Note:
869: 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
870: defined by $B^T$. This is why this has the funny form that it computes $B^T * A^T$
872: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplyTranspose()`, `PCApplyBAorAB()`
873: @*/
874: PetscErrorCode PCApplyBAorABTranspose(PC pc, PCSide side, Vec x, Vec y, Vec work)
875: {
876: PetscFunctionBegin;
881: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
882: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 3, PETSC_TRUE));
883: if (pc->ops->applyBAtranspose) {
884: PetscUseTypeMethod(pc, applyBAtranspose, side, x, y, work);
885: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 4, PETSC_FALSE));
886: PetscFunctionReturn(PETSC_SUCCESS);
887: }
888: PetscCheck(side == PC_LEFT || side == PC_RIGHT, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Side must be right or left");
890: PetscCall(PCSetUp(pc));
891: if (side == PC_RIGHT) {
892: PetscCall(PCApplyTranspose(pc, x, work));
893: PetscCall(MatMultTranspose(pc->mat, work, y));
894: } else if (side == PC_LEFT) {
895: PetscCall(MatMultTranspose(pc->mat, x, work));
896: PetscCall(PCApplyTranspose(pc, work, y));
897: }
898: /* add support for PC_SYMMETRIC */
899: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 4, PETSC_FALSE));
900: PetscFunctionReturn(PETSC_SUCCESS);
901: }
903: /*@
904: PCApplyRichardsonExists - Determines whether a particular preconditioner has a
905: built-in fast application of Richardson's method.
907: Not Collective
909: Input Parameter:
910: . pc - the preconditioner
912: Output Parameter:
913: . exists - `PETSC_TRUE` or `PETSC_FALSE`
915: Level: developer
917: .seealso: [](ch_ksp), `PC`, `KSPRICHARDSON`, `PCApplyRichardson()`
918: @*/
919: PetscErrorCode PCApplyRichardsonExists(PC pc, PetscBool *exists)
920: {
921: PetscFunctionBegin;
923: PetscAssertPointer(exists, 2);
924: if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
925: else *exists = PETSC_FALSE;
926: PetscFunctionReturn(PETSC_SUCCESS);
927: }
929: /*@
930: PCApplyRichardson - Applies several steps of Richardson iteration with
931: the particular preconditioner. This routine is usually used by the
932: Krylov solvers and not the application code directly.
934: Collective
936: Input Parameters:
937: + pc - the `PC` preconditioner context
938: . b - the right-hand side
939: . w - one work vector
940: . rtol - relative decrease in residual norm convergence criteria
941: . abstol - absolute residual norm convergence criteria
942: . dtol - divergence residual norm increase criteria
943: . its - the number of iterations to apply.
944: - guesszero - if the input x contains nonzero initial guess
946: Output Parameters:
947: + outits - number of iterations actually used (for SOR this always equals its)
948: . reason - the reason the apply terminated
949: - y - the solution (also contains initial guess if guesszero is `PETSC_FALSE`
951: Level: developer
953: Notes:
954: Most preconditioners do not support this function. Use the command
955: `PCApplyRichardsonExists()` to determine if one does.
957: Except for the `PCMG` this routine ignores the convergence tolerances
958: and always runs for the number of iterations
960: .seealso: [](ch_ksp), `PC`, `PCApplyRichardsonExists()`
961: @*/
962: PetscErrorCode PCApplyRichardson(PC pc, Vec b, Vec y, Vec w, PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt its, PetscBool guesszero, PetscInt *outits, PCRichardsonConvergedReason *reason)
963: {
964: PetscFunctionBegin;
969: PetscCheck(b != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "b and y must be different vectors");
970: PetscCall(PCSetUp(pc));
971: PetscUseTypeMethod(pc, applyrichardson, b, y, w, rtol, abstol, dtol, its, guesszero, outits, reason);
972: PetscFunctionReturn(PETSC_SUCCESS);
973: }
975: /*@
976: PCSetFailedReason - Sets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail
978: Logically Collective
980: Input Parameters:
981: + pc - the `PC` preconditioner context
982: - reason - the reason it failed
984: Level: advanced
986: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCFailedReason`
987: @*/
988: PetscErrorCode PCSetFailedReason(PC pc, PCFailedReason reason)
989: {
990: PetscFunctionBegin;
992: pc->failedreason = reason;
993: PetscFunctionReturn(PETSC_SUCCESS);
994: }
996: /*@
997: PCGetFailedReason - Gets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail
999: Not Collective
1001: Input Parameter:
1002: . pc - the `PC` preconditioner context
1004: Output Parameter:
1005: . reason - the reason it failed
1007: Level: advanced
1009: Note:
1010: After a call to `KSPCheckDot()` or `KSPCheckNorm()` inside a `KSPSolve()` or a call to `PCReduceFailedReason()`
1011: this is the maximum reason over all MPI processes in the `PC` communicator and hence logically collective.
1012: Otherwise it returns the local value.
1014: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCSetFailedReason()`, `PCFailedReason`
1015: @*/
1016: PetscErrorCode PCGetFailedReason(PC pc, PCFailedReason *reason)
1017: {
1018: PetscFunctionBegin;
1020: if (pc->setupcalled < 0) *reason = (PCFailedReason)pc->setupcalled;
1021: else *reason = pc->failedreason;
1022: PetscFunctionReturn(PETSC_SUCCESS);
1023: }
1025: /*@
1026: PCReduceFailedReason - Reduce the failed reason among the MPI processes that share the `PC`
1028: Collective
1030: Input Parameter:
1031: . pc - the `PC` preconditioner context
1033: Level: advanced
1035: Note:
1036: Different MPI processes may have different reasons or no reason, see `PCGetFailedReason()`. This routine
1037: makes them have a common value (failure if any MPI process had a failure).
1039: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCGetFailedReason()`, `PCSetFailedReason()`, `PCFailedReason`
1040: @*/
1041: PetscErrorCode PCReduceFailedReason(PC pc)
1042: {
1043: PetscInt buf;
1045: PetscFunctionBegin;
1047: buf = (PetscInt)pc->failedreason;
1048: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &buf, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)pc)));
1049: pc->failedreason = (PCFailedReason)buf;
1050: PetscFunctionReturn(PETSC_SUCCESS);
1051: }
1053: /*
1054: a setupcall of 0 indicates never setup,
1055: 1 indicates has been previously setup
1056: -1 indicates a PCSetUp() was attempted and failed
1057: */
1058: /*@
1059: PCSetUp - Prepares for the use of a preconditioner. Performs all the one-time operations needed before the preconditioner
1060: can be used with `PCApply()`
1062: Collective
1064: Input Parameter:
1065: . pc - the `PC` preconditioner context
1067: Level: developer
1069: Notes:
1070: For example, for `PCLU` this will compute the factorization.
1072: This is called automatically by `KSPSetUp()` or `PCApply()` so rarely needs to be called directly.
1074: For nested preconditioners, such as `PCFIELDSPLIT` or `PCBJACOBI` this may not finish the construction of the preconditioner
1075: on the inner levels, the routine `PCSetUpOnBlocks()` may compute more of the preconditioner in those situations.
1077: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `KSPSetUp()`, `PCSetUpOnBlocks()`
1078: @*/
1079: PetscErrorCode PCSetUp(PC pc)
1080: {
1081: const char *def;
1082: PetscObjectState matstate, matnonzerostate;
1084: PetscFunctionBegin;
1086: PetscCheck(pc->mat, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Matrix must be set first");
1088: if (pc->setupcalled && pc->reusepreconditioner) {
1089: PetscCall(PetscInfo(pc, "Leaving PC with identical preconditioner since reuse preconditioner is set\n"));
1090: PetscFunctionReturn(PETSC_SUCCESS);
1091: }
1093: PetscCall(PetscObjectStateGet((PetscObject)pc->pmat, &matstate));
1094: PetscCall(MatGetNonzeroState(pc->pmat, &matnonzerostate));
1095: if (!pc->setupcalled) {
1096: //PetscCall(PetscInfo(pc, "Setting up PC for first time\n"));
1097: pc->flag = DIFFERENT_NONZERO_PATTERN;
1098: } else if (matstate == pc->matstate) PetscFunctionReturn(PETSC_SUCCESS);
1099: else {
1100: if (matnonzerostate != pc->matnonzerostate) {
1101: PetscCall(PetscInfo(pc, "Setting up PC with different nonzero pattern\n"));
1102: pc->flag = DIFFERENT_NONZERO_PATTERN;
1103: } else {
1104: //PetscCall(PetscInfo(pc, "Setting up PC with same nonzero pattern\n"));
1105: pc->flag = SAME_NONZERO_PATTERN;
1106: }
1107: }
1108: pc->matstate = matstate;
1109: pc->matnonzerostate = matnonzerostate;
1111: if (!((PetscObject)pc)->type_name) {
1112: PetscCall(PCGetDefaultType_Private(pc, &def));
1113: PetscCall(PCSetType(pc, def));
1114: }
1116: PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
1117: PetscCall(MatSetErrorIfFailure(pc->mat, pc->erroriffailure));
1118: PetscCall(PetscLogEventBegin(PC_SetUp, pc, 0, 0, 0));
1119: if (pc->ops->setup) {
1120: PetscCall(PCLogEventsDeactivatePush());
1121: PetscUseTypeMethod(pc, setup);
1122: PetscCall(PCLogEventsDeactivatePop());
1123: }
1124: PetscCall(PetscLogEventEnd(PC_SetUp, pc, 0, 0, 0));
1125: if (!pc->setupcalled) pc->setupcalled = 1;
1126: PetscFunctionReturn(PETSC_SUCCESS);
1127: }
1129: /*@
1130: PCSetUpOnBlocks - Sets up the preconditioner for each block in
1131: the block Jacobi, overlapping Schwarz, and fieldsplit methods.
1133: Collective
1135: Input Parameter:
1136: . pc - the `PC` preconditioner context
1138: Level: developer
1140: Notes:
1141: For nested preconditioners such as `PCBJACOBI`, `PCSetUp()` is not called on each sub-`KSP` when `PCSetUp()` is
1142: called on the outer `PC`, this routine ensures it is called.
1144: It calls `PCSetUp()` if not yet called.
1146: .seealso: [](ch_ksp), `PC`, `PCSetUp()`, `PCCreate()`, `PCApply()`, `PCDestroy()`
1147: @*/
1148: PetscErrorCode PCSetUpOnBlocks(PC pc)
1149: {
1150: PetscFunctionBegin;
1152: if (!pc->setupcalled) PetscCall(PCSetUp(pc)); /* "if" to prevent -info extra prints */
1153: if (!pc->ops->setuponblocks) PetscFunctionReturn(PETSC_SUCCESS);
1154: PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
1155: PetscCall(PetscLogEventBegin(PC_SetUpOnBlocks, pc, 0, 0, 0));
1156: PetscCall(PCLogEventsDeactivatePush());
1157: PetscUseTypeMethod(pc, setuponblocks);
1158: PetscCall(PCLogEventsDeactivatePop());
1159: PetscCall(PetscLogEventEnd(PC_SetUpOnBlocks, pc, 0, 0, 0));
1160: PetscFunctionReturn(PETSC_SUCCESS);
1161: }
1163: /*@C
1164: PCSetModifySubMatrices - Sets a user-defined routine for modifying the
1165: submatrices that arise within certain subdomain-based preconditioners such as `PCASM`
1167: Logically Collective
1169: Input Parameters:
1170: + pc - the `PC` preconditioner context
1171: . func - routine for modifying the submatrices
1172: - ctx - optional user-defined context (may be `NULL`)
1174: Calling sequence of `func`:
1175: + pc - the `PC` preconditioner context
1176: . nsub - number of index sets
1177: . row - an array of index sets that contain the global row numbers
1178: that comprise each local submatrix
1179: . col - an array of index sets that contain the global column numbers
1180: that comprise each local submatrix
1181: . submat - array of local submatrices
1182: - ctx - optional user-defined context for private data for the
1183: user-defined func routine (may be `NULL`)
1185: Level: advanced
1187: Notes:
1188: The basic submatrices are extracted from the preconditioner matrix as
1189: usual; the user can then alter these (for example, to set different boundary
1190: conditions for each submatrix) before they are used for the local solves.
1192: `PCSetModifySubMatrices()` MUST be called before `KSPSetUp()` and
1193: `KSPSolve()`.
1195: A routine set by `PCSetModifySubMatrices()` is currently called within
1196: the block Jacobi (`PCBJACOBI`) and additive Schwarz (`PCASM`)
1197: preconditioners. All other preconditioners ignore this routine.
1199: .seealso: [](ch_ksp), `PC`, `PCBJACOBI`, `PCASM`, `PCModifySubMatrices()`
1200: @*/
1201: PetscErrorCode PCSetModifySubMatrices(PC pc, PetscErrorCode (*func)(PC pc, PetscInt nsub, const IS row[], const IS col[], Mat submat[], void *ctx), void *ctx)
1202: {
1203: PetscFunctionBegin;
1205: pc->modifysubmatrices = func;
1206: pc->modifysubmatricesP = ctx;
1207: PetscFunctionReturn(PETSC_SUCCESS);
1208: }
1210: /*@C
1211: PCModifySubMatrices - Calls an optional user-defined routine within
1212: certain preconditioners if one has been set with `PCSetModifySubMatrices()`.
1214: Collective
1216: Input Parameters:
1217: + pc - the `PC` preconditioner context
1218: . nsub - the number of local submatrices
1219: . row - an array of index sets that contain the global row numbers
1220: that comprise each local submatrix
1221: . col - an array of index sets that contain the global column numbers
1222: that comprise each local submatrix
1223: . submat - array of local submatrices
1224: - ctx - optional user-defined context for private data for the
1225: user-defined routine (may be `NULL`)
1227: Output Parameter:
1228: . submat - array of local submatrices (the entries of which may
1229: have been modified)
1231: Level: developer
1233: Note:
1234: The user should NOT generally call this routine, as it will
1235: automatically be called within certain preconditioners.
1237: .seealso: [](ch_ksp), `PC`, `PCSetModifySubMatrices()`
1238: @*/
1239: PetscErrorCode PCModifySubMatrices(PC pc, PetscInt nsub, const IS row[], const IS col[], Mat submat[], void *ctx)
1240: {
1241: PetscFunctionBegin;
1243: if (!pc->modifysubmatrices) PetscFunctionReturn(PETSC_SUCCESS);
1244: PetscCall(PetscLogEventBegin(PC_ModifySubMatrices, pc, 0, 0, 0));
1245: PetscCall((*pc->modifysubmatrices)(pc, nsub, row, col, submat, ctx));
1246: PetscCall(PetscLogEventEnd(PC_ModifySubMatrices, pc, 0, 0, 0));
1247: PetscFunctionReturn(PETSC_SUCCESS);
1248: }
1250: /*@
1251: PCSetOperators - Sets the matrix associated with the linear system and
1252: a (possibly) different one associated with the preconditioner.
1254: Logically Collective
1256: Input Parameters:
1257: + pc - the `PC` preconditioner context
1258: . Amat - the matrix that defines the linear system
1259: - Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
1261: Level: intermediate
1263: Notes:
1264: Passing a `NULL` for `Amat` or `Pmat` removes the matrix that is currently used.
1266: If you wish to replace either `Amat` or `Pmat` but leave the other one untouched then
1267: first call `KSPGetOperators()` to get the one you wish to keep, call `PetscObjectReference()`
1268: on it and then pass it back in in your call to `KSPSetOperators()`.
1270: More Notes about Repeated Solution of Linear Systems:
1271: PETSc does NOT reset the matrix entries of either `Amat` or `Pmat`
1272: to zero after a linear solve; the user is completely responsible for
1273: matrix assembly. See the routine `MatZeroEntries()` if desiring to
1274: zero all elements of a matrix.
1276: .seealso: [](ch_ksp), `PC`, `PCGetOperators()`, `MatZeroEntries()`
1277: @*/
1278: PetscErrorCode PCSetOperators(PC pc, Mat Amat, Mat Pmat)
1279: {
1280: PetscInt m1, n1, m2, n2;
1282: PetscFunctionBegin;
1286: if (Amat) PetscCheckSameComm(pc, 1, Amat, 2);
1287: if (Pmat) PetscCheckSameComm(pc, 1, Pmat, 3);
1288: if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
1289: PetscCall(MatGetLocalSize(Amat, &m1, &n1));
1290: PetscCall(MatGetLocalSize(pc->mat, &m2, &n2));
1291: 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);
1292: PetscCall(MatGetLocalSize(Pmat, &m1, &n1));
1293: PetscCall(MatGetLocalSize(pc->pmat, &m2, &n2));
1294: 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);
1295: }
1297: if (Pmat != pc->pmat) {
1298: /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
1299: pc->matnonzerostate = -1;
1300: pc->matstate = -1;
1301: }
1303: /* reference first in case the matrices are the same */
1304: if (Amat) PetscCall(PetscObjectReference((PetscObject)Amat));
1305: PetscCall(MatDestroy(&pc->mat));
1306: if (Pmat) PetscCall(PetscObjectReference((PetscObject)Pmat));
1307: PetscCall(MatDestroy(&pc->pmat));
1308: pc->mat = Amat;
1309: pc->pmat = Pmat;
1310: PetscFunctionReturn(PETSC_SUCCESS);
1311: }
1313: /*@
1314: PCSetReusePreconditioner - reuse the current preconditioner even if the operator in the preconditioner `PC` has changed.
1316: Logically Collective
1318: Input Parameters:
1319: + pc - the `PC` preconditioner context
1320: - flag - `PETSC_TRUE` do not compute a new preconditioner, `PETSC_FALSE` do compute a new preconditioner
1322: Level: intermediate
1324: Note:
1325: Normally if a matrix inside a `PC` changes the `PC` automatically updates itself using information from the changed matrix. This option
1326: prevents this.
1328: .seealso: [](ch_ksp), `PC`, `PCGetOperators()`, `MatZeroEntries()`, `PCGetReusePreconditioner()`, `KSPSetReusePreconditioner()`
1329: @*/
1330: PetscErrorCode PCSetReusePreconditioner(PC pc, PetscBool flag)
1331: {
1332: PetscFunctionBegin;
1335: pc->reusepreconditioner = flag;
1336: PetscFunctionReturn(PETSC_SUCCESS);
1337: }
1339: /*@
1340: PCGetReusePreconditioner - Determines if the `PC` reuses the current preconditioner even if the operator in the preconditioner has changed.
1342: Not Collective
1344: Input Parameter:
1345: . pc - the `PC` preconditioner context
1347: Output Parameter:
1348: . flag - `PETSC_TRUE` do not compute a new preconditioner, `PETSC_FALSE` do compute a new preconditioner
1350: Level: intermediate
1352: .seealso: [](ch_ksp), `PC`, `PCGetOperators()`, `MatZeroEntries()`, `PCSetReusePreconditioner()`
1353: @*/
1354: PetscErrorCode PCGetReusePreconditioner(PC pc, PetscBool *flag)
1355: {
1356: PetscFunctionBegin;
1358: PetscAssertPointer(flag, 2);
1359: *flag = pc->reusepreconditioner;
1360: PetscFunctionReturn(PETSC_SUCCESS);
1361: }
1363: /*@
1364: PCGetOperators - Gets the matrix associated with the linear system and
1365: possibly a different one which is used to construct the preconditioner.
1367: Not Collective, though parallel `Mat`s are returned if `pc` is parallel
1369: Input Parameter:
1370: . pc - the `PC` preconditioner context
1372: Output Parameters:
1373: + Amat - the matrix defining the linear system
1374: - Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.
1376: Level: intermediate
1378: Note:
1379: Does not increase the reference count of the matrices, so you should not destroy them
1381: Alternative usage: If the operators have NOT been set with `KSPSetOperators()`/`PCSetOperators()` then the operators
1382: are created in `PC` and returned to the user. In this case, if both operators
1383: mat and pmat are requested, two DIFFERENT operators will be returned. If
1384: only one is requested both operators in the PC will be the same (i.e. as
1385: if one had called `KSPSetOperators()`/`PCSetOperators()` with the same argument for both Mats).
1386: The user must set the sizes of the returned matrices and their type etc just
1387: as if the user created them with `MatCreate()`. For example,
1389: .vb
1390: KSP/PCGetOperators(ksp/pc,&Amat,NULL); is equivalent to
1391: set size, type, etc of Amat
1393: MatCreate(comm,&mat);
1394: KSP/PCSetOperators(ksp/pc,Amat,Amat);
1395: PetscObjectDereference((PetscObject)mat);
1396: set size, type, etc of Amat
1397: .ve
1399: and
1401: .vb
1402: KSP/PCGetOperators(ksp/pc,&Amat,&Pmat); is equivalent to
1403: set size, type, etc of Amat and Pmat
1405: MatCreate(comm,&Amat);
1406: MatCreate(comm,&Pmat);
1407: KSP/PCSetOperators(ksp/pc,Amat,Pmat);
1408: PetscObjectDereference((PetscObject)Amat);
1409: PetscObjectDereference((PetscObject)Pmat);
1410: set size, type, etc of Amat and Pmat
1411: .ve
1413: The rationale for this support is so that when creating a `TS`, `SNES`, or `KSP` the hierarchy
1414: of underlying objects (i.e. `SNES`, `KSP`, `PC`, `Mat`) and their lifespans can be completely
1415: managed by the top most level object (i.e. the `TS`, `SNES`, or `KSP`). Another way to look
1416: at this is when you create a `SNES` you do not NEED to create a `KSP` and attach it to
1417: the `SNES` object (the `SNES` object manages it for you). Similarly when you create a KSP
1418: you do not need to attach a `PC` to it (the `KSP` object manages the `PC` object for you).
1419: Thus, why should YOU have to create the `Mat` and attach it to the `SNES`/`KSP`/`PC`, when
1420: it can be created for you?
1422: .seealso: [](ch_ksp), `PC`, `PCSetOperators()`, `KSPGetOperators()`, `KSPSetOperators()`, `PCGetOperatorsSet()`
1423: @*/
1424: PetscErrorCode PCGetOperators(PC pc, Mat *Amat, Mat *Pmat)
1425: {
1426: PetscFunctionBegin;
1428: if (Amat) {
1429: if (!pc->mat) {
1430: if (pc->pmat && !Pmat) { /* Pmat has been set, but user did not request it, so use for Amat */
1431: pc->mat = pc->pmat;
1432: PetscCall(PetscObjectReference((PetscObject)pc->mat));
1433: } else { /* both Amat and Pmat are empty */
1434: PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &pc->mat));
1435: if (!Pmat) { /* user did NOT request Pmat, so make same as Amat */
1436: pc->pmat = pc->mat;
1437: PetscCall(PetscObjectReference((PetscObject)pc->pmat));
1438: }
1439: }
1440: }
1441: *Amat = pc->mat;
1442: }
1443: if (Pmat) {
1444: if (!pc->pmat) {
1445: if (pc->mat && !Amat) { /* Amat has been set but was not requested, so use for pmat */
1446: pc->pmat = pc->mat;
1447: PetscCall(PetscObjectReference((PetscObject)pc->pmat));
1448: } else {
1449: PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &pc->pmat));
1450: if (!Amat) { /* user did NOT request Amat, so make same as Pmat */
1451: pc->mat = pc->pmat;
1452: PetscCall(PetscObjectReference((PetscObject)pc->mat));
1453: }
1454: }
1455: }
1456: *Pmat = pc->pmat;
1457: }
1458: PetscFunctionReturn(PETSC_SUCCESS);
1459: }
1461: /*@
1462: PCGetOperatorsSet - Determines if the matrix associated with the linear system and
1463: possibly a different one associated with the preconditioner have been set in the `PC`.
1465: Not Collective, though the results on all processes should be the same
1467: Input Parameter:
1468: . pc - the `PC` preconditioner context
1470: Output Parameters:
1471: + mat - the matrix associated with the linear system was set
1472: - pmat - matrix associated with the preconditioner was set, usually the same
1474: Level: intermediate
1476: .seealso: [](ch_ksp), `PC`, `PCSetOperators()`, `KSPGetOperators()`, `KSPSetOperators()`, `PCGetOperators()`
1477: @*/
1478: PetscErrorCode PCGetOperatorsSet(PC pc, PetscBool *mat, PetscBool *pmat)
1479: {
1480: PetscFunctionBegin;
1482: if (mat) *mat = (pc->mat) ? PETSC_TRUE : PETSC_FALSE;
1483: if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1484: PetscFunctionReturn(PETSC_SUCCESS);
1485: }
1487: /*@
1488: PCFactorGetMatrix - Gets the factored matrix from the
1489: preconditioner context. This routine is valid only for the `PCLU`,
1490: `PCILU`, `PCCHOLESKY`, and `PCICC` methods.
1492: Not Collective though `mat` is parallel if `pc` is parallel
1494: Input Parameter:
1495: . pc - the `PC` preconditioner context
1497: Output Parameters:
1498: . mat - the factored matrix
1500: Level: advanced
1502: Note:
1503: Does not increase the reference count for `mat` so DO NOT destroy it
1505: .seealso: [](ch_ksp), `PC`, `PCLU`, `PCILU`, `PCCHOLESKY`, `PCICC`
1506: @*/
1507: PetscErrorCode PCFactorGetMatrix(PC pc, Mat *mat)
1508: {
1509: PetscFunctionBegin;
1511: PetscAssertPointer(mat, 2);
1512: PetscCall(PCFactorSetUpMatSolverType(pc));
1513: PetscUseTypeMethod(pc, getfactoredmatrix, mat);
1514: PetscFunctionReturn(PETSC_SUCCESS);
1515: }
1517: /*@
1518: PCSetOptionsPrefix - Sets the prefix used for searching for all
1519: `PC` options in the database.
1521: Logically Collective
1523: Input Parameters:
1524: + pc - the `PC` preconditioner context
1525: - prefix - the prefix string to prepend to all `PC` option requests
1527: Note:
1528: A hyphen (-) must NOT be given at the beginning of the prefix name.
1529: The first character of all runtime options is AUTOMATICALLY the
1530: hyphen.
1532: Level: advanced
1534: .seealso: [](ch_ksp), `PC`, `PCSetFromOptions`, `PCAppendOptionsPrefix()`, `PCGetOptionsPrefix()`
1535: @*/
1536: PetscErrorCode PCSetOptionsPrefix(PC pc, const char prefix[])
1537: {
1538: PetscFunctionBegin;
1540: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc, prefix));
1541: PetscFunctionReturn(PETSC_SUCCESS);
1542: }
1544: /*@
1545: PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1546: `PC` options in the database.
1548: Logically Collective
1550: Input Parameters:
1551: + pc - the `PC` preconditioner context
1552: - prefix - the prefix string to prepend to all `PC` option requests
1554: Note:
1555: A hyphen (-) must NOT be given at the beginning of the prefix name.
1556: The first character of all runtime options is AUTOMATICALLY the
1557: hyphen.
1559: Level: advanced
1561: .seealso: [](ch_ksp), `PC`, `PCSetFromOptions`, `PCSetOptionsPrefix()`, `PCGetOptionsPrefix()`
1562: @*/
1563: PetscErrorCode PCAppendOptionsPrefix(PC pc, const char prefix[])
1564: {
1565: PetscFunctionBegin;
1567: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)pc, prefix));
1568: PetscFunctionReturn(PETSC_SUCCESS);
1569: }
1571: /*@
1572: PCGetOptionsPrefix - Gets the prefix used for searching for all
1573: PC options in the database.
1575: Not Collective
1577: Input Parameter:
1578: . pc - the `PC` preconditioner context
1580: Output Parameter:
1581: . prefix - pointer to the prefix string used, is returned
1583: Level: advanced
1585: .seealso: [](ch_ksp), `PC`, `PCSetFromOptions`, `PCSetOptionsPrefix()`, `PCAppendOptionsPrefix()`
1586: @*/
1587: PetscErrorCode PCGetOptionsPrefix(PC pc, const char *prefix[])
1588: {
1589: PetscFunctionBegin;
1591: PetscAssertPointer(prefix, 2);
1592: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, prefix));
1593: PetscFunctionReturn(PETSC_SUCCESS);
1594: }
1596: /*
1597: Indicates the right-hand side will be changed by KSPSolve(), this occurs for a few
1598: preconditioners including BDDC and Eisentat that transform the equations before applying
1599: the Krylov methods
1600: */
1601: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC pc, PetscBool *change)
1602: {
1603: PetscFunctionBegin;
1605: PetscAssertPointer(change, 2);
1606: *change = PETSC_FALSE;
1607: PetscTryMethod(pc, "PCPreSolveChangeRHS_C", (PC, PetscBool *), (pc, change));
1608: PetscFunctionReturn(PETSC_SUCCESS);
1609: }
1611: /*@
1612: PCPreSolve - Optional pre-solve phase, intended for any preconditioner-specific actions that must be performed before
1613: the iterative solve itself. Used in conjunction with `PCPostSolve()`
1615: Collective
1617: Input Parameters:
1618: + pc - the `PC` preconditioner context
1619: - ksp - the Krylov subspace context
1621: Level: developer
1623: Notes:
1624: `KSPSolve()` calls this directly, so is rarely called by the user.
1626: Certain preconditioners, such as the `PCType` of `PCEISENSTAT`, change the formulation of the linear system to be solved iteratively.
1627: This function performs that transformation. `PCPostSolve()` then transforms the system back to its original form after the solve.
1628: `PCPostSolve()` also transforms the resulting solution of the transformed system to the solution of the original problem.
1630: `KSPSetPreSolve()` and `KSPSetPostSolve()` provide an alternative way to provide such transformations.
1632: .seealso: [](ch_ksp), `PC`, `PCPostSolve()`, `KSP`, `PCSetPreSolve()`, `KSPSetPreSolve()`, `KSPSetPostSolve()`
1633: @*/
1634: PetscErrorCode PCPreSolve(PC pc, KSP ksp)
1635: {
1636: Vec x, rhs;
1638: PetscFunctionBegin;
1641: pc->presolvedone++;
1642: PetscCheck(pc->presolvedone <= 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot embed PCPreSolve() more than twice");
1643: PetscCall(KSPGetSolution(ksp, &x));
1644: PetscCall(KSPGetRhs(ksp, &rhs));
1646: if (pc->ops->presolve) PetscUseTypeMethod(pc, presolve, ksp, rhs, x);
1647: else if (pc->presolve) PetscCall(pc->presolve(pc, ksp));
1648: PetscFunctionReturn(PETSC_SUCCESS);
1649: }
1651: /*@C
1652: PCSetPreSolve - Sets function used by `PCPreSolve()` which is intended for any
1653: preconditioner-specific actions that must be performed before
1654: the iterative solve itself.
1656: Logically Collective
1658: Input Parameters:
1659: + pc - the preconditioner object
1660: - presolve - the function to call before the solve
1662: Calling sequence of `presolve`:
1663: + pc - the `PC` context
1664: - ksp - the `KSP` context
1666: Level: developer
1668: .seealso: [](ch_ksp), `PC`, `PCSetUp()`, `PCPreSolve()`
1669: @*/
1670: PetscErrorCode PCSetPreSolve(PC pc, PetscErrorCode (*presolve)(PC pc, KSP ksp))
1671: {
1672: PetscFunctionBegin;
1674: pc->presolve = presolve;
1675: PetscFunctionReturn(PETSC_SUCCESS);
1676: }
1678: /*@
1679: PCPostSolve - Optional post-solve phase, intended for any
1680: preconditioner-specific actions that must be performed after
1681: the iterative solve itself.
1683: Collective
1685: Input Parameters:
1686: + pc - the `PC` preconditioner context
1687: - ksp - the `KSP` Krylov subspace context
1689: Example Usage:
1690: .vb
1691: PCPreSolve(pc,ksp);
1692: KSPSolve(ksp,b,x);
1693: PCPostSolve(pc,ksp);
1694: .ve
1696: Level: developer
1698: Note:
1699: `KSPSolve()` calls this routine directly, so it is rarely called by the user.
1701: .seealso: [](ch_ksp), `PC`, `PCSetPreSolve()`, `KSPSetPostSolve()`, `KSPSetPreSolve()`, `PCPreSolve()`, `KSPSolve()`
1702: @*/
1703: PetscErrorCode PCPostSolve(PC pc, KSP ksp)
1704: {
1705: Vec x, rhs;
1707: PetscFunctionBegin;
1710: pc->presolvedone--;
1711: PetscCall(KSPGetSolution(ksp, &x));
1712: PetscCall(KSPGetRhs(ksp, &rhs));
1713: PetscTryTypeMethod(pc, postsolve, ksp, rhs, x);
1714: PetscFunctionReturn(PETSC_SUCCESS);
1715: }
1717: /*@
1718: PCLoad - Loads a `PC` that has been stored in binary with `PCView()`.
1720: Collective
1722: Input Parameters:
1723: + newdm - the newly loaded `PC`, this needs to have been created with `PCCreate()` or
1724: some related function before a call to `PCLoad()`.
1725: - viewer - binary file viewer `PETSCVIEWERBINARY`, obtained from `PetscViewerBinaryOpen()`
1727: Level: intermediate
1729: Note:
1730: The type is determined by the data in the file, any `PCType` set into the `PC` before this call is ignored.
1732: .seealso: [](ch_ksp), `PC`, `PetscViewerBinaryOpen()`, `PCView()`, `MatLoad()`, `VecLoad()`, `PETSCVIEWERBINARY`
1733: @*/
1734: PetscErrorCode PCLoad(PC newdm, PetscViewer viewer)
1735: {
1736: PetscBool isbinary;
1737: PetscInt classid;
1738: char type[256];
1740: PetscFunctionBegin;
1743: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1744: PetscCheck(isbinary, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1746: PetscCall(PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT));
1747: PetscCheck(classid == PC_FILE_CLASSID, PetscObjectComm((PetscObject)newdm), PETSC_ERR_ARG_WRONG, "Not PC next in file");
1748: PetscCall(PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR));
1749: PetscCall(PCSetType(newdm, type));
1750: PetscTryTypeMethod(newdm, load, viewer);
1751: PetscFunctionReturn(PETSC_SUCCESS);
1752: }
1754: #include <petscdraw.h>
1755: #if defined(PETSC_HAVE_SAWS)
1756: #include <petscviewersaws.h>
1757: #endif
1759: /*@
1760: PCViewFromOptions - View (print or provide information about) the `PC`, based on options in the options database
1762: Collective
1764: Input Parameters:
1765: + A - the `PC` context
1766: . obj - Optional object that provides the options prefix
1767: - name - command line option name
1769: Level: developer
1771: .seealso: [](ch_ksp), `PC`, `PCView`, `PetscObjectViewFromOptions()`, `PCCreate()`
1772: @*/
1773: PetscErrorCode PCViewFromOptions(PC A, PetscObject obj, const char name[])
1774: {
1775: PetscFunctionBegin;
1777: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
1778: PetscFunctionReturn(PETSC_SUCCESS);
1779: }
1781: /*@
1782: PCView - Prints information about the `PC`
1784: Collective
1786: Input Parameters:
1787: + pc - the `PC` preconditioner context
1788: - viewer - optional `PetscViewer` visualization context
1790: Level: intermediate
1792: Notes:
1793: The available visualization contexts include
1794: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
1795: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
1796: output where only the first processor opens
1797: the file. All other processors send their
1798: data to the first processor to print.
1800: The user can open an alternative visualization contexts with
1801: `PetscViewerASCIIOpen()` (output to a specified file).
1803: .seealso: [](ch_ksp), `PC`, `PetscViewer`, `PetscViewerType`, `KSPView()`, `PetscViewerASCIIOpen()`
1804: @*/
1805: PetscErrorCode PCView(PC pc, PetscViewer viewer)
1806: {
1807: PCType cstr;
1808: PetscViewerFormat format;
1809: PetscBool iascii, isstring, isbinary, isdraw, pop = PETSC_FALSE;
1810: #if defined(PETSC_HAVE_SAWS)
1811: PetscBool issaws;
1812: #endif
1814: PetscFunctionBegin;
1816: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc), &viewer));
1818: PetscCheckSameComm(pc, 1, viewer, 2);
1820: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1821: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
1822: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1823: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1824: #if defined(PETSC_HAVE_SAWS)
1825: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws));
1826: #endif
1828: if (iascii) {
1829: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)pc, viewer));
1830: if (!pc->setupcalled) PetscCall(PetscViewerASCIIPrintf(viewer, " PC has not been set up so information may be incomplete\n"));
1831: PetscCall(PetscViewerASCIIPushTab(viewer));
1832: PetscTryTypeMethod(pc, view, viewer);
1833: PetscCall(PetscViewerASCIIPopTab(viewer));
1834: if (pc->mat) {
1835: PetscCall(PetscViewerGetFormat(viewer, &format));
1836: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
1837: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
1838: pop = PETSC_TRUE;
1839: }
1840: if (pc->pmat == pc->mat) {
1841: PetscCall(PetscViewerASCIIPrintf(viewer, " linear system matrix = precond matrix:\n"));
1842: PetscCall(PetscViewerASCIIPushTab(viewer));
1843: PetscCall(MatView(pc->mat, viewer));
1844: PetscCall(PetscViewerASCIIPopTab(viewer));
1845: } else {
1846: if (pc->pmat) {
1847: PetscCall(PetscViewerASCIIPrintf(viewer, " linear system matrix followed by preconditioner matrix:\n"));
1848: } else {
1849: PetscCall(PetscViewerASCIIPrintf(viewer, " linear system matrix:\n"));
1850: }
1851: PetscCall(PetscViewerASCIIPushTab(viewer));
1852: PetscCall(MatView(pc->mat, viewer));
1853: if (pc->pmat) PetscCall(MatView(pc->pmat, viewer));
1854: PetscCall(PetscViewerASCIIPopTab(viewer));
1855: }
1856: if (pop) PetscCall(PetscViewerPopFormat(viewer));
1857: }
1858: } else if (isstring) {
1859: PetscCall(PCGetType(pc, &cstr));
1860: PetscCall(PetscViewerStringSPrintf(viewer, " PCType: %-7.7s", cstr));
1861: PetscTryTypeMethod(pc, view, viewer);
1862: if (pc->mat) PetscCall(MatView(pc->mat, viewer));
1863: if (pc->pmat && pc->pmat != pc->mat) PetscCall(MatView(pc->pmat, viewer));
1864: } else if (isbinary) {
1865: PetscInt classid = PC_FILE_CLASSID;
1866: MPI_Comm comm;
1867: PetscMPIInt rank;
1868: char type[256];
1870: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
1871: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1872: if (rank == 0) {
1873: PetscCall(PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT));
1874: PetscCall(PetscStrncpy(type, ((PetscObject)pc)->type_name, 256));
1875: PetscCall(PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR));
1876: }
1877: PetscTryTypeMethod(pc, view, viewer);
1878: } else if (isdraw) {
1879: PetscDraw draw;
1880: char str[25];
1881: PetscReal x, y, bottom, h;
1882: PetscInt n;
1884: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1885: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
1886: if (pc->mat) {
1887: PetscCall(MatGetSize(pc->mat, &n, NULL));
1888: PetscCall(PetscSNPrintf(str, 25, "PC: %s (%" PetscInt_FMT ")", ((PetscObject)pc)->type_name, n));
1889: } else {
1890: PetscCall(PetscSNPrintf(str, 25, "PC: %s", ((PetscObject)pc)->type_name));
1891: }
1892: PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h));
1893: bottom = y - h;
1894: PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
1895: PetscTryTypeMethod(pc, view, viewer);
1896: PetscCall(PetscDrawPopCurrentPoint(draw));
1897: #if defined(PETSC_HAVE_SAWS)
1898: } else if (issaws) {
1899: PetscMPIInt rank;
1901: PetscCall(PetscObjectName((PetscObject)pc));
1902: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1903: if (!((PetscObject)pc)->amsmem && rank == 0) PetscCall(PetscObjectViewSAWs((PetscObject)pc, viewer));
1904: if (pc->mat) PetscCall(MatView(pc->mat, viewer));
1905: if (pc->pmat && pc->pmat != pc->mat) PetscCall(MatView(pc->pmat, viewer));
1906: #endif
1907: }
1908: PetscFunctionReturn(PETSC_SUCCESS);
1909: }
1911: /*@C
1912: PCRegister - Adds a method (`PCType`) to the PETSc preconditioner package.
1914: Not collective. No Fortran Support
1916: Input Parameters:
1917: + sname - name of a new user-defined solver
1918: - function - routine to create the method context which will be stored in a `PC` when `PCSetType()` is called
1920: Example Usage:
1921: .vb
1922: PCRegister("my_solver", MySolverCreate);
1923: .ve
1925: Then, your solver can be chosen with the procedural interface via
1926: $ PCSetType(pc, "my_solver")
1927: or at runtime via the option
1928: $ -pc_type my_solver
1930: Level: advanced
1932: Note:
1933: A simpler alternative to using `PCRegister()` for an application specific preconditioner is to use a `PC` of `PCType` `PCSHELL` and
1934: provide your customizations with `PCShellSetContext()` and `PCShellSetApply()`
1936: `PCRegister()` may be called multiple times to add several user-defined preconditioners.
1938: .seealso: [](ch_ksp), `PC`, `PCType`, `PCRegisterAll()`, `PCSetType()`, `PCShellSetContext()`, `PCShellSetApply()`, `PCSHELL`
1939: @*/
1940: PetscErrorCode PCRegister(const char sname[], PetscErrorCode (*function)(PC))
1941: {
1942: PetscFunctionBegin;
1943: PetscCall(PCInitializePackage());
1944: PetscCall(PetscFunctionListAdd(&PCList, sname, function));
1945: PetscFunctionReturn(PETSC_SUCCESS);
1946: }
1948: static PetscErrorCode MatMult_PC(Mat A, Vec X, Vec Y)
1949: {
1950: PC pc;
1952: PetscFunctionBegin;
1953: PetscCall(MatShellGetContext(A, &pc));
1954: PetscCall(PCApply(pc, X, Y));
1955: PetscFunctionReturn(PETSC_SUCCESS);
1956: }
1958: /*@
1959: PCComputeOperator - Computes the explicit preconditioned operator as a matrix `Mat`.
1961: Collective
1963: Input Parameters:
1964: + pc - the `PC` preconditioner object
1965: - mattype - the `MatType` to be used for the operator
1967: Output Parameter:
1968: . mat - the explicit preconditioned operator
1970: Level: advanced
1972: Note:
1973: This computation is done by applying the operators to columns of the identity matrix.
1974: This routine is costly in general, and is recommended for use only with relatively small systems.
1975: Currently, this routine uses a dense matrix format when `mattype` == `NULL`
1977: Developer Note:
1978: This should be called `PCCreateExplicitOperator()`
1980: .seealso: [](ch_ksp), `PC`, `KSPComputeOperator()`, `MatType`
1981: @*/
1982: PetscErrorCode PCComputeOperator(PC pc, MatType mattype, Mat *mat)
1983: {
1984: PetscInt N, M, m, n;
1985: Mat A, Apc;
1987: PetscFunctionBegin;
1989: PetscAssertPointer(mat, 3);
1990: PetscCall(PCGetOperators(pc, &A, NULL));
1991: PetscCall(MatGetLocalSize(A, &m, &n));
1992: PetscCall(MatGetSize(A, &M, &N));
1993: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)pc), m, n, M, N, pc, &Apc));
1994: PetscCall(MatShellSetOperation(Apc, MATOP_MULT, (void (*)(void))MatMult_PC));
1995: PetscCall(MatComputeOperator(Apc, mattype, mat));
1996: PetscCall(MatDestroy(&Apc));
1997: PetscFunctionReturn(PETSC_SUCCESS);
1998: }
2000: /*@
2001: PCSetCoordinates - sets the coordinates of all the nodes (degrees of freedom in the vector) on the local process
2003: Collective
2005: Input Parameters:
2006: + pc - the `PC` preconditioner context
2007: . dim - the dimension of the coordinates 1, 2, or 3
2008: . nloc - the blocked size of the coordinates array
2009: - coords - the coordinates array
2011: Level: intermediate
2013: Notes:
2014: `coords` is an array of the dim coordinates for the nodes on
2015: the local processor, of size `dim`*`nloc`.
2016: If there are 108 equations (dofs) on a processor
2017: for a 3d displacement finite element discretization of elasticity (so
2018: that there are nloc = 36 = 108/3 nodes) then the array must have 108
2019: double precision values (ie, 3 * 36). These x y z coordinates
2020: should be ordered for nodes 0 to N-1 like so: [ 0.x, 0.y, 0.z, 1.x,
2021: ... , N-1.z ].
2023: The information provided here can be used by some preconditioners, such as `PCGAMG`, to produce a better preconditioner.
2024: See also `MatSetNearNullSpace()`.
2026: .seealso: [](ch_ksp), `PC`, `MatSetNearNullSpace()`
2027: @*/
2028: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[])
2029: {
2030: PetscFunctionBegin;
2033: PetscTryMethod(pc, "PCSetCoordinates_C", (PC, PetscInt, PetscInt, PetscReal[]), (pc, dim, nloc, coords));
2034: PetscFunctionReturn(PETSC_SUCCESS);
2035: }
2037: /*@
2038: PCGetInterpolations - Gets interpolation matrices for all levels (except level 0)
2040: Logically Collective
2042: Input Parameter:
2043: . pc - the precondition context
2045: Output Parameters:
2046: + num_levels - the number of levels
2047: - interpolations - the interpolation matrices (size of `num_levels`-1)
2049: Level: advanced
2051: Developer Note:
2052: Why is this here instead of in `PCMG` etc?
2054: .seealso: [](ch_ksp), `PC`, `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetInterpolation()`, `PCGetCoarseOperators()`
2055: @*/
2056: PetscErrorCode PCGetInterpolations(PC pc, PetscInt *num_levels, Mat *interpolations[])
2057: {
2058: PetscFunctionBegin;
2060: PetscAssertPointer(num_levels, 2);
2061: PetscAssertPointer(interpolations, 3);
2062: PetscUseMethod(pc, "PCGetInterpolations_C", (PC, PetscInt *, Mat *[]), (pc, num_levels, interpolations));
2063: PetscFunctionReturn(PETSC_SUCCESS);
2064: }
2066: /*@
2067: PCGetCoarseOperators - Gets coarse operator matrices for all levels (except the finest level)
2069: Logically Collective
2071: Input Parameter:
2072: . pc - the precondition context
2074: Output Parameters:
2075: + num_levels - the number of levels
2076: - coarseOperators - the coarse operator matrices (size of `num_levels`-1)
2078: Level: advanced
2080: Developer Note:
2081: Why is this here instead of in `PCMG` etc?
2083: .seealso: [](ch_ksp), `PC`, `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetRScale()`, `PCMGGetInterpolation()`, `PCGetInterpolations()`
2084: @*/
2085: PetscErrorCode PCGetCoarseOperators(PC pc, PetscInt *num_levels, Mat *coarseOperators[])
2086: {
2087: PetscFunctionBegin;
2089: PetscAssertPointer(num_levels, 2);
2090: PetscAssertPointer(coarseOperators, 3);
2091: PetscUseMethod(pc, "PCGetCoarseOperators_C", (PC, PetscInt *, Mat *[]), (pc, num_levels, coarseOperators));
2092: PetscFunctionReturn(PETSC_SUCCESS);
2093: }