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:   MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
 21:   if (pc->pmat) {
 22:     MatHasOperation(pc->pmat, MATOP_GET_DIAGONAL_BLOCK, &hasop);
 23:     if (size == 1) {
 24:       MatGetFactorAvailable(pc->pmat, "petsc", MAT_FACTOR_ICC, &flg1);
 25:       MatGetFactorAvailable(pc->pmat, "petsc", MAT_FACTOR_ILU, &flg2);
 26:       MatIsSymmetricKnown(pc->pmat, &set, &flg3);
 27:       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:   return 0;
 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: {
 74:   PetscTryTypeMethod(pc, reset);
 75:   VecDestroy(&pc->diagonalscaleright);
 76:   VecDestroy(&pc->diagonalscaleleft);
 77:   MatDestroy(&pc->pmat);
 78:   MatDestroy(&pc->mat);

 80:   pc->setupcalled = 0;
 81:   return 0;
 82: }

 84: /*@C
 85:    PCDestroy - Destroys `PC` context that was created with `PCCreate()`.

 87:    Collective

 89:    Input Parameter:
 90: .  pc - the preconditioner context

 92:    Level: developer

 94: .seealso: `PC`, `PCCreate()`, `PCSetUp()`
 95: @*/
 96: PetscErrorCode PCDestroy(PC *pc)
 97: {
 98:   if (!*pc) return 0;
100:   if (--((PetscObject)(*pc))->refct > 0) {
101:     *pc = NULL;
102:     return 0;
103:   }

105:   PCReset(*pc);

107:   /* if memory was published with SAWs then destroy it */
108:   PetscObjectSAWsViewOff((PetscObject)*pc);
109:   PetscTryTypeMethod((*pc), destroy);
110:   DMDestroy(&(*pc)->dm);
111:   PetscHeaderDestroy(pc);
112:   return 0;
113: }

115: /*@C
116:    PCGetDiagonalScale - Indicates if the preconditioner applies an additional left and right
117:       scaling as needed by certain time-stepping codes.

119:    Logically Collective

121:    Input Parameter:
122: .  pc - the preconditioner context

124:    Output Parameter:
125: .  flag - `PETSC_TRUE` if it applies the scaling

127:    Level: developer

129:    Note:
130:     If this returns `PETSC_TRUE` then the system solved via the Krylov method is
131: .vb
132:       D M A D^{-1} y = D M b  for left preconditioning or
133:       D A M D^{-1} z = D b for right preconditioning
134: .ve

136: .seealso: `PC`, `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCDiagonalScaleRight()`, `PCSetDiagonalScale()`
137: @*/
138: PetscErrorCode PCGetDiagonalScale(PC pc, PetscBool *flag)
139: {
142:   *flag = pc->diagonalscale;
143:   return 0;
144: }

146: /*@
147:    PCSetDiagonalScale - Indicates the left scaling to use to apply an additional left and right
148:       scaling as needed by certain time-stepping codes.

150:    Logically Collective

152:    Input Parameters:
153: +  pc - the preconditioner context
154: -  s - scaling vector

156:    Level: intermediate

158:    Notes:
159:     The system solved via the Krylov method is
160: .vb
161:            D M A D^{-1} y = D M b  for left preconditioning or
162:            D A M D^{-1} z = D b for right preconditioning
163: .ve

165:    `PCDiagonalScaleLeft()` scales a vector by D. `PCDiagonalScaleRight()` scales a vector by D^{-1}.

167: .seealso: `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCDiagonalScaleRight()`, `PCGetDiagonalScale()`
168: @*/
169: PetscErrorCode PCSetDiagonalScale(PC pc, Vec s)
170: {
173:   pc->diagonalscale = PETSC_TRUE;

175:   PetscObjectReference((PetscObject)s);
176:   VecDestroy(&pc->diagonalscaleleft);

178:   pc->diagonalscaleleft = s;

180:   VecDuplicate(s, &pc->diagonalscaleright);
181:   VecCopy(s, pc->diagonalscaleright);
182:   VecReciprocal(pc->diagonalscaleright);
183:   return 0;
184: }

186: /*@
187:    PCDiagonalScaleLeft - Scales a vector by the left scaling as needed by certain time-stepping codes.

189:    Logically Collective

191:    Input Parameters:
192: +  pc - the preconditioner context
193: .  in - input vector
194: -  out - scaled vector (maybe the same as in)

196:    Level: intermediate

198:    Notes:
199:     The system solved via the Krylov method is
200: .vb
201:         D M A D^{-1} y = D M b  for left preconditioning or
202:         D A M D^{-1} z = D b for right preconditioning
203: .ve

205:    `PCDiagonalScaleLeft()` scales a vector by D. `PCDiagonalScaleRight()` scales a vector by D^{-1}.

207:    If diagonal scaling is turned off and in is not out then in is copied to out

209: .seealso: `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleSet()`, `PCDiagonalScaleRight()`, `PCDiagonalScale()`
210: @*/
211: PetscErrorCode PCDiagonalScaleLeft(PC pc, Vec in, Vec out)
212: {
216:   if (pc->diagonalscale) {
217:     VecPointwiseMult(out, pc->diagonalscaleleft, in);
218:   } else if (in != out) {
219:     VecCopy(in, out);
220:   }
221:   return 0;
222: }

224: /*@
225:    PCDiagonalScaleRight - Scales a vector by the right scaling as needed by certain time-stepping codes.

227:    Logically Collective

229:    Input Parameters:
230: +  pc - the preconditioner context
231: .  in - input vector
232: -  out - scaled vector (maybe the same as in)

234:    Level: intermediate

236:    Notes:
237:     The system solved via the Krylov method is
238: .vb
239:         D M A D^{-1} y = D M b  for left preconditioning or
240:         D A M D^{-1} z = D b for right preconditioning
241: .ve

243:    `PCDiagonalScaleLeft()` scales a vector by D. `PCDiagonalScaleRight()` scales a vector by D^{-1}.

245:    If diagonal scaling is turned off and in is not out then in is copied to out

247: .seealso: `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCDiagonalScaleSet()`, `PCDiagonalScale()`
248: @*/
249: PetscErrorCode PCDiagonalScaleRight(PC pc, Vec in, Vec out)
250: {
254:   if (pc->diagonalscale) {
255:     VecPointwiseMult(out, pc->diagonalscaleright, in);
256:   } else if (in != out) {
257:     VecCopy(in, out);
258:   }
259:   return 0;
260: }

262: /*@
263:    PCSetUseAmat - Sets a flag to indicate that when the preconditioner needs to apply (part of) the
264:    operator during the preconditioning process it applies the Amat provided to `TSSetRHSJacobian()`,
265:    `TSSetIJacobian()`, `SNESSetJacobian()`, `KSPSetOperators()` or `PCSetOperators()` not the Pmat.

267:    Logically Collective

269:    Input Parameters:
270: +  pc - the preconditioner context
271: -  flg - `PETSC_TRUE` to use the Amat, `PETSC_FALSE` to use the Pmat (default is false)

273:    Options Database Key:
274: .  -pc_use_amat <true,false> - use the amat to apply the operator

276:    Note:
277:    For the common case in which the linear system matrix and the matrix used to construct the
278:    preconditioner are identical, this routine is does nothing.

280:    Level: intermediate

282: .seealso: `PC`, `PCGetUseAmat()`, `PCBJACOBI`, `PGMG`, `PCFIELDSPLIT`, `PCCOMPOSITE`
283: @*/
284: PetscErrorCode PCSetUseAmat(PC pc, PetscBool flg)
285: {
287:   pc->useAmat = flg;
288:   return 0;
289: }

291: /*@
292:    PCSetErrorIfFailure - Causes `PC` to generate an error if a FPE, for example a zero pivot, is detected.

294:    Logically Collective

296:    Input Parameters:
297: +  pc - iterative context obtained from PCCreate()
298: -  flg - `PETSC_TRUE` indicates you want the error generated

300:    Level: advanced

302:    Notes:
303:     Normally PETSc continues if a linear solver fails due to a failed setup of a preconditioner, you can call `KSPGetConvergedReason()` after a `KSPSolve()`
304:     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
305:     detected.

307:     This is propagated into KSPs used by this PC, which then propagate it into PCs used by those KSPs

309: .seealso: `PC`, `KSPSetErrorIfNotConverged()`, `PCGetInitialGuessNonzero()`, `PCSetInitialGuessKnoll()`, `PCGetInitialGuessKnoll()`
310: @*/
311: PetscErrorCode PCSetErrorIfFailure(PC pc, PetscBool flg)
312: {
315:   pc->erroriffailure = flg;
316:   return 0;
317: }

319: /*@
320:    PCGetUseAmat - Gets a flag to indicate that when the preconditioner needs to apply (part of) the
321:    operator during the preconditioning process it applies the Amat provided to `TSSetRHSJacobian()`,
322:    `TSSetIJacobian()`, `SNESSetJacobian()`, `KSPSetOperators()` or `PCSetOperators()` not the Pmat.

324:    Logically Collective

326:    Input Parameter:
327: .  pc - the preconditioner context

329:    Output Parameter:
330: .  flg - `PETSC_TRUE` to use the Amat, `PETSC_FALSE` to use the Pmat (default is false)

332:    Note:
333:    For the common case in which the linear system matrix and the matrix used to construct the
334:    preconditioner are identical, this routine is does nothing.

336:    Level: intermediate

338: .seealso: `PC`, `PCSetUseAmat()`, `PCBJACOBI`, `PGMG`, `PCFIELDSPLIT`, `PCCOMPOSITE`
339: @*/
340: PetscErrorCode PCGetUseAmat(PC pc, PetscBool *flg)
341: {
343:   *flg = pc->useAmat;
344:   return 0;
345: }

347: /*@
348:    PCCreate - Creates a preconditioner context, `PC`

350:    Collective

352:    Input Parameter:
353: .  comm - MPI communicator

355:    Output Parameter:
356: .  pc - location to put the preconditioner context

358:    Note:
359:    The default preconditioner for sparse matrices is `PCILU` or `PCICC` with 0 fill on one process and block Jacobi (`PCBJACOBI`) with `PCILU` or `PCICC`
360:    in parallel. For dense matrices it is always `PCNONE`.

362:    Level: developer

364: .seealso: `PC`, `PCSetUp()`, `PCApply()`, `PCDestroy()`
365: @*/
366: PetscErrorCode PCCreate(MPI_Comm comm, PC *newpc)
367: {
368:   PC pc;

371:   *newpc = NULL;
372:   PCInitializePackage();

374:   PetscHeaderCreate(pc, PC_CLASSID, "PC", "Preconditioner", "PC", comm, PCDestroy, PCView);

376:   pc->mat                  = NULL;
377:   pc->pmat                 = NULL;
378:   pc->setupcalled          = 0;
379:   pc->setfromoptionscalled = 0;
380:   pc->data                 = NULL;
381:   pc->diagonalscale        = PETSC_FALSE;
382:   pc->diagonalscaleleft    = NULL;
383:   pc->diagonalscaleright   = NULL;

385:   pc->modifysubmatrices  = NULL;
386:   pc->modifysubmatricesP = NULL;

388:   *newpc = pc;
389:   return 0;
390: }

392: /*@
393:    PCApply - Applies the preconditioner to a vector.

395:    Collective

397:    Input Parameters:
398: +  pc - the preconditioner context
399: -  x - input vector

401:    Output Parameter:
402: .  y - output vector

404:    Level: developer

406: .seealso: `PC`, `PCApplyTranspose()`, `PCApplyBAorAB()`
407: @*/
408: PetscErrorCode PCApply(PC pc, Vec x, Vec y)
409: {
410:   PetscInt m, n, mv, nv;

416:   if (pc->erroriffailure) VecValidValues_Internal(x, 2, PETSC_TRUE);
417:   /* use pmat to check vector sizes since for KSPLSQR the pmat may be of a different size than mat */
418:   MatGetLocalSize(pc->pmat, &m, &n);
419:   VecGetLocalSize(x, &mv);
420:   VecGetLocalSize(y, &nv);
421:   /* check pmat * y = x is feasible */
424:   VecSetErrorIfLocked(y, 3);

426:   PCSetUp(pc);
427:   VecLockReadPush(x);
428:   PetscLogEventBegin(PC_Apply, pc, x, y, 0);
429:   PetscUseTypeMethod(pc, apply, x, y);
430:   PetscLogEventEnd(PC_Apply, pc, x, y, 0);
431:   if (pc->erroriffailure) VecValidValues_Internal(y, 3, PETSC_FALSE);
432:   VecLockReadPop(x);
433:   return 0;
434: }

436: /*@
437:    PCMatApply - Applies the preconditioner to multiple vectors stored as a `MATDENSE`. Like `PCApply()`, Y and X must be different matrices.

439:    Collective

441:    Input Parameters:
442: +  pc - the preconditioner context
443: -  X - block of input vectors

445:    Output Parameter:
446: .  Y - block of output vectors

448:    Level: developer

450: .seealso: `PC`, `PCApply()`, `KSPMatSolve()`
451: @*/
452: PetscErrorCode PCMatApply(PC pc, Mat X, Mat Y)
453: {
454:   Mat       A;
455:   Vec       cy, cx;
456:   PetscInt  m1, M1, m2, M2, n1, N1, n2, N2, m3, M3, n3, N3;
457:   PetscBool match;

465:   PCGetOperators(pc, NULL, &A);
466:   MatGetLocalSize(A, &m3, &n3);
467:   MatGetLocalSize(X, &m2, &n2);
468:   MatGetLocalSize(Y, &m1, &n1);
469:   MatGetSize(A, &M3, &N3);
470:   MatGetSize(X, &M2, &N2);
471:   MatGetSize(Y, &M1, &N1);
475:   PetscObjectBaseTypeCompareAny((PetscObject)Y, &match, MATSEQDENSE, MATMPIDENSE, "");
477:   PetscObjectBaseTypeCompareAny((PetscObject)X, &match, MATSEQDENSE, MATMPIDENSE, "");
479:   PCSetUp(pc);
480:   if (pc->ops->matapply) {
481:     PetscLogEventBegin(PC_MatApply, pc, X, Y, 0);
482:     PetscUseTypeMethod(pc, matapply, X, Y);
483:     PetscLogEventEnd(PC_MatApply, pc, X, Y, 0);
484:   } else {
485:     PetscInfo(pc, "PC type %s applying column by column\n", ((PetscObject)pc)->type_name);
486:     for (n1 = 0; n1 < N1; ++n1) {
487:       MatDenseGetColumnVecRead(X, n1, &cx);
488:       MatDenseGetColumnVecWrite(Y, n1, &cy);
489:       PCApply(pc, cx, cy);
490:       MatDenseRestoreColumnVecWrite(Y, n1, &cy);
491:       MatDenseRestoreColumnVecRead(X, n1, &cx);
492:     }
493:   }
494:   return 0;
495: }

497: /*@
498:    PCApplySymmetricLeft - Applies the left part of a symmetric preconditioner to a vector.

500:    Collective

502:    Input Parameters:
503: +  pc - the preconditioner context
504: -  x - input vector

506:    Output Parameter:
507: .  y - output vector

509:    Note:
510:    Currently, this routine is implemented only for `PCICC` and `PCJACOBI` preconditioners.

512:    Level: developer

514: .seealso: `PC`, `PCApply()`, `PCApplySymmetricRight()`
515: @*/
516: PetscErrorCode PCApplySymmetricLeft(PC pc, Vec x, Vec y)
517: {
522:   if (pc->erroriffailure) VecValidValues_Internal(x, 2, PETSC_TRUE);
523:   PCSetUp(pc);
524:   VecLockReadPush(x);
525:   PetscLogEventBegin(PC_ApplySymmetricLeft, pc, x, y, 0);
526:   PetscUseTypeMethod(pc, applysymmetricleft, x, y);
527:   PetscLogEventEnd(PC_ApplySymmetricLeft, pc, x, y, 0);
528:   VecLockReadPop(x);
529:   if (pc->erroriffailure) VecValidValues_Internal(y, 3, PETSC_FALSE);
530:   return 0;
531: }

533: /*@
534:    PCApplySymmetricRight - Applies the right part of a symmetric preconditioner to a vector.

536:    Collective

538:    Input Parameters:
539: +  pc - the preconditioner context
540: -  x - input vector

542:    Output Parameter:
543: .  y - output vector

545:    Level: developer

547:    Note:
548:    Currently, this routine is implemented only for `PCICC` and `PCJACOBI` preconditioners.

550: .seealso: `PC`, `PCApply()`, `PCApplySymmetricLeft()`
551: @*/
552: PetscErrorCode PCApplySymmetricRight(PC pc, Vec x, Vec y)
553: {
558:   if (pc->erroriffailure) VecValidValues_Internal(x, 2, PETSC_TRUE);
559:   PCSetUp(pc);
560:   VecLockReadPush(x);
561:   PetscLogEventBegin(PC_ApplySymmetricRight, pc, x, y, 0);
562:   PetscUseTypeMethod(pc, applysymmetricright, x, y);
563:   PetscLogEventEnd(PC_ApplySymmetricRight, pc, x, y, 0);
564:   VecLockReadPop(x);
565:   if (pc->erroriffailure) VecValidValues_Internal(y, 3, PETSC_FALSE);
566:   return 0;
567: }

569: /*@
570:    PCApplyTranspose - Applies the transpose of preconditioner to a vector.

572:    Collective

574:    Input Parameters:
575: +  pc - the preconditioner context
576: -  x - input vector

578:    Output Parameter:
579: .  y - output vector

581:    Note:
582:     For complex numbers this applies the non-Hermitian transpose.

584:    Developer Note:
585:     We need to implement a `PCApplyHermitianTranspose()`

587:    Level: developer

589: .seealso: `PC`, `PCApply()`, `PCApplyBAorAB()`, `PCApplyBAorABTranspose()`, `PCApplyTransposeExists()`
590: @*/
591: PetscErrorCode PCApplyTranspose(PC pc, Vec x, Vec y)
592: {
597:   if (pc->erroriffailure) VecValidValues_Internal(x, 2, PETSC_TRUE);
598:   PCSetUp(pc);
599:   VecLockReadPush(x);
600:   PetscLogEventBegin(PC_Apply, pc, x, y, 0);
601:   PetscUseTypeMethod(pc, applytranspose, x, y);
602:   PetscLogEventEnd(PC_Apply, pc, x, y, 0);
603:   VecLockReadPop(x);
604:   if (pc->erroriffailure) VecValidValues_Internal(y, 3, PETSC_FALSE);
605:   return 0;
606: }

608: /*@
609:    PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation

611:    Collective

613:    Input Parameter:
614: .  pc - the preconditioner context

616:    Output Parameter:
617: .  flg - `PETSC_TRUE` if a transpose operation is defined

619:    Level: developer

621: .seealso: `PC`, `PCApplyTranspose()`
622: @*/
623: PetscErrorCode PCApplyTransposeExists(PC pc, PetscBool *flg)
624: {
627:   if (pc->ops->applytranspose) *flg = PETSC_TRUE;
628:   else *flg = PETSC_FALSE;
629:   return 0;
630: }

632: /*@
633:    PCApplyBAorAB - Applies the preconditioner and operator to a vector. y = B*A*x or y = A*B*x.

635:    Collective

637:    Input Parameters:
638: +  pc - the preconditioner context
639: .  side - indicates the preconditioner side, one of `PC_LEFT`, `PC_RIGHT`, or `PC_SYMMETRIC`
640: .  x - input vector
641: -  work - work vector

643:    Output Parameter:
644: .  y - output vector

646:    Level: developer

648:    Note:
649:     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
650:    specific `KSPSolve()` method must also be written to handle the post-solve "correction" for the diagonal scaling.

652: .seealso: `PC`, `PCApply()`, `PCApplyTranspose()`, `PCApplyBAorABTranspose()`
653: @*/
654: PetscErrorCode PCApplyBAorAB(PC pc, PCSide side, Vec x, Vec y, Vec work)
655: {
667:   if (pc->erroriffailure) VecValidValues_Internal(x, 3, PETSC_TRUE);

669:   PCSetUp(pc);
670:   if (pc->diagonalscale) {
671:     if (pc->ops->applyBA) {
672:       Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
673:       VecDuplicate(x, &work2);
674:       PCDiagonalScaleRight(pc, x, work2);
675:       PetscUseTypeMethod(pc, applyBA, side, work2, y, work);
676:       PCDiagonalScaleLeft(pc, y, y);
677:       VecDestroy(&work2);
678:     } else if (side == PC_RIGHT) {
679:       PCDiagonalScaleRight(pc, x, y);
680:       PCApply(pc, y, work);
681:       MatMult(pc->mat, work, y);
682:       PCDiagonalScaleLeft(pc, y, y);
683:     } else if (side == PC_LEFT) {
684:       PCDiagonalScaleRight(pc, x, y);
685:       MatMult(pc->mat, y, work);
686:       PCApply(pc, work, y);
687:       PCDiagonalScaleLeft(pc, y, y);
689:   } else {
690:     if (pc->ops->applyBA) {
691:       PetscUseTypeMethod(pc, applyBA, side, x, y, work);
692:     } else if (side == PC_RIGHT) {
693:       PCApply(pc, x, work);
694:       MatMult(pc->mat, work, y);
695:     } else if (side == PC_LEFT) {
696:       MatMult(pc->mat, x, work);
697:       PCApply(pc, work, y);
698:     } else if (side == PC_SYMMETRIC) {
699:       /* There's an extra copy here; maybe should provide 2 work vectors instead? */
700:       PCApplySymmetricRight(pc, x, work);
701:       MatMult(pc->mat, work, y);
702:       VecCopy(y, work);
703:       PCApplySymmetricLeft(pc, work, y);
704:     }
705:   }
706:   if (pc->erroriffailure) VecValidValues_Internal(y, 4, PETSC_FALSE);
707:   return 0;
708: }

710: /*@
711:    PCApplyBAorABTranspose - Applies the transpose of the preconditioner
712:    and operator to a vector. That is, applies tr(B) * tr(A) with left preconditioning,
713:    NOT tr(B*A) = tr(A)*tr(B).

715:    Collective

717:    Input Parameters:
718: +  pc - the preconditioner context
719: .  side - indicates the preconditioner side, one of `PC_LEFT`, `PC_RIGHT`, or `PC_SYMMETRIC`
720: .  x - input vector
721: -  work - work vector

723:    Output Parameter:
724: .  y - output vector

726:    Note:
727:     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
728:       defined by B'. This is why this has the funny form that it computes tr(B) * tr(A)

730:     Level: developer

732: .seealso: `PC`, `PCApply()`, `PCApplyTranspose()`, `PCApplyBAorAB()`
733: @*/
734: PetscErrorCode PCApplyBAorABTranspose(PC pc, PCSide side, Vec x, Vec y, Vec work)
735: {
741:   if (pc->erroriffailure) VecValidValues_Internal(x, 3, PETSC_TRUE);
742:   if (pc->ops->applyBAtranspose) {
743:     PetscUseTypeMethod(pc, applyBAtranspose, side, x, y, work);
744:     if (pc->erroriffailure) VecValidValues_Internal(y, 4, PETSC_FALSE);
745:     return 0;
746:   }

749:   PCSetUp(pc);
750:   if (side == PC_RIGHT) {
751:     PCApplyTranspose(pc, x, work);
752:     MatMultTranspose(pc->mat, work, y);
753:   } else if (side == PC_LEFT) {
754:     MatMultTranspose(pc->mat, x, work);
755:     PCApplyTranspose(pc, work, y);
756:   }
757:   /* add support for PC_SYMMETRIC */
758:   if (pc->erroriffailure) VecValidValues_Internal(y, 4, PETSC_FALSE);
759:   return 0;
760: }

762: /*@
763:    PCApplyRichardsonExists - Determines whether a particular preconditioner has a
764:    built-in fast application of Richardson's method.

766:    Not Collective

768:    Input Parameter:
769: .  pc - the preconditioner

771:    Output Parameter:
772: .  exists - `PETSC_TRUE` or `PETSC_FALSE`

774:    Level: developer

776: .seealso: `PC`, `PCRICHARDSON`, `PCApplyRichardson()`
777: @*/
778: PetscErrorCode PCApplyRichardsonExists(PC pc, PetscBool *exists)
779: {
782:   if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
783:   else *exists = PETSC_FALSE;
784:   return 0;
785: }

787: /*@
788:    PCApplyRichardson - Applies several steps of Richardson iteration with
789:    the particular preconditioner. This routine is usually used by the
790:    Krylov solvers and not the application code directly.

792:    Collective

794:    Input Parameters:
795: +  pc  - the preconditioner context
796: .  b   - the right hand side
797: .  w   - one work vector
798: .  rtol - relative decrease in residual norm convergence criteria
799: .  abstol - absolute residual norm convergence criteria
800: .  dtol - divergence residual norm increase criteria
801: .  its - the number of iterations to apply.
802: -  guesszero - if the input x contains nonzero initial guess

804:    Output Parameters:
805: +  outits - number of iterations actually used (for SOR this always equals its)
806: .  reason - the reason the apply terminated
807: -  y - the solution (also contains initial guess if guesszero is `PETSC_FALSE`

809:    Notes:
810:    Most preconditioners do not support this function. Use the command
811:    `PCApplyRichardsonExists()` to determine if one does.

813:    Except for the `PCMG` this routine ignores the convergence tolerances
814:    and always runs for the number of iterations

816:    Level: developer

818: .seealso: `PC`, `PCApplyRichardsonExists()`
819: @*/
820: PetscErrorCode PCApplyRichardson(PC pc, Vec b, Vec y, Vec w, PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt its, PetscBool guesszero, PetscInt *outits, PCRichardsonConvergedReason *reason)
821: {
827:   PCSetUp(pc);
828:   PetscUseTypeMethod(pc, applyrichardson, b, y, w, rtol, abstol, dtol, its, guesszero, outits, reason);
829:   return 0;
830: }

832: /*@
833:    PCSetFailedReason - Sets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail

835:    Logically Collective

837:    Input Parameters:
838: +  pc - the preconditioner context
839: -  reason - the reason it failedx

841:    Level: advanced

843: .seealso: `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCFailedReason`
844: @*/
845: PetscErrorCode PCSetFailedReason(PC pc, PCFailedReason reason)
846: {
847:   pc->failedreason = reason;
848:   return 0;
849: }

851: /*@
852:    PCGetFailedReason - Gets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail

854:    Logically Collective

856:    Input Parameter:
857: .  pc - the preconditioner context

859:    Output Parameter:
860: .  reason - the reason it failed

862:    Level: advanced

864:    Note:
865:    This is the maximum over reason over all ranks in the PC communicator. It is only valid after
866:    a call `KSPCheckDot()` or  `KSPCheckNorm()` inside a `KSPSolve()`. It is not valid immediately after a `PCSetUp()`
867:    or `PCApply()`, then use `PCGetFailedReasonRank()`

869: .seealso: PC`, ``PCCreate()`, `PCApply()`, `PCDestroy()`, `PCGetFailedReasonRank()`, `PCSetFailedReason()`
870: @*/
871: PetscErrorCode PCGetFailedReason(PC pc, PCFailedReason *reason)
872: {
873:   if (pc->setupcalled < 0) *reason = (PCFailedReason)pc->setupcalled;
874:   else *reason = pc->failedreason;
875:   return 0;
876: }

878: /*@
879:    PCGetFailedReasonRank - Gets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail on this MPI rank

881:    Not Collective

883:    Input Parameter:
884: .  pc - the preconditioner context

886:    Output Parameter:
887: .  reason - the reason it failed

889:    Note:
890:      Different ranks may have different reasons or no reason, see `PCGetFailedReason()`

892:    Level: advanced

894: .seealso: `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCGetFailedReason()`, `PCSetFailedReason()`
895: @*/
896: PetscErrorCode PCGetFailedReasonRank(PC pc, PCFailedReason *reason)
897: {
898:   if (pc->setupcalled < 0) *reason = (PCFailedReason)pc->setupcalled;
899:   else *reason = pc->failedreason;
900:   return 0;
901: }

903: /*  Next line needed to deactivate KSP_Solve logging */
904: #include <petsc/private/kspimpl.h>

906: /*
907:       a setupcall of 0 indicates never setup,
908:                      1 indicates has been previously setup
909:                     -1 indicates a PCSetUp() was attempted and failed
910: */
911: /*@
912:    PCSetUp - Prepares for the use of a preconditioner.

914:    Collective

916:    Input Parameter:
917: .  pc - the preconditioner context

919:    Level: developer

921: .seealso: `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`
922: @*/
923: PetscErrorCode PCSetUp(PC pc)
924: {
925:   const char      *def;
926:   PetscObjectState matstate, matnonzerostate;


931:   if (pc->setupcalled && pc->reusepreconditioner) {
932:     PetscInfo(pc, "Leaving PC with identical preconditioner since reuse preconditioner is set\n");
933:     return 0;
934:   }

936:   PetscObjectStateGet((PetscObject)pc->pmat, &matstate);
937:   MatGetNonzeroState(pc->pmat, &matnonzerostate);
938:   if (!pc->setupcalled) {
939:     PetscInfo(pc, "Setting up PC for first time\n");
940:     pc->flag = DIFFERENT_NONZERO_PATTERN;
941:   } else if (matstate == pc->matstate) {
942:     PetscInfo(pc, "Leaving PC with identical preconditioner since operator is unchanged\n");
943:     return 0;
944:   } else {
945:     if (matnonzerostate > pc->matnonzerostate) {
946:       PetscInfo(pc, "Setting up PC with different nonzero pattern\n");
947:       pc->flag = DIFFERENT_NONZERO_PATTERN;
948:     } else {
949:       PetscInfo(pc, "Setting up PC with same nonzero pattern\n");
950:       pc->flag = SAME_NONZERO_PATTERN;
951:     }
952:   }
953:   pc->matstate        = matstate;
954:   pc->matnonzerostate = matnonzerostate;

956:   if (!((PetscObject)pc)->type_name) {
957:     PCGetDefaultType_Private(pc, &def);
958:     PCSetType(pc, def);
959:   }

961:   MatSetErrorIfFailure(pc->pmat, pc->erroriffailure);
962:   MatSetErrorIfFailure(pc->mat, pc->erroriffailure);
963:   PetscLogEventBegin(PC_SetUp, pc, 0, 0, 0);
964:   if (pc->ops->setup) {
965:     /* do not log solves and applications of preconditioners while constructing preconditioners; perhaps they should be logged separately from the regular solves */
966:     KSPInitializePackage();
967:     PetscLogEventDeactivatePush(KSP_Solve);
968:     PetscLogEventDeactivatePush(PC_Apply);
969:     PetscUseTypeMethod(pc, setup);
970:     PetscLogEventDeactivatePop(KSP_Solve);
971:     PetscLogEventDeactivatePop(PC_Apply);
972:   }
973:   PetscLogEventEnd(PC_SetUp, pc, 0, 0, 0);
974:   if (!pc->setupcalled) pc->setupcalled = 1;
975:   return 0;
976: }

978: /*@
979:    PCSetUpOnBlocks - Sets up the preconditioner for each block in
980:    the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
981:    methods.

983:    Collective

985:    Input Parameter:
986: .  pc - the preconditioner context

988:    Level: developer

990:    Note:
991:    For nested preconditioners such as `PCBJACOBI` `PCSetUp()` is not called on each sub-`KSP` when `PCSetUp()` is
992:    called on the outer `PC`, this routine ensures it is called.

994: .seealso: `PC`, `PCSetUp()`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCSetUp()`
995: @*/
996: PetscErrorCode PCSetUpOnBlocks(PC pc)
997: {
999:   if (!pc->ops->setuponblocks) return 0;
1000:   PetscLogEventBegin(PC_SetUpOnBlocks, pc, 0, 0, 0);
1001:   PetscUseTypeMethod(pc, setuponblocks);
1002:   PetscLogEventEnd(PC_SetUpOnBlocks, pc, 0, 0, 0);
1003:   return 0;
1004: }

1006: /*@C
1007:    PCSetModifySubMatrices - Sets a user-defined routine for modifying the
1008:    submatrices that arise within certain subdomain-based preconditioners.
1009:    The basic submatrices are extracted from the preconditioner matrix as
1010:    usual; the user can then alter these (for example, to set different boundary
1011:    conditions for each submatrix) before they are used for the local solves.

1013:    Logically Collective

1015:    Input Parameters:
1016: +  pc - the preconditioner context
1017: .  func - routine for modifying the submatrices
1018: -  ctx - optional user-defined context (may be null)

1020:    Calling sequence of func:
1021: $     func (PC pc,PetscInt nsub,IS *row,IS *col,Mat *submat,void *ctx);

1023: +  row - an array of index sets that contain the global row numbers
1024:          that comprise each local submatrix
1025: .  col - an array of index sets that contain the global column numbers
1026:          that comprise each local submatrix
1027: .  submat - array of local submatrices
1028: -  ctx - optional user-defined context for private data for the
1029:          user-defined func routine (may be null)

1031:    Notes:
1032:    `PCSetModifySubMatrices()` MUST be called before `KSPSetUp()` and
1033:    `KSPSolve()`.

1035:    A routine set by `PCSetModifySubMatrices()` is currently called within
1036:    the block Jacobi (`PCBJACOBI`) and additive Schwarz (`PCASM`)
1037:    preconditioners.  All other preconditioners ignore this routine.

1039:    Level: advanced

1041: .seealso: `PC`, `PCBJACOBI`, `PCASM`, `PCModifySubMatrices()`
1042: @*/
1043: PetscErrorCode PCSetModifySubMatrices(PC pc, PetscErrorCode (*func)(PC, PetscInt, const IS[], const IS[], Mat[], void *), void *ctx)
1044: {
1046:   pc->modifysubmatrices  = func;
1047:   pc->modifysubmatricesP = ctx;
1048:   return 0;
1049: }

1051: /*@C
1052:    PCModifySubMatrices - Calls an optional user-defined routine within
1053:    certain preconditioners if one has been set with `PCSetModifySubMatrices()`.

1055:    Collective

1057:    Input Parameters:
1058: +  pc - the preconditioner context
1059: .  nsub - the number of local submatrices
1060: .  row - an array of index sets that contain the global row numbers
1061:          that comprise each local submatrix
1062: .  col - an array of index sets that contain the global column numbers
1063:          that comprise each local submatrix
1064: .  submat - array of local submatrices
1065: -  ctx - optional user-defined context for private data for the
1066:          user-defined routine (may be null)

1068:    Output Parameter:
1069: .  submat - array of local submatrices (the entries of which may
1070:             have been modified)

1072:    Notes:
1073:    The user should NOT generally call this routine, as it will
1074:    automatically be called within certain preconditioners (currently
1075:    block Jacobi, additive Schwarz) if set.

1077:    The basic submatrices are extracted from the preconditioner matrix
1078:    as usual; the user can then alter these (for example, to set different
1079:    boundary conditions for each submatrix) before they are used for the
1080:    local solves.

1082:    Level: developer

1084: .seealso: `PC`, `PCSetModifySubMatrices()`
1085: @*/
1086: PetscErrorCode PCModifySubMatrices(PC pc, PetscInt nsub, const IS row[], const IS col[], Mat submat[], void *ctx)
1087: {
1089:   if (!pc->modifysubmatrices) return 0;
1090:   PetscLogEventBegin(PC_ModifySubMatrices, pc, 0, 0, 0);
1091:   (*pc->modifysubmatrices)(pc, nsub, row, col, submat, ctx);
1092:   PetscLogEventEnd(PC_ModifySubMatrices, pc, 0, 0, 0);
1093:   return 0;
1094: }

1096: /*@
1097:    PCSetOperators - Sets the matrix associated with the linear system and
1098:    a (possibly) different one associated with the preconditioner.

1100:    Logically Collective

1102:    Input Parameters:
1103: +  pc - the preconditioner context
1104: .  Amat - the matrix that defines the linear system
1105: -  Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.

1107:    Notes:
1108:     Passing a NULL for Amat or Pmat removes the matrix that is currently used.

1110:     If you wish to replace either Amat or Pmat but leave the other one untouched then
1111:     first call `KSPGetOperators()` to get the one you wish to keep, call `PetscObjectReference()`
1112:     on it and then pass it back in in your call to `KSPSetOperators()`.

1114:    More Notes about Repeated Solution of Linear Systems:
1115:    PETSc does NOT reset the matrix entries of either Amat or Pmat
1116:    to zero after a linear solve; the user is completely responsible for
1117:    matrix assembly.  See the routine `MatZeroEntries()` if desiring to
1118:    zero all elements of a matrix.

1120:    Level: intermediate

1122: .seealso: `PC`, `PCGetOperators()`, `MatZeroEntries()`
1123:  @*/
1124: PetscErrorCode PCSetOperators(PC pc, Mat Amat, Mat Pmat)
1125: {
1126:   PetscInt m1, n1, m2, n2;

1133:   if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
1134:     MatGetLocalSize(Amat, &m1, &n1);
1135:     MatGetLocalSize(pc->mat, &m2, &n2);
1137:     MatGetLocalSize(Pmat, &m1, &n1);
1138:     MatGetLocalSize(pc->pmat, &m2, &n2);
1140:   }

1142:   if (Pmat != pc->pmat) {
1143:     /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
1144:     pc->matnonzerostate = -1;
1145:     pc->matstate        = -1;
1146:   }

1148:   /* reference first in case the matrices are the same */
1149:   if (Amat) PetscObjectReference((PetscObject)Amat);
1150:   MatDestroy(&pc->mat);
1151:   if (Pmat) PetscObjectReference((PetscObject)Pmat);
1152:   MatDestroy(&pc->pmat);
1153:   pc->mat  = Amat;
1154:   pc->pmat = Pmat;
1155:   return 0;
1156: }

1158: /*@
1159:    PCSetReusePreconditioner - reuse the current preconditioner even if the operator in the preconditioner has changed.

1161:    Logically Collective

1163:    Input Parameters:
1164: +  pc - the preconditioner context
1165: -  flag - `PETSC_TRUE` do not compute a new preconditioner, `PETSC_FALSE` do compute a new preconditioner

1167:     Level: intermediate

1169:    Note:
1170:    Normally if a matrix inside a `PC` changes the `PC` automatically updates itself using information from the changed matrix. This option
1171:    prevents this.

1173: .seealso: `PC`, `PCGetOperators()`, `MatZeroEntries()`, `PCGetReusePreconditioner()`, `KSPSetReusePreconditioner()`
1174:  @*/
1175: PetscErrorCode PCSetReusePreconditioner(PC pc, PetscBool flag)
1176: {
1179:   pc->reusepreconditioner = flag;
1180:   return 0;
1181: }

1183: /*@
1184:    PCGetReusePreconditioner - Determines if the `PC` reuses the current preconditioner even if the operator in the preconditioner has changed.

1186:    Not Collective

1188:    Input Parameter:
1189: .  pc - the preconditioner context

1191:    Output Parameter:
1192: .  flag - `PETSC_TRUE` do not compute a new preconditioner, `PETSC_FALSE` do compute a new preconditioner

1194:    Level: intermediate

1196: .seealso: `PC`, `PCGetOperators()`, `MatZeroEntries()`, `PCSetReusePreconditioner()`
1197:  @*/
1198: PetscErrorCode PCGetReusePreconditioner(PC pc, PetscBool *flag)
1199: {
1202:   *flag = pc->reusepreconditioner;
1203:   return 0;
1204: }

1206: /*@
1207:    PCGetOperators - Gets the matrix associated with the linear system and
1208:    possibly a different one associated with the preconditioner.

1210:    Not collective, though parallel mats are returned if pc is parallel

1212:    Input Parameter:
1213: .  pc - the preconditioner context

1215:    Output Parameters:
1216: +  Amat - the matrix defining the linear system
1217: -  Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.

1219:    Level: intermediate

1221:    Note:
1222:     Does not increase the reference count of the matrices, so you should not destroy them

1224:    Alternative usage: If the operators have NOT been set with `KSPSetOperators()`/`PCSetOperators()` then the operators
1225:       are created in `PC` and returned to the user. In this case, if both operators
1226:       mat and pmat are requested, two DIFFERENT operators will be returned. If
1227:       only one is requested both operators in the PC will be the same (i.e. as
1228:       if one had called `KSPSetOperators()`/`PCSetOperators()` with the same argument for both Mats).
1229:       The user must set the sizes of the returned matrices and their type etc just
1230:       as if the user created them with `MatCreate()`. For example,

1232: .vb
1233:          KSP/PCGetOperators(ksp/pc,&Amat,NULL); is equivalent to
1234:            set size, type, etc of Amat

1236:          MatCreate(comm,&mat);
1237:          KSP/PCSetOperators(ksp/pc,Amat,Amat);
1238:          PetscObjectDereference((PetscObject)mat);
1239:            set size, type, etc of Amat
1240: .ve

1242:      and

1244: .vb
1245:          KSP/PCGetOperators(ksp/pc,&Amat,&Pmat); is equivalent to
1246:            set size, type, etc of Amat and Pmat

1248:          MatCreate(comm,&Amat);
1249:          MatCreate(comm,&Pmat);
1250:          KSP/PCSetOperators(ksp/pc,Amat,Pmat);
1251:          PetscObjectDereference((PetscObject)Amat);
1252:          PetscObjectDereference((PetscObject)Pmat);
1253:            set size, type, etc of Amat and Pmat
1254: .ve

1256:     The rationale for this support is so that when creating a `TS`, `SNES`, or `KSP` the hierarchy
1257:     of underlying objects (i.e. `SNES`, `KSP`, `PC`, `Mat`) and their lifespans can be completely
1258:     managed by the top most level object (i.e. the `TS`, `SNES`, or `KSP`). Another way to look
1259:     at this is when you create a `SNES` you do not NEED to create a `KSP` and attach it to
1260:     the `SNES` object (the `SNES` object manages it for you). Similarly when you create a KSP
1261:     you do not need to attach a `PC` to it (the `KSP` object manages the `PC` object for you).
1262:     Thus, why should YOU have to create the `Mat` and attach it to the `SNES`/`KSP`/`PC`, when
1263:     it can be created for you?

1265: .seealso: `PC`, `PCSetOperators()`, `KSPGetOperators()`, `KSPSetOperators()`, `PCGetOperatorsSet()`
1266: @*/
1267: PetscErrorCode PCGetOperators(PC pc, Mat *Amat, Mat *Pmat)
1268: {
1270:   if (Amat) {
1271:     if (!pc->mat) {
1272:       if (pc->pmat && !Pmat) { /* Apmat has been set, but user did not request it, so use for Amat */
1273:         pc->mat = pc->pmat;
1274:         PetscObjectReference((PetscObject)pc->mat);
1275:       } else { /* both Amat and Pmat are empty */
1276:         MatCreate(PetscObjectComm((PetscObject)pc), &pc->mat);
1277:         if (!Pmat) { /* user did NOT request Pmat, so make same as Amat */
1278:           pc->pmat = pc->mat;
1279:           PetscObjectReference((PetscObject)pc->pmat);
1280:         }
1281:       }
1282:     }
1283:     *Amat = pc->mat;
1284:   }
1285:   if (Pmat) {
1286:     if (!pc->pmat) {
1287:       if (pc->mat && !Amat) { /* Amat has been set but was not requested, so use for pmat */
1288:         pc->pmat = pc->mat;
1289:         PetscObjectReference((PetscObject)pc->pmat);
1290:       } else {
1291:         MatCreate(PetscObjectComm((PetscObject)pc), &pc->pmat);
1292:         if (!Amat) { /* user did NOT request Amat, so make same as Pmat */
1293:           pc->mat = pc->pmat;
1294:           PetscObjectReference((PetscObject)pc->mat);
1295:         }
1296:       }
1297:     }
1298:     *Pmat = pc->pmat;
1299:   }
1300:   return 0;
1301: }

1303: /*@C
1304:    PCGetOperatorsSet - Determines if the matrix associated with the linear system and
1305:    possibly a different one associated with the preconditioner have been set in the `PC`.

1307:    Not collective, though the results on all processes should be the same

1309:    Input Parameter:
1310: .  pc - the preconditioner context

1312:    Output Parameters:
1313: +  mat - the matrix associated with the linear system was set
1314: -  pmat - matrix associated with the preconditioner was set, usually the same

1316:    Level: intermediate

1318: .seealso: `PC`, `PCSetOperators()`, `KSPGetOperators()`, `KSPSetOperators()`, `PCGetOperators()`
1319: @*/
1320: PetscErrorCode PCGetOperatorsSet(PC pc, PetscBool *mat, PetscBool *pmat)
1321: {
1323:   if (mat) *mat = (pc->mat) ? PETSC_TRUE : PETSC_FALSE;
1324:   if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1325:   return 0;
1326: }

1328: /*@
1329:    PCFactorGetMatrix - Gets the factored matrix from the
1330:    preconditioner context.  This routine is valid only for the `PCLU`,
1331:    `PCILU`, `PCCHOLESKY`, and `PCICC` methods.

1333:    Not Collective though mat is parallel if pc is parallel

1335:    Input Parameter:
1336: .  pc - the preconditioner context

1338:    Output parameters:
1339: .  mat - the factored matrix

1341:    Level: advanced

1343:    Note:
1344:     Does not increase the reference count for the matrix so DO NOT destroy it

1346: .seealso: `PC`, `PCLU`, `PCILU`, `PCCHOLESKY`, `PCICC`
1347: @*/
1348: PetscErrorCode PCFactorGetMatrix(PC pc, Mat *mat)
1349: {
1352:   PetscUseTypeMethod(pc, getfactoredmatrix, mat);
1353:   return 0;
1354: }

1356: /*@C
1357:    PCSetOptionsPrefix - Sets the prefix used for searching for all
1358:    `PC` options in the database.

1360:    Logically Collective

1362:    Input Parameters:
1363: +  pc - the preconditioner context
1364: -  prefix - the prefix string to prepend to all `PC` option requests

1366:    Note:
1367:    A hyphen (-) must NOT be given at the beginning of the prefix name.
1368:    The first character of all runtime options is AUTOMATICALLY the
1369:    hyphen.

1371:    Level: advanced

1373: .seealso: `PC`, `PCSetFromOptions`, `PCAppendOptionsPrefix()`, `PCGetOptionsPrefix()`
1374: @*/
1375: PetscErrorCode PCSetOptionsPrefix(PC pc, const char prefix[])
1376: {
1378:   PetscObjectSetOptionsPrefix((PetscObject)pc, prefix);
1379:   return 0;
1380: }

1382: /*@C
1383:    PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1384:    `PC` options in the database.

1386:    Logically Collective

1388:    Input Parameters:
1389: +  pc - the preconditioner context
1390: -  prefix - the prefix string to prepend to all `PC` option requests

1392:    Note:
1393:    A hyphen (-) must NOT be given at the beginning of the prefix name.
1394:    The first character of all runtime options is AUTOMATICALLY the
1395:    hyphen.

1397:    Level: advanced

1399: .seealso: `PC`, `PCSetFromOptions`, `PCSetOptionsPrefix()`, `PCGetOptionsPrefix()`
1400: @*/
1401: PetscErrorCode PCAppendOptionsPrefix(PC pc, const char prefix[])
1402: {
1404:   PetscObjectAppendOptionsPrefix((PetscObject)pc, prefix);
1405:   return 0;
1406: }

1408: /*@C
1409:    PCGetOptionsPrefix - Gets the prefix used for searching for all
1410:    PC options in the database.

1412:    Not Collective

1414:    Input Parameter:
1415: .  pc - the preconditioner context

1417:    Output Parameter:
1418: .  prefix - pointer to the prefix string used, is returned

1420:    Fortran Note:
1421:    The user should pass in a string 'prefix' of
1422:    sufficient length to hold the prefix.

1424:    Level: advanced

1426: .seealso: `PC`, `PCSetFromOptions`, `PCSetOptionsPrefix()`, `PCAppendOptionsPrefix()`
1427: @*/
1428: PetscErrorCode PCGetOptionsPrefix(PC pc, const char *prefix[])
1429: {
1432:   PetscObjectGetOptionsPrefix((PetscObject)pc, prefix);
1433:   return 0;
1434: }

1436: /*
1437:    Indicates the right hand side will be changed by KSPSolve(), this occurs for a few
1438:   preconditioners including BDDC and Eisentat that transform the equations before applying
1439:   the Krylov methods
1440: */
1441: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC pc, PetscBool *change)
1442: {
1445:   *change = PETSC_FALSE;
1446:   PetscTryMethod(pc, "PCPreSolveChangeRHS_C", (PC, PetscBool *), (pc, change));
1447:   return 0;
1448: }

1450: /*@
1451:    PCPreSolve - Optional pre-solve phase, intended for any
1452:    preconditioner-specific actions that must be performed before
1453:    the iterative solve itself.

1455:    Collective

1457:    Input Parameters:
1458: +  pc - the preconditioner context
1459: -  ksp - the Krylov subspace context

1461:    Level: developer

1463:    Sample of Usage:
1464: .vb
1465:     PCPreSolve(pc,ksp);
1466:     KSPSolve(ksp,b,x);
1467:     PCPostSolve(pc,ksp);
1468: .ve

1470:    Notes:
1471:    The pre-solve phase is distinct from the `PCSetUp()` phase.

1473:    `KSPSolve()` calls this directly, so is rarely called by the user.

1475: .seealso: `PC`, `PCPostSolve()`
1476: @*/
1477: PetscErrorCode PCPreSolve(PC pc, KSP ksp)
1478: {
1479:   Vec x, rhs;

1483:   pc->presolvedone++;
1485:   KSPGetSolution(ksp, &x);
1486:   KSPGetRhs(ksp, &rhs);

1488:   if (pc->ops->presolve) PetscUseTypeMethod(pc, presolve, ksp, rhs, x);
1489:   else if (pc->presolve) (pc->presolve)(pc, ksp);
1490:   return 0;
1491: }

1493: /*@C
1494:    PCSetPreSolve - Sets function used by `PCPreSolve()` which is intended for any
1495:    preconditioner-specific actions that must be performed before
1496:    the iterative solve itself.

1498:    Logically Collective

1500:    Input Parameters:
1501: +   pc - the preconditioner object
1502: -   presolve - the function to call before the solve

1504:    Calling sequence of presolve:
1505: $  func(PC pc,KSP ksp)

1507: +  pc - the PC context
1508: -  ksp - the KSP context

1510:    Level: developer

1512: .seealso: `PC`, `PCSetUp()`, `PCPreSolve()`
1513: @*/
1514: PetscErrorCode PCSetPreSolve(PC pc, PetscErrorCode (*presolve)(PC, KSP))
1515: {
1517:   pc->presolve = presolve;
1518:   return 0;
1519: }

1521: /*@
1522:    PCPostSolve - Optional post-solve phase, intended for any
1523:    preconditioner-specific actions that must be performed after
1524:    the iterative solve itself.

1526:    Collective

1528:    Input Parameters:
1529: +  pc - the preconditioner context
1530: -  ksp - the Krylov subspace context

1532:    Sample of Usage:
1533: .vb
1534:     PCPreSolve(pc,ksp);
1535:     KSPSolve(ksp,b,x);
1536:     PCPostSolve(pc,ksp);
1537: .ve

1539:    Note:
1540:    `KSPSolve()` calls this routine directly, so it is rarely called by the user.

1542:    Level: developer

1544: .seealso: `PC`, `PCSetPostSolve()`, `PCSetPresolve()`, `PCPreSolve()`, `KSPSolve()`
1545: @*/
1546: PetscErrorCode PCPostSolve(PC pc, KSP ksp)
1547: {
1548:   Vec x, rhs;

1552:   pc->presolvedone--;
1553:   KSPGetSolution(ksp, &x);
1554:   KSPGetRhs(ksp, &rhs);
1555:   PetscTryTypeMethod(pc, postsolve, ksp, rhs, x);
1556:   return 0;
1557: }

1559: /*@C
1560:   PCLoad - Loads a `PC` that has been stored in binary  with `PCView()`.

1562:   Collective

1564:   Input Parameters:
1565: + newdm - the newly loaded `PC`, this needs to have been created with `PCCreate()` or
1566:            some related function before a call to `PCLoad()`.
1567: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()`

1569:    Level: intermediate

1571:   Note:
1572:    The type is determined by the data in the file, any type set into the PC before this call is ignored.

1574: .seealso: `PC`, `PetscViewerBinaryOpen()`, `PCView()`, `MatLoad()`, `VecLoad()`
1575: @*/
1576: PetscErrorCode PCLoad(PC newdm, PetscViewer viewer)
1577: {
1578:   PetscBool isbinary;
1579:   PetscInt  classid;
1580:   char      type[256];

1584:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary);

1587:   PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT);
1589:   PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR);
1590:   PCSetType(newdm, type);
1591:   PetscTryTypeMethod(newdm, load, viewer);
1592:   return 0;
1593: }

