Actual source code: precon.c

petsc-3.5.3 2015-01-31
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:    PCGetUseAmat - Gets a flag to indicate that when the preconditioner needs to apply (part of) the
324:    operator during the preconditioning process it applies the Amat provided to TSSetRHSJacobian(),
325:    TSSetIJacobian(), SNESSetJacobian(), KSPSetOperator() or PCSetOperator() not the Pmat.

327:    Logically Collective on PC

329:    Input Parameter:
330: .  pc - the preconditioner context

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

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

339:    Level: intermediate

341: .seealso: PCSetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
342: @*/
343: PetscErrorCode  PCGetUseAmat(PC pc,PetscBool *flg)
344: {
347:   *flg = pc->useAmat;
348:   return(0);
349: }

353: /*@
354:    PCCreate - Creates a preconditioner context.

356:    Collective on MPI_Comm

358:    Input Parameter:
359: .  comm - MPI communicator

361:    Output Parameter:
362: .  pc - location to put the preconditioner context

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

368:    Level: developer

370: .keywords: PC, create, context

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

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

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

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

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

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

401: }

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

407: /*@
408:    PCApply - Applies the preconditioner to a vector.

410:    Collective on PC and Vec

412:    Input Parameters:
413: +  pc - the preconditioner context
414: -  x - input vector

416:    Output Parameter:
417: .  y - output vector

419:    Level: developer

421: .keywords: PC, apply

423: .seealso: PCApplyTranspose(), PCApplyBAorAB()
424: @*/
425: PetscErrorCode  PCApply(PC pc,Vec x,Vec y)
426: {

433:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
434:   VecValidValues(x,2,PETSC_TRUE);
435:   if (pc->setupcalled < 2) {
436:     PCSetUp(pc);
437:   }
438:   if (!pc->ops->apply) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply");
439:   PetscLogEventBegin(PC_Apply,pc,x,y,0);
440:   (*pc->ops->apply)(pc,x,y);
441:   PetscLogEventEnd(PC_Apply,pc,x,y,0);
442:   VecValidValues(y,3,PETSC_FALSE);
443:   return(0);
444: }

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

451:    Collective on PC and Vec

453:    Input Parameters:
454: +  pc - the preconditioner context
455: -  x - input vector

457:    Output Parameter:
458: .  y - output vector

460:    Notes:
461:    Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.

463:    Level: developer

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

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

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

492: /*@
493:    PCApplySymmetricRight - Applies the right 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:    Level: developer

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

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

511: .seealso: PCApply(), PCApplySymmetricLeft()
512: @*/
513: PetscErrorCode  PCApplySymmetricRight(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:   VecValidValues(x,2,PETSC_TRUE);
523:   if (pc->setupcalled < 2) {
524:     PCSetUp(pc);
525:   }
526:   if (!pc->ops->applysymmetricright) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have left symmetric apply");
527:   PetscLogEventBegin(PC_ApplySymmetricRight,pc,x,y,0);
528:   (*pc->ops->applysymmetricright)(pc,x,y);
529:   PetscLogEventEnd(PC_ApplySymmetricRight,pc,x,y,0);
530:   VecValidValues(y,3,PETSC_FALSE);
531:   return(0);
532: }

