Actual source code: precon.c

petsc-3.11.4 2019-09-28
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: .keywords: PC, destroy

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

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

 85:   pc->setupcalled = 0;
 86:   return(0);
 87: }

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

 92:    Collective on PC

 94:    Input Parameter:
 95: .  pc - the preconditioner context

 97:    Level: developer

 99: .keywords: PC, destroy

101: .seealso: PCCreate(), PCSetUp()
102: @*/
103: PetscErrorCode  PCDestroy(PC *pc)
104: {

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

112:   PCReset(*pc);

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

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

126:    Logically Collective on PC

128:    Input Parameter:
129: .  pc - the preconditioner context

131:    Output Parameter:
132: .  flag - PETSC_TRUE if it applies the scaling

134:    Level: developer

136:    Notes:
137:     If this returns PETSC_TRUE then the system solved via the Krylov method is
138: $           D M A D^{-1} y = D M b  for left preconditioning or
139: $           D A M D^{-1} z = D b for right preconditioning

141: .keywords: PC

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

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

158:    Logically Collective on PC

160:    Input Parameters:
161: +  pc - the preconditioner context
162: -  s - scaling vector

164:    Level: intermediate

166:    Notes:
167:     The system solved via the Krylov method is
168: $           D M A D^{-1} y = D M b  for left preconditioning or
169: $           D A M D^{-1} z = D b for right preconditioning

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

173: .keywords: PC

175: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCGetDiagonalScale()
176: @*/
177: PetscErrorCode  PCSetDiagonalScale(PC pc,Vec s)
178: {

184:   pc->diagonalscale     = PETSC_TRUE;

186:   PetscObjectReference((PetscObject)s);
187:   VecDestroy(&pc->diagonalscaleleft);

189:   pc->diagonalscaleleft = s;

191:   VecDuplicate(s,&pc->diagonalscaleright);
192:   VecCopy(s,pc->diagonalscaleright);
193:   VecReciprocal(pc->diagonalscaleright);
194:   return(0);
195: }

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

200:    Logically Collective on PC

202:    Input Parameters:
203: +  pc - the preconditioner context
204: .  in - input vector
205: +  out - scaled vector (maybe the same as in)

207:    Level: intermediate

209:    Notes:
210:     The system solved via the Krylov method is
211: $           D M A D^{-1} y = D M b  for left preconditioning or
212: $           D A M D^{-1} z = D b for right preconditioning

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

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

218: .keywords: PC

220: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleSet(), PCDiagonalScaleRight(), PCDiagonalScale()
221: @*/
222: PetscErrorCode  PCDiagonalScaleLeft(PC pc,Vec in,Vec out)
223: {

230:   if (pc->diagonalscale) {
231:     VecPointwiseMult(out,pc->diagonalscaleleft,in);
232:   } else if (in != out) {
233:     VecCopy(in,out);
234:   }
235:   return(0);
236: }

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

241:    Logically Collective on PC

243:    Input Parameters:
244: +  pc - the preconditioner context
245: .  in - input vector
246: +  out - scaled vector (maybe the same as in)

248:    Level: intermediate

250:    Notes:
251:     The system solved via the Krylov method is
252: $           D M A D^{-1} y = D M b  for left preconditioning or
253: $           D A M D^{-1} z = D b for right preconditioning

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

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

259: .keywords: PC

261: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleSet(), PCDiagonalScale()
262: @*/
263: PetscErrorCode  PCDiagonalScaleRight(PC pc,Vec in,Vec out)
264: {

271:   if (pc->diagonalscale) {
272:     VecPointwiseMult(out,pc->diagonalscaleright,in);
273:   } else if (in != out) {
274:     VecCopy(in,out);
275:   }
276:   return(0);
277: }

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

284:    Logically Collective on PC

286:    Input Parameters:
287: +  pc - the preconditioner context
288: -  flg - PETSC_TRUE to use the Amat, PETSC_FALSE to use the Pmat (default is false)

290:    Options Database Key:
291: .  -pc_use_amat <true,false>

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

297:    Level: intermediate

299: .seealso: PCGetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
300: @*/
301: PetscErrorCode  PCSetUseAmat(PC pc,PetscBool flg)
302: {
305:   pc->useAmat = flg;
306:   return(0);
307: }

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

312:    Logically Collective on PC

314:    Input Parameters:
315: +  pc - iterative context obtained from PCCreate()
316: -  flg - PETSC_TRUE indicates you want the error generated

318:    Level: advanced

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

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

327: .keywords: PC, set, initial guess, nonzero

329: .seealso: PCGetInitialGuessNonzero(), PCSetInitialGuessKnoll(), PCGetInitialGuessKnoll()
330: @*/
331: PetscErrorCode  PCSetErrorIfFailure(PC pc,PetscBool flg)
332: {
336:   pc->erroriffailure = flg;
337:   return(0);
338: }

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

345:    Logically Collective on PC

347:    Input Parameter:
348: .  pc - the preconditioner context

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

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

357:    Level: intermediate

359: .seealso: PCSetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
360: @*/
361: PetscErrorCode  PCGetUseAmat(PC pc,PetscBool *flg)
362: {
365:   *flg = pc->useAmat;
366:   return(0);
367: }

369: /*@
370:    PCCreate - Creates a preconditioner context.

372:    Collective on MPI_Comm

374:    Input Parameter:
375: .  comm - MPI communicator

377:    Output Parameter:
378: .  pc - location to put the preconditioner context

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

384:    Level: developer

386: .keywords: PC, create, context

388: .seealso: PCSetUp(), PCApply(), PCDestroy()
389: @*/
390: PetscErrorCode  PCCreate(MPI_Comm comm,PC *newpc)
391: {
392:   PC             pc;

397:   *newpc = 0;
398:   PCInitializePackage();

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

402:   pc->mat                  = 0;
403:   pc->pmat                 = 0;
404:   pc->setupcalled          = 0;
405:   pc->setfromoptionscalled = 0;
406:   pc->data                 = 0;
407:   pc->diagonalscale        = PETSC_FALSE;
408:   pc->diagonalscaleleft    = 0;
409:   pc->diagonalscaleright   = 0;

411:   pc->modifysubmatrices  = 0;
412:   pc->modifysubmatricesP = 0;

414:   *newpc = pc;
415:   return(0);

417: }

419: /* -------------------------------------------------------------------------------*/

421: /*@
422:    PCApply - Applies the preconditioner to a vector.

424:    Collective on PC and Vec

426:    Input Parameters:
427: +  pc - the preconditioner context
428: -  x - input vector

430:    Output Parameter:
431: .  y - output vector

433:    Level: developer

435: .keywords: PC, apply

437: .seealso: PCApplyTranspose(), PCApplyBAorAB()
438: @*/
439: PetscErrorCode  PCApply(PC pc,Vec x,Vec y)
440: {
442:   PetscInt       m,n,mv,nv;

448:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
449:   if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
450:   /* use pmat to check vector sizes since for KSPLQR the pmat may be of a different size than mat */
451:   MatGetLocalSize(pc->pmat,&m,&n);
452:   VecGetLocalSize(x,&nv);
453:   VecGetLocalSize(y,&mv);
454:   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);
455:   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);
456:   VecSetErrorIfLocked(y,3);

