Actual source code: precon.c

petsc-3.7.7 2017-09-25
Report Typos and Errors
  2: /*
  3:     The PC (preconditioner) interface routines, callable by users.
  4: */
  5: #include <petsc/private/pcimpl.h>            /*I "petscksp.h" I*/
  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;

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

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

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

 62:    Collective on PC

 64:    Input Parameter:
 65: .  pc - the preconditioner context

 67:    Level: developer

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

 71: .keywords: PC, destroy

 73: .seealso: PCCreate(), PCSetUp()
 74: @*/
 75: PetscErrorCode  PCReset(PC pc)
 76: {

 81:   if (pc->ops->reset) {
 82:     (*pc->ops->reset)(pc);
 83:   }
 84:   VecDestroy(&pc->diagonalscaleright);
 85:   VecDestroy(&pc->diagonalscaleleft);
 86:   MatDestroy(&pc->pmat);
 87:   MatDestroy(&pc->mat);

 89:   pc->setupcalled = 0;
 90:   return(0);
 91: }

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

 98:    Collective on PC

100:    Input Parameter:
101: .  pc - the preconditioner context

103:    Level: developer

105: .keywords: PC, destroy

107: .seealso: PCCreate(), PCSetUp()
108: @*/
109: PetscErrorCode  PCDestroy(PC *pc)
110: {

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

118:   PCReset(*pc);

120:   /* if memory was published with SAWs then destroy it */
121:   PetscObjectSAWsViewOff((PetscObject)*pc);
122:   if ((*pc)->ops->destroy) {(*(*pc)->ops->destroy)((*pc));}
123:   DMDestroy(&(*pc)->dm);
124:   PetscHeaderDestroy(pc);
125:   return(0);
126: }

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

134:    Logically Collective on PC

136:    Input Parameter:
137: .  pc - the preconditioner context

139:    Output Parameter:
140: .  flag - PETSC_TRUE if it applies the scaling

142:    Level: developer

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

148: .keywords: PC

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

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

167:    Logically Collective on PC

169:    Input Parameters:
170: +  pc - the preconditioner context
171: -  s - scaling vector

173:    Level: intermediate

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

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

181: .keywords: PC

183: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCGetDiagonalScale()
184: @*/
185: PetscErrorCode  PCSetDiagonalScale(PC pc,Vec s)
186: {

192:   pc->diagonalscale     = PETSC_TRUE;

194:   PetscObjectReference((PetscObject)s);
195:   VecDestroy(&pc->diagonalscaleleft);

197:   pc->diagonalscaleleft = s;

199:   VecDuplicate(s,&pc->diagonalscaleright);
200:   VecCopy(s,pc->diagonalscaleright);
201:   VecReciprocal(pc->diagonalscaleright);
202:   return(0);
203: }

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

210:    Logically Collective on PC

212:    Input Parameters:
213: +  pc - the preconditioner context
214: .  in - input vector
215: +  out - scaled vector (maybe the same as in)

217:    Level: intermediate

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

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

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

227: .keywords: PC

229: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleSet(), PCDiagonalScaleRight(), PCDiagonalScale()
230: @*/
231: PetscErrorCode  PCDiagonalScaleLeft(PC pc,Vec in,Vec out)
232: {

239:   if (pc->diagonalscale) {
240:     VecPointwiseMult(out,pc->diagonalscaleleft,in);
241:   } else if (in != out) {
242:     VecCopy(in,out);
243:   }
244:   return(0);
245: }

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

252:    Logically Collective on PC

254:    Input Parameters:
255: +  pc - the preconditioner context
256: .  in - input vector
257: +  out - scaled vector (maybe the same as in)

259:    Level: intermediate

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

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

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

269: .keywords: PC

271: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleSet(), PCDiagonalScale()
272: @*/
273: PetscErrorCode  PCDiagonalScaleRight(PC pc,Vec in,Vec out)
274: {

281:   if (pc->diagonalscale) {
282:     VecPointwiseMult(out,pc->diagonalscaleright,in);
283:   } else if (in != out) {
284:     VecCopy(in,out);
285:   }
286:   return(0);
287: }

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

296:    Logically Collective on PC

298:    Input Parameters:
299: +  pc - the preconditioner context
300: -  flg - PETSC_TRUE to use the Amat, PETSC_FALSE to use the Pmat (default is false)

302:    Options Database Key:
303: .  -pc_use_amat <true,false>

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

309:    Level: intermediate

311: .seealso: PCGetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
312: @*/
313: PetscErrorCode  PCSetUseAmat(PC pc,PetscBool flg)
314: {
317:   pc->useAmat = flg;
318:   return(0);
319: }

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

326:    Logically Collective on PC

328:    Input Parameters:
329: +  pc - iterative context obtained from PCCreate()
330: -  flg - PETSC_TRUE indicates you want the error generated

332:    Level: advanced

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

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

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

343: .seealso: PCGetInitialGuessNonzero(), PCSetInitialGuessKnoll(), PCGetInitialGuessKnoll()
344: @*/
345: PetscErrorCode  PCSetErrorIfFailure(PC pc,PetscBool flg)
346: {
350:   pc->erroriffailure = flg;
351:   return(0);
352: }

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

361:    Logically Collective on PC

363:    Input Parameter:
364: .  pc - the preconditioner context

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

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

