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: }