458:   PCSetUp(pc);
459:   if (!pc->ops->apply) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply");
460:   VecLockReadPush(x);
461:   PetscLogEventBegin(PC_Apply,pc,x,y,0);
462:   (*pc->ops->apply)(pc,x,y);
463:   PetscLogEventEnd(PC_Apply,pc,x,y,0);
464:   if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
465:   VecLockReadPop(x);
466:   return(0);
467: }

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

472:    Collective on PC and Vec

474:    Input Parameters:
475: +  pc - the preconditioner context
476: -  x - input vector

478:    Output Parameter:
479: .  y - output vector

481:    Notes:
482:    Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.

484:    Level: developer

486: .keywords: PC, apply, symmetric, left

488: .seealso: PCApply(), PCApplySymmetricRight()
489: @*/
490: PetscErrorCode  PCApplySymmetricLeft(PC pc,Vec x,Vec y)
491: {

498:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
499:   if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
500:   PCSetUp(pc);
501:   if (!pc->ops->applysymmetricleft) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have left symmetric apply");
502:   VecLockReadPush(x);
503:   PetscLogEventBegin(PC_ApplySymmetricLeft,pc,x,y,0);
504:   (*pc->ops->applysymmetricleft)(pc,x,y);
505:   PetscLogEventEnd(PC_ApplySymmetricLeft,pc,x,y,0);
506:   VecLockReadPop(x);
507:   if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
508:   return(0);
509: }

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

514:    Collective on PC and Vec

516:    Input Parameters:
517: +  pc - the preconditioner context
518: -  x - input vector

520:    Output Parameter:
521: .  y - output vector

523:    Level: developer

525:    Notes:
526:    Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.

528: .keywords: PC, apply, symmetric, right

530: .seealso: PCApply(), PCApplySymmetricLeft()
531: @*/
532: PetscErrorCode  PCApplySymmetricRight(PC pc,Vec x,Vec y)
533: {

540:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
541:   if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
542:   PCSetUp(pc);
543:   if (!pc->ops->applysymmetricright) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have left symmetric apply");
544:   VecLockReadPush(x);
545:   PetscLogEventBegin(PC_ApplySymmetricRight,pc,x,y,0);
546:   (*pc->ops->applysymmetricright)(pc,x,y);
547:   PetscLogEventEnd(PC_ApplySymmetricRight,pc,x,y,0);
548:   VecLockReadPop(x);
549:   if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
550:   return(0);
551: }

553: /*@
554:    PCApplyTranspose - Applies the transpose of preconditioner to a vector.

556:    Collective on PC and Vec

558:    Input Parameters:
559: +  pc - the preconditioner context
560: -  x - input vector

562:    Output Parameter:
563: .  y - output vector

565:    Notes:
566:     For complex numbers this applies the non-Hermitian transpose.

568:    Developer Notes:
569:     We need to implement a PCApplyHermitianTranspose()

571:    Level: developer

573: .keywords: PC, apply, transpose

575: .seealso: PCApply(), PCApplyBAorAB(), PCApplyBAorABTranspose(), PCApplyTransposeExists()
576: @*/
577: PetscErrorCode  PCApplyTranspose(PC pc,Vec x,Vec y)
578: {

585:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
586:   if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
587:   PCSetUp(pc);
588:   if (!pc->ops->applytranspose) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply transpose");
589:   VecLockReadPush(x);
590:   PetscLogEventBegin(PC_Apply,pc,x,y,0);
591:   (*pc->ops->applytranspose)(pc,x,y);
592:   PetscLogEventEnd(PC_Apply,pc,x,y,0);
593:   VecLockReadPop(x);
594:   if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
595:   return(0);
596: }

598: /*@
599:    PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation

601:    Collective on PC and Vec

603:    Input Parameters:
604: .  pc - the preconditioner context

606:    Output Parameter:
607: .  flg - PETSC_TRUE if a transpose operation is defined

609:    Level: developer

611: .keywords: PC, apply, transpose

613: .seealso: PCApplyTranspose()
614: @*/
615: PetscErrorCode  PCApplyTransposeExists(PC pc,PetscBool  *flg)
616: {
620:   if (pc->ops->applytranspose) *flg = PETSC_TRUE;
621:   else *flg = PETSC_FALSE;
622:   return(0);
623: }

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

628:    Collective on PC and Vec