1595: #include <petscdraw.h>
1596: #if defined(PETSC_HAVE_SAWS)
1597: #include <petscviewersaws.h>
1598: #endif

1600: /*@C
1601:    PCViewFromOptions - View from the `PC` based on options in the database

1603:    Collective

1605:    Input Parameters:
1606: +  A - the PC context
1607: .  obj - Optional object that provides the options prefix
1608: -  name - command line option

1610:    Level: intermediate

1612: .seealso: `PC`, `PCView`, `PetscObjectViewFromOptions()`, `PCCreate()`
1613: @*/
1614: PetscErrorCode PCViewFromOptions(PC A, PetscObject obj, const char name[])
1615: {
1617:   PetscObjectViewFromOptions((PetscObject)A, obj, name);
1618:   return 0;
1619: }

1621: /*@C
1622:    PCView - Prints information about the `PC`

1624:    Collective

1626:    Input Parameters:
1627: +  PC - the `PC` context
1628: -  viewer - optional visualization context

1630:    Notes:
1631:    The available visualization contexts include
1632: +     `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
1633: -     `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
1634:          output where only the first processor opens
1635:          the file.  All other processors send their
1636:          data to the first processor to print.

1638:    The user can open an alternative visualization contexts with
1639:    `PetscViewerASCIIOpen()` (output to a specified file).

1641:    Level: developer

1643: .seealso: `PC`, `PetscViewer`, `KSPView()`, `PetscViewerASCIIOpen()`
1644: @*/
1645: PetscErrorCode PCView(PC pc, PetscViewer viewer)
1646: {
1647:   PCType    cstr;
1648:   PetscBool iascii, isstring, isbinary, isdraw;
1649: #if defined(PETSC_HAVE_SAWS)
1650:   PetscBool issaws;
1651: #endif

1654:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc), &viewer);