373:    Level: intermediate

375: .seealso: PCSetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
376: @*/
377: PetscErrorCode  PCGetUseAmat(PC pc,PetscBool *flg)
378: {
381:   *flg = pc->useAmat;
382:   return(0);
383: }

387: /*@
388:    PCCreate - Creates a preconditioner context.

390:    Collective on MPI_Comm

392:    Input Parameter:
393: .  comm - MPI communicator

395:    Output Parameter:
396: .  pc - location to put the preconditioner context

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

402:    Level: developer

404: .keywords: PC, create, context

406: .seealso: PCSetUp(), PCApply(), PCDestroy()
407: @*/
408: PetscErrorCode  PCCreate(MPI_Comm comm,PC *newpc)
409: {
410:   PC             pc;

415:   *newpc = 0;
416:   PCInitializePackage();

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

420:   pc->mat                  = 0;
421:   pc->pmat                 = 0;
422:   pc->setupcalled          = 0;
423:   pc->setfromoptionscalled = 0;
424:   pc->data                 = 0;
425:   pc->diagonalscale        = PETSC_FALSE;
426:   pc->diagonalscaleleft    = 0;
427:   pc->diagonalscaleright   = 0;

429:   pc->modifysubmatrices  = 0;
430:   pc->modifysubmatricesP = 0;

432:   *newpc = pc;
433:   return(0);

435: }

437: /* -------------------------------------------------------------------------------*/

441: /*@
442:    PCApply - Applies the preconditioner to a vector.

444:    Collective on PC and Vec

446:    Input Parameters:
447: +  pc - the preconditioner context
448: -  x - input vector

450:    Output Parameter:
451: .  y - output vector

453:    Level: developer

455: .keywords: PC, apply

457: .seealso: PCApplyTranspose(), PCApplyBAorAB()
458: @*/
459: PetscErrorCode  PCApply(PC pc,Vec x,Vec y)
460: {
462:   PetscInt       m,n,mv,nv;

468:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
469:   if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
470:   /* use pmat to check vector sizes since for KSPLQR the pmat may be of a different size than mat */
471:   MatGetLocalSize(pc->pmat,&m,&n);
472:   VecGetLocalSize(x,&nv);
473:   VecGetLocalSize(y,&mv);
474:   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);
475:   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);
476:   VecLocked(y,3);

478:   PCSetUp(pc);
479:   if (!pc->ops->apply) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply");
480:   VecLockPush(x);
481:   PetscLogEventBegin(PC_Apply,pc,x,y,0);
482:   (*pc->ops->apply)(pc,x,y);
483:   PetscLogEventEnd(PC_Apply,pc,x,y,0);
484:   if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
485:   VecLockPop(x);
486:   return(0);
487: }

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

494:    Collective on PC and Vec

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

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

503:    Notes:
504:    Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.

506:    Level: developer

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

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

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

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

538:    Collective on PC and Vec

540:    Input Parameters:
541: +  pc - the preconditioner context
542: -  x - input vector

544:    Output Parameter:
545: .  y - output vector

547:    Level: developer

549:    Notes:
550:    Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.

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

554: .seealso: PCApply(), PCApplySymmetricLeft()
555: @*/
556: PetscErrorCode  PCApplySymmetricRight(PC pc,Vec x,Vec y)
557: {

564:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
565:   if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
566:   PCSetUp(pc);
567:   if (!pc->ops->applysymmetricright) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have left symmetric apply");
568:   VecLockPush(x);
569:   PetscLogEventBegin(PC_ApplySymmetricRight,pc,x,y,0);
570:   (*pc->ops->applysymmetricright)(pc,x,y);
571:   PetscLogEventEnd(PC_ApplySymmetricRight,pc,x,y,0);
572:   VecLockPop(x);
573:   if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
574:   return(0);
575: }

579: /*@
580:    PCApplyTranspose - Applies the transpose of preconditioner to a vector.

582:    Collective on PC and Vec

584:    Input Parameters:
585: +  pc - the preconditioner context
586: -  x - input vector

588:    Output Parameter:
589: .  y - output vector

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

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

595:    Level: developer

597: .keywords: PC, apply, transpose

599: .seealso: PCApply(), PCApplyBAorAB(), PCApplyBAorABTranspose(), PCApplyTransposeExists()
600: @*/
601: PetscErrorCode  PCApplyTranspose(PC pc,Vec x,Vec y)
602: {

609:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
610:   if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
611:   PCSetUp(pc);
612:   if (!pc->ops->applytranspose) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply transpose");
613:   VecLockPush(x);
614:   PetscLogEventBegin(PC_Apply,pc,x,y,0);
615:   (*pc->ops->applytranspose)(pc,x,y);
616:   PetscLogEventEnd(PC_Apply,pc,x,y,0);
617:   VecLockPop(x);
618:   if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
619:   return(0);
620: }

624: /*@
625:    PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation

627:    Collective on PC and Vec

629:    Input Parameters:
630: .  pc - the preconditioner context

632:    Output Parameter:
633: .  flg - PETSC_TRUE if a transpose operation is defined

635:    Level: developer

637: .keywords: PC, apply, transpose

639: .seealso: PCApplyTranspose()
640: @*/
641: PetscErrorCode  PCApplyTransposeExists(PC pc,PetscBool  *flg)
642: {
646:   if (pc->ops->applytranspose) *flg = PETSC_TRUE;
647:   else *flg = PETSC_FALSE;
648:   return(0);
649: }

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