536: /*@
537:    PCApplyTranspose - Applies the transpose of preconditioner to a vector.

539:    Collective on PC and Vec

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

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

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

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

552:    Level: developer

554: .keywords: PC, apply, transpose

556: .seealso: PCApply(), PCApplyBAorAB(), PCApplyBAorABTranspose(), PCApplyTransposeExists()
557: @*/
558: PetscErrorCode  PCApplyTranspose(PC pc,Vec x,Vec y)
559: {

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

581: /*@
582:    PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation

584:    Collective on PC and Vec

586:    Input Parameters:
587: .  pc - the preconditioner context

589:    Output Parameter:
590: .  flg - PETSC_TRUE if a transpose operation is defined

592:    Level: developer

594: .keywords: PC, apply, transpose

596: .seealso: PCApplyTranspose()
597: @*/
598: PetscErrorCode  PCApplyTransposeExists(PC pc,PetscBool  *flg)
599: {
603:   if (pc->ops->applytranspose) *flg = PETSC_TRUE;
604:   else *flg = PETSC_FALSE;
605:   return(0);
606: }

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

613:    Collective on PC and Vec

615:    Input Parameters:
616: +  pc - the preconditioner context
617: .  side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
618: .  x - input vector
619: -  work - work vector

621:    Output Parameter:
622: .  y - output vector

624:    Level: developer

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

629: .keywords: PC, apply, operator

631: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorABTranspose()
632: @*/
633: PetscErrorCode  PCApplyBAorAB(PC pc,PCSide side,Vec x,Vec y,Vec work)
634: {

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

647:   if (pc->setupcalled < 2) {
648:     PCSetUp(pc);
649:   }

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

693: /*@
694:    PCApplyBAorABTranspose - Applies the transpose of the preconditioner
695:    and operator to a vector. That is, applies tr(B) * tr(A) with left preconditioning,
696:    NOT tr(B*A) = tr(A)*tr(B).

698:    Collective on PC and Vec

700:    Input Parameters:
701: +  pc - the preconditioner context
702: .  side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
703: .  x - input vector
704: -  work - work vector

706:    Output Parameter:
707: .  y - output vector


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

713:     Level: developer

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

717: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorAB()
718: @*/
719: PetscErrorCode  PCApplyBAorABTranspose(PC pc,PCSide side,Vec x,Vec y,Vec work)
720: {

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

737:   if (pc->setupcalled < 2) {
738:     PCSetUp(pc);
739:   }

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

753: /* -------------------------------------------------------------------------------*/

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

761:    Not Collective

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

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

769:    Level: developer

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

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

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

792:    Collective on PC

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

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

809:    Notes:
810:    Most preconditioners do not support this function. Use the command
811:    PCApplyRichardsonExists() to determine if one does.

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

816:    Level: developer

818: .keywords: PC, apply, Richardson

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

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

840: /*
841:       a setupcall of 0 indicates never setup,
842:                      1 indicates has been previously setup
843: */
846: /*@
847:    PCSetUp - Prepares for the use of a preconditioner.

849:    Collective on PC

851:    Input Parameter:
852: .  pc - the preconditioner context

854:    Level: developer

856: .keywords: PC, setup

858: .seealso: PCCreate(), PCApply(), PCDestroy()
859: @*/
860: PetscErrorCode  PCSetUp(PC pc)
861: {
862:   PetscErrorCode   ierr;
863:   const char       *def;
864:   PetscObjectState matstate, matnonzerostate;

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

870:   if (pc->setupcalled && pc->reusepreconditioner) {
871:     PetscInfo(pc,"Leaving PC with identical preconditioner since reuse preconditioner is set\n");
872:     return(0);
873:   }

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

895:   if (!((PetscObject)pc)->type_name) {
896:     PCGetDefaultType_Private(pc,&def);
897:     PCSetType(pc,def);
898:   }

900:   PetscLogEventBegin(PC_SetUp,pc,0,0,0);
901:   if (pc->ops->setup) {
902:     (*pc->ops->setup)(pc);
903:   }
904:   PetscLogEventEnd(PC_SetUp,pc,0,0,0);
905:   pc->setupcalled = 1;
906:   return(0);
907: }

911: /*@
912:    PCSetUpOnBlocks - Sets up the preconditioner for each block in
913:    the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
914:    methods.

916:    Collective on PC

918:    Input Parameters:
919: .  pc - the preconditioner context

921:    Level: developer

923: .keywords: PC, setup, blocks

925: .seealso: PCCreate(), PCApply(), PCDestroy(), PCSetUp()
926: @*/
927: PetscErrorCode  PCSetUpOnBlocks(PC pc)
928: {

933:   if (!pc->ops->setuponblocks) return(0);
934:   PetscLogEventBegin(PC_SetUpOnBlocks,pc,0,0,0);
935:   (*pc->ops->setuponblocks)(pc);
936:   PetscLogEventEnd(PC_SetUpOnBlocks,pc,0,0,0);
937:   return(0);
938: }

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

949:    Logically Collective on PC

951:    Input Parameters:
952: +  pc - the preconditioner context
953: .  func - routine for modifying the submatrices
954: -  ctx - optional user-defined context (may be null)

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

959: .  row - an array of index sets that contain the global row numbers
960:          that comprise each local submatrix
961: .  col - an array of index sets that contain the global column numbers
962:          that comprise each local submatrix
963: .  submat - array of local submatrices
964: -  ctx - optional user-defined context for private data for the
965:          user-defined func routine (may be null)

967:    Notes:
968:    PCSetModifySubMatrices() MUST be called before KSPSetUp() and
969:    KSPSolve().

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

975:    Level: advanced

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

979: .seealso: PCModifySubMatrices(), PCASMGetSubMatrices()
980: @*/
981: PetscErrorCode  PCSetModifySubMatrices(PC pc,PetscErrorCode (*func)(PC,PetscInt,const IS[],const IS[],Mat[],void*),void *ctx)
982: {
985:   pc->modifysubmatrices  = func;
986:   pc->modifysubmatricesP = ctx;
987:   return(0);
988: }

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

996:    Collective on PC

998:    Input Parameters:
999: +  pc - the preconditioner context
1000: .  nsub - the number of local submatrices
1001: .  row - an array of index sets that contain the global row numbers
1002:          that comprise each local submatrix
1003: .  col - an array of index sets that contain the global column numbers
1004:          that comprise each local submatrix
1005: .  submat - array of local submatrices
1006: -  ctx - optional user-defined context for private data for the
1007:          user-defined routine (may be null)

1009:    Output Parameter:
1010: .  submat - array of local submatrices (the entries of which may
1011:             have been modified)

1013:    Notes:
1014:    The user should NOT generally call this routine, as it will
1015:    automatically be called within certain preconditioners (currently
1016:    block Jacobi, additive Schwarz) if set.

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

1023:    Level: developer

1025: .keywords: PC, modify, submatrices

1027: .seealso: PCSetModifySubMatrices()
1028: @*/
1029: PetscErrorCode  PCModifySubMatrices(PC pc,PetscInt nsub,const IS row[],const IS col[],Mat submat[],void *ctx)
1030: {

1035:   if (!pc->modifysubmatrices) return(0);
1036:   PetscLogEventBegin(PC_ModifySubMatrices,pc,0,0,0);
1037:   (*pc->modifysubmatrices)(pc,nsub,row,col,submat,ctx);
1038:   PetscLogEventEnd(PC_ModifySubMatrices,pc,0,0,0);
1039:   return(0);
1040: }

1044: /*@
1045:    PCSetOperators - Sets the matrix associated with the linear system and
1046:    a (possibly) different one associated with the preconditioner.

1048:    Logically Collective on PC and Mat

1050:    Input Parameters:
1051: +  pc - the preconditioner context
1052: .  Amat - the matrix that defines the linear system
1053: -  Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.

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

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

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

1068:    Level: intermediate

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

1072: .seealso: PCGetOperators(), MatZeroEntries()
1073:  @*/
1074: PetscErrorCode  PCSetOperators(PC pc,Mat Amat,Mat Pmat)
1075: {
1076:   PetscErrorCode   ierr;
1077:   PetscInt         m1,n1,m2,n2;

1085:   if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
1086:     MatGetLocalSize(Amat,&m1,&n1);
1087:     MatGetLocalSize(pc->mat,&m2,&n2);
1088:     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);
1089:     MatGetLocalSize(Pmat,&m1,&n1);
1090:     MatGetLocalSize(pc->pmat,&m2,&n2);
1091:     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);
1092:   }

1094:   if (Pmat != pc->pmat) {
1095:     /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
1096:     pc->matnonzerostate = -1;
1097:     pc->matstate        = -1;
1098:   }

1100:   /* reference first in case the matrices are the same */
1101:   if (Amat) {PetscObjectReference((PetscObject)Amat);}
1102:   MatDestroy(&pc->mat);
1103:   if (Pmat) {PetscObjectReference((PetscObject)Pmat);}
1104:   MatDestroy(&pc->pmat);
1105:   pc->mat  = Amat;
1106:   pc->pmat = Pmat;
1107:   return(0);
1108: }

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