1658:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1659:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring);
1660:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary);
1661:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
1662: #if defined(PETSC_HAVE_SAWS)
1663:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws);
1664: #endif

1666:   if (iascii) {
1667:     PetscObjectPrintClassNamePrefixType((PetscObject)pc, viewer);
1668:     if (!pc->setupcalled) PetscViewerASCIIPrintf(viewer, "  PC has not been set up so information may be incomplete\n");
1669:     PetscViewerASCIIPushTab(viewer);
1670:     PetscTryTypeMethod(pc, view, viewer);
1671:     PetscViewerASCIIPopTab(viewer);
1672:     if (pc->mat) {
1673:       PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);
1674:       if (pc->pmat == pc->mat) {
1675:         PetscViewerASCIIPrintf(viewer, "  linear system matrix = precond matrix:\n");
1676:         PetscViewerASCIIPushTab(viewer);
1677:         MatView(pc->mat, viewer);
1678:         PetscViewerASCIIPopTab(viewer);
1679:       } else {
1680:         if (pc->pmat) {
1681:           PetscViewerASCIIPrintf(viewer, "  linear system matrix followed by preconditioner matrix:\n");
1682:         } else {
1683:           PetscViewerASCIIPrintf(viewer, "  linear system matrix:\n");
1684:         }
1685:         PetscViewerASCIIPushTab(viewer);
1686:         MatView(pc->mat, viewer);
1687:         if (pc->pmat) MatView(pc->pmat, viewer);
1688:         PetscViewerASCIIPopTab(viewer);
1689:       }
1690:       PetscViewerPopFormat(viewer);
1691:     }
1692:   } else if (isstring) {
1693:     PCGetType(pc, &cstr);
1694:     PetscViewerStringSPrintf(viewer, " PCType: %-7.7s", cstr);
1695:     PetscTryTypeMethod(pc, view, viewer);
1696:     if (pc->mat) MatView(pc->mat, viewer);
1697:     if (pc->pmat && pc->pmat != pc->mat) MatView(pc->pmat, viewer);
1698:   } else if (isbinary) {
1699:     PetscInt    classid = PC_FILE_CLASSID;
1700:     MPI_Comm    comm;
1701:     PetscMPIInt rank;
1702:     char        type[256];

1704:     PetscObjectGetComm((PetscObject)pc, &comm);
1705:     MPI_Comm_rank(comm, &rank);
1706:     if (rank == 0) {
1707:       PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT);
1708:       PetscStrncpy(type, ((PetscObject)pc)->type_name, 256);
1709:       PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR);
1710:     }
1711:     PetscTryTypeMethod(pc, view, viewer);
1712:   } else if (isdraw) {
1713:     PetscDraw draw;
1714:     char      str[25];
1715:     PetscReal x, y, bottom, h;
1716:     PetscInt  n;

1718:     PetscViewerDrawGetDraw(viewer, 0, &draw);
1719:     PetscDrawGetCurrentPoint(draw, &x, &y);
1720:     if (pc->mat) {
1721:       MatGetSize(pc->mat, &n, NULL);
1722:       PetscSNPrintf(str, 25, "PC: %s (%" PetscInt_FMT ")", ((PetscObject)pc)->type_name, n);
1723:     } else {
1724:       PetscSNPrintf(str, 25, "PC: %s", ((PetscObject)pc)->type_name);
1725:     }
1726:     PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h);
1727:     bottom = y - h;
1728:     PetscDrawPushCurrentPoint(draw, x, bottom);
1729:     PetscTryTypeMethod(pc, view, viewer);
1730:     PetscDrawPopCurrentPoint(draw);
1731: #if defined(PETSC_HAVE_SAWS)
1732:   } else if (issaws) {
1733:     PetscMPIInt rank;

1735:     PetscObjectName((PetscObject)pc);
1736:     MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
1737:     if (!((PetscObject)pc)->amsmem && rank == 0) PetscObjectViewSAWs((PetscObject)pc, viewer);
1738:     if (pc->mat) MatView(pc->mat, viewer);
1739:     if (pc->pmat && pc->pmat != pc->mat) MatView(pc->pmat, viewer);
1740: #endif
1741:   }
1742:   return 0;
1743: }

