Actual source code: precon.c

petsc-3.6.3 2015-12-03
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;

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

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

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

 61:    Collective on PC

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

 66:    Level: developer

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

 70: .keywords: PC, destroy

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

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

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

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

 97:    Collective on PC

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

102:    Level: developer

104: .keywords: PC, destroy

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

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

117:   PCReset(*pc);

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

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

133:    Logically Collective on PC

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

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

141:    Level: developer

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

147: .keywords: PC

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

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

166:    Logically Collective on PC

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

172:    Level: intermediate

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

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

180: .keywords: PC

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

191:   pc->diagonalscale     = PETSC_TRUE;

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

196:   pc->diagonalscaleleft = s;

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

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

209:    Logically Collective on PC

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

216:    Level: intermediate

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

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

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

226: .keywords: PC

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

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

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

251:    Logically Collective on PC

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

258:    Level: intermediate

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

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

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

268: .keywords: PC

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

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

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

295:    Logically Collective on PC

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

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

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

308:    Level: intermediate

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

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

325:    Logically Collective on PC

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

331:    Level: advanced

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

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

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

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

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

360:    Logically Collective on PC

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

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

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

372:    Level: intermediate

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

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

389:    Collective on MPI_Comm

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

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

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

401:    Level: developer

403: .keywords: PC, create, context

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

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

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

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

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

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

434: }

436: /* -------------------------------------------------------------------------------*/

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

443:    Collective on PC and Vec

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

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

452:    Level: developer

454: .keywords: PC, apply

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

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

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

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

495:    Collective on PC and Vec

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

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

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

507:    Level: developer

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

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

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

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

541:    Collective on PC and Vec

543:    Input Parameters:
544: +  pc - the preconditioner context
545: -  x - input vector

547:    Output Parameter:
548: .  y - output vector

550:    Level: developer

552:    Notes:
553:    Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.

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