1115:    Logically Collective on PC

1117:    Input Parameters:
1118: +  pc - the preconditioner context
1119: -  flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner

1121: .seealso: PCGetOperators(), MatZeroEntries(), PCGetReusePreconditioner(), KSPSetReusePreconditioner()
1122:  @*/
1123: PetscErrorCode  PCSetReusePreconditioner(PC pc,PetscBool flag)
1124: {
1127:   pc->reusepreconditioner = flag;
1128:   return(0);
1129: }

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

1136:    Not Collective

1138:    Input Parameter:
1139: .  pc - the preconditioner context

1141:    Output Parameter:
1142: .  flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner

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

1156: /*@C
1157:    PCGetOperators - Gets the matrix associated with the linear system and
1158:    possibly a different one associated with the preconditioner.

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

1162:    Input Parameter:
1163: .  pc - the preconditioner context

1165:    Output Parameters:
1166: +  Amat - the matrix defining the linear system
1167: -  Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.

1169:    Level: intermediate

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

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

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

1184: $         MatCreate(comm,&mat);
1185: $         KSP/PCSetOperators(ksp/pc,Amat,Amat);
1186: $         PetscObjectDereference((PetscObject)mat);
1187: $           set size, type, etc of Amat

1189:      and

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

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

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


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

1213: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperatorsSet()
1214: @*/
1215: PetscErrorCode  PCGetOperators(PC pc,Mat *Amat,Mat *Pmat)
1216: {

1221:   if (Amat) {
1222:     if (!pc->mat) {
1223:       if (pc->pmat && !Pmat) {  /* Apmat has been set, but user did not request it, so use for Amat */
1224:         pc->mat = pc->pmat;
1225:         PetscObjectReference((PetscObject)pc->mat);
1226:       } else {                  /* both Amat and Pmat are empty */
1227:         MatCreate(PetscObjectComm((PetscObject)pc),&pc->mat);
1228:         if (!Pmat) { /* user did NOT request Pmat, so make same as Amat */
1229:           pc->pmat = pc->mat;
1230:           PetscObjectReference((PetscObject)pc->pmat);
1231:         }
1232:       }
1233:     }
1234:     *Amat = pc->mat;
1235:   }
1236:   if (Pmat) {
1237:     if (!pc->pmat) {
1238:       if (pc->mat && !Amat) {    /* Amat has been set but was not requested, so use for pmat */
1239:         pc->pmat = pc->mat;
1240:         PetscObjectReference((PetscObject)pc->pmat);
1241:       } else {
1242:         MatCreate(PetscObjectComm((PetscObject)pc),&pc->pmat);
1243:         if (!Amat) { /* user did NOT request Amat, so make same as Pmat */
1244:           pc->mat = pc->pmat;
1245:           PetscObjectReference((PetscObject)pc->mat);
1246:         }
1247:       }
1248:     }
1249:     *Pmat = pc->pmat;
1250:   }
1251:   return(0);
1252: }

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

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