1745: /*@C
1746:   PCRegister -  Adds a method to the preconditioner package.

1748:    Not collective

1750:    Input Parameters:
1751: +  name_solver - name of a new user-defined solver
1752: -  routine_create - routine to create method context

1754:    Note:
1755:    `PCRegister()` may be called multiple times to add several user-defined preconditioners.

1757:    Sample usage:
1758: .vb
1759:    PCRegister("my_solver", MySolverCreate);
1760: .ve

1762:    Then, your solver can be chosen with the procedural interface via
1763: $     PCSetType(pc,"my_solver")
1764:    or at runtime via the option
1765: $     -pc_type my_solver

1767:    Level: advanced

1769: .seealso: `PC`, `PCType`, `PCRegisterAll()`
1770: @*/
1771: PetscErrorCode PCRegister(const char sname[], PetscErrorCode (*function)(PC))
1772: {
1773:   PCInitializePackage();
1774:   PetscFunctionListAdd(&PCList, sname, function);
1775:   return 0;
1776: }

1778: static PetscErrorCode MatMult_PC(Mat A, Vec X, Vec Y)
1779: {
1780:   PC pc;

1782:   MatShellGetContext(A, &pc);
1783:   PCApply(pc, X, Y);
1784:   return 0;
1785: }

1787: /*@
1788:     PCComputeOperator - Computes the explicit preconditioned operator.

1790:     Collective

1792:     Input Parameters:
1793: +   pc - the preconditioner object
1794: -   mattype - the matrix type to be used for the operator

1796:     Output Parameter:
1797: .   mat - the explicit preconditioned operator

1799:     Note:
1800:     This computation is done by applying the operators to columns of the identity matrix.
1801:     This routine is costly in general, and is recommended for use only with relatively small systems.
1802:     Currently, this routine uses a dense matrix format when mattype == NULL

1804:     Level: advanced

1806: .seealso: `PC`, `KSPComputeOperator()`, `MatType`

1808: @*/
1809: PetscErrorCode PCComputeOperator(PC pc, MatType mattype, Mat *mat)
1810: {
1811:   PetscInt N, M, m, n;
1812:   Mat      A, Apc;

1816:   PCGetOperators(pc, &A, NULL);
1817:   MatGetLocalSize(A, &m, &n);
1818:   MatGetSize(A, &M, &N);
1819:   MatCreateShell(PetscObjectComm((PetscObject)pc), m, n, M, N, pc, &Apc);
1820:   MatShellSetOperation(Apc, MATOP_MULT, (void (*)(void))MatMult_PC);
1821:   MatComputeOperator(Apc, mattype, mat);
1822:   MatDestroy(&Apc);
1823:   return 0;
1824: }