557: .seealso: PCApply(), PCApplySymmetricLeft()
558: @*/
559: PetscErrorCode  PCApplySymmetricRight(PC pc,Vec x,Vec y)
560: {

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

584: /*@
585:    PCApplyTranspose - Applies the transpose of preconditioner to a vector.

587:    Collective on PC and Vec

589:    Input Parameters:
590: +  pc - the preconditioner context
591: -  x - input vector

593:    Output Parameter:
594: .  y - output vector

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

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

600:    Level: developer

602: .keywords: PC, apply, transpose

604: .seealso: PCApply(), PCApplyBAorAB(), PCApplyBAorABTranspose(), PCApplyTransposeExists()
605: @*/
606: PetscErrorCode  PCApplyTranspose(PC pc,Vec x,Vec y)
607: {

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

631: /*@
632:    PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation

634:    Collective on PC and Vec

636:    Input Parameters:
637: .  pc - the preconditioner context

639:    Output Parameter:
640: .  flg - PETSC_TRUE if a transpose operation is defined

642:    Level: developer

644: .keywords: PC, apply, transpose

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

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

663:    Collective on PC and Vec

665:    Input Parameters:
666: +  pc - the preconditioner context
667: .  side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
668: .  x - input vector
669: -  work - work vector

671:    Output Parameter:
672: .  y - output vector

674:    Level: developer

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

679: .keywords: PC, apply, operator

681: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorABTranspose()
682: @*/
683: PetscErrorCode  PCApplyBAorAB(PC pc,PCSide side,Vec x,Vec y,Vec work)
684: {

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

697:   if (pc->setupcalled < 2) {
698:     PCSetUp(pc);
699:   }

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

743: /*@
744:    PCApplyBAorABTranspose - Applies the transpose of the preconditioner
745:    and operator to a vector. That is, applies tr(B) * tr(A) with left preconditioning,
746:    NOT tr(B*A) = tr(A)*tr(B).

748:    Collective on PC and Vec

750:    Input Parameters:
751: +  pc - the preconditioner context
752: .  side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
753: .  x - input vector
754: -  work - work vector

756:    Output Parameter:
757: .  y - output vector


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

763:     Level: developer

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

767: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorAB()
768: @*/
769: PetscErrorCode  PCApplyBAorABTranspose(PC pc,PCSide side,Vec x,Vec y,Vec work)
770: {

778:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
779:   if (pc->erroriffailure) {VecValidValues(x,3,PETSC_TRUE);}
780:   if (pc->ops->applyBAtranspose) {
781:     (*pc->ops->applyBAtranspose)(pc,side,x,y,work);
782:     if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
783:     return(0);
784:   }
785:   if (side != PC_LEFT && side != PC_RIGHT) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Side must be right or left");

787:   if (pc->setupcalled < 2) {
788:     PCSetUp(pc);
789:   }

791:   if (side == PC_RIGHT) {
792:     PCApplyTranspose(pc,x,work);
793:     MatMultTranspose(pc->mat,work,y);
794:   } else if (side == PC_LEFT) {
795:     MatMultTranspose(pc->mat,x,work);
796:     PCApplyTranspose(pc,work,y);
797:   }
798:   /* add support for PC_SYMMETRIC */
799:   if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
800:   return(0);
801: }

803: /* -------------------------------------------------------------------------------*/

807: /*@
808:    PCApplyRichardsonExists - Determines whether a particular preconditioner has a
809:    built-in fast application of Richardson's method.

811:    Not Collective

813:    Input Parameter:
814: .  pc - the preconditioner

816:    Output Parameter:
817: .  exists - PETSC_TRUE or PETSC_FALSE

819:    Level: developer

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

823: .seealso: PCApplyRichardson()
824: @*/
825: PetscErrorCode  PCApplyRichardsonExists(PC pc,PetscBool  *exists)
826: {
830:   if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
831:   else *exists = PETSC_FALSE;
832:   return(0);
833: }

837: /*@
838:    PCApplyRichardson - Applies several steps of Richardson iteration with
839:    the particular preconditioner. This routine is usually used by the
840:    Krylov solvers and not the application code directly.

842:    Collective on PC

844:    Input Parameters:
845: +  pc  - the preconditioner context
846: .  b   - the right hand side
847: .  w   - one work vector
848: .  rtol - relative decrease in residual norm convergence criteria
849: .  abstol - absolute residual norm convergence criteria
850: .  dtol - divergence residual norm increase criteria
851: .  its - the number of iterations to apply.
852: -  guesszero - if the input x contains nonzero initial guess

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

859:    Notes:
860:    Most preconditioners do not support this function. Use the command
861:    PCApplyRichardsonExists() to determine if one does.

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

866:    Level: developer

868: .keywords: PC, apply, Richardson

870: .seealso: PCApplyRichardsonExists()
871: @*/
872: PetscErrorCode  PCApplyRichardson(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
873: {

881:   if (b == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"b and y must be different vectors");
882:   if (pc->setupcalled < 2) {
883:     PCSetUp(pc);
884:   }
885:   if (!pc->ops->applyrichardson) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply richardson");
886:   (*pc->ops->applyrichardson)(pc,b,y,w,rtol,abstol,dtol,its,guesszero,outits,reason);
887:   return(0);
888: }

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

895:    Logically Collective on PC

897:    Input Parameter:
898: .  pc - the preconditioner context

900:    Output Parameter:
901: .  reason - the reason it failed, currently only -1 

903:    Level: advanced

905: .keywords: PC, setup

907: .seealso: PCCreate(), PCApply(), PCDestroy()
908: @*/
909: PetscErrorCode  PCGetSetUpFailedReason(PC pc,PetscInt *reason)
910: {
912:   if (pc->setupcalled < 0) *reason = pc->setupcalled;
913:   else *reason = 0;
914:   return(0);
915: }


918: /*
919:       a setupcall of 0 indicates never setup,
920:                      1 indicates has been previously setup
921:                     -1 indicates a PCSetUp() was attempted and failed
922: */
925: /*@
926:    PCSetUp - Prepares for the use of a preconditioner.

928:    Collective on PC

930:    Input Parameter:
931: .  pc - the preconditioner context

933:    Level: developer

935: .keywords: PC, setup

937: .seealso: PCCreate(), PCApply(), PCDestroy()
938: @*/
939: PetscErrorCode  PCSetUp(PC pc)
940: {
941:   PetscErrorCode   ierr;
942:   const char       *def;
943:   PetscObjectState matstate, matnonzerostate;

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

949:   if (pc->setupcalled && pc->reusepreconditioner) {
950:     PetscInfo(pc,"Leaving PC with identical preconditioner since reuse preconditioner is set\n");
951:     return(0);
952:   }

954:   PetscObjectStateGet((PetscObject)pc->pmat,&matstate);
955:   MatGetNonzeroState(pc->pmat,&matnonzerostate);
956:   if (!pc->setupcalled) {
957:     PetscInfo(pc,"Setting up PC for first time\n");
958:     pc->flag        = DIFFERENT_NONZERO_PATTERN;
959:   } else if (matstate == pc->matstate) {
960:     PetscInfo(pc,"Leaving PC with identical preconditioner since operator is unchanged\n");
961:     return(0);
962:   } else {
963:     if (matnonzerostate > pc->matnonzerostate) {
964:        PetscInfo(pc,"Setting up PC with different nonzero pattern\n");
965:        pc->flag            = DIFFERENT_NONZERO_PATTERN;
966:     } else {
967:       PetscInfo(pc,"Setting up PC with same nonzero pattern\n");
968:       pc->flag            = SAME_NONZERO_PATTERN;
969:     }
970:   }
971:   pc->matstate        = matstate;
972:   pc->matnonzerostate = matnonzerostate;

974:   if (!((PetscObject)pc)->type_name) {
975:     PCGetDefaultType_Private(pc,&def);
976:     PCSetType(pc,def);
977:   }

979:   MatSetErrorIfFPE(pc->pmat,pc->erroriffailure);
980:   MatSetErrorIfFPE(pc->mat,pc->erroriffailure);
981:   PetscLogEventBegin(PC_SetUp,pc,0,0,0);
982:   if (pc->ops->setup) {
983:     (*pc->ops->setup)(pc);
984:   }
985:   PetscLogEventEnd(PC_SetUp,pc,0,0,0);
986:   if (!pc->setupcalled) pc->setupcalled = 1;
987:   return(0);
988: }

992: /*@
993:    PCSetUpOnBlocks - Sets up the preconditioner for each block in
994:    the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
995:    methods.

997:    Collective on PC

999:    Input Parameters:
1000: .  pc - the preconditioner context

1002:    Level: developer

1004: .keywords: PC, setup, blocks

1006: .seealso: PCCreate(), PCApply(), PCDestroy(), PCSetUp()
1007: @*/
1008: PetscErrorCode  PCSetUpOnBlocks(PC pc)
1009: {

1014:   if (!pc->ops->setuponblocks) return(0);
1015:   PetscLogEventBegin(PC_SetUpOnBlocks,pc,0,0,0);
1016:   (*pc->ops->setuponblocks)(pc);
1017:   PetscLogEventEnd(PC_SetUpOnBlocks,pc,0,0,0);
1018:   return(0);
1019: }

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

1030:    Logically Collective on PC

1032:    Input Parameters:
1033: +  pc - the preconditioner context
1034: .  func - routine for modifying the submatrices
1035: -  ctx - optional user-defined context (may be null)

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

1040: .  row - an array of index sets that contain the global row numbers
1041:          that comprise each local submatrix
1042: .  col - an array of index sets that contain the global column numbers
1043:          that comprise each local submatrix
1044: .  submat - array of local submatrices
1045: -  ctx - optional user-defined context for private data for the
1046:          user-defined func routine (may be null)

1048:    Notes:
1049:    PCSetModifySubMatrices() MUST be called before KSPSetUp() and
1050:    KSPSolve().

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

1056:    Level: advanced

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

1060: .seealso: PCModifySubMatrices(), PCASMGetSubMatrices()
1061: @*/
1062: PetscErrorCode  PCSetModifySubMatrices(PC pc,PetscErrorCode (*func)(PC,PetscInt,const IS[],const IS[],Mat[],void*),void *ctx)
1063: {
1066:   pc->modifysubmatrices  = func;
1067:   pc->modifysubmatricesP = ctx;
1068:   return(0);
1069: }

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

1077:    Collective on PC

1079:    Input Parameters:
1080: +  pc - the preconditioner context
1081: .  nsub - the number of local submatrices
1082: .  row - an array of index sets that contain the global row numbers
1083:          that comprise each local submatrix
1084: .  col - an array of index sets that contain the global column numbers
1085:          that comprise each local submatrix
1086: .  submat - array of local submatrices
1087: -  ctx - optional user-defined context for private data for the
1088:          user-defined routine (may be null)

1090:    Output Parameter:
1091: .  submat - array of local submatrices (the entries of which may
1092:             have been modified)

1094:    Notes:
1095:    The user should NOT generally call this routine, as it will
1096:    automatically be called within certain preconditioners (currently
1097:    block Jacobi, additive Schwarz) if set.

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

1104:    Level: developer

1106: .keywords: PC, modify, submatrices

1108: .seealso: PCSetModifySubMatrices()
1109: @*/
1110: PetscErrorCode  PCModifySubMatrices(PC pc,PetscInt nsub,const IS row[],const IS col[],Mat submat[],void *ctx)
1111: {

1116:   if (!pc->modifysubmatrices) return(0);
1117:   PetscLogEventBegin(PC_ModifySubMatrices,pc,0,0,0);
1118:   (*pc->modifysubmatrices)(pc,nsub,row,col,submat,ctx);
1119:   PetscLogEventEnd(PC_ModifySubMatrices,pc,0,0,0);
1120:   return(0);
1121: }

1125: /*@
1126:    PCSetOperators - Sets the matrix associated with the linear system and
1127:    a (possibly) different one associated with the preconditioner.

1129:    Logically Collective on PC and Mat

1131:    Input Parameters:
1132: +  pc - the preconditioner context
1133: .  Amat - the matrix that defines the linear system
1134: -  Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.

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

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

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

1149:    Level: intermediate

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

1153: .seealso: PCGetOperators(), MatZeroEntries()
1154:  @*/
1155: PetscErrorCode  PCSetOperators(PC pc,Mat Amat,Mat Pmat)
1156: {
1157:   PetscErrorCode   ierr;
1158:   PetscInt         m1,n1,m2,n2;

1166:   if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
1167:     MatGetLocalSize(Amat,&m1,&n1);
1168:     MatGetLocalSize(pc->mat,&m2,&n2);
1169:     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);
1170:     MatGetLocalSize(Pmat,&m1,&n1);
1171:     MatGetLocalSize(pc->pmat,&m2,&n2);
1172:     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);
1173:   }

1175:   if (Pmat != pc->pmat) {
1176:     /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
1177:     pc->matnonzerostate = -1;
1178:     pc->matstate        = -1;
1179:   }

1181:   /* reference first in case the matrices are the same */
1182:   if (Amat) {PetscObjectReference((PetscObject)Amat);}
1183:   MatDestroy(&pc->mat);
1184:   if (Pmat) {PetscObjectReference((PetscObject)Pmat);}
1185:   MatDestroy(&pc->pmat);
1186:   pc->mat  = Amat;
1187:   pc->pmat = Pmat;
1188:   return(0);
1189: }

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

1196:    Logically Collective on PC

1198:    Input Parameters:
1199: +  pc - the preconditioner context
1200: -  flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner

1202:     Level: intermediate

1204: .seealso: PCGetOperators(), MatZeroEntries(), PCGetReusePreconditioner(), KSPSetReusePreconditioner()
1205:  @*/
1206: PetscErrorCode  PCSetReusePreconditioner(PC pc,PetscBool flag)
1207: {
1210:   pc->reusepreconditioner = flag;
1211:   return(0);
1212: }

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

1219:    Not Collective

1221:    Input Parameter:
1222: .  pc - the preconditioner context

1224:    Output Parameter:
1225: .  flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner

1227:    Level: intermediate

1229: .seealso: PCGetOperators(), MatZeroEntries(), PCSetReusePreconditioner()
1230:  @*/
1231: PetscErrorCode  PCGetReusePreconditioner(PC pc,PetscBool *flag)
1232: {
1235:   *flag = pc->reusepreconditioner;
1236:   return(0);
1237: }

1241: /*@C
1242:    PCGetOperators - Gets the matrix associated with the linear system and
1243:    possibly a different one associated with the preconditioner.

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

1247:    Input Parameter:
1248: .  pc - the preconditioner context

1250:    Output Parameters:
1251: +  Amat - the matrix defining the linear system
1252: -  Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.

1254:    Level: intermediate

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

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

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

1269: $         MatCreate(comm,&mat);
1270: $         KSP/PCSetOperators(ksp/pc,Amat,Amat);
1271: $         PetscObjectDereference((PetscObject)mat);
1272: $           set size, type, etc of Amat

1274:      and

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

1279: $         MatCreate(comm,&Amat);
1280: $         MatCreate(comm,&Pmat);
1281: $         KSP/PCSetOperators(ksp/pc,Amat,Pmat);
1282: $         PetscObjectDereference((PetscObject)Amat);
1283: $         PetscObjectDereference((PetscObject)Pmat);
1284: $           set size, type, etc of Amat and Pmat

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


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

1298: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperatorsSet()
1299: @*/
1300: PetscErrorCode  PCGetOperators(PC pc,Mat *Amat,Mat *Pmat)
1301: {

1306:   if (Amat) {
1307:     if (!pc->mat) {
1308:       if (pc->pmat && !Pmat) {  /* Apmat has been set, but user did not request it, so use for Amat */
1309:         pc->mat = pc->pmat;
1310:         PetscObjectReference((PetscObject)pc->mat);
1311:       } else {                  /* both Amat and Pmat are empty */
1312:         MatCreate(PetscObjectComm((PetscObject)pc),&pc->mat);
1313:         if (!Pmat) { /* user did NOT request Pmat, so make same as Amat */
1314:           pc->pmat = pc->mat;
1315:           PetscObjectReference((PetscObject)pc->pmat);
1316:         }
1317:       }
1318:     }
1319:     *Amat = pc->mat;
1320:   }
1321:   if (Pmat) {
1322:     if (!pc->pmat) {
1323:       if (pc->mat && !Amat) {    /* Amat has been set but was not requested, so use for pmat */
1324:         pc->pmat = pc->mat;
1325:         PetscObjectReference((PetscObject)pc->pmat);
1326:       } else {
1327:         MatCreate(PetscObjectComm((PetscObject)pc),&pc->pmat);
1328:         if (!Amat) { /* user did NOT request Amat, so make same as Pmat */
1329:           pc->mat = pc->pmat;
1330:           PetscObjectReference((PetscObject)pc->mat);
1331:         }
1332:       }
1333:     }
1334:     *Pmat = pc->pmat;
1335:   }
1336:   return(0);
1337: }

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

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

1347:    Input Parameter:
1348: .  pc - the preconditioner context

1350:    Output Parameters:
1351: +  mat - the matrix associated with the linear system was set
1352: -  pmat - matrix associated with the preconditioner was set, usually the same

1354:    Level: intermediate

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

1358: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperators()
1359: @*/
1360: PetscErrorCode  PCGetOperatorsSet(PC pc,PetscBool  *mat,PetscBool  *pmat)
1361: {
1364:   if (mat) *mat = (pc->mat)  ? PETSC_TRUE : PETSC_FALSE;
1365:   if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1366:   return(0);
1367: }

1371: /*@
1372:    PCFactorGetMatrix - Gets the factored matrix from the
1373:    preconditioner context.  This routine is valid only for the LU,
1374:    incomplete LU, Cholesky, and incomplete Cholesky methods.

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

1378:    Input Parameters:
1379: .  pc - the preconditioner context

1381:    Output parameters:
1382: .  mat - the factored matrix

1384:    Level: advanced

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

1388: .keywords: PC, get, factored, matrix
1389: @*/
1390: PetscErrorCode  PCFactorGetMatrix(PC pc,Mat *mat)
1391: {

1397:   if (pc->ops->getfactoredmatrix) {
1398:     (*pc->ops->getfactoredmatrix)(pc,mat);
1399:   } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC type does not support getting factor matrix");
1400:   return(0);
1401: }

1405: /*@C
1406:    PCSetOptionsPrefix - Sets the prefix used for searching for all
1407:    PC options in the database.

1409:    Logically Collective on PC

1411:    Input Parameters:
1412: +  pc - the preconditioner context
1413: -  prefix - the prefix string to prepend to all PC option requests

1415:    Notes:
1416:    A hyphen (-) must NOT be given at the beginning of the prefix name.
1417:    The first character of all runtime options is AUTOMATICALLY the
1418:    hyphen.

1420:    Level: advanced

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

1424: .seealso: PCAppendOptionsPrefix(), PCGetOptionsPrefix()
1425: @*/
1426: PetscErrorCode  PCSetOptionsPrefix(PC pc,const char prefix[])
1427: {

1432:   PetscObjectSetOptionsPrefix((PetscObject)pc,prefix);
1433:   return(0);
1434: }

1438: /*@C
1439:    PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1440:    PC options in the database.

1442:    Logically Collective on PC

1444:    Input Parameters:
1445: +  pc - the preconditioner context
1446: -  prefix - the prefix string to prepend to all PC option requests

1448:    Notes:
1449:    A hyphen (-) must NOT be given at the beginning of the prefix name.
1450:    The first character of all runtime options is AUTOMATICALLY the
1451:    hyphen.

1453:    Level: advanced

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

1457: .seealso: PCSetOptionsPrefix(), PCGetOptionsPrefix()
1458: @*/
1459: PetscErrorCode  PCAppendOptionsPrefix(PC pc,const char prefix[])
1460: {

1465:   PetscObjectAppendOptionsPrefix((PetscObject)pc,prefix);
1466:   return(0);
1467: }

1471: /*@C
1472:    PCGetOptionsPrefix - Gets the prefix used for searching for all
1473:    PC options in the database.

1475:    Not Collective

1477:    Input Parameters:
1478: .  pc - the preconditioner context

1480:    Output Parameters:
1481: .  prefix - pointer to the prefix string used, is returned

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

1486:    Level: advanced

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

1490: .seealso: PCSetOptionsPrefix(), PCAppendOptionsPrefix()
1491: @*/
1492: PetscErrorCode  PCGetOptionsPrefix(PC pc,const char *prefix[])
1493: {

1499:   PetscObjectGetOptionsPrefix((PetscObject)pc,prefix);
1500:   return(0);
1501: }

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

1510:    Collective on PC

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

1516:    Level: developer

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

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

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

1530: .keywords: PC, pre-solve

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

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

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

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

1560:    Collective on PC

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

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

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

1576:    Level: developer

1578: .keywords: PC, post-solve

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

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

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

1604:   Collective on PetscViewer

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

1611:    Level: intermediate

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

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

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

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

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

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

1659:    Collective on PC

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

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

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

1676:    Level: developer

1678: .keywords: PC, view

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

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

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

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

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

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

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


1799: /*@
1800:    PCSetInitialGuessNonzero - Tells the iterative solver that the
1801:    initial guess is nonzero; otherwise PC assumes the initial guess
1802:    is to be zero (and thus zeros it out before solving).

1804:    Logically Collective on PC

1806:    Input Parameters:
1807: +  pc - iterative context obtained from PCCreate()
1808: -  flg - PETSC_TRUE indicates the guess is non-zero, PETSC_FALSE indicates the guess is zero

1810:    Level: Developer

1812:    Notes:
1813:     This is a weird function. Since PC's are linear operators on the right hand side they
1814:     CANNOT use an initial guess. This function is for the "pass-through" preconditioners
1815:     PCKSP and PCREDUNDANT  and causes the inner KSP object to use the nonzero
1816:     initial guess. Not currently working for PCREDUNDANT, that has to be rewritten to use KSP.


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

1821: .seealso: PCGetInitialGuessNonzero(), PCSetInitialGuessKnoll(), PCGetInitialGuessKnoll()
1822: @*/
1823: PetscErrorCode  PCSetInitialGuessNonzero(PC pc,PetscBool flg)
1824: {
1827:   pc->nonzero_guess = flg;
1828:   return(0);
1829: }

1833: /*@
1834:    PCGetInitialGuessNonzero - Determines if the iterative solver assumes that the
1835:    initial guess is nonzero; otherwise PC assumes the initial guess
1836:    is to be zero (and thus zeros it out before solving).

1838:    Logically Collective on PC

1840:    Input Parameter:
1841: .   pc - iterative context obtained from PCCreate()

1843:    Output Parameter:
1844: .  flg - PETSC_TRUE indicates the guess is non-zero, PETSC_FALSE indicates the guess is zero

1846:    Level: Developer

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

1850: .seealso: PCGetInitialGuessNonzero(), PCSetInitialGuessKnoll(), PCGetInitialGuessKnoll(), PCSetInitialGuessNonzero()
1851: @*/
1852: PetscErrorCode  PCGetInitialGuessNonzero(PC pc,PetscBool *flg)
1853: {
1855:   *flg = pc->nonzero_guess;
1856:   return(0);
1857: }

1861: /*@C
1862:   PCRegister -  Adds a method to the preconditioner package.

1864:    Not collective

1866:    Input Parameters:
1867: +  name_solver - name of a new user-defined solver
1868: -  routine_create - routine to create method context

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

1873:    Sample usage:
1874: .vb
1875:    PCRegister("my_solver", MySolverCreate);
1876: .ve

1878:    Then, your solver can be chosen with the procedural interface via
1879: $     PCSetType(pc,"my_solver")
1880:    or at runtime via the option
1881: $     -pc_type my_solver

1883:    Level: advanced

1885: .keywords: PC, register

1887: .seealso: PCRegisterAll(), PCRegisterDestroy()
1888: @*/
1889: PetscErrorCode  PCRegister(const char sname[],PetscErrorCode (*function)(PC))
1890: {

1894:   PetscFunctionListAdd(&PCList,sname,function);
1895:   return(0);
1896: }

1900: /*@
1901:     PCComputeExplicitOperator - Computes the explicit preconditioned operator.

1903:     Collective on PC

1905:     Input Parameter:
1906: .   pc - the preconditioner object

1908:     Output Parameter:
1909: .   mat - the explict preconditioned operator

1911:     Notes:
1912:     This computation is done by applying the operators to columns of the
1913:     identity matrix.

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

1919:     Level: advanced

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

1923: .seealso: KSPComputeExplicitOperator()

1925: @*/
1926: PetscErrorCode  PCComputeExplicitOperator(PC pc,Mat *mat)
1927: {
1928:   Vec            in,out;
1930:   PetscInt       i,M,m,*rows,start,end;
1931:   PetscMPIInt    size;
1932:   MPI_Comm       comm;
1933:   PetscScalar    *array,one = 1.0;


1939:   PetscObjectGetComm((PetscObject)pc,&comm);
1940:   MPI_Comm_size(comm,&size);

1942:   if (!pc->pmat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"You must call KSPSetOperators() or PCSetOperators() before this call");
1943:   MatCreateVecs(pc->pmat,&in,0);
1944:   VecDuplicate(in,&out);
1945:   VecGetOwnershipRange(in,&start,&end);
1946:   VecGetSize(in,&M);
1947:   VecGetLocalSize(in,&m);
1948:   PetscMalloc1(m+1,&rows);
1949:   for (i=0; i<m; i++) rows[i] = start + i;

1951:   MatCreate(comm,mat);
1952:   MatSetSizes(*mat,m,m,M,M);
1953:   if (size == 1) {
1954:     MatSetType(*mat,MATSEQDENSE);
1955:     MatSeqDenseSetPreallocation(*mat,NULL);
1956:   } else {
1957:     MatSetType(*mat,MATMPIAIJ);
1958:     MatMPIAIJSetPreallocation(*mat,0,NULL,0,NULL);
1959:   }
1960:   MatSetOption(*mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);

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

1964:     VecSet(in,0.0);
1965:     VecSetValues(in,1,&i,&one,INSERT_VALUES);
1966:     VecAssemblyBegin(in);
1967:     VecAssemblyEnd(in);

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

1972:     VecGetArray(out,&array);
1973:     MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
1974:     VecRestoreArray(out,&array);

1976:   }
1977:   PetscFree(rows);
1978:   VecDestroy(&out);
1979:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
1980:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
1981:   return(0);
1982: }

1986: /*@
1987:    PCSetCoordinates - sets the coordinates of all the nodes on the local process

1989:    Collective on PC

1991:    Input Parameters:
1992: +  pc - the solver context
1993: .  dim - the dimension of the coordinates 1, 2, or 3
1994: -  coords - the coordinates

1996:    Level: intermediate

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

2006: .seealso: MatSetNearNullSpace
2007: @*/
2008: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
2009: {

2013:   PetscTryMethod(pc,"PCSetCoordinates_C",(PC,PetscInt,PetscInt,PetscReal*),(pc,dim,nloc,coords));
2014:   return(0);
2015: }