656:    Collective on PC and Vec

658:    Input Parameters:
659: +  pc - the preconditioner context
660: .  side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
661: .  x - input vector
662: -  work - work vector

664:    Output Parameter:
665: .  y - output vector

667:    Level: developer

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

672: .keywords: PC, apply, operator

674: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorABTranspose()
675: @*/
676: PetscErrorCode  PCApplyBAorAB(PC pc,PCSide side,Vec x,Vec y,Vec work)
677: {

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

690:   PCSetUp(pc);
691:   if (pc->diagonalscale) {
692:     if (pc->ops->applyBA) {
693:       Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
694:       VecDuplicate(x,&work2);
695:       PCDiagonalScaleRight(pc,x,work2);
696:       (*pc->ops->applyBA)(pc,side,work2,y,work);
697:       PCDiagonalScaleLeft(pc,y,y);
698:       VecDestroy(&work2);
699:     } else if (side == PC_RIGHT) {
700:       PCDiagonalScaleRight(pc,x,y);
701:       PCApply(pc,y,work);
702:       MatMult(pc->mat,work,y);
703:       PCDiagonalScaleLeft(pc,y,y);
704:     } else if (side == PC_LEFT) {
705:       PCDiagonalScaleRight(pc,x,y);
706:       MatMult(pc->mat,y,work);
707:       PCApply(pc,work,y);
708:       PCDiagonalScaleLeft(pc,y,y);
709:     } else if (side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot provide diagonal scaling with symmetric application of preconditioner");
710:   } else {
711:     if (pc->ops->applyBA) {
712:       (*pc->ops->applyBA)(pc,side,x,y,work);
713:     } else if (side == PC_RIGHT) {
714:       PCApply(pc,x,work);
715:       MatMult(pc->mat,work,y);
716:     } else if (side == PC_LEFT) {
717:       MatMult(pc->mat,x,work);
718:       PCApply(pc,work,y);
719:     } else if (side == PC_SYMMETRIC) {
720:       /* There's an extra copy here; maybe should provide 2 work vectors instead? */
721:       PCApplySymmetricRight(pc,x,work);
722:       MatMult(pc->mat,work,y);
723:       VecCopy(y,work);
724:       PCApplySymmetricLeft(pc,work,y);
725:     }
726:   }
727:   if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
728:   return(0);
729: }

733: /*@
734:    PCApplyBAorABTranspose - Applies the transpose of the preconditioner
735:    and operator to a vector. That is, applies tr(B) * tr(A) with left preconditioning,
736:    NOT tr(B*A) = tr(A)*tr(B).

738:    Collective on PC and Vec

740:    Input Parameters:
741: +  pc - the preconditioner context
742: .  side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
743: .  x - input vector
744: -  work - work vector

746:    Output Parameter:
747: .  y - output vector


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

753:     Level: developer

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

757: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorAB()
758: @*/
759: PetscErrorCode  PCApplyBAorABTranspose(PC pc,PCSide side,Vec x,Vec y,Vec work)
760: {

768:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
769:   if (pc->erroriffailure) {VecValidValues(x,3,PETSC_TRUE);}
770:   if (pc->ops->applyBAtranspose) {
771:     (*pc->ops->applyBAtranspose)(pc,side,x,y,work);
772:     if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
773:     return(0);
774:   }
775:   if (side != PC_LEFT && side != PC_RIGHT) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Side must be right or left");

777:   PCSetUp(pc);
778:   if (side == PC_RIGHT) {
779:     PCApplyTranspose(pc,x,work);
780:     MatMultTranspose(pc->mat,work,y);
781:   } else if (side == PC_LEFT) {
782:     MatMultTranspose(pc->mat,x,work);
783:     PCApplyTranspose(pc,work,y);
784:   }
785:   /* add support for PC_SYMMETRIC */
786:   if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
787:   return(0);
788: }

790: /* -------------------------------------------------------------------------------*/

794: /*@
795:    PCApplyRichardsonExists - Determines whether a particular preconditioner has a
796:    built-in fast application of Richardson's method.

798:    Not Collective

800:    Input Parameter:
801: .  pc - the preconditioner

803:    Output Parameter:
804: .  exists - PETSC_TRUE or PETSC_FALSE

806:    Level: developer

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

810: .seealso: PCApplyRichardson()
811: @*/
812: PetscErrorCode  PCApplyRichardsonExists(PC pc,PetscBool  *exists)
813: {
817:   if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
818:   else *exists = PETSC_FALSE;
819:   return(0);
820: }

824: /*@
825:    PCApplyRichardson - Applies several steps of Richardson iteration with
826:    the particular preconditioner. This routine is usually used by the
827:    Krylov solvers and not the application code directly.

829:    Collective on PC

831:    Input Parameters:
832: +  pc  - the preconditioner context
833: .  b   - the right hand side
834: .  w   - one work vector
835: .  rtol - relative decrease in residual norm convergence criteria
836: .  abstol - absolute residual norm convergence criteria
837: .  dtol - divergence residual norm increase criteria
838: .  its - the number of iterations to apply.
839: -  guesszero - if the input x contains nonzero initial guess

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

846:    Notes:
847:    Most preconditioners do not support this function. Use the command
848:    PCApplyRichardsonExists() to determine if one does.

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

853:    Level: developer

855: .keywords: PC, apply, Richardson

857: .seealso: PCApplyRichardsonExists()
858: @*/
859: PetscErrorCode  PCApplyRichardson(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
860: {

868:   if (b == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"b and y must be different vectors");
869:   PCSetUp(pc);
870:   if (!pc->ops->applyrichardson) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply richardson");
871:   (*pc->ops->applyrichardson)(pc,b,y,w,rtol,abstol,dtol,its,guesszero,outits,reason);
872:   return(0);
873: }

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

880:    Logically Collective on PC

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

885:    Output Parameter:
886: .  reason - the reason it failed, currently only -1 

888:    Level: advanced

890: .keywords: PC, setup

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


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

913:    Collective on PC

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

918:    Level: developer

920: .keywords: PC, setup

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

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

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

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

959:   if (!((PetscObject)pc)->type_name) {
960:     PCGetDefaultType_Private(pc,&def);
961:     PCSetType(pc,def);
962:   }

964:   MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
965:   MatSetErrorIfFailure(pc->mat,pc->erroriffailure);
966:   PetscLogEventBegin(PC_SetUp,pc,0,0,0);
967:   if (pc->ops->setup) {
968:     (*pc->ops->setup)(pc);
969:   }
970:   PetscLogEventEnd(PC_SetUp,pc,0,0,0);
971:   if (!pc->setupcalled) pc->setupcalled = 1;
972:   return(0);
973: }

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

982:    Collective on PC

984:    Input Parameters:
985: .  pc - the preconditioner context

987:    Level: developer

989: .keywords: PC, setup, blocks

991: .seealso: PCCreate(), PCApply(), PCDestroy(), PCSetUp()
992: @*/
993: PetscErrorCode  PCSetUpOnBlocks(PC pc)
994: {

999:   if (!pc->ops->setuponblocks) return(0);
1000:   PetscLogEventBegin(PC_SetUpOnBlocks,pc,0,0,0);
1001:   (*pc->ops->setuponblocks)(pc);
1002:   PetscLogEventEnd(PC_SetUpOnBlocks,pc,0,0,0);
1003:   return(0);
1004: }

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

1015:    Logically Collective on PC

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

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

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 func routine (may be null)

1033:    Notes:
1034:    PCSetModifySubMatrices() MUST be called before KSPSetUp() and
1035:    KSPSolve().

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

1041:    Level: advanced

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

1045: .seealso: PCModifySubMatrices(), PCASMGetSubMatrices()
1046: @*/
1047: PetscErrorCode  PCSetModifySubMatrices(PC pc,PetscErrorCode (*func)(PC,PetscInt,const IS[],const IS[],Mat[],void*),void *ctx)
1048: {
1051:   pc->modifysubmatrices  = func;
1052:   pc->modifysubmatricesP = ctx;
1053:   return(0);
1054: }

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

1062:    Collective on PC

1064:    Input Parameters:
1065: +  pc - the preconditioner context
1066: .  nsub - the number of local submatrices
1067: .  row - an array of index sets that contain the global row numbers
1068:          that comprise each local submatrix
1069: .  col - an array of index sets that contain the global column numbers
1070:          that comprise each local submatrix
1071: .  submat - array of local submatrices
1072: -  ctx - optional user-defined context for private data for the
1073:          user-defined routine (may be null)

1075:    Output Parameter:
1076: .  submat - array of local submatrices (the entries of which may
1077:             have been modified)

1079:    Notes:
1080:    The user should NOT generally call this routine, as it will
1081:    automatically be called within certain preconditioners (currently
1082:    block Jacobi, additive Schwarz) if set.

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

1089:    Level: developer

1091: .keywords: PC, modify, submatrices

1093: .seealso: PCSetModifySubMatrices()
1094: @*/
1095: PetscErrorCode  PCModifySubMatrices(PC pc,PetscInt nsub,const IS row[],const IS col[],Mat submat[],void *ctx)
1096: {

1101:   if (!pc->modifysubmatrices) return(0);
1102:   PetscLogEventBegin(PC_ModifySubMatrices,pc,0,0,0);
1103:   (*pc->modifysubmatrices)(pc,nsub,row,col,submat,ctx);
1104:   PetscLogEventEnd(PC_ModifySubMatrices,pc,0,0,0);
1105:   return(0);
1106: }

1110: /*@
1111:    PCSetOperators - Sets the matrix associated with the linear system and
1112:    a (possibly) different one associated with the preconditioner.

1114:    Logically Collective on PC and Mat

1116:    Input Parameters:
1117: +  pc - the preconditioner context
1118: .  Amat - the matrix that defines the linear system
1119: -  Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.

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

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

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

1134:    Level: intermediate

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

1138: .seealso: PCGetOperators(), MatZeroEntries()
1139:  @*/
1140: PetscErrorCode  PCSetOperators(PC pc,Mat Amat,Mat Pmat)
1141: {
1142:   PetscErrorCode   ierr;
1143:   PetscInt         m1,n1,m2,n2;

1151:   if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
1152:     MatGetLocalSize(Amat,&m1,&n1);
1153:     MatGetLocalSize(pc->mat,&m2,&n2);
1154:     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);
1155:     MatGetLocalSize(Pmat,&m1,&n1);
1156:     MatGetLocalSize(pc->pmat,&m2,&n2);
1157:     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);
1158:   }