1826: /*@
1827:    PCSetCoordinates - sets the coordinates of all the nodes on the local process

1829:    Collective

1831:    Input Parameters:
1832: +  pc - the solver context
1833: .  dim - the dimension of the coordinates 1, 2, or 3
1834: .  nloc - the blocked size of the coordinates array
1835: -  coords - the coordinates array

1837:    Level: intermediate

1839:    Note:
1840:    coords is an array of the dim coordinates for the nodes on
1841:    the local processor, of size dim*nloc.
1842:    If there are 108 equation on a processor
1843:    for a displacement finite element discretization of elasticity (so
1844:    that there are nloc = 36 = 108/3 nodes) then the array must have 108
1845:    double precision values (ie, 3 * 36).  These x y z coordinates
1846:    should be ordered for nodes 0 to N-1 like so: [ 0.x, 0.y, 0.z, 1.x,
1847:    ... , N-1.z ].

1849: .seealso: `PC`, `MatSetNearNullSpace()`
1850: @*/
1851: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[])
1852: {
1855:   PetscTryMethod(pc, "PCSetCoordinates_C", (PC, PetscInt, PetscInt, PetscReal *), (pc, dim, nloc, coords));
1856:   return 0;
1857: }

1859: /*@
1860:    PCGetInterpolations - Gets interpolation matrices for all levels (except level 0)

1862:    Logically Collective

1864:    Input Parameter:
1865: .  pc - the precondition context

1867:    Output Parameters:
1868: +  num_levels - the number of levels
1869: -  interpolations - the interpolation matrices (size of num_levels-1)

1871:    Level: advanced

1873:    Developer Note:
1874:    Why is this here instead of in `PCMG` etc?

1876: .seealso: `PC`, `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetInterpolation()`, `PCGetCoarseOperators()`
1877: @*/
1878: PetscErrorCode PCGetInterpolations(PC pc, PetscInt *num_levels, Mat *interpolations[])
1879: {
1883:   PetscUseMethod(pc, "PCGetInterpolations_C", (PC, PetscInt *, Mat *[]), (pc, num_levels, interpolations));
1884:   return 0;
1885: }

1887: /*@
1888:    PCGetCoarseOperators - Gets coarse operator matrices for all levels (except the finest level)

1890:    Logically Collective

1892:    Input Parameter:
1893: .  pc - the precondition context

1895:    Output Parameters:
1896: +  num_levels - the number of levels
1897: -  coarseOperators - the coarse operator matrices (size of num_levels-1)

1899:    Level: advanced

1901:    Developer Note:
1902:    Why is this here instead of in `PCMG` etc?

1904: .seealso: `PC`, `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetRScale()`, `PCMGGetInterpolation()`, `PCGetInterpolations()`
1905: @*/
1906: PetscErrorCode PCGetCoarseOperators(PC pc, PetscInt *num_levels, Mat *coarseOperators[])
1907: {
1911:   PetscUseMethod(pc, "PCGetCoarseOperators_C", (PC, PetscInt *, Mat *[]), (pc, num_levels, coarseOperators));
1912:   return 0;
1913: }