630:    Input Parameters:
631: +  pc - the preconditioner context
632: .  side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
633: .  x - input vector
634: -  work - work vector

636:    Output Parameter:
637: .  y - output vector

639:    Level: developer

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

645: .keywords: PC, apply, operator

647: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorABTranspose()
648: @*/
649: PetscErrorCode  PCApplyBAorAB(PC pc,PCSide side,Vec x,Vec y,Vec work)
650: {

658:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
659:   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");
660:   if (pc->diagonalscale && side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot include diagonal scaling with symmetric preconditioner application");
661:   if (pc->erroriffailure) {VecValidValues(x,3,PETSC_TRUE);}

663:   PCSetUp(pc);
664:   if (pc->diagonalscale) {
665:     if (pc->ops->applyBA) {
666:       Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
667:       VecDuplicate(x,&work2);
668:       PCDiagonalScaleRight(pc,x,work2);
669:       (*pc->ops->applyBA)(pc,side,work2,y,work);
670:       PCDiagonalScaleLeft(pc,y,y);
671:       VecDestroy(&work2);
672:     } else if (side == PC_RIGHT) {
673:       PCDiagonalScaleRight(pc,x,y);
674:       PCApply(pc,y,work);
675:       MatMult(pc->mat,work,y);
676:       PCDiagonalScaleLeft(pc,y,y);
677:     } else if (side == PC_LEFT) {
678:       PCDiagonalScaleRight(pc,x,y);
679:       MatMult(pc->mat,y,work);
680:       PCApply(pc,work,y);
681:       PCDiagonalScaleLeft(pc,y,y);
682:     } else if (side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot provide diagonal scaling with symmetric application of preconditioner");
683:   } else {
684:     if (pc->ops->applyBA) {
685:       (*pc->ops->applyBA)(pc,side,x,y,work);
686:     } else if (side == PC_RIGHT) {
687:       PCApply(pc,x,work);
688:       MatMult(pc->mat,work,y);
689:     } else if (side == PC_LEFT) {
690:       MatMult(pc->mat,x,work);
691:       PCApply(pc,work,y);
692:     } else if (side == PC_SYMMETRIC) {
693:       /* There's an extra copy here; maybe should provide 2 work vectors instead? */
694:       PCApplySymmetricRight(pc,x,work);
695:       MatMult(pc->mat,work,y);
696:       VecCopy(y,work);
697:       PCApplySymmetricLeft(pc,work,y);
698:     }
699:   }
700:   if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
701:   return(0);
702: }

704: /*@
705:    PCApplyBAorABTranspose - Applies the transpose of the preconditioner
706:    and operator to a vector. That is, applies tr(B) * tr(A) with left preconditioning,
707:    NOT tr(B*A) = tr(A)*tr(B).

709:    Collective on PC and Vec

711:    Input Parameters:
712: +  pc - the preconditioner context
713: .  side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
714: .  x - input vector
715: -  work - work vector

717:    Output Parameter:
718: .  y - output vector


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

725:     Level: developer

727: .keywords: PC, apply, operator, transpose

729: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorAB()
730: @*/
731: PetscErrorCode  PCApplyBAorABTranspose(PC pc,PCSide side,Vec x,Vec y,Vec work)
732: {

740:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
741:   if (pc->erroriffailure) {VecValidValues(x,3,PETSC_TRUE);}
742:   if (pc->ops->applyBAtranspose) {
743:     (*pc->ops->applyBAtranspose)(pc,side,x,y,work);
744:     if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
745:     return(0);
746:   }
747:   if (side != PC_LEFT && side != PC_RIGHT) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Side must be right or left");

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(y,4,PETSC_FALSE);}
759:   return(0);
760: }

762: /* -------------------------------------------------------------------------------*/

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

768:    Not Collective

770:    Input Parameter:
771: .  pc - the preconditioner

773:    Output Parameter:
774: .  exists - PETSC_TRUE or PETSC_FALSE

776:    Level: developer

778: .keywords: PC, apply, Richardson, exists

780: .seealso: PCApplyRichardson()
781: @*/
782: PetscErrorCode  PCApplyRichardsonExists(PC pc,PetscBool  *exists)
783: {
787:   if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
788:   else *exists = PETSC_FALSE;
789:   return(0);
790: }

792: /*@
793:    PCApplyRichardson - Applies several steps of Richardson iteration with
794:    the particular preconditioner. This routine is usually used by the
795:    Krylov solvers and not the application code directly.

797:    Collective on PC

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

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

814:    Notes:
815:    Most preconditioners do not support this function. Use the command
816:    PCApplyRichardsonExists() to determine if one does.

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

821:    Level: developer

823: .keywords: PC, apply, Richardson

825: .seealso: PCApplyRichardsonExists()
826: @*/
827: PetscErrorCode  PCApplyRichardson(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
828: {

836:   if (b == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"b and y must be different vectors");
837:   PCSetUp(pc);
838:   if (!pc->ops->applyrichardson) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply richardson");
839:   (*pc->ops->applyrichardson)(pc,b,y,w,rtol,abstol,dtol,its,guesszero,outits,reason);
840:   return(0);
841: }

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

846:    Logically Collective on PC

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

851:    Output Parameter:
852: .  reason - the reason it failed, currently only -1

854:    Level: advanced

856: .keywords: PC, setup

858: .seealso: PCCreate(), PCApply(), PCDestroy()
859: @*/
860: PetscErrorCode PCGetFailedReason(PC pc,PCFailedReason *reason)
861: {
863:   if (pc->setupcalled < 0) *reason = (PCFailedReason)pc->setupcalled;
864:   else *reason = pc->failedreason;
865:   return(0);
866: }


869: /*
870:       a setupcall of 0 indicates never setup,
871:                      1 indicates has been previously setup
872:                     -1 indicates a PCSetUp() was attempted and failed
873: */
874: /*@
875:    PCSetUp - Prepares for the use of a preconditioner.

877:    Collective on PC

879:    Input Parameter:
880: .  pc - the preconditioner context

882:    Level: developer

884: .keywords: PC, setup

886: .seealso: PCCreate(), PCApply(), PCDestroy()
887: @*/
888: PetscErrorCode  PCSetUp(PC pc)
889: {
890:   PetscErrorCode   ierr;
891:   const char       *def;
892:   PetscObjectState matstate, matnonzerostate;

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

898:   if (pc->setupcalled && pc->reusepreconditioner) {
899:     PetscInfo(pc,"Leaving PC with identical preconditioner since reuse preconditioner is set\n");
900:     return(0);
901:   }

903:   PetscObjectStateGet((PetscObject)pc->pmat,&matstate);
904:   MatGetNonzeroState(pc->pmat,&matnonzerostate);
905:   if (!pc->setupcalled) {
906:     PetscInfo(pc,"Setting up PC for first time\n");
907:     pc->flag        = DIFFERENT_NONZERO_PATTERN;
908:   } else if (matstate == pc->matstate) {
909:     PetscInfo(pc,"Leaving PC with identical preconditioner since operator is unchanged\n");
910:     return(0);
911:   } else {
912:     if (matnonzerostate > pc->matnonzerostate) {
913:        PetscInfo(pc,"Setting up PC with different nonzero pattern\n");
914:        pc->flag            = DIFFERENT_NONZERO_PATTERN;
915:     } else {
916:       PetscInfo(pc,"Setting up PC with same nonzero pattern\n");
917:       pc->flag            = SAME_NONZERO_PATTERN;
918:     }
919:   }
920:   pc->matstate        = matstate;
921:   pc->matnonzerostate = matnonzerostate;

923:   if (!((PetscObject)pc)->type_name) {
924:     PCGetDefaultType_Private(pc,&def);
925:     PCSetType(pc,def);
926:   }

928:   MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
929:   MatSetErrorIfFailure(pc->mat,pc->erroriffailure);
930:   PetscLogEventBegin(PC_SetUp,pc,0,0,0);
931:   if (pc->ops->setup) {
932:     (*pc->ops->setup)(pc);
933:   }
934:   PetscLogEventEnd(PC_SetUp,pc,0,0,0);
935:   if (!pc->setupcalled) pc->setupcalled = 1;
936:   return(0);
937: }

939: /*@
940:    PCSetUpOnBlocks - Sets up the preconditioner for each block in
941:    the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
942:    methods.

944:    Collective on PC

946:    Input Parameters:
947: .  pc - the preconditioner context

949:    Level: developer

951: .keywords: PC, setup, blocks

953: .seealso: PCCreate(), PCApply(), PCDestroy(), PCSetUp()
954: @*/
955: PetscErrorCode  PCSetUpOnBlocks(PC pc)
956: {

961:   if (!pc->ops->setuponblocks) return(0);
962:   PetscLogEventBegin(PC_SetUpOnBlocks,pc,0,0,0);
963:   (*pc->ops->setuponblocks)(pc);
964:   PetscLogEventEnd(PC_SetUpOnBlocks,pc,0,0,0);
965:   return(0);
966: }

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

975:    Logically Collective on PC

977:    Input Parameters:
978: +  pc - the preconditioner context
979: .  func - routine for modifying the submatrices
980: -  ctx - optional user-defined context (may be null)

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

985: .  row - an array of index sets that contain the global row numbers
986:          that comprise each local submatrix
987: .  col - an array of index sets that contain the global column numbers
988:          that comprise each local submatrix
989: .  submat - array of local submatrices
990: -  ctx - optional user-defined context for private data for the
991:          user-defined func routine (may be null)

993:    Notes:
994:    PCSetModifySubMatrices() MUST be called before KSPSetUp() and
995:    KSPSolve().

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

1001:    Level: advanced

1003: .keywords: PC, set, modify, submatrices

1005: .seealso: PCModifySubMatrices(), PCASMGetSubMatrices()
1006: @*/
1007: PetscErrorCode  PCSetModifySubMatrices(PC pc,PetscErrorCode (*func)(PC,PetscInt,const IS[],const IS[],Mat[],void*),void *ctx)
1008: {
1011:   pc->modifysubmatrices  = func;
1012:   pc->modifysubmatricesP = ctx;
1013:   return(0);
1014: }

1016: /*@C
1017:    PCModifySubMatrices - Calls an optional user-defined routine within
1018:    certain preconditioners if one has been set with PCSetModifySubMarices().

1020:    Collective on PC

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

1033:    Output Parameter:
1034: .  submat - array of local submatrices (the entries of which may
1035:             have been modified)

1037:    Notes:
1038:    The user should NOT generally call this routine, as it will
1039:    automatically be called within certain preconditioners (currently
1040:    block Jacobi, additive Schwarz) if set.

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

1047:    Level: developer

1049: .keywords: PC, modify, submatrices

1051: .seealso: PCSetModifySubMatrices()
1052: @*/
1053: PetscErrorCode  PCModifySubMatrices(PC pc,PetscInt nsub,const IS row[],const IS col[],Mat submat[],void *ctx)
1054: {

1059:   if (!pc->modifysubmatrices) return(0);
1060:   PetscLogEventBegin(PC_ModifySubMatrices,pc,0,0,0);
1061:   (*pc->modifysubmatrices)(pc,nsub,row,col,submat,ctx);
1062:   PetscLogEventEnd(PC_ModifySubMatrices,pc,0,0,0);
1063:   return(0);
1064: }

1066: /*@
1067:    PCSetOperators - Sets the matrix associated with the linear system and
1068:    a (possibly) different one associated with the preconditioner.

1070:    Logically Collective on PC and Mat

1072:    Input Parameters:
1073: +  pc - the preconditioner context
1074: .  Amat - the matrix that defines the linear system
1075: -  Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.

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

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

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

1090:    Level: intermediate

1092: .keywords: PC, set, operators, matrix, linear system

1094: .seealso: PCGetOperators(), MatZeroEntries()
1095:  @*/
1096: PetscErrorCode  PCSetOperators(PC pc,Mat Amat,Mat Pmat)
1097: {
1098:   PetscErrorCode   ierr;
1099:   PetscInt         m1,n1,m2,n2;

1107:   if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
1108:     MatGetLocalSize(Amat,&m1,&n1);
1109:     MatGetLocalSize(pc->mat,&m2,&n2);
1110:     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);
1111:     MatGetLocalSize(Pmat,&m1,&n1);
1112:     MatGetLocalSize(pc->pmat,&m2,&n2);
1113:     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);
1114:   }

