Actual source code: precon.c

petsc-3.8.3 2017-12-09
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, PC_ApplyOnMproc;
 12: PetscInt      PetscMGLevelId;

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

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

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

 58:    Collective on PC

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

 63:    Level: developer

 65:    Notes: 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: If this returns PETSC_TRUE then the system solved via the Krylov method is
137: $           D M A D^{-1} y = D M b  for left preconditioning or
138: $           D A M D^{-1} z = D b for right preconditioning

140: .keywords: PC

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

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

157:    Logically Collective on PC

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

163:    Level: intermediate

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

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

171: .keywords: PC

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

182:   pc->diagonalscale     = PETSC_TRUE;

184:   PetscObjectReference((PetscObject)s);
185:   VecDestroy(&pc->diagonalscaleleft);

187:   pc->diagonalscaleleft = s;

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

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

198:    Logically Collective on PC

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

205:    Level: intermediate

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

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

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

215: .keywords: PC

217: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleSet(), PCDiagonalScaleRight(), PCDiagonalScale()
218: @*/
219: PetscErrorCode  PCDiagonalScaleLeft(PC pc,Vec in,Vec out)
220: {

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

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

238:    Logically Collective on PC

240:    Input Parameters:
241: +  pc - the preconditioner context
242: .  in - input vector
243: +  out - scaled vector (maybe the same as in)

245:    Level: intermediate

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

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

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

255: .keywords: PC

257: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleSet(), PCDiagonalScale()
258: @*/
259: PetscErrorCode  PCDiagonalScaleRight(PC pc,Vec in,Vec out)
260: {

267:   if (pc->diagonalscale) {
268:     VecPointwiseMult(out,pc->diagonalscaleright,in);
269:   } else if (in != out) {
270:     VecCopy(in,out);
271:   }
272:   return(0);
273: }

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

280:    Logically Collective on PC

282:    Input Parameters:
283: +  pc - the preconditioner context
284: -  flg - PETSC_TRUE to use the Amat, PETSC_FALSE to use the Pmat (default is false)

286:    Options Database Key:
287: .  -pc_use_amat <true,false>

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

293:    Level: intermediate

295: .seealso: PCGetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
296: @*/
297: PetscErrorCode  PCSetUseAmat(PC pc,PetscBool flg)
298: {
301:   pc->useAmat = flg;
302:   return(0);
303: }

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

308:    Logically Collective on PC

310:    Input Parameters:
311: +  pc - iterative context obtained from PCCreate()
312: -  flg - PETSC_TRUE indicates you want the error generated

314:    Level: advanced

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

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

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

325: .seealso: PCGetInitialGuessNonzero(), PCSetInitialGuessKnoll(), PCGetInitialGuessKnoll()
326: @*/
327: PetscErrorCode  PCSetErrorIfFailure(PC pc,PetscBool flg)
328: {
332:   pc->erroriffailure = flg;
333:   return(0);
334: }

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

341:    Logically Collective on PC

343:    Input Parameter:
344: .  pc - the preconditioner context

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

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

353:    Level: intermediate

355: .seealso: PCSetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
356: @*/
357: PetscErrorCode  PCGetUseAmat(PC pc,PetscBool *flg)
358: {
361:   *flg = pc->useAmat;
362:   return(0);
363: }

365: /*@
366:    PCCreate - Creates a preconditioner context.

368:    Collective on MPI_Comm

370:    Input Parameter:
371: .  comm - MPI communicator

373:    Output Parameter:
374: .  pc - location to put the preconditioner context

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

380:    Level: developer

382: .keywords: PC, create, context

384: .seealso: PCSetUp(), PCApply(), PCDestroy()
385: @*/
386: PetscErrorCode  PCCreate(MPI_Comm comm,PC *newpc)
387: {
388:   PC             pc;

393:   *newpc = 0;
394:   PCInitializePackage();

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

398:   pc->mat                  = 0;
399:   pc->pmat                 = 0;
400:   pc->setupcalled          = 0;
401:   pc->setfromoptionscalled = 0;
402:   pc->data                 = 0;
403:   pc->diagonalscale        = PETSC_FALSE;
404:   pc->diagonalscaleleft    = 0;
405:   pc->diagonalscaleright   = 0;

407:   pc->modifysubmatrices  = 0;
408:   pc->modifysubmatricesP = 0;

410:   *newpc = pc;
411:   return(0);

413: }

415: /* -------------------------------------------------------------------------------*/

417: /*@
418:    PCApply - Applies the preconditioner to a vector.

420:    Collective on PC and Vec

422:    Input Parameters:
423: +  pc - the preconditioner context
424: -  x - input vector

426:    Output Parameter:
427: .  y - output vector

429:    Level: developer

431: .keywords: PC, apply

433: .seealso: PCApplyTranspose(), PCApplyBAorAB()
434: @*/
435: PetscErrorCode  PCApply(PC pc,Vec x,Vec y)
436: {
438:   PetscInt       m,n,mv,nv;

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

454:   PCSetUp(pc);
455:   if (!pc->ops->apply) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply");
456:   VecLockPush(x);
457:   PetscLogEventBegin(PC_Apply,pc,x,y,0);
458:   (*pc->ops->apply)(pc,x,y);
459:   PetscLogEventEnd(PC_Apply,pc,x,y,0);
460:   if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
461:   VecLockPop(x);
462:   return(0);
463: }

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

468:    Collective on PC and Vec

470:    Input Parameters:
471: +  pc - the preconditioner context
472: -  x - input vector

474:    Output Parameter:
475: .  y - output vector

477:    Notes:
478:    Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.

480:    Level: developer

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

484: .seealso: PCApply(), PCApplySymmetricRight()
485: @*/
486: PetscErrorCode  PCApplySymmetricLeft(PC pc,Vec x,Vec y)
487: {

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

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

510:    Collective on PC and Vec

512:    Input Parameters:
513: +  pc - the preconditioner context
514: -  x - input vector

516:    Output Parameter:
517: .  y - output vector

519:    Level: developer

521:    Notes:
522:    Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.

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

526: .seealso: PCApply(), PCApplySymmetricLeft()
527: @*/
528: PetscErrorCode  PCApplySymmetricRight(PC pc,Vec x,Vec y)
529: {

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

549: /*@
550:    PCApplyTranspose - Applies the transpose of preconditioner to a vector.

552:    Collective on PC and Vec

554:    Input Parameters:
555: +  pc - the preconditioner context
556: -  x - input vector

558:    Output Parameter:
559: .  y - output vector

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

563:    Developer Notes: We need to implement a PCApplyHermitianTranspose()

565:    Level: developer

567: .keywords: PC, apply, transpose

569: .seealso: PCApply(), PCApplyBAorAB(), PCApplyBAorABTranspose(), PCApplyTransposeExists()
570: @*/
571: PetscErrorCode  PCApplyTranspose(PC pc,Vec x,Vec y)
572: {

579:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
580:   if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
581:   PCSetUp(pc);
582:   if (!pc->ops->applytranspose) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply transpose");
583:   VecLockPush(x);
584:   PetscLogEventBegin(PC_Apply,pc,x,y,0);
585:   (*pc->ops->applytranspose)(pc,x,y);
586:   PetscLogEventEnd(PC_Apply,pc,x,y,0);
587:   VecLockPop(x);
588:   if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
589:   return(0);
590: }

592: /*@
593:    PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation

595:    Collective on PC and Vec

597:    Input Parameters:
598: .  pc - the preconditioner context

600:    Output Parameter:
601: .  flg - PETSC_TRUE if a transpose operation is defined

603:    Level: developer

605: .keywords: PC, apply, transpose

607: .seealso: PCApplyTranspose()
608: @*/
609: PetscErrorCode  PCApplyTransposeExists(PC pc,PetscBool  *flg)
610: {
614:   if (pc->ops->applytranspose) *flg = PETSC_TRUE;
615:   else *flg = PETSC_FALSE;
616:   return(0);
617: }

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

622:    Collective on PC and Vec

624:    Input Parameters:
625: +  pc - the preconditioner context
626: .  side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
627: .  x - input vector
628: -  work - work vector

630:    Output Parameter:
631: .  y - output vector

633:    Level: developer

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

638: .keywords: PC, apply, operator

640: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorABTranspose()
641: @*/
642: PetscErrorCode  PCApplyBAorAB(PC pc,PCSide side,Vec x,Vec y,Vec work)
643: {

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

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

697: /*@
698:    PCApplyBAorABTranspose - Applies the transpose of the preconditioner
699:    and operator to a vector. That is, applies tr(B) * tr(A) with left preconditioning,
700:    NOT tr(B*A) = tr(A)*tr(B).

702:    Collective on PC and Vec

704:    Input Parameters:
705: +  pc - the preconditioner context
706: .  side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
707: .  x - input vector
708: -  work - work vector

710:    Output Parameter:
711: .  y - output vector


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

717:     Level: developer

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

721: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorAB()
722: @*/
723: PetscErrorCode  PCApplyBAorABTranspose(PC pc,PCSide side,Vec x,Vec y,Vec work)
724: {

732:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
733:   if (pc->erroriffailure) {VecValidValues(x,3,PETSC_TRUE);}
734:   if (pc->ops->applyBAtranspose) {
735:     (*pc->ops->applyBAtranspose)(pc,side,x,y,work);
736:     if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
737:     return(0);
738:   }
739:   if (side != PC_LEFT && side != PC_RIGHT) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Side must be right or left");

741:   PCSetUp(pc);
742:   if (side == PC_RIGHT) {
743:     PCApplyTranspose(pc,x,work);
744:     MatMultTranspose(pc->mat,work,y);
745:   } else if (side == PC_LEFT) {
746:     MatMultTranspose(pc->mat,x,work);
747:     PCApplyTranspose(pc,work,y);
748:   }
749:   /* add support for PC_SYMMETRIC */
750:   if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
751:   return(0);
752: }

754: /* -------------------------------------------------------------------------------*/

756: /*@
757:    PCApplyRichardsonExists - Determines whether a particular preconditioner has a
758:    built-in fast application of Richardson's method.

760:    Not Collective

762:    Input Parameter:
763: .  pc - the preconditioner

765:    Output Parameter:
766: .  exists - PETSC_TRUE or PETSC_FALSE

768:    Level: developer

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

772: .seealso: PCApplyRichardson()
773: @*/
774: PetscErrorCode  PCApplyRichardsonExists(PC pc,PetscBool  *exists)
775: {
779:   if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
780:   else *exists = PETSC_FALSE;
781:   return(0);
782: }

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

789:    Collective on PC

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

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

806:    Notes:
807:    Most preconditioners do not support this function. Use the command
808:    PCApplyRichardsonExists() to determine if one does.

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

813:    Level: developer

815: .keywords: PC, apply, Richardson

817: .seealso: PCApplyRichardsonExists()
818: @*/
819: PetscErrorCode  PCApplyRichardson(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
820: {

828:   if (b == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"b and y must be different vectors");
829:   PCSetUp(pc);
830:   if (!pc->ops->applyrichardson) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply richardson");
831:   (*pc->ops->applyrichardson)(pc,b,y,w,rtol,abstol,dtol,its,guesszero,outits,reason);
832:   return(0);
833: }

835: /*@
836:    PCGetSetUpFailedReason - Gets the reason a PCSetUp() failed or 0 if it did not fail

838:    Logically Collective on PC

840:    Input Parameter:
841: .  pc - the preconditioner context

843:    Output Parameter:
844: .  reason - the reason it failed, currently only -1 

846:    Level: advanced

848: .keywords: PC, setup

850: .seealso: PCCreate(), PCApply(), PCDestroy()
851: @*/
852: PetscErrorCode PCGetSetUpFailedReason(PC pc,PCFailedReason *reason)
853: {
855:   if (pc->setupcalled < 0) *reason = (PCFailedReason)pc->setupcalled;
856:   else *reason = pc->failedreason;
857:   return(0);
858: }


861: /*
862:       a setupcall of 0 indicates never setup,
863:                      1 indicates has been previously setup
864:                     -1 indicates a PCSetUp() was attempted and failed
865: */
866: /*@
867:    PCSetUp - Prepares for the use of a preconditioner.

869:    Collective on PC

871:    Input Parameter:
872: .  pc - the preconditioner context

874:    Level: developer

876: .keywords: PC, setup

878: .seealso: PCCreate(), PCApply(), PCDestroy()
879: @*/
880: PetscErrorCode  PCSetUp(PC pc)
881: {
882:   PetscErrorCode   ierr;
883:   const char       *def;
884:   PetscObjectState matstate, matnonzerostate;

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

890:   if (pc->setupcalled && pc->reusepreconditioner) {
891:     PetscInfo(pc,"Leaving PC with identical preconditioner since reuse preconditioner is set\n");
892:     return(0);
893:   }

895:   PetscObjectStateGet((PetscObject)pc->pmat,&matstate);
896:   MatGetNonzeroState(pc->pmat,&matnonzerostate);
897:   if (!pc->setupcalled) {
898:     PetscInfo(pc,"Setting up PC for first time\n");
899:     pc->flag        = DIFFERENT_NONZERO_PATTERN;
900:   } else if (matstate == pc->matstate) {
901:     PetscInfo(pc,"Leaving PC with identical preconditioner since operator is unchanged\n");
902:     return(0);
903:   } else {
904:     if (matnonzerostate > pc->matnonzerostate) {
905:        PetscInfo(pc,"Setting up PC with different nonzero pattern\n");
906:        pc->flag            = DIFFERENT_NONZERO_PATTERN;
907:     } else {
908:       PetscInfo(pc,"Setting up PC with same nonzero pattern\n");
909:       pc->flag            = SAME_NONZERO_PATTERN;
910:     }
911:   }
912:   pc->matstate        = matstate;
913:   pc->matnonzerostate = matnonzerostate;

915:   if (!((PetscObject)pc)->type_name) {
916:     PCGetDefaultType_Private(pc,&def);
917:     PCSetType(pc,def);
918:   }

920:   MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
921:   MatSetErrorIfFailure(pc->mat,pc->erroriffailure);
922:   PetscLogEventBegin(PC_SetUp,pc,0,0,0);
923:   if (pc->ops->setup) {
924:     (*pc->ops->setup)(pc);
925:   }
926:   PetscLogEventEnd(PC_SetUp,pc,0,0,0);
927:   if (!pc->setupcalled) pc->setupcalled = 1;
928:   return(0);
929: }

931: /*@
932:    PCSetUpOnBlocks - Sets up the preconditioner for each block in
933:    the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
934:    methods.

936:    Collective on PC

938:    Input Parameters:
939: .  pc - the preconditioner context

941:    Level: developer

943: .keywords: PC, setup, blocks

945: .seealso: PCCreate(), PCApply(), PCDestroy(), PCSetUp()
946: @*/
947: PetscErrorCode  PCSetUpOnBlocks(PC pc)
948: {

953:   if (!pc->ops->setuponblocks) return(0);
954:   PetscLogEventBegin(PC_SetUpOnBlocks,pc,0,0,0);
955:   (*pc->ops->setuponblocks)(pc);
956:   PetscLogEventEnd(PC_SetUpOnBlocks,pc,0,0,0);
957:   return(0);
958: }

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

967:    Logically Collective on PC

969:    Input Parameters:
970: +  pc - the preconditioner context
971: .  func - routine for modifying the submatrices
972: -  ctx - optional user-defined context (may be null)

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

977: .  row - an array of index sets that contain the global row numbers
978:          that comprise each local submatrix
979: .  col - an array of index sets that contain the global column numbers
980:          that comprise each local submatrix
981: .  submat - array of local submatrices
982: -  ctx - optional user-defined context for private data for the
983:          user-defined func routine (may be null)

985:    Notes:
986:    PCSetModifySubMatrices() MUST be called before KSPSetUp() and
987:    KSPSolve().

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

993:    Level: advanced

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

997: .seealso: PCModifySubMatrices(), PCASMGetSubMatrices()
998: @*/
999: PetscErrorCode  PCSetModifySubMatrices(PC pc,PetscErrorCode (*func)(PC,PetscInt,const IS[],const IS[],Mat[],void*),void *ctx)
1000: {
1003:   pc->modifysubmatrices  = func;
1004:   pc->modifysubmatricesP = ctx;
1005:   return(0);
1006: }

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

1012:    Collective on PC

1014:    Input Parameters:
1015: +  pc - the preconditioner context
1016: .  nsub - the number of local submatrices
1017: .  row - an array of index sets that contain the global row numbers
1018:          that comprise each local submatrix
1019: .  col - an array of index sets that contain the global column numbers
1020:          that comprise each local submatrix
1021: .  submat - array of local submatrices
1022: -  ctx - optional user-defined context for private data for the
1023:          user-defined routine (may be null)

1025:    Output Parameter:
1026: .  submat - array of local submatrices (the entries of which may
1027:             have been modified)

1029:    Notes:
1030:    The user should NOT generally call this routine, as it will
1031:    automatically be called within certain preconditioners (currently
1032:    block Jacobi, additive Schwarz) if set.

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

1039:    Level: developer

1041: .keywords: PC, modify, submatrices

1043: .seealso: PCSetModifySubMatrices()
1044: @*/
1045: PetscErrorCode  PCModifySubMatrices(PC pc,PetscInt nsub,const IS row[],const IS col[],Mat submat[],void *ctx)
1046: {

1051:   if (!pc->modifysubmatrices) return(0);
1052:   PetscLogEventBegin(PC_ModifySubMatrices,pc,0,0,0);
1053:   (*pc->modifysubmatrices)(pc,nsub,row,col,submat,ctx);
1054:   PetscLogEventEnd(PC_ModifySubMatrices,pc,0,0,0);
1055:   return(0);
1056: }

1058: /*@
1059:    PCSetOperators - Sets the matrix associated with the linear system and
1060:    a (possibly) different one associated with the preconditioner.

1062:    Logically Collective on PC and Mat

1064:    Input Parameters:
1065: +  pc - the preconditioner context
1066: .  Amat - the matrix that defines the linear system
1067: -  Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.

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

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

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

1082:    Level: intermediate

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

1086: .seealso: PCGetOperators(), MatZeroEntries()
1087:  @*/
1088: PetscErrorCode  PCSetOperators(PC pc,Mat Amat,Mat Pmat)
1089: {
1090:   PetscErrorCode   ierr;
1091:   PetscInt         m1,n1,m2,n2;

1099:   if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
1100:     MatGetLocalSize(Amat,&m1,&n1);
1101:     MatGetLocalSize(pc->mat,&m2,&n2);
1102:     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);
1103:     MatGetLocalSize(Pmat,&m1,&n1);
1104:     MatGetLocalSize(pc->pmat,&m2,&n2);
1105:     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);
1106:   }

1108:   if (Pmat != pc->pmat) {
1109:     /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
1110:     pc->matnonzerostate = -1;
1111:     pc->matstate        = -1;
1112:   }

1114:   /* reference first in case the matrices are the same */
1115:   if (Amat) {PetscObjectReference((PetscObject)Amat);}
1116:   MatDestroy(&pc->mat);
1117:   if (Pmat) {PetscObjectReference((PetscObject)Pmat);}
1118:   MatDestroy(&pc->pmat);
1119:   pc->mat  = Amat;
1120:   pc->pmat = Pmat;
1121:   return(0);
1122: }

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

1127:    Logically Collective on PC

1129:    Input Parameters:
1130: +  pc - the preconditioner context
1131: -  flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner

1133:     Level: intermediate

1135: .seealso: PCGetOperators(), MatZeroEntries(), PCGetReusePreconditioner(), KSPSetReusePreconditioner()
1136:  @*/
1137: PetscErrorCode  PCSetReusePreconditioner(PC pc,PetscBool flag)
1138: {
1141:   pc->reusepreconditioner = flag;
1142:   return(0);
1143: }

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

1148:    Not Collective

1150:    Input Parameter:
1151: .  pc - the preconditioner context

1153:    Output Parameter:
1154: .  flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner

1156:    Level: intermediate

1158: .seealso: PCGetOperators(), MatZeroEntries(), PCSetReusePreconditioner()
1159:  @*/
1160: PetscErrorCode  PCGetReusePreconditioner(PC pc,PetscBool *flag)
1161: {
1164:   *flag = pc->reusepreconditioner;
1165:   return(0);
1166: }

1168: /*@C
1169:    PCGetOperators - Gets the matrix associated with the linear system and
1170:    possibly a different one associated with the preconditioner.

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

1174:    Input Parameter:
1175: .  pc - the preconditioner context

1177:    Output Parameters:
1178: +  Amat - the matrix defining the linear system
1179: -  Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.

1181:    Level: intermediate

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

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

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

1196: $         MatCreate(comm,&mat);
1197: $         KSP/PCSetOperators(ksp/pc,Amat,Amat);
1198: $         PetscObjectDereference((PetscObject)mat);
1199: $           set size, type, etc of Amat

1201:      and

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

1206: $         MatCreate(comm,&Amat);
1207: $         MatCreate(comm,&Pmat);
1208: $         KSP/PCSetOperators(ksp/pc,Amat,Pmat);
1209: $         PetscObjectDereference((PetscObject)Amat);
1210: $         PetscObjectDereference((PetscObject)Pmat);
1211: $           set size, type, etc of Amat and Pmat

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


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

1225: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperatorsSet()
1226: @*/
1227: PetscErrorCode  PCGetOperators(PC pc,Mat *Amat,Mat *Pmat)
1228: {

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

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

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

1272:    Input Parameter:
1273: .  pc - the preconditioner context

1275:    Output Parameters:
1276: +  mat - the matrix associated with the linear system was set
1277: -  pmat - matrix associated with the preconditioner was set, usually the same

1279:    Level: intermediate

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

1283: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperators()
1284: @*/
1285: PetscErrorCode  PCGetOperatorsSet(PC pc,PetscBool  *mat,PetscBool  *pmat)
1286: {
1289:   if (mat) *mat = (pc->mat)  ? PETSC_TRUE : PETSC_FALSE;
1290:   if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1291:   return(0);
1292: }

1294: /*@
1295:    PCFactorGetMatrix - Gets the factored matrix from the
1296:    preconditioner context.  This routine is valid only for the LU,
1297:    incomplete LU, Cholesky, and incomplete Cholesky methods.

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

1301:    Input Parameters:
1302: .  pc - the preconditioner context

1304:    Output parameters:
1305: .  mat - the factored matrix

1307:    Level: advanced

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

1311: .keywords: PC, get, factored, matrix
1312: @*/
1313: PetscErrorCode  PCFactorGetMatrix(PC pc,Mat *mat)
1314: {

1320:   if (pc->ops->getfactoredmatrix) {
1321:     (*pc->ops->getfactoredmatrix)(pc,mat);
1322:   } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC type does not support getting factor matrix");
1323:   return(0);
1324: }

1326: /*@C
1327:    PCSetOptionsPrefix - Sets the prefix used for searching for all
1328:    PC options in the database.

1330:    Logically Collective on PC

1332:    Input Parameters:
1333: +  pc - the preconditioner context
1334: -  prefix - the prefix string to prepend to all PC option requests

1336:    Notes:
1337:    A hyphen (-) must NOT be given at the beginning of the prefix name.
1338:    The first character of all runtime options is AUTOMATICALLY the
1339:    hyphen.

1341:    Level: advanced

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

1345: .seealso: PCAppendOptionsPrefix(), PCGetOptionsPrefix()
1346: @*/
1347: PetscErrorCode  PCSetOptionsPrefix(PC pc,const char prefix[])
1348: {

1353:   PetscObjectSetOptionsPrefix((PetscObject)pc,prefix);
1354:   return(0);
1355: }

1357: /*@C
1358:    PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1359:    PC options in the database.

1361:    Logically Collective on PC

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

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

1372:    Level: advanced

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

1376: .seealso: PCSetOptionsPrefix(), PCGetOptionsPrefix()
1377: @*/
1378: PetscErrorCode  PCAppendOptionsPrefix(PC pc,const char prefix[])
1379: {

1384:   PetscObjectAppendOptionsPrefix((PetscObject)pc,prefix);
1385:   return(0);
1386: }

1388: /*@C
1389:    PCGetOptionsPrefix - Gets the prefix used for searching for all
1390:    PC options in the database.

1392:    Not Collective

1394:    Input Parameters:
1395: .  pc - the preconditioner context

1397:    Output Parameters:
1398: .  prefix - pointer to the prefix string used, is returned

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

1403:    Level: advanced

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

1407: .seealso: PCSetOptionsPrefix(), PCAppendOptionsPrefix()
1408: @*/
1409: PetscErrorCode  PCGetOptionsPrefix(PC pc,const char *prefix[])
1410: {

1416:   PetscObjectGetOptionsPrefix((PetscObject)pc,prefix);
1417:   return(0);
1418: }

1420: /*
1421:    Indicates the right hand side will be changed by KSPSolve(), this occurs for a few
1422:   preconditioners including BDDC and Eisentat that transform the equations before applying
1423:   the Krylov methods
1424: */
1425: PETSC_INTERN PetscErrorCode  PCPreSolveChangeRHS(PC pc,PetscBool *change)
1426: {

1432:   *change = PETSC_FALSE;
1433:   PetscTryMethod(pc,"PCPreSolveChangeRHS_C",(PC,PetscBool*),(pc,change));
1434:   return(0);
1435: }

1437: /*@
1438:    PCPreSolve - Optional pre-solve phase, intended for any
1439:    preconditioner-specific actions that must be performed before
1440:    the iterative solve itself.

1442:    Collective on PC

1444:    Input Parameters:
1445: +  pc - the preconditioner context
1446: -  ksp - the Krylov subspace context

1448:    Level: developer

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

1457:    Notes:
1458:    The pre-solve phase is distinct from the PCSetUp() phase.

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

1462: .keywords: PC, pre-solve

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

1474:   pc->presolvedone++;
1475:   if (pc->presolvedone > 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot embed PCPreSolve() more than twice");
1476:   KSPGetSolution(ksp,&x);
1477:   KSPGetRhs(ksp,&rhs);

1479:   if (pc->ops->presolve) {
1480:     (*pc->ops->presolve)(pc,ksp,rhs,x);
1481:   }
1482:   return(0);
1483: }

1485: /*@
1486:    PCPostSolve - Optional post-solve phase, intended for any
1487:    preconditioner-specific actions that must be performed after
1488:    the iterative solve itself.

1490:    Collective on PC

1492:    Input Parameters:
1493: +  pc - the preconditioner context
1494: -  ksp - the Krylov subspace context

1496:    Sample of Usage:
1497: .vb
1498:     PCPreSolve(pc,ksp);
1499:     KSPSolve(ksp,b,x);
1500:     PCPostSolve(pc,ksp);
1501: .ve

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

1506:    Level: developer

1508: .keywords: PC, post-solve

1510: .seealso: PCPreSolve(), KSPSolve()
1511: @*/
1512: PetscErrorCode  PCPostSolve(PC pc,KSP ksp)
1513: {
1515:   Vec            x,rhs;

1520:   pc->presolvedone--;
1521:   KSPGetSolution(ksp,&x);
1522:   KSPGetRhs(ksp,&rhs);
1523:   if (pc->ops->postsolve) {
1524:      (*pc->ops->postsolve)(pc,ksp,rhs,x);
1525:   }
1526:   return(0);
1527: }

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

1532:   Collective on PetscViewer

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

1539:    Level: intermediate

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

1544:   Notes for advanced users:
1545:   Most users should not need to know the details of the binary storage
1546:   format, since PCLoad() and PCView() completely hide these details.
1547:   But for anyone who's interested, the standard binary matrix storage
1548:   format is
1549: .vb
1550:      has not yet been determined
1551: .ve

1553: .seealso: PetscViewerBinaryOpen(), PCView(), MatLoad(), VecLoad()
1554: @*/
1555: PetscErrorCode  PCLoad(PC newdm, PetscViewer viewer)
1556: {
1558:   PetscBool      isbinary;
1559:   PetscInt       classid;
1560:   char           type[256];

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

1568:   PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
1569:   if (classid != PC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not PC next in file");
1570:   PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
1571:   PCSetType(newdm, type);
1572:   if (newdm->ops->load) {
1573:     (*newdm->ops->load)(newdm,viewer);
1574:   }
1575:   return(0);
1576: }

1578:  #include <petscdraw.h>
1579: #if defined(PETSC_HAVE_SAWS)
1580:  #include <petscviewersaws.h>
1581: #endif
1582: /*@C
1583:    PCView - Prints the PC data structure.

1585:    Collective on PC

1587:    Input Parameters:
1588: +  PC - the PC context
1589: -  viewer - optional visualization context

1591:    Note:
1592:    The available visualization contexts include
1593: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
1594: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1595:          output where only the first processor opens
1596:          the file.  All other processors send their
1597:          data to the first processor to print.

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

1602:    Level: developer

1604: .keywords: PC, view

1606: .seealso: KSPView(), PetscViewerASCIIOpen()
1607: @*/
1608: PetscErrorCode  PCView(PC pc,PetscViewer viewer)
1609: {
1610:   PCType            cstr;
1611:   PetscErrorCode    ierr;
1612:   PetscBool         iascii,isstring,isbinary,isdraw;
1613:   PetscViewerFormat format;
1614: #if defined(PETSC_HAVE_SAWS)
1615:   PetscBool         issaws;
1616: #endif

1620:   if (!viewer) {
1621:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);
1622:   }

1626:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1627:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
1628:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1629:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1630: #if defined(PETSC_HAVE_SAWS)
1631:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws);
1632: #endif

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

1675:     PetscObjectGetComm((PetscObject)pc,&comm);
1676:     MPI_Comm_rank(comm,&rank);
1677:     if (!rank) {
1678:       PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
1679:       PetscStrncpy(type,((PetscObject)pc)->type_name,256);
1680:       PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);
1681:     }
1682:     if (pc->ops->view) {
1683:       (*pc->ops->view)(pc,viewer);
1684:     }
1685:   } else if (isdraw) {
1686:     PetscDraw draw;
1687:     char      str[25];
1688:     PetscReal x,y,bottom,h;
1689:     PetscInt  n;

1691:     PetscViewerDrawGetDraw(viewer,0,&draw);
1692:     PetscDrawGetCurrentPoint(draw,&x,&y);
1693:     if (pc->mat) {
1694:       MatGetSize(pc->mat,&n,NULL);
1695:       PetscSNPrintf(str,25,"PC: %s (%D)",((PetscObject)pc)->type_name,n);
1696:     } else {
1697:       PetscSNPrintf(str,25,"PC: %s",((PetscObject)pc)->type_name);
1698:     }
1699:     PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
1700:     bottom = y - h;
1701:     PetscDrawPushCurrentPoint(draw,x,bottom);
1702:     if (pc->ops->view) {
1703:       (*pc->ops->view)(pc,viewer);
1704:     }
1705:     PetscDrawPopCurrentPoint(draw);
1706: #if defined(PETSC_HAVE_SAWS)
1707:   } else if (issaws) {
1708:     PetscMPIInt rank;

1710:     PetscObjectName((PetscObject)pc);
1711:     MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
1712:     if (!((PetscObject)pc)->amsmem && !rank) {
1713:       PetscObjectViewSAWs((PetscObject)pc,viewer);
1714:     }
1715:     if (pc->mat) {MatView(pc->mat,viewer);}
1716:     if (pc->pmat && pc->pmat != pc->mat) {MatView(pc->pmat,viewer);}
1717: #endif
1718:   }
1719:   return(0);
1720: }