1160:   if (Pmat != pc->pmat) {
1161:     /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
1162:     pc->matnonzerostate = -1;
1163:     pc->matstate        = -1;
1164:   }

1166:   /* reference first in case the matrices are the same */
1167:   if (Amat) {PetscObjectReference((PetscObject)Amat);}
1168:   MatDestroy(&pc->mat);
1169:   if (Pmat) {PetscObjectReference((PetscObject)Pmat);}
1170:   MatDestroy(&pc->pmat);
1171:   pc->mat  = Amat;
1172:   pc->pmat = Pmat;
1173:   return(0);
1174: }

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

1181:    Logically Collective on PC

1183:    Input Parameters:
1184: +  pc - the preconditioner context
1185: -  flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner

1187:     Level: intermediate

1189: .seealso: PCGetOperators(), MatZeroEntries(), PCGetReusePreconditioner(), KSPSetReusePreconditioner()
1190:  @*/
1191: PetscErrorCode  PCSetReusePreconditioner(PC pc,PetscBool flag)
1192: {
1195:   pc->reusepreconditioner = flag;
1196:   return(0);
1197: }

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

1204:    Not Collective

1206:    Input Parameter:
1207: .  pc - the preconditioner context

1209:    Output Parameter:
1210: .  flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner

1212:    Level: intermediate

1214: .seealso: PCGetOperators(), MatZeroEntries(), PCSetReusePreconditioner()
1215:  @*/
1216: PetscErrorCode  PCGetReusePreconditioner(PC pc,PetscBool *flag)
1217: {
1220:   *flag = pc->reusepreconditioner;
1221:   return(0);
1222: }

1226: /*@C
1227:    PCGetOperators - Gets the matrix associated with the linear system and
1228:    possibly a different one associated with the preconditioner.

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

1232:    Input Parameter:
1233: .  pc - the preconditioner context

1235:    Output Parameters:
1236: +  Amat - the matrix defining the linear system
1237: -  Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.

1239:    Level: intermediate

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

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

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

1254: $         MatCreate(comm,&mat);
1255: $         KSP/PCSetOperators(ksp/pc,Amat,Amat);
1256: $         PetscObjectDereference((PetscObject)mat);
1257: $           set size, type, etc of Amat

1259:      and

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

1264: $         MatCreate(comm,&Amat);
1265: $         MatCreate(comm,&Pmat);
1266: $         KSP/PCSetOperators(ksp/pc,Amat,Pmat);
1267: $         PetscObjectDereference((PetscObject)Amat);
1268: $         PetscObjectDereference((PetscObject)Pmat);
1269: $           set size, type, etc of Amat and Pmat

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


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

1283: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperatorsSet()
1284: @*/
1285: PetscErrorCode  PCGetOperators(PC pc,Mat *Amat,Mat *Pmat)
1286: {

1291:   if (Amat) {
1292:     if (!pc->mat) {
1293:       if (pc->pmat && !Pmat) {  /* Apmat has been set, but user did not request it, so use for Amat */
1294:         pc->mat = pc->pmat;
1295:         PetscObjectReference((PetscObject)pc->mat);
1296:       } else {                  /* both Amat and Pmat are empty */
1297:         MatCreate(PetscObjectComm((PetscObject)pc),&pc->mat);
1298:         if (!Pmat) { /* user did NOT request Pmat, so make same as Amat */
1299:           pc->pmat = pc->mat;
1300:           PetscObjectReference((PetscObject)pc->pmat);
1301:         }
1302:       }
1303:     }
1304:     *Amat = pc->mat;
1305:   }
1306:   if (Pmat) {
1307:     if (!pc->pmat) {
1308:       if (pc->mat && !Amat) {    /* Amat has been set but was not requested, so use for pmat */
1309:         pc->pmat = pc->mat;
1310:         PetscObjectReference((PetscObject)pc->pmat);
1311:       } else {
1312:         MatCreate(PetscObjectComm((PetscObject)pc),&pc->pmat);
1313:         if (!Amat) { /* user did NOT request Amat, so make same as Pmat */
1314:           pc->mat = pc->pmat;
1315:           PetscObjectReference((PetscObject)pc->mat);
1316:         }
1317:       }
1318:     }
1319:     *Pmat = pc->pmat;
1320:   }
1321:   return(0);
1322: }

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

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

1332:    Input Parameter:
1333: .  pc - the preconditioner context

1335:    Output Parameters:
1336: +  mat - the matrix associated with the linear system was set
1337: -  pmat - matrix associated with the preconditioner was set, usually the same

1339:    Level: intermediate

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

1343: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperators()
1344: @*/
1345: PetscErrorCode  PCGetOperatorsSet(PC pc,PetscBool  *mat,PetscBool  *pmat)
1346: {
1349:   if (mat) *mat = (pc->mat)  ? PETSC_TRUE : PETSC_FALSE;
1350:   if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1351:   return(0);
1352: }