1262:    Input Parameter:
1263: .  pc - the preconditioner context

1265:    Output Parameters:
1266: +  mat - the matrix associated with the linear system was set
1267: -  pmat - matrix associated with the preconditioner was set, usually the same

1269:    Level: intermediate

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

1273: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperators()
1274: @*/
1275: PetscErrorCode  PCGetOperatorsSet(PC pc,PetscBool  *mat,PetscBool  *pmat)
1276: {
1279:   if (mat) *mat = (pc->mat)  ? PETSC_TRUE : PETSC_FALSE;
1280:   if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1281:   return(0);
1282: }

1286: /*@
1287:    PCFactorGetMatrix - Gets the factored matrix from the
1288:    preconditioner context.  This routine is valid only for the LU,
1289:    incomplete LU, Cholesky, and incomplete Cholesky methods.

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

1293:    Input Parameters:
1294: .  pc - the preconditioner context

1296:    Output parameters:
1297: .  mat - the factored matrix

1299:    Level: advanced

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

1303: .keywords: PC, get, factored, matrix
1304: @*/
1305: PetscErrorCode  PCFactorGetMatrix(PC pc,Mat *mat)
1306: {

1312:   if (pc->ops->getfactoredmatrix) {
1313:     (*pc->ops->getfactoredmatrix)(pc,mat);
1314:   } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC type does not support getting factor matrix");
1315:   return(0);
1316: }

1320: /*@C
1321:    PCSetOptionsPrefix - Sets the prefix used for searching for all
1322:    PC options in the database.

1324:    Logically Collective on PC

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

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

1335:    Level: advanced

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

1339: .seealso: PCAppendOptionsPrefix(), PCGetOptionsPrefix()
1340: @*/
1341: PetscErrorCode  PCSetOptionsPrefix(PC pc,const char prefix[])
1342: {

1347:   PetscObjectSetOptionsPrefix((PetscObject)pc,prefix);
1348:   return(0);
1349: }

