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: PetscErrorCode PCGetDefaultType_Private(PC pc, const char *type[])
15: {
16: PetscMPIInt size;
17: PetscBool hasop, 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, &hasop));
23: if (size == 1) {
24: PetscCall(MatGetFactorAvailable(pc->pmat, "petsc", MAT_FACTOR_ICC, &flg1));
25: PetscCall(MatGetFactorAvailable(pc->pmat, "petsc", MAT_FACTOR_ILU, &flg2));
26: PetscCall(MatIsSymmetricKnown(pc->pmat, &set, &flg3));
27: PetscCall(PetscObjectTypeCompareAny((PetscObject)pc->pmat, &isnormal, MATNORMAL, MATNORMALHERMITIAN, NULL));
28: if (flg1 && (!flg2 || (set && flg3))) {
29: *type = PCICC;
30: } else if (flg2) {
31: *type = PCILU;
32: } else if (isnormal) {
33: *type = PCNONE;
34: } else if (hasop) { /* likely is a parallel matrix run on one processor */
35: *type = PCBJACOBI;
36: } else {
37: *type = PCNONE;
38: }
39: } else {
40: if (hasop) {
41: *type = PCBJACOBI;
42: } else {
43: *type = PCNONE;
44: }
45: }
46: } else {
47: if (size == 1) {
48: *type = PCILU;
49: } else {
50: *type = PCBJACOBI;
51: }
52: }
53: PetscFunctionReturn(PETSC_SUCCESS);
54: }
56: /*@
57: PCReset - Resets a `PC` context to the pcsetupcalled = 0 state and removes any allocated `Vec`s and `Mat`s
59: Collective
61: Input Parameter:
62: . pc - the preconditioner context
64: Level: developer
66: Note:
67: This allows a `PC` to be reused for a different sized linear system but using the same options that have been previously set in `pc`
69: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetUp()`
70: @*/
71: PetscErrorCode PCReset(PC pc)
72: {
73: PetscFunctionBegin;
75: PetscTryTypeMethod(pc, reset);
76: PetscCall(VecDestroy(&pc->diagonalscaleright));
77: PetscCall(VecDestroy(&pc->diagonalscaleleft));
78: PetscCall(MatDestroy(&pc->pmat));
79: PetscCall(MatDestroy(&pc->mat));
81: pc->setupcalled = 0;
82: PetscFunctionReturn(PETSC_SUCCESS);
83: }
85: /*@C
86: PCDestroy - Destroys `PC` context that was created with `PCCreate()`.
88: Collective
90: Input Parameter:
91: . pc - the preconditioner context
93: Level: developer
95: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetUp()`
96: @*/
97: PetscErrorCode PCDestroy(PC *pc)
98: {
99: PetscFunctionBegin;
100: if (!*pc) PetscFunctionReturn(PETSC_SUCCESS);
102: if (--((PetscObject)(*pc))->refct > 0) {
103: *pc = NULL;
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
107: PetscCall(PCReset(*pc));
109: /* if memory was published with SAWs then destroy it */
110: PetscCall(PetscObjectSAWsViewOff((PetscObject)*pc));
111: PetscTryTypeMethod((*pc), destroy);
112: PetscCall(DMDestroy(&(*pc)->dm));
113: PetscCall(PetscHeaderDestroy(pc));
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }
117: /*@C
118: PCGetDiagonalScale - Indicates if the preconditioner applies an additional left and right
119: scaling as needed by certain time-stepping codes.
121: Logically Collective
123: Input Parameter:
124: . pc - the preconditioner context
126: Output Parameter:
127: . flag - `PETSC_TRUE` if it applies the scaling
129: Level: developer
131: Note:
132: If this returns `PETSC_TRUE` then the system solved via the Krylov method is, for left and right preconditioning,
134: $$
135: \begin{align*}
136: D M A D^{-1} y = D M b \\
137: D A M D^{-1} z = D b.
138: \end{align*}
139: $$
141: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCDiagonalScaleRight()`, `PCSetDiagonalScale()`
142: @*/
143: PetscErrorCode PCGetDiagonalScale(PC pc, PetscBool *flag)
144: {
145: PetscFunctionBegin;
147: PetscAssertPointer(flag, 2);
148: *flag = pc->diagonalscale;
149: PetscFunctionReturn(PETSC_SUCCESS);
150: }
152: /*@
153: PCSetDiagonalScale - Indicates the left scaling to use to apply an additional left and right
154: scaling as needed by certain time-stepping codes.
156: Logically Collective
158: Input Parameters:
159: + pc - the preconditioner context
160: - s - scaling vector
162: Level: intermediate
164: Notes:
165: The system solved via the Krylov method is, for left and right preconditioning,
166: $$
167: \begin{align*}
168: D M A D^{-1} y = D M b \\
169: D A M D^{-1} z = D b.
170: \end{align*}
171: $$
173: `PCDiagonalScaleLeft()` scales a vector by $D$. `PCDiagonalScaleRight()` scales a vector by $D^{-1}$.
175: .seealso: [](ch_ksp), `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCDiagonalScaleRight()`, `PCGetDiagonalScale()`
176: @*/
177: PetscErrorCode PCSetDiagonalScale(PC pc, Vec s)
178: {
179: PetscFunctionBegin;
182: pc->diagonalscale = PETSC_TRUE;
184: PetscCall(PetscObjectReference((PetscObject)s));
185: PetscCall(VecDestroy(&pc->diagonalscaleleft));
187: pc->diagonalscaleleft = s;
189: PetscCall(VecDuplicate(s, &pc->diagonalscaleright));
190: PetscCall(VecCopy(s, pc->diagonalscaleright));
191: PetscCall(VecReciprocal(pc->diagonalscaleright));
192: PetscFunctionReturn(PETSC_SUCCESS);
193: }
195: /*@
196: PCDiagonalScaleLeft - Scales a vector by the left scaling as needed by certain time-stepping codes.
198: Logically Collective
200: Input Parameters:
201: + pc - the preconditioner context
202: . in - input vector
203: - out - scaled vector (maybe the same as in)
205: Level: intermediate
207: Notes:
208: The system solved via the Krylov method is, for left and right preconditioning,
210: $$
211: \begin{align*}
212: D M A D^{-1} y = D M b \\
213: D A M D^{-1} z = D b.
214: \end{align*}
215: $$
217: `PCDiagonalScaleLeft()` scales a vector by $D$. `PCDiagonalScaleRight()` scales a vector by $D^{-1}$.
219: If diagonal scaling is turned off and `in` is not `out` then `in` is copied to `out`
221: .seealso: [](ch_ksp), `PCCreate()`, `PCSetUp()`, `PCSetDiagonalScale()`, `PCDiagonalScaleRight()`, `PCDiagonalScale()`
222: @*/
223: PetscErrorCode PCDiagonalScaleLeft(PC pc, Vec in, Vec out)
224: {
225: PetscFunctionBegin;
229: if (pc->diagonalscale) {
230: PetscCall(VecPointwiseMult(out, pc->diagonalscaleleft, in));
231: } else if (in != out) {
232: PetscCall(VecCopy(in, out));
233: }
234: PetscFunctionReturn(PETSC_SUCCESS);
235: }
237: /*@
238: PCDiagonalScaleRight - Scales a vector by the right scaling as needed by certain time-stepping codes.
240: Logically Collective
242: Input Parameters:
243: + pc - the preconditioner context
244: . in - input vector
245: - out - scaled vector (maybe the same as in)
247: Level: intermediate
249: Notes:
250: The system solved via the Krylov method is, for left and right preconditioning,
252: $$
253: \begin{align*}
254: D M A D^{-1} y = D M b \\
255: D A M D^{-1} z = D b.
256: \end{align*}
257: $$
259: `PCDiagonalScaleLeft()` scales a vector by $D$. `PCDiagonalScaleRight()` scales a vector by $D^{-1}$.
261: If diagonal scaling is turned off and `in` is not `out` then `in` is copied to `out`
263: .seealso: [](ch_ksp), `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCSetDiagonalScale()`, `PCDiagonalScale()`
264: @*/
265: PetscErrorCode PCDiagonalScaleRight(PC pc, Vec in, Vec out)
266: {
267: PetscFunctionBegin;
271: if (pc->diagonalscale) {
272: PetscCall(VecPointwiseMult(out, pc->diagonalscaleright, in));
273: } else if (in != out) {
274: PetscCall(VecCopy(in, out));
275: }
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: /*@
280: PCSetUseAmat - Sets a flag to indicate that when the preconditioner needs to apply (part of) the
281: operator during the preconditioning process it applies the Amat provided to `TSSetRHSJacobian()`,
282: `TSSetIJacobian()`, `SNESSetJacobian()`, `KSPSetOperators()` or `PCSetOperators()` not the Pmat.
284: Logically Collective
286: Input Parameters:
287: + pc - the preconditioner context
288: - flg - `PETSC_TRUE` to use the Amat, `PETSC_FALSE` to use the Pmat (default is false)
290: Options Database Key:
291: . -pc_use_amat <true,false> - use the amat argument to `KSPSetOperators()` or `PCSetOperators()` to apply the operator
293: Level: intermediate
295: Note:
296: For the common case in which the linear system matrix and the matrix used to construct the
297: preconditioner are identical, this routine has no affect.
299: .seealso: [](ch_ksp), `PC`, `PCGetUseAmat()`, `PCBJACOBI`, `PGMG`, `PCFIELDSPLIT`, `PCCOMPOSITE`,
300: `KSPSetOperators()`, `PCSetOperators()`
301: @*/
302: PetscErrorCode PCSetUseAmat(PC pc, PetscBool flg)
303: {
304: PetscFunctionBegin;
306: pc->useAmat = flg;
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: /*@
311: PCSetErrorIfFailure - Causes `PC` to generate an error if a floating point exception, for example a zero pivot, is detected.
313: Logically Collective
315: Input Parameters:
316: + pc - iterative context obtained from `PCCreate()`
317: - flg - `PETSC_TRUE` indicates you want the error generated
319: Level: advanced
321: Notes:
322: Normally PETSc continues if a linear solver fails due to a failed setup of a preconditioner, you can call `KSPGetConvergedReason()` after a `KSPSolve()`
323: 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
324: detected.
326: This is propagated into `KSP`s used by this `PC`, which then propagate it into `PC`s used by those `KSP`s
328: .seealso: [](ch_ksp), `PC`, `KSPSetErrorIfNotConverged()`, `PCGetInitialGuessNonzero()`, `PCSetInitialGuessKnoll()`, `PCGetInitialGuessKnoll()`
329: @*/
330: PetscErrorCode PCSetErrorIfFailure(PC pc, PetscBool flg)
331: {
332: PetscFunctionBegin;
335: pc->erroriffailure = flg;
336: PetscFunctionReturn(PETSC_SUCCESS);
337: }
339: /*@
340: PCGetUseAmat - Gets a flag to indicate that when the preconditioner needs to apply (part of) the
341: operator during the preconditioning process it applies the Amat provided to `TSSetRHSJacobian()`,
342: `TSSetIJacobian()`, `SNESSetJacobian()`, `KSPSetOperators()` or `PCSetOperators()` not the Pmat.
344: Logically Collective
346: Input Parameter:
347: . pc - the preconditioner context
349: Output Parameter:
350: . flg - `PETSC_TRUE` to use the Amat, `PETSC_FALSE` to use the Pmat (default is false)
352: Level: intermediate
354: Note:
355: For the common case in which the linear system matrix and the matrix used to construct the
356: preconditioner are identical, this routine is does nothing.
358: .seealso: [](ch_ksp), `PC`, `PCSetUseAmat()`, `PCBJACOBI`, `PGMG`, `PCFIELDSPLIT`, `PCCOMPOSITE`
359: @*/
360: PetscErrorCode PCGetUseAmat(PC pc, PetscBool *flg)
361: {
362: PetscFunctionBegin;
364: *flg = pc->useAmat;
365: PetscFunctionReturn(PETSC_SUCCESS);
366: }
368: /*@
369: PCSetKSPNestLevel - sets the amount of nesting the `KSP` that contains this `PC` has
371: Collective
373: Input Parameters:
374: + pc - the `PC`
375: - level - the nest level
377: Level: developer
379: .seealso: [](ch_ksp), `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSP`, `KSPGMRES`, `KSPType`, `KSPGetNestLevel()`, `PCGetKSPNestLevel()`, `KSPSetNestLevel()`
380: @*/
381: PetscErrorCode PCSetKSPNestLevel(PC pc, PetscInt level)
382: {
383: PetscFunctionBegin;
386: pc->kspnestlevel = level;
387: PetscFunctionReturn(PETSC_SUCCESS);
388: }
390: /*@
391: PCGetKSPNestLevel - gets the amount of nesting the `KSP` that contains this `PC` has
393: Not Collective
395: Input Parameter:
396: . pc - the `PC`
398: Output Parameter:
399: . level - the nest level
401: Level: developer
403: .seealso: [](ch_ksp), `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSP`, `KSPGMRES`, `KSPType`, `KSPSetNestLevel()`, `PCSetKSPNestLevel()`, `KSPGetNestLevel()`
404: @*/
405: PetscErrorCode PCGetKSPNestLevel(PC pc, PetscInt *level)
406: {
407: PetscFunctionBegin;
409: PetscAssertPointer(level, 2);
410: *level = pc->kspnestlevel;
411: PetscFunctionReturn(PETSC_SUCCESS);
412: }
414: /*@
415: PCCreate - Creates a preconditioner context, `PC`
417: Collective
419: Input Parameter:
420: . comm - MPI communicator
422: Output Parameter:
423: . newpc - location to put the preconditioner context
425: Level: developer
427: Note:
428: The default preconditioner for sparse matrices is `PCILU` or `PCICC` with 0 fill on one process and block Jacobi (`PCBJACOBI`) with `PCILU` or `PCICC`
429: in parallel. For dense matrices it is always `PCNONE`.
431: .seealso: [](ch_ksp), `PC`, `PCSetUp()`, `PCApply()`, `PCDestroy()`
432: @*/
433: PetscErrorCode PCCreate(MPI_Comm comm, PC *newpc)
434: {
435: PC pc;
437: PetscFunctionBegin;
438: PetscAssertPointer(newpc, 2);
439: *newpc = NULL;
440: PetscCall(PCInitializePackage());
442: PetscCall(PetscHeaderCreate(pc, PC_CLASSID, "PC", "Preconditioner", "PC", comm, PCDestroy, PCView));
444: pc->mat = NULL;
445: pc->pmat = NULL;
446: pc->setupcalled = 0;
447: pc->setfromoptionscalled = 0;
448: pc->data = NULL;
449: pc->diagonalscale = PETSC_FALSE;
450: pc->diagonalscaleleft = NULL;
451: pc->diagonalscaleright = NULL;
453: pc->modifysubmatrices = NULL;
454: pc->modifysubmatricesP = NULL;
456: *newpc = pc;
457: PetscFunctionReturn(PETSC_SUCCESS);
458: }
460: /*@
461: PCApply - Applies the preconditioner to a vector.
463: Collective
465: Input Parameters:
466: + pc - the preconditioner context
467: - x - input vector
469: Output Parameter:
470: . y - output vector
472: Level: developer
474: .seealso: [](ch_ksp), `PC`, `PCApplyTranspose()`, `PCApplyBAorAB()`
475: @*/
476: PetscErrorCode PCApply(PC pc, Vec x, Vec y)
477: {
478: PetscInt m, n, mv, nv;
480: PetscFunctionBegin;
484: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
485: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
486: /* use pmat to check vector sizes since for KSPLSQR the pmat may be of a different size than mat */
487: PetscCall(MatGetLocalSize(pc->pmat, &m, &n));
488: PetscCall(VecGetLocalSize(x, &mv));
489: PetscCall(VecGetLocalSize(y, &nv));
490: /* check pmat * y = x is feasible */
491: 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);
492: 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);
493: PetscCall(VecSetErrorIfLocked(y, 3));
495: PetscCall(PCSetUp(pc));
496: PetscCall(VecLockReadPush(x));
497: PetscCall(PetscLogEventBegin(PC_Apply, pc, x, y, 0));
498: PetscUseTypeMethod(pc, apply, x, y);
499: PetscCall(PetscLogEventEnd(PC_Apply, pc, x, y, 0));
500: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
501: PetscCall(VecLockReadPop(x));
502: PetscFunctionReturn(PETSC_SUCCESS);
503: }
505: /*@
506: PCMatApply - Applies the preconditioner to multiple vectors stored as a `MATDENSE`. Like `PCApply()`, `Y` and `X` must be different matrices.
508: Collective
510: Input Parameters:
511: + pc - the preconditioner context
512: - X - block of input vectors
514: Output Parameter:
515: . Y - block of output vectors
517: Level: developer
519: .seealso: [](ch_ksp), `PC`, `PCApply()`, `KSPMatSolve()`
520: @*/
521: PetscErrorCode PCMatApply(PC pc, Mat X, Mat Y)
522: {
523: Mat A;
524: Vec cy, cx;
525: PetscInt m1, M1, m2, M2, n1, N1, n2, N2, m3, M3, n3, N3;
526: PetscBool match;
528: PetscFunctionBegin;
532: PetscCheckSameComm(pc, 1, X, 2);
533: PetscCheckSameComm(pc, 1, Y, 3);
534: PetscCheck(Y != X, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "Y and X must be different matrices");
535: PetscCall(PCGetOperators(pc, NULL, &A));
536: PetscCall(MatGetLocalSize(A, &m3, &n3));
537: PetscCall(MatGetLocalSize(X, &m2, &n2));
538: PetscCall(MatGetLocalSize(Y, &m1, &n1));
539: PetscCall(MatGetSize(A, &M3, &N3));
540: PetscCall(MatGetSize(X, &M2, &N2));
541: PetscCall(MatGetSize(Y, &M1, &N1));
542: 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);
543: 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);
544: 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);
545: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)Y, &match, MATSEQDENSE, MATMPIDENSE, ""));
546: PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided block of output vectors not stored in a dense Mat");
547: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)X, &match, MATSEQDENSE, MATMPIDENSE, ""));
548: PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided block of input vectors not stored in a dense Mat");
549: PetscCall(PCSetUp(pc));
550: if (pc->ops->matapply) {
551: PetscCall(PetscLogEventBegin(PC_MatApply, pc, X, Y, 0));
552: PetscUseTypeMethod(pc, matapply, X, Y);
553: PetscCall(PetscLogEventEnd(PC_MatApply, pc, X, Y, 0));
554: } else {
555: PetscCall(PetscInfo(pc, "PC type %s applying column by column\n", ((PetscObject)pc)->type_name));
556: for (n1 = 0; n1 < N1; ++n1) {
557: PetscCall(MatDenseGetColumnVecRead(X, n1, &cx));
558: PetscCall(MatDenseGetColumnVecWrite(Y, n1, &cy));
559: PetscCall(PCApply(pc, cx, cy));
560: PetscCall(MatDenseRestoreColumnVecWrite(Y, n1, &cy));
561: PetscCall(MatDenseRestoreColumnVecRead(X, n1, &cx));
562: }
563: }
564: PetscFunctionReturn(PETSC_SUCCESS);
565: }
567: /*@
568: PCApplySymmetricLeft - Applies the left part of a symmetric preconditioner to a vector.
570: Collective
572: Input Parameters:
573: + pc - the preconditioner context
574: - x - input vector
576: Output Parameter:
577: . y - output vector
579: Level: developer
581: Note:
582: Currently, this routine is implemented only for `PCICC` and `PCJACOBI` preconditioners.
584: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplySymmetricRight()`
585: @*/
586: PetscErrorCode PCApplySymmetricLeft(PC pc, Vec x, Vec y)
587: {
588: PetscFunctionBegin;
592: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
593: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
594: PetscCall(PCSetUp(pc));
595: PetscCall(VecLockReadPush(x));
596: PetscCall(PetscLogEventBegin(PC_ApplySymmetricLeft, pc, x, y, 0));
597: PetscUseTypeMethod(pc, applysymmetricleft, x, y);
598: PetscCall(PetscLogEventEnd(PC_ApplySymmetricLeft, pc, x, y, 0));
599: PetscCall(VecLockReadPop(x));
600: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
601: PetscFunctionReturn(PETSC_SUCCESS);
602: }
604: /*@
605: PCApplySymmetricRight - Applies the right part of a symmetric preconditioner to a vector.
607: Collective
609: Input Parameters:
610: + pc - the preconditioner context
611: - x - input vector
613: Output Parameter:
614: . y - output vector
616: Level: developer
618: Note:
619: Currently, this routine is implemented only for `PCICC` and `PCJACOBI` preconditioners.
621: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplySymmetricLeft()`
622: @*/
623: PetscErrorCode PCApplySymmetricRight(PC pc, Vec x, Vec y)
624: {
625: PetscFunctionBegin;
629: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
630: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
631: PetscCall(PCSetUp(pc));
632: PetscCall(VecLockReadPush(x));
633: PetscCall(PetscLogEventBegin(PC_ApplySymmetricRight, pc, x, y, 0));
634: PetscUseTypeMethod(pc, applysymmetricright, x, y);
635: PetscCall(PetscLogEventEnd(PC_ApplySymmetricRight, pc, x, y, 0));
636: PetscCall(VecLockReadPop(x));
637: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
638: PetscFunctionReturn(PETSC_SUCCESS);
639: }
641: /*@
642: PCApplyTranspose - Applies the transpose of preconditioner to a vector.
644: Collective
646: Input Parameters:
647: + pc - the preconditioner context
648: - x - input vector
650: Output Parameter:
651: . y - output vector
653: Level: developer
655: Note:
656: For complex numbers this applies the non-Hermitian transpose.
658: Developer Note:
659: We need to implement a `PCApplyHermitianTranspose()`
661: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplyBAorAB()`, `PCApplyBAorABTranspose()`, `PCApplyTransposeExists()`
662: @*/
663: PetscErrorCode PCApplyTranspose(PC pc, Vec x, Vec y)
664: {
665: PetscFunctionBegin;
669: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
670: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
671: PetscCall(PCSetUp(pc));
672: PetscCall(VecLockReadPush(x));
673: PetscCall(PetscLogEventBegin(PC_Apply, pc, x, y, 0));
674: PetscUseTypeMethod(pc, applytranspose, x, y);
675: PetscCall(PetscLogEventEnd(PC_Apply, pc, x, y, 0));
676: PetscCall(VecLockReadPop(x));
677: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
678: PetscFunctionReturn(PETSC_SUCCESS);
679: }
681: /*@
682: PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation
684: Collective
686: Input Parameter:
687: . pc - the preconditioner context
689: Output Parameter:
690: . flg - `PETSC_TRUE` if a transpose operation is defined
692: Level: developer
694: .seealso: [](ch_ksp), `PC`, `PCApplyTranspose()`
695: @*/
696: PetscErrorCode PCApplyTransposeExists(PC pc, PetscBool *flg)
697: {
698: PetscFunctionBegin;
700: PetscAssertPointer(flg, 2);
701: if (pc->ops->applytranspose) *flg = PETSC_TRUE;
702: else *flg = PETSC_FALSE;
703: PetscFunctionReturn(PETSC_SUCCESS);
704: }
706: /*@
707: PCApplyBAorAB - Applies the preconditioner and operator to a vector. $y = B*A*x $ or $ y = A*B*x$.
709: Collective
711: Input Parameters:
712: + pc - the preconditioner context
713: . side - indicates the preconditioner side, one of `PC_LEFT`, `PC_RIGHT`, or `PC_SYMMETRIC`
714: . x - input vector
715: - work - work vector
717: Output Parameter:
718: . y - output vector
720: Level: developer
722: Note:
723: 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.
724: The specific `KSPSolve()` method must also be written to handle the post-solve "correction" for the diagonal scaling.
726: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplyTranspose()`, `PCApplyBAorABTranspose()`
727: @*/
728: PetscErrorCode PCApplyBAorAB(PC pc, PCSide side, Vec x, Vec y, Vec work)
729: {
730: PetscFunctionBegin;
736: PetscCheckSameComm(pc, 1, x, 3);
737: PetscCheckSameComm(pc, 1, y, 4);
738: PetscCheckSameComm(pc, 1, work, 5);
739: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
740: PetscCheck(side == PC_LEFT || side == PC_SYMMETRIC || side == PC_RIGHT, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Side must be right, left, or symmetric");
741: PetscCheck(!pc->diagonalscale || side != PC_SYMMETRIC, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot include diagonal scaling with symmetric preconditioner application");
742: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 3, PETSC_TRUE));
744: PetscCall(PCSetUp(pc));
745: if (pc->diagonalscale) {
746: if (pc->ops->applyBA) {
747: Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
748: PetscCall(VecDuplicate(x, &work2));
749: PetscCall(PCDiagonalScaleRight(pc, x, work2));
750: PetscUseTypeMethod(pc, applyBA, side, work2, y, work);
751: PetscCall(PCDiagonalScaleLeft(pc, y, y));
752: PetscCall(VecDestroy(&work2));
753: } else if (side == PC_RIGHT) {
754: PetscCall(PCDiagonalScaleRight(pc, x, y));
755: PetscCall(PCApply(pc, y, work));
756: PetscCall(MatMult(pc->mat, work, y));
757: PetscCall(PCDiagonalScaleLeft(pc, y, y));
758: } else if (side == PC_LEFT) {
759: PetscCall(PCDiagonalScaleRight(pc, x, y));
760: PetscCall(MatMult(pc->mat, y, work));
761: PetscCall(PCApply(pc, work, y));
762: PetscCall(PCDiagonalScaleLeft(pc, y, y));
763: } else PetscCheck(side != PC_SYMMETRIC, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot provide diagonal scaling with symmetric application of preconditioner");
764: } else {
765: if (pc->ops->applyBA) {
766: PetscUseTypeMethod(pc, applyBA, side, x, y, work);
767: } else if (side == PC_RIGHT) {
768: PetscCall(PCApply(pc, x, work));
769: PetscCall(MatMult(pc->mat, work, y));
770: } else if (side == PC_LEFT) {
771: PetscCall(MatMult(pc->mat, x, work));
772: PetscCall(PCApply(pc, work, y));
773: } else if (side == PC_SYMMETRIC) {
774: /* There's an extra copy here; maybe should provide 2 work vectors instead? */
775: PetscCall(PCApplySymmetricRight(pc, x, work));
776: PetscCall(MatMult(pc->mat, work, y));
777: PetscCall(VecCopy(y, work));
778: PetscCall(PCApplySymmetricLeft(pc, work, y));
779: }
780: }
781: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 4, PETSC_FALSE));
782: PetscFunctionReturn(PETSC_SUCCESS);
783: }
785: /*@
786: PCApplyBAorABTranspose - Applies the transpose of the preconditioner
787: and operator to a vector. That is, applies $B^T * A^T$ with left preconditioning,
788: NOT $(B*A)^T = A^T*B^T$.
790: Collective
792: Input Parameters:
793: + pc - the preconditioner context
794: . side - indicates the preconditioner side, one of `PC_LEFT`, `PC_RIGHT`, or `PC_SYMMETRIC`
795: . x - input vector
796: - work - work vector
798: Output Parameter:
799: . y - output vector
801: Level: developer
803: Note:
804: 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
805: defined by $B^T$. This is why this has the funny form that it computes $B^T * A^T$
807: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplyTranspose()`, `PCApplyBAorAB()`
808: @*/
809: PetscErrorCode PCApplyBAorABTranspose(PC pc, PCSide side, Vec x, Vec y, Vec work)
810: {
811: PetscFunctionBegin;
816: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
817: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 3, PETSC_TRUE));
818: if (pc->ops->applyBAtranspose) {
819: PetscUseTypeMethod(pc, applyBAtranspose, side, x, y, work);
820: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 4, PETSC_FALSE));
821: PetscFunctionReturn(PETSC_SUCCESS);
822: }
823: PetscCheck(side == PC_LEFT || side == PC_RIGHT, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Side must be right or left");
825: PetscCall(PCSetUp(pc));
826: if (side == PC_RIGHT) {
827: PetscCall(PCApplyTranspose(pc, x, work));
828: PetscCall(MatMultTranspose(pc->mat, work, y));
829: } else if (side == PC_LEFT) {
830: PetscCall(MatMultTranspose(pc->mat, x, work));
831: PetscCall(PCApplyTranspose(pc, work, y));
832: }
833: /* add support for PC_SYMMETRIC */
834: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 4, PETSC_FALSE));
835: PetscFunctionReturn(PETSC_SUCCESS);
836: }
838: /*@
839: PCApplyRichardsonExists - Determines whether a particular preconditioner has a
840: built-in fast application of Richardson's method.
842: Not Collective
844: Input Parameter:
845: . pc - the preconditioner
847: Output Parameter:
848: . exists - `PETSC_TRUE` or `PETSC_FALSE`
850: Level: developer
852: .seealso: [](ch_ksp), `PC`, `PCRICHARDSON`, `PCApplyRichardson()`
853: @*/
854: PetscErrorCode PCApplyRichardsonExists(PC pc, PetscBool *exists)
855: {
856: PetscFunctionBegin;
858: PetscAssertPointer(exists, 2);
859: if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
860: else *exists = PETSC_FALSE;
861: PetscFunctionReturn(PETSC_SUCCESS);
862: }
864: /*@
865: PCApplyRichardson - Applies several steps of Richardson iteration with
866: the particular preconditioner. This routine is usually used by the
867: Krylov solvers and not the application code directly.
869: Collective
871: Input Parameters:
872: + pc - the preconditioner context
873: . b - the right hand side
874: . w - one work vector
875: . rtol - relative decrease in residual norm convergence criteria
876: . abstol - absolute residual norm convergence criteria
877: . dtol - divergence residual norm increase criteria
878: . its - the number of iterations to apply.
879: - guesszero - if the input x contains nonzero initial guess
881: Output Parameters:
882: + outits - number of iterations actually used (for SOR this always equals its)
883: . reason - the reason the apply terminated
884: - y - the solution (also contains initial guess if guesszero is `PETSC_FALSE`
886: Level: developer
888: Notes:
889: Most preconditioners do not support this function. Use the command
890: `PCApplyRichardsonExists()` to determine if one does.
892: Except for the `PCMG` this routine ignores the convergence tolerances
893: and always runs for the number of iterations
895: .seealso: [](ch_ksp), `PC`, `PCApplyRichardsonExists()`
896: @*/
897: PetscErrorCode PCApplyRichardson(PC pc, Vec b, Vec y, Vec w, PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt its, PetscBool guesszero, PetscInt *outits, PCRichardsonConvergedReason *reason)
898: {
899: PetscFunctionBegin;
904: PetscCheck(b != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "b and y must be different vectors");
905: PetscCall(PCSetUp(pc));
906: PetscUseTypeMethod(pc, applyrichardson, b, y, w, rtol, abstol, dtol, its, guesszero, outits, reason);
907: PetscFunctionReturn(PETSC_SUCCESS);
908: }
910: /*@
911: PCSetFailedReason - Sets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail
913: Logically Collective
915: Input Parameters:
916: + pc - the preconditioner context
917: - reason - the reason it failedx
919: Level: advanced
921: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCFailedReason`
922: @*/
923: PetscErrorCode PCSetFailedReason(PC pc, PCFailedReason reason)
924: {
925: PetscFunctionBegin;
927: pc->failedreason = reason;
928: PetscFunctionReturn(PETSC_SUCCESS);
929: }
931: /*@
932: PCGetFailedReason - Gets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail
934: Logically Collective
936: Input Parameter:
937: . pc - the preconditioner context
939: Output Parameter:
940: . reason - the reason it failed
942: Level: advanced
944: Note:
945: This is the maximum over reason over all ranks in the PC communicator. It is only valid after
946: a call `KSPCheckDot()` or `KSPCheckNorm()` inside a `KSPSolve()` or `PCReduceFailedReason()`.
947: It is not valid immediately after a `PCSetUp()` or `PCApply()`, then use `PCGetFailedReasonRank()`
949: .seealso: [](ch_ksp), `PC`, ``PCCreate()`, `PCApply()`, `PCDestroy()`, `PCGetFailedReasonRank()`, `PCSetFailedReason()`
950: @*/
951: PetscErrorCode PCGetFailedReason(PC pc, PCFailedReason *reason)
952: {
953: PetscFunctionBegin;
955: if (pc->setupcalled < 0) *reason = (PCFailedReason)pc->setupcalled;
956: else *reason = pc->failedreason;
957: PetscFunctionReturn(PETSC_SUCCESS);
958: }
960: /*@
961: PCGetFailedReasonRank - Gets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail on this MPI rank
963: Not Collective
965: Input Parameter:
966: . pc - the preconditioner context
968: Output Parameter:
969: . reason - the reason it failed
971: Level: advanced
973: Note:
974: Different processes may have different reasons or no reason, see `PCGetFailedReason()`
976: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCGetFailedReason()`, `PCSetFailedReason()`, `PCReduceFailedReason()`
977: @*/
978: PetscErrorCode PCGetFailedReasonRank(PC pc, PCFailedReason *reason)
979: {
980: PetscFunctionBegin;
982: if (pc->setupcalled < 0) *reason = (PCFailedReason)pc->setupcalled;
983: else *reason = pc->failedreason;
984: PetscFunctionReturn(PETSC_SUCCESS);
985: }
987: /*@
988: PCReduceFailedReason - Reduce the failed reason among the MPI processes that share the `PC`
990: Collective
992: Input Parameter:
993: . pc - the preconditioner context
995: Level: advanced
997: Note:
998: Different MPI processes may have different reasons or no reason, see `PCGetFailedReason()`. This routine
999: makes them have a common value (failure if any MPI process had a failure).
1001: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCGetFailedReason()`, `PCSetFailedReason()`
1002: @*/
1003: PetscErrorCode PCReduceFailedReason(PC pc)
1004: {
1005: PetscInt buf;
1007: PetscFunctionBegin;
1009: buf = (PetscInt)pc->failedreason;
1010: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &buf, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)pc)));
1011: pc->failedreason = (PCFailedReason)buf;
1012: PetscFunctionReturn(PETSC_SUCCESS);
1013: }
1015: /* Next line needed to deactivate KSP_Solve logging */
1016: #include <petsc/private/kspimpl.h>
1018: /*
1019: a setupcall of 0 indicates never setup,
1020: 1 indicates has been previously setup
1021: -1 indicates a PCSetUp() was attempted and failed
1022: */
1023: /*@
1024: PCSetUp - Prepares for the use of a preconditioner.
1026: Collective
1028: Input Parameter:
1029: . pc - the preconditioner context
1031: Level: developer
1033: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`
1034: @*/
1035: PetscErrorCode PCSetUp(PC pc)
1036: {
1037: const char *def;
1038: PetscObjectState matstate, matnonzerostate;
1040: PetscFunctionBegin;
1042: PetscCheck(pc->mat, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Matrix must be set first");
1044: if (pc->setupcalled && pc->reusepreconditioner) {
1045: PetscCall(PetscInfo(pc, "Leaving PC with identical preconditioner since reuse preconditioner is set\n"));
1046: PetscFunctionReturn(PETSC_SUCCESS);
1047: }
1049: PetscCall(PetscObjectStateGet((PetscObject)pc->pmat, &matstate));
1050: PetscCall(MatGetNonzeroState(pc->pmat, &matnonzerostate));
1051: if (!pc->setupcalled) {
1052: PetscCall(PetscInfo(pc, "Setting up PC for first time\n"));
1053: pc->flag = DIFFERENT_NONZERO_PATTERN;
1054: } else if (matstate == pc->matstate) PetscFunctionReturn(PETSC_SUCCESS);
1055: else {
1056: if (matnonzerostate != pc->matnonzerostate) {
1057: PetscCall(PetscInfo(pc, "Setting up PC with different nonzero pattern\n"));
1058: pc->flag = DIFFERENT_NONZERO_PATTERN;
1059: } else {
1060: PetscCall(PetscInfo(pc, "Setting up PC with same nonzero pattern\n"));
1061: pc->flag = SAME_NONZERO_PATTERN;
1062: }
1063: }
1064: pc->matstate = matstate;
1065: pc->matnonzerostate = matnonzerostate;
1067: if (!((PetscObject)pc)->type_name) {
1068: PetscCall(PCGetDefaultType_Private(pc, &def));
1069: PetscCall(PCSetType(pc, def));
1070: }
1072: PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
1073: PetscCall(MatSetErrorIfFailure(pc->mat, pc->erroriffailure));
1074: PetscCall(PetscLogEventBegin(PC_SetUp, pc, 0, 0, 0));
1075: if (pc->ops->setup) {
1076: /* do not log solves and applications of preconditioners while constructing preconditioners; perhaps they should be logged separately from the regular solves */
1077: PetscCall(KSPInitializePackage());
1078: PetscCall(PetscLogEventDeactivatePush(KSP_Solve));
1079: PetscCall(PetscLogEventDeactivatePush(PC_Apply));
1080: PetscUseTypeMethod(pc, setup);
1081: PetscCall(PetscLogEventDeactivatePop(KSP_Solve));
1082: PetscCall(PetscLogEventDeactivatePop(PC_Apply));
1083: }
1084: PetscCall(PetscLogEventEnd(PC_SetUp, pc, 0, 0, 0));
1085: if (!pc->setupcalled) pc->setupcalled = 1;
1086: PetscFunctionReturn(PETSC_SUCCESS);
1087: }
1089: /*@
1090: PCSetUpOnBlocks - Sets up the preconditioner for each block in
1091: the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
1092: methods.
1094: Collective
1096: Input Parameter:
1097: . pc - the preconditioner context
1099: Level: developer
1101: Note:
1102: For nested preconditioners such as `PCBJACOBI` `PCSetUp()` is not called on each sub-`KSP` when `PCSetUp()` is
1103: called on the outer `PC`, this routine ensures it is called.
1105: .seealso: [](ch_ksp), `PC`, `PCSetUp()`, `PCCreate()`, `PCApply()`, `PCDestroy()`
1106: @*/
1107: PetscErrorCode PCSetUpOnBlocks(PC pc)
1108: {
1109: PetscFunctionBegin;
1111: if (!pc->ops->setuponblocks) PetscFunctionReturn(PETSC_SUCCESS);
1112: PetscCall(PetscLogEventBegin(PC_SetUpOnBlocks, pc, 0, 0, 0));
1113: PetscUseTypeMethod(pc, setuponblocks);
1114: PetscCall(PetscLogEventEnd(PC_SetUpOnBlocks, pc, 0, 0, 0));
1115: PetscFunctionReturn(PETSC_SUCCESS);
1116: }
1118: /*@C
1119: PCSetModifySubMatrices - Sets a user-defined routine for modifying the
1120: submatrices that arise within certain subdomain-based preconditioners such as `PCASM`
1122: Logically Collective
1124: Input Parameters:
1125: + pc - the preconditioner context
1126: . func - routine for modifying the submatrices
1127: - ctx - optional user-defined context (may be `NULL`)
1129: Calling sequence of `func`:
1130: + pc - the preconditioner context
1131: . nsub - number of index sets
1132: . row - an array of index sets that contain the global row numbers
1133: that comprise each local submatrix
1134: . col - an array of index sets that contain the global column numbers
1135: that comprise each local submatrix
1136: . submat - array of local submatrices
1137: - ctx - optional user-defined context for private data for the
1138: user-defined func routine (may be `NULL`)
1140: Level: advanced
1142: Notes:
1143: The basic submatrices are extracted from the preconditioner matrix as
1144: usual; the user can then alter these (for example, to set different boundary
1145: conditions for each submatrix) before they are used for the local solves.
1147: `PCSetModifySubMatrices()` MUST be called before `KSPSetUp()` and
1148: `KSPSolve()`.
1150: A routine set by `PCSetModifySubMatrices()` is currently called within
1151: the block Jacobi (`PCBJACOBI`) and additive Schwarz (`PCASM`)
1152: preconditioners. All other preconditioners ignore this routine.
1154: .seealso: [](ch_ksp), `PC`, `PCBJACOBI`, `PCASM`, `PCModifySubMatrices()`
1155: @*/
1156: PetscErrorCode PCSetModifySubMatrices(PC pc, PetscErrorCode (*func)(PC pc, PetscInt nsub, const IS row[], const IS col[], Mat submat[], void *ctx), void *ctx)
1157: {
1158: PetscFunctionBegin;
1160: pc->modifysubmatrices = func;
1161: pc->modifysubmatricesP = ctx;
1162: PetscFunctionReturn(PETSC_SUCCESS);
1163: }
1165: /*@C
1166: PCModifySubMatrices - Calls an optional user-defined routine within
1167: certain preconditioners if one has been set with `PCSetModifySubMatrices()`.
1169: Collective
1171: Input Parameters:
1172: + pc - the preconditioner context
1173: . nsub - the number of local submatrices
1174: . row - an array of index sets that contain the global row numbers
1175: that comprise each local submatrix
1176: . col - an array of index sets that contain the global column numbers
1177: that comprise each local submatrix
1178: . submat - array of local submatrices
1179: - ctx - optional user-defined context for private data for the
1180: user-defined routine (may be `NULL`)
1182: Output Parameter:
1183: . submat - array of local submatrices (the entries of which may
1184: have been modified)
1186: Level: developer
1188: Note:
1189: The user should NOT generally call this routine, as it will
1190: automatically be called within certain preconditioners.
1192: .seealso: [](ch_ksp), `PC`, `PCSetModifySubMatrices()`
1193: @*/
1194: PetscErrorCode PCModifySubMatrices(PC pc, PetscInt nsub, const IS row[], const IS col[], Mat submat[], void *ctx)
1195: {
1196: PetscFunctionBegin;
1198: if (!pc->modifysubmatrices) PetscFunctionReturn(PETSC_SUCCESS);
1199: PetscCall(PetscLogEventBegin(PC_ModifySubMatrices, pc, 0, 0, 0));
1200: PetscCall((*pc->modifysubmatrices)(pc, nsub, row, col, submat, ctx));
1201: PetscCall(PetscLogEventEnd(PC_ModifySubMatrices, pc, 0, 0, 0));
1202: PetscFunctionReturn(PETSC_SUCCESS);
1203: }
1205: /*@
1206: PCSetOperators - Sets the matrix associated with the linear system and
1207: a (possibly) different one associated with the preconditioner.
1209: Logically Collective
1211: Input Parameters:
1212: + pc - the preconditioner context
1213: . Amat - the matrix that defines the linear system
1214: - Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
1216: Level: intermediate
1218: Notes:
1219: Passing a `NULL` for `Amat` or `Pmat` removes the matrix that is currently used.
1221: If you wish to replace either `Amat` or `Pmat` but leave the other one untouched then
1222: first call `KSPGetOperators()` to get the one you wish to keep, call `PetscObjectReference()`
1223: on it and then pass it back in in your call to `KSPSetOperators()`.
1225: More Notes about Repeated Solution of Linear Systems:
1226: PETSc does NOT reset the matrix entries of either `Amat` or `Pmat`
1227: to zero after a linear solve; the user is completely responsible for
1228: matrix assembly. See the routine `MatZeroEntries()` if desiring to
1229: zero all elements of a matrix.
1231: .seealso: [](ch_ksp), `PC`, `PCGetOperators()`, `MatZeroEntries()`
1232: @*/
1233: PetscErrorCode PCSetOperators(PC pc, Mat Amat, Mat Pmat)
1234: {
1235: PetscInt m1, n1, m2, n2;
1237: PetscFunctionBegin;
1241: if (Amat) PetscCheckSameComm(pc, 1, Amat, 2);
1242: if (Pmat) PetscCheckSameComm(pc, 1, Pmat, 3);
1243: if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
1244: PetscCall(MatGetLocalSize(Amat, &m1, &n1));
1245: PetscCall(MatGetLocalSize(pc->mat, &m2, &n2));
1246: 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);
1247: PetscCall(MatGetLocalSize(Pmat, &m1, &n1));
1248: PetscCall(MatGetLocalSize(pc->pmat, &m2, &n2));
1249: 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);
1250: }
1252: if (Pmat != pc->pmat) {
1253: /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
1254: pc->matnonzerostate = -1;
1255: pc->matstate = -1;
1256: }
1258: /* reference first in case the matrices are the same */
1259: if (Amat) PetscCall(PetscObjectReference((PetscObject)Amat));
1260: PetscCall(MatDestroy(&pc->mat));
1261: if (Pmat) PetscCall(PetscObjectReference((PetscObject)Pmat));
1262: PetscCall(MatDestroy(&pc->pmat));
1263: pc->mat = Amat;
1264: pc->pmat = Pmat;
1265: PetscFunctionReturn(PETSC_SUCCESS);
1266: }
1268: /*@
1269: PCSetReusePreconditioner - reuse the current preconditioner even if the operator in the preconditioner has changed.
1271: Logically Collective
1273: Input Parameters:
1274: + pc - the preconditioner context
1275: - flag - `PETSC_TRUE` do not compute a new preconditioner, `PETSC_FALSE` do compute a new preconditioner
1277: Level: intermediate
1279: Note:
1280: Normally if a matrix inside a `PC` changes the `PC` automatically updates itself using information from the changed matrix. This option
1281: prevents this.
1283: .seealso: [](ch_ksp), `PC`, `PCGetOperators()`, `MatZeroEntries()`, `PCGetReusePreconditioner()`, `KSPSetReusePreconditioner()`
1284: @*/
1285: PetscErrorCode PCSetReusePreconditioner(PC pc, PetscBool flag)
1286: {
1287: PetscFunctionBegin;
1290: pc->reusepreconditioner = flag;
1291: PetscFunctionReturn(PETSC_SUCCESS);
1292: }
1294: /*@
1295: PCGetReusePreconditioner - Determines if the `PC` reuses the current preconditioner even if the operator in the preconditioner has changed.
1297: Not Collective
1299: Input Parameter:
1300: . pc - the preconditioner context
1302: Output Parameter:
1303: . flag - `PETSC_TRUE` do not compute a new preconditioner, `PETSC_FALSE` do compute a new preconditioner
1305: Level: intermediate
1307: .seealso: [](ch_ksp), `PC`, `PCGetOperators()`, `MatZeroEntries()`, `PCSetReusePreconditioner()`
1308: @*/
1309: PetscErrorCode PCGetReusePreconditioner(PC pc, PetscBool *flag)
1310: {
1311: PetscFunctionBegin;
1313: PetscAssertPointer(flag, 2);
1314: *flag = pc->reusepreconditioner;
1315: PetscFunctionReturn(PETSC_SUCCESS);
1316: }
1318: /*@
1319: PCGetOperators - Gets the matrix associated with the linear system and
1320: possibly a different one associated with the preconditioner.
1322: Not Collective, though parallel `Mat`s are returned if `pc` is parallel
1324: Input Parameter:
1325: . pc - the preconditioner context
1327: Output Parameters:
1328: + Amat - the matrix defining the linear system
1329: - Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.
1331: Level: intermediate
1333: Note:
1334: Does not increase the reference count of the matrices, so you should not destroy them
1336: Alternative usage: If the operators have NOT been set with `KSPSetOperators()`/`PCSetOperators()` then the operators
1337: are created in `PC` and returned to the user. In this case, if both operators
1338: mat and pmat are requested, two DIFFERENT operators will be returned. If
1339: only one is requested both operators in the PC will be the same (i.e. as
1340: if one had called `KSPSetOperators()`/`PCSetOperators()` with the same argument for both Mats).
1341: The user must set the sizes of the returned matrices and their type etc just
1342: as if the user created them with `MatCreate()`. For example,
1344: .vb
1345: KSP/PCGetOperators(ksp/pc,&Amat,NULL); is equivalent to
1346: set size, type, etc of Amat
1348: MatCreate(comm,&mat);
1349: KSP/PCSetOperators(ksp/pc,Amat,Amat);
1350: PetscObjectDereference((PetscObject)mat);
1351: set size, type, etc of Amat
1352: .ve
1354: and
1356: .vb
1357: KSP/PCGetOperators(ksp/pc,&Amat,&Pmat); is equivalent to
1358: set size, type, etc of Amat and Pmat
1360: MatCreate(comm,&Amat);
1361: MatCreate(comm,&Pmat);
1362: KSP/PCSetOperators(ksp/pc,Amat,Pmat);
1363: PetscObjectDereference((PetscObject)Amat);
1364: PetscObjectDereference((PetscObject)Pmat);
1365: set size, type, etc of Amat and Pmat
1366: .ve
1368: The rationale for this support is so that when creating a `TS`, `SNES`, or `KSP` the hierarchy
1369: of underlying objects (i.e. `SNES`, `KSP`, `PC`, `Mat`) and their lifespans can be completely
1370: managed by the top most level object (i.e. the `TS`, `SNES`, or `KSP`). Another way to look
1371: at this is when you create a `SNES` you do not NEED to create a `KSP` and attach it to
1372: the `SNES` object (the `SNES` object manages it for you). Similarly when you create a KSP
1373: you do not need to attach a `PC` to it (the `KSP` object manages the `PC` object for you).
1374: Thus, why should YOU have to create the `Mat` and attach it to the `SNES`/`KSP`/`PC`, when
1375: it can be created for you?
1377: .seealso: [](ch_ksp), `PC`, `PCSetOperators()`, `KSPGetOperators()`, `KSPSetOperators()`, `PCGetOperatorsSet()`
1378: @*/
1379: PetscErrorCode PCGetOperators(PC pc, Mat *Amat, Mat *Pmat)
1380: {
1381: PetscFunctionBegin;
1383: if (Amat) {
1384: if (!pc->mat) {
1385: if (pc->pmat && !Pmat) { /* Pmat has been set, but user did not request it, so use for Amat */
1386: pc->mat = pc->pmat;
1387: PetscCall(PetscObjectReference((PetscObject)pc->mat));
1388: } else { /* both Amat and Pmat are empty */
1389: PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &pc->mat));
1390: if (!Pmat) { /* user did NOT request Pmat, so make same as Amat */
1391: pc->pmat = pc->mat;
1392: PetscCall(PetscObjectReference((PetscObject)pc->pmat));
1393: }
1394: }
1395: }
1396: *Amat = pc->mat;
1397: }
1398: if (Pmat) {
1399: if (!pc->pmat) {
1400: if (pc->mat && !Amat) { /* Amat has been set but was not requested, so use for pmat */
1401: pc->pmat = pc->mat;
1402: PetscCall(PetscObjectReference((PetscObject)pc->pmat));
1403: } else {
1404: PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &pc->pmat));
1405: if (!Amat) { /* user did NOT request Amat, so make same as Pmat */
1406: pc->mat = pc->pmat;
1407: PetscCall(PetscObjectReference((PetscObject)pc->mat));
1408: }
1409: }
1410: }
1411: *Pmat = pc->pmat;
1412: }
1413: PetscFunctionReturn(PETSC_SUCCESS);
1414: }
1416: /*@C
1417: PCGetOperatorsSet - Determines if the matrix associated with the linear system and
1418: possibly a different one associated with the preconditioner have been set in the `PC`.
1420: Not Collective, though the results on all processes should be the same
1422: Input Parameter:
1423: . pc - the preconditioner context
1425: Output Parameters:
1426: + mat - the matrix associated with the linear system was set
1427: - pmat - matrix associated with the preconditioner was set, usually the same
1429: Level: intermediate
1431: .seealso: [](ch_ksp), `PC`, `PCSetOperators()`, `KSPGetOperators()`, `KSPSetOperators()`, `PCGetOperators()`
1432: @*/
1433: PetscErrorCode PCGetOperatorsSet(PC pc, PetscBool *mat, PetscBool *pmat)
1434: {
1435: PetscFunctionBegin;
1437: if (mat) *mat = (pc->mat) ? PETSC_TRUE : PETSC_FALSE;
1438: if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1439: PetscFunctionReturn(PETSC_SUCCESS);
1440: }
1442: /*@
1443: PCFactorGetMatrix - Gets the factored matrix from the
1444: preconditioner context. This routine is valid only for the `PCLU`,
1445: `PCILU`, `PCCHOLESKY`, and `PCICC` methods.
1447: Not Collective though `mat` is parallel if `pc` is parallel
1449: Input Parameter:
1450: . pc - the preconditioner context
1452: Output Parameters:
1453: . mat - the factored matrix
1455: Level: advanced
1457: Note:
1458: Does not increase the reference count for `mat` so DO NOT destroy it
1460: .seealso: [](ch_ksp), `PC`, `PCLU`, `PCILU`, `PCCHOLESKY`, `PCICC`
1461: @*/
1462: PetscErrorCode PCFactorGetMatrix(PC pc, Mat *mat)
1463: {
1464: PetscFunctionBegin;
1466: PetscAssertPointer(mat, 2);
1467: PetscCall(PCFactorSetUpMatSolverType(pc));
1468: PetscUseTypeMethod(pc, getfactoredmatrix, mat);
1469: PetscFunctionReturn(PETSC_SUCCESS);
1470: }
1472: /*@C
1473: PCSetOptionsPrefix - Sets the prefix used for searching for all
1474: `PC` options in the database.
1476: Logically Collective
1478: Input Parameters:
1479: + pc - the preconditioner context
1480: - prefix - the prefix string to prepend to all `PC` option requests
1482: Note:
1483: A hyphen (-) must NOT be given at the beginning of the prefix name.
1484: The first character of all runtime options is AUTOMATICALLY the
1485: hyphen.
1487: Level: advanced
1489: .seealso: [](ch_ksp), `PC`, `PCSetFromOptions`, `PCAppendOptionsPrefix()`, `PCGetOptionsPrefix()`
1490: @*/
1491: PetscErrorCode PCSetOptionsPrefix(PC pc, const char prefix[])
1492: {
1493: PetscFunctionBegin;
1495: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc, prefix));
1496: PetscFunctionReturn(PETSC_SUCCESS);
1497: }
1499: /*@C
1500: PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1501: `PC` options in the database.
1503: Logically Collective
1505: Input Parameters:
1506: + pc - the preconditioner context
1507: - prefix - the prefix string to prepend to all `PC` option requests
1509: Note:
1510: A hyphen (-) must NOT be given at the beginning of the prefix name.
1511: The first character of all runtime options is AUTOMATICALLY the
1512: hyphen.
1514: Level: advanced
1516: .seealso: [](ch_ksp), `PC`, `PCSetFromOptions`, `PCSetOptionsPrefix()`, `PCGetOptionsPrefix()`
1517: @*/
1518: PetscErrorCode PCAppendOptionsPrefix(PC pc, const char prefix[])
1519: {
1520: PetscFunctionBegin;
1522: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)pc, prefix));
1523: PetscFunctionReturn(PETSC_SUCCESS);
1524: }
1526: /*@C
1527: PCGetOptionsPrefix - Gets the prefix used for searching for all
1528: PC options in the database.
1530: Not Collective
1532: Input Parameter:
1533: . pc - the preconditioner context
1535: Output Parameter:
1536: . prefix - pointer to the prefix string used, is returned
1538: Level: advanced
1540: Fortran Note:
1541: The user should pass in a string `prefix` of
1542: sufficient length to hold the prefix.
1544: .seealso: [](ch_ksp), `PC`, `PCSetFromOptions`, `PCSetOptionsPrefix()`, `PCAppendOptionsPrefix()`
1545: @*/
1546: PetscErrorCode PCGetOptionsPrefix(PC pc, const char *prefix[])
1547: {
1548: PetscFunctionBegin;
1550: PetscAssertPointer(prefix, 2);
1551: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, prefix));
1552: PetscFunctionReturn(PETSC_SUCCESS);
1553: }
1555: /*
1556: Indicates the right hand side will be changed by KSPSolve(), this occurs for a few
1557: preconditioners including BDDC and Eisentat that transform the equations before applying
1558: the Krylov methods
1559: */
1560: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC pc, PetscBool *change)
1561: {
1562: PetscFunctionBegin;
1564: PetscAssertPointer(change, 2);
1565: *change = PETSC_FALSE;
1566: PetscTryMethod(pc, "PCPreSolveChangeRHS_C", (PC, PetscBool *), (pc, change));
1567: PetscFunctionReturn(PETSC_SUCCESS);
1568: }
1570: /*@
1571: PCPreSolve - Optional pre-solve phase, intended for any
1572: preconditioner-specific actions that must be performed before
1573: the iterative solve itself.
1575: Collective
1577: Input Parameters:
1578: + pc - the preconditioner context
1579: - ksp - the Krylov subspace context
1581: Level: developer
1583: Example Usage:
1584: .vb
1585: PCPreSolve(pc,ksp);
1586: KSPSolve(ksp,b,x);
1587: PCPostSolve(pc,ksp);
1588: .ve
1590: Notes:
1591: The pre-solve phase is distinct from the `PCSetUp()` phase.
1593: `KSPSolve()` calls this directly, so is rarely called by the user.
1595: .seealso: [](ch_ksp), `PC`, `PCPostSolve()`
1596: @*/
1597: PetscErrorCode PCPreSolve(PC pc, KSP ksp)
1598: {
1599: Vec x, rhs;
1601: PetscFunctionBegin;
1604: pc->presolvedone++;
1605: PetscCheck(pc->presolvedone <= 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot embed PCPreSolve() more than twice");
1606: PetscCall(KSPGetSolution(ksp, &x));
1607: PetscCall(KSPGetRhs(ksp, &rhs));
1609: if (pc->ops->presolve) PetscUseTypeMethod(pc, presolve, ksp, rhs, x);
1610: else if (pc->presolve) PetscCall((pc->presolve)(pc, ksp));
1611: PetscFunctionReturn(PETSC_SUCCESS);
1612: }
1614: /*@C
1615: PCSetPreSolve - Sets function used by `PCPreSolve()` which is intended for any
1616: preconditioner-specific actions that must be performed before
1617: the iterative solve itself.
1619: Logically Collective
1621: Input Parameters:
1622: + pc - the preconditioner object
1623: - presolve - the function to call before the solve
1625: Calling sequence of `presolve`:
1626: + pc - the `PC` context
1627: - ksp - the `KSP` context
1629: Level: developer
1631: .seealso: [](ch_ksp), `PC`, `PCSetUp()`, `PCPreSolve()`
1632: @*/
1633: PetscErrorCode PCSetPreSolve(PC pc, PetscErrorCode (*presolve)(PC pc, KSP ksp))
1634: {
1635: PetscFunctionBegin;
1637: pc->presolve = presolve;
1638: PetscFunctionReturn(PETSC_SUCCESS);
1639: }
1641: /*@
1642: PCPostSolve - Optional post-solve phase, intended for any
1643: preconditioner-specific actions that must be performed after
1644: the iterative solve itself.
1646: Collective
1648: Input Parameters:
1649: + pc - the preconditioner context
1650: - ksp - the Krylov subspace context
1652: Example Usage:
1653: .vb
1654: PCPreSolve(pc,ksp);
1655: KSPSolve(ksp,b,x);
1656: PCPostSolve(pc,ksp);
1657: .ve
1659: Level: developer
1661: Note:
1662: `KSPSolve()` calls this routine directly, so it is rarely called by the user.
1664: .seealso: [](ch_ksp), `PC`, `PCSetPostSolve()`, `PCSetPresolve()`, `PCPreSolve()`, `KSPSolve()`
1665: @*/
1666: PetscErrorCode PCPostSolve(PC pc, KSP ksp)
1667: {
1668: Vec x, rhs;
1670: PetscFunctionBegin;
1673: pc->presolvedone--;
1674: PetscCall(KSPGetSolution(ksp, &x));
1675: PetscCall(KSPGetRhs(ksp, &rhs));
1676: PetscTryTypeMethod(pc, postsolve, ksp, rhs, x);
1677: PetscFunctionReturn(PETSC_SUCCESS);
1678: }
1680: /*@C
1681: PCLoad - Loads a `PC` that has been stored in binary with `PCView()`.
1683: Collective
1685: Input Parameters:
1686: + newdm - the newly loaded `PC`, this needs to have been created with `PCCreate()` or
1687: some related function before a call to `PCLoad()`.
1688: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()`
1690: Level: intermediate
1692: Note:
1693: The type is determined by the data in the file, any `PCType` set into the `PC` before this call is ignored.
1695: .seealso: [](ch_ksp), `PC`, `PetscViewerBinaryOpen()`, `PCView()`, `MatLoad()`, `VecLoad()`
1696: @*/
1697: PetscErrorCode PCLoad(PC newdm, PetscViewer viewer)
1698: {
1699: PetscBool isbinary;
1700: PetscInt classid;
1701: char type[256];
1703: PetscFunctionBegin;
1706: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1707: PetscCheck(isbinary, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1709: PetscCall(PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT));
1710: PetscCheck(classid == PC_FILE_CLASSID, PetscObjectComm((PetscObject)newdm), PETSC_ERR_ARG_WRONG, "Not PC next in file");
1711: PetscCall(PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR));
1712: PetscCall(PCSetType(newdm, type));
1713: PetscTryTypeMethod(newdm, load, viewer);
1714: PetscFunctionReturn(PETSC_SUCCESS);
1715: }
1717: #include <petscdraw.h>
1718: #if defined(PETSC_HAVE_SAWS)
1719: #include <petscviewersaws.h>
1720: #endif
1722: /*@C
1723: PCViewFromOptions - View from the `PC` based on options in the options database
1725: Collective
1727: Input Parameters:
1728: + A - the `PC` context
1729: . obj - Optional object that provides the options prefix
1730: - name - command line option
1732: Level: intermediate
1734: .seealso: [](ch_ksp), `PC`, `PCView`, `PetscObjectViewFromOptions()`, `PCCreate()`
1735: @*/
1736: PetscErrorCode PCViewFromOptions(PC A, PetscObject obj, const char name[])
1737: {
1738: PetscFunctionBegin;
1740: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
1741: PetscFunctionReturn(PETSC_SUCCESS);
1742: }
1744: /*@C
1745: PCView - Prints information about the `PC`
1747: Collective
1749: Input Parameters:
1750: + pc - the `PC` context
1751: - viewer - optional visualization context
1753: Level: developer
1755: Notes:
1756: The available visualization contexts include
1757: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
1758: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
1759: output where only the first processor opens
1760: the file. All other processors send their
1761: data to the first processor to print.
1763: The user can open an alternative visualization contexts with
1764: `PetscViewerASCIIOpen()` (output to a specified file).
1766: .seealso: [](ch_ksp), `PC`, `PetscViewer`, `KSPView()`, `PetscViewerASCIIOpen()`
1767: @*/
1768: PetscErrorCode PCView(PC pc, PetscViewer viewer)
1769: {
1770: PCType cstr;
1771: PetscViewerFormat format;
1772: PetscBool iascii, isstring, isbinary, isdraw, pop = PETSC_FALSE;
1773: #if defined(PETSC_HAVE_SAWS)
1774: PetscBool issaws;
1775: #endif
1777: PetscFunctionBegin;
1779: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc), &viewer));
1781: PetscCheckSameComm(pc, 1, viewer, 2);
1783: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1784: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
1785: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1786: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1787: #if defined(PETSC_HAVE_SAWS)
1788: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws));
1789: #endif
1791: if (iascii) {
1792: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)pc, viewer));
1793: if (!pc->setupcalled) PetscCall(PetscViewerASCIIPrintf(viewer, " PC has not been set up so information may be incomplete\n"));
1794: PetscCall(PetscViewerASCIIPushTab(viewer));
1795: PetscTryTypeMethod(pc, view, viewer);
1796: PetscCall(PetscViewerASCIIPopTab(viewer));
1797: if (pc->mat) {
1798: PetscCall(PetscViewerGetFormat(viewer, &format));
1799: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
1800: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
1801: pop = PETSC_TRUE;
1802: }
1803: if (pc->pmat == pc->mat) {
1804: PetscCall(PetscViewerASCIIPrintf(viewer, " linear system matrix = precond matrix:\n"));
1805: PetscCall(PetscViewerASCIIPushTab(viewer));
1806: PetscCall(MatView(pc->mat, viewer));
1807: PetscCall(PetscViewerASCIIPopTab(viewer));
1808: } else {
1809: if (pc->pmat) {
1810: PetscCall(PetscViewerASCIIPrintf(viewer, " linear system matrix followed by preconditioner matrix:\n"));
1811: } else {
1812: PetscCall(PetscViewerASCIIPrintf(viewer, " linear system matrix:\n"));
1813: }
1814: PetscCall(PetscViewerASCIIPushTab(viewer));
1815: PetscCall(MatView(pc->mat, viewer));
1816: if (pc->pmat) PetscCall(MatView(pc->pmat, viewer));
1817: PetscCall(PetscViewerASCIIPopTab(viewer));
1818: }
1819: if (pop) PetscCall(PetscViewerPopFormat(viewer));
1820: }
1821: } else if (isstring) {
1822: PetscCall(PCGetType(pc, &cstr));
1823: PetscCall(PetscViewerStringSPrintf(viewer, " PCType: %-7.7s", cstr));
1824: PetscTryTypeMethod(pc, view, viewer);
1825: if (pc->mat) PetscCall(MatView(pc->mat, viewer));
1826: if (pc->pmat && pc->pmat != pc->mat) PetscCall(MatView(pc->pmat, viewer));
1827: } else if (isbinary) {
1828: PetscInt classid = PC_FILE_CLASSID;
1829: MPI_Comm comm;
1830: PetscMPIInt rank;
1831: char type[256];
1833: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
1834: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1835: if (rank == 0) {
1836: PetscCall(PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT));
1837: PetscCall(PetscStrncpy(type, ((PetscObject)pc)->type_name, 256));
1838: PetscCall(PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR));
1839: }
1840: PetscTryTypeMethod(pc, view, viewer);
1841: } else if (isdraw) {
1842: PetscDraw draw;
1843: char str[25];
1844: PetscReal x, y, bottom, h;
1845: PetscInt n;
1847: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1848: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
1849: if (pc->mat) {
1850: PetscCall(MatGetSize(pc->mat, &n, NULL));
1851: PetscCall(PetscSNPrintf(str, 25, "PC: %s (%" PetscInt_FMT ")", ((PetscObject)pc)->type_name, n));
1852: } else {
1853: PetscCall(PetscSNPrintf(str, 25, "PC: %s", ((PetscObject)pc)->type_name));
1854: }
1855: PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h));
1856: bottom = y - h;
1857: PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
1858: PetscTryTypeMethod(pc, view, viewer);
1859: PetscCall(PetscDrawPopCurrentPoint(draw));
1860: #if defined(PETSC_HAVE_SAWS)
1861: } else if (issaws) {
1862: PetscMPIInt rank;
1864: PetscCall(PetscObjectName((PetscObject)pc));
1865: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1866: if (!((PetscObject)pc)->amsmem && rank == 0) PetscCall(PetscObjectViewSAWs((PetscObject)pc, viewer));
1867: if (pc->mat) PetscCall(MatView(pc->mat, viewer));
1868: if (pc->pmat && pc->pmat != pc->mat) PetscCall(MatView(pc->pmat, viewer));
1869: #endif
1870: }
1871: PetscFunctionReturn(PETSC_SUCCESS);
1872: }
1874: /*@C
1875: PCRegister - Adds a method (`PCType`) to the preconditioner package.
1877: Not collective
1879: Input Parameters:
1880: + sname - name of a new user-defined solver
1881: - function - routine to create method context
1883: Example Usage:
1884: .vb
1885: PCRegister("my_solver", MySolverCreate);
1886: .ve
1888: Then, your solver can be chosen with the procedural interface via
1889: $ PCSetType(pc, "my_solver")
1890: or at runtime via the option
1891: $ -pc_type my_solver
1893: Level: advanced
1895: Note:
1896: `PCRegister()` may be called multiple times to add several user-defined preconditioners.
1898: .seealso: [](ch_ksp), `PC`, `PCType`, `PCRegisterAll()`
1899: @*/
1900: PetscErrorCode PCRegister(const char sname[], PetscErrorCode (*function)(PC))
1901: {
1902: PetscFunctionBegin;
1903: PetscCall(PCInitializePackage());
1904: PetscCall(PetscFunctionListAdd(&PCList, sname, function));
1905: PetscFunctionReturn(PETSC_SUCCESS);
1906: }
1908: static PetscErrorCode MatMult_PC(Mat A, Vec X, Vec Y)
1909: {
1910: PC pc;
1912: PetscFunctionBegin;
1913: PetscCall(MatShellGetContext(A, &pc));
1914: PetscCall(PCApply(pc, X, Y));
1915: PetscFunctionReturn(PETSC_SUCCESS);
1916: }
1918: /*@C
1919: PCComputeOperator - Computes the explicit preconditioned operator.
1921: Collective
1923: Input Parameters:
1924: + pc - the preconditioner object
1925: - mattype - the `MatType` to be used for the operator
1927: Output Parameter:
1928: . mat - the explicit preconditioned operator
1930: Level: advanced
1932: Note:
1933: This computation is done by applying the operators to columns of the identity matrix.
1934: This routine is costly in general, and is recommended for use only with relatively small systems.
1935: Currently, this routine uses a dense matrix format when `mattype` == `NULL`
1937: .seealso: [](ch_ksp), `PC`, `KSPComputeOperator()`, `MatType`
1938: @*/
1939: PetscErrorCode PCComputeOperator(PC pc, MatType mattype, Mat *mat)
1940: {
1941: PetscInt N, M, m, n;
1942: Mat A, Apc;
1944: PetscFunctionBegin;
1946: PetscAssertPointer(mat, 3);
1947: PetscCall(PCGetOperators(pc, &A, NULL));
1948: PetscCall(MatGetLocalSize(A, &m, &n));
1949: PetscCall(MatGetSize(A, &M, &N));
1950: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)pc), m, n, M, N, pc, &Apc));
1951: PetscCall(MatShellSetOperation(Apc, MATOP_MULT, (void (*)(void))MatMult_PC));
1952: PetscCall(MatComputeOperator(Apc, mattype, mat));
1953: PetscCall(MatDestroy(&Apc));
1954: PetscFunctionReturn(PETSC_SUCCESS);
1955: }
1957: /*@
1958: PCSetCoordinates - sets the coordinates of all the nodes on the local process
1960: Collective
1962: Input Parameters:
1963: + pc - the solver context
1964: . dim - the dimension of the coordinates 1, 2, or 3
1965: . nloc - the blocked size of the coordinates array
1966: - coords - the coordinates array
1968: Level: intermediate
1970: Note:
1971: `coords` is an array of the dim coordinates for the nodes on
1972: the local processor, of size `dim`*`nloc`.
1973: If there are 108 equation on a processor
1974: for a displacement finite element discretization of elasticity (so
1975: that there are nloc = 36 = 108/3 nodes) then the array must have 108
1976: double precision values (ie, 3 * 36). These x y z coordinates
1977: should be ordered for nodes 0 to N-1 like so: [ 0.x, 0.y, 0.z, 1.x,
1978: ... , N-1.z ].
1980: .seealso: [](ch_ksp), `PC`, `MatSetNearNullSpace()`
1981: @*/
1982: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[])
1983: {
1984: PetscFunctionBegin;
1987: PetscTryMethod(pc, "PCSetCoordinates_C", (PC, PetscInt, PetscInt, PetscReal[]), (pc, dim, nloc, coords));
1988: PetscFunctionReturn(PETSC_SUCCESS);
1989: }
1991: /*@
1992: PCGetInterpolations - Gets interpolation matrices for all levels (except level 0)
1994: Logically Collective
1996: Input Parameter:
1997: . pc - the precondition context
1999: Output Parameters:
2000: + num_levels - the number of levels
2001: - interpolations - the interpolation matrices (size of `num_levels`-1)
2003: Level: advanced
2005: Developer Note:
2006: Why is this here instead of in `PCMG` etc?
2008: .seealso: [](ch_ksp), `PC`, `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetInterpolation()`, `PCGetCoarseOperators()`
2009: @*/
2010: PetscErrorCode PCGetInterpolations(PC pc, PetscInt *num_levels, Mat *interpolations[])
2011: {
2012: PetscFunctionBegin;
2014: PetscAssertPointer(num_levels, 2);
2015: PetscAssertPointer(interpolations, 3);
2016: PetscUseMethod(pc, "PCGetInterpolations_C", (PC, PetscInt *, Mat *[]), (pc, num_levels, interpolations));
2017: PetscFunctionReturn(PETSC_SUCCESS);
2018: }
2020: /*@
2021: PCGetCoarseOperators - Gets coarse operator matrices for all levels (except the finest level)
2023: Logically Collective
2025: Input Parameter:
2026: . pc - the precondition context
2028: Output Parameters:
2029: + num_levels - the number of levels
2030: - coarseOperators - the coarse operator matrices (size of `num_levels`-1)
2032: Level: advanced
2034: Developer Note:
2035: Why is this here instead of in `PCMG` etc?
2037: .seealso: [](ch_ksp), `PC`, `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetRScale()`, `PCMGGetInterpolation()`, `PCGetInterpolations()`
2038: @*/
2039: PetscErrorCode PCGetCoarseOperators(PC pc, PetscInt *num_levels, Mat *coarseOperators[])
2040: {
2041: PetscFunctionBegin;
2043: PetscAssertPointer(num_levels, 2);
2044: PetscAssertPointer(coarseOperators, 3);
2045: PetscUseMethod(pc, "PCGetCoarseOperators_C", (PC, PetscInt *, Mat *[]), (pc, num_levels, coarseOperators));
2046: PetscFunctionReturn(PETSC_SUCCESS);
2047: }