Actual source code: precon.c

petsc-3.13.5 2020-09-01
Report Typos and Errors

  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_ApplyCoarse, PC_ApplyMultiple, PC_ApplySymmetricLeft;
 11: PetscLogEvent PC_ApplySymmetricRight, PC_ModifySubMatrices, PC_ApplyOnBlocks, PC_ApplyTransposeOnBlocks;
 12: PetscInt      PetscMGLevelId;

 14: PetscErrorCode PCGetDefaultType_Private(PC pc,const char *type[])
 15: {
 17:   PetscMPIInt    size;
 18:   PetscBool      hasop,flg1,flg2,set,flg3;

 21:   MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
 22:   if (pc->pmat) {
 23:     MatHasOperation(pc->pmat,MATOP_GET_DIAGONAL_BLOCK,&hasop);
 24:     if (size == 1) {
 25:       MatGetFactorAvailable(pc->pmat,"petsc",MAT_FACTOR_ICC,&flg1);
 26:       MatGetFactorAvailable(pc->pmat,"petsc",MAT_FACTOR_ILU,&flg2);
 27:       MatIsSymmetricKnown(pc->pmat,&set,&flg3);
 28:       if (flg1 && (!flg2 || (set && flg3))) {
 29:         *type = PCICC;
 30:       } else if (flg2) {
 31:         *type = PCILU;
 32:       } else if (hasop) { /* likely is a parallel matrix run on one processor */
 33:         *type = PCBJACOBI;
 34:       } else {
 35:         *type = PCNONE;
 36:       }
 37:     } else {
 38:        if (hasop) {
 39:         *type = PCBJACOBI;
 40:       } else {
 41:         *type = PCNONE;
 42:       }
 43:     }
 44:   } else {
 45:     if (size == 1) {
 46:       *type = PCILU;
 47:     } else {
 48:       *type = PCBJACOBI;
 49:     }
 50:   }
 51:   return(0);
 52: }

 54: /*@
 55:    PCReset - Resets a PC context to the pcsetupcalled = 0 state and removes any allocated Vecs and Mats

 57:    Collective on PC

 59:    Input Parameter:
 60: .  pc - the preconditioner context

 62:    Level: developer

 64:    Notes:
 65:     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

 67: .seealso: PCCreate(), PCSetUp()
 68: @*/
 69: PetscErrorCode  PCReset(PC pc)
 70: {

 75:   if (pc->ops->reset) {
 76:     (*pc->ops->reset)(pc);
 77:   }
 78:   VecDestroy(&pc->diagonalscaleright);
 79:   VecDestroy(&pc->diagonalscaleleft);
 80:   MatDestroy(&pc->pmat);
 81:   MatDestroy(&pc->mat);

 83:   pc->setupcalled = 0;
 84:   return(0);
 85: }

 87: /*@
 88:    PCDestroy - Destroys PC context that was created with PCCreate().

 90:    Collective on PC

 92:    Input Parameter:
 93: .  pc - the preconditioner context

 95:    Level: developer

 97: .seealso: PCCreate(), PCSetUp()
 98: @*/
 99: PetscErrorCode  PCDestroy(PC *pc)
100: {

104:   if (!*pc) return(0);
106:   if (--((PetscObject)(*pc))->refct > 0) {*pc = 0; return(0);}

108:   PCReset(*pc);

110:   /* if memory was published with SAWs then destroy it */
111:   PetscObjectSAWsViewOff((PetscObject)*pc);
112:   if ((*pc)->ops->destroy) {(*(*pc)->ops->destroy)((*pc));}
113:   DMDestroy(&(*pc)->dm);
114:   PetscHeaderDestroy(pc);
115:   return(0);
116: }

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

122:    Logically Collective on PC

124:    Input Parameter:
125: .  pc - the preconditioner context

127:    Output Parameter:
128: .  flag - PETSC_TRUE if it applies the scaling

130:    Level: developer

132:    Notes:
133:     If this returns PETSC_TRUE then the system solved via the Krylov method is
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

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

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

152:    Logically Collective on PC

154:    Input Parameters:
155: +  pc - the preconditioner context
156: -  s - scaling vector

158:    Level: intermediate

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

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

176:   pc->diagonalscale     = PETSC_TRUE;

178:   PetscObjectReference((PetscObject)s);
179:   VecDestroy(&pc->diagonalscaleleft);

181:   pc->diagonalscaleleft = s;

183:   VecDuplicate(s,&pc->diagonalscaleright);
184:   VecCopy(s,pc->diagonalscaleright);
185:   VecReciprocal(pc->diagonalscaleright);
186:   return(0);
187: }

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

192:    Logically Collective on PC

194:    Input Parameters:
195: +  pc - the preconditioner context
196: .  in - input vector
197: -  out - scaled vector (maybe the same as in)

199:    Level: intermediate

201:    Notes:
202:     The system solved via the Krylov method is
203: $           D M A D^{-1} y = D M b  for left preconditioning or
204: $           D A M D^{-1} z = D b for right preconditioning

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

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

210: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleSet(), PCDiagonalScaleRight(), PCDiagonalScale()
211: @*/
212: PetscErrorCode  PCDiagonalScaleLeft(PC pc,Vec in,Vec out)
213: {

220:   if (pc->diagonalscale) {
221:     VecPointwiseMult(out,pc->diagonalscaleleft,in);
222:   } else if (in != out) {
223:     VecCopy(in,out);
224:   }
225:   return(0);
226: }

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

231:    Logically Collective on PC

233:    Input Parameters:
234: +  pc - the preconditioner context
235: .  in - input vector
236: -  out - scaled vector (maybe the same as in)

238:    Level: intermediate

240:    Notes:
241:     The system solved via the Krylov method is
242: $           D M A D^{-1} y = D M b  for left preconditioning or
243: $           D A M D^{-1} z = D b for right preconditioning

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

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

249: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleSet(), PCDiagonalScale()
250: @*/
251: PetscErrorCode  PCDiagonalScaleRight(PC pc,Vec in,Vec out)
252: {

259:   if (pc->diagonalscale) {
260:     VecPointwiseMult(out,pc->diagonalscaleright,in);
261:   } else if (in != out) {
262:     VecCopy(in,out);
263:   }
264:   return(0);
265: }

267: /*@
268:    PCSetUseAmat - Sets a flag to indicate that when the preconditioner needs to apply (part of) the
269:    operator during the preconditioning process it applies the Amat provided to TSSetRHSJacobian(),
270:    TSSetIJacobian(), SNESSetJacobian(), KSPSetOperator() or PCSetOperator() not the Pmat.

272:    Logically Collective on PC

274:    Input Parameters:
275: +  pc - the preconditioner context
276: -  flg - PETSC_TRUE to use the Amat, PETSC_FALSE to use the Pmat (default is false)

278:    Options Database Key:
279: .  -pc_use_amat <true,false>

281:    Notes:
282:    For the common case in which the linear system matrix and the matrix used to construct the
283:    preconditioner are identical, this routine is does nothing.

285:    Level: intermediate

287: .seealso: PCGetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
288: @*/
289: PetscErrorCode  PCSetUseAmat(PC pc,PetscBool flg)
290: {
293:   pc->useAmat = flg;
294:   return(0);
295: }

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