1353: /*@C
1354:    PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1355:    PC options in the database.

1357:    Logically Collective on PC

1359:    Input Parameters:
1360: +  pc - the preconditioner context
1361: -  prefix - the prefix string to prepend to all PC option requests

1363:    Notes:
1364:    A hyphen (-) must NOT be given at the beginning of the prefix name.
1365:    The first character of all runtime options is AUTOMATICALLY the
1366:    hyphen.

1368:    Level: advanced

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

1372: .seealso: PCSetOptionsPrefix(), PCGetOptionsPrefix()
1373: @*/
1374: PetscErrorCode  PCAppendOptionsPrefix(PC pc,const char prefix[])
1375: {

1380:   PetscObjectAppendOptionsPrefix((PetscObject)pc,prefix);
1381:   return(0);
1382: }

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

1390:    Not Collective

1392:    Input Parameters:
1393: .  pc - the preconditioner context

1395:    Output Parameters:
1396: .  prefix - pointer to the prefix string used, is returned

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

1401:    Level: advanced

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

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

1414:   PetscObjectGetOptionsPrefix((PetscObject)pc,prefix);
1415:   return(0);
1416: }

1420: /*@
1421:    PCPreSolve - Optional pre-solve phase, intended for any
1422:    preconditioner-specific actions that must be performed before
1423:    the iterative solve itself.

1425:    Collective on PC

1427:    Input Parameters:
1428: +  pc - the preconditioner context
1429: -  ksp - the Krylov subspace context

1431:    Level: developer

1433:    Sample of Usage:
1434: .vb
1435:     PCPreSolve(pc,ksp);
1436:     KSPSolve(ksp,b,x);
1437:     PCPostSolve(pc,ksp);
1438: .ve

1440:    Notes:
1441:    The pre-solve phase is distinct from the PCSetUp() phase.

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

1445: .keywords: PC, pre-solve

1447: .seealso: PCPostSolve()
1448: @*/
1449: PetscErrorCode  PCPreSolve(PC pc,KSP ksp)
1450: {
1452:   Vec            x,rhs;

1457:   pc->presolvedone++;
1458:   if (pc->presolvedone > 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot embed PCPreSolve() more than twice");
1459:   KSPGetSolution(ksp,&x);
1460:   KSPGetRhs(ksp,&rhs);

1462:   if (pc->ops->presolve) {
1463:     (*pc->ops->presolve)(pc,ksp,rhs,x);
1464:   }
1465:   return(0);
1466: }

1470: /*@
1471:    PCPostSolve - Optional post-solve phase, intended for any
1472:    preconditioner-specific actions that must be performed after
1473:    the iterative solve itself.

1475:    Collective on PC

1477:    Input Parameters:
1478: +  pc - the preconditioner context
1479: -  ksp - the Krylov subspace context

1481:    Sample of Usage:
1482: .vb
1483:     PCPreSolve(pc,ksp);
1484:     KSPSolve(ksp,b,x);
1485:     PCPostSolve(pc,ksp);
1486: .ve

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

1491:    Level: developer

1493: .keywords: PC, post-solve

1495: .seealso: PCPreSolve(), KSPSolve()
1496: @*/
1497: PetscErrorCode  PCPostSolve(PC pc,KSP ksp)
1498: {
1500:   Vec            x,rhs;

1505:   pc->presolvedone--;
1506:   KSPGetSolution(ksp,&x);
1507:   KSPGetRhs(ksp,&rhs);
1508:   if (pc->ops->postsolve) {
1509:      (*pc->ops->postsolve)(pc,ksp,rhs,x);
1510:   }
1511:   return(0);
1512: }

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

1519:   Collective on PetscViewer

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

1526:    Level: intermediate

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

1531:   Notes for advanced users:
1532:   Most users should not need to know the details of the binary storage
1533:   format, since PCLoad() and PCView() completely hide these details.
1534:   But for anyone who's interested, the standard binary matrix storage
1535:   format is
1536: .vb
1537:      has not yet been determined
1538: .ve

1540: .seealso: PetscViewerBinaryOpen(), PCView(), MatLoad(), VecLoad()
1541: @*/
1542: PetscErrorCode  PCLoad(PC newdm, PetscViewer viewer)
1543: {
1545:   PetscBool      isbinary;
1546:   PetscInt       classid;
1547:   char           type[256];

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

1555:   PetscViewerBinaryRead(viewer,&classid,1,PETSC_INT);
1556:   if (classid != PC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not PC next in file");
1557:   PetscViewerBinaryRead(viewer,type,256,PETSC_CHAR);
1558:   PCSetType(newdm, type);
1559:   if (newdm->ops->load) {
1560:     (*newdm->ops->load)(newdm,viewer);
1561:   }
1562:   return(0);
1563: }

1565: #include <petscdraw.h>
1566: #if defined(PETSC_HAVE_SAWS)
1567: #include <petscviewersaws.h>
1568: #endif
1571: /*@C
1572:    PCView - Prints the PC data structure.

1574:    Collective on PC

1576:    Input Parameters:
1577: +  PC - the PC context
1578: -  viewer - optional visualization context

1580:    Note:
1581:    The available visualization contexts include
1582: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
1583: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1584:          output where only the first processor opens
1585:          the file.  All other processors send their
1586:          data to the first processor to print.

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

1591:    Level: developer

1593: .keywords: PC, view

1595: .seealso: KSPView(), PetscViewerASCIIOpen()
1596: @*/
1597: PetscErrorCode  PCView(PC pc,PetscViewer viewer)
1598: {
1599:   PCType            cstr;
1600:   PetscErrorCode    ierr;
1601:   PetscBool         iascii,isstring,isbinary,isdraw;
1602:   PetscViewerFormat format;
1603: #if defined(PETSC_HAVE_SAWS)
1604:   PetscBool         isams;
1605: #endif

1609:   if (!viewer) {
1610:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);
1611:   }

1615:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1616:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
1617:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1618:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1619: #if defined(PETSC_HAVE_SAWS)
1620:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&isams);
1621: #endif

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

1664:     PetscObjectGetComm((PetscObject)pc,&comm);
1665:     MPI_Comm_rank(comm,&rank);
1666:     if (!rank) {
1667:       PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
1668:       PetscStrncpy(type,((PetscObject)pc)->type_name,256);
1669:       PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);
1670:     }
1671:     if (pc->ops->view) {
1672:       (*pc->ops->view)(pc,viewer);
1673:     }
1674:   } else if (isdraw) {
1675:     PetscDraw draw;
1676:     char      str[25];
1677:     PetscReal x,y,bottom,h;
1678:     PetscInt  n;

1680:     PetscViewerDrawGetDraw(viewer,0,&draw);
1681:     PetscDrawGetCurrentPoint(draw,&x,&y);
1682:     if (pc->mat) {
1683:       MatGetSize(pc->mat,&n,NULL);
1684:       PetscSNPrintf(str,25,"PC: %s (%D)",((PetscObject)pc)->type_name,n);
1685:     } else {
1686:       PetscSNPrintf(str,25,"PC: %s",((PetscObject)pc)->type_name);
1687:     }
1688:     PetscDrawBoxedString(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
1689:     bottom = y - h;
1690:     PetscDrawPushCurrentPoint(draw,x,bottom);
1691:     if (pc->ops->view) {
1692:       (*pc->ops->view)(pc,viewer);
1693:     }
1694:     PetscDrawPopCurrentPoint(draw);
1695: #if defined(PETSC_HAVE_SAWS)
1696:   } else if (isams) {
1697:     PetscMPIInt rank;

1699:     PetscObjectName((PetscObject)pc);
1700:     MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
1701:     if (!((PetscObject)pc)->amsmem && !rank) {
1702:       PetscObjectViewSAWs((PetscObject)pc,viewer);
1703:     }
1704:     if (pc->mat) {MatView(pc->mat,viewer);}
1705:     if (pc->pmat && pc->pmat != pc->mat) {MatView(pc->pmat,viewer);}
1706: #endif
1707:   }
1708:   return(0);
1709: }


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