1116:   if (Pmat != pc->pmat) {
1117:     /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
1118:     pc->matnonzerostate = -1;
1119:     pc->matstate        = -1;
1120:   }

1122:   /* reference first in case the matrices are the same */
1123:   if (Amat) {PetscObjectReference((PetscObject)Amat);}
1124:   MatDestroy(&pc->mat);
1125:   if (Pmat) {PetscObjectReference((PetscObject)Pmat);}
1126:   MatDestroy(&pc->pmat);
1127:   pc->mat  = Amat;
1128:   pc->pmat = Pmat;
1129:   return(0);
1130: }

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

1135:    Logically Collective on PC

1137:    Input Parameters:
1138: +  pc - the preconditioner context
1139: -  flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner

1141:     Level: intermediate

1143: .seealso: PCGetOperators(), MatZeroEntries(), PCGetReusePreconditioner(), KSPSetReusePreconditioner()
1144:  @*/
1145: PetscErrorCode  PCSetReusePreconditioner(PC pc,PetscBool flag)
1146: {
1149:   pc->reusepreconditioner = flag;
1150:   return(0);
1151: }

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

1156:    Not Collective

1158:    Input Parameter:
1159: .  pc - the preconditioner context

1161:    Output Parameter:
1162: .  flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner

1164:    Level: intermediate

1166: .seealso: PCGetOperators(), MatZeroEntries(), PCSetReusePreconditioner()
1167:  @*/
1168: PetscErrorCode  PCGetReusePreconditioner(PC pc,PetscBool *flag)
1169: {
1172:   *flag = pc->reusepreconditioner;
1173:   return(0);
1174: }

1176: /*@
1177:    PCGetOperators - Gets the matrix associated with the linear system and
1178:    possibly a different one associated with the preconditioner.

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

1182:    Input Parameter:
1183: .  pc - the preconditioner context

1185:    Output Parameters:
1186: +  Amat - the matrix defining the linear system
1187: -  Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.

1189:    Level: intermediate

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

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

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

1205: $         MatCreate(comm,&mat);
1206: $         KSP/PCSetOperators(ksp/pc,Amat,Amat);
1207: $         PetscObjectDereference((PetscObject)mat);
1208: $           set size, type, etc of Amat

1210:      and

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

1215: $         MatCreate(comm,&Amat);
1216: $         MatCreate(comm,&Pmat);
1217: $         KSP/PCSetOperators(ksp/pc,Amat,Pmat);
1218: $         PetscObjectDereference((PetscObject)Amat);
1219: $         PetscObjectDereference((PetscObject)Pmat);
1220: $           set size, type, etc of Amat and Pmat

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


1232: .keywords: PC, get, operators, matrix, linear system

1234: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperatorsSet()
1235: @*/
1236: PetscErrorCode  PCGetOperators(PC pc,Mat *Amat,Mat *Pmat)
1237: {

1242:   if (Amat) {
1243:     if (!pc->mat) {
1244:       if (pc->pmat && !Pmat) {  /* Apmat has been set, but user did not request it, so use for Amat */
1245:         pc->mat = pc->pmat;
1246:         PetscObjectReference((PetscObject)pc->mat);
1247:       } else {                  /* both Amat and Pmat are empty */
1248:         MatCreate(PetscObjectComm((PetscObject)pc),&pc->mat);
1249:         if (!Pmat) { /* user did NOT request Pmat, so make same as Amat */
1250:           pc->pmat = pc->mat;
1251:           PetscObjectReference((PetscObject)pc->pmat);
1252:         }
1253:       }
1254:     }
1255:     *Amat = pc->mat;
1256:   }
1257:   if (Pmat) {
1258:     if (!pc->pmat) {
1259:       if (pc->mat && !Amat) {    /* Amat has been set but was not requested, so use for pmat */
1260:         pc->pmat = pc->mat;
1261:         PetscObjectReference((PetscObject)pc->pmat);
1262:       } else {
1263:         MatCreate(PetscObjectComm((PetscObject)pc),&pc->pmat);
1264:         if (!Amat) { /* user did NOT request Amat, so make same as Pmat */
1265:           pc->mat = pc->pmat;
1266:           PetscObjectReference((PetscObject)pc->mat);
1267:         }
1268:       }
1269:     }
1270:     *Pmat = pc->pmat;
1271:   }
1272:   return(0);
1273: }

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

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

1281:    Input Parameter:
1282: .  pc - the preconditioner context

1284:    Output Parameters:
1285: +  mat - the matrix associated with the linear system was set
1286: -  pmat - matrix associated with the preconditioner was set, usually the same

1288:    Level: intermediate

1290: .keywords: PC, get, operators, matrix, linear system