300:    Logically Collective on PC

302:    Input Parameters:
303: +  pc - iterative context obtained from PCCreate()
304: -  flg - PETSC_TRUE indicates you want the error generated

306:    Level: advanced

308:    Notes:
309:     Normally PETSc continues if a linear solver fails due to a failed setup of a preconditioner, you can call KSPGetConvergedReason() after a KSPSolve()
310:     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
311:     detected.

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

315: .seealso: PCGetInitialGuessNonzero(), PCSetInitialGuessKnoll(), PCGetInitialGuessKnoll()
316: @*/
317: PetscErrorCode  PCSetErrorIfFailure(PC pc,PetscBool flg)
318: {
322:   pc->erroriffailure = flg;
323:   return(0);
324: }

326: /*@
327:    PCGetUseAmat - Gets a flag to indicate that when the preconditioner needs to apply (part of) the
328:    operator during the preconditioning process it applies the Amat provided to TSSetRHSJacobian(),
329:    TSSetIJacobian(), SNESSetJacobian(), KSPSetOperator() or PCSetOperator() not the Pmat.

331:    Logically Collective on PC

333:    Input Parameter:
334: .  pc - the preconditioner context

336:    Output Parameter:
337: .  flg - PETSC_TRUE to use the Amat, PETSC_FALSE to use the Pmat (default is false)

339:    Notes:
340:    For the common case in which the linear system matrix and the matrix used to construct the
341:    preconditioner are identical, this routine is does nothing.

343:    Level: intermediate

345: .seealso: PCSetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
346: @*/
347: PetscErrorCode  PCGetUseAmat(PC pc,PetscBool *flg)
348: {
351:   *flg = pc->useAmat;
352:   return(0);
353: }

355: /*@
356:    PCCreate - Creates a preconditioner context.

358:    Collective

360:    Input Parameter:
361: .  comm - MPI communicator

363:    Output Parameter:
364: .  pc - location to put the preconditioner context

366:    Notes:
367:    The default preconditioner for sparse matrices is PCILU or PCICC with 0 fill on one process and block Jacobi with PCILU or PCICC
368:    in parallel. For dense matrices it is always PCNONE.

370:    Level: developer

372: .seealso: PCSetUp(), PCApply(), PCDestroy()
373: @*/
374: PetscErrorCode  PCCreate(MPI_Comm comm,PC *newpc)
375: {
376:   PC             pc;

381:   *newpc = 0;
382:   PCInitializePackage();

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

386:   pc->mat                  = 0;
387:   pc->pmat                 = 0;
388:   pc->setupcalled          = 0;
389:   pc->setfromoptionscalled = 0;
390:   pc->data                 = 0;
391:   pc->diagonalscale        = PETSC_FALSE;
392:   pc->diagonalscaleleft    = 0;
393:   pc->diagonalscaleright   = 0;

395:   pc->modifysubmatrices  = 0;
396:   pc->modifysubmatricesP = 0;

398:   *newpc = pc;
399:   return(0);

401: }

403: /* -------------------------------------------------------------------------------*/

