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