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 the PC
69: .seealso: `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: `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
133: .vb
134: D M A D^{-1} y = D M b for left preconditioning or
135: D A M D^{-1} z = D b for right preconditioning
136: .ve
138: .seealso: `PC`, `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCDiagonalScaleRight()`, `PCSetDiagonalScale()`
139: @*/
140: PetscErrorCode PCGetDiagonalScale(PC pc, PetscBool *flag)
141: {
142: PetscFunctionBegin;
144: PetscAssertPointer(flag, 2);
145: *flag = pc->diagonalscale;
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: /*@
150: PCSetDiagonalScale - Indicates the left scaling to use to apply an additional left and right
151: scaling as needed by certain time-stepping codes.
153: Logically Collective
155: Input Parameters:
156: + pc - the preconditioner context
157: - s - scaling vector
159: Level: intermediate
161: Notes:
162: The system solved via the Krylov method is
163: .vb
164: D M A D^{-1} y = D M b for left preconditioning or
165: D A M D^{-1} z = D b for right preconditioning
166: .ve
168: `PCDiagonalScaleLeft()` scales a vector by D. `PCDiagonalScaleRight()` scales a vector by D^{-1}.
170: .seealso: `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCDiagonalScaleRight()`, `PCGetDiagonalScale()`
171: @*/
172: PetscErrorCode PCSetDiagonalScale(PC pc, Vec s)
173: {
174: PetscFunctionBegin;
177: pc->diagonalscale = PETSC_TRUE;
179: PetscCall(PetscObjectReference((PetscObject)s));
180: PetscCall(VecDestroy(&pc->diagonalscaleleft));
182: pc->diagonalscaleleft = s;
184: PetscCall(VecDuplicate(s, &pc->diagonalscaleright));
185: PetscCall(VecCopy(s, pc->diagonalscaleright));
186: PetscCall(VecReciprocal(pc->diagonalscaleright));
187: PetscFunctionReturn(PETSC_SUCCESS);
188: }
190: /*@
191: PCDiagonalScaleLeft - Scales a vector by the left scaling as needed by certain time-stepping codes.
193: Logically Collective
195: Input Parameters:
196: + pc - the preconditioner context
197: . in - input vector
198: - out - scaled vector (maybe the same as in)
200: Level: intermediate
202: Notes:
203: The system solved via the Krylov method is
204: .vb
205: D M A D^{-1} y = D M b for left preconditioning or
206: D A M D^{-1} z = D b for right preconditioning
207: .ve
209: `PCDiagonalScaleLeft()` scales a vector by D. `PCDiagonalScaleRight()` scales a vector by D^{-1}.
211: If diagonal scaling is turned off and in is not out then in is copied to out
213: .seealso: `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleSet()`, `PCDiagonalScaleRight()`, `PCDiagonalScale()`
214: @*/
215: PetscErrorCode PCDiagonalScaleLeft(PC pc, Vec in, Vec out)
216: {
217: PetscFunctionBegin;
221: if (pc->diagonalscale) {
222: PetscCall(VecPointwiseMult(out, pc->diagonalscaleleft, in));
223: } else if (in != out) {
224: PetscCall(VecCopy(in, out));
225: }
226: PetscFunctionReturn(PETSC_SUCCESS);
227: }
229: /*@
230: PCDiagonalScaleRight - Scales a vector by the right scaling as needed by certain time-stepping codes.
232: Logically Collective
234: Input Parameters:
235: + pc - the preconditioner context
236: . in - input vector
237: - out - scaled vector (maybe the same as in)
239: Level: intermediate
241: Notes:
242: The system solved via the Krylov method is
243: .vb
244: D M A D^{-1} y = D M b for left preconditioning or
245: D A M D^{-1} z = D b for right preconditioning
246: .ve
248: `PCDiagonalScaleLeft()` scales a vector by D. `PCDiagonalScaleRight()` scales a vector by D^{-1}.
250: If diagonal scaling is turned off and in is not out then in is copied to out
252: .seealso: `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCDiagonalScaleSet()`, `PCDiagonalScale()`
253: @*/
254: PetscErrorCode PCDiagonalScaleRight(PC pc, Vec in, Vec out)
255: {
256: PetscFunctionBegin;
260: if (pc->diagonalscale) {
261: PetscCall(VecPointwiseMult(out, pc->diagonalscaleright, in));
262: } else if (in != out) {
263: PetscCall(VecCopy(in, out));
264: }
265: PetscFunctionReturn(PETSC_SUCCESS);
266: }
268: /*@
269: PCSetUseAmat - Sets a flag to indicate that when the preconditioner needs to apply (part of) the
270: operator during the preconditioning process it applies the Amat provided to `TSSetRHSJacobian()`,
271: `TSSetIJacobian()`, `SNESSetJacobian()`, `KSPSetOperators()` or `PCSetOperators()` not the Pmat.
273: Logically Collective
275: Input Parameters:
276: + pc - the preconditioner context
277: - flg - `PETSC_TRUE` to use the Amat, `PETSC_FALSE` to use the Pmat (default is false)
279: Options Database Key:
280: . -pc_use_amat <true,false> - use the amat to apply the operator
282: Level: intermediate
284: Note:
285: For the common case in which the linear system matrix and the matrix used to construct the
286: preconditioner are identical, this routine is does nothing.
288: .seealso: `PC`, `PCGetUseAmat()`, `PCBJACOBI`, `PGMG`, `PCFIELDSPLIT`, `PCCOMPOSITE`
289: @*/
290: PetscErrorCode PCSetUseAmat(PC pc, PetscBool flg)
291: {
292: PetscFunctionBegin;
294: pc->useAmat = flg;
295: PetscFunctionReturn(PETSC_SUCCESS);
296: }
298: /*@
299: PCSetErrorIfFailure - Causes `PC` to generate an error if a FPE, for example a zero pivot, is detected.
301: Logically Collective
303: Input Parameters:
304: + pc - iterative context obtained from PCCreate()
305: - flg - `PETSC_TRUE` indicates you want the error generated
307: Level: advanced
309: Notes:
310: Normally PETSc continues if a linear solver fails due to a failed setup of a preconditioner, you can call `KSPGetConvergedReason()` after a `KSPSolve()`
311: 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
312: detected.
314: This is propagated into KSPs used by this PC, which then propagate it into PCs used by those KSPs
316: .seealso: `PC`, `KSPSetErrorIfNotConverged()`, `PCGetInitialGuessNonzero()`, `PCSetInitialGuessKnoll()`, `PCGetInitialGuessKnoll()`
317: @*/
318: PetscErrorCode PCSetErrorIfFailure(PC pc, PetscBool flg)
319: {
320: PetscFunctionBegin;
323: pc->erroriffailure = flg;
324: PetscFunctionReturn(PETSC_SUCCESS);
325: }
327: /*@
328: PCGetUseAmat - Gets a flag to indicate that when the preconditioner needs to apply (part of) the
329: operator during the preconditioning process it applies the Amat provided to `TSSetRHSJacobian()`,
330: `TSSetIJacobian()`, `SNESSetJacobian()`, `KSPSetOperators()` or `PCSetOperators()` not the Pmat.
332: Logically Collective
334: Input Parameter:
335: . pc - the preconditioner context
337: Output Parameter:
338: . flg - `PETSC_TRUE` to use the Amat, `PETSC_FALSE` to use the Pmat (default is false)
340: Level: intermediate
342: Note:
343: For the common case in which the linear system matrix and the matrix used to construct the
344: preconditioner are identical, this routine is does nothing.
346: .seealso: `PC`, `PCSetUseAmat()`, `PCBJACOBI`, `PGMG`, `PCFIELDSPLIT`, `PCCOMPOSITE`
347: @*/
348: PetscErrorCode PCGetUseAmat(PC pc, PetscBool *flg)
349: {
350: PetscFunctionBegin;
352: *flg = pc->useAmat;
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: /*@
357: PCSetKSPNestLevel - sets the amount of nesting the `KSP` that contains this `PC` has
359: Collective
361: Input Parameters:
362: + pc - the `PC`
363: - level - the nest level
365: Level: developer
367: .seealso: [](ch_ksp), `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSP`, `KSPGMRES`, `KSPType`, `KSPGetNestLevel()`, `PCGetKSPNestLevel()`, `KSPSetNestLevel()`
368: @*/
369: PetscErrorCode PCSetKSPNestLevel(PC pc, PetscInt level)
370: {
371: PetscFunctionBegin;
374: pc->kspnestlevel = level;
375: PetscFunctionReturn(PETSC_SUCCESS);
376: }
378: /*@
379: PCGetKSPNestLevel - gets the amount of nesting the `KSP` that contains this `PC` has
381: Not Collective
383: Input Parameter:
384: . pc - the `PC`
386: Output Parameter:
387: . level - the nest level
389: Level: developer
391: .seealso: [](ch_ksp), `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSP`, `KSPGMRES`, `KSPType`, `KSPSetNestLevel()`, `PCSetKSPNestLevel()`, `KSPGetNestLevel()`
392: @*/
393: PetscErrorCode PCGetKSPNestLevel(PC pc, PetscInt *level)
394: {
395: PetscFunctionBegin;
397: PetscAssertPointer(level, 2);
398: *level = pc->kspnestlevel;
399: PetscFunctionReturn(PETSC_SUCCESS);
400: }
402: /*@
403: PCCreate - Creates a preconditioner context, `PC`
405: Collective
407: Input Parameter:
408: . comm - MPI communicator
410: Output Parameter:
411: . newpc - location to put the preconditioner context
413: Level: developer
415: Note:
416: The default preconditioner for sparse matrices is `PCILU` or `PCICC` with 0 fill on one process and block Jacobi (`PCBJACOBI`) with `PCILU` or `PCICC`
417: in parallel. For dense matrices it is always `PCNONE`.
419: .seealso: `PC`, `PCSetUp()`, `PCApply()`, `PCDestroy()`
420: @*/
421: PetscErrorCode PCCreate(MPI_Comm comm, PC *newpc)
422: {
423: PC pc;
425: PetscFunctionBegin;
426: PetscAssertPointer(newpc, 2);
427: *newpc = NULL;
428: PetscCall(PCInitializePackage());
430: PetscCall(PetscHeaderCreate(pc, PC_CLASSID, "PC", "Preconditioner", "PC", comm, PCDestroy, PCView));
432: pc->mat = NULL;
433: pc->pmat = NULL;
434: pc->setupcalled = 0;
435: pc->setfromoptionscalled = 0;
436: pc->data = NULL;
437: pc->diagonalscale = PETSC_FALSE;
438: pc->diagonalscaleleft = NULL;
439: pc->diagonalscaleright = NULL;
441: pc->modifysubmatrices = NULL;
442: pc->modifysubmatricesP = NULL;
444: *newpc = pc;
445: PetscFunctionReturn(PETSC_SUCCESS);
446: }
448: /*@
449: PCApply - Applies the preconditioner to a vector.
451: Collective
453: Input Parameters:
454: + pc - the preconditioner context
455: - x - input vector
457: Output Parameter:
458: . y - output vector
460: Level: developer
462: .seealso: `PC`, `PCApplyTranspose()`, `PCApplyBAorAB()`
463: @*/
464: PetscErrorCode PCApply(PC pc, Vec x, Vec y)
465: {
466: PetscInt m, n, mv, nv;
468: PetscFunctionBegin;
472: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
473: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
474: /* use pmat to check vector sizes since for KSPLSQR the pmat may be of a different size than mat */
475: PetscCall(MatGetLocalSize(pc->pmat, &m, &n));
476: PetscCall(VecGetLocalSize(x, &mv));
477: PetscCall(VecGetLocalSize(y, &nv));
478: /* check pmat * y = x is feasible */
479: 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);
480: 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);
481: PetscCall(VecSetErrorIfLocked(y, 3));
483: PetscCall(PCSetUp(pc));
484: PetscCall(VecLockReadPush(x));
485: PetscCall(PetscLogEventBegin(PC_Apply, pc, x, y, 0));
486: PetscUseTypeMethod(pc, apply, x, y);
487: PetscCall(PetscLogEventEnd(PC_Apply, pc, x, y, 0));
488: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
489: PetscCall(VecLockReadPop(x));
490: PetscFunctionReturn(PETSC_SUCCESS);
491: }
493: /*@
494: PCMatApply - Applies the preconditioner to multiple vectors stored as a `MATDENSE`. Like `PCApply()`, Y and X must be different matrices.
496: Collective
498: Input Parameters:
499: + pc - the preconditioner context
500: - X - block of input vectors
502: Output Parameter:
503: . Y - block of output vectors
505: Level: developer
507: .seealso: `PC`, `PCApply()`, `KSPMatSolve()`
508: @*/
509: PetscErrorCode PCMatApply(PC pc, Mat X, Mat Y)
510: {
511: Mat A;
512: Vec cy, cx;
513: PetscInt m1, M1, m2, M2, n1, N1, n2, N2, m3, M3, n3, N3;
514: PetscBool match;
516: PetscFunctionBegin;
520: PetscCheckSameComm(pc, 1, X, 2);
521: PetscCheckSameComm(pc, 1, Y, 3);
522: PetscCheck(Y != X, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "Y and X must be different matrices");
523: PetscCall(PCGetOperators(pc, NULL, &A));
524: PetscCall(MatGetLocalSize(A, &m3, &n3));
525: PetscCall(MatGetLocalSize(X, &m2, &n2));
526: PetscCall(MatGetLocalSize(Y, &m1, &n1));
527: PetscCall(MatGetSize(A, &M3, &N3));
528: PetscCall(MatGetSize(X, &M2, &N2));
529: PetscCall(MatGetSize(Y, &M1, &N1));
530: 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);
531: 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);
532: 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);
533: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)Y, &match, MATSEQDENSE, MATMPIDENSE, ""));
534: PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided block of output vectors not stored in a dense Mat");
535: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)X, &match, MATSEQDENSE, MATMPIDENSE, ""));
536: PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided block of input vectors not stored in a dense Mat");
537: PetscCall(PCSetUp(pc));
538: if (pc->ops->matapply) {
539: PetscCall(PetscLogEventBegin(PC_MatApply, pc, X, Y, 0));
540: PetscUseTypeMethod(pc, matapply, X, Y);
541: PetscCall(PetscLogEventEnd(PC_MatApply, pc, X, Y, 0));
542: } else {
543: PetscCall(PetscInfo(pc, "PC type %s applying column by column\n", ((PetscObject)pc)->type_name));
544: for (n1 = 0; n1 < N1; ++n1) {
545: PetscCall(MatDenseGetColumnVecRead(X, n1, &cx));
546: PetscCall(MatDenseGetColumnVecWrite(Y, n1, &cy));
547: PetscCall(PCApply(pc, cx, cy));
548: PetscCall(MatDenseRestoreColumnVecWrite(Y, n1, &cy));
549: PetscCall(MatDenseRestoreColumnVecRead(X, n1, &cx));
550: }
551: }
552: PetscFunctionReturn(PETSC_SUCCESS);
553: }
555: /*@
556: PCApplySymmetricLeft - Applies the left part of a symmetric preconditioner to a vector.
558: Collective
560: Input Parameters:
561: + pc - the preconditioner context
562: - x - input vector
564: Output Parameter:
565: . y - output vector
567: Level: developer
569: Note:
570: Currently, this routine is implemented only for `PCICC` and `PCJACOBI` preconditioners.
572: .seealso: `PC`, `PCApply()`, `PCApplySymmetricRight()`
573: @*/
574: PetscErrorCode PCApplySymmetricLeft(PC pc, Vec x, Vec y)
575: {
576: PetscFunctionBegin;
580: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
581: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
582: PetscCall(PCSetUp(pc));
583: PetscCall(VecLockReadPush(x));
584: PetscCall(PetscLogEventBegin(PC_ApplySymmetricLeft, pc, x, y, 0));
585: PetscUseTypeMethod(pc, applysymmetricleft, x, y);
586: PetscCall(PetscLogEventEnd(PC_ApplySymmetricLeft, pc, x, y, 0));
587: PetscCall(VecLockReadPop(x));
588: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
589: PetscFunctionReturn(PETSC_SUCCESS);
590: }
592: /*@
593: PCApplySymmetricRight - Applies the right part of a symmetric preconditioner to a vector.
595: Collective
597: Input Parameters:
598: + pc - the preconditioner context
599: - x - input vector
601: Output Parameter:
602: . y - output vector
604: Level: developer
606: Note:
607: Currently, this routine is implemented only for `PCICC` and `PCJACOBI` preconditioners.
609: .seealso: `PC`, `PCApply()`, `PCApplySymmetricLeft()`
610: @*/
611: PetscErrorCode PCApplySymmetricRight(PC pc, Vec x, Vec y)
612: {
613: PetscFunctionBegin;
617: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
618: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
619: PetscCall(PCSetUp(pc));
620: PetscCall(VecLockReadPush(x));
621: PetscCall(PetscLogEventBegin(PC_ApplySymmetricRight, pc, x, y, 0));
622: PetscUseTypeMethod(pc, applysymmetricright, x, y);
623: PetscCall(PetscLogEventEnd(PC_ApplySymmetricRight, pc, x, y, 0));
624: PetscCall(VecLockReadPop(x));
625: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
626: PetscFunctionReturn(PETSC_SUCCESS);
627: }
629: /*@
630: PCApplyTranspose - Applies the transpose of preconditioner to a vector.
632: Collective
634: Input Parameters:
635: + pc - the preconditioner context
636: - x - input vector
638: Output Parameter:
639: . y - output vector
641: Level: developer
643: Note:
644: For complex numbers this applies the non-Hermitian transpose.
646: Developer Notes:
647: We need to implement a `PCApplyHermitianTranspose()`
649: .seealso: `PC`, `PCApply()`, `PCApplyBAorAB()`, `PCApplyBAorABTranspose()`, `PCApplyTransposeExists()`
650: @*/
651: PetscErrorCode PCApplyTranspose(PC pc, Vec x, Vec y)
652: {
653: PetscFunctionBegin;
657: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
658: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
659: PetscCall(PCSetUp(pc));
660: PetscCall(VecLockReadPush(x));
661: PetscCall(PetscLogEventBegin(PC_Apply, pc, x, y, 0));
662: PetscUseTypeMethod(pc, applytranspose, x, y);
663: PetscCall(PetscLogEventEnd(PC_Apply, pc, x, y, 0));
664: PetscCall(VecLockReadPop(x));
665: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
666: PetscFunctionReturn(PETSC_SUCCESS);
667: }
669: /*@
670: PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation
672: Collective
674: Input Parameter:
675: . pc - the preconditioner context
677: Output Parameter:
678: . flg - `PETSC_TRUE` if a transpose operation is defined
680: Level: developer
682: .seealso: `PC`, `PCApplyTranspose()`
683: @*/
684: PetscErrorCode PCApplyTransposeExists(PC pc, PetscBool *flg)
685: {
686: PetscFunctionBegin;
688: PetscAssertPointer(flg, 2);
689: if (pc->ops->applytranspose) *flg = PETSC_TRUE;
690: else *flg = PETSC_FALSE;
691: PetscFunctionReturn(PETSC_SUCCESS);
692: }
694: /*@
695: PCApplyBAorAB - Applies the preconditioner and operator to a vector. y = B*A*x or y = A*B*x.
697: Collective
699: Input Parameters:
700: + pc - the preconditioner context
701: . side - indicates the preconditioner side, one of `PC_LEFT`, `PC_RIGHT`, or `PC_SYMMETRIC`
702: . x - input vector
703: - work - work vector
705: Output Parameter:
706: . y - output vector
708: Level: developer
710: Note:
711: 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
712: specific `KSPSolve()` method must also be written to handle the post-solve "correction" for the diagonal scaling.
714: .seealso: `PC`, `PCApply()`, `PCApplyTranspose()`, `PCApplyBAorABTranspose()`
715: @*/
716: PetscErrorCode PCApplyBAorAB(PC pc, PCSide side, Vec x, Vec y, Vec work)
717: {
718: PetscFunctionBegin;
724: PetscCheckSameComm(pc, 1, x, 3);
725: PetscCheckSameComm(pc, 1, y, 4);
726: PetscCheckSameComm(pc, 1, work, 5);
727: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
728: PetscCheck(side == PC_LEFT || side == PC_SYMMETRIC || side == PC_RIGHT, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Side must be right, left, or symmetric");
729: PetscCheck(!pc->diagonalscale || side != PC_SYMMETRIC, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot include diagonal scaling with symmetric preconditioner application");
730: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 3, PETSC_TRUE));
732: PetscCall(PCSetUp(pc));
733: if (pc->diagonalscale) {
734: if (pc->ops->applyBA) {
735: Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
736: PetscCall(VecDuplicate(x, &work2));
737: PetscCall(PCDiagonalScaleRight(pc, x, work2));
738: PetscUseTypeMethod(pc, applyBA, side, work2, y, work);
739: PetscCall(PCDiagonalScaleLeft(pc, y, y));
740: PetscCall(VecDestroy(&work2));
741: } else if (side == PC_RIGHT) {
742: PetscCall(PCDiagonalScaleRight(pc, x, y));
743: PetscCall(PCApply(pc, y, work));
744: PetscCall(MatMult(pc->mat, work, y));
745: PetscCall(PCDiagonalScaleLeft(pc, y, y));
746: } else if (side == PC_LEFT) {
747: PetscCall(PCDiagonalScaleRight(pc, x, y));
748: PetscCall(MatMult(pc->mat, y, work));
749: PetscCall(PCApply(pc, work, y));
750: PetscCall(PCDiagonalScaleLeft(pc, y, y));
751: } else PetscCheck(side != PC_SYMMETRIC, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot provide diagonal scaling with symmetric application of preconditioner");
752: } else {
753: if (pc->ops->applyBA) {
754: PetscUseTypeMethod(pc, applyBA, side, x, y, work);
755: } else if (side == PC_RIGHT) {
756: PetscCall(PCApply(pc, x, work));
757: PetscCall(MatMult(pc->mat, work, y));
758: } else if (side == PC_LEFT) {
759: PetscCall(MatMult(pc->mat, x, work));
760: PetscCall(PCApply(pc, work, y));
761: } else if (side == PC_SYMMETRIC) {
762: /* There's an extra copy here; maybe should provide 2 work vectors instead? */
763: PetscCall(PCApplySymmetricRight(pc, x, work));
764: PetscCall(MatMult(pc->mat, work, y));
765: PetscCall(VecCopy(y, work));
766: PetscCall(PCApplySymmetricLeft(pc, work, y));
767: }
768: }
769: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 4, PETSC_FALSE));
770: PetscFunctionReturn(PETSC_SUCCESS);
771: }
773: /*@
774: PCApplyBAorABTranspose - Applies the transpose of the preconditioner
775: and operator to a vector. That is, applies tr(B) * tr(A) with left preconditioning,
776: NOT tr(B*A) = tr(A)*tr(B).
778: Collective
780: Input Parameters:
781: + pc - the preconditioner context
782: . side - indicates the preconditioner side, one of `PC_LEFT`, `PC_RIGHT`, or `PC_SYMMETRIC`
783: . x - input vector
784: - work - work vector
786: Output Parameter:
787: . y - output vector
789: Level: developer
791: Note:
792: 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
793: defined by B'. This is why this has the funny form that it computes tr(B) * tr(A)
795: .seealso: `PC`, `PCApply()`, `PCApplyTranspose()`, `PCApplyBAorAB()`
796: @*/
797: PetscErrorCode PCApplyBAorABTranspose(PC pc, PCSide side, Vec x, Vec y, Vec work)
798: {
799: PetscFunctionBegin;
804: PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
805: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 3, PETSC_TRUE));
806: if (pc->ops->applyBAtranspose) {
807: PetscUseTypeMethod(pc, applyBAtranspose, side, x, y, work);
808: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 4, PETSC_FALSE));
809: PetscFunctionReturn(PETSC_SUCCESS);
810: }
811: PetscCheck(side == PC_LEFT || side == PC_RIGHT, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Side must be right or left");
813: PetscCall(PCSetUp(pc));
814: if (side == PC_RIGHT) {
815: PetscCall(PCApplyTranspose(pc, x, work));
816: PetscCall(MatMultTranspose(pc->mat, work, y));
817: } else if (side == PC_LEFT) {
818: PetscCall(MatMultTranspose(pc->mat, x, work));
819: PetscCall(PCApplyTranspose(pc, work, y));
820: }
821: /* add support for PC_SYMMETRIC */
822: if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 4, PETSC_FALSE));
823: PetscFunctionReturn(PETSC_SUCCESS);
824: }
826: /*@
827: PCApplyRichardsonExists - Determines whether a particular preconditioner has a
828: built-in fast application of Richardson's method.
830: Not Collective
832: Input Parameter:
833: . pc - the preconditioner
835: Output Parameter:
836: . exists - `PETSC_TRUE` or `PETSC_FALSE`
838: Level: developer
840: .seealso: `PC`, `PCRICHARDSON`, `PCApplyRichardson()`
841: @*/
842: PetscErrorCode PCApplyRichardsonExists(PC pc, PetscBool *exists)
843: {
844: PetscFunctionBegin;
846: PetscAssertPointer(exists, 2);
847: if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
848: else *exists = PETSC_FALSE;
849: PetscFunctionReturn(PETSC_SUCCESS);
850: }
852: /*@
853: PCApplyRichardson - Applies several steps of Richardson iteration with
854: the particular preconditioner. This routine is usually used by the
855: Krylov solvers and not the application code directly.
857: Collective
859: Input Parameters:
860: + pc - the preconditioner context
861: . b - the right hand side
862: . w - one work vector
863: . rtol - relative decrease in residual norm convergence criteria
864: . abstol - absolute residual norm convergence criteria
865: . dtol - divergence residual norm increase criteria
866: . its - the number of iterations to apply.
867: - guesszero - if the input x contains nonzero initial guess
869: Output Parameters:
870: + outits - number of iterations actually used (for SOR this always equals its)
871: . reason - the reason the apply terminated
872: - y - the solution (also contains initial guess if guesszero is `PETSC_FALSE`
874: Level: developer
876: Notes:
877: Most preconditioners do not support this function. Use the command
878: `PCApplyRichardsonExists()` to determine if one does.
880: Except for the `PCMG` this routine ignores the convergence tolerances
881: and always runs for the number of iterations
883: .seealso: `PC`, `PCApplyRichardsonExists()`
884: @*/
885: PetscErrorCode PCApplyRichardson(PC pc, Vec b, Vec y, Vec w, PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt its, PetscBool guesszero, PetscInt *outits, PCRichardsonConvergedReason *reason)
886: {
887: PetscFunctionBegin;
892: PetscCheck(b != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "b and y must be different vectors");
893: PetscCall(PCSetUp(pc));
894: PetscUseTypeMethod(pc, applyrichardson, b, y, w, rtol, abstol, dtol, its, guesszero, outits, reason);
895: PetscFunctionReturn(PETSC_SUCCESS);
896: }
898: /*@
899: PCSetFailedReason - Sets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail
901: Logically Collective
903: Input Parameters:
904: + pc - the preconditioner context
905: - reason - the reason it failedx
907: Level: advanced
909: .seealso: `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCFailedReason`
910: @*/
911: PetscErrorCode PCSetFailedReason(PC pc, PCFailedReason reason)
912: {
913: PetscFunctionBegin;
915: pc->failedreason = reason;
916: PetscFunctionReturn(PETSC_SUCCESS);
917: }
919: /*@
920: PCGetFailedReason - Gets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail
922: Logically Collective
924: Input Parameter:
925: . pc - the preconditioner context
927: Output Parameter:
928: . reason - the reason it failed
930: Level: advanced
932: Note:
933: This is the maximum over reason over all ranks in the PC communicator. It is only valid after
934: a call `KSPCheckDot()` or `KSPCheckNorm()` inside a `KSPSolve()` or `PCReduceFailedReason()`.
935: It is not valid immediately after a `PCSetUp()` or `PCApply()`, then use `PCGetFailedReasonRank()`
937: .seealso: `PC`, ``PCCreate()`, `PCApply()`, `PCDestroy()`, `PCGetFailedReasonRank()`, `PCSetFailedReason()`
938: @*/
939: PetscErrorCode PCGetFailedReason(PC pc, PCFailedReason *reason)
940: {
941: PetscFunctionBegin;
943: if (pc->setupcalled < 0) *reason = (PCFailedReason)pc->setupcalled;
944: else *reason = pc->failedreason;
945: PetscFunctionReturn(PETSC_SUCCESS);
946: }
948: /*@
949: PCGetFailedReasonRank - Gets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail on this MPI rank
951: Not Collective
953: Input Parameter:
954: . pc - the preconditioner context
956: Output Parameter:
957: . reason - the reason it failed
959: Level: advanced
961: Note:
962: Different ranks may have different reasons or no reason, see `PCGetFailedReason()`
964: .seealso: `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCGetFailedReason()`, `PCSetFailedReason()`, `PCReduceFailedReason()`
965: @*/
966: PetscErrorCode PCGetFailedReasonRank(PC pc, PCFailedReason *reason)
967: {
968: PetscFunctionBegin;
970: if (pc->setupcalled < 0) *reason = (PCFailedReason)pc->setupcalled;
971: else *reason = pc->failedreason;
972: PetscFunctionReturn(PETSC_SUCCESS);
973: }
975: /*@
976: PCReduceFailedReason - Reduce the failed reason among the MPI processes that share the `PC`
978: Collective
980: Input Parameter:
981: . pc - the preconditioner context
983: Level: advanced
985: Note:
986: Different MPI ranks may have different reasons or no reason, see `PCGetFailedReason()`. This routine
987: makes them have a common value (failure if any MPI process had a failure).
989: .seealso: `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCGetFailedReason()`, `PCSetFailedReason()`
990: @*/
991: PetscErrorCode PCReduceFailedReason(PC pc)
992: {
993: PetscInt buf;
995: PetscFunctionBegin;
997: buf = (PetscInt)pc->failedreason;
998: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &buf, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)pc)));
999: pc->failedreason = (PCFailedReason)buf;
1000: PetscFunctionReturn(PETSC_SUCCESS);
1001: }
1003: /* Next line needed to deactivate KSP_Solve logging */
1004: #include <petsc/private/kspimpl.h>
1006: /*
1007: a setupcall of 0 indicates never setup,
1008: 1 indicates has been previously setup
1009: -1 indicates a PCSetUp() was attempted and failed
1010: */
1011: /*@
1012: PCSetUp - Prepares for the use of a preconditioner.
1014: Collective
1016: Input Parameter:
1017: . pc - the preconditioner context
1019: Level: developer
1021: .seealso: `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`
1022: @*/
1023: PetscErrorCode PCSetUp(PC pc)
1024: {
1025: const char *def;
1026: PetscObjectState matstate, matnonzerostate;
1028: PetscFunctionBegin;
1030: PetscCheck(pc->mat, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Matrix must be set first");
1032: if (pc->setupcalled && pc->reusepreconditioner) {
1033: PetscCall(PetscInfo(pc, "Leaving PC with identical preconditioner since reuse preconditioner is set\n"));
1034: PetscFunctionReturn(PETSC_SUCCESS);
1035: }
1037: PetscCall(PetscObjectStateGet((PetscObject)pc->pmat, &matstate));
1038: PetscCall(MatGetNonzeroState(pc->pmat, &matnonzerostate));
1039: if (!pc->setupcalled) {
1040: PetscCall(PetscInfo(pc, "Setting up PC for first time\n"));
1041: pc->flag = DIFFERENT_NONZERO_PATTERN;
1042: } else if (matstate == pc->matstate) PetscFunctionReturn(PETSC_SUCCESS);
1043: else {
1044: if (matnonzerostate != pc->matnonzerostate) {
1045: PetscCall(PetscInfo(pc, "Setting up PC with different nonzero pattern\n"));
1046: pc->flag = DIFFERENT_NONZERO_PATTERN;
1047: } else {
1048: PetscCall(PetscInfo(pc, "Setting up PC with same nonzero pattern\n"));
1049: pc->flag = SAME_NONZERO_PATTERN;
1050: }
1051: }
1052: pc->matstate = matstate;
1053: pc->matnonzerostate = matnonzerostate;
1055: if (!((PetscObject)pc)->type_name) {
1056: PetscCall(PCGetDefaultType_Private(pc, &def));
1057: PetscCall(PCSetType(pc, def));
1058: }
1060: PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
1061: PetscCall(MatSetErrorIfFailure(pc->mat, pc->erroriffailure));
1062: PetscCall(PetscLogEventBegin(PC_SetUp, pc, 0, 0, 0));
1063: if (pc->ops->setup) {
1064: /* do not log solves and applications of preconditioners while constructing preconditioners; perhaps they should be logged separately from the regular solves */
1065: PetscCall(KSPInitializePackage());
1066: PetscCall(PetscLogEventDeactivatePush(KSP_Solve));
1067: PetscCall(PetscLogEventDeactivatePush(PC_Apply));
1068: PetscUseTypeMethod(pc, setup);
1069: PetscCall(PetscLogEventDeactivatePop(KSP_Solve));
1070: PetscCall(PetscLogEventDeactivatePop(PC_Apply));
1071: }
1072: PetscCall(PetscLogEventEnd(PC_SetUp, pc, 0, 0, 0));
1073: if (!pc->setupcalled) pc->setupcalled = 1;
1074: PetscFunctionReturn(PETSC_SUCCESS);
1075: }
1077: /*@
1078: PCSetUpOnBlocks - Sets up the preconditioner for each block in
1079: the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
1080: methods.
1082: Collective
1084: Input Parameter:
1085: . pc - the preconditioner context
1087: Level: developer
1089: Note:
1090: For nested preconditioners such as `PCBJACOBI` `PCSetUp()` is not called on each sub-`KSP` when `PCSetUp()` is
1091: called on the outer `PC`, this routine ensures it is called.
1093: .seealso: `PC`, `PCSetUp()`, `PCCreate()`, `PCApply()`, `PCDestroy()`
1094: @*/
1095: PetscErrorCode PCSetUpOnBlocks(PC pc)
1096: {
1097: PetscFunctionBegin;
1099: if (!pc->ops->setuponblocks) PetscFunctionReturn(PETSC_SUCCESS);
1100: PetscCall(PetscLogEventBegin(PC_SetUpOnBlocks, pc, 0, 0, 0));
1101: PetscUseTypeMethod(pc, setuponblocks);
1102: PetscCall(PetscLogEventEnd(PC_SetUpOnBlocks, pc, 0, 0, 0));
1103: PetscFunctionReturn(PETSC_SUCCESS);
1104: }
1106: /*@C
1107: PCSetModifySubMatrices - Sets a user-defined routine for modifying the
1108: submatrices that arise within certain subdomain-based preconditioners such as `PCASM`
1110: Logically Collective
1112: Input Parameters:
1113: + pc - the preconditioner context
1114: . func - routine for modifying the submatrices
1115: - ctx - optional user-defined context (may be `NULL`)
1117: Calling sequence of `func`:
1118: + pc - the preconditioner context
1119: . nsub - number of index sets
1120: . row - an array of index sets that contain the global row numbers
1121: that comprise each local submatrix
1122: . col - an array of index sets that contain the global column numbers
1123: that comprise each local submatrix
1124: . submat - array of local submatrices
1125: - ctx - optional user-defined context for private data for the
1126: user-defined func routine (may be `NULL`)
1128: Level: advanced
1130: Notes:
1131: The basic submatrices are extracted from the preconditioner matrix as
1132: usual; the user can then alter these (for example, to set different boundary
1133: conditions for each submatrix) before they are used for the local solves.
1135: `PCSetModifySubMatrices()` MUST be called before `KSPSetUp()` and
1136: `KSPSolve()`.
1138: A routine set by `PCSetModifySubMatrices()` is currently called within
1139: the block Jacobi (`PCBJACOBI`) and additive Schwarz (`PCASM`)
1140: preconditioners. All other preconditioners ignore this routine.
1142: .seealso: `PC`, `PCBJACOBI`, `PCASM`, `PCModifySubMatrices()`
1143: @*/
1144: PetscErrorCode PCSetModifySubMatrices(PC pc, PetscErrorCode (*func)(PC pc, PetscInt nsub, const IS row[], const IS col[], Mat submat[], void *ctx), void *ctx)
1145: {
1146: PetscFunctionBegin;
1148: pc->modifysubmatrices = func;
1149: pc->modifysubmatricesP = ctx;
1150: PetscFunctionReturn(PETSC_SUCCESS);
1151: }
1153: /*@C
1154: PCModifySubMatrices - Calls an optional user-defined routine within
1155: certain preconditioners if one has been set with `PCSetModifySubMatrices()`.
1157: Collective
1159: Input Parameters:
1160: + pc - the preconditioner context
1161: . nsub - the number of local submatrices
1162: . row - an array of index sets that contain the global row numbers
1163: that comprise each local submatrix
1164: . col - an array of index sets that contain the global column numbers
1165: that comprise each local submatrix
1166: . submat - array of local submatrices
1167: - ctx - optional user-defined context for private data for the
1168: user-defined routine (may be null)
1170: Output Parameter:
1171: . submat - array of local submatrices (the entries of which may
1172: have been modified)
1174: Level: developer
1176: Note:
1177: The user should NOT generally call this routine, as it will
1178: automatically be called within certain preconditioners.
1180: .seealso: `PC`, `PCSetModifySubMatrices()`
1181: @*/
1182: PetscErrorCode PCModifySubMatrices(PC pc, PetscInt nsub, const IS row[], const IS col[], Mat submat[], void *ctx)
1183: {
1184: PetscFunctionBegin;
1186: if (!pc->modifysubmatrices) PetscFunctionReturn(PETSC_SUCCESS);
1187: PetscCall(PetscLogEventBegin(PC_ModifySubMatrices, pc, 0, 0, 0));
1188: PetscCall((*pc->modifysubmatrices)(pc, nsub, row, col, submat, ctx));
1189: PetscCall(PetscLogEventEnd(PC_ModifySubMatrices, pc, 0, 0, 0));
1190: PetscFunctionReturn(PETSC_SUCCESS);
1191: }
1193: /*@
1194: PCSetOperators - Sets the matrix associated with the linear system and
1195: a (possibly) different one associated with the preconditioner.
1197: Logically Collective
1199: Input Parameters:
1200: + pc - the preconditioner context
1201: . Amat - the matrix that defines the linear system
1202: - Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
1204: Level: intermediate
1206: Notes:
1207: Passing a `NULL` for `Amat` or `Pmat` removes the matrix that is currently used.
1209: If you wish to replace either `Amat` or `Pmat` but leave the other one untouched then
1210: first call `KSPGetOperators()` to get the one you wish to keep, call `PetscObjectReference()`
1211: on it and then pass it back in in your call to `KSPSetOperators()`.
1213: More Notes about Repeated Solution of Linear Systems:
1214: PETSc does NOT reset the matrix entries of either `Amat` or `Pmat`
1215: to zero after a linear solve; the user is completely responsible for
1216: matrix assembly. See the routine `MatZeroEntries()` if desiring to
1217: zero all elements of a matrix.
1219: .seealso: `PC`, `PCGetOperators()`, `MatZeroEntries()`
1220: @*/
1221: PetscErrorCode PCSetOperators(PC pc, Mat Amat, Mat Pmat)
1222: {
1223: PetscInt m1, n1, m2, n2;
1225: PetscFunctionBegin;
1229: if (Amat) PetscCheckSameComm(pc, 1, Amat, 2);
1230: if (Pmat) PetscCheckSameComm(pc, 1, Pmat, 3);
1231: if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
1232: PetscCall(MatGetLocalSize(Amat, &m1, &n1));
1233: PetscCall(MatGetLocalSize(pc->mat, &m2, &n2));
1234: 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);
1235: PetscCall(MatGetLocalSize(Pmat, &m1, &n1));
1236: PetscCall(MatGetLocalSize(pc->pmat, &m2, &n2));
1237: 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);
1238: }
1240: if (Pmat != pc->pmat) {
1241: /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
1242: pc->matnonzerostate = -1;
1243: pc->matstate = -1;
1244: }
1246: /* reference first in case the matrices are the same */
1247: if (Amat) PetscCall(PetscObjectReference((PetscObject)Amat));
1248: PetscCall(MatDestroy(&pc->mat));
1249: if (Pmat) PetscCall(PetscObjectReference((PetscObject)Pmat));
1250: PetscCall(MatDestroy(&pc->pmat));
1251: pc->mat = Amat;
1252: pc->pmat = Pmat;
1253: PetscFunctionReturn(PETSC_SUCCESS);
1254: }
1256: /*@
1257: PCSetReusePreconditioner - reuse the current preconditioner even if the operator in the preconditioner has changed.
1259: Logically Collective
1261: Input Parameters:
1262: + pc - the preconditioner context
1263: - flag - `PETSC_TRUE` do not compute a new preconditioner, `PETSC_FALSE` do compute a new preconditioner
1265: Level: intermediate
1267: Note:
1268: Normally if a matrix inside a `PC` changes the `PC` automatically updates itself using information from the changed matrix. This option
1269: prevents this.
1271: .seealso: `PC`, `PCGetOperators()`, `MatZeroEntries()`, `PCGetReusePreconditioner()`, `KSPSetReusePreconditioner()`
1272: @*/
1273: PetscErrorCode PCSetReusePreconditioner(PC pc, PetscBool flag)
1274: {
1275: PetscFunctionBegin;
1278: pc->reusepreconditioner = flag;
1279: PetscFunctionReturn(PETSC_SUCCESS);
1280: }
1282: /*@
1283: PCGetReusePreconditioner - Determines if the `PC` reuses the current preconditioner even if the operator in the preconditioner has changed.
1285: Not Collective
1287: Input Parameter:
1288: . pc - the preconditioner context
1290: Output Parameter:
1291: . flag - `PETSC_TRUE` do not compute a new preconditioner, `PETSC_FALSE` do compute a new preconditioner
1293: Level: intermediate
1295: .seealso: `PC`, `PCGetOperators()`, `MatZeroEntries()`, `PCSetReusePreconditioner()`
1296: @*/
1297: PetscErrorCode PCGetReusePreconditioner(PC pc, PetscBool *flag)
1298: {
1299: PetscFunctionBegin;
1301: PetscAssertPointer(flag, 2);
1302: *flag = pc->reusepreconditioner;
1303: PetscFunctionReturn(PETSC_SUCCESS);
1304: }
1306: /*@
1307: PCGetOperators - Gets the matrix associated with the linear system and
1308: possibly a different one associated with the preconditioner.
1310: Not Collective, though parallel mats are returned if pc is parallel
1312: Input Parameter:
1313: . pc - the preconditioner context
1315: Output Parameters:
1316: + Amat - the matrix defining the linear system
1317: - Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.
1319: Level: intermediate
1321: Note:
1322: Does not increase the reference count of the matrices, so you should not destroy them
1324: Alternative usage: If the operators have NOT been set with `KSPSetOperators()`/`PCSetOperators()` then the operators
1325: are created in `PC` and returned to the user. In this case, if both operators
1326: mat and pmat are requested, two DIFFERENT operators will be returned. If
1327: only one is requested both operators in the PC will be the same (i.e. as
1328: if one had called `KSPSetOperators()`/`PCSetOperators()` with the same argument for both Mats).
1329: The user must set the sizes of the returned matrices and their type etc just
1330: as if the user created them with `MatCreate()`. For example,
1332: .vb
1333: KSP/PCGetOperators(ksp/pc,&Amat,NULL); is equivalent to
1334: set size, type, etc of Amat
1336: MatCreate(comm,&mat);
1337: KSP/PCSetOperators(ksp/pc,Amat,Amat);
1338: PetscObjectDereference((PetscObject)mat);
1339: set size, type, etc of Amat
1340: .ve
1342: and
1344: .vb
1345: KSP/PCGetOperators(ksp/pc,&Amat,&Pmat); is equivalent to
1346: set size, type, etc of Amat and Pmat
1348: MatCreate(comm,&Amat);
1349: MatCreate(comm,&Pmat);
1350: KSP/PCSetOperators(ksp/pc,Amat,Pmat);
1351: PetscObjectDereference((PetscObject)Amat);
1352: PetscObjectDereference((PetscObject)Pmat);
1353: set size, type, etc of Amat and Pmat
1354: .ve
1356: The rationale for this support is so that when creating a `TS`, `SNES`, or `KSP` the hierarchy
1357: of underlying objects (i.e. `SNES`, `KSP`, `PC`, `Mat`) and their lifespans can be completely
1358: managed by the top most level object (i.e. the `TS`, `SNES`, or `KSP`). Another way to look
1359: at this is when you create a `SNES` you do not NEED to create a `KSP` and attach it to
1360: the `SNES` object (the `SNES` object manages it for you). Similarly when you create a KSP
1361: you do not need to attach a `PC` to it (the `KSP` object manages the `PC` object for you).
1362: Thus, why should YOU have to create the `Mat` and attach it to the `SNES`/`KSP`/`PC`, when
1363: it can be created for you?
1365: .seealso: `PC`, `PCSetOperators()`, `KSPGetOperators()`, `KSPSetOperators()`, `PCGetOperatorsSet()`
1366: @*/
1367: PetscErrorCode PCGetOperators(PC pc, Mat *Amat, Mat *Pmat)
1368: {
1369: PetscFunctionBegin;
1371: if (Amat) {
1372: if (!pc->mat) {
1373: if (pc->pmat && !Pmat) { /* Pmat has been set, but user did not request it, so use for Amat */
1374: pc->mat = pc->pmat;
1375: PetscCall(PetscObjectReference((PetscObject)pc->mat));
1376: } else { /* both Amat and Pmat are empty */
1377: PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &pc->mat));
1378: if (!Pmat) { /* user did NOT request Pmat, so make same as Amat */
1379: pc->pmat = pc->mat;
1380: PetscCall(PetscObjectReference((PetscObject)pc->pmat));
1381: }
1382: }
1383: }
1384: *Amat = pc->mat;
1385: }
1386: if (Pmat) {
1387: if (!pc->pmat) {
1388: if (pc->mat && !Amat) { /* Amat has been set but was not requested, so use for pmat */
1389: pc->pmat = pc->mat;
1390: PetscCall(PetscObjectReference((PetscObject)pc->pmat));
1391: } else {
1392: PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &pc->pmat));
1393: if (!Amat) { /* user did NOT request Amat, so make same as Pmat */
1394: pc->mat = pc->pmat;
1395: PetscCall(PetscObjectReference((PetscObject)pc->mat));
1396: }
1397: }
1398: }
1399: *Pmat = pc->pmat;
1400: }
1401: PetscFunctionReturn(PETSC_SUCCESS);
1402: }
1404: /*@C
1405: PCGetOperatorsSet - Determines if the matrix associated with the linear system and
1406: possibly a different one associated with the preconditioner have been set in the `PC`.
1408: Not Collective, though the results on all processes should be the same
1410: Input Parameter:
1411: . pc - the preconditioner context
1413: Output Parameters:
1414: + mat - the matrix associated with the linear system was set
1415: - pmat - matrix associated with the preconditioner was set, usually the same
1417: Level: intermediate
1419: .seealso: `PC`, `PCSetOperators()`, `KSPGetOperators()`, `KSPSetOperators()`, `PCGetOperators()`
1420: @*/
1421: PetscErrorCode PCGetOperatorsSet(PC pc, PetscBool *mat, PetscBool *pmat)
1422: {
1423: PetscFunctionBegin;
1425: if (mat) *mat = (pc->mat) ? PETSC_TRUE : PETSC_FALSE;
1426: if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1427: PetscFunctionReturn(PETSC_SUCCESS);
1428: }
1430: /*@
1431: PCFactorGetMatrix - Gets the factored matrix from the
1432: preconditioner context. This routine is valid only for the `PCLU`,
1433: `PCILU`, `PCCHOLESKY`, and `PCICC` methods.
1435: Not Collective though mat is parallel if pc is parallel
1437: Input Parameter:
1438: . pc - the preconditioner context
1440: Output Parameters:
1441: . mat - the factored matrix
1443: Level: advanced
1445: Note:
1446: Does not increase the reference count for the matrix so DO NOT destroy it
1448: .seealso: `PC`, `PCLU`, `PCILU`, `PCCHOLESKY`, `PCICC`
1449: @*/
1450: PetscErrorCode PCFactorGetMatrix(PC pc, Mat *mat)
1451: {
1452: PetscFunctionBegin;
1454: PetscAssertPointer(mat, 2);
1455: PetscCall(PCFactorSetUpMatSolverType(pc));
1456: PetscUseTypeMethod(pc, getfactoredmatrix, mat);
1457: PetscFunctionReturn(PETSC_SUCCESS);
1458: }
1460: /*@C
1461: PCSetOptionsPrefix - Sets the prefix used for searching for all
1462: `PC` options in the database.
1464: Logically Collective
1466: Input Parameters:
1467: + pc - the preconditioner context
1468: - prefix - the prefix string to prepend to all `PC` option requests
1470: Note:
1471: A hyphen (-) must NOT be given at the beginning of the prefix name.
1472: The first character of all runtime options is AUTOMATICALLY the
1473: hyphen.
1475: Level: advanced
1477: .seealso: `PC`, `PCSetFromOptions`, `PCAppendOptionsPrefix()`, `PCGetOptionsPrefix()`
1478: @*/
1479: PetscErrorCode PCSetOptionsPrefix(PC pc, const char prefix[])
1480: {
1481: PetscFunctionBegin;
1483: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc, prefix));
1484: PetscFunctionReturn(PETSC_SUCCESS);
1485: }
1487: /*@C
1488: PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1489: `PC` options in the database.
1491: Logically Collective
1493: Input Parameters:
1494: + pc - the preconditioner context
1495: - prefix - the prefix string to prepend to all `PC` option requests
1497: Note:
1498: A hyphen (-) must NOT be given at the beginning of the prefix name.
1499: The first character of all runtime options is AUTOMATICALLY the
1500: hyphen.
1502: Level: advanced
1504: .seealso: `PC`, `PCSetFromOptions`, `PCSetOptionsPrefix()`, `PCGetOptionsPrefix()`
1505: @*/
1506: PetscErrorCode PCAppendOptionsPrefix(PC pc, const char prefix[])
1507: {
1508: PetscFunctionBegin;
1510: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)pc, prefix));
1511: PetscFunctionReturn(PETSC_SUCCESS);
1512: }
1514: /*@C
1515: PCGetOptionsPrefix - Gets the prefix used for searching for all
1516: PC options in the database.
1518: Not Collective
1520: Input Parameter:
1521: . pc - the preconditioner context
1523: Output Parameter:
1524: . prefix - pointer to the prefix string used, is returned
1526: Fortran Notes:
1527: The user should pass in a string 'prefix' of
1528: sufficient length to hold the prefix.
1530: Level: advanced
1532: .seealso: `PC`, `PCSetFromOptions`, `PCSetOptionsPrefix()`, `PCAppendOptionsPrefix()`
1533: @*/
1534: PetscErrorCode PCGetOptionsPrefix(PC pc, const char *prefix[])
1535: {
1536: PetscFunctionBegin;
1538: PetscAssertPointer(prefix, 2);
1539: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, prefix));
1540: PetscFunctionReturn(PETSC_SUCCESS);
1541: }
1543: /*
1544: Indicates the right hand side will be changed by KSPSolve(), this occurs for a few
1545: preconditioners including BDDC and Eisentat that transform the equations before applying
1546: the Krylov methods
1547: */
1548: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC pc, PetscBool *change)
1549: {
1550: PetscFunctionBegin;
1552: PetscAssertPointer(change, 2);
1553: *change = PETSC_FALSE;
1554: PetscTryMethod(pc, "PCPreSolveChangeRHS_C", (PC, PetscBool *), (pc, change));
1555: PetscFunctionReturn(PETSC_SUCCESS);
1556: }
1558: /*@
1559: PCPreSolve - Optional pre-solve phase, intended for any
1560: preconditioner-specific actions that must be performed before
1561: the iterative solve itself.
1563: Collective
1565: Input Parameters:
1566: + pc - the preconditioner context
1567: - ksp - the Krylov subspace context
1569: Level: developer
1571: Example Usage:
1572: .vb
1573: PCPreSolve(pc,ksp);
1574: KSPSolve(ksp,b,x);
1575: PCPostSolve(pc,ksp);
1576: .ve
1578: Notes:
1579: The pre-solve phase is distinct from the `PCSetUp()` phase.
1581: `KSPSolve()` calls this directly, so is rarely called by the user.
1583: .seealso: `PC`, `PCPostSolve()`
1584: @*/
1585: PetscErrorCode PCPreSolve(PC pc, KSP ksp)
1586: {
1587: Vec x, rhs;
1589: PetscFunctionBegin;
1592: pc->presolvedone++;
1593: PetscCheck(pc->presolvedone <= 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot embed PCPreSolve() more than twice");
1594: PetscCall(KSPGetSolution(ksp, &x));
1595: PetscCall(KSPGetRhs(ksp, &rhs));
1597: if (pc->ops->presolve) PetscUseTypeMethod(pc, presolve, ksp, rhs, x);
1598: else if (pc->presolve) PetscCall((pc->presolve)(pc, ksp));
1599: PetscFunctionReturn(PETSC_SUCCESS);
1600: }
1602: /*@C
1603: PCSetPreSolve - Sets function used by `PCPreSolve()` which is intended for any
1604: preconditioner-specific actions that must be performed before
1605: the iterative solve itself.
1607: Logically Collective
1609: Input Parameters:
1610: + pc - the preconditioner object
1611: - presolve - the function to call before the solve
1613: Calling sequence of `presolve`:
1614: + pc - the `PC` context
1615: - ksp - the `KSP` context
1617: Level: developer
1619: .seealso: `PC`, `PCSetUp()`, `PCPreSolve()`
1620: @*/
1621: PetscErrorCode PCSetPreSolve(PC pc, PetscErrorCode (*presolve)(PC pc, KSP ksp))
1622: {
1623: PetscFunctionBegin;
1625: pc->presolve = presolve;
1626: PetscFunctionReturn(PETSC_SUCCESS);
1627: }
1629: /*@
1630: PCPostSolve - Optional post-solve phase, intended for any
1631: preconditioner-specific actions that must be performed after
1632: the iterative solve itself.
1634: Collective
1636: Input Parameters:
1637: + pc - the preconditioner context
1638: - ksp - the Krylov subspace context
1640: Example Usage:
1641: .vb
1642: PCPreSolve(pc,ksp);
1643: KSPSolve(ksp,b,x);
1644: PCPostSolve(pc,ksp);
1645: .ve
1647: Note:
1648: `KSPSolve()` calls this routine directly, so it is rarely called by the user.
1650: Level: developer
1652: .seealso: `PC`, `PCSetPostSolve()`, `PCSetPresolve()`, `PCPreSolve()`, `KSPSolve()`
1653: @*/
1654: PetscErrorCode PCPostSolve(PC pc, KSP ksp)
1655: {
1656: Vec x, rhs;
1658: PetscFunctionBegin;
1661: pc->presolvedone--;
1662: PetscCall(KSPGetSolution(ksp, &x));
1663: PetscCall(KSPGetRhs(ksp, &rhs));
1664: PetscTryTypeMethod(pc, postsolve, ksp, rhs, x);
1665: PetscFunctionReturn(PETSC_SUCCESS);
1666: }
1668: /*@C
1669: PCLoad - Loads a `PC` that has been stored in binary with `PCView()`.
1671: Collective
1673: Input Parameters:
1674: + newdm - the newly loaded `PC`, this needs to have been created with `PCCreate()` or
1675: some related function before a call to `PCLoad()`.
1676: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()`
1678: Level: intermediate
1680: Note:
1681: The type is determined by the data in the file, any type set into the PC before this call is ignored.
1683: .seealso: `PC`, `PetscViewerBinaryOpen()`, `PCView()`, `MatLoad()`, `VecLoad()`
1684: @*/
1685: PetscErrorCode PCLoad(PC newdm, PetscViewer viewer)
1686: {
1687: PetscBool isbinary;
1688: PetscInt classid;
1689: char type[256];
1691: PetscFunctionBegin;
1694: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1695: PetscCheck(isbinary, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1697: PetscCall(PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT));
1698: PetscCheck(classid == PC_FILE_CLASSID, PetscObjectComm((PetscObject)newdm), PETSC_ERR_ARG_WRONG, "Not PC next in file");
1699: PetscCall(PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR));
1700: PetscCall(PCSetType(newdm, type));
1701: PetscTryTypeMethod(newdm, load, viewer);
1702: PetscFunctionReturn(PETSC_SUCCESS);
1703: }
1705: #include <petscdraw.h>
1706: #if defined(PETSC_HAVE_SAWS)
1707: #include <petscviewersaws.h>
1708: #endif
1710: /*@C
1711: PCViewFromOptions - View from the `PC` based on options in the database
1713: Collective
1715: Input Parameters:
1716: + A - the `PC` context
1717: . obj - Optional object that provides the options prefix
1718: - name - command line option
1720: Level: intermediate
1722: .seealso: `PC`, `PCView`, `PetscObjectViewFromOptions()`, `PCCreate()`
1723: @*/
1724: PetscErrorCode PCViewFromOptions(PC A, PetscObject obj, const char name[])
1725: {
1726: PetscFunctionBegin;
1728: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
1729: PetscFunctionReturn(PETSC_SUCCESS);
1730: }
1732: /*@C
1733: PCView - Prints information about the `PC`
1735: Collective
1737: Input Parameters:
1738: + pc - the `PC` context
1739: - viewer - optional visualization context
1741: Level: developer
1743: Notes:
1744: The available visualization contexts include
1745: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
1746: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
1747: output where only the first processor opens
1748: the file. All other processors send their
1749: data to the first processor to print.
1751: The user can open an alternative visualization contexts with
1752: `PetscViewerASCIIOpen()` (output to a specified file).
1754: .seealso: `PC`, `PetscViewer`, `KSPView()`, `PetscViewerASCIIOpen()`
1755: @*/
1756: PetscErrorCode PCView(PC pc, PetscViewer viewer)
1757: {
1758: PCType cstr;
1759: PetscViewerFormat format;
1760: PetscBool iascii, isstring, isbinary, isdraw, pop = PETSC_FALSE;
1761: #if defined(PETSC_HAVE_SAWS)
1762: PetscBool issaws;
1763: #endif
1765: PetscFunctionBegin;
1767: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc), &viewer));
1769: PetscCheckSameComm(pc, 1, viewer, 2);
1771: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1772: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
1773: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1774: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1775: #if defined(PETSC_HAVE_SAWS)
1776: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws));
1777: #endif
1779: if (iascii) {
1780: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)pc, viewer));
1781: if (!pc->setupcalled) PetscCall(PetscViewerASCIIPrintf(viewer, " PC has not been set up so information may be incomplete\n"));
1782: PetscCall(PetscViewerASCIIPushTab(viewer));
1783: PetscTryTypeMethod(pc, view, viewer);
1784: PetscCall(PetscViewerASCIIPopTab(viewer));
1785: if (pc->mat) {
1786: PetscCall(PetscViewerGetFormat(viewer, &format));
1787: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
1788: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
1789: pop = PETSC_TRUE;
1790: }
1791: if (pc->pmat == pc->mat) {
1792: PetscCall(PetscViewerASCIIPrintf(viewer, " linear system matrix = precond matrix:\n"));
1793: PetscCall(PetscViewerASCIIPushTab(viewer));
1794: PetscCall(MatView(pc->mat, viewer));
1795: PetscCall(PetscViewerASCIIPopTab(viewer));
1796: } else {
1797: if (pc->pmat) {
1798: PetscCall(PetscViewerASCIIPrintf(viewer, " linear system matrix followed by preconditioner matrix:\n"));
1799: } else {
1800: PetscCall(PetscViewerASCIIPrintf(viewer, " linear system matrix:\n"));
1801: }
1802: PetscCall(PetscViewerASCIIPushTab(viewer));
1803: PetscCall(MatView(pc->mat, viewer));
1804: if (pc->pmat) PetscCall(MatView(pc->pmat, viewer));
1805: PetscCall(PetscViewerASCIIPopTab(viewer));
1806: }
1807: if (pop) PetscCall(PetscViewerPopFormat(viewer));
1808: }
1809: } else if (isstring) {
1810: PetscCall(PCGetType(pc, &cstr));
1811: PetscCall(PetscViewerStringSPrintf(viewer, " PCType: %-7.7s", cstr));
1812: PetscTryTypeMethod(pc, view, viewer);
1813: if (pc->mat) PetscCall(MatView(pc->mat, viewer));
1814: if (pc->pmat && pc->pmat != pc->mat) PetscCall(MatView(pc->pmat, viewer));
1815: } else if (isbinary) {
1816: PetscInt classid = PC_FILE_CLASSID;
1817: MPI_Comm comm;
1818: PetscMPIInt rank;
1819: char type[256];
1821: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
1822: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1823: if (rank == 0) {
1824: PetscCall(PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT));
1825: PetscCall(PetscStrncpy(type, ((PetscObject)pc)->type_name, 256));
1826: PetscCall(PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR));
1827: }
1828: PetscTryTypeMethod(pc, view, viewer);
1829: } else if (isdraw) {
1830: PetscDraw draw;
1831: char str[25];
1832: PetscReal x, y, bottom, h;
1833: PetscInt n;
1835: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1836: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
1837: if (pc->mat) {
1838: PetscCall(MatGetSize(pc->mat, &n, NULL));
1839: PetscCall(PetscSNPrintf(str, 25, "PC: %s (%" PetscInt_FMT ")", ((PetscObject)pc)->type_name, n));
1840: } else {
1841: PetscCall(PetscSNPrintf(str, 25, "PC: %s", ((PetscObject)pc)->type_name));
1842: }
1843: PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h));
1844: bottom = y - h;
1845: PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
1846: PetscTryTypeMethod(pc, view, viewer);
1847: PetscCall(PetscDrawPopCurrentPoint(draw));
1848: #if defined(PETSC_HAVE_SAWS)
1849: } else if (issaws) {
1850: PetscMPIInt rank;
1852: PetscCall(PetscObjectName((PetscObject)pc));
1853: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1854: if (!((PetscObject)pc)->amsmem && rank == 0) PetscCall(PetscObjectViewSAWs((PetscObject)pc, viewer));
1855: if (pc->mat) PetscCall(MatView(pc->mat, viewer));
1856: if (pc->pmat && pc->pmat != pc->mat) PetscCall(MatView(pc->pmat, viewer));
1857: #endif
1858: }
1859: PetscFunctionReturn(PETSC_SUCCESS);
1860: }
1862: /*@C
1863: PCRegister - Adds a method to the preconditioner package.
1865: Not collective
1867: Input Parameters:
1868: + sname - name of a new user-defined solver
1869: - function - routine to create method context
1871: Example Usage:
1872: .vb
1873: PCRegister("my_solver", MySolverCreate);
1874: .ve
1876: Then, your solver can be chosen with the procedural interface via
1877: $ PCSetType(pc, "my_solver")
1878: or at runtime via the option
1879: $ -pc_type my_solver
1881: Level: advanced
1883: Note:
1884: `PCRegister()` may be called multiple times to add several user-defined preconditioners.
1886: .seealso: `PC`, `PCType`, `PCRegisterAll()`
1887: @*/
1888: PetscErrorCode PCRegister(const char sname[], PetscErrorCode (*function)(PC))
1889: {
1890: PetscFunctionBegin;
1891: PetscCall(PCInitializePackage());
1892: PetscCall(PetscFunctionListAdd(&PCList, sname, function));
1893: PetscFunctionReturn(PETSC_SUCCESS);
1894: }
1896: static PetscErrorCode MatMult_PC(Mat A, Vec X, Vec Y)
1897: {
1898: PC pc;
1900: PetscFunctionBegin;
1901: PetscCall(MatShellGetContext(A, &pc));
1902: PetscCall(PCApply(pc, X, Y));
1903: PetscFunctionReturn(PETSC_SUCCESS);
1904: }
1906: /*@C
1907: PCComputeOperator - Computes the explicit preconditioned operator.
1909: Collective
1911: Input Parameters:
1912: + pc - the preconditioner object
1913: - mattype - the matrix type to be used for the operator
1915: Output Parameter:
1916: . mat - the explicit preconditioned operator
1918: Level: advanced
1920: Note:
1921: This computation is done by applying the operators to columns of the identity matrix.
1922: This routine is costly in general, and is recommended for use only with relatively small systems.
1923: Currently, this routine uses a dense matrix format when mattype == NULL
1925: .seealso: `PC`, `KSPComputeOperator()`, `MatType`
1926: @*/
1927: PetscErrorCode PCComputeOperator(PC pc, MatType mattype, Mat *mat)
1928: {
1929: PetscInt N, M, m, n;
1930: Mat A, Apc;
1932: PetscFunctionBegin;
1934: PetscAssertPointer(mat, 3);
1935: PetscCall(PCGetOperators(pc, &A, NULL));
1936: PetscCall(MatGetLocalSize(A, &m, &n));
1937: PetscCall(MatGetSize(A, &M, &N));
1938: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)pc), m, n, M, N, pc, &Apc));
1939: PetscCall(MatShellSetOperation(Apc, MATOP_MULT, (void (*)(void))MatMult_PC));
1940: PetscCall(MatComputeOperator(Apc, mattype, mat));
1941: PetscCall(MatDestroy(&Apc));
1942: PetscFunctionReturn(PETSC_SUCCESS);
1943: }
1945: /*@
1946: PCSetCoordinates - sets the coordinates of all the nodes on the local process
1948: Collective
1950: Input Parameters:
1951: + pc - the solver context
1952: . dim - the dimension of the coordinates 1, 2, or 3
1953: . nloc - the blocked size of the coordinates array
1954: - coords - the coordinates array
1956: Level: intermediate
1958: Note:
1959: `coords` is an array of the dim coordinates for the nodes on
1960: the local processor, of size `dim`*`nloc`.
1961: If there are 108 equation on a processor
1962: for a displacement finite element discretization of elasticity (so
1963: that there are nloc = 36 = 108/3 nodes) then the array must have 108
1964: double precision values (ie, 3 * 36). These x y z coordinates
1965: should be ordered for nodes 0 to N-1 like so: [ 0.x, 0.y, 0.z, 1.x,
1966: ... , N-1.z ].
1968: .seealso: `PC`, `MatSetNearNullSpace()`
1969: @*/
1970: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[])
1971: {
1972: PetscFunctionBegin;
1975: PetscTryMethod(pc, "PCSetCoordinates_C", (PC, PetscInt, PetscInt, PetscReal[]), (pc, dim, nloc, coords));
1976: PetscFunctionReturn(PETSC_SUCCESS);
1977: }
1979: /*@
1980: PCGetInterpolations - Gets interpolation matrices for all levels (except level 0)
1982: Logically Collective
1984: Input Parameter:
1985: . pc - the precondition context
1987: Output Parameters:
1988: + num_levels - the number of levels
1989: - interpolations - the interpolation matrices (size of num_levels-1)
1991: Level: advanced
1993: Developer Notes:
1994: Why is this here instead of in `PCMG` etc?
1996: .seealso: `PC`, `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetInterpolation()`, `PCGetCoarseOperators()`
1997: @*/
1998: PetscErrorCode PCGetInterpolations(PC pc, PetscInt *num_levels, Mat *interpolations[])
1999: {
2000: PetscFunctionBegin;
2002: PetscAssertPointer(num_levels, 2);
2003: PetscAssertPointer(interpolations, 3);
2004: PetscUseMethod(pc, "PCGetInterpolations_C", (PC, PetscInt *, Mat *[]), (pc, num_levels, interpolations));
2005: PetscFunctionReturn(PETSC_SUCCESS);
2006: }
2008: /*@
2009: PCGetCoarseOperators - Gets coarse operator matrices for all levels (except the finest level)
2011: Logically Collective
2013: Input Parameter:
2014: . pc - the precondition context
2016: Output Parameters:
2017: + num_levels - the number of levels
2018: - coarseOperators - the coarse operator matrices (size of num_levels-1)
2020: Level: advanced
2022: Developer Notes:
2023: Why is this here instead of in `PCMG` etc?
2025: .seealso: `PC`, `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetRScale()`, `PCMGGetInterpolation()`, `PCGetInterpolations()`
2026: @*/
2027: PetscErrorCode PCGetCoarseOperators(PC pc, PetscInt *num_levels, Mat *coarseOperators[])
2028: {
2029: PetscFunctionBegin;
2031: PetscAssertPointer(num_levels, 2);
2032: PetscAssertPointer(coarseOperators, 3);
2033: PetscUseMethod(pc, "PCGetCoarseOperators_C", (PC, PetscInt *, Mat *[]), (pc, num_levels, coarseOperators));
2034: PetscFunctionReturn(PETSC_SUCCESS);
2035: }