Actual source code: precon.c
2: /*
3: The PC (preconditioner) interface routines, callable by users.
4: */
5: #include <petsc/private/pcimpl.h>
6: #include <petscdm.h>
8: /* Logging support */
9: PetscClassId PC_CLASSID;
10: PetscLogEvent PC_SetUp, PC_SetUpOnBlocks, PC_Apply, PC_MatApply, PC_ApplyCoarse, PC_ApplyMultiple, PC_ApplySymmetricLeft;
11: PetscLogEvent PC_ApplySymmetricRight, PC_ModifySubMatrices, PC_ApplyOnBlocks, PC_ApplyTransposeOnBlocks;
12: PetscInt PetscMGLevelId;
13: PetscLogStage PCMPIStage;
15: PetscErrorCode PCGetDefaultType_Private(PC pc, const char *type[])
16: {
17: PetscMPIInt size;
18: PetscBool hasop, flg1, flg2, set, flg3, isnormal;
20: PetscFunctionBegin;
21: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
22: if (pc->pmat) {
23: PetscCall(MatHasOperation(pc->pmat, MATOP_GET_DIAGONAL_BLOCK, &hasop));
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 (hasop) { /* likely is a parallel matrix run on one processor */
36: *type = PCBJACOBI;
37: } else {
38: *type = PCNONE;
39: }
40: } else {
41: if (hasop) {
42: *type = PCBJACOBI;
43: } else {
44: *type = PCNONE;
45: }
46: }
47: } else {
48: if (size == 1) {
49: *type = PCILU;
50: } else {
51: *type = PCBJACOBI;
52: }
53: }
54: PetscFunctionReturn(PETSC_SUCCESS);
55: }
57: /*@
58: PCReset - Resets a PC context to the pcsetupcalled = 0 state and removes any allocated `Vec`s and `Mat`s
60: Collective
62: Input Parameter:
63: . pc - the preconditioner context
65: Level: developer
67: Note:
68: This allows a `PC` to be reused for a different sized linear system but using the same options that have been previously set in the PC
70: .seealso: `PC`, `PCCreate()`, `PCSetUp()`
71: @*/
72: PetscErrorCode PCReset(PC pc)
73: {
74: PetscFunctionBegin;
76: PetscTryTypeMethod(pc, reset);
77: PetscCall(VecDestroy(&pc->diagonalscaleright));
78: PetscCall(VecDestroy(&pc->diagonalscaleleft));
79: PetscCall(MatDestroy(&pc->pmat));
80: PetscCall(MatDestroy(&pc->mat));
82: pc->setupcalled = 0;
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: /*@C
87: PCDestroy - Destroys `PC` context that was created with `PCCreate()`.
89: Collective
91: Input Parameter:
92: . pc - the preconditioner context
94: Level: developer
96: .seealso: `PC`, `PCCreate()`, `PCSetUp()`
97: @*/
98: PetscErrorCode PCDestroy(PC *pc)
99: {
100: PetscFunctionBegin;
101: if (!*pc) PetscFunctionReturn(PETSC_SUCCESS);
103: if (--((PetscObject)(*pc))->refct > 0) {
104: *pc = NULL;
105: PetscFunctionReturn(PETSC_SUCCESS);
106: }
108: PetscCall(PCReset(*pc));
110: /* if memory was published with SAWs then destroy it */
111: PetscCall(PetscObjectSAWsViewOff((PetscObject)*pc));
112: PetscTryTypeMethod((*pc), destroy);
113: PetscCall(DMDestroy(&(*pc)->dm));
114: PetscCall(PetscHeaderDestroy(pc));
115: PetscFunctionReturn(PETSC_SUCCESS);
116: }
118: /*@C
119: PCGetDiagonalScale - Indicates if the preconditioner applies an additional left and right
120: scaling as needed by certain time-stepping codes.
122: Logically Collective
124: Input Parameter:
125: . pc - the preconditioner context
127: Output Parameter:
128: . flag - `PETSC_TRUE` if it applies the scaling
130: Level: developer
132: Note:
133: If this returns `PETSC_TRUE` then the system solved via the Krylov method is
134: .vb
135: D M A D^{-1} y = D M b for left preconditioning or
136: D A M D^{-1} z = D b for right preconditioning
137: .ve
139: .seealso: `PC`, `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCDiagonalScaleRight()`, `PCSetDiagonalScale()`
140: @*/
141: PetscErrorCode PCGetDiagonalScale(PC pc, PetscBool *flag)
142: {
143: PetscFunctionBegin;
146: *flag = pc->diagonalscale;
147: PetscFunctionReturn(PETSC_SUCCESS);
148: }
150: /*@
151: PCSetDiagonalScale - Indicates the left scaling to use to apply an additional left and right
152: scaling as needed by certain time-stepping codes.
154: Logically Collective
156: Input Parameters:
157: + pc - the preconditioner context
158: - s - scaling vector
160: Level: intermediate
162: Notes:
163: The system solved via the Krylov method is
164: .vb
165: D M A D^{-1} y = D M b for left preconditioning or
166: D A M D^{-1} z = D b for right preconditioning
167: .ve
169: `PCDiagonalScaleLeft()` scales a vector by D. `PCDiagonalScaleRight()` scales a vector by D^{-1}.
171: .seealso: `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCDiagonalScaleRight()`, `PCGetDiagonalScale()`
172: @*/
173: PetscErrorCode PCSetDiagonalScale(PC pc, Vec s)
174: {
175: PetscFunctionBegin;
178: pc->diagonalscale = PETSC_TRUE;
180: PetscCall(PetscObjectReference((PetscObject)s));
181: PetscCall(VecDestroy(&pc->diagonalscaleleft));
183: pc->diagonalscaleleft = s;
185: PetscCall(VecDuplicate(s, &pc->diagonalscaleright));
186: PetscCall(VecCopy(s, pc->diagonalscaleright));
187: PetscCall(VecReciprocal(pc->diagonalscaleright));
188: PetscFunctionReturn(PETSC_SUCCESS);
189: }
191: /*@
192: PCDiagonalScaleLeft - Scales a vector by the left scaling as needed by certain time-stepping codes.
194: Logically Collective
196: Input Parameters:
197: + pc - the preconditioner context
198: . in - input vector
199: - out - scaled vector (maybe the same as in)
201: Level: intermediate
203: Notes:
204: The system solved via the Krylov method is
205: .vb
206: D M A D^{-1} y = D M b for left preconditioning or
207: D A M D^{-1} z = D b for right preconditioning
208: .ve
210: `PCDiagonalScaleLeft()` scales a vector by D. `PCDiagonalScaleRight()` scales a vector by D^{-1}.
212: If diagonal scaling is turned off and in is not out then in is copied to out
214: .seealso: `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleSet()`, `PCDiagonalScaleRight()`, `PCDiagonalScale()`
215: @*/
216: PetscErrorCode PCDiagonalScaleLeft(PC pc, Vec in, Vec out)
217: {
218: PetscFunctionBegin;
222: if (pc->diagonalscale) {
223: PetscCall(VecPointwiseMult(out, pc->diagonalscaleleft, in));
224: } else if (in != out) {
225: PetscCall(VecCopy(in, out));
226: }
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
230: /*@
231: PCDiagonalScaleRight - Scales a vector by the right scaling as needed by certain time-stepping codes.
233: Logically Collective
235: Input Parameters:
236: + pc - the preconditioner context
237: . in - input vector
238: - out - scaled vector (maybe the same as in)
240: Level: intermediate
242: Notes:
243: The system solved via the Krylov method is
244: .vb
245: D M A D^{-1} y = D M b for left preconditioning or
246: D A M D^{-1} z = D b for right preconditioning
247: .ve
249: `PCDiagonalScaleLeft()` scales a vector by D. `PCDiagonalScaleRight()` scales a vector by D^{-1}.
251: If diagonal scaling is turned off and in is not out then in is copied to out
253: .seealso: `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCDiagonalScaleSet()`, `PCDiagonalScale()`
254: @*/
255: PetscErrorCode PCDiagonalScaleRight(PC pc, Vec in, Vec out)
256: {
257: PetscFunctionBegin;
261: if (pc->diagonalscale) {
262: PetscCall(VecPointwiseMult(out, pc->diagonalscaleright, in));
263: } else if (in != out) {
264: PetscCall(VecCopy(in, out));
265: }
266: PetscFunctionReturn(PETSC_SUCCESS);
267: }
269: /*@
270: PCSetUseAmat - Sets a flag to indicate that when the preconditioner needs to apply (part of) the
271: operator during the preconditioning process it applies the Amat provided to `TSSetRHSJacobian()`,
272: `TSSetIJacobian()`, `SNESSetJacobian()`, `KSPSetOperators()` or `PCSetOperators()` not the Pmat.
274: Logically Collective
276: Input Parameters:
277: + pc - the preconditioner context
278: - flg - `PETSC_TRUE` to use the Amat, `PETSC_FALSE` to use the Pmat (default is false)
280: Options Database Key:
281: . -pc_use_amat <true,false> - use the amat to apply the operator
283: Level: intermediate
285: Note:
286: For the common case in which the linear system matrix and the matrix used to construct the
287: preconditioner are identical, this routine is does nothing.
289: .seealso: `PC`, `PCGetUseAmat()`, `PCBJACOBI`, `PGMG`, `PCFIELDSPLIT`, `PCCOMPOSITE`
290: @*/
291: PetscErrorCode PCSetUseAmat(PC pc, PetscBool flg)
292: {
293: PetscFunctionBegin;
295: pc->useAmat = flg;
296: PetscFunctionReturn(PETSC_SUCCESS);
297: }
299: /*@
300: PCSetErrorIfFailure - Causes `PC` to generate an error if a FPE, for example a zero pivot, is detected.
302: Logically Collective
304: Input Parameters:
305: + pc - iterative context obtained from PCCreate()
306: - flg - `PETSC_TRUE` indicates you want the error generated
308: Level: advanced
310: Notes:
311: Normally PETSc continues if a linear solver fails due to a failed setup of a preconditioner, you can call `KSPGetConvergedReason()` after a `KSPSolve()`
312: 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
313: detected.
315: This is propagated into KSPs used by this PC, which then propagate it into PCs used by those KSPs
317: .seealso: `PC`, `KSPSetErrorIfNotConverged()`, `PCGetInitialGuessNonzero()`, `PCSetInitialGuessKnoll()`, `PCGetInitialGuessKnoll()`
318: @*/
319: PetscErrorCode PCSetErrorIfFailure(PC pc, PetscBool flg)
320: {
321: PetscFunctionBegin;
324: pc->erroriffailure = flg;
325: PetscFunctionReturn(PETSC_SUCCESS);
326: }
328: /*@
329: PCGetUseAmat - Gets a flag to indicate that when the preconditioner needs to apply (part of) the
330: operator during the preconditioning process it applies the Amat provided to `TSSetRHSJacobian()`,
331: `TSSetIJacobian()`, `SNESSetJacobian()`, `KSPSetOperators()` or `PCSetOperators()` not the Pmat.
333: Logically Collective
335: Input Parameter:
336: . pc - the preconditioner context
338: Output Parameter:
339: . flg - `PETSC_TRUE` to use the Amat, `PETSC_FALSE` to use the Pmat (default is false)
341: Level: intermediate
343: Note:
344: For the common case in which the linear system matrix and the matrix used to construct the
345: preconditioner are identical, this routine is does nothing.
347: .seealso: `PC`, `PCSetUseAmat()`, `PCBJACOBI`, `PGMG`, `PCFIELDSPLIT`, `PCCOMPOSITE`
348: @*/
349: PetscErrorCode PCGetUseAmat(PC pc, PetscBool *flg)
350: {
351: PetscFunctionBegin;
353: *flg = pc->useAmat;
354: PetscFunctionReturn(PETSC_SUCCESS);
355: }
357: /*@
358: PCCreate - Creates a preconditioner context, `PC`
360: Collective
362: Input Parameter:
363: . comm - MPI communicator
365: Output Parameter:
366: . pc - location to put the preconditioner context
368: Level: developer
370: Note:
371: The default preconditioner for sparse matrices is `PCILU` or `PCICC` with 0 fill on one process and block Jacobi (`PCBJACOBI`) with `PCILU` or `PCICC`
372: in parallel. For dense matrices it is always `PCNONE`.
374: .seealso: `PC`, `PCSetUp()`, `PCApply()`, `PCDestroy()`
375: @*/
376: PetscErrorCode PCCreate(MPI_Comm comm, PC *newpc)
377: {
378: PC pc;
380: PetscFunctionBegin;
382: *newpc = NULL;
383: PetscCall(PCInitializePackage());
385: PetscCall(PetscHeaderCreate(pc, PC_CLASSID, "PC", "Preconditioner", "PC", comm, PCDestroy, PCView));
387: pc->mat = NULL;
388: pc->pmat = NULL;
389: pc->setupcalled = 0;
390: pc->setfromoptionscalled = 0;
391: pc->data = NULL;
392: pc->diagonalscale = PETSC_FALSE;
393: pc->diagonalscaleleft = NULL;
394: pc->diagonalscaleright = NULL;
396: pc->modifysubmatrices = NULL;
397: pc->modifysubmatricesP = NULL;
399: *newpc = pc;
400: PetscFunctionReturn(PETSC_SUCCESS);
401: }
403: /*@
404: PCApply - Applies the preconditioner to a vector.
406: Collective
408: Input Parameters:
409: + pc - the preconditioner context
410: - x - input vector
412: Output Parameter:
413: . y - output vector
415: Level: developer
417: .seealso: `PC`, `PCApplyTranspose()`, `PCApplyBAorAB()`
418: @*/
419: PetscErrorCode PCApply(PC pc, Vec x, Vec y)
420: {
421: PetscInt m, n, mv, nv;
423: PetscFunctionBegin;
427: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
428: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
429: /* use pmat to check vector sizes since for KSPLSQR the pmat may be of a different size than mat */
430: PetscCall(MatGetLocalSize(pc->pmat, &m, &n));
431: PetscCall(VecGetLocalSize(x, &mv));
432: PetscCall(VecGetLocalSize(y, &nv));
433: /* check pmat * y = x is feasible */
434: 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);
435: 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);
436: PetscCall(VecSetErrorIfLocked(y, 3));
438: PetscCall(PCSetUp(pc));
439: PetscCall(VecLockReadPush(x));
440: PetscCall(PetscLogEventBegin(PC_Apply, pc, x, y, 0));
441: PetscUseTypeMethod(pc, apply, x, y);
442: PetscCall(PetscLogEventEnd(PC_Apply, pc, x, y, 0));
443: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
444: PetscCall(VecLockReadPop(x));
445: PetscFunctionReturn(PETSC_SUCCESS);
446: }
448: /*@
449: PCMatApply - Applies the preconditioner to multiple vectors stored as a `MATDENSE`. Like `PCApply()`, Y and X must be different matrices.
451: Collective
453: Input Parameters:
454: + pc - the preconditioner context
455: - X - block of input vectors
457: Output Parameter:
458: . Y - block of output vectors
460: Level: developer
462: .seealso: `PC`, `PCApply()`, `KSPMatSolve()`
463: @*/
464: PetscErrorCode PCMatApply(PC pc, Mat X, Mat Y)
465: {
466: Mat A;
467: Vec cy, cx;
468: PetscInt m1, M1, m2, M2, n1, N1, n2, N2, m3, M3, n3, N3;
469: PetscBool match;
471: PetscFunctionBegin;
475: PetscCheckSameComm(pc, 1, X, 2);
476: PetscCheckSameComm(pc, 1, Y, 3);
477: PetscCheck(Y != X, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "Y and X must be different matrices");
478: PetscCall(PCGetOperators(pc, NULL, &A));
479: PetscCall(MatGetLocalSize(A, &m3, &n3));
480: PetscCall(MatGetLocalSize(X, &m2, &n2));
481: PetscCall(MatGetLocalSize(Y, &m1, &n1));
482: PetscCall(MatGetSize(A, &M3, &N3));
483: PetscCall(MatGetSize(X, &M2, &N2));
484: PetscCall(MatGetSize(Y, &M1, &N1));
485: 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);
486: 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);
487: 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);
488: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)Y, &match, MATSEQDENSE, MATMPIDENSE, ""));
489: PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided block of output vectors not stored in a dense Mat");
490: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)X, &match, MATSEQDENSE, MATMPIDENSE, ""));
491: PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided block of input vectors not stored in a dense Mat");
492: PetscCall(PCSetUp(pc));
493: if (pc->ops->matapply) {
494: PetscCall(PetscLogEventBegin(PC_MatApply, pc, X, Y, 0));
495: PetscUseTypeMethod(pc, matapply, X, Y);
496: PetscCall(PetscLogEventEnd(PC_MatApply, pc, X, Y, 0));
497: } else {
498: PetscCall(PetscInfo(pc, "PC type %s applying column by column\n", ((PetscObject)pc)->type_name));
499: for (n1 = 0; n1 < N1; ++n1) {
500: PetscCall(MatDenseGetColumnVecRead(X, n1, &cx));
501: PetscCall(MatDenseGetColumnVecWrite(Y, n1, &cy));
502: PetscCall(PCApply(pc, cx, cy));
503: PetscCall(MatDenseRestoreColumnVecWrite(Y, n1, &cy));
504: PetscCall(MatDenseRestoreColumnVecRead(X, n1, &cx));
505: }
506: }
507: PetscFunctionReturn(PETSC_SUCCESS);
508: }
510: /*@
511: PCApplySymmetricLeft - Applies the left part of a symmetric preconditioner to a vector.
513: Collective
515: Input Parameters:
516: + pc - the preconditioner context
517: - x - input vector
519: Output Parameter:
520: . y - output vector
522: Level: developer
524: Note:
525: Currently, this routine is implemented only for `PCICC` and `PCJACOBI` preconditioners.
527: .seealso: `PC`, `PCApply()`, `PCApplySymmetricRight()`
528: @*/
529: PetscErrorCode PCApplySymmetricLeft(PC pc, Vec x, Vec y)
530: {
531: PetscFunctionBegin;
535: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
536: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
537: PetscCall(PCSetUp(pc));
538: PetscCall(VecLockReadPush(x));
539: PetscCall(PetscLogEventBegin(PC_ApplySymmetricLeft, pc, x, y, 0));
540: PetscUseTypeMethod(pc, applysymmetricleft, x, y);
541: PetscCall(PetscLogEventEnd(PC_ApplySymmetricLeft, pc, x, y, 0));
542: PetscCall(VecLockReadPop(x));
543: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
544: PetscFunctionReturn(PETSC_SUCCESS);
545: }
547: /*@
548: PCApplySymmetricRight - Applies the right part of a symmetric preconditioner to a vector.
550: Collective
552: Input Parameters:
553: + pc - the preconditioner context
554: - x - input vector
556: Output Parameter:
557: . y - output vector
559: Level: developer
561: Note:
562: Currently, this routine is implemented only for `PCICC` and `PCJACOBI` preconditioners.
564: .seealso: `PC`, `PCApply()`, `PCApplySymmetricLeft()`
565: @*/
566: PetscErrorCode PCApplySymmetricRight(PC pc, Vec x, Vec y)
567: {
568: PetscFunctionBegin;
572: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
573: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
574: PetscCall(PCSetUp(pc));
575: PetscCall(VecLockReadPush(x));
576: PetscCall(PetscLogEventBegin(PC_ApplySymmetricRight, pc, x, y, 0));
577: PetscUseTypeMethod(pc, applysymmetricright, x, y);
578: PetscCall(PetscLogEventEnd(PC_ApplySymmetricRight, pc, x, y, 0));
579: PetscCall(VecLockReadPop(x));
580: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
581: PetscFunctionReturn(PETSC_SUCCESS);
582: }
584: /*@
585: PCApplyTranspose - Applies the transpose of preconditioner to a vector.
587: Collective
589: Input Parameters:
590: + pc - the preconditioner context
591: - x - input vector
593: Output Parameter:
594: . y - output vector
596: Level: developer
598: Note:
599: For complex numbers this applies the non-Hermitian transpose.
601: Developer Note:
602: We need to implement a `PCApplyHermitianTranspose()`
604: .seealso: `PC`, `PCApply()`, `PCApplyBAorAB()`, `PCApplyBAorABTranspose()`, `PCApplyTransposeExists()`
605: @*/
606: PetscErrorCode PCApplyTranspose(PC pc, Vec x, Vec y)
607: {
608: PetscFunctionBegin;
612: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
613: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
614: PetscCall(PCSetUp(pc));
615: PetscCall(VecLockReadPush(x));
616: PetscCall(PetscLogEventBegin(PC_Apply, pc, x, y, 0));
617: PetscUseTypeMethod(pc, applytranspose, x, y);
618: PetscCall(PetscLogEventEnd(PC_Apply, pc, x, y, 0));
619: PetscCall(VecLockReadPop(x));
620: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
621: PetscFunctionReturn(PETSC_SUCCESS);
622: }
624: /*@
625: PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation
627: Collective
629: Input Parameter:
630: . pc - the preconditioner context
632: Output Parameter:
633: . flg - `PETSC_TRUE` if a transpose operation is defined
635: Level: developer
637: .seealso: `PC`, `PCApplyTranspose()`
638: @*/
639: PetscErrorCode PCApplyTransposeExists(PC pc, PetscBool *flg)
640: {
641: PetscFunctionBegin;
644: if (pc->ops->applytranspose) *flg = PETSC_TRUE;
645: else *flg = PETSC_FALSE;
646: PetscFunctionReturn(PETSC_SUCCESS);
647: }
649: /*@
650: PCApplyBAorAB - Applies the preconditioner and operator to a vector. y = B*A*x or y = A*B*x.
652: Collective
654: Input Parameters:
655: + pc - the preconditioner context
656: . side - indicates the preconditioner side, one of `PC_LEFT`, `PC_RIGHT`, or `PC_SYMMETRIC`
657: . x - input vector
658: - work - work vector
660: Output Parameter:
661: . y - output vector
663: Level: developer
665: Note:
666: 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. Note that the
667: specific `KSPSolve()` method must also be written to handle the post-solve "correction" for the diagonal scaling.
669: .seealso: `PC`, `PCApply()`, `PCApplyTranspose()`, `PCApplyBAorABTranspose()`
670: @*/
671: PetscErrorCode PCApplyBAorAB(PC pc, PCSide side, Vec x, Vec y, Vec work)
672: {
673: PetscFunctionBegin;
679: PetscCheckSameComm(pc, 1, x, 3);
680: PetscCheckSameComm(pc, 1, y, 4);
681: PetscCheckSameComm(pc, 1, work, 5);
682: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
683: PetscCheck(side == PC_LEFT || side == PC_SYMMETRIC || side == PC_RIGHT, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Side must be right, left, or symmetric");
684: PetscCheck(!pc->diagonalscale || side != PC_SYMMETRIC, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot include diagonal scaling with symmetric preconditioner application");
685: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 3, PETSC_TRUE));
687: PetscCall(PCSetUp(pc));
688: if (pc->diagonalscale) {
689: if (pc->ops->applyBA) {
690: Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
691: PetscCall(VecDuplicate(x, &work2));
692: PetscCall(PCDiagonalScaleRight(pc, x, work2));
693: PetscUseTypeMethod(pc, applyBA, side, work2, y, work);
694: PetscCall(PCDiagonalScaleLeft(pc, y, y));
695: PetscCall(VecDestroy(&work2));
696: } else if (side == PC_RIGHT) {
697: PetscCall(PCDiagonalScaleRight(pc, x, y));
698: PetscCall(PCApply(pc, y, work));
699: PetscCall(MatMult(pc->mat, work, y));
700: PetscCall(PCDiagonalScaleLeft(pc, y, y));
701: } else if (side == PC_LEFT) {
702: PetscCall(PCDiagonalScaleRight(pc, x, y));
703: PetscCall(MatMult(pc->mat, y, work));
704: PetscCall(PCApply(pc, work, y));
705: PetscCall(PCDiagonalScaleLeft(pc, y, y));
706: } else PetscCheck(side != PC_SYMMETRIC, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot provide diagonal scaling with symmetric application of preconditioner");
707: } else {
708: if (pc->ops->applyBA) {
709: PetscUseTypeMethod(pc, applyBA, side, x, y, work);
710: } else if (side == PC_RIGHT) {
711: PetscCall(PCApply(pc, x, work));
712: PetscCall(MatMult(pc->mat, work, y));
713: } else if (side == PC_LEFT) {
714: PetscCall(MatMult(pc->mat, x, work));
715: PetscCall(PCApply(pc, work, y));
716: } else if (side == PC_SYMMETRIC) {
717: /* There's an extra copy here; maybe should provide 2 work vectors instead? */
718: PetscCall(PCApplySymmetricRight(pc, x, work));
719: PetscCall(MatMult(pc->mat, work, y));
720: PetscCall(VecCopy(y, work));
721: PetscCall(PCApplySymmetricLeft(pc, work, y));
722: }
723: }
724: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 4, PETSC_FALSE));
725: PetscFunctionReturn(PETSC_SUCCESS);
726: }
728: /*@
729: PCApplyBAorABTranspose - Applies the transpose of the preconditioner
730: and operator to a vector. That is, applies tr(B) * tr(A) with left preconditioning,
731: NOT tr(B*A) = tr(A)*tr(B).
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: This routine is used internally so that the same Krylov code can be used to solve A x = b and A' x = b, with a preconditioner
748: defined by B'. This is why this has the funny form that it computes tr(B) * tr(A)
750: .seealso: `PC`, `PCApply()`, `PCApplyTranspose()`, `PCApplyBAorAB()`
751: @*/
752: PetscErrorCode PCApplyBAorABTranspose(PC pc, PCSide side, Vec x, Vec y, Vec work)
753: {
754: PetscFunctionBegin;
759: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
760: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 3, PETSC_TRUE));
761: if (pc->ops->applyBAtranspose) {
762: PetscUseTypeMethod(pc, applyBAtranspose, side, x, y, work);
763: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 4, PETSC_FALSE));
764: PetscFunctionReturn(PETSC_SUCCESS);
765: }
766: PetscCheck(side == PC_LEFT || side == PC_RIGHT, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Side must be right or left");
768: PetscCall(PCSetUp(pc));
769: if (side == PC_RIGHT) {
770: PetscCall(PCApplyTranspose(pc, x, work));
771: PetscCall(MatMultTranspose(pc->mat, work, y));
772: } else if (side == PC_LEFT) {
773: PetscCall(MatMultTranspose(pc->mat, x, work));
774: PetscCall(PCApplyTranspose(pc, work, y));
775: }
776: /* add support for PC_SYMMETRIC */
777: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 4, PETSC_FALSE));
778: PetscFunctionReturn(PETSC_SUCCESS);
779: }
781: /*@
782: PCApplyRichardsonExists - Determines whether a particular preconditioner has a
783: built-in fast application of Richardson's method.
785: Not Collective
787: Input Parameter:
788: . pc - the preconditioner
790: Output Parameter:
791: . exists - `PETSC_TRUE` or `PETSC_FALSE`
793: Level: developer
795: .seealso: `PC`, `PCRICHARDSON`, `PCApplyRichardson()`
796: @*/
797: PetscErrorCode PCApplyRichardsonExists(PC pc, PetscBool *exists)
798: {
799: PetscFunctionBegin;
802: if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
803: else *exists = PETSC_FALSE;
804: PetscFunctionReturn(PETSC_SUCCESS);
805: }
807: /*@
808: PCApplyRichardson - Applies several steps of Richardson iteration with
809: the particular preconditioner. This routine is usually used by the
810: Krylov solvers and not the application code directly.
812: Collective
814: Input Parameters:
815: + pc - the preconditioner context
816: . b - the right hand side
817: . w - one work vector
818: . rtol - relative decrease in residual norm convergence criteria
819: . abstol - absolute residual norm convergence criteria
820: . dtol - divergence residual norm increase criteria
821: . its - the number of iterations to apply.
822: - guesszero - if the input x contains nonzero initial guess
824: Output Parameters:
825: + outits - number of iterations actually used (for SOR this always equals its)
826: . reason - the reason the apply terminated
827: - y - the solution (also contains initial guess if guesszero is `PETSC_FALSE`
829: Level: developer
831: Notes:
832: Most preconditioners do not support this function. Use the command
833: `PCApplyRichardsonExists()` to determine if one does.
835: Except for the `PCMG` this routine ignores the convergence tolerances
836: and always runs for the number of iterations
838: .seealso: `PC`, `PCApplyRichardsonExists()`
839: @*/
840: PetscErrorCode PCApplyRichardson(PC pc, Vec b, Vec y, Vec w, PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt its, PetscBool guesszero, PetscInt *outits, PCRichardsonConvergedReason *reason)
841: {
842: PetscFunctionBegin;
847: PetscCheck(b != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "b and y must be different vectors");
848: PetscCall(PCSetUp(pc));
849: PetscUseTypeMethod(pc, applyrichardson, b, y, w, rtol, abstol, dtol, its, guesszero, outits, reason);
850: PetscFunctionReturn(PETSC_SUCCESS);
851: }
853: /*@
854: PCSetFailedReason - Sets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail
856: Logically Collective
858: Input Parameters:
859: + pc - the preconditioner context
860: - reason - the reason it failedx
862: Level: advanced
864: .seealso: `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCFailedReason`
865: @*/
866: PetscErrorCode PCSetFailedReason(PC pc, PCFailedReason reason)
867: {
868: PetscFunctionBegin;
869: pc->failedreason = reason;
870: PetscFunctionReturn(PETSC_SUCCESS);
871: }
873: /*@
874: PCGetFailedReason - Gets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail
876: Logically Collective
878: Input Parameter:
879: . pc - the preconditioner context
881: Output Parameter:
882: . reason - the reason it failed
884: Level: advanced
886: Note:
887: This is the maximum over reason over all ranks in the PC communicator. It is only valid after
888: a call `KSPCheckDot()` or `KSPCheckNorm()` inside a `KSPSolve()`. It is not valid immediately after a `PCSetUp()`
889: or `PCApply()`, then use `PCGetFailedReasonRank()`
891: .seealso: PC`, ``PCCreate()`, `PCApply()`, `PCDestroy()`, `PCGetFailedReasonRank()`, `PCSetFailedReason()`
892: @*/
893: PetscErrorCode PCGetFailedReason(PC pc, PCFailedReason *reason)
894: {
895: PetscFunctionBegin;
896: if (pc->setupcalled < 0) *reason = (PCFailedReason)pc->setupcalled;
897: else *reason = pc->failedreason;
898: PetscFunctionReturn(PETSC_SUCCESS);
899: }
901: /*@
902: PCGetFailedReasonRank - Gets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail on this MPI rank
904: Not Collective
906: Input Parameter:
907: . pc - the preconditioner context
909: Output Parameter:
910: . reason - the reason it failed
912: Level: advanced
914: Note:
915: Different ranks may have different reasons or no reason, see `PCGetFailedReason()`
917: .seealso: `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCGetFailedReason()`, `PCSetFailedReason()`
918: @*/
919: PetscErrorCode PCGetFailedReasonRank(PC pc, PCFailedReason *reason)
920: {
921: PetscFunctionBegin;
922: if (pc->setupcalled < 0) *reason = (PCFailedReason)pc->setupcalled;
923: else *reason = pc->failedreason;
924: PetscFunctionReturn(PETSC_SUCCESS);
925: }
927: /* Next line needed to deactivate KSP_Solve logging */
928: #include <petsc/private/kspimpl.h>
930: /*
931: a setupcall of 0 indicates never setup,
932: 1 indicates has been previously setup
933: -1 indicates a PCSetUp() was attempted and failed
934: */
935: /*@
936: PCSetUp - Prepares for the use of a preconditioner.
938: Collective
940: Input Parameter:
941: . pc - the preconditioner context
943: Level: developer
945: .seealso: `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`
946: @*/
947: PetscErrorCode PCSetUp(PC pc)
948: {
949: const char *def;
950: PetscObjectState matstate, matnonzerostate;
952: PetscFunctionBegin;
954: PetscCheck(pc->mat, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Matrix must be set first");
956: if (pc->setupcalled && pc->reusepreconditioner) {
957: PetscCall(PetscInfo(pc, "Leaving PC with identical preconditioner since reuse preconditioner is set\n"));
958: PetscFunctionReturn(PETSC_SUCCESS);
959: }
961: PetscCall(PetscObjectStateGet((PetscObject)pc->pmat, &matstate));
962: PetscCall(MatGetNonzeroState(pc->pmat, &matnonzerostate));
963: if (!pc->setupcalled) {
964: PetscCall(PetscInfo(pc, "Setting up PC for first time\n"));
965: pc->flag = DIFFERENT_NONZERO_PATTERN;
966: } else if (matstate == pc->matstate) {
967: PetscCall(PetscInfo(pc, "Leaving PC with identical preconditioner since operator is unchanged\n"));
968: PetscFunctionReturn(PETSC_SUCCESS);
969: } else {
970: if (matnonzerostate > pc->matnonzerostate) {
971: PetscCall(PetscInfo(pc, "Setting up PC with different nonzero pattern\n"));
972: pc->flag = DIFFERENT_NONZERO_PATTERN;
973: } else {
974: PetscCall(PetscInfo(pc, "Setting up PC with same nonzero pattern\n"));
975: pc->flag = SAME_NONZERO_PATTERN;
976: }
977: }
978: pc->matstate = matstate;
979: pc->matnonzerostate = matnonzerostate;
981: if (!((PetscObject)pc)->type_name) {
982: PetscCall(PCGetDefaultType_Private(pc, &def));
983: PetscCall(PCSetType(pc, def));
984: }
986: PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
987: PetscCall(MatSetErrorIfFailure(pc->mat, pc->erroriffailure));
988: PetscCall(PetscLogEventBegin(PC_SetUp, pc, 0, 0, 0));
989: if (pc->ops->setup) {
990: /* do not log solves and applications of preconditioners while constructing preconditioners; perhaps they should be logged separately from the regular solves */
991: PetscCall(KSPInitializePackage());
992: PetscCall(PetscLogEventDeactivatePush(KSP_Solve));
993: PetscCall(PetscLogEventDeactivatePush(PC_Apply));
994: PetscUseTypeMethod(pc, setup);
995: PetscCall(PetscLogEventDeactivatePop(KSP_Solve));
996: PetscCall(PetscLogEventDeactivatePop(PC_Apply));
997: }
998: PetscCall(PetscLogEventEnd(PC_SetUp, pc, 0, 0, 0));
999: if (!pc->setupcalled) pc->setupcalled = 1;
1000: PetscFunctionReturn(PETSC_SUCCESS);
1001: }
1003: /*@
1004: PCSetUpOnBlocks - Sets up the preconditioner for each block in
1005: the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
1006: methods.
1008: Collective
1010: Input Parameter:
1011: . pc - the preconditioner context
1013: Level: developer
1015: Note:
1016: For nested preconditioners such as `PCBJACOBI` `PCSetUp()` is not called on each sub-`KSP` when `PCSetUp()` is
1017: called on the outer `PC`, this routine ensures it is called.
1019: .seealso: `PC`, `PCSetUp()`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCSetUp()`
1020: @*/
1021: PetscErrorCode PCSetUpOnBlocks(PC pc)
1022: {
1023: PetscFunctionBegin;
1025: if (!pc->ops->setuponblocks) PetscFunctionReturn(PETSC_SUCCESS);
1026: PetscCall(PetscLogEventBegin(PC_SetUpOnBlocks, pc, 0, 0, 0));
1027: PetscUseTypeMethod(pc, setuponblocks);
1028: PetscCall(PetscLogEventEnd(PC_SetUpOnBlocks, pc, 0, 0, 0));
1029: PetscFunctionReturn(PETSC_SUCCESS);
1030: }
1032: /*@C
1033: PCSetModifySubMatrices - Sets a user-defined routine for modifying the
1034: submatrices that arise within certain subdomain-based preconditioners.
1035: The basic submatrices are extracted from the preconditioner matrix as
1036: usual; the user can then alter these (for example, to set different boundary
1037: conditions for each submatrix) before they are used for the local solves.
1039: Logically Collective
1041: Input Parameters:
1042: + pc - the preconditioner context
1043: . func - routine for modifying the submatrices
1044: - ctx - optional user-defined context (may be null)
1046: Calling sequence of `func`:
1047: $ PetscErrorCode func(PC pc, PetscInt nsub, IS *row, IS *col, Mat *submat, void *ctx);
1048: + pc - the preconditioner context
1049: . row - an array of index sets that contain the global row numbers
1050: that comprise each local submatrix
1051: . col - an array of index sets that contain the global column numbers
1052: that comprise each local submatrix
1053: . submat - array of local submatrices
1054: - ctx - optional user-defined context for private data for the
1055: user-defined func routine (may be null)
1057: Level: advanced
1059: Notes:
1060: `PCSetModifySubMatrices()` MUST be called before `KSPSetUp()` and
1061: `KSPSolve()`.
1063: A routine set by `PCSetModifySubMatrices()` is currently called within
1064: the block Jacobi (`PCBJACOBI`) and additive Schwarz (`PCASM`)
1065: preconditioners. All other preconditioners ignore this routine.
1067: .seealso: `PC`, `PCBJACOBI`, `PCASM`, `PCModifySubMatrices()`
1068: @*/
1069: PetscErrorCode PCSetModifySubMatrices(PC pc, PetscErrorCode (*func)(PC, PetscInt, const IS[], const IS[], Mat[], void *), void *ctx)
1070: {
1071: PetscFunctionBegin;
1073: pc->modifysubmatrices = func;
1074: pc->modifysubmatricesP = ctx;
1075: PetscFunctionReturn(PETSC_SUCCESS);
1076: }
1078: /*@C
1079: PCModifySubMatrices - Calls an optional user-defined routine within
1080: certain preconditioners if one has been set with `PCSetModifySubMatrices()`.
1082: Collective
1084: Input Parameters:
1085: + pc - the preconditioner context
1086: . nsub - the number of local submatrices
1087: . row - an array of index sets that contain the global row numbers
1088: that comprise each local submatrix
1089: . col - an array of index sets that contain the global column numbers
1090: that comprise each local submatrix
1091: . submat - array of local submatrices
1092: - ctx - optional user-defined context for private data for the
1093: user-defined routine (may be null)
1095: Output Parameter:
1096: . submat - array of local submatrices (the entries of which may
1097: have been modified)
1099: Level: developer
1101: Notes:
1102: The user should NOT generally call this routine, as it will
1103: automatically be called within certain preconditioners (currently
1104: block Jacobi, additive Schwarz) if set.
1106: The basic submatrices are extracted from the preconditioner matrix
1107: as usual; the user can then alter these (for example, to set different
1108: boundary conditions for each submatrix) before they are used for the
1109: local solves.
1111: .seealso: `PC`, `PCSetModifySubMatrices()`
1112: @*/
1113: PetscErrorCode PCModifySubMatrices(PC pc, PetscInt nsub, const IS row[], const IS col[], Mat submat[], void *ctx)
1114: {
1115: PetscFunctionBegin;
1117: if (!pc->modifysubmatrices) PetscFunctionReturn(PETSC_SUCCESS);
1118: PetscCall(PetscLogEventBegin(PC_ModifySubMatrices, pc, 0, 0, 0));
1119: PetscCall((*pc->modifysubmatrices)(pc, nsub, row, col, submat, ctx));
1120: PetscCall(PetscLogEventEnd(PC_ModifySubMatrices, pc, 0, 0, 0));
1121: PetscFunctionReturn(PETSC_SUCCESS);
1122: }
1124: /*@
1125: PCSetOperators - Sets the matrix associated with the linear system and
1126: a (possibly) different one associated with the preconditioner.
1128: Logically Collective
1130: Input Parameters:
1131: + pc - the preconditioner context
1132: . Amat - the matrix that defines the linear system
1133: - Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
1135: Level: intermediate
1137: Notes:
1138: Passing a `NULL` for `Amat` or `Pmat` removes the matrix that is currently used.
1140: If you wish to replace either `Amat` or `Pmat` but leave the other one untouched then
1141: first call `KSPGetOperators()` to get the one you wish to keep, call `PetscObjectReference()`
1142: on it and then pass it back in in your call to `KSPSetOperators()`.
1144: More Notes about Repeated Solution of Linear Systems:
1145: PETSc does NOT reset the matrix entries of either `Amat` or `Pmat`
1146: to zero after a linear solve; the user is completely responsible for
1147: matrix assembly. See the routine `MatZeroEntries()` if desiring to
1148: zero all elements of a matrix.
1150: .seealso: `PC`, `PCGetOperators()`, `MatZeroEntries()`
1151: @*/
1152: PetscErrorCode PCSetOperators(PC pc, Mat Amat, Mat Pmat)
1153: {
1154: PetscInt m1, n1, m2, n2;
1156: PetscFunctionBegin;
1160: if (Amat) PetscCheckSameComm(pc, 1, Amat, 2);
1161: if (Pmat) PetscCheckSameComm(pc, 1, Pmat, 3);
1162: if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
1163: PetscCall(MatGetLocalSize(Amat, &m1, &n1));
1164: PetscCall(MatGetLocalSize(pc->mat, &m2, &n2));
1165: 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);
1166: PetscCall(MatGetLocalSize(Pmat, &m1, &n1));
1167: PetscCall(MatGetLocalSize(pc->pmat, &m2, &n2));
1168: 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);
1169: }
1171: if (Pmat != pc->pmat) {
1172: /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
1173: pc->matnonzerostate = -1;
1174: pc->matstate = -1;
1175: }
1177: /* reference first in case the matrices are the same */
1178: if (Amat) PetscCall(PetscObjectReference((PetscObject)Amat));
1179: PetscCall(MatDestroy(&pc->mat));
1180: if (Pmat) PetscCall(PetscObjectReference((PetscObject)Pmat));
1181: PetscCall(MatDestroy(&pc->pmat));
1182: pc->mat = Amat;
1183: pc->pmat = Pmat;
1184: PetscFunctionReturn(PETSC_SUCCESS);
1185: }
1187: /*@
1188: PCSetReusePreconditioner - reuse the current preconditioner even if the operator in the preconditioner has changed.
1190: Logically Collective
1192: Input Parameters:
1193: + pc - the preconditioner context
1194: - flag - `PETSC_TRUE` do not compute a new preconditioner, `PETSC_FALSE` do compute a new preconditioner
1196: Level: intermediate
1198: Note:
1199: Normally if a matrix inside a `PC` changes the `PC` automatically updates itself using information from the changed matrix. This option
1200: prevents this.
1202: .seealso: `PC`, `PCGetOperators()`, `MatZeroEntries()`, `PCGetReusePreconditioner()`, `KSPSetReusePreconditioner()`
1203: @*/
1204: PetscErrorCode PCSetReusePreconditioner(PC pc, PetscBool flag)
1205: {
1206: PetscFunctionBegin;
1209: pc->reusepreconditioner = flag;
1210: PetscFunctionReturn(PETSC_SUCCESS);
1211: }
1213: /*@
1214: PCGetReusePreconditioner - Determines if the `PC` reuses the current preconditioner even if the operator in the preconditioner has changed.
1216: Not Collective
1218: Input Parameter:
1219: . pc - the preconditioner context
1221: Output Parameter:
1222: . flag - `PETSC_TRUE` do not compute a new preconditioner, `PETSC_FALSE` do compute a new preconditioner
1224: Level: intermediate
1226: .seealso: `PC`, `PCGetOperators()`, `MatZeroEntries()`, `PCSetReusePreconditioner()`
1227: @*/
1228: PetscErrorCode PCGetReusePreconditioner(PC pc, PetscBool *flag)
1229: {
1230: PetscFunctionBegin;
1233: *flag = pc->reusepreconditioner;
1234: PetscFunctionReturn(PETSC_SUCCESS);
1235: }
1237: /*@
1238: PCGetOperators - Gets the matrix associated with the linear system and
1239: possibly a different one associated with the preconditioner.
1241: Not Collective, though parallel mats are returned if pc is parallel
1243: Input Parameter:
1244: . pc - the preconditioner context
1246: Output Parameters:
1247: + Amat - the matrix defining the linear system
1248: - Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.
1250: Level: intermediate
1252: Note:
1253: Does not increase the reference count of the matrices, so you should not destroy them
1255: Alternative usage: If the operators have NOT been set with `KSPSetOperators()`/`PCSetOperators()` then the operators
1256: are created in `PC` and returned to the user. In this case, if both operators
1257: mat and pmat are requested, two DIFFERENT operators will be returned. If
1258: only one is requested both operators in the PC will be the same (i.e. as
1259: if one had called `KSPSetOperators()`/`PCSetOperators()` with the same argument for both Mats).
1260: The user must set the sizes of the returned matrices and their type etc just
1261: as if the user created them with `MatCreate()`. For example,
1263: .vb
1264: KSP/PCGetOperators(ksp/pc,&Amat,NULL); is equivalent to
1265: set size, type, etc of Amat
1267: MatCreate(comm,&mat);
1268: KSP/PCSetOperators(ksp/pc,Amat,Amat);
1269: PetscObjectDereference((PetscObject)mat);
1270: set size, type, etc of Amat
1271: .ve
1273: and
1275: .vb
1276: KSP/PCGetOperators(ksp/pc,&Amat,&Pmat); is equivalent to
1277: set size, type, etc of Amat and Pmat
1279: MatCreate(comm,&Amat);
1280: MatCreate(comm,&Pmat);
1281: KSP/PCSetOperators(ksp/pc,Amat,Pmat);
1282: PetscObjectDereference((PetscObject)Amat);
1283: PetscObjectDereference((PetscObject)Pmat);
1284: set size, type, etc of Amat and Pmat
1285: .ve
1287: The rationale for this support is so that when creating a `TS`, `SNES`, or `KSP` the hierarchy
1288: of underlying objects (i.e. `SNES`, `KSP`, `PC`, `Mat`) and their lifespans can be completely
1289: managed by the top most level object (i.e. the `TS`, `SNES`, or `KSP`). Another way to look
1290: at this is when you create a `SNES` you do not NEED to create a `KSP` and attach it to
1291: the `SNES` object (the `SNES` object manages it for you). Similarly when you create a KSP
1292: you do not need to attach a `PC` to it (the `KSP` object manages the `PC` object for you).
1293: Thus, why should YOU have to create the `Mat` and attach it to the `SNES`/`KSP`/`PC`, when
1294: it can be created for you?
1296: .seealso: `PC`, `PCSetOperators()`, `KSPGetOperators()`, `KSPSetOperators()`, `PCGetOperatorsSet()`
1297: @*/
1298: PetscErrorCode PCGetOperators(PC pc, Mat *Amat, Mat *Pmat)
1299: {
1300: PetscFunctionBegin;
1302: if (Amat) {
1303: if (!pc->mat) {
1304: if (pc->pmat && !Pmat) { /* Pmat has been set, but user did not request it, so use for Amat */
1305: pc->mat = pc->pmat;
1306: PetscCall(PetscObjectReference((PetscObject)pc->mat));
1307: } else { /* both Amat and Pmat are empty */
1308: PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &pc->mat));
1309: if (!Pmat) { /* user did NOT request Pmat, so make same as Amat */
1310: pc->pmat = pc->mat;
1311: PetscCall(PetscObjectReference((PetscObject)pc->pmat));
1312: }
1313: }
1314: }
1315: *Amat = pc->mat;
1316: }
1317: if (Pmat) {
1318: if (!pc->pmat) {
1319: if (pc->mat && !Amat) { /* Amat has been set but was not requested, so use for pmat */
1320: pc->pmat = pc->mat;
1321: PetscCall(PetscObjectReference((PetscObject)pc->pmat));
1322: } else {
1323: PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &pc->pmat));
1324: if (!Amat) { /* user did NOT request Amat, so make same as Pmat */
1325: pc->mat = pc->pmat;
1326: PetscCall(PetscObjectReference((PetscObject)pc->mat));
1327: }
1328: }
1329: }
1330: *Pmat = pc->pmat;
1331: }
1332: PetscFunctionReturn(PETSC_SUCCESS);
1333: }
1335: /*@C
1336: PCGetOperatorsSet - Determines if the matrix associated with the linear system and
1337: possibly a different one associated with the preconditioner have been set in the `PC`.
1339: Not Collective, though the results on all processes should be the same
1341: Input Parameter:
1342: . pc - the preconditioner context
1344: Output Parameters:
1345: + mat - the matrix associated with the linear system was set
1346: - pmat - matrix associated with the preconditioner was set, usually the same
1348: Level: intermediate
1350: .seealso: `PC`, `PCSetOperators()`, `KSPGetOperators()`, `KSPSetOperators()`, `PCGetOperators()`
1351: @*/
1352: PetscErrorCode PCGetOperatorsSet(PC pc, PetscBool *mat, PetscBool *pmat)
1353: {
1354: PetscFunctionBegin;
1356: if (mat) *mat = (pc->mat) ? PETSC_TRUE : PETSC_FALSE;
1357: if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1358: PetscFunctionReturn(PETSC_SUCCESS);
1359: }
1361: /*@
1362: PCFactorGetMatrix - Gets the factored matrix from the
1363: preconditioner context. This routine is valid only for the `PCLU`,
1364: `PCILU`, `PCCHOLESKY`, and `PCICC` methods.
1366: Not Collective though mat is parallel if pc is parallel
1368: Input Parameter:
1369: . pc - the preconditioner context
1371: Output parameters:
1372: . mat - the factored matrix
1374: Level: advanced
1376: Note:
1377: Does not increase the reference count for the matrix so DO NOT destroy it
1379: .seealso: `PC`, `PCLU`, `PCILU`, `PCCHOLESKY`, `PCICC`
1380: @*/
1381: PetscErrorCode PCFactorGetMatrix(PC pc, Mat *mat)
1382: {
1383: PetscFunctionBegin;
1386: PetscUseTypeMethod(pc, getfactoredmatrix, mat);
1387: PetscFunctionReturn(PETSC_SUCCESS);
1388: }
1390: /*@C
1391: PCSetOptionsPrefix - Sets the prefix used for searching for all
1392: `PC` options in the database.
1394: Logically Collective
1396: Input Parameters:
1397: + pc - the preconditioner context
1398: - prefix - the prefix string to prepend to all `PC` option requests
1400: Note:
1401: A hyphen (-) must NOT be given at the beginning of the prefix name.
1402: The first character of all runtime options is AUTOMATICALLY the
1403: hyphen.
1405: Level: advanced
1407: .seealso: `PC`, `PCSetFromOptions`, `PCAppendOptionsPrefix()`, `PCGetOptionsPrefix()`
1408: @*/
1409: PetscErrorCode PCSetOptionsPrefix(PC pc, const char prefix[])
1410: {
1411: PetscFunctionBegin;
1413: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc, prefix));
1414: PetscFunctionReturn(PETSC_SUCCESS);
1415: }
1417: /*@C
1418: PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1419: `PC` options in the database.
1421: Logically Collective
1423: Input Parameters:
1424: + pc - the preconditioner context
1425: - prefix - the prefix string to prepend to all `PC` option requests
1427: Note:
1428: A hyphen (-) must NOT be given at the beginning of the prefix name.
1429: The first character of all runtime options is AUTOMATICALLY the
1430: hyphen.
1432: Level: advanced
1434: .seealso: `PC`, `PCSetFromOptions`, `PCSetOptionsPrefix()`, `PCGetOptionsPrefix()`
1435: @*/
1436: PetscErrorCode PCAppendOptionsPrefix(PC pc, const char prefix[])
1437: {
1438: PetscFunctionBegin;
1440: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)pc, prefix));
1441: PetscFunctionReturn(PETSC_SUCCESS);
1442: }
1444: /*@C
1445: PCGetOptionsPrefix - Gets the prefix used for searching for all
1446: PC options in the database.
1448: Not Collective
1450: Input Parameter:
1451: . pc - the preconditioner context
1453: Output Parameter:
1454: . prefix - pointer to the prefix string used, is returned
1456: Fortran Note:
1457: The user should pass in a string 'prefix' of
1458: sufficient length to hold the prefix.
1460: Level: advanced
1462: .seealso: `PC`, `PCSetFromOptions`, `PCSetOptionsPrefix()`, `PCAppendOptionsPrefix()`
1463: @*/
1464: PetscErrorCode PCGetOptionsPrefix(PC pc, const char *prefix[])
1465: {
1466: PetscFunctionBegin;
1469: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, prefix));
1470: PetscFunctionReturn(PETSC_SUCCESS);
1471: }
1473: /*
1474: Indicates the right hand side will be changed by KSPSolve(), this occurs for a few
1475: preconditioners including BDDC and Eisentat that transform the equations before applying
1476: the Krylov methods
1477: */
1478: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC pc, PetscBool *change)
1479: {
1480: PetscFunctionBegin;
1483: *change = PETSC_FALSE;
1484: PetscTryMethod(pc, "PCPreSolveChangeRHS_C", (PC, PetscBool *), (pc, change));
1485: PetscFunctionReturn(PETSC_SUCCESS);
1486: }
1488: /*@
1489: PCPreSolve - Optional pre-solve phase, intended for any
1490: preconditioner-specific actions that must be performed before
1491: the iterative solve itself.
1493: Collective
1495: Input Parameters:
1496: + pc - the preconditioner context
1497: - ksp - the Krylov subspace context
1499: Level: developer
1501: Sample of Usage:
1502: .vb
1503: PCPreSolve(pc,ksp);
1504: KSPSolve(ksp,b,x);
1505: PCPostSolve(pc,ksp);
1506: .ve
1508: Notes:
1509: The pre-solve phase is distinct from the `PCSetUp()` phase.
1511: `KSPSolve()` calls this directly, so is rarely called by the user.
1513: .seealso: `PC`, `PCPostSolve()`
1514: @*/
1515: PetscErrorCode PCPreSolve(PC pc, KSP ksp)
1516: {
1517: Vec x, rhs;
1519: PetscFunctionBegin;
1522: pc->presolvedone++;
1523: PetscCheck(pc->presolvedone <= 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot embed PCPreSolve() more than twice");
1524: PetscCall(KSPGetSolution(ksp, &x));
1525: PetscCall(KSPGetRhs(ksp, &rhs));
1527: if (pc->ops->presolve) PetscUseTypeMethod(pc, presolve, ksp, rhs, x);
1528: else if (pc->presolve) PetscCall((pc->presolve)(pc, ksp));
1529: PetscFunctionReturn(PETSC_SUCCESS);
1530: }
1532: /*@C
1533: PCSetPreSolve - Sets function used by `PCPreSolve()` which is intended for any
1534: preconditioner-specific actions that must be performed before
1535: the iterative solve itself.
1537: Logically Collective
1539: Input Parameters:
1540: + pc - the preconditioner object
1541: - presolve - the function to call before the solve
1543: Calling sequence of `presolve`:
1544: $ PetscErrorCode presolve(PC pc, KSP ksp)
1545: + pc - the `PC` context
1546: - ksp - the `KSP` context
1548: Level: developer
1550: .seealso: `PC`, `PCSetUp()`, `PCPreSolve()`
1551: @*/
1552: PetscErrorCode PCSetPreSolve(PC pc, PetscErrorCode (*presolve)(PC, KSP))
1553: {
1554: PetscFunctionBegin;
1556: pc->presolve = presolve;
1557: PetscFunctionReturn(PETSC_SUCCESS);
1558: }
1560: /*@
1561: PCPostSolve - Optional post-solve phase, intended for any
1562: preconditioner-specific actions that must be performed after
1563: the iterative solve itself.
1565: Collective
1567: Input Parameters:
1568: + pc - the preconditioner context
1569: - ksp - the Krylov subspace context
1571: Sample of Usage:
1572: .vb
1573: PCPreSolve(pc,ksp);
1574: KSPSolve(ksp,b,x);
1575: PCPostSolve(pc,ksp);
1576: .ve
1578: Note:
1579: `KSPSolve()` calls this routine directly, so it is rarely called by the user.
1581: Level: developer
1583: .seealso: `PC`, `PCSetPostSolve()`, `PCSetPresolve()`, `PCPreSolve()`, `KSPSolve()`
1584: @*/
1585: PetscErrorCode PCPostSolve(PC pc, KSP ksp)
1586: {
1587: Vec x, rhs;
1589: PetscFunctionBegin;
1592: pc->presolvedone--;
1593: PetscCall(KSPGetSolution(ksp, &x));
1594: PetscCall(KSPGetRhs(ksp, &rhs));
1595: PetscTryTypeMethod(pc, postsolve, ksp, rhs, x);
1596: PetscFunctionReturn(PETSC_SUCCESS);
1597: }
1599: /*@C
1600: PCLoad - Loads a `PC` that has been stored in binary with `PCView()`.
1602: Collective
1604: Input Parameters:
1605: + newdm - the newly loaded `PC`, this needs to have been created with `PCCreate()` or
1606: some related function before a call to `PCLoad()`.
1607: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()`
1609: Level: intermediate
1611: Note:
1612: The type is determined by the data in the file, any type set into the PC before this call is ignored.
1614: .seealso: `PC`, `PetscViewerBinaryOpen()`, `PCView()`, `MatLoad()`, `VecLoad()`
1615: @*/
1616: PetscErrorCode PCLoad(PC newdm, PetscViewer viewer)
1617: {
1618: PetscBool isbinary;
1619: PetscInt classid;
1620: char type[256];
1622: PetscFunctionBegin;
1625: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1626: PetscCheck(isbinary, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1628: PetscCall(PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT));
1629: PetscCheck(classid == PC_FILE_CLASSID, PetscObjectComm((PetscObject)newdm), PETSC_ERR_ARG_WRONG, "Not PC next in file");
1630: PetscCall(PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR));
1631: PetscCall(PCSetType(newdm, type));
1632: PetscTryTypeMethod(newdm, load, viewer);
1633: PetscFunctionReturn(PETSC_SUCCESS);
1634: }
1636: #include <petscdraw.h>
1637: #if defined(PETSC_HAVE_SAWS)
1638: #include <petscviewersaws.h>
1639: #endif
1641: /*@C
1642: PCViewFromOptions - View from the `PC` based on options in the database
1644: Collective
1646: Input Parameters:
1647: + A - the `PC` context
1648: . obj - Optional object that provides the options prefix
1649: - name - command line option
1651: Level: intermediate
1653: .seealso: `PC`, `PCView`, `PetscObjectViewFromOptions()`, `PCCreate()`
1654: @*/
1655: PetscErrorCode PCViewFromOptions(PC A, PetscObject obj, const char name[])
1656: {
1657: PetscFunctionBegin;
1659: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
1660: PetscFunctionReturn(PETSC_SUCCESS);
1661: }
1663: /*@C
1664: PCView - Prints information about the `PC`
1666: Collective
1668: Input Parameters:
1669: + PC - the `PC` context
1670: - viewer - optional visualization context
1672: Level: developer
1674: Notes:
1675: The available visualization contexts include
1676: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
1677: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
1678: output where only the first processor opens
1679: the file. All other processors send their
1680: data to the first processor to print.
1682: The user can open an alternative visualization contexts with
1683: `PetscViewerASCIIOpen()` (output to a specified file).
1685: .seealso: `PC`, `PetscViewer`, `KSPView()`, `PetscViewerASCIIOpen()`
1686: @*/
1687: PetscErrorCode PCView(PC pc, PetscViewer viewer)
1688: {
1689: PCType cstr;
1690: PetscBool iascii, isstring, isbinary, isdraw;
1691: #if defined(PETSC_HAVE_SAWS)
1692: PetscBool issaws;
1693: #endif
1695: PetscFunctionBegin;
1697: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc), &viewer));
1699: PetscCheckSameComm(pc, 1, viewer, 2);
1701: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1702: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
1703: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1704: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1705: #if defined(PETSC_HAVE_SAWS)
1706: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws));
1707: #endif
1709: if (iascii) {
1710: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)pc, viewer));
1711: if (!pc->setupcalled) PetscCall(PetscViewerASCIIPrintf(viewer, " PC has not been set up so information may be incomplete\n"));
1712: PetscCall(PetscViewerASCIIPushTab(viewer));
1713: PetscTryTypeMethod(pc, view, viewer);
1714: PetscCall(PetscViewerASCIIPopTab(viewer));
1715: if (pc->mat) {
1716: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
1717: if (pc->pmat == pc->mat) {
1718: PetscCall(PetscViewerASCIIPrintf(viewer, " linear system matrix = precond matrix:\n"));
1719: PetscCall(PetscViewerASCIIPushTab(viewer));
1720: PetscCall(MatView(pc->mat, viewer));
1721: PetscCall(PetscViewerASCIIPopTab(viewer));
1722: } else {
1723: if (pc->pmat) {
1724: PetscCall(PetscViewerASCIIPrintf(viewer, " linear system matrix followed by preconditioner matrix:\n"));
1725: } else {
1726: PetscCall(PetscViewerASCIIPrintf(viewer, " linear system matrix:\n"));
1727: }
1728: PetscCall(PetscViewerASCIIPushTab(viewer));
1729: PetscCall(MatView(pc->mat, viewer));
1730: if (pc->pmat) PetscCall(MatView(pc->pmat, viewer));
1731: PetscCall(PetscViewerASCIIPopTab(viewer));
1732: }
1733: PetscCall(PetscViewerPopFormat(viewer));
1734: }
1735: } else if (isstring) {
1736: PetscCall(PCGetType(pc, &cstr));
1737: PetscCall(PetscViewerStringSPrintf(viewer, " PCType: %-7.7s", cstr));
1738: PetscTryTypeMethod(pc, view, viewer);
1739: if (pc->mat) PetscCall(MatView(pc->mat, viewer));
1740: if (pc->pmat && pc->pmat != pc->mat) PetscCall(MatView(pc->pmat, viewer));
1741: } else if (isbinary) {
1742: PetscInt classid = PC_FILE_CLASSID;
1743: MPI_Comm comm;
1744: PetscMPIInt rank;
1745: char type[256];
1747: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
1748: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1749: if (rank == 0) {
1750: PetscCall(PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT));
1751: PetscCall(PetscStrncpy(type, ((PetscObject)pc)->type_name, 256));
1752: PetscCall(PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR));
1753: }
1754: PetscTryTypeMethod(pc, view, viewer);
1755: } else if (isdraw) {
1756: PetscDraw draw;
1757: char str[25];
1758: PetscReal x, y, bottom, h;
1759: PetscInt n;
1761: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1762: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
1763: if (pc->mat) {
1764: PetscCall(MatGetSize(pc->mat, &n, NULL));
1765: PetscCall(PetscSNPrintf(str, 25, "PC: %s (%" PetscInt_FMT ")", ((PetscObject)pc)->type_name, n));
1766: } else {
1767: PetscCall(PetscSNPrintf(str, 25, "PC: %s", ((PetscObject)pc)->type_name));
1768: }
1769: PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h));
1770: bottom = y - h;
1771: PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
1772: PetscTryTypeMethod(pc, view, viewer);
1773: PetscCall(PetscDrawPopCurrentPoint(draw));
1774: #if defined(PETSC_HAVE_SAWS)
1775: } else if (issaws) {
1776: PetscMPIInt rank;
1778: PetscCall(PetscObjectName((PetscObject)pc));
1779: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1780: if (!((PetscObject)pc)->amsmem && rank == 0) PetscCall(PetscObjectViewSAWs((PetscObject)pc, viewer));
1781: if (pc->mat) PetscCall(MatView(pc->mat, viewer));
1782: if (pc->pmat && pc->pmat != pc->mat) PetscCall(MatView(pc->pmat, viewer));
1783: #endif
1784: }
1785: PetscFunctionReturn(PETSC_SUCCESS);
1786: }
1788: /*@C
1789: PCRegister - Adds a method to the preconditioner package.
1791: Not collective
1793: Input Parameters:
1794: + sname - name of a new user-defined solver
1795: - function - routine to create method context
1797: Sample usage:
1798: .vb
1799: PCRegister("my_solver", MySolverCreate);
1800: .ve
1802: Then, your solver can be chosen with the procedural interface via
1803: $ PCSetType(pc, "my_solver")
1804: or at runtime via the option
1805: $ -pc_type my_solver
1807: Level: advanced
1809: Note:
1810: `PCRegister()` may be called multiple times to add several user-defined preconditioners.
1812: .seealso: `PC`, `PCType`, `PCRegisterAll()`
1813: @*/
1814: PetscErrorCode PCRegister(const char sname[], PetscErrorCode (*function)(PC))
1815: {
1816: PetscFunctionBegin;
1817: PetscCall(PCInitializePackage());
1818: PetscCall(PetscFunctionListAdd(&PCList, sname, function));
1819: PetscFunctionReturn(PETSC_SUCCESS);
1820: }
1822: static PetscErrorCode MatMult_PC(Mat A, Vec X, Vec Y)
1823: {
1824: PC pc;
1826: PetscFunctionBegin;
1827: PetscCall(MatShellGetContext(A, &pc));
1828: PetscCall(PCApply(pc, X, Y));
1829: PetscFunctionReturn(PETSC_SUCCESS);
1830: }
1832: /*@
1833: PCComputeOperator - Computes the explicit preconditioned operator.
1835: Collective
1837: Input Parameters:
1838: + pc - the preconditioner object
1839: - mattype - the matrix type to be used for the operator
1841: Output Parameter:
1842: . mat - the explicit preconditioned operator
1844: Level: advanced
1846: Note:
1847: This computation is done by applying the operators to columns of the identity matrix.
1848: This routine is costly in general, and is recommended for use only with relatively small systems.
1849: Currently, this routine uses a dense matrix format when mattype == NULL
1851: .seealso: `PC`, `KSPComputeOperator()`, `MatType`
1852: @*/
1853: PetscErrorCode PCComputeOperator(PC pc, MatType mattype, Mat *mat)
1854: {
1855: PetscInt N, M, m, n;
1856: Mat A, Apc;
1858: PetscFunctionBegin;
1861: PetscCall(PCGetOperators(pc, &A, NULL));
1862: PetscCall(MatGetLocalSize(A, &m, &n));
1863: PetscCall(MatGetSize(A, &M, &N));
1864: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)pc), m, n, M, N, pc, &Apc));
1865: PetscCall(MatShellSetOperation(Apc, MATOP_MULT, (void (*)(void))MatMult_PC));
1866: PetscCall(MatComputeOperator(Apc, mattype, mat));
1867: PetscCall(MatDestroy(&Apc));
1868: PetscFunctionReturn(PETSC_SUCCESS);
1869: }
1871: /*@
1872: PCSetCoordinates - sets the coordinates of all the nodes on the local process
1874: Collective
1876: Input Parameters:
1877: + pc - the solver context
1878: . dim - the dimension of the coordinates 1, 2, or 3
1879: . nloc - the blocked size of the coordinates array
1880: - coords - the coordinates array
1882: Level: intermediate
1884: Note:
1885: `coords` is an array of the dim coordinates for the nodes on
1886: the local processor, of size `dim`*`nloc`.
1887: If there are 108 equation on a processor
1888: for a displacement finite element discretization of elasticity (so
1889: that there are nloc = 36 = 108/3 nodes) then the array must have 108
1890: double precision values (ie, 3 * 36). These x y z coordinates
1891: should be ordered for nodes 0 to N-1 like so: [ 0.x, 0.y, 0.z, 1.x,
1892: ... , N-1.z ].
1894: .seealso: `PC`, `MatSetNearNullSpace()`
1895: @*/
1896: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[])
1897: {
1898: PetscFunctionBegin;
1901: PetscTryMethod(pc, "PCSetCoordinates_C", (PC, PetscInt, PetscInt, PetscReal *), (pc, dim, nloc, coords));
1902: PetscFunctionReturn(PETSC_SUCCESS);
1903: }
1905: /*@
1906: PCGetInterpolations - Gets interpolation matrices for all levels (except level 0)
1908: Logically Collective
1910: Input Parameter:
1911: . pc - the precondition context
1913: Output Parameters:
1914: + num_levels - the number of levels
1915: - interpolations - the interpolation matrices (size of num_levels-1)
1917: Level: advanced
1919: Developer Note:
1920: Why is this here instead of in `PCMG` etc?
1922: .seealso: `PC`, `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetInterpolation()`, `PCGetCoarseOperators()`
1923: @*/
1924: PetscErrorCode PCGetInterpolations(PC pc, PetscInt *num_levels, Mat *interpolations[])
1925: {
1926: PetscFunctionBegin;
1930: PetscUseMethod(pc, "PCGetInterpolations_C", (PC, PetscInt *, Mat *[]), (pc, num_levels, interpolations));
1931: PetscFunctionReturn(PETSC_SUCCESS);
1932: }
1934: /*@
1935: PCGetCoarseOperators - Gets coarse operator matrices for all levels (except the finest level)
1937: Logically Collective
1939: Input Parameter:
1940: . pc - the precondition context
1942: Output Parameters:
1943: + num_levels - the number of levels
1944: - coarseOperators - the coarse operator matrices (size of num_levels-1)
1946: Level: advanced
1948: Developer Note:
1949: Why is this here instead of in `PCMG` etc?
1951: .seealso: `PC`, `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetRScale()`, `PCMGGetInterpolation()`, `PCGetInterpolations()`
1952: @*/
1953: PetscErrorCode PCGetCoarseOperators(PC pc, PetscInt *num_levels, Mat *coarseOperators[])
1954: {
1955: PetscFunctionBegin;
1959: PetscUseMethod(pc, "PCGetCoarseOperators_C", (PC, PetscInt *, Mat *[]), (pc, num_levels, coarseOperators));
1960: PetscFunctionReturn(PETSC_SUCCESS);
1961: }