1356: /*@
1357:    PCFactorGetMatrix - Gets the factored matrix from the
1358:    preconditioner context.  This routine is valid only for the LU,
1359:    incomplete LU, Cholesky, and incomplete Cholesky methods.

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

1363:    Input Parameters:
1364: .  pc - the preconditioner context

1366:    Output parameters:
1367: .  mat - the factored matrix

1369:    Level: advanced

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

1373: .keywords: PC, get, factored, matrix
1374: @*/
1375: PetscErrorCode  PCFactorGetMatrix(PC pc,Mat *mat)
1376: {

1382:   if (pc->ops->getfactoredmatrix) {
1383:     (*pc->ops->getfactoredmatrix)(pc,mat);
1384:   } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC type does not support getting factor matrix");
1385:   return(0);
1386: }

1390: /*@C
1391:    PCSetOptionsPrefix - Sets the prefix used for searching for all
1392:    PC options in the database.

1394:    Logically Collective on PC

1396:    Input Parameters:
1397: +  pc - the preconditioner context
1398: -  prefix - the prefix string to prepend to all PC option requests

1400:    Notes:
1401:    A hyphen (-) must NOT be given at the beginning of the prefix name.
1402:    The first character of all runtime options is AUTOMATICALLY the
1403:    hyphen.

1405:    Level: advanced

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

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

1417:   PetscObjectSetOptionsPrefix((PetscObject)pc,prefix);
1418:   return(0);
1419: }

1423: /*@C
1424:    PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1425:    PC options in the database.

1427:    Logically Collective on PC

1429:    Input Parameters:
1430: +  pc - the preconditioner context
1431: -  prefix - the prefix string to prepend to all PC option requests

1433:    Notes:
1434:    A hyphen (-) must NOT be given at the beginning of the prefix name.
1435:    The first character of all runtime options is AUTOMATICALLY the
1436:    hyphen.

1438:    Level: advanced

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

1442: .seealso: PCSetOptionsPrefix(), PCGetOptionsPrefix()
1443: @*/
1444: PetscErrorCode  PCAppendOptionsPrefix(PC pc,const char prefix[])
1445: {

1450:   PetscObjectAppendOptionsPrefix((PetscObject)pc,prefix);
1451:   return(0);
1452: }

1456: /*@C
1457:    PCGetOptionsPrefix - Gets the prefix used for searching for all
1458:    PC options in the database.

1460:    Not Collective

1462:    Input Parameters:
1463: .  pc - the preconditioner context

1465:    Output Parameters:
1466: .  prefix - pointer to the prefix string used, is returned

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

1471:    Level: advanced

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

1475: .seealso: PCSetOptionsPrefix(), PCAppendOptionsPrefix()
1476: @*/
1477: PetscErrorCode  PCGetOptionsPrefix(PC pc,const char *prefix[])
1478: {

1484:   PetscObjectGetOptionsPrefix((PetscObject)pc,prefix);
1485:   return(0);
1486: }

1490: PETSC_INTERN PetscErrorCode  PCPreSolveChangeRHS(PC pc,PetscBool *change)
1491: {

1497:   *change = PETSC_FALSE;
1498:   PetscTryMethod(pc,"PCPreSolveChangeRHS_C",(PC,PetscBool*),(pc,change));
1499:   return(0);
1500: }