405: /*@
406:    PCApply - Applies the preconditioner to a vector.

408:    Collective on PC

410:    Input Parameters:
411: +  pc - the preconditioner context
412: -  x - input vector

414:    Output Parameter:
415: .  y - output vector

417:    Level: developer

419: .seealso: PCApplyTranspose(), PCApplyBAorAB()
420: @*/
421: PetscErrorCode  PCApply(PC pc,Vec x,Vec y)
422: {
424:   PetscInt       m,n,mv,nv;

430:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
431:   if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
432:   /* use pmat to check vector sizes since for KSPLQR the pmat may be of a different size than mat */
433:   MatGetLocalSize(pc->pmat,&m,&n);
434:   VecGetLocalSize(x,&nv);
435:   VecGetLocalSize(y,&mv);
436:   if (mv != m) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Preconditioner number of local rows %D does not equal resulting vector number of rows %D",m,mv);
437:   if (nv != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Preconditioner number of local columns %D does not equal resulting vector number of rows %D",n,nv);
438:   VecSetErrorIfLocked(y,3);

440:   PCSetUp(pc);
441:   if (!pc->ops->apply) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply");
442:   VecLockReadPush(x);
443:   PetscLogEventBegin(PC_Apply,pc,x,y,0);
444:   (*pc->ops->apply)(pc,x,y);
445:   PetscLogEventEnd(PC_Apply,pc,x,y,0);
446:   if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
447:   VecLockReadPop(x);
448:   return(0);
449: }

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

454:    Collective on PC

456:    Input Parameters:
457: +  pc - the preconditioner context
458: -  x - input vector

460:    Output Parameter:
461: .  y - output vector

463:    Notes:
464:    Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.

466:    Level: developer

468: .seealso: PCApply(), PCApplySymmetricRight()
469: @*/
470: PetscErrorCode  PCApplySymmetricLeft(PC pc,Vec x,Vec y)
471: {

478:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
479:   if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
480:   PCSetUp(pc);
481:   if (!pc->ops->applysymmetricleft) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have left symmetric apply");
482:   VecLockReadPush(x);
483:   PetscLogEventBegin(PC_ApplySymmetricLeft,pc,x,y,0);
484:   (*pc->ops->applysymmetricleft)(pc,x,y);
485:   PetscLogEventEnd(PC_ApplySymmetricLeft,pc,x,y,0);
486:   VecLockReadPop(x);
487:   if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
488:   return(0);
489: }

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

494:    Collective on PC

496:    Input Parameters:
497: +  pc - the preconditioner context
498: -  x - input vector

500:    Output Parameter:
501: .  y - output vector

503:    Level: developer

505:    Notes:
506:    Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.

508: .seealso: PCApply(), PCApplySymmetricLeft()
509: @*/
510: PetscErrorCode  PCApplySymmetricRight(PC pc,Vec x,Vec y)
511: {

518:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
519:   if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
520:   PCSetUp(pc);
521:   if (!pc->ops->applysymmetricright) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have left symmetric apply");
522:   VecLockReadPush(x);
523:   PetscLogEventBegin(PC_ApplySymmetricRight,pc,x,y,0);
524:   (*pc->ops->applysymmetricright)(pc,x,y);
525:   PetscLogEventEnd(PC_ApplySymmetricRight,pc,x,y,0);
526:   VecLockReadPop(x);
527:   if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
528:   return(0);
529: }

531: /*@
532:    PCApplyTranspose - Applies the transpose of preconditioner to a vector.

534:    Collective on PC

536:    Input Parameters:
537: +  pc - the preconditioner context
538: -  x - input vector

540:    Output Parameter:
541: .  y - output vector

543:    Notes:
544:     For complex numbers this applies the non-Hermitian transpose.

546:    Developer Notes:
547:     We need to implement a PCApplyHermitianTranspose()

549:    Level: developer

551: .seealso: PCApply(), PCApplyBAorAB(), PCApplyBAorABTranspose(), PCApplyTransposeExists()
552: @*/
553: PetscErrorCode  PCApplyTranspose(PC pc,Vec x,Vec y)
554: {

561:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
562:   if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
563:   PCSetUp(pc);
564:   if (!pc->ops->applytranspose) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply transpose");
565:   VecLockReadPush(x);
566:   PetscLogEventBegin(PC_Apply,pc,x,y,0);
567:   (*pc->ops->applytranspose)(pc,x,y);
568:   PetscLogEventEnd(PC_Apply,pc,x,y,0);
569:   VecLockReadPop(x);
570:   if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
571:   return(0);
572: }

574: /*@
575:    PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation

577:    Collective on PC

579:    Input Parameters:
580: .  pc - the preconditioner context

582:    Output Parameter:
583: .  flg - PETSC_TRUE if a transpose operation is defined

585:    Level: developer

587: .seealso: PCApplyTranspose()
588: @*/
589: PetscErrorCode  PCApplyTransposeExists(PC pc,PetscBool  *flg)
590: {
594:   if (pc->ops->applytranspose) *flg = PETSC_TRUE;
595:   else *flg = PETSC_FALSE;
596:   return(0);
597: }

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

602:    Collective on PC

604:    Input Parameters:
605: +  pc - the preconditioner context
606: .  side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
607: .  x - input vector
608: -  work - work vector

610:    Output Parameter:
611: .  y - output vector

613:    Level: developer

615:    Notes:
616:     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
617:    specific KSPSolve() method must also be written to handle the post-solve "correction" for the diagonal scaling.

619: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorABTranspose()
620: @*/
621: PetscErrorCode  PCApplyBAorAB(PC pc,PCSide side,Vec x,Vec y,Vec work)
622: {

634:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
635:   if (side != PC_LEFT && side != PC_SYMMETRIC && side != PC_RIGHT) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Side must be right, left, or symmetric");
636:   if (pc->diagonalscale && side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot include diagonal scaling with symmetric preconditioner Section 1.5 Writing Application Codes with PETSc");
637:   if (pc->erroriffailure) {VecValidValues(x,3,PETSC_TRUE);}

639:   PCSetUp(pc);
640:   if (pc->diagonalscale) {
641:     if (pc->ops->applyBA) {
642:       Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
643:       VecDuplicate(x,&work2);
644:       PCDiagonalScaleRight(pc,x,work2);
645:       (*pc->ops->applyBA)(pc,side,work2,y,work);
646:       PCDiagonalScaleLeft(pc,y,y);
647:       VecDestroy(&work2);
648:     } else if (side == PC_RIGHT) {
649:       PCDiagonalScaleRight(pc,x,y);
650:       PCApply(pc,y,work);
651:       MatMult(pc->mat,work,y);
652:       PCDiagonalScaleLeft(pc,y,y);
653:     } else if (side == PC_LEFT) {
654:       PCDiagonalScaleRight(pc,x,y);
655:       MatMult(pc->mat,y,work);
656:       PCApply(pc,work,y);
657:       PCDiagonalScaleLeft(pc,y,y);
658:     } else if (side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot provide diagonal scaling with symmetric Section 1.5 Writing Application Codes with PETSc of preconditioner");
659:   } else {
660:     if (pc->ops->applyBA) {
661:       (*pc->ops->applyBA)(pc,side,x,y,work);
662:     } else if (side == PC_RIGHT) {
663:       PCApply(pc,x,work);
664:       MatMult(pc->mat,work,y);
665:     } else if (side == PC_LEFT) {
666:       MatMult(pc->mat,x,work);
667:       PCApply(pc,work,y);
668:     } else if (side == PC_SYMMETRIC) {
669:       /* There's an extra copy here; maybe should provide 2 work vectors instead? */
670:       PCApplySymmetricRight(pc,x,work);
671:       MatMult(pc->mat,work,y);
672:       VecCopy(y,work);
673:       PCApplySymmetricLeft(pc,work,y);
674:     }
675:   }
676:   if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
677:   return(0);
678: }

680: /*@
681:    PCApplyBAorABTranspose - Applies the transpose of the preconditioner
682:    and operator to a vector. That is, applies tr(B) * tr(A) with left preconditioning,
683:    NOT tr(B*A) = tr(A)*tr(B).

685:    Collective on PC

687:    Input Parameters:
688: +  pc - the preconditioner context
689: .  side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
690: .  x - input vector
691: -  work - work vector

693:    Output Parameter:
694: .  y - output vector


697:    Notes:
698:     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
699:       defined by B'. This is why this has the funny form that it computes tr(B) * tr(A)

701:     Level: developer

703: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorAB()
704: @*/
705: PetscErrorCode  PCApplyBAorABTranspose(PC pc,PCSide side,Vec x,Vec y,Vec work)
706: {

714:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
715:   if (pc->erroriffailure) {VecValidValues(x,3,PETSC_TRUE);}
716:   if (pc->ops->applyBAtranspose) {
717:     (*pc->ops->applyBAtranspose)(pc,side,x,y,work);
718:     if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
719:     return(0);
720:   }
721:   if (side != PC_LEFT && side != PC_RIGHT) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Side must be right or left");

723:   PCSetUp(pc);
724:   if (side == PC_RIGHT) {
725:     PCApplyTranspose(pc,x,work);
726:     MatMultTranspose(pc->mat,work,y);
727:   } else if (side == PC_LEFT) {
728:     MatMultTranspose(pc->mat,x,work);
729:     PCApplyTranspose(pc,work,y);
730:   }
731:   /* add support for PC_SYMMETRIC */
732:   if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
733:   return(0);
734: }

736: /* -------------------------------------------------------------------------------*/

738: /*@
739:    PCApplyRichardsonExists - Determines whether a particular preconditioner has a
740:    built-in fast Section 1.5 Writing Application Codes with PETSc of Richardson's method.

742:    Not Collective

744:    Input Parameter:
745: .  pc - the preconditioner

747:    Output Parameter:
748: .  exists - PETSC_TRUE or PETSC_FALSE

750:    Level: developer

752: .seealso: PCApplyRichardson()
753: @*/
754: PetscErrorCode  PCApplyRichardsonExists(PC pc,PetscBool  *exists)
755: {
759:   if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
760:   else *exists = PETSC_FALSE;
761:   return(0);
762: }

764: /*@
765:    PCApplyRichardson - Applies several steps of Richardson iteration with
766:    the particular preconditioner. This routine is usually used by the
767:    Krylov solvers and not the Section 1.5 Writing Application Codes with PETSc code directly.

769:    Collective on PC

771:    Input Parameters:
772: +  pc  - the preconditioner context
773: .  b   - the right hand side
774: .  w   - one work vector
775: .  rtol - relative decrease in residual norm convergence criteria
776: .  abstol - absolute residual norm convergence criteria
777: .  dtol - divergence residual norm increase criteria
778: .  its - the number of iterations to apply.
779: -  guesszero - if the input x contains nonzero initial guess

781:    Output Parameter:
782: +  outits - number of iterations actually used (for SOR this always equals its)
783: .  reason - the reason the apply terminated
784: -  y - the solution (also contains initial guess if guesszero is PETSC_FALSE

786:    Notes:
787:    Most preconditioners do not support this function. Use the command
788:    PCApplyRichardsonExists() to determine if one does.

790:    Except for the multigrid PC this routine ignores the convergence tolerances
791:    and always runs for the number of iterations

793:    Level: developer

795: .seealso: PCApplyRichardsonExists()
796: @*/
797: PetscErrorCode  PCApplyRichardson(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
798: {

806:   if (b == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"b and y must be different vectors");
807:   PCSetUp(pc);
808:   if (!pc->ops->applyrichardson) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply richardson");
809:   (*pc->ops->applyrichardson)(pc,b,y,w,rtol,abstol,dtol,its,guesszero,outits,reason);
810:   return(0);
811: }

813: /*@
814:    PCGetFailedReason - Gets the reason a PCSetUp() failed or 0 if it did not fail

816:    Logically Collective on PC

818:    Input Parameter:
819: .  pc - the preconditioner context

821:    Output Parameter:
822: .  reason - the reason it failed, currently only -1

824:    Level: advanced

826: .seealso: PCCreate(), PCApply(), PCDestroy()
827: @*/
828: PetscErrorCode PCGetFailedReason(PC pc,PCFailedReason *reason)
829: {
831:   if (pc->setupcalled < 0) *reason = (PCFailedReason)pc->setupcalled;
832:   else *reason = pc->failedreason;
833:   return(0);
834: }


837: /*
838:       a setupcall of 0 indicates never setup,
839:                      1 indicates has been previously setup
840:                     -1 indicates a PCSetUp() was attempted and failed
841: */
842: /*@
843:    PCSetUp - Prepares for the use of a preconditioner.

845:    Collective on PC

847:    Input Parameter:
848: .  pc - the preconditioner context

850:    Level: developer

852: .seealso: PCCreate(), PCApply(), PCDestroy()
853: @*/
854: PetscErrorCode  PCSetUp(PC pc)
855: {
856:   PetscErrorCode   ierr;
857:   const char       *def;
858:   PetscObjectState matstate, matnonzerostate;

862:   if (!pc->mat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Matrix must be set first");

864:   if (pc->setupcalled && pc->reusepreconditioner) {
865:     PetscInfo(pc,"Leaving PC with identical preconditioner since reuse preconditioner is set\n");
866:     return(0);
867:   }

869:   PetscObjectStateGet((PetscObject)pc->pmat,&matstate);
870:   MatGetNonzeroState(pc->pmat,&matnonzerostate);
871:   if (!pc->setupcalled) {
872:     PetscInfo(pc,"Setting up PC for first time\n");
873:     pc->flag = DIFFERENT_NONZERO_PATTERN;
874:   } else if (matstate == pc->matstate) {
875:     PetscInfo(pc,"Leaving PC with identical preconditioner since operator is unchanged\n");
876:     return(0);
877:   } else {
878:     if (matnonzerostate > pc->matnonzerostate) {
879:        PetscInfo(pc,"Setting up PC with different nonzero pattern\n");
880:        pc->flag = DIFFERENT_NONZERO_PATTERN;
881:     } else {
882:       PetscInfo(pc,"Setting up PC with same nonzero pattern\n");
883:       pc->flag = SAME_NONZERO_PATTERN;
884:     }
885:   }
886:   pc->matstate        = matstate;
887:   pc->matnonzerostate = matnonzerostate;

889:   if (!((PetscObject)pc)->type_name) {
890:     PCGetDefaultType_Private(pc,&def);
891:     PCSetType(pc,def);
892:   }

894:   MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
895:   MatSetErrorIfFailure(pc->mat,pc->erroriffailure);
896:   PetscLogEventBegin(PC_SetUp,pc,0,0,0);
897:   if (pc->ops->setup) {
898:     (*pc->ops->setup)(pc);
899:   }
900:   PetscLogEventEnd(PC_SetUp,pc,0,0,0);
901:   if (!pc->setupcalled) pc->setupcalled = 1;
902:   return(0);
903: }

905: /*@
906:    PCSetUpOnBlocks - Sets up the preconditioner for each block in
907:    the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
908:    methods.

910:    Collective on PC

912:    Input Parameters:
913: .  pc - the preconditioner context

915:    Level: developer

917: .seealso: PCCreate(), PCApply(), PCDestroy(), PCSetUp()
918: @*/
919: PetscErrorCode  PCSetUpOnBlocks(PC pc)
920: {

925:   if (!pc->ops->setuponblocks) return(0);
926:   PetscLogEventBegin(PC_SetUpOnBlocks,pc,0,0,0);
927:   (*pc->ops->setuponblocks)(pc);
928:   PetscLogEventEnd(PC_SetUpOnBlocks,pc,0,0,0);
929:   return(0);
930: }

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

939:    Logically Collective on PC

941:    Input Parameters:
942: +  pc - the preconditioner context
943: .  func - routine for modifying the submatrices
944: -  ctx - optional user-defined context (may be null)

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

949: +  row - an array of index sets that contain the global row numbers
950:          that comprise each local submatrix
951: .  col - an array of index sets that contain the global column numbers
952:          that comprise each local submatrix
953: .  submat - array of local submatrices
954: -  ctx - optional user-defined context for private data for the
955:          user-defined func routine (may be null)

957:    Notes:
958:    PCSetModifySubMatrices() MUST be called before KSPSetUp() and
959:    KSPSolve().

961:    A routine set by PCSetModifySubMatrices() is currently called within
962:    the block Jacobi (PCBJACOBI) and additive Schwarz (PCASM)
963:    preconditioners.  All other preconditioners ignore this routine.

965:    Level: advanced

967: .seealso: PCModifySubMatrices()
968: @*/
969: PetscErrorCode  PCSetModifySubMatrices(PC pc,PetscErrorCode (*func)(PC,PetscInt,const IS[],const IS[],Mat[],void*),void *ctx)
970: {
973:   pc->modifysubmatrices  = func;
974:   pc->modifysubmatricesP = ctx;
975:   return(0);
976: }

978: /*@C
979:    PCModifySubMatrices - Calls an optional user-defined routine within
980:    certain preconditioners if one has been set with PCSetModifySubMatrices().

982:    Collective on PC

984:    Input Parameters:
985: +  pc - the preconditioner context
986: .  nsub - the number of local submatrices
987: .  row - an array of index sets that contain the global row numbers
988:          that comprise each local submatrix
989: .  col - an array of index sets that contain the global column numbers
990:          that comprise each local submatrix
991: .  submat - array of local submatrices
992: -  ctx - optional user-defined context for private data for the
993:          user-defined routine (may be null)

995:    Output Parameter:
996: .  submat - array of local submatrices (the entries of which may
997:             have been modified)

999:    Notes:
1000:    The user should NOT generally call this routine, as it will
1001:    automatically be called within certain preconditioners (currently
1002:    block Jacobi, additive Schwarz) if set.

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

1009:    Level: developer

1011: .seealso: PCSetModifySubMatrices()
1012: @*/
1013: PetscErrorCode  PCModifySubMatrices(PC pc,PetscInt nsub,const IS row[],const IS col[],Mat submat[],void *ctx)
1014: {

1019:   if (!pc->modifysubmatrices) return(0);
1020:   PetscLogEventBegin(PC_ModifySubMatrices,pc,0,0,0);
1021:   (*pc->modifysubmatrices)(pc,nsub,row,col,submat,ctx);
1022:   PetscLogEventEnd(PC_ModifySubMatrices,pc,0,0,0);
1023:   return(0);
1024: }

1026: /*@
1027:    PCSetOperators - Sets the matrix associated with the linear system and
1028:    a (possibly) different one associated with the preconditioner.

1030:    Logically Collective on PC

1032:    Input Parameters:
1033: +  pc - the preconditioner context
1034: .  Amat - the matrix that defines the linear system
1035: -  Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.

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

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

1044:    More Notes about Repeated Solution of Linear Systems:
1045:    PETSc does NOT reset the matrix entries of either Amat or Pmat
1046:    to zero after a linear solve; the user is completely responsible for
1047:    matrix assembly.  See the routine MatZeroEntries() if desiring to
1048:    zero all elements of a matrix.

1050:    Level: intermediate

1052: .seealso: PCGetOperators(), MatZeroEntries()
1053:  @*/
1054: PetscErrorCode  PCSetOperators(PC pc,Mat Amat,Mat Pmat)
1055: {
1056:   PetscErrorCode   ierr;
1057:   PetscInt         m1,n1,m2,n2;

1065:   if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
1066:     MatGetLocalSize(Amat,&m1,&n1);
1067:     MatGetLocalSize(pc->mat,&m2,&n2);
1068:     if (m1 != m2 || n1 != n2) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot change local size of Amat after use old sizes %D %D new sizes %D %D",m2,n2,m1,n1);
1069:     MatGetLocalSize(Pmat,&m1,&n1);
1070:     MatGetLocalSize(pc->pmat,&m2,&n2);
1071:     if (m1 != m2 || n1 != n2) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot change local size of Pmat after use old sizes %D %D new sizes %D %D",m2,n2,m1,n1);
1072:   }

1074:   if (Pmat != pc->pmat) {
1075:     /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
1076:     pc->matnonzerostate = -1;
1077:     pc->matstate        = -1;
1078:   }

1080:   /* reference first in case the matrices are the same */
1081:   if (Amat) {PetscObjectReference((PetscObject)Amat);}
1082:   MatDestroy(&pc->mat);
1083:   if (Pmat) {PetscObjectReference((PetscObject)Pmat);}
1084:   MatDestroy(&pc->pmat);
1085:   pc->mat  = Amat;
1086:   pc->pmat = Pmat;
1087:   return(0);
1088: }

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

1093:    Logically Collective on PC

1095:    Input Parameters:
1096: +  pc - the preconditioner context
1097: -  flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner

1099:     Level: intermediate

1101: .seealso: PCGetOperators(), MatZeroEntries(), PCGetReusePreconditioner(), KSPSetReusePreconditioner()
1102:  @*/
1103: PetscErrorCode  PCSetReusePreconditioner(PC pc,PetscBool flag)
1104: {
1107:   pc->reusepreconditioner = flag;
1108:   return(0);
1109: }

1111: /*@
1112:    PCGetReusePreconditioner - Determines if the PC reuses the current preconditioner even if the operator in the preconditioner has changed.

1114:    Not Collective

1116:    Input Parameter:
1117: .  pc - the preconditioner context

1119:    Output Parameter:
1120: .  flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner

1122:    Level: intermediate

1124: .seealso: PCGetOperators(), MatZeroEntries(), PCSetReusePreconditioner()
1125:  @*/
1126: PetscErrorCode  PCGetReusePreconditioner(PC pc,PetscBool *flag)
1127: {
1130:   *flag = pc->reusepreconditioner;
1131:   return(0);
1132: }

1134: /*@
1135:    PCGetOperators - Gets the matrix associated with the linear system and
1136:    possibly a different one associated with the preconditioner.

1138:    Not collective, though parallel Mats are returned if the PC is parallel

1140:    Input Parameter:
1141: .  pc - the preconditioner context

1143:    Output Parameters:
1144: +  Amat - the matrix defining the linear system
1145: -  Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.

1147:    Level: intermediate

1149:    Notes:
1150:     Does not increase the reference count of the matrices, so you should not destroy them

1152:    Alternative usage: If the operators have NOT been set with KSP/PCSetOperators() then the operators
1153:       are created in PC and returned to the user. In this case, if both operators
1154:       mat and pmat are requested, two DIFFERENT operators will be returned. If
1155:       only one is requested both operators in the PC will be the same (i.e. as
1156:       if one had called KSP/PCSetOperators() with the same argument for both Mats).
1157:       The user must set the sizes of the returned matrices and their type etc just
1158:       as if the user created them with MatCreate(). For example,

1160: $         KSP/PCGetOperators(ksp/pc,&Amat,NULL); is equivalent to
1161: $           set size, type, etc of Amat

1163: $         MatCreate(comm,&mat);
1164: $         KSP/PCSetOperators(ksp/pc,Amat,Amat);
1165: $         PetscObjectDereference((PetscObject)mat);
1166: $           set size, type, etc of Amat

1168:      and

1170: $         KSP/PCGetOperators(ksp/pc,&Amat,&Pmat); is equivalent to
1171: $           set size, type, etc of Amat and Pmat

1173: $         MatCreate(comm,&Amat);
1174: $         MatCreate(comm,&Pmat);
1175: $         KSP/PCSetOperators(ksp/pc,Amat,Pmat);
1176: $         PetscObjectDereference((PetscObject)Amat);
1177: $         PetscObjectDereference((PetscObject)Pmat);
1178: $           set size, type, etc of Amat and Pmat

1180:     The rational for this support is so that when creating a TS, SNES, or KSP the hierarchy
1181:     of underlying objects (i.e. SNES, KSP, PC, Mat) and their livespans can be completely
1182:     managed by the top most level object (i.e. the TS, SNES, or KSP). Another way to look
1183:     at this is when you create a SNES you do not NEED to create a KSP and attach it to
1184:     the SNES object (the SNES object manages it for you). Similarly when you create a KSP
1185:     you do not need to attach a PC to it (the KSP object manages the PC object for you).
1186:     Thus, why should YOU have to create the Mat and attach it to the SNES/KSP/PC, when
1187:     it can be created for you?


1190: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperatorsSet()
1191: @*/
1192: PetscErrorCode  PCGetOperators(PC pc,Mat *Amat,Mat *Pmat)
1193: {

1198:   if (Amat) {
1199:     if (!pc->mat) {
1200:       if (pc->pmat && !Pmat) {  /* Apmat has been set, but user did not request it, so use for Amat */
1201:         pc->mat = pc->pmat;
1202:         PetscObjectReference((PetscObject)pc->mat);
1203:       } else {                  /* both Amat and Pmat are empty */
1204:         MatCreate(PetscObjectComm((PetscObject)pc),&pc->mat);
1205:         if (!Pmat) { /* user did NOT request Pmat, so make same as Amat */
1206:           pc->pmat = pc->mat;
1207:           PetscObjectReference((PetscObject)pc->pmat);
1208:         }
1209:       }
1210:     }
1211:     *Amat = pc->mat;
1212:   }
1213:   if (Pmat) {
1214:     if (!pc->pmat) {
1215:       if (pc->mat && !Amat) {    /* Amat has been set but was not requested, so use for pmat */
1216:         pc->pmat = pc->mat;
1217:         PetscObjectReference((PetscObject)pc->pmat);
1218:       } else {
1219:         MatCreate(PetscObjectComm((PetscObject)pc),&pc->pmat);
1220:         if (!Amat) { /* user did NOT request Amat, so make same as Pmat */
1221:           pc->mat = pc->pmat;
1222:           PetscObjectReference((PetscObject)pc->mat);
1223:         }
1224:       }
1225:     }
1226:     *Pmat = pc->pmat;
1227:   }
1228:   return(0);
1229: }

1231: /*@C
1232:    PCGetOperatorsSet - Determines if the matrix associated with the linear system and
1233:    possibly a different one associated with the preconditioner have been set in the PC.

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

1237:    Input Parameter:
1238: .  pc - the preconditioner context

1240:    Output Parameters:
1241: +  mat - the matrix associated with the linear system was set
1242: -  pmat - matrix associated with the preconditioner was set, usually the same

1244:    Level: intermediate

1246: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperators()
1247: @*/
1248: PetscErrorCode  PCGetOperatorsSet(PC pc,PetscBool  *mat,PetscBool  *pmat)
1249: {
1252:   if (mat) *mat = (pc->mat)  ? PETSC_TRUE : PETSC_FALSE;
1253:   if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1254:   return(0);
1255: }

1257: /*@
1258:    PCFactorGetMatrix - Gets the factored matrix from the
1259:    preconditioner context.  This routine is valid only for the LU,
1260:    incomplete LU, Cholesky, and incomplete Cholesky methods.

1262:    Not Collective on PC though Mat is parallel if PC is parallel

1264:    Input Parameters:
1265: .  pc - the preconditioner context

1267:    Output parameters:
1268: .  mat - the factored matrix

1270:    Level: advanced

1272:    Notes:
1273:     Does not increase the reference count for the matrix so DO NOT destroy it

1275: @*/
1276: PetscErrorCode  PCFactorGetMatrix(PC pc,Mat *mat)
1277: {

1283:   if (pc->ops->getfactoredmatrix) {
1284:     (*pc->ops->getfactoredmatrix)(pc,mat);
1285:   } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC type does not support getting factor matrix");
1286:   return(0);
1287: }

1289: /*@C
1290:    PCSetOptionsPrefix - Sets the prefix used for searching for all
1291:    PC options in the database.

1293:    Logically Collective on PC

1295:    Input Parameters:
1296: +  pc - the preconditioner context
1297: -  prefix - the prefix string to prepend to all PC option requests

1299:    Notes:
1300:    A hyphen (-) must NOT be given at the beginning of the prefix name.
1301:    The first character of all runtime options is AUTOMATICALLY the
1302:    hyphen.

1304:    Level: advanced

1306: .seealso: PCAppendOptionsPrefix(), PCGetOptionsPrefix()
1307: @*/
1308: PetscErrorCode  PCSetOptionsPrefix(PC pc,const char prefix[])
1309: {

1314:   PetscObjectSetOptionsPrefix((PetscObject)pc,prefix);
1315:   return(0);
1316: }

1318: /*@C
1319:    PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1320:    PC options in the database.

1322:    Logically Collective on PC

1324:    Input Parameters:
1325: +  pc - the preconditioner context
1326: -  prefix - the prefix string to prepend to all PC option requests

1328:    Notes:
1329:    A hyphen (-) must NOT be given at the beginning of the prefix name.
1330:    The first character of all runtime options is AUTOMATICALLY the
1331:    hyphen.

1333:    Level: advanced

1335: .seealso: PCSetOptionsPrefix(), PCGetOptionsPrefix()
1336: @*/
1337: PetscErrorCode  PCAppendOptionsPrefix(PC pc,const char prefix[])
1338: {

1343:   PetscObjectAppendOptionsPrefix((PetscObject)pc,prefix);
1344:   return(0);
1345: }

1347: /*@C
1348:    PCGetOptionsPrefix - Gets the prefix used for searching for all
1349:    PC options in the database.

1351:    Not Collective

1353:    Input Parameters:
1354: .  pc - the preconditioner context

1356:    Output Parameters:
1357: .  prefix - pointer to the prefix string used, is returned

1359:    Notes:
1360:     On the fortran side, the user should pass in a string 'prifix' of
1361:    sufficient length to hold the prefix.

1363:    Level: advanced

1365: .seealso: PCSetOptionsPrefix(), PCAppendOptionsPrefix()
1366: @*/
1367: PetscErrorCode  PCGetOptionsPrefix(PC pc,const char *prefix[])
1368: {

1374:   PetscObjectGetOptionsPrefix((PetscObject)pc,prefix);
1375:   return(0);
1376: }

1378: /*
1379:    Indicates the right hand side will be changed by KSPSolve(), this occurs for a few
1380:   preconditioners including BDDC and Eisentat that transform the equations before applying
1381:   the Krylov methods
1382: */
1383: PETSC_INTERN PetscErrorCode  PCPreSolveChangeRHS(PC pc,PetscBool *change)
1384: {

1390:   *change = PETSC_FALSE;
1391:   PetscTryMethod(pc,"PCPreSolveChangeRHS_C",(PC,PetscBool*),(pc,change));
1392:   return(0);
1393: }

1395: /*@
1396:    PCPreSolve - Optional pre-solve phase, intended for any
1397:    preconditioner-specific actions that must be performed before
1398:    the iterative solve itself.

1400:    Collective on PC

1402:    Input Parameters:
1403: +  pc - the preconditioner context
1404: -  ksp - the Krylov subspace context

1406:    Level: developer

1408:    Sample of Usage:
1409: .vb
1410:     PCPreSolve(pc,ksp);
1411:     KSPSolve(ksp,b,x);
1412:     PCPostSolve(pc,ksp);
1413: .ve

1415:    Notes:
1416:    The pre-solve phase is distinct from the PCSetUp() phase.

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

1420: .seealso: PCPostSolve()
1421: @*/
1422: PetscErrorCode  PCPreSolve(PC pc,KSP ksp)
1423: {
1425:   Vec            x,rhs;

1430:   pc->presolvedone++;
1431:   if (pc->presolvedone > 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot embed PCPreSolve() more than twice");
1432:   KSPGetSolution(ksp,&x);
1433:   KSPGetRhs(ksp,&rhs);

1435:   if (pc->ops->presolve) {
1436:     (*pc->ops->presolve)(pc,ksp,rhs,x);
1437:   }
1438:   return(0);
1439: }

1441: /*@
1442:    PCPostSolve - Optional post-solve phase, intended for any
1443:    preconditioner-specific actions that must be performed after
1444:    the iterative solve itself.

1446:    Collective on PC

1448:    Input Parameters:
1449: +  pc - the preconditioner context
1450: -  ksp - the Krylov subspace context

1452:    Sample of Usage:
1453: .vb
1454:     PCPreSolve(pc,ksp);
1455:     KSPSolve(ksp,b,x);
1456:     PCPostSolve(pc,ksp);
1457: .ve

1459:    Note:
1460:    KSPSolve() calls this routine directly, so it is rarely called by the user.

1462:    Level: developer

1464: .seealso: PCPreSolve(), KSPSolve()
1465: @*/
1466: PetscErrorCode  PCPostSolve(PC pc,KSP ksp)
1467: {
1469:   Vec            x,rhs;

1474:   pc->presolvedone--;
1475:   KSPGetSolution(ksp,&x);
1476:   KSPGetRhs(ksp,&rhs);
1477:   if (pc->ops->postsolve) {
1478:      (*pc->ops->postsolve)(pc,ksp,rhs,x);
1479:   }
1480:   return(0);
1481: }

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

1486:   Collective on PetscViewer

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

1493:    Level: intermediate

1495:   Notes:
1496:    The type is determined by the data in the file, any type set into the PC before this call is ignored.

1498:   Notes for advanced users:
1499:   Most users should not need to know the details of the binary storage
1500:   format, since PCLoad() and PCView() completely hide these details.
1501:   But for anyone who's interested, the standard binary matrix storage
1502:   format is
1503: .vb
1504:      has not yet been determined
1505: .ve

1507: .seealso: PetscViewerBinaryOpen(), PCView(), MatLoad(), VecLoad()
1508: @*/
1509: PetscErrorCode  PCLoad(PC newdm, PetscViewer viewer)
1510: {
1512:   PetscBool      isbinary;
1513:   PetscInt       classid;
1514:   char           type[256];

1519:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1520:   if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");

1522:   PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
1523:   if (classid != PC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not PC next in file");
1524:   PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
1525:   PCSetType(newdm, type);
1526:   if (newdm->ops->load) {
1527:     (*newdm->ops->load)(newdm,viewer);
1528:   }
1529:   return(0);
1530: }

1532:  #include <petscdraw.h>
1533: #if defined(PETSC_HAVE_SAWS)
1534:  #include <petscviewersaws.h>
1535: #endif

1537: /*@C
1538:    PCViewFromOptions - View from Options

1540:    Collective on PC

1542:    Input Parameters:
1543: +  A - the PC context
1544: .  obj - Optional object
1545: -  name - command line option

1547:    Level: intermediate
1548: .seealso:  PC, PCView, PetscObjectViewFromOptions(), PCCreate()
1549: @*/
1550: PetscErrorCode  PCViewFromOptions(PC A,PetscObject obj,const char name[])
1551: {

1556:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
1557:   return(0);
1558: }

1560: /*@C
1561:    PCView - Prints the PC data structure.

1563:    Collective on PC

1565:    Input Parameters:
1566: +  PC - the PC context
1567: -  viewer - optional visualization context

1569:    Note:
1570:    The available visualization contexts include
1571: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
1572: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1573:          output where only the first processor opens
1574:          the file.  All other processors send their
1575:          data to the first processor to print.

1577:    The user can open an alternative visualization contexts with
1578:    PetscViewerASCIIOpen() (output to a specified file).

1580:    Level: developer

1582: .seealso: KSPView(), PetscViewerASCIIOpen()
1583: @*/
1584: PetscErrorCode  PCView(PC pc,PetscViewer viewer)
1585: {
1586:   PCType         cstr;
1588:   PetscBool      iascii,isstring,isbinary,isdraw;
1589: #if defined(PETSC_HAVE_SAWS)
1590:   PetscBool      issaws;
1591: #endif

1595:   if (!viewer) {
1596:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);
1597:   }

1601:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1602:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
1603:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1604:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1605: #if defined(PETSC_HAVE_SAWS)
1606:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws);
1607: #endif

1609:   if (iascii) {
1610:     PetscObjectPrintClassNamePrefixType((PetscObject)pc,viewer);
1611:     if (!pc->setupcalled) {
1612:       PetscViewerASCIIPrintf(viewer,"  PC has not been set up so information may be incomplete\n");
1613:     }
1614:     if (pc->ops->view) {
1615:       PetscViewerASCIIPushTab(viewer);
1616:       (*pc->ops->view)(pc,viewer);
1617:       PetscViewerASCIIPopTab(viewer);
1618:     }
1619:     if (pc->mat) {
1620:       PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
1621:       if (pc->pmat == pc->mat) {
1622:         PetscViewerASCIIPrintf(viewer,"  linear system matrix = precond matrix:\n");
1623:         PetscViewerASCIIPushTab(viewer);
1624:         MatView(pc->mat,viewer);
1625:         PetscViewerASCIIPopTab(viewer);
1626:       } else {
1627:         if (pc->pmat) {
1628:           PetscViewerASCIIPrintf(viewer,"  linear system matrix followed by preconditioner matrix:\n");
1629:         } else {
1630:           PetscViewerASCIIPrintf(viewer,"  linear system matrix:\n");
1631:         }
1632:         PetscViewerASCIIPushTab(viewer);
1633:         MatView(pc->mat,viewer);
1634:         if (pc->pmat) {MatView(pc->pmat,viewer);}
1635:         PetscViewerASCIIPopTab(viewer);
1636:       }
1637:       PetscViewerPopFormat(viewer);
1638:     }
1639:   } else if (isstring) {
1640:     PCGetType(pc,&cstr);
1641:     PetscViewerStringSPrintf(viewer," PCType: %-7.7s",cstr);
1642:     if (pc->ops->view) {(*pc->ops->view)(pc,viewer);}
1643:     if (pc->mat) {MatView(pc->mat,viewer);}
1644:     if (pc->pmat && pc->pmat != pc->mat) {MatView(pc->pmat,viewer);}
1645:   } else if (isbinary) {
1646:     PetscInt    classid = PC_FILE_CLASSID;
1647:     MPI_Comm    comm;
1648:     PetscMPIInt rank;
1649:     char        type[256];

1651:     PetscObjectGetComm((PetscObject)pc,&comm);
1652:     MPI_Comm_rank(comm,&rank);
1653:     if (!rank) {
1654:       PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT);
1655:       PetscStrncpy(type,((PetscObject)pc)->type_name,256);
1656:       PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR);
1657:     }
1658:     if (pc->ops->view) {
1659:       (*pc->ops->view)(pc,viewer);
1660:     }
1661:   } else if (isdraw) {
1662:     PetscDraw draw;
1663:     char      str[25];
1664:     PetscReal x,y,bottom,h;
1665:     PetscInt  n;

1667:     PetscViewerDrawGetDraw(viewer,0,&draw);
1668:     PetscDrawGetCurrentPoint(draw,&x,&y);
1669:     if (pc->mat) {
1670:       MatGetSize(pc->mat,&n,NULL);
1671:       PetscSNPrintf(str,25,"PC: %s (%D)",((PetscObject)pc)->type_name,n);
1672:     } else {
1673:       PetscSNPrintf(str,25,"PC: %s",((PetscObject)pc)->type_name);
1674:     }
1675:     PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
1676:     bottom = y - h;
1677:     PetscDrawPushCurrentPoint(draw,x,bottom);
1678:     if (pc->ops->view) {
1679:       (*pc->ops->view)(pc,viewer);
1680:     }
1681:     PetscDrawPopCurrentPoint(draw);
1682: #if defined(PETSC_HAVE_SAWS)
1683:   } else if (issaws) {
1684:     PetscMPIInt rank;

1686:     PetscObjectName((PetscObject)pc);
1687:     MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
1688:     if (!((PetscObject)pc)->amsmem && !rank) {
1689:       PetscObjectViewSAWs((PetscObject)pc,viewer);
1690:     }
1691:     if (pc->mat) {MatView(pc->mat,viewer);}
1692:     if (pc->pmat && pc->pmat != pc->mat) {MatView(pc->pmat,viewer);}
1693: #endif
1694:   }
1695:   return(0);
1696: }

1698: /*@C
1699:   PCRegister -  Adds a method to the preconditioner package.

1701:    Not collective

1703:    Input Parameters:
1704: +  name_solver - name of a new user-defined solver
1705: -  routine_create - routine to create method context

1707:    Notes:
1708:    PCRegister() may be called multiple times to add several user-defined preconditioners.

1710:    Sample usage:
1711: .vb
1712:    PCRegister("my_solver", MySolverCreate);
1713: .ve

1715:    Then, your solver can be chosen with the procedural interface via
1716: $     PCSetType(pc,"my_solver")
1717:    or at runtime via the option
1718: $     -pc_type my_solver

1720:    Level: advanced

1722: .seealso: PCRegisterAll()
1723: @*/
1724: PetscErrorCode  PCRegister(const char sname[],PetscErrorCode (*function)(PC))
1725: {

1729:   PCInitializePackage();
1730:   PetscFunctionListAdd(&PCList,sname,function);
1731:   return(0);
1732: }

1734: static PetscErrorCode MatMult_PC(Mat A,Vec X,Vec Y)
1735: {
1736:   PC             pc;

1740:   MatShellGetContext(A,&pc);
1741:   PCApply(pc,X,Y);
1742:   return(0);
1743: }

1745: /*@
1746:     PCComputeOperator - Computes the explicit preconditioned operator.

1748:     Collective on PC

1750:     Input Parameter:
1751: +   pc - the preconditioner object
1752: -   mattype - the matrix type to be used for the operator

1754:     Output Parameter:
1755: .   mat - the explict preconditioned operator

1757:     Notes:
1758:     This computation is done by applying the operators to columns of the identity matrix.
1759:     This routine is costly in general, and is recommended for use only with relatively small systems.
1760:     Currently, this routine uses a dense matrix format when mattype == NULL

1762:     Level: advanced

1764: .seealso: KSPComputeOperator(), MatType

1766: @*/
1767: PetscErrorCode  PCComputeOperator(PC pc,MatType mattype,Mat *mat)
1768: {
1770:   PetscInt       N,M,m,n;
1771:   Mat            A,Apc;

1776:   PCGetOperators(pc,&A,NULL);
1777:   MatGetLocalSize(A,&m,&n);
1778:   MatGetSize(A,&M,&N);
1779:   MatCreateShell(PetscObjectComm((PetscObject)pc),m,n,M,N,pc,&Apc);
1780:   MatShellSetOperation(Apc,MATOP_MULT,(void (*)(void))MatMult_PC);
1781:   MatComputeOperator(Apc,mattype,mat);
1782:   MatDestroy(&Apc);
1783:   return(0);
1784: }

1786: /*@
1787:    PCSetCoordinates - sets the coordinates of all the nodes on the local process

1789:    Collective on PC

1791:    Input Parameters:
1792: +  pc - the solver context
1793: .  dim - the dimension of the coordinates 1, 2, or 3
1794: .  nloc - the blocked size of the coordinates array
1795: -  coords - the coordinates array

1797:    Level: intermediate

1799:    Notes:
1800:    coords is an array of the dim coordinates for the nodes on
1801:    the local processor, of size dim*nloc.
1802:    If there are 108 equation on a processor
1803:    for a displacement finite element discretization of elasticity (so
1804:    that there are nloc = 36 = 108/3 nodes) then the array must have 108
1805:    double precision values (ie, 3 * 36).  These x y z coordinates
1806:    should be ordered for nodes 0 to N-1 like so: [ 0.x, 0.y, 0.z, 1.x,
1807:    ... , N-1.z ].

1809: .seealso: MatSetNearNullSpace()
1810: @*/
1811: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[])
1812: {

1818:   PetscTryMethod(pc,"PCSetCoordinates_C",(PC,PetscInt,PetscInt,PetscReal*),(pc,dim,nloc,coords));
1819:   return(0);
1820: }

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

1825:    Logically Collective on PC

1827:    Input Parameters:
1828: +  pc - the precondition context

1830:    Output Parameter:
1831: -  num_levels - the number of levels
1832: .  interpolations - the interpolation matrices (size of num_levels-1)

1834:    Level: advanced

1836: .keywords: MG, GAMG, BoomerAMG, multigrid, interpolation, level

1838: .seealso: PCMGGetRestriction(), PCMGSetInterpolation(), PCMGGetInterpolation(), PCGetCoarseOperators()
1839: @*/
1840: PetscErrorCode PCGetInterpolations(PC pc,PetscInt *num_levels,Mat *interpolations[])
1841: {

1848:   PetscUseMethod(pc,"PCGetInterpolations_C",(PC,PetscInt*,Mat*[]),(pc,num_levels,interpolations));
1849:   return(0);
1850: }

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

1855:    Logically Collective on PC

1857:    Input Parameters:
1858: +  pc - the precondition context

1860:    Output Parameter:
1861: -  num_levels - the number of levels
1862: .  coarseOperators - the coarse operator matrices (size of num_levels-1)

1864:    Level: advanced

1866: .keywords: MG, GAMG, BoomerAMG, get, multigrid, interpolation, level

1868: .seealso: PCMGGetRestriction(), PCMGSetInterpolation(), PCMGGetRScale(), PCMGGetInterpolation(), PCGetInterpolations()
1869: @*/
1870: PetscErrorCode PCGetCoarseOperators(PC pc,PetscInt *num_levels,Mat *coarseOperators[])
1871: {

1878:   PetscUseMethod(pc,"PCGetCoarseOperators_C",(PC,PetscInt*,Mat*[]),(pc,num_levels,coarseOperators));
1879:   return(0);
1880: }