1292: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperators()
1293: @*/
1294: PetscErrorCode  PCGetOperatorsSet(PC pc,PetscBool  *mat,PetscBool  *pmat)
1295: {
1298:   if (mat) *mat = (pc->mat)  ? PETSC_TRUE : PETSC_FALSE;
1299:   if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1300:   return(0);
1301: }

1303: /*@
1304:    PCFactorGetMatrix - Gets the factored matrix from the
1305:    preconditioner context.  This routine is valid only for the LU,
1306:    incomplete LU, Cholesky, and incomplete Cholesky methods.

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

1310:    Input Parameters:
1311: .  pc - the preconditioner context

1313:    Output parameters:
1314: .  mat - the factored matrix

1316:    Level: advanced

1318:    Notes:
1319:     Does not increase the reference count for the matrix so DO NOT destroy it

1321: .keywords: PC, get, factored, matrix
1322: @*/
1323: PetscErrorCode  PCFactorGetMatrix(PC pc,Mat *mat)
1324: {

1330:   if (pc->ops->getfactoredmatrix) {
1331:     (*pc->ops->getfactoredmatrix)(pc,mat);
1332:   } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC type does not support getting factor matrix");
1333:   return(0);
1334: }

1336: /*@C
1337:    PCSetOptionsPrefix - Sets the prefix used for searching for all
1338:    PC options in the database.

1340:    Logically Collective on PC

1342:    Input Parameters:
1343: +  pc - the preconditioner context
1344: -  prefix - the prefix string to prepend to all PC option requests

1346:    Notes:
1347:    A hyphen (-) must NOT be given at the beginning of the prefix name.
1348:    The first character of all runtime options is AUTOMATICALLY the
1349:    hyphen.

1351:    Level: advanced

1353: .keywords: PC, set, options, prefix, database

1355: .seealso: PCAppendOptionsPrefix(), PCGetOptionsPrefix()
1356: @*/
1357: PetscErrorCode  PCSetOptionsPrefix(PC pc,const char prefix[])
1358: {

1363:   PetscObjectSetOptionsPrefix((PetscObject)pc,prefix);
1364:   return(0);
1365: }

1367: /*@C
1368:    PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1369:    PC options in the database.

1371:    Logically Collective on PC

1373:    Input Parameters:
1374: +  pc - the preconditioner context
1375: -  prefix - the prefix string to prepend to all PC option requests

1377:    Notes:
1378:    A hyphen (-) must NOT be given at the beginning of the prefix name.
1379:    The first character of all runtime options is AUTOMATICALLY the
1380:    hyphen.

1382:    Level: advanced

1384: .keywords: PC, append, options, prefix, database

1386: .seealso: PCSetOptionsPrefix(), PCGetOptionsPrefix()
1387: @*/
1388: PetscErrorCode  PCAppendOptionsPrefix(PC pc,const char prefix[])
1389: {

1394:   PetscObjectAppendOptionsPrefix((PetscObject)pc,prefix);
1395:   return(0);
1396: }

1398: /*@C
1399:    PCGetOptionsPrefix - Gets the prefix used for searching for all
1400:    PC options in the database.

1402:    Not Collective

1404:    Input Parameters:
1405: .  pc - the preconditioner context

1407:    Output Parameters:
1408: .  prefix - pointer to the prefix string used, is returned

1410:    Notes:
1411:     On the fortran side, the user should pass in a string 'prifix' of
1412:    sufficient length to hold the prefix.

1414:    Level: advanced

1416: .keywords: PC, get, options, prefix, database

1418: .seealso: PCSetOptionsPrefix(), PCAppendOptionsPrefix()
1419: @*/
1420: PetscErrorCode  PCGetOptionsPrefix(PC pc,const char *prefix[])
1421: {

1427:   PetscObjectGetOptionsPrefix((PetscObject)pc,prefix);
1428:   return(0);
1429: }

1431: /*
1432:    Indicates the right hand side will be changed by KSPSolve(), this occurs for a few
1433:   preconditioners including BDDC and Eisentat that transform the equations before applying
1434:   the Krylov methods
1435: */
1436: PETSC_INTERN PetscErrorCode  PCPreSolveChangeRHS(PC pc,PetscBool *change)
1437: {

1443:   *change = PETSC_FALSE;
1444:   PetscTryMethod(pc,"PCPreSolveChangeRHS_C",(PC,PetscBool*),(pc,change));
1445:   return(0);
1446: }

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

1453:    Collective on PC

1455:    Input Parameters:
1456: +  pc - the preconditioner context
1457: -  ksp - the Krylov subspace context

1459:    Level: developer

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

1468:    Notes:
1469:    The pre-solve phase is distinct from the PCSetUp() phase.

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

1473: .keywords: PC, pre-solve

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

1485:   pc->presolvedone++;
1486:   if (pc->presolvedone > 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot embed PCPreSolve() more than twice");
1487:   KSPGetSolution(ksp,&x);
1488:   KSPGetRhs(ksp,&rhs);

1490:   if (pc->ops->presolve) {
1491:     (*pc->ops->presolve)(pc,ksp,rhs,x);
1492:   }
1493:   return(0);
1494: }

1496: /*@
1497:    PCPostSolve - Optional post-solve phase, intended for any
1498:    preconditioner-specific actions that must be performed after
1499:    the iterative solve itself.

1501:    Collective on PC

1503:    Input Parameters:
1504: +  pc - the preconditioner context
1505: -  ksp - the Krylov subspace context

1507:    Sample of Usage:
1508: .vb
1509:     PCPreSolve(pc,ksp);
1510:     KSPSolve(ksp,b,x);
1511:     PCPostSolve(pc,ksp);
1512: .ve

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

1517:    Level: developer

1519: .keywords: PC, post-solve

1521: .seealso: PCPreSolve(), KSPSolve()
1522: @*/
1523: PetscErrorCode  PCPostSolve(PC pc,KSP ksp)
1524: {
1526:   Vec            x,rhs;

1531:   pc->presolvedone--;
1532:   KSPGetSolution(ksp,&x);
1533:   KSPGetRhs(ksp,&rhs);
1534:   if (pc->ops->postsolve) {
1535:      (*pc->ops->postsolve)(pc,ksp,rhs,x);
1536:   }
1537:   return(0);
1538: }

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