1504: /*@
1505:    PCPreSolve - Optional pre-solve phase, intended for any
1506:    preconditioner-specific actions that must be performed before
1507:    the iterative solve itself.

1509:    Collective on PC

1511:    Input Parameters:
1512: +  pc - the preconditioner context
1513: -  ksp - the Krylov subspace context

1515:    Level: developer

1517:    Sample of Usage:
1518: .vb
1519:     PCPreSolve(pc,ksp);
1520:     KSPSolve(ksp,b,x);
1521:     PCPostSolve(pc,ksp);
1522: .ve

1524:    Notes:
1525:    The pre-solve phase is distinct from the PCSetUp() phase.

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

1529: .keywords: PC, pre-solve

1531: .seealso: PCPostSolve()
1532: @*/
1533: PetscErrorCode  PCPreSolve(PC pc,KSP ksp)
1534: {
1536:   Vec            x,rhs;

1541:   pc->presolvedone++;
1542:   if (pc->presolvedone > 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot embed PCPreSolve() more than twice");
1543:   KSPGetSolution(ksp,&x);
1544:   KSPGetRhs(ksp,&rhs);

1546:   if (pc->ops->presolve) {
1547:     (*pc->ops->presolve)(pc,ksp,rhs,x);
1548:   }
1549:   return(0);
1550: }

1554: /*@
1555:    PCPostSolve - Optional post-solve phase, intended for any
1556:    preconditioner-specific actions that must be performed after
1557:    the iterative solve itself.

1559:    Collective on PC

1561:    Input Parameters:
1562: +  pc - the preconditioner context
1563: -  ksp - the Krylov subspace context

1565:    Sample of Usage:
1566: .vb
1567:     PCPreSolve(pc,ksp);
1568:     KSPSolve(ksp,b,x);
1569:     PCPostSolve(pc,ksp);
1570: .ve

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

1575:    Level: developer

1577: .keywords: PC, post-solve

1579: .seealso: PCPreSolve(), KSPSolve()
1580: @*/
1581: PetscErrorCode  PCPostSolve(PC pc,KSP ksp)
1582: {
1584:   Vec            x,rhs;

1589:   pc->presolvedone--;
1590:   KSPGetSolution(ksp,&x);
1591:   KSPGetRhs(ksp,&rhs);
1592:   if (pc->ops->postsolve) {
1593:      (*pc->ops->postsolve)(pc,ksp,rhs,x);
1594:   }
1595:   return(0);
1596: }

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

1603:   Collective on PetscViewer

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

1610:    Level: intermediate

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

1615:   Notes for advanced users:
1616:   Most users should not need to know the details of the binary storage
1617:   format, since PCLoad() and PCView() completely hide these details.
1618:   But for anyone who's interested, the standard binary matrix storage
1619:   format is
1620: .vb
1621:      has not yet been determined
1622: .ve

1624: .seealso: PetscViewerBinaryOpen(), PCView(), MatLoad(), VecLoad()
1625: @*/
1626: PetscErrorCode  PCLoad(PC newdm, PetscViewer viewer)
1627: {
1629:   PetscBool      isbinary;
1630:   PetscInt       classid;
1631:   char           type[256];

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

1639:   PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
1640:   if (classid != PC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not PC next in file");
1641:   PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
1642:   PCSetType(newdm, type);
1643:   if (newdm->ops->load) {
1644:     (*newdm->ops->load)(newdm,viewer);
1645:   }
1646:   return(0);
1647: }

1649: #include <petscdraw.h>
1650: #if defined(PETSC_HAVE_SAWS)
1651: #include <petscviewersaws.h>
1652: #endif
1655: /*@C
1656:    PCView - Prints the PC data structure.

1658:    Collective on PC

1660:    Input Parameters:
1661: +  PC - the PC context
1662: -  viewer - optional visualization context

1664:    Note:
1665:    The available visualization contexts include
1666: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
1667: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1668:          output where only the first processor opens
1669:          the file.  All other processors send their
1670:          data to the first processor to print.

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

1675:    Level: developer

1677: .keywords: PC, view

1679: .seealso: KSPView(), PetscViewerASCIIOpen()
1680: @*/
1681: PetscErrorCode  PCView(PC pc,PetscViewer viewer)
1682: {
1683:   PCType            cstr;
1684:   PetscErrorCode    ierr;
1685:   PetscBool         iascii,isstring,isbinary,isdraw;
1686:   PetscViewerFormat format;
1687: #if defined(PETSC_HAVE_SAWS)
1688:   PetscBool         issaws;
1689: #endif

1693:   if (!viewer) {
1694:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);
1695:   }

1699:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1700:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
1701:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1702:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1703: #if defined(PETSC_HAVE_SAWS)
1704:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws);
1705: #endif

1707:   if (iascii) {
1708:     PetscViewerGetFormat(viewer,&format);
1709:     PetscObjectPrintClassNamePrefixType((PetscObject)pc,viewer);
1710:     if (!pc->setupcalled) {
1711:       PetscViewerASCIIPrintf(viewer,"  PC has not been set up so information may be incomplete\n");
1712:     }
1713:     if (pc->ops->view) {
1714:       PetscViewerASCIIPushTab(viewer);
1715:       (*pc->ops->view)(pc,viewer);
1716:       PetscViewerASCIIPopTab(viewer);
1717:     }
1718:     if (pc->mat) {
1719:       PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
1720:       if (pc->pmat == pc->mat) {
1721:         PetscViewerASCIIPrintf(viewer,"  linear system matrix = precond matrix:\n");
1722:         PetscViewerASCIIPushTab(viewer);
1723:         MatView(pc->mat,viewer);
1724:         PetscViewerASCIIPopTab(viewer);
1725:       } else {
1726:         if (pc->pmat) {
1727:           PetscViewerASCIIPrintf(viewer,"  linear system matrix followed by preconditioner matrix:\n");
1728:         } else {
1729:           PetscViewerASCIIPrintf(viewer,"  linear system matrix:\n");
1730:         }
1731:         PetscViewerASCIIPushTab(viewer);
1732:         MatView(pc->mat,viewer);
1733:         if (pc->pmat) {MatView(pc->pmat,viewer);}
1734:         PetscViewerASCIIPopTab(viewer);
1735:       }
1736:       PetscViewerPopFormat(viewer);
1737:     }
1738:   } else if (isstring) {
1739:     PCGetType(pc,&cstr);
1740:     PetscViewerStringSPrintf(viewer," %-7.7s",cstr);
1741:     if (pc->ops->view) {(*pc->ops->view)(pc,viewer);}
1742:   } else if (isbinary) {
1743:     PetscInt    classid = PC_FILE_CLASSID;
1744:     MPI_Comm    comm;
1745:     PetscMPIInt rank;
1746:     char        type[256];

1748:     PetscObjectGetComm((PetscObject)pc,&comm);
1749:     MPI_Comm_rank(comm,&rank);
1750:     if (!rank) {
1751:       PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
1752:       PetscStrncpy(type,((PetscObject)pc)->type_name,256);
1753:       PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);
1754:     }
1755:     if (pc->ops->view) {
1756:       (*pc->ops->view)(pc,viewer);
1757:     }
1758:   } else if (isdraw) {
1759:     PetscDraw draw;
1760:     char      str[25];
1761:     PetscReal x,y,bottom,h;
1762:     PetscInt  n;

1764:     PetscViewerDrawGetDraw(viewer,0,&draw);
1765:     PetscDrawGetCurrentPoint(draw,&x,&y);
1766:     if (pc->mat) {
1767:       MatGetSize(pc->mat,&n,NULL);
1768:       PetscSNPrintf(str,25,"PC: %s (%D)",((PetscObject)pc)->type_name,n);
1769:     } else {
1770:       PetscSNPrintf(str,25,"PC: %s",((PetscObject)pc)->type_name);
1771:     }
1772:     PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
1773:     bottom = y - h;
1774:     PetscDrawPushCurrentPoint(draw,x,bottom);
1775:     if (pc->ops->view) {
1776:       (*pc->ops->view)(pc,viewer);
1777:     }
1778:     PetscDrawPopCurrentPoint(draw);
1779: #if defined(PETSC_HAVE_SAWS)
1780:   } else if (issaws) {
1781:     PetscMPIInt rank;

1783:     PetscObjectName((PetscObject)pc);
1784:     MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
1785:     if (!((PetscObject)pc)->amsmem && !rank) {
1786:       PetscObjectViewSAWs((PetscObject)pc,viewer);
1787:     }
1788:     if (pc->mat) {MatView(pc->mat,viewer);}
1789:     if (pc->pmat && pc->pmat != pc->mat) {MatView(pc->pmat,viewer);}
1790: #endif
1791:   }
1792:   return(0);
1793: }