1719:    Logically Collective on PC

1721:    Input Parameters:
1722: +  pc - iterative context obtained from PCCreate()
1723: -  flg - PETSC_TRUE indicates the guess is non-zero, PETSC_FALSE indicates the guess is zero

1725:    Level: Developer

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


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

1736: .seealso: PCGetInitialGuessNonzero(), PCSetInitialGuessKnoll(), PCGetInitialGuessKnoll()
1737: @*/
1738: PetscErrorCode  PCSetInitialGuessNonzero(PC pc,PetscBool flg)
1739: {
1742:   pc->nonzero_guess = flg;
1743:   return(0);
1744: }

1748: /*@C
1749:   PCRegister -  Adds a method to the preconditioner package.

1751:    Not collective

1753:    Input Parameters:
1754: +  name_solver - name of a new user-defined solver
1755: -  routine_create - routine to create method context

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

1760:    Sample usage:
1761: .vb
1762:    PCRegister("my_solver", MySolverCreate);
1763: .ve

1765:    Then, your solver can be chosen with the procedural interface via
1766: $     PCSetType(pc,"my_solver")
1767:    or at runtime via the option
1768: $     -pc_type my_solver

1770:    Level: advanced

1772: .keywords: PC, register

1774: .seealso: PCRegisterAll(), PCRegisterDestroy()
1775: @*/
1776: PetscErrorCode  PCRegister(const char sname[],PetscErrorCode (*function)(PC))
1777: {

1781:   PetscFunctionListAdd(&PCList,sname,function);
1782:   return(0);
1783: }