1543:   Collective on PetscViewer

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

1550:    Level: intermediate

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

1555:   Notes for advanced users:
1556:   Most users should not need to know the details of the binary storage
1557:   format, since PCLoad() and PCView() completely hide these details.
1558:   But for anyone who's interested, the standard binary matrix storage
1559:   format is
1560: .vb
1561:      has not yet been determined
1562: .ve

1564: .seealso: PetscViewerBinaryOpen(), PCView(), MatLoad(), VecLoad()
1565: @*/
1566: PetscErrorCode  PCLoad(PC newdm, PetscViewer viewer)
1567: {
1569:   PetscBool      isbinary;
1570:   PetscInt       classid;
1571:   char           type[256];

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

1579:   PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
1580:   if (classid != PC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not PC next in file");
1581:   PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
1582:   PCSetType(newdm, type);
1583:   if (newdm->ops->load) {
1584:     (*newdm->ops->load)(newdm,viewer);
1585:   }
1586:   return(0);
1587: }

1589:  #include <petscdraw.h>
1590: #if defined(PETSC_HAVE_SAWS)
1591:  #include <petscviewersaws.h>
1592: #endif
1593: /*@C
1594:    PCView - Prints the PC data structure.

1596:    Collective on PC

1598:    Input Parameters:
1599: +  PC - the PC context
1600: -  viewer - optional visualization context

1602:    Note:
1603:    The available visualization contexts include
1604: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
1605: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1606:          output where only the first processor opens
1607:          the file.  All other processors send their
1608:          data to the first processor to print.

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

1613:    Level: developer

1615: .keywords: PC, view

1617: .seealso: KSPView(), PetscViewerASCIIOpen()
1618: @*/
1619: PetscErrorCode  PCView(PC pc,PetscViewer viewer)
1620: {
1621:   PCType         cstr;
1623:   PetscBool      iascii,isstring,isbinary,isdraw;
1624: #if defined(PETSC_HAVE_SAWS)
1625:   PetscBool      issaws;
1626: #endif

1630:   if (!viewer) {
1631:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);
1632:   }

1636:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1637:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
1638:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1639:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1640: #if defined(PETSC_HAVE_SAWS)
1641:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws);
1642: #endif

1644:   if (iascii) {
1645:     PetscObjectPrintClassNamePrefixType((PetscObject)pc,viewer);
1646:     if (!pc->setupcalled) {
1647:       PetscViewerASCIIPrintf(viewer,"  PC has not been set up so information may be incomplete\n");
1648:     }
1649:     if (pc->ops->view) {
1650:       PetscViewerASCIIPushTab(viewer);
1651:       (*pc->ops->view)(pc,viewer);
1652:       PetscViewerASCIIPopTab(viewer);
1653:     }
1654:     if (pc->mat) {
1655:       PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
1656:       if (pc->pmat == pc->mat) {
1657:         PetscViewerASCIIPrintf(viewer,"  linear system matrix = precond matrix:\n");
1658:         PetscViewerASCIIPushTab(viewer);
1659:         MatView(pc->mat,viewer);
1660:         PetscViewerASCIIPopTab(viewer);
1661:       } else {
1662:         if (pc->pmat) {
1663:           PetscViewerASCIIPrintf(viewer,"  linear system matrix followed by preconditioner matrix:\n");
1664:         } else {
1665:           PetscViewerASCIIPrintf(viewer,"  linear system matrix:\n");
1666:         }
1667:         PetscViewerASCIIPushTab(viewer);
1668:         MatView(pc->mat,viewer);
1669:         if (pc->pmat) {MatView(pc->pmat,viewer);}
1670:         PetscViewerASCIIPopTab(viewer);
1671:       }
1672:       PetscViewerPopFormat(viewer);
1673:     }
1674:   } else if (isstring) {
1675:     PCGetType(pc,&cstr);
1676:     PetscViewerStringSPrintf(viewer," %-7.7s",cstr);
1677:     if (pc->ops->view) {(*pc->ops->view)(pc,viewer);}
1678:   } else if (isbinary) {
1679:     PetscInt    classid = PC_FILE_CLASSID;
1680:     MPI_Comm    comm;
1681:     PetscMPIInt rank;
1682:     char        type[256];

1684:     PetscObjectGetComm((PetscObject)pc,&comm);
1685:     MPI_Comm_rank(comm,&rank);
1686:     if (!rank) {
1687:       PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
1688:       PetscStrncpy(type,((PetscObject)pc)->type_name,256);
1689:       PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);
1690:     }
1691:     if (pc->ops->view) {
1692:       (*pc->ops->view)(pc,viewer);
1693:     }
1694:   } else if (isdraw) {
1695:     PetscDraw draw;
1696:     char      str[25];
1697:     PetscReal x,y,bottom,h;
1698:     PetscInt  n;

1700:     PetscViewerDrawGetDraw(viewer,0,&draw);
1701:     PetscDrawGetCurrentPoint(draw,&x,&y);
1702:     if (pc->mat) {
1703:       MatGetSize(pc->mat,&n,NULL);
1704:       PetscSNPrintf(str,25,"PC: %s (%D)",((PetscObject)pc)->type_name,n);
1705:     } else {
1706:       PetscSNPrintf(str,25,"PC: %s",((PetscObject)pc)->type_name);
1707:     }
1708:     PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
1709:     bottom = y - h;
1710:     PetscDrawPushCurrentPoint(draw,x,bottom);
1711:     if (pc->ops->view) {
1712:       (*pc->ops->view)(pc,viewer);
1713:     }
1714:     PetscDrawPopCurrentPoint(draw);
1715: #if defined(PETSC_HAVE_SAWS)
1716:   } else if (issaws) {
1717:     PetscMPIInt rank;

1719:     PetscObjectName((PetscObject)pc);
1720:     MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
1721:     if (!((PetscObject)pc)->amsmem && !rank) {
1722:       PetscObjectViewSAWs((PetscObject)pc,viewer);
1723:     }
1724:     if (pc->mat) {MatView(pc->mat,viewer);}
1725:     if (pc->pmat && pc->pmat != pc->mat) {MatView(pc->pmat,viewer);}
1726: #endif
1727:   }
1728:   return(0);
1729: }