1797: /*@C
1798:   PCRegister -  Adds a method to the preconditioner package.

1800:    Not collective

1802:    Input Parameters:
1803: +  name_solver - name of a new user-defined solver
1804: -  routine_create - routine to create method context

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

1809:    Sample usage:
1810: .vb
1811:    PCRegister("my_solver", MySolverCreate);
1812: .ve

1814:    Then, your solver can be chosen with the procedural interface via
1815: $     PCSetType(pc,"my_solver")
1816:    or at runtime via the option
1817: $     -pc_type my_solver

1819:    Level: advanced

1821: .keywords: PC, register

1823: .seealso: PCRegisterAll(), PCRegisterDestroy()
1824: @*/
1825: PetscErrorCode  PCRegister(const char sname[],PetscErrorCode (*function)(PC))
1826: {

1830:   PetscFunctionListAdd(&PCList,sname,function);
1831:   return(0);
1832: }

1836: /*@
1837:     PCComputeExplicitOperator - Computes the explicit preconditioned operator.

1839:     Collective on PC

1841:     Input Parameter:
1842: .   pc - the preconditioner object

1844:     Output Parameter:
1845: .   mat - the explict preconditioned operator

1847:     Notes:
1848:     This computation is done by applying the operators to columns of the
1849:     identity matrix.

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

1855:     Level: advanced

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

1859: .seealso: KSPComputeExplicitOperator()

1861: @*/
1862: PetscErrorCode  PCComputeExplicitOperator(PC pc,Mat *mat)
1863: {
1864:   Vec            in,out;
1866:   PetscInt       i,M,m,*rows,start,end;
1867:   PetscMPIInt    size;
1868:   MPI_Comm       comm;
1869:   PetscScalar    *array,one = 1.0;


1875:   PetscObjectGetComm((PetscObject)pc,&comm);
1876:   MPI_Comm_size(comm,&size);

1878:   if (!pc->pmat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"You must call KSPSetOperators() or PCSetOperators() before this call");
1879:   MatCreateVecs(pc->pmat,&in,0);
1880:   VecDuplicate(in,&out);
1881:   VecGetOwnershipRange(in,&start,&end);
1882:   VecGetSize(in,&M);
1883:   VecGetLocalSize(in,&m);
1884:   PetscMalloc1(m+1,&rows);
1885:   for (i=0; i<m; i++) rows[i] = start + i;

1887:   MatCreate(comm,mat);
1888:   MatSetSizes(*mat,m,m,M,M);
1889:   if (size == 1) {
1890:     MatSetType(*mat,MATSEQDENSE);
1891:     MatSeqDenseSetPreallocation(*mat,NULL);
1892:   } else {
1893:     MatSetType(*mat,MATMPIAIJ);
1894:     MatMPIAIJSetPreallocation(*mat,0,NULL,0,NULL);
1895:   }
1896:   MatSetOption(*mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);

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

1900:     VecSet(in,0.0);
1901:     VecSetValues(in,1,&i,&one,INSERT_VALUES);
1902:     VecAssemblyBegin(in);
1903:     VecAssemblyEnd(in);

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

1908:     VecGetArray(out,&array);
1909:     MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
1910:     VecRestoreArray(out,&array);

1912:   }
1913:   PetscFree(rows);
1914:   VecDestroy(&out);
1915:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
1916:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
1917:   return(0);
1918: }

1922: /*@
1923:    PCSetCoordinates - sets the coordinates of all the nodes on the local process

1925:    Collective on PC

1927:    Input Parameters:
1928: +  pc - the solver context
1929: .  dim - the dimension of the coordinates 1, 2, or 3
1930: -  coords - the coordinates

1932:    Level: intermediate

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

1942: .seealso: MatSetNearNullSpace
1943: @*/
1944: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
1945: {

1949:   PetscTryMethod(pc,"PCSetCoordinates_C",(PC,PetscInt,PetscInt,PetscReal*),(pc,dim,nloc,coords));
1950:   return(0);
1951: }