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_ApplyMultiple, 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 pcsetupcalled = 0 state and removes any allocated `Vec`s and `Mat`s
85: Collective
87: Input Parameter:
88: . pc - the preconditioner context
90: Level: developer
92: Note:
93: 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`
95: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetUp()`
96: @*/
97: PetscErrorCode PCReset(PC pc)
98: {
99: PetscFunctionBegin;
101: PetscTryTypeMethod(pc, reset);
102: PetscCall(VecDestroy(&pc->diagonalscaleright));
103: PetscCall(VecDestroy(&pc->diagonalscaleleft));
104: PetscCall(MatDestroy(&pc->pmat));
105: PetscCall(MatDestroy(&pc->mat));
107: pc->setupcalled = 0;
108: PetscFunctionReturn(PETSC_SUCCESS);
109: }
111: /*@
112: PCDestroy - Destroys `PC` context that was created with `PCCreate()`.
114: Collective
116: Input Parameter:
117: . pc - the preconditioner context
119: Level: developer
121: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetUp()`
122: @*/
123: PetscErrorCode PCDestroy(PC *pc)
124: {
125: PetscFunctionBegin;
126: if (!*pc) PetscFunctionReturn(PETSC_SUCCESS);
128: if (--((PetscObject)*pc)->refct > 0) {
129: *pc = NULL;
130: PetscFunctionReturn(PETSC_SUCCESS);
131: }
133: PetscCall(PCReset(*pc));
135: /* if memory was published with SAWs then destroy it */
136: PetscCall(PetscObjectSAWsViewOff((PetscObject)*pc));
137: PetscTryTypeMethod(*pc, destroy);
138: PetscCall(DMDestroy(&(*pc)->dm));
139: PetscCall(PetscHeaderDestroy(pc));
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
143: /*@
144: PCGetDiagonalScale - Indicates if the preconditioner applies an additional left and right
145: scaling as needed by certain time-stepping codes.
147: Logically Collective
149: Input Parameter:
150: . pc - the preconditioner context
152: Output Parameter:
153: . flag - `PETSC_TRUE` if it applies the scaling
155: Level: developer
157: Note:
158: If this returns `PETSC_TRUE` then the system solved via the Krylov method is, for left and right preconditioning,
160: $$
161: \begin{align*}
162: D M A D^{-1} y = D M b \\
163: D A M D^{-1} z = D b.
164: \end{align*}
165: $$
167: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCDiagonalScaleRight()`, `PCSetDiagonalScale()`
168: @*/
169: PetscErrorCode PCGetDiagonalScale(PC pc, PetscBool *flag)
170: {
171: PetscFunctionBegin;
173: PetscAssertPointer(flag, 2);
174: *flag = pc->diagonalscale;
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: /*@
179: PCSetDiagonalScale - Indicates the left scaling to use to apply an additional left and right
180: scaling as needed by certain time-stepping codes.
182: Logically Collective
184: Input Parameters:
185: + pc - the preconditioner context
186: - s - scaling vector
188: Level: intermediate
190: Notes:
191: The system solved via the Krylov method is, for left and right preconditioning,
192: $$
193: \begin{align*}
194: D M A D^{-1} y = D M b \\
195: D A M D^{-1} z = D b.
196: \end{align*}
197: $$
199: `PCDiagonalScaleLeft()` scales a vector by $D$. `PCDiagonalScaleRight()` scales a vector by $D^{-1}$.
201: .seealso: [](ch_ksp), `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCDiagonalScaleRight()`, `PCGetDiagonalScale()`
202: @*/
203: PetscErrorCode PCSetDiagonalScale(PC pc, Vec s)
204: {
205: PetscFunctionBegin;
208: pc->diagonalscale = PETSC_TRUE;
210: PetscCall(PetscObjectReference((PetscObject)s));
211: PetscCall(VecDestroy(&pc->diagonalscaleleft));
213: pc->diagonalscaleleft = s;
215: PetscCall(VecDuplicate(s, &pc->diagonalscaleright));
216: PetscCall(VecCopy(s, pc->diagonalscaleright));
217: PetscCall(VecReciprocal(pc->diagonalscaleright));
218: PetscFunctionReturn(PETSC_SUCCESS);
219: }
221: /*@
222: PCDiagonalScaleLeft - Scales a vector by the left scaling as needed by certain time-stepping codes.
224: Logically Collective
226: Input Parameters:
227: + pc - the preconditioner context
228: . in - input vector
229: - out - scaled vector (maybe the same as in)
231: Level: intermediate
233: Notes:
234: The system solved via the Krylov method is, for left and right preconditioning,
236: $$
237: \begin{align*}
238: D M A D^{-1} y = D M b \\
239: D A M D^{-1} z = D b.
240: \end{align*}
241: $$
243: `PCDiagonalScaleLeft()` scales a vector by $D$. `PCDiagonalScaleRight()` scales a vector by $D^{-1}$.
245: If diagonal scaling is turned off and `in` is not `out` then `in` is copied to `out`
247: .seealso: [](ch_ksp), `PCCreate()`, `PCSetUp()`, `PCSetDiagonalScale()`, `PCDiagonalScaleRight()`, `MatDiagonalScale()`
248: @*/
249: PetscErrorCode PCDiagonalScaleLeft(PC pc, Vec in, Vec out)
250: {
251: PetscFunctionBegin;
255: if (pc->diagonalscale) {
256: PetscCall(VecPointwiseMult(out, pc->diagonalscaleleft, in));
257: } else if (in != out) {
258: PetscCall(VecCopy(in, out));
259: }
260: PetscFunctionReturn(PETSC_SUCCESS);
261: }
263: /*@
264: PCDiagonalScaleRight - Scales a vector by the right scaling as needed by certain time-stepping codes.
266: Logically Collective
268: Input Parameters:
269: + pc - the preconditioner context
270: . in - input vector
271: - out - scaled vector (maybe the same as in)
273: Level: intermediate
275: Notes:
276: The system solved via the Krylov method is, for left and right preconditioning,
278: $$
279: \begin{align*}
280: D M A D^{-1} y = D M b \\
281: D A M D^{-1} z = D b.
282: \end{align*}
283: $$
285: `PCDiagonalScaleLeft()` scales a vector by $D$. `PCDiagonalScaleRight()` scales a vector by $D^{-1}$.
287: If diagonal scaling is turned off and `in` is not `out` then `in` is copied to `out`
289: .seealso: [](ch_ksp), `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCSetDiagonalScale()`, `MatDiagonalScale()`
290: @*/
291: PetscErrorCode PCDiagonalScaleRight(PC pc, Vec in, Vec out)
292: {
293: PetscFunctionBegin;
297: if (pc->diagonalscale) {
298: PetscCall(VecPointwiseMult(out, pc->diagonalscaleright, in));
299: } else if (in != out) {
300: PetscCall(VecCopy(in, out));
301: }
302: PetscFunctionReturn(PETSC_SUCCESS);
303: }
305: /*@
306: PCSetUseAmat - Sets a flag to indicate that when the preconditioner needs to apply (part of) the
307: operator during the preconditioning process it applies the Amat provided to `TSSetRHSJacobian()`,
308: `TSSetIJacobian()`, `SNESSetJacobian()`, `KSPSetOperators()` or `PCSetOperators()` not the Pmat.
310: Logically Collective
312: Input Parameters:
313: + pc - the preconditioner context
314: - flg - `PETSC_TRUE` to use the Amat, `PETSC_FALSE` to use the Pmat (default is false)
316: Options Database Key:
317: . -pc_use_amat <true,false> - use the amat argument to `KSPSetOperators()` or `PCSetOperators()` to apply the operator
319: Level: intermediate
321: Note:
322: For the common case in which the linear system matrix and the matrix used to construct the
323: preconditioner are identical, this routine has no affect.
325: .seealso: [](ch_ksp), `PC`, `PCGetUseAmat()`, `PCBJACOBI`, `PCMG`, `PCFIELDSPLIT`, `PCCOMPOSITE`,
326: `KSPSetOperators()`, `PCSetOperators()`
327: @*/
328: PetscErrorCode PCSetUseAmat(PC pc, PetscBool flg)
329: {
330: PetscFunctionBegin;
332: pc->useAmat = flg;
333: PetscFunctionReturn(PETSC_SUCCESS);
334: }
336: /*@
337: PCSetErrorIfFailure - Causes `PC` to generate an error if a floating point exception, for example a zero pivot, is detected.
339: Logically Collective
341: Input Parameters:
342: + pc - iterative context obtained from `PCCreate()`
343: - flg - `PETSC_TRUE` indicates you want the error generated
345: Level: advanced
347: Notes:
348: Normally PETSc continues if a linear solver fails due to a failed setup of a preconditioner, you can call `KSPGetConvergedReason()` after a `KSPSolve()`
349: 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
350: detected.
352: This is propagated into `KSP`s used by this `PC`, which then propagate it into `PC`s used by those `KSP`s
354: .seealso: [](ch_ksp), `PC`, `KSPSetErrorIfNotConverged()`, `PCGetInitialGuessNonzero()`, `PCSetInitialGuessKnoll()`, `PCGetInitialGuessKnoll()`
355: @*/
356: PetscErrorCode PCSetErrorIfFailure(PC pc, PetscBool flg)
357: {
358: PetscFunctionBegin;
361: pc->erroriffailure = flg;
362: PetscFunctionReturn(PETSC_SUCCESS);
363: }
365: /*@
366: PCGetUseAmat - Gets a flag to indicate that when the preconditioner needs to apply (part of) the
367: operator during the preconditioning process it applies the Amat provided to `TSSetRHSJacobian()`,
368: `TSSetIJacobian()`, `SNESSetJacobian()`, `KSPSetOperators()` or `PCSetOperators()` not the Pmat.
370: Logically Collective
372: Input Parameter:
373: . pc - the preconditioner context
375: Output Parameter:
376: . flg - `PETSC_TRUE` to use the Amat, `PETSC_FALSE` to use the Pmat (default is false)
378: Level: intermediate
380: Note:
381: For the common case in which the linear system matrix and the matrix used to construct the
382: preconditioner are identical, this routine is does nothing.
384: .seealso: [](ch_ksp), `PC`, `PCSetUseAmat()`, `PCBJACOBI`, `PCMG`, `PCFIELDSPLIT`, `PCCOMPOSITE`
385: @*/
386: PetscErrorCode PCGetUseAmat(PC pc, PetscBool *flg)
387: {
388: PetscFunctionBegin;
390: *flg = pc->useAmat;
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
394: /*@
395: PCSetKSPNestLevel - sets the amount of nesting the `KSP` that contains this `PC` has
397: Collective
399: Input Parameters:
400: + pc - the `PC`
401: - level - the nest level
403: Level: developer
405: .seealso: [](ch_ksp), `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSP`, `KSPGMRES`, `KSPType`, `KSPGetNestLevel()`, `PCGetKSPNestLevel()`, `KSPSetNestLevel()`
406: @*/
407: PetscErrorCode PCSetKSPNestLevel(PC pc, PetscInt level)
408: {
409: PetscFunctionBegin;
412: pc->kspnestlevel = level;
413: PetscFunctionReturn(PETSC_SUCCESS);
414: }
416: /*@
417: PCGetKSPNestLevel - gets the amount of nesting the `KSP` that contains this `PC` has
419: Not Collective
421: Input Parameter:
422: . pc - the `PC`
424: Output Parameter:
425: . level - the nest level
427: Level: developer
429: .seealso: [](ch_ksp), `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSP`, `KSPGMRES`, `KSPType`, `KSPSetNestLevel()`, `PCSetKSPNestLevel()`, `KSPGetNestLevel()`
430: @*/
431: PetscErrorCode PCGetKSPNestLevel(PC pc, PetscInt *level)
432: {
433: PetscFunctionBegin;
435: PetscAssertPointer(level, 2);
436: *level = pc->kspnestlevel;
437: PetscFunctionReturn(PETSC_SUCCESS);
438: }
440: /*@
441: PCCreate - Creates a preconditioner context, `PC`
443: Collective
445: Input Parameter:
446: . comm - MPI communicator
448: Output Parameter:
449: . newpc - location to put the preconditioner context
451: Level: developer
453: Note:
454: The default preconditioner for sparse matrices is `PCILU` or `PCICC` with 0 fill on one process and block Jacobi (`PCBJACOBI`) with `PCILU` or `PCICC`
455: in parallel. For dense matrices it is always `PCNONE`.
457: .seealso: [](ch_ksp), `PC`, `PCSetUp()`, `PCApply()`, `PCDestroy()`
458: @*/
459: PetscErrorCode PCCreate(MPI_Comm comm, PC *newpc)
460: {
461: PC pc;
463: PetscFunctionBegin;
464: PetscAssertPointer(newpc, 2);
465: PetscCall(PCInitializePackage());
467: PetscCall(PetscHeaderCreate(pc, PC_CLASSID, "PC", "Preconditioner", "PC", comm, PCDestroy, PCView));
468: pc->mat = NULL;
469: pc->pmat = NULL;
470: pc->setupcalled = 0;
471: pc->setfromoptionscalled = 0;
472: pc->data = NULL;
473: pc->diagonalscale = PETSC_FALSE;
474: pc->diagonalscaleleft = NULL;
475: pc->diagonalscaleright = NULL;
477: pc->modifysubmatrices = NULL;
478: pc->modifysubmatricesP = NULL;
480: *newpc = pc;
481: PetscFunctionReturn(PETSC_SUCCESS);
482: }
484: /*@
485: PCApply - Applies the preconditioner to a vector.
487: Collective
489: Input Parameters:
490: + pc - the preconditioner context
491: - x - input vector
493: Output Parameter:
494: . y - output vector
496: Level: developer
498: .seealso: [](ch_ksp), `PC`, `PCApplyTranspose()`, `PCApplyBAorAB()`
499: @*/
500: PetscErrorCode PCApply(PC pc, Vec x, Vec y)
501: {
502: PetscInt m, n, mv, nv;
504: PetscFunctionBegin;
508: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
509: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
510: /* use pmat to check vector sizes since for KSPLSQR the pmat may be of a different size than mat */
511: PetscCall(MatGetLocalSize(pc->pmat, &m, &n));
512: PetscCall(VecGetLocalSize(x, &mv));
513: PetscCall(VecGetLocalSize(y, &nv));
514: /* check pmat * y = x is feasible */
515: 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);
516: 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);
517: PetscCall(VecSetErrorIfLocked(y, 3));
519: PetscCall(PCSetUp(pc));
520: PetscCall(VecLockReadPush(x));
521: PetscCall(PetscLogEventBegin(PC_Apply, pc, x, y, 0));
522: PetscUseTypeMethod(pc, apply, x, y);
523: PetscCall(PetscLogEventEnd(PC_Apply, pc, x, y, 0));
524: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
525: PetscCall(VecLockReadPop(x));
526: PetscFunctionReturn(PETSC_SUCCESS);
527: }
529: /*@
530: PCMatApply - Applies the preconditioner to multiple vectors stored as a `MATDENSE`. Like `PCApply()`, `Y` and `X` must be different matrices.
532: Collective
534: Input Parameters:
535: + pc - the preconditioner context
536: - X - block of input vectors
538: Output Parameter:
539: . Y - block of output vectors
541: Level: developer
543: .seealso: [](ch_ksp), `PC`, `PCApply()`, `KSPMatSolve()`
544: @*/
545: PetscErrorCode PCMatApply(PC pc, Mat X, Mat Y)
546: {
547: Mat A;
548: Vec cy, cx;
549: PetscInt m1, M1, m2, M2, n1, N1, n2, N2, m3, M3, n3, N3;
550: PetscBool match;
552: PetscFunctionBegin;
556: PetscCheckSameComm(pc, 1, X, 2);
557: PetscCheckSameComm(pc, 1, Y, 3);
558: PetscCheck(Y != X, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "Y and X must be different matrices");
559: PetscCall(PCGetOperators(pc, NULL, &A));
560: PetscCall(MatGetLocalSize(A, &m3, &n3));
561: PetscCall(MatGetLocalSize(X, &m2, &n2));
562: PetscCall(MatGetLocalSize(Y, &m1, &n1));
563: PetscCall(MatGetSize(A, &M3, &N3));
564: PetscCall(MatGetSize(X, &M2, &N2));
565: PetscCall(MatGetSize(Y, &M1, &N1));
566: 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);
567: 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);
568: 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);
569: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)Y, &match, MATSEQDENSE, MATMPIDENSE, ""));
570: PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided block of output vectors not stored in a dense Mat");
571: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)X, &match, MATSEQDENSE, MATMPIDENSE, ""));
572: PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided block of input vectors not stored in a dense Mat");
573: PetscCall(PCSetUp(pc));
574: if (pc->ops->matapply) {
575: PetscCall(PetscLogEventBegin(PC_MatApply, pc, X, Y, 0));
576: PetscUseTypeMethod(pc, matapply, X, Y);
577: PetscCall(PetscLogEventEnd(PC_MatApply, pc, X, Y, 0));
578: } else {
579: PetscCall(PetscInfo(pc, "PC type %s applying column by column\n", ((PetscObject)pc)->type_name));
580: for (n1 = 0; n1 < N1; ++n1) {
581: PetscCall(MatDenseGetColumnVecRead(X, n1, &cx));
582: PetscCall(MatDenseGetColumnVecWrite(Y, n1, &cy));
583: PetscCall(PCApply(pc, cx, cy));
584: PetscCall(MatDenseRestoreColumnVecWrite(Y, n1, &cy));
585: PetscCall(MatDenseRestoreColumnVecRead(X, n1, &cx));
586: }
587: }
588: PetscFunctionReturn(PETSC_SUCCESS);
589: }
591: /*@
592: PCApplySymmetricLeft - Applies the left part of a symmetric preconditioner to a vector.
594: Collective
596: Input Parameters:
597: + pc - the preconditioner context
598: - x - input vector
600: Output Parameter:
601: . y - output vector
603: Level: developer
605: Note:
606: Currently, this routine is implemented only for `PCICC` and `PCJACOBI` preconditioners.
608: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplySymmetricRight()`
609: @*/
610: PetscErrorCode PCApplySymmetricLeft(PC pc, Vec x, Vec y)
611: {
612: PetscFunctionBegin;
616: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
617: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
618: PetscCall(PCSetUp(pc));
619: PetscCall(VecLockReadPush(x));
620: PetscCall(PetscLogEventBegin(PC_ApplySymmetricLeft, pc, x, y, 0));
621: PetscUseTypeMethod(pc, applysymmetricleft, x, y);
622: PetscCall(PetscLogEventEnd(PC_ApplySymmetricLeft, pc, x, y, 0));
623: PetscCall(VecLockReadPop(x));
624: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
625: PetscFunctionReturn(PETSC_SUCCESS);
626: }
628: /*@
629: PCApplySymmetricRight - Applies the right part of a symmetric preconditioner to a vector.
631: Collective
633: Input Parameters:
634: + pc - the preconditioner context
635: - x - input vector
637: Output Parameter:
638: . y - output vector
640: Level: developer
642: Note:
643: Currently, this routine is implemented only for `PCICC` and `PCJACOBI` preconditioners.
645: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplySymmetricLeft()`
646: @*/
647: PetscErrorCode PCApplySymmetricRight(PC pc, Vec x, Vec y)
648: {
649: PetscFunctionBegin;
653: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
654: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
655: PetscCall(PCSetUp(pc));
656: PetscCall(VecLockReadPush(x));
657: PetscCall(PetscLogEventBegin(PC_ApplySymmetricRight, pc, x, y, 0));
658: PetscUseTypeMethod(pc, applysymmetricright, x, y);
659: PetscCall(PetscLogEventEnd(PC_ApplySymmetricRight, pc, x, y, 0));
660: PetscCall(VecLockReadPop(x));
661: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
662: PetscFunctionReturn(PETSC_SUCCESS);
663: }
665: /*@
666: PCApplyTranspose - Applies the transpose of preconditioner to a vector.
668: Collective
670: Input Parameters:
671: + pc - the preconditioner context
672: - x - input vector
674: Output Parameter:
675: . y - output vector
677: Level: developer
679: Note:
680: For complex numbers this applies the non-Hermitian transpose.
682: Developer Note:
683: We need to implement a `PCApplyHermitianTranspose()`
685: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplyBAorAB()`, `PCApplyBAorABTranspose()`, `PCApplyTransposeExists()`
686: @*/
687: PetscErrorCode PCApplyTranspose(PC pc, Vec x, Vec y)
688: {
689: PetscFunctionBegin;
693: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
694: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
695: PetscCall(PCSetUp(pc));
696: PetscCall(VecLockReadPush(x));
697: PetscCall(PetscLogEventBegin(PC_Apply, pc, x, y, 0));
698: PetscUseTypeMethod(pc, applytranspose, x, y);
699: PetscCall(PetscLogEventEnd(PC_Apply, pc, x, y, 0));
700: PetscCall(VecLockReadPop(x));
701: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
702: PetscFunctionReturn(PETSC_SUCCESS);
703: }
705: /*@
706: PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation
708: Collective
710: Input Parameter:
711: . pc - the preconditioner context
713: Output Parameter:
714: . flg - `PETSC_TRUE` if a transpose operation is defined
716: Level: developer
718: .seealso: [](ch_ksp), `PC`, `PCApplyTranspose()`
719: @*/
720: PetscErrorCode PCApplyTransposeExists(PC pc, PetscBool *flg)
721: {
722: PetscFunctionBegin;
724: PetscAssertPointer(flg, 2);
725: if (pc->ops->applytranspose) *flg = PETSC_TRUE;
726: else *flg = PETSC_FALSE;
727: PetscFunctionReturn(PETSC_SUCCESS);
728: }
730: /*@
731: PCApplyBAorAB - Applies the preconditioner and operator to a vector. $y = B*A*x $ or $ y = A*B*x$.
733: Collective
735: Input Parameters:
736: + pc - the preconditioner context
737: . side - indicates the preconditioner side, one of `PC_LEFT`, `PC_RIGHT`, or `PC_SYMMETRIC`
738: . x - input vector
739: - work - work vector
741: Output Parameter:
742: . y - output vector
744: Level: developer
746: Note:
747: 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.
748: The specific `KSPSolve()` method must also be written to handle the post-solve "correction" for the diagonal scaling.
750: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplyTranspose()`, `PCApplyBAorABTranspose()`
751: @*/
752: PetscErrorCode PCApplyBAorAB(PC pc, PCSide side, Vec x, Vec y, Vec work)
753: {
754: PetscFunctionBegin;
760: PetscCheckSameComm(pc, 1, x, 3);
761: PetscCheckSameComm(pc, 1, y, 4);
762: PetscCheckSameComm(pc, 1, work, 5);
763: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
764: PetscCheck(side == PC_LEFT || side == PC_SYMMETRIC || side == PC_RIGHT, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Side must be right, left, or symmetric");
765: PetscCheck(!pc->diagonalscale || side != PC_SYMMETRIC, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot include diagonal scaling with symmetric preconditioner application");
766: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 3, PETSC_TRUE));
768: PetscCall(PCSetUp(pc));
769: if (pc->diagonalscale) {
770: if (pc->ops->applyBA) {
771: Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
772: PetscCall(VecDuplicate(x, &work2));
773: PetscCall(PCDiagonalScaleRight(pc, x, work2));
774: PetscUseTypeMethod(pc, applyBA, side, work2, y, work);
775: PetscCall(PCDiagonalScaleLeft(pc, y, y));
776: PetscCall(VecDestroy(&work2));
777: } else if (side == PC_RIGHT) {
778: PetscCall(PCDiagonalScaleRight(pc, x, y));
779: PetscCall(PCApply(pc, y, work));
780: PetscCall(MatMult(pc->mat, work, y));
781: PetscCall(PCDiagonalScaleLeft(pc, y, y));
782: } else if (side == PC_LEFT) {
783: PetscCall(PCDiagonalScaleRight(pc, x, y));
784: PetscCall(MatMult(pc->mat, y, work));
785: PetscCall(PCApply(pc, work, y));
786: PetscCall(PCDiagonalScaleLeft(pc, y, y));
787: } else PetscCheck(side != PC_SYMMETRIC, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot provide diagonal scaling with symmetric application of preconditioner");
788: } else {
789: if (pc->ops->applyBA) {
790: PetscUseTypeMethod(pc, applyBA, side, x, y, work);
791: } else if (side == PC_RIGHT) {
792: PetscCall(PCApply(pc, x, work));
793: PetscCall(MatMult(pc->mat, work, y));
794: } else if (side == PC_LEFT) {
795: PetscCall(MatMult(pc->mat, x, work));
796: PetscCall(PCApply(pc, work, y));
797: } else if (side == PC_SYMMETRIC) {
798: /* There's an extra copy here; maybe should provide 2 work vectors instead? */
799: PetscCall(PCApplySymmetricRight(pc, x, work));
800: PetscCall(MatMult(pc->mat, work, y));
801: PetscCall(VecCopy(y, work));
802: PetscCall(PCApplySymmetricLeft(pc, work, y));
803: }
804: }
805: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 4, PETSC_FALSE));
806: PetscFunctionReturn(PETSC_SUCCESS);
807: }
809: /*@
810: PCApplyBAorABTranspose - Applies the transpose of the preconditioner
811: and operator to a vector. That is, applies $B^T * A^T$ with left preconditioning,
812: NOT $(B*A)^T = A^T*B^T$.
814: Collective
816: Input Parameters:
817: + pc - the preconditioner context
818: . side - indicates the preconditioner side, one of `PC_LEFT`, `PC_RIGHT`, or `PC_SYMMETRIC`
819: . x - input vector
820: - work - work vector
822: Output Parameter:
823: . y - output vector
825: Level: developer
827: Note:
828: 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
829: defined by $B^T$. This is why this has the funny form that it computes $B^T * A^T$
831: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplyTranspose()`, `PCApplyBAorAB()`
832: @*/
833: PetscErrorCode PCApplyBAorABTranspose(PC pc, PCSide side, Vec x, Vec y, Vec work)
834: {
835: PetscFunctionBegin;
840: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
841: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 3, PETSC_TRUE));
842: if (pc->ops->applyBAtranspose) {
843: PetscUseTypeMethod(pc, applyBAtranspose, side, x, y, work);
844: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 4, PETSC_FALSE));
845: PetscFunctionReturn(PETSC_SUCCESS);
846: }
847: PetscCheck(side == PC_LEFT || side == PC_RIGHT, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Side must be right or left");
849: PetscCall(PCSetUp(pc));
850: if (side == PC_RIGHT) {
851: PetscCall(PCApplyTranspose(pc, x, work));
852: PetscCall(MatMultTranspose(pc->mat, work, y));
853: } else if (side == PC_LEFT) {
854: PetscCall(MatMultTranspose(pc->mat, x, work));
855: PetscCall(PCApplyTranspose(pc, work, y));
856: }
857: /* add support for PC_SYMMETRIC */
858: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 4, PETSC_FALSE));
859: PetscFunctionReturn(PETSC_SUCCESS);
860: }
862: /*@
863: PCApplyRichardsonExists - Determines whether a particular preconditioner has a
864: built-in fast application of Richardson's method.
866: Not Collective
868: Input Parameter:
869: . pc - the preconditioner
871: Output Parameter:
872: . exists - `PETSC_TRUE` or `PETSC_FALSE`
874: Level: developer
876: .seealso: [](ch_ksp), `PC`, `KSPRICHARDSON`, `PCApplyRichardson()`
877: @*/
878: PetscErrorCode PCApplyRichardsonExists(PC pc, PetscBool *exists)
879: {
880: PetscFunctionBegin;
882: PetscAssertPointer(exists, 2);
883: if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
884: else *exists = PETSC_FALSE;
885: PetscFunctionReturn(PETSC_SUCCESS);
886: }
888: /*@
889: PCApplyRichardson - Applies several steps of Richardson iteration with
890: the particular preconditioner. This routine is usually used by the
891: Krylov solvers and not the application code directly.
893: Collective
895: Input Parameters:
896: + pc - the preconditioner context
897: . b - the right-hand side
898: . w - one work vector
899: . rtol - relative decrease in residual norm convergence criteria
900: . abstol - absolute residual norm convergence criteria
901: . dtol - divergence residual norm increase criteria
902: . its - the number of iterations to apply.
903: - guesszero - if the input x contains nonzero initial guess
905: Output Parameters:
906: + outits - number of iterations actually used (for SOR this always equals its)
907: . reason - the reason the apply terminated
908: - y - the solution (also contains initial guess if guesszero is `PETSC_FALSE`
910: Level: developer
912: Notes:
913: Most preconditioners do not support this function. Use the command
914: `PCApplyRichardsonExists()` to determine if one does.
916: Except for the `PCMG` this routine ignores the convergence tolerances
917: and always runs for the number of iterations
919: .seealso: [](ch_ksp), `PC`, `PCApplyRichardsonExists()`
920: @*/
921: PetscErrorCode PCApplyRichardson(PC pc, Vec b, Vec y, Vec w, PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt its, PetscBool guesszero, PetscInt *outits, PCRichardsonConvergedReason *reason)
922: {
923: PetscFunctionBegin;
928: PetscCheck(b != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "b and y must be different vectors");
929: PetscCall(PCSetUp(pc));
930: PetscUseTypeMethod(pc, applyrichardson, b, y, w, rtol, abstol, dtol, its, guesszero, outits, reason);
931: PetscFunctionReturn(PETSC_SUCCESS);
932: }
934: /*@
935: PCSetFailedReason - Sets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail
937: Logically Collective
939: Input Parameters:
940: + pc - the preconditioner context
941: - reason - the reason it failedx
943: Level: advanced
945: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCFailedReason`
946: @*/
947: PetscErrorCode PCSetFailedReason(PC pc, PCFailedReason reason)
948: {
949: PetscFunctionBegin;
951: pc->failedreason = reason;
952: PetscFunctionReturn(PETSC_SUCCESS);
953: }
955: /*@
956: PCGetFailedReason - Gets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail
958: Not Collective
960: Input Parameter:
961: . pc - the preconditioner context
963: Output Parameter:
964: . reason - the reason it failed
966: Level: advanced
968: Note:
969: After call `KSPCheckDot()` or `KSPCheckNorm()` inside a `KSPSolve()` or a call to `PCReduceFailedReason()`
970: this is the maximum over reason over all ranks in the `PC` communicator and hence logically collective.
971: Otherwise it returns the local value.
973: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCSetFailedReason()`
974: @*/
975: PetscErrorCode PCGetFailedReason(PC pc, PCFailedReason *reason)
976: {
977: PetscFunctionBegin;
979: if (pc->setupcalled < 0) *reason = (PCFailedReason)pc->setupcalled;
980: else *reason = pc->failedreason;
981: PetscFunctionReturn(PETSC_SUCCESS);
982: }
984: /*@
985: PCReduceFailedReason - Reduce the failed reason among the MPI processes that share the `PC`
987: Collective
989: Input Parameter:
990: . pc - the preconditioner context
992: Level: advanced
994: Note:
995: Different MPI processes may have different reasons or no reason, see `PCGetFailedReason()`. This routine
996: makes them have a common value (failure if any MPI process had a failure).
998: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCGetFailedReason()`, `PCSetFailedReason()`
999: @*/
1000: PetscErrorCode PCReduceFailedReason(PC pc)
1001: {
1002: PetscInt buf;
1004: PetscFunctionBegin;
1006: buf = (PetscInt)pc->failedreason;
1007: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &buf, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)pc)));
1008: pc->failedreason = (PCFailedReason)buf;
1009: PetscFunctionReturn(PETSC_SUCCESS);
1010: }
1012: /*
1013: a setupcall of 0 indicates never setup,
1014: 1 indicates has been previously setup
1015: -1 indicates a PCSetUp() was attempted and failed
1016: */
1017: /*@
1018: PCSetUp - Prepares for the use of a preconditioner.
1020: Collective
1022: Input Parameter:
1023: . pc - the preconditioner context
1025: Level: developer
1027: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`
1028: @*/
1029: PetscErrorCode PCSetUp(PC pc)
1030: {
1031: const char *def;
1032: PetscObjectState matstate, matnonzerostate;
1034: PetscFunctionBegin;
1036: PetscCheck(pc->mat, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Matrix must be set first");
1038: if (pc->setupcalled && pc->reusepreconditioner) {
1039: PetscCall(PetscInfo(pc, "Leaving PC with identical preconditioner since reuse preconditioner is set\n"));
1040: PetscFunctionReturn(PETSC_SUCCESS);
1041: }
1043: PetscCall(PetscObjectStateGet((PetscObject)pc->pmat, &matstate));
1044: PetscCall(MatGetNonzeroState(pc->pmat, &matnonzerostate));
1045: if (!pc->setupcalled) {
1046: //PetscCall(PetscInfo(pc, "Setting up PC for first time\n"));
1047: pc->flag = DIFFERENT_NONZERO_PATTERN;
1048: } else if (matstate == pc->matstate) PetscFunctionReturn(PETSC_SUCCESS);
1049: else {
1050: if (matnonzerostate != pc->matnonzerostate) {
1051: PetscCall(PetscInfo(pc, "Setting up PC with different nonzero pattern\n"));
1052: pc->flag = DIFFERENT_NONZERO_PATTERN;
1053: } else {
1054: //PetscCall(PetscInfo(pc, "Setting up PC with same nonzero pattern\n"));
1055: pc->flag = SAME_NONZERO_PATTERN;
1056: }
1057: }
1058: pc->matstate = matstate;
1059: pc->matnonzerostate = matnonzerostate;
1061: if (!((PetscObject)pc)->type_name) {
1062: PetscCall(PCGetDefaultType_Private(pc, &def));
1063: PetscCall(PCSetType(pc, def));
1064: }
1066: PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
1067: PetscCall(MatSetErrorIfFailure(pc->mat, pc->erroriffailure));
1068: PetscCall(PetscLogEventBegin(PC_SetUp, pc, 0, 0, 0));
1069: if (pc->ops->setup) {
1070: PetscCall(PCLogEventsDeactivatePush());
1071: PetscUseTypeMethod(pc, setup);
1072: PetscCall(PCLogEventsDeactivatePop());
1073: }
1074: PetscCall(PetscLogEventEnd(PC_SetUp, pc, 0, 0, 0));
1075: if (!pc->setupcalled) pc->setupcalled = 1;
1076: PetscFunctionReturn(PETSC_SUCCESS);
1077: }
1079: /*@
1080: PCSetUpOnBlocks - Sets up the preconditioner for each block in
1081: the block Jacobi, overlapping Schwarz, and fieldsplit methods.
1083: Collective
1085: Input Parameter:
1086: . pc - the preconditioner context
1088: Level: developer
1090: Notes:
1091: For nested preconditioners such as `PCBJACOBI`, `PCSetUp()` is not called on each sub-`KSP` when `PCSetUp()` is
1092: called on the outer `PC`, this routine ensures it is called.
1094: It calls `PCSetUp()` if not yet called.
1096: .seealso: [](ch_ksp), `PC`, `PCSetUp()`, `PCCreate()`, `PCApply()`, `PCDestroy()`
1097: @*/
1098: PetscErrorCode PCSetUpOnBlocks(PC pc)
1099: {
1100: PetscFunctionBegin;
1102: if (!pc->setupcalled) PetscCall(PCSetUp(pc)); /* "if" to prevent -info extra prints */
1103: if (!pc->ops->setuponblocks) PetscFunctionReturn(PETSC_SUCCESS);
1104: PetscCall(PetscLogEventBegin(PC_SetUpOnBlocks, pc, 0, 0, 0));
1105: PetscCall(PCLogEventsDeactivatePush());
1106: PetscUseTypeMethod(pc, setuponblocks);
1107: PetscCall(PCLogEventsDeactivatePop());
1108: PetscCall(PetscLogEventEnd(PC_SetUpOnBlocks, pc, 0, 0, 0));
1109: PetscFunctionReturn(PETSC_SUCCESS);
1110: }
1112: /*@C
1113: PCSetModifySubMatrices - Sets a user-defined routine for modifying the
1114: submatrices that arise within certain subdomain-based preconditioners such as `PCASM`
1116: Logically Collective
1118: Input Parameters:
1119: + pc - the preconditioner context
1120: . func - routine for modifying the submatrices
1121: - ctx - optional user-defined context (may be `NULL`)
1123: Calling sequence of `func`:
1124: + pc - the preconditioner context
1125: . nsub - number of index sets
1126: . row - an array of index sets that contain the global row numbers
1127: that comprise each local submatrix
1128: . col - an array of index sets that contain the global column numbers
1129: that comprise each local submatrix
1130: . submat - array of local submatrices
1131: - ctx - optional user-defined context for private data for the
1132: user-defined func routine (may be `NULL`)
1134: Level: advanced
1136: Notes:
1137: The basic submatrices are extracted from the preconditioner matrix as
1138: usual; the user can then alter these (for example, to set different boundary
1139: conditions for each submatrix) before they are used for the local solves.
1141: `PCSetModifySubMatrices()` MUST be called before `KSPSetUp()` and
1142: `KSPSolve()`.
1144: A routine set by `PCSetModifySubMatrices()` is currently called within
1145: the block Jacobi (`PCBJACOBI`) and additive Schwarz (`PCASM`)
1146: preconditioners. All other preconditioners ignore this routine.
1148: .seealso: [](ch_ksp), `PC`, `PCBJACOBI`, `PCASM`, `PCModifySubMatrices()`
1149: @*/
1150: PetscErrorCode PCSetModifySubMatrices(PC pc, PetscErrorCode (*func)(PC pc, PetscInt nsub, const IS row[], const IS col[], Mat submat[], void *ctx), void *ctx)
1151: {
1152: PetscFunctionBegin;
1154: pc->modifysubmatrices = func;
1155: pc->modifysubmatricesP = ctx;
1156: PetscFunctionReturn(PETSC_SUCCESS);
1157: }
1159: /*@C
1160: PCModifySubMatrices - Calls an optional user-defined routine within
1161: certain preconditioners if one has been set with `PCSetModifySubMatrices()`.
1163: Collective
1165: Input Parameters:
1166: + pc - the preconditioner context
1167: . nsub - the number of local submatrices
1168: . row - an array of index sets that contain the global row numbers
1169: that comprise each local submatrix
1170: . col - an array of index sets that contain the global column numbers
1171: that comprise each local submatrix
1172: . submat - array of local submatrices
1173: - ctx - optional user-defined context for private data for the
1174: user-defined routine (may be `NULL`)
1176: Output Parameter:
1177: . submat - array of local submatrices (the entries of which may
1178: have been modified)
1180: Level: developer
1182: Note:
1183: The user should NOT generally call this routine, as it will
1184: automatically be called within certain preconditioners.
1186: .seealso: [](ch_ksp), `PC`, `PCSetModifySubMatrices()`
1187: @*/
1188: PetscErrorCode PCModifySubMatrices(PC pc, PetscInt nsub, const IS row[], const IS col[], Mat submat[], void *ctx)
1189: {
1190: PetscFunctionBegin;
1192: if (!pc->modifysubmatrices) PetscFunctionReturn(PETSC_SUCCESS);
1193: PetscCall(PetscLogEventBegin(PC_ModifySubMatrices, pc, 0, 0, 0));
1194: PetscCall((*pc->modifysubmatrices)(pc, nsub, row, col, submat, ctx));
1195: PetscCall(PetscLogEventEnd(PC_ModifySubMatrices, pc, 0, 0, 0));
1196: PetscFunctionReturn(PETSC_SUCCESS);
1197: }
1199: /*@
1200: PCSetOperators - Sets the matrix associated with the linear system and
1201: a (possibly) different one associated with the preconditioner.
1203: Logically Collective
1205: Input Parameters:
1206: + pc - the preconditioner context
1207: . Amat - the matrix that defines the linear system
1208: - Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
1210: Level: intermediate
1212: Notes:
1213: Passing a `NULL` for `Amat` or `Pmat` removes the matrix that is currently used.
1215: If you wish to replace either `Amat` or `Pmat` but leave the other one untouched then
1216: first call `KSPGetOperators()` to get the one you wish to keep, call `PetscObjectReference()`
1217: on it and then pass it back in in your call to `KSPSetOperators()`.
1219: More Notes about Repeated Solution of Linear Systems:
1220: PETSc does NOT reset the matrix entries of either `Amat` or `Pmat`
1221: to zero after a linear solve; the user is completely responsible for
1222: matrix assembly. See the routine `MatZeroEntries()` if desiring to
1223: zero all elements of a matrix.
1225: .seealso: [](ch_ksp), `PC`, `PCGetOperators()`, `MatZeroEntries()`
1226: @*/
1227: PetscErrorCode PCSetOperators(PC pc, Mat Amat, Mat Pmat)
1228: {
1229: PetscInt m1, n1, m2, n2;
1231: PetscFunctionBegin;
1235: if (Amat) PetscCheckSameComm(pc, 1, Amat, 2);
1236: if (Pmat) PetscCheckSameComm(pc, 1, Pmat, 3);
1237: if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
1238: PetscCall(MatGetLocalSize(Amat, &m1, &n1));
1239: PetscCall(MatGetLocalSize(pc->mat, &m2, &n2));
1240: 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);
1241: PetscCall(MatGetLocalSize(Pmat, &m1, &n1));
1242: PetscCall(MatGetLocalSize(pc->pmat, &m2, &n2));
1243: 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);
1244: }
1246: if (Pmat != pc->pmat) {
1247: /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
1248: pc->matnonzerostate = -1;
1249: pc->matstate = -1;
1250: }
1252: /* reference first in case the matrices are the same */
1253: if (Amat) PetscCall(PetscObjectReference((PetscObject)Amat));
1254: PetscCall(MatDestroy(&pc->mat));
1255: if (Pmat) PetscCall(PetscObjectReference((PetscObject)Pmat));
1256: PetscCall(MatDestroy(&pc->pmat));
1257: pc->mat = Amat;
1258: pc->pmat = Pmat;
1259: PetscFunctionReturn(PETSC_SUCCESS);
1260: }
1262: /*@
1263: PCSetReusePreconditioner - reuse the current preconditioner even if the operator in the preconditioner has changed.
1265: Logically Collective
1267: Input Parameters:
1268: + pc - the preconditioner context
1269: - flag - `PETSC_TRUE` do not compute a new preconditioner, `PETSC_FALSE` do compute a new preconditioner
1271: Level: intermediate
1273: Note:
1274: Normally if a matrix inside a `PC` changes the `PC` automatically updates itself using information from the changed matrix. This option
1275: prevents this.
1277: .seealso: [](ch_ksp), `PC`, `PCGetOperators()`, `MatZeroEntries()`, `PCGetReusePreconditioner()`, `KSPSetReusePreconditioner()`
1278: @*/
1279: PetscErrorCode PCSetReusePreconditioner(PC pc, PetscBool flag)
1280: {
1281: PetscFunctionBegin;
1284: pc->reusepreconditioner = flag;
1285: PetscFunctionReturn(PETSC_SUCCESS);
1286: }
1288: /*@
1289: PCGetReusePreconditioner - Determines if the `PC` reuses the current preconditioner even if the operator in the preconditioner has changed.
1291: Not Collective
1293: Input Parameter:
1294: . pc - the preconditioner context
1296: Output Parameter:
1297: . flag - `PETSC_TRUE` do not compute a new preconditioner, `PETSC_FALSE` do compute a new preconditioner
1299: Level: intermediate
1301: .seealso: [](ch_ksp), `PC`, `PCGetOperators()`, `MatZeroEntries()`, `PCSetReusePreconditioner()`
1302: @*/
1303: PetscErrorCode PCGetReusePreconditioner(PC pc, PetscBool *flag)
1304: {
1305: PetscFunctionBegin;
1307: PetscAssertPointer(flag, 2);
1308: *flag = pc->reusepreconditioner;
1309: PetscFunctionReturn(PETSC_SUCCESS);
1310: }
1312: /*@
1313: PCGetOperators - Gets the matrix associated with the linear system and
1314: possibly a different one associated with the preconditioner.
1316: Not Collective, though parallel `Mat`s are returned if `pc` is parallel
1318: Input Parameter:
1319: . pc - the preconditioner context
1321: Output Parameters:
1322: + Amat - the matrix defining the linear system
1323: - Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.
1325: Level: intermediate
1327: Note:
1328: Does not increase the reference count of the matrices, so you should not destroy them
1330: Alternative usage: If the operators have NOT been set with `KSPSetOperators()`/`PCSetOperators()` then the operators
1331: are created in `PC` and returned to the user. In this case, if both operators
1332: mat and pmat are requested, two DIFFERENT operators will be returned. If
1333: only one is requested both operators in the PC will be the same (i.e. as
1334: if one had called `KSPSetOperators()`/`PCSetOperators()` with the same argument for both Mats).
1335: The user must set the sizes of the returned matrices and their type etc just
1336: as if the user created them with `MatCreate()`. For example,
1338: .vb
1339: KSP/PCGetOperators(ksp/pc,&Amat,NULL); is equivalent to
1340: set size, type, etc of Amat
1342: MatCreate(comm,&mat);
1343: KSP/PCSetOperators(ksp/pc,Amat,Amat);
1344: PetscObjectDereference((PetscObject)mat);
1345: set size, type, etc of Amat
1346: .ve
1348: and
1350: .vb
1351: KSP/PCGetOperators(ksp/pc,&Amat,&Pmat); is equivalent to
1352: set size, type, etc of Amat and Pmat
1354: MatCreate(comm,&Amat);
1355: MatCreate(comm,&Pmat);
1356: KSP/PCSetOperators(ksp/pc,Amat,Pmat);
1357: PetscObjectDereference((PetscObject)Amat);
1358: PetscObjectDereference((PetscObject)Pmat);
1359: set size, type, etc of Amat and Pmat
1360: .ve
1362: The rationale for this support is so that when creating a `TS`, `SNES`, or `KSP` the hierarchy
1363: of underlying objects (i.e. `SNES`, `KSP`, `PC`, `Mat`) and their lifespans can be completely
1364: managed by the top most level object (i.e. the `TS`, `SNES`, or `KSP`). Another way to look
1365: at this is when you create a `SNES` you do not NEED to create a `KSP` and attach it to
1366: the `SNES` object (the `SNES` object manages it for you). Similarly when you create a KSP
1367: you do not need to attach a `PC` to it (the `KSP` object manages the `PC` object for you).
1368: Thus, why should YOU have to create the `Mat` and attach it to the `SNES`/`KSP`/`PC`, when
1369: it can be created for you?
1371: .seealso: [](ch_ksp), `PC`, `PCSetOperators()`, `KSPGetOperators()`, `KSPSetOperators()`, `PCGetOperatorsSet()`
1372: @*/
1373: PetscErrorCode PCGetOperators(PC pc, Mat *Amat, Mat *Pmat)
1374: {
1375: PetscFunctionBegin;
1377: if (Amat) {
1378: if (!pc->mat) {
1379: if (pc->pmat && !Pmat) { /* Pmat has been set, but user did not request it, so use for Amat */
1380: pc->mat = pc->pmat;
1381: PetscCall(PetscObjectReference((PetscObject)pc->mat));
1382: } else { /* both Amat and Pmat are empty */
1383: PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &pc->mat));
1384: if (!Pmat) { /* user did NOT request Pmat, so make same as Amat */
1385: pc->pmat = pc->mat;
1386: PetscCall(PetscObjectReference((PetscObject)pc->pmat));
1387: }
1388: }
1389: }
1390: *Amat = pc->mat;
1391: }
1392: if (Pmat) {
1393: if (!pc->pmat) {
1394: if (pc->mat && !Amat) { /* Amat has been set but was not requested, so use for pmat */
1395: pc->pmat = pc->mat;
1396: PetscCall(PetscObjectReference((PetscObject)pc->pmat));
1397: } else {
1398: PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &pc->pmat));
1399: if (!Amat) { /* user did NOT request Amat, so make same as Pmat */
1400: pc->mat = pc->pmat;
1401: PetscCall(PetscObjectReference((PetscObject)pc->mat));
1402: }
1403: }
1404: }
1405: *Pmat = pc->pmat;
1406: }
1407: PetscFunctionReturn(PETSC_SUCCESS);
1408: }
1410: /*@
1411: PCGetOperatorsSet - Determines if the matrix associated with the linear system and
1412: possibly a different one associated with the preconditioner have been set in the `PC`.
1414: Not Collective, though the results on all processes should be the same
1416: Input Parameter:
1417: . pc - the preconditioner context
1419: Output Parameters:
1420: + mat - the matrix associated with the linear system was set
1421: - pmat - matrix associated with the preconditioner was set, usually the same
1423: Level: intermediate
1425: .seealso: [](ch_ksp), `PC`, `PCSetOperators()`, `KSPGetOperators()`, `KSPSetOperators()`, `PCGetOperators()`
1426: @*/
1427: PetscErrorCode PCGetOperatorsSet(PC pc, PetscBool *mat, PetscBool *pmat)
1428: {
1429: PetscFunctionBegin;
1431: if (mat) *mat = (pc->mat) ? PETSC_TRUE : PETSC_FALSE;
1432: if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1433: PetscFunctionReturn(PETSC_SUCCESS);
1434: }
1436: /*@
1437: PCFactorGetMatrix - Gets the factored matrix from the
1438: preconditioner context. This routine is valid only for the `PCLU`,
1439: `PCILU`, `PCCHOLESKY`, and `PCICC` methods.
1441: Not Collective though `mat` is parallel if `pc` is parallel
1443: Input Parameter:
1444: . pc - the preconditioner context
1446: Output Parameters:
1447: . mat - the factored matrix
1449: Level: advanced
1451: Note:
1452: Does not increase the reference count for `mat` so DO NOT destroy it
1454: .seealso: [](ch_ksp), `PC`, `PCLU`, `PCILU`, `PCCHOLESKY`, `PCICC`
1455: @*/
1456: PetscErrorCode PCFactorGetMatrix(PC pc, Mat *mat)
1457: {
1458: PetscFunctionBegin;
1460: PetscAssertPointer(mat, 2);
1461: PetscCall(PCFactorSetUpMatSolverType(pc));
1462: PetscUseTypeMethod(pc, getfactoredmatrix, mat);
1463: PetscFunctionReturn(PETSC_SUCCESS);
1464: }
1466: /*@
1467: PCSetOptionsPrefix - Sets the prefix used for searching for all
1468: `PC` options in the database.
1470: Logically Collective
1472: Input Parameters:
1473: + pc - the preconditioner context
1474: - prefix - the prefix string to prepend to all `PC` option requests
1476: Note:
1477: A hyphen (-) must NOT be given at the beginning of the prefix name.
1478: The first character of all runtime options is AUTOMATICALLY the
1479: hyphen.
1481: Level: advanced
1483: .seealso: [](ch_ksp), `PC`, `PCSetFromOptions`, `PCAppendOptionsPrefix()`, `PCGetOptionsPrefix()`
1484: @*/
1485: PetscErrorCode PCSetOptionsPrefix(PC pc, const char prefix[])
1486: {
1487: PetscFunctionBegin;
1489: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc, prefix));
1490: PetscFunctionReturn(PETSC_SUCCESS);
1491: }
1493: /*@
1494: PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1495: `PC` options in the database.
1497: Logically Collective
1499: Input Parameters:
1500: + pc - the preconditioner context
1501: - prefix - the prefix string to prepend to all `PC` option requests
1503: Note:
1504: A hyphen (-) must NOT be given at the beginning of the prefix name.
1505: The first character of all runtime options is AUTOMATICALLY the
1506: hyphen.
1508: Level: advanced
1510: .seealso: [](ch_ksp), `PC`, `PCSetFromOptions`, `PCSetOptionsPrefix()`, `PCGetOptionsPrefix()`
1511: @*/
1512: PetscErrorCode PCAppendOptionsPrefix(PC pc, const char prefix[])
1513: {
1514: PetscFunctionBegin;
1516: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)pc, prefix));
1517: PetscFunctionReturn(PETSC_SUCCESS);
1518: }
1520: /*@
1521: PCGetOptionsPrefix - Gets the prefix used for searching for all
1522: PC options in the database.
1524: Not Collective
1526: Input Parameter:
1527: . pc - the preconditioner context
1529: Output Parameter:
1530: . prefix - pointer to the prefix string used, is returned
1532: Level: advanced
1534: Fortran Note:
1535: The user should pass in a string `prefix` of
1536: sufficient length to hold the prefix.
1538: .seealso: [](ch_ksp), `PC`, `PCSetFromOptions`, `PCSetOptionsPrefix()`, `PCAppendOptionsPrefix()`
1539: @*/
1540: PetscErrorCode PCGetOptionsPrefix(PC pc, const char *prefix[])
1541: {
1542: PetscFunctionBegin;
1544: PetscAssertPointer(prefix, 2);
1545: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, prefix));
1546: PetscFunctionReturn(PETSC_SUCCESS);
1547: }
1549: /*
1550: Indicates the right-hand side will be changed by KSPSolve(), this occurs for a few
1551: preconditioners including BDDC and Eisentat that transform the equations before applying
1552: the Krylov methods
1553: */
1554: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC pc, PetscBool *change)
1555: {
1556: PetscFunctionBegin;
1558: PetscAssertPointer(change, 2);
1559: *change = PETSC_FALSE;
1560: PetscTryMethod(pc, "PCPreSolveChangeRHS_C", (PC, PetscBool *), (pc, change));
1561: PetscFunctionReturn(PETSC_SUCCESS);
1562: }
1564: /*@
1565: PCPreSolve - Optional pre-solve phase, intended for any
1566: preconditioner-specific actions that must be performed before
1567: the iterative solve itself.
1569: Collective
1571: Input Parameters:
1572: + pc - the preconditioner context
1573: - ksp - the Krylov subspace context
1575: Level: developer
1577: Example Usage:
1578: .vb
1579: PCPreSolve(pc,ksp);
1580: KSPSolve(ksp,b,x);
1581: PCPostSolve(pc,ksp);
1582: .ve
1584: Notes:
1585: The pre-solve phase is distinct from the `PCSetUp()` phase.
1587: `KSPSolve()` calls this directly, so is rarely called by the user.
1589: .seealso: [](ch_ksp), `PC`, `PCPostSolve()`
1590: @*/
1591: PetscErrorCode PCPreSolve(PC pc, KSP ksp)
1592: {
1593: Vec x, rhs;
1595: PetscFunctionBegin;
1598: pc->presolvedone++;
1599: PetscCheck(pc->presolvedone <= 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot embed PCPreSolve() more than twice");
1600: PetscCall(KSPGetSolution(ksp, &x));
1601: PetscCall(KSPGetRhs(ksp, &rhs));
1603: if (pc->ops->presolve) PetscUseTypeMethod(pc, presolve, ksp, rhs, x);
1604: else if (pc->presolve) PetscCall(pc->presolve(pc, ksp));
1605: PetscFunctionReturn(PETSC_SUCCESS);
1606: }
1608: /*@C
1609: PCSetPreSolve - Sets function used by `PCPreSolve()` which is intended for any
1610: preconditioner-specific actions that must be performed before
1611: the iterative solve itself.
1613: Logically Collective
1615: Input Parameters:
1616: + pc - the preconditioner object
1617: - presolve - the function to call before the solve
1619: Calling sequence of `presolve`:
1620: + pc - the `PC` context
1621: - ksp - the `KSP` context
1623: Level: developer
1625: .seealso: [](ch_ksp), `PC`, `PCSetUp()`, `PCPreSolve()`
1626: @*/
1627: PetscErrorCode PCSetPreSolve(PC pc, PetscErrorCode (*presolve)(PC pc, KSP ksp))
1628: {
1629: PetscFunctionBegin;
1631: pc->presolve = presolve;
1632: PetscFunctionReturn(PETSC_SUCCESS);
1633: }
1635: /*@
1636: PCPostSolve - Optional post-solve phase, intended for any
1637: preconditioner-specific actions that must be performed after
1638: the iterative solve itself.
1640: Collective
1642: Input Parameters:
1643: + pc - the preconditioner context
1644: - ksp - the Krylov subspace context
1646: Example Usage:
1647: .vb
1648: PCPreSolve(pc,ksp);
1649: KSPSolve(ksp,b,x);
1650: PCPostSolve(pc,ksp);
1651: .ve
1653: Level: developer
1655: Note:
1656: `KSPSolve()` calls this routine directly, so it is rarely called by the user.
1658: .seealso: [](ch_ksp), `PC`, `PCSetPostSolve()`, `PCSetPresolve()`, `PCPreSolve()`, `KSPSolve()`
1659: @*/
1660: PetscErrorCode PCPostSolve(PC pc, KSP ksp)
1661: {
1662: Vec x, rhs;
1664: PetscFunctionBegin;
1667: pc->presolvedone--;
1668: PetscCall(KSPGetSolution(ksp, &x));
1669: PetscCall(KSPGetRhs(ksp, &rhs));
1670: PetscTryTypeMethod(pc, postsolve, ksp, rhs, x);
1671: PetscFunctionReturn(PETSC_SUCCESS);
1672: }
1674: /*@
1675: PCLoad - Loads a `PC` that has been stored in binary with `PCView()`.
1677: Collective
1679: Input Parameters:
1680: + newdm - the newly loaded `PC`, this needs to have been created with `PCCreate()` or
1681: some related function before a call to `PCLoad()`.
1682: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()`
1684: Level: intermediate
1686: Note:
1687: The type is determined by the data in the file, any `PCType` set into the `PC` before this call is ignored.
1689: .seealso: [](ch_ksp), `PC`, `PetscViewerBinaryOpen()`, `PCView()`, `MatLoad()`, `VecLoad()`
1690: @*/
1691: PetscErrorCode PCLoad(PC newdm, PetscViewer viewer)
1692: {
1693: PetscBool isbinary;
1694: PetscInt classid;
1695: char type[256];
1697: PetscFunctionBegin;
1700: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1701: PetscCheck(isbinary, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1703: PetscCall(PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT));
1704: PetscCheck(classid == PC_FILE_CLASSID, PetscObjectComm((PetscObject)newdm), PETSC_ERR_ARG_WRONG, "Not PC next in file");
1705: PetscCall(PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR));
1706: PetscCall(PCSetType(newdm, type));
1707: PetscTryTypeMethod(newdm, load, viewer);
1708: PetscFunctionReturn(PETSC_SUCCESS);
1709: }
1711: #include <petscdraw.h>
1712: #if defined(PETSC_HAVE_SAWS)
1713: #include <petscviewersaws.h>
1714: #endif
1716: /*@
1717: PCViewFromOptions - View from the `PC` based on options in the options database
1719: Collective
1721: Input Parameters:
1722: + A - the `PC` context
1723: . obj - Optional object that provides the options prefix
1724: - name - command line option
1726: Level: intermediate
1728: .seealso: [](ch_ksp), `PC`, `PCView`, `PetscObjectViewFromOptions()`, `PCCreate()`
1729: @*/
1730: PetscErrorCode PCViewFromOptions(PC A, PetscObject obj, const char name[])
1731: {
1732: PetscFunctionBegin;
1734: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
1735: PetscFunctionReturn(PETSC_SUCCESS);
1736: }
1738: /*@
1739: PCView - Prints information about the `PC`
1741: Collective
1743: Input Parameters:
1744: + pc - the `PC` context
1745: - viewer - optional visualization context
1747: Level: developer
1749: Notes:
1750: The available visualization contexts include
1751: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
1752: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
1753: output where only the first processor opens
1754: the file. All other processors send their
1755: data to the first processor to print.
1757: The user can open an alternative visualization contexts with
1758: `PetscViewerASCIIOpen()` (output to a specified file).
1760: .seealso: [](ch_ksp), `PC`, `PetscViewer`, `KSPView()`, `PetscViewerASCIIOpen()`
1761: @*/
1762: PetscErrorCode PCView(PC pc, PetscViewer viewer)
1763: {
1764: PCType cstr;
1765: PetscViewerFormat format;
1766: PetscBool iascii, isstring, isbinary, isdraw, pop = PETSC_FALSE;
1767: #if defined(PETSC_HAVE_SAWS)
1768: PetscBool issaws;
1769: #endif
1771: PetscFunctionBegin;
1773: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc), &viewer));
1775: PetscCheckSameComm(pc, 1, viewer, 2);
1777: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1778: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
1779: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1780: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1781: #if defined(PETSC_HAVE_SAWS)
1782: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws));
1783: #endif
1785: if (iascii) {
1786: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)pc, viewer));
1787: if (!pc->setupcalled) PetscCall(PetscViewerASCIIPrintf(viewer, " PC has not been set up so information may be incomplete\n"));
1788: PetscCall(PetscViewerASCIIPushTab(viewer));
1789: PetscTryTypeMethod(pc, view, viewer);
1790: PetscCall(PetscViewerASCIIPopTab(viewer));
1791: if (pc->mat) {
1792: PetscCall(PetscViewerGetFormat(viewer, &format));
1793: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
1794: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
1795: pop = PETSC_TRUE;
1796: }
1797: if (pc->pmat == pc->mat) {
1798: PetscCall(PetscViewerASCIIPrintf(viewer, " linear system matrix = precond matrix:\n"));
1799: PetscCall(PetscViewerASCIIPushTab(viewer));
1800: PetscCall(MatView(pc->mat, viewer));
1801: PetscCall(PetscViewerASCIIPopTab(viewer));
1802: } else {
1803: if (pc->pmat) {
1804: PetscCall(PetscViewerASCIIPrintf(viewer, " linear system matrix followed by preconditioner matrix:\n"));
1805: } else {
1806: PetscCall(PetscViewerASCIIPrintf(viewer, " linear system matrix:\n"));
1807: }
1808: PetscCall(PetscViewerASCIIPushTab(viewer));
1809: PetscCall(MatView(pc->mat, viewer));
1810: if (pc->pmat) PetscCall(MatView(pc->pmat, viewer));
1811: PetscCall(PetscViewerASCIIPopTab(viewer));
1812: }
1813: if (pop) PetscCall(PetscViewerPopFormat(viewer));
1814: }
1815: } else if (isstring) {
1816: PetscCall(PCGetType(pc, &cstr));
1817: PetscCall(PetscViewerStringSPrintf(viewer, " PCType: %-7.7s", cstr));
1818: PetscTryTypeMethod(pc, view, viewer);
1819: if (pc->mat) PetscCall(MatView(pc->mat, viewer));
1820: if (pc->pmat && pc->pmat != pc->mat) PetscCall(MatView(pc->pmat, viewer));
1821: } else if (isbinary) {
1822: PetscInt classid = PC_FILE_CLASSID;
1823: MPI_Comm comm;
1824: PetscMPIInt rank;
1825: char type[256];
1827: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
1828: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1829: if (rank == 0) {
1830: PetscCall(PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT));
1831: PetscCall(PetscStrncpy(type, ((PetscObject)pc)->type_name, 256));
1832: PetscCall(PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR));
1833: }
1834: PetscTryTypeMethod(pc, view, viewer);
1835: } else if (isdraw) {
1836: PetscDraw draw;
1837: char str[25];
1838: PetscReal x, y, bottom, h;
1839: PetscInt n;
1841: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1842: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
1843: if (pc->mat) {
1844: PetscCall(MatGetSize(pc->mat, &n, NULL));
1845: PetscCall(PetscSNPrintf(str, 25, "PC: %s (%" PetscInt_FMT ")", ((PetscObject)pc)->type_name, n));
1846: } else {
1847: PetscCall(PetscSNPrintf(str, 25, "PC: %s", ((PetscObject)pc)->type_name));
1848: }
1849: PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h));
1850: bottom = y - h;
1851: PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
1852: PetscTryTypeMethod(pc, view, viewer);
1853: PetscCall(PetscDrawPopCurrentPoint(draw));
1854: #if defined(PETSC_HAVE_SAWS)
1855: } else if (issaws) {
1856: PetscMPIInt rank;
1858: PetscCall(PetscObjectName((PetscObject)pc));
1859: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1860: if (!((PetscObject)pc)->amsmem && rank == 0) PetscCall(PetscObjectViewSAWs((PetscObject)pc, viewer));
1861: if (pc->mat) PetscCall(MatView(pc->mat, viewer));
1862: if (pc->pmat && pc->pmat != pc->mat) PetscCall(MatView(pc->pmat, viewer));
1863: #endif
1864: }
1865: PetscFunctionReturn(PETSC_SUCCESS);
1866: }
1868: /*@C
1869: PCRegister - Adds a method (`PCType`) to the preconditioner package.
1871: Not collective. No Fortran Support
1873: Input Parameters:
1874: + sname - name of a new user-defined solver
1875: - function - routine to create method context
1877: Example Usage:
1878: .vb
1879: PCRegister("my_solver", MySolverCreate);
1880: .ve
1882: Then, your solver can be chosen with the procedural interface via
1883: $ PCSetType(pc, "my_solver")
1884: or at runtime via the option
1885: $ -pc_type my_solver
1887: Level: advanced
1889: Note:
1890: `PCRegister()` may be called multiple times to add several user-defined preconditioners.
1892: .seealso: [](ch_ksp), `PC`, `PCType`, `PCRegisterAll()`
1893: @*/
1894: PetscErrorCode PCRegister(const char sname[], PetscErrorCode (*function)(PC))
1895: {
1896: PetscFunctionBegin;
1897: PetscCall(PCInitializePackage());
1898: PetscCall(PetscFunctionListAdd(&PCList, sname, function));
1899: PetscFunctionReturn(PETSC_SUCCESS);
1900: }
1902: static PetscErrorCode MatMult_PC(Mat A, Vec X, Vec Y)
1903: {
1904: PC pc;
1906: PetscFunctionBegin;
1907: PetscCall(MatShellGetContext(A, &pc));
1908: PetscCall(PCApply(pc, X, Y));
1909: PetscFunctionReturn(PETSC_SUCCESS);
1910: }
1912: /*@
1913: PCComputeOperator - Computes the explicit preconditioned operator.
1915: Collective
1917: Input Parameters:
1918: + pc - the preconditioner object
1919: - mattype - the `MatType` to be used for the operator
1921: Output Parameter:
1922: . mat - the explicit preconditioned operator
1924: Level: advanced
1926: Note:
1927: This computation is done by applying the operators to columns of the identity matrix.
1928: This routine is costly in general, and is recommended for use only with relatively small systems.
1929: Currently, this routine uses a dense matrix format when `mattype` == `NULL`
1931: .seealso: [](ch_ksp), `PC`, `KSPComputeOperator()`, `MatType`
1932: @*/
1933: PetscErrorCode PCComputeOperator(PC pc, MatType mattype, Mat *mat)
1934: {
1935: PetscInt N, M, m, n;
1936: Mat A, Apc;
1938: PetscFunctionBegin;
1940: PetscAssertPointer(mat, 3);
1941: PetscCall(PCGetOperators(pc, &A, NULL));
1942: PetscCall(MatGetLocalSize(A, &m, &n));
1943: PetscCall(MatGetSize(A, &M, &N));
1944: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)pc), m, n, M, N, pc, &Apc));
1945: PetscCall(MatShellSetOperation(Apc, MATOP_MULT, (void (*)(void))MatMult_PC));
1946: PetscCall(MatComputeOperator(Apc, mattype, mat));
1947: PetscCall(MatDestroy(&Apc));
1948: PetscFunctionReturn(PETSC_SUCCESS);
1949: }
1951: /*@
1952: PCSetCoordinates - sets the coordinates of all the nodes on the local process
1954: Collective
1956: Input Parameters:
1957: + pc - the solver context
1958: . dim - the dimension of the coordinates 1, 2, or 3
1959: . nloc - the blocked size of the coordinates array
1960: - coords - the coordinates array
1962: Level: intermediate
1964: Note:
1965: `coords` is an array of the dim coordinates for the nodes on
1966: the local processor, of size `dim`*`nloc`.
1967: If there are 108 equation on a processor
1968: for a displacement finite element discretization of elasticity (so
1969: that there are nloc = 36 = 108/3 nodes) then the array must have 108
1970: double precision values (ie, 3 * 36). These x y z coordinates
1971: should be ordered for nodes 0 to N-1 like so: [ 0.x, 0.y, 0.z, 1.x,
1972: ... , N-1.z ].
1974: .seealso: [](ch_ksp), `PC`, `MatSetNearNullSpace()`
1975: @*/
1976: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[])
1977: {
1978: PetscFunctionBegin;
1981: PetscTryMethod(pc, "PCSetCoordinates_C", (PC, PetscInt, PetscInt, PetscReal[]), (pc, dim, nloc, coords));
1982: PetscFunctionReturn(PETSC_SUCCESS);
1983: }
1985: /*@
1986: PCGetInterpolations - Gets interpolation matrices for all levels (except level 0)
1988: Logically Collective
1990: Input Parameter:
1991: . pc - the precondition context
1993: Output Parameters:
1994: + num_levels - the number of levels
1995: - interpolations - the interpolation matrices (size of `num_levels`-1)
1997: Level: advanced
1999: Developer Note:
2000: Why is this here instead of in `PCMG` etc?
2002: .seealso: [](ch_ksp), `PC`, `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetInterpolation()`, `PCGetCoarseOperators()`
2003: @*/
2004: PetscErrorCode PCGetInterpolations(PC pc, PetscInt *num_levels, Mat *interpolations[])
2005: {
2006: PetscFunctionBegin;
2008: PetscAssertPointer(num_levels, 2);
2009: PetscAssertPointer(interpolations, 3);
2010: PetscUseMethod(pc, "PCGetInterpolations_C", (PC, PetscInt *, Mat *[]), (pc, num_levels, interpolations));
2011: PetscFunctionReturn(PETSC_SUCCESS);
2012: }
2014: /*@
2015: PCGetCoarseOperators - Gets coarse operator matrices for all levels (except the finest level)
2017: Logically Collective
2019: Input Parameter:
2020: . pc - the precondition context
2022: Output Parameters:
2023: + num_levels - the number of levels
2024: - coarseOperators - the coarse operator matrices (size of `num_levels`-1)
2026: Level: advanced
2028: Developer Note:
2029: Why is this here instead of in `PCMG` etc?
2031: .seealso: [](ch_ksp), `PC`, `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetRScale()`, `PCMGGetInterpolation()`, `PCGetInterpolations()`
2032: @*/
2033: PetscErrorCode PCGetCoarseOperators(PC pc, PetscInt *num_levels, Mat *coarseOperators[])
2034: {
2035: PetscFunctionBegin;
2037: PetscAssertPointer(num_levels, 2);
2038: PetscAssertPointer(coarseOperators, 3);
2039: PetscUseMethod(pc, "PCGetCoarseOperators_C", (PC, PetscInt *, Mat *[]), (pc, num_levels, coarseOperators));
2040: PetscFunctionReturn(PETSC_SUCCESS);
2041: }