1722: /*@C
1723:   PCRegister -  Adds a method to the preconditioner package.

1725:    Not collective

1727:    Input Parameters:
1728: +  name_solver - name of a new user-defined solver
1729: -  routine_create - routine to create method context

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

1734:    Sample usage:
1735: .vb
1736:    PCRegister("my_solver", MySolverCreate);
1737: .ve

1739:    Then, your solver can be chosen with the procedural interface via
1740: $     PCSetType(pc,"my_solver")
1741:    or at runtime via the option
1742: $     -pc_type my_solver

1744:    Level: advanced

1746: .keywords: PC, register

1748: .seealso: PCRegisterAll(), PCRegisterDestroy()
1749: @*/
1750: PetscErrorCode  PCRegister(const char sname[],PetscErrorCode (*function)(PC))
1751: {

1755:   PetscFunctionListAdd(&PCList,sname,function);
1756:   return(0);
1757: }

1759: /*@
1760:     PCComputeExplicitOperator - Computes the explicit preconditioned operator.

1762:     Collective on PC

1764:     Input Parameter:
1765: .   pc - the preconditioner object

1767:     Output Parameter:
1768: .   mat - the explict preconditioned operator

1770:     Notes:
1771:     This computation is done by applying the operators to columns of the
1772:     identity matrix.

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

1778:     Level: advanced

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

1782: .seealso: KSPComputeExplicitOperator()

1784: @*/
1785: PetscErrorCode  PCComputeExplicitOperator(PC pc,Mat *mat)
1786: {
1787:   Vec            in,out;
1789:   PetscInt       i,M,m,*rows,start,end;
1790:   PetscMPIInt    size;
1791:   MPI_Comm       comm;
1792:   PetscScalar    *array,one = 1.0;


1798:   PetscObjectGetComm((PetscObject)pc,&comm);
1799:   MPI_Comm_size(comm,&size);

1801:   if (!pc->pmat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"You must call KSPSetOperators() or PCSetOperators() before this call");
1802:   MatCreateVecs(pc->pmat,&in,0);
1803:   VecDuplicate(in,&out);
1804:   VecGetOwnershipRange(in,&start,&end);
1805:   VecGetSize(in,&M);
1806:   VecGetLocalSize(in,&m);
1807:   PetscMalloc1(m+1,&rows);
1808:   for (i=0; i<m; i++) rows[i] = start + i;

1810:   MatCreate(comm,mat);
1811:   MatSetSizes(*mat,m,m,M,M);
1812:   if (size == 1) {
1813:     MatSetType(*mat,MATSEQDENSE);
1814:     MatSeqDenseSetPreallocation(*mat,NULL);
1815:   } else {
1816:     MatSetType(*mat,MATMPIAIJ);
1817:     MatMPIAIJSetPreallocation(*mat,0,NULL,0,NULL);
1818:   }
1819:   MatSetOption(*mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);

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

1823:     VecSet(in,0.0);
1824:     VecSetValues(in,1,&i,&one,INSERT_VALUES);
1825:     VecAssemblyBegin(in);
1826:     VecAssemblyEnd(in);

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

1831:     VecGetArray(out,&array);
1832:     MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
1833:     VecRestoreArray(out,&array);

1835:   }
1836:   PetscFree(rows);
1837:   VecDestroy(&out);
1838:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
1839:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
1840:   return(0);
1841: }

1843: /*@
1844:    PCSetCoordinates - sets the coordinates of all the nodes on the local process

1846:    Collective on PC

1848:    Input Parameters:
1849: +  pc - the solver context
1850: .  dim - the dimension of the coordinates 1, 2, or 3
1851: -  coords - the coordinates

1853:    Level: intermediate

1855:    Notes: coords is an array of the 3D coordinates for the nodes on
1856:    the local processor.  So if there are 108 equation on a processor
1857:    for a displacement finite element discretization of elasticity (so
1858:    that there are 36 = 108/3 nodes) then the array must have 108
1859:    double precision values (ie, 3 * 36).  These x y z coordinates
1860:    should be ordered for nodes 0 to N-1 like so: [ 0.x, 0.y, 0.z, 1.x,
1861:    ... , N-1.z ].

1863: .seealso: MatSetNearNullSpace
1864: @*/
1865: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
1866: {

1870:   PetscTryMethod(pc,"PCSetCoordinates_C",(PC,PetscInt,PetscInt,PetscReal*),(pc,dim,nloc,coords));
1871:   return(0);
1872: }