1731: /*@C
1732:   PCRegister -  Adds a method to the preconditioner package.

1734:    Not collective

1736:    Input Parameters:
1737: +  name_solver - name of a new user-defined solver
1738: -  routine_create - routine to create method context

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

1743:    Sample usage:
1744: .vb
1745:    PCRegister("my_solver", MySolverCreate);
1746: .ve

1748:    Then, your solver can be chosen with the procedural interface via
1749: $     PCSetType(pc,"my_solver")
1750:    or at runtime via the option
1751: $     -pc_type my_solver

1753:    Level: advanced

1755: .keywords: PC, register

1757: .seealso: PCRegisterAll(), PCRegisterDestroy()
1758: @*/
1759: PetscErrorCode  PCRegister(const char sname[],PetscErrorCode (*function)(PC))
1760: {

1764:   PCInitializePackage();
1765:   PetscFunctionListAdd(&PCList,sname,function);
1766:   return(0);
1767: }

1769: /*@
1770:     PCComputeExplicitOperator - Computes the explicit preconditioned operator.

1772:     Collective on PC

1774:     Input Parameter:
1775: .   pc - the preconditioner object

1777:     Output Parameter:
1778: .   mat - the explict preconditioned operator

1780:     Notes:
1781:     This computation is done by applying the operators to columns of the
1782:     identity matrix.

1784:     Currently, this routine uses a dense matrix format when 1 processor
1785:     is used and a sparse format otherwise.  This routine is costly in general,
1786:     and is recommended for use only with relatively small systems.

1788:     Level: advanced

1790: .keywords: PC, compute, explicit, operator

1792: .seealso: KSPComputeExplicitOperator()

1794: @*/
1795: PetscErrorCode  PCComputeExplicitOperator(PC pc,Mat *mat)
1796: {
1797:   Vec            in,out;
1799:   PetscInt       i,M,m,*rows,start,end;
1800:   PetscMPIInt    size;
1801:   MPI_Comm       comm;
1802:   PetscScalar    *array,one = 1.0;


1808:   PetscObjectGetComm((PetscObject)pc,&comm);
1809:   MPI_Comm_size(comm,&size);

1811:   if (!pc->pmat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"You must call KSPSetOperators() or PCSetOperators() before this call");
1812:   MatCreateVecs(pc->pmat,&in,0);
1813:   VecDuplicate(in,&out);
1814:   VecGetOwnershipRange(in,&start,&end);
1815:   VecGetSize(in,&M);
1816:   VecGetLocalSize(in,&m);
1817:   PetscMalloc1(m+1,&rows);
1818:   for (i=0; i<m; i++) rows[i] = start + i;

1820:   MatCreate(comm,mat);
1821:   MatSetSizes(*mat,m,m,M,M);
1822:   if (size == 1) {
1823:     MatSetType(*mat,MATSEQDENSE);
1824:     MatSeqDenseSetPreallocation(*mat,NULL);
1825:   } else {
1826:     MatSetType(*mat,MATMPIAIJ);
1827:     MatMPIAIJSetPreallocation(*mat,0,NULL,0,NULL);
1828:   }
1829:   MatSetOption(*mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);

1831:   for (i=0; i<M; i++) {

1833:     VecSet(in,0.0);
1834:     VecSetValues(in,1,&i,&one,INSERT_VALUES);
1835:     VecAssemblyBegin(in);
1836:     VecAssemblyEnd(in);

1838:     /* should fix, allowing user to choose side */
1839:     PCApply(pc,in,out);

1841:     VecGetArray(out,&array);
1842:     MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
1843:     VecRestoreArray(out,&array);

1845:   }
1846:   PetscFree(rows);
1847:   VecDestroy(&out);
1848:   VecDestroy(&in);
1849:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
1850:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
1851:   return(0);
1852: }

1854: /*@
1855:    PCSetCoordinates - sets the coordinates of all the nodes on the local process

1857:    Collective on PC

1859:    Input Parameters:
1860: +  pc - the solver context
1861: .  dim - the dimension of the coordinates 1, 2, or 3
1862: .  nloc - the blocked size of the coordinates array
1863: -  coords - the coordinates array

1865:    Level: intermediate

1867:    Notes:
1868:    coords is an array of the dim coordinates for the nodes on
1869:    the local processor, of size dim*nloc.
1870:    If there are 108 equation on a processor
1871:    for a displacement finite element discretization of elasticity (so
1872:    that there are nloc = 36 = 108/3 nodes) then the array must have 108
1873:    double precision values (ie, 3 * 36).  These x y z coordinates
1874:    should be ordered for nodes 0 to N-1 like so: [ 0.x, 0.y, 0.z, 1.x,
1875:    ... , N-1.z ].

1877: .seealso: MatSetNearNullSpace()
1878: @*/
1879: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
1880: {

1886:   PetscTryMethod(pc,"PCSetCoordinates_C",(PC,PetscInt,PetscInt,PetscReal*),(pc,dim,nloc,coords));
1887:   return(0);
1888: }