1787: /*@
1788:     PCComputeExplicitOperator - Computes the explicit preconditioned operator.

1790:     Collective on PC

1792:     Input Parameter:
1793: .   pc - the preconditioner object

1795:     Output Parameter:
1796: .   mat - the explict preconditioned operator

1798:     Notes:
1799:     This computation is done by applying the operators to columns of the
1800:     identity matrix.

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

1806:     Level: advanced

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

1810: .seealso: KSPComputeExplicitOperator()

1812: @*/
1813: PetscErrorCode  PCComputeExplicitOperator(PC pc,Mat *mat)
1814: {
1815:   Vec            in,out;
1817:   PetscInt       i,M,m,*rows,start,end;
1818:   PetscMPIInt    size;
1819:   MPI_Comm       comm;
1820:   PetscScalar    *array,one = 1.0;


1826:   PetscObjectGetComm((PetscObject)pc,&comm);
1827:   MPI_Comm_size(comm,&size);

1829:   if (!pc->pmat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"You must call KSPSetOperators() or PCSetOperators() before this call");
1830:   MatGetVecs(pc->pmat,&in,0);
1831:   VecDuplicate(in,&out);
1832:   VecGetOwnershipRange(in,&start,&end);
1833:   VecGetSize(in,&M);
1834:   VecGetLocalSize(in,&m);
1835:   PetscMalloc1((m+1),&rows);
1836:   for (i=0; i<m; i++) rows[i] = start + i;

1838:   MatCreate(comm,mat);
1839:   MatSetSizes(*mat,m,m,M,M);
1840:   if (size == 1) {
1841:     MatSetType(*mat,MATSEQDENSE);
1842:     MatSeqDenseSetPreallocation(*mat,NULL);
1843:   } else {
1844:     MatSetType(*mat,MATMPIAIJ);
1845:     MatMPIAIJSetPreallocation(*mat,0,NULL,0,NULL);
1846:   }
1847:   MatSetOption(*mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);

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

1851:     VecSet(in,0.0);
1852:     VecSetValues(in,1,&i,&one,INSERT_VALUES);
1853:     VecAssemblyBegin(in);
1854:     VecAssemblyEnd(in);

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

1859:     VecGetArray(out,&array);
1860:     MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
1861:     VecRestoreArray(out,&array);

1863:   }
1864:   PetscFree(rows);
1865:   VecDestroy(&out);
1866:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
1867:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
1868:   return(0);
1869: }

1873: /*@
1874:    PCSetCoordinates - sets the coordinates of all the nodes on the local process

1876:    Collective on PC

1878:    Input Parameters:
1879: +  pc - the solver context
1880: .  dim - the dimension of the coordinates 1, 2, or 3
1881: -  coords - the coordinates

1883:    Level: intermediate

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

1893: .seealso: MatSetNearNullSpace
1894: @*/
1895: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
1896: {

1900:   PetscTryMethod(pc,"PCSetCoordinates_C",(PC,PetscInt,PetscInt,PetscReal*),(pc,dim,nloc,coords));
1901:   return(0);
1902: }