Actual source code: precon.c

      1: #define PETSCKSP_DLL

      3: /*
      4:     The PC (preconditioner) interface routines, callable by users.
      5: */
     6:  #include src/ksp/pc/pcimpl.h

      8: /* Logging support */
      9: PetscCookie PC_COOKIE = 0;
     10: PetscEvent  PC_SetUp = 0, PC_SetUpOnBlocks = 0, PC_Apply = 0, PC_ApplyCoarse = 0, PC_ApplyMultiple = 0, PC_ApplySymmetricLeft = 0;
     11: PetscEvent  PC_ApplySymmetricRight = 0, PC_ModifySubMatrices = 0;

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

     22:   MPI_Comm_size(pc->comm,&size);
     23:   if (pc->pmat) {
     24:     PetscErrorCode (*f)(Mat,PetscTruth*,MatReuse,Mat*);
     25:     PetscObjectQueryFunction((PetscObject)pc->pmat,"MatGetDiagonalBlock_C",(void (**)(void))&f);
     26:     if (size == 1) {
     27:       MatHasOperation(pc->pmat,MATOP_ICCFACTOR_SYMBOLIC,&flg1);
     28:       MatHasOperation(pc->pmat,MATOP_ILUFACTOR_SYMBOLIC,&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: /*@C
     59:    PCDestroy - Destroys PC context that was created with PCCreate().

     61:    Collective on PC

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

     66:    Level: developer

     68: .keywords: PC, destroy

     70: .seealso: PCCreate(), PCSetUp()
     71: @*/
     72: PetscErrorCode PETSCKSP_DLLEXPORT PCDestroy(PC pc)
     73: {

     78:   if (--pc->refct > 0) return(0);

     80:   /* if memory was published with AMS then destroy it */
     81:   PetscObjectDepublish(pc);

     83:   if (pc->ops->destroy)       { (*pc->ops->destroy)(pc);}
     84:   if (pc->diagonalscaleright) {VecDestroy(pc->diagonalscaleright);}
     85:   if (pc->diagonalscaleleft)  {VecDestroy(pc->diagonalscaleleft);}

     87:   PetscHeaderDestroy(pc);
     88:   return(0);
     89: }

     93: /*@C
     94:    PCDiagonalScale - Indicates if the preconditioner applies an additional left and right
     95:       scaling as needed by certain time-stepping codes.

     97:    Collective on PC

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

    102:    Output Parameter:
    103: .  flag - PETSC_TRUE if it applies the scaling

    105:    Level: developer

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

    111: .keywords: PC

    113: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCDiagonalScaleSet()
    114: @*/
    115: PetscErrorCode PETSCKSP_DLLEXPORT PCDiagonalScale(PC pc,PetscTruth *flag)
    116: {
    120:   *flag = pc->diagonalscale;
    121:   return(0);
    122: }

    126: /*@
    127:    PCDiagonalScaleSet - Indicates the left scaling to use to apply an additional left and right
    128:       scaling as needed by certain time-stepping codes.

    130:    Collective on PC

    132:    Input Parameters:
    133: +  pc - the preconditioner context
    134: -  s - scaling vector

    136:    Level: intermediate

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

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

    144: .keywords: PC

    146: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCDiagonalScale()
    147: @*/
    148: PetscErrorCode PETSCKSP_DLLEXPORT PCDiagonalScaleSet(PC pc,Vec s)
    149: {

    155:   pc->diagonalscale     = PETSC_TRUE;
    156:   if (pc->diagonalscaleleft) {
    157:     VecDestroy(pc->diagonalscaleleft);
    158:   }
    159:   pc->diagonalscaleleft = s;
    160:   PetscObjectReference((PetscObject)s);
    161:   if (!pc->diagonalscaleright) {
    162:     VecDuplicate(s,&pc->diagonalscaleright);
    163:   }
    164:   VecCopy(s,pc->diagonalscaleright);
    165:   VecReciprocal(pc->diagonalscaleright);
    166:   return(0);
    167: }

    171: /*@C
    172:    PCDiagonalScaleLeft - Indicates the left scaling to use to apply an additional left and right
    173:       scaling as needed by certain time-stepping codes.

    175:    Collective on PC

    177:    Input Parameters:
    178: +  pc - the preconditioner context
    179: .  in - input vector
    180: +  out - scaled vector (maybe the same as in)

    182:    Level: intermediate

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

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

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

    192: .keywords: PC

    194: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleSet(), PCDiagonalScaleRight(), PCDiagonalScale()
    195: @*/
    196: PetscErrorCode PETSCKSP_DLLEXPORT PCDiagonalScaleLeft(PC pc,Vec in,Vec out)
    197: {

    204:   if (pc->diagonalscale) {
    205:     VecPointwiseMult(out,pc->diagonalscaleleft,in);
    206:   } else if (in != out) {
    207:     VecCopy(in,out);
    208:   }
    209:   return(0);
    210: }

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

    217:    Collective on PC

    219:    Input Parameters:
    220: +  pc - the preconditioner context
    221: .  in - input vector
    222: +  out - scaled vector (maybe the same as in)

    224:    Level: intermediate

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

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

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

    234: .keywords: PC

    236: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleSet(), PCDiagonalScale()
    237: @*/
    238: PetscErrorCode PETSCKSP_DLLEXPORT PCDiagonalScaleRight(PC pc,Vec in,Vec out)
    239: {

    246:   if (pc->diagonalscale) {
    247:     VecPointwiseMult(out,pc->diagonalscaleright,in);
    248:   } else if (in != out) {
    249:     VecCopy(in,out);
    250:   }
    251:   return(0);
    252: }

    256: static PetscErrorCode PCPublish_Petsc(PetscObject obj)
    257: {
    259:   return(0);
    260: }

    264: /*@C
    265:    PCCreate - Creates a preconditioner context.

    267:    Collective on MPI_Comm

    269:    Input Parameter:
    270: .  comm - MPI communicator

    272:    Output Parameter:
    273: .  pc - location to put the preconditioner context

    275:    Notes:
    276:    The default preconditioner on one processor is PCILU with 0 fill on more
    277:    then one it is PCBJACOBI with ILU() on each processor.

    279:    Level: developer

    281: .keywords: PC, create, context

    283: .seealso: PCSetUp(), PCApply(), PCDestroy()
    284: @*/
    285: PetscErrorCode PETSCKSP_DLLEXPORT PCCreate(MPI_Comm comm,PC *newpc)
    286: {
    287:   PC             pc;

    292:   *newpc = 0;
    293: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
    294:   PCInitializePackage(PETSC_NULL);
    295: #endif

    297:   PetscHeaderCreate(pc,_p_PC,struct _PCOps,PC_COOKIE,-1,"PC",comm,PCDestroy,PCView);
    298:   pc->bops->publish      = PCPublish_Petsc;
    299:   pc->mat                = 0;
    300:   pc->pmat               = 0;
    301:   pc->setupcalled        = 0;
    302:   pc->data               = 0;
    303:   pc->diagonalscale      = PETSC_FALSE;
    304:   pc->diagonalscaleleft  = 0;
    305:   pc->diagonalscaleright = 0;

    307:   pc->ops->destroy             = 0;
    308:   pc->ops->apply               = 0;
    309:   pc->ops->applytranspose      = 0;
    310:   pc->ops->applyBA             = 0;
    311:   pc->ops->applyBAtranspose    = 0;
    312:   pc->ops->applyrichardson     = 0;
    313:   pc->ops->view                = 0;
    314:   pc->ops->getfactoredmatrix   = 0;
    315:   pc->ops->applysymmetricright = 0;
    316:   pc->ops->applysymmetricleft  = 0;
    317:   pc->ops->setuponblocks       = 0;

    319:   pc->modifysubmatrices   = 0;
    320:   pc->modifysubmatricesP  = 0;
    321:   *newpc                  = pc;
    322:   PetscPublishAll(pc);
    323:   return(0);

    325: }

    327: /* -------------------------------------------------------------------------------*/

    331: /*@
    332:    PCApply - Applies the preconditioner to a vector.

    334:    Collective on PC and Vec

    336:    Input Parameters:
    337: +  pc - the preconditioner context
    338: -  x - input vector

    340:    Output Parameter:
    341: .  y - output vector

    343:    Level: developer

    345: .keywords: PC, apply

    347: .seealso: PCApplyTranspose(), PCApplyBAorAB()
    348: @*/
    349: PetscErrorCode PETSCKSP_DLLEXPORT PCApply(PC pc,Vec x,Vec y)
    350: {

    357:   if (x == y) SETERRQ(PETSC_ERR_ARG_IDN,"x and y must be different vectors");

    359:   if (pc->setupcalled < 2) {
    360:     PCSetUp(pc);
    361:   }
    362:   PetscLogEventBegin(PC_Apply,pc,x,y,0);
    363:   (*pc->ops->apply)(pc,x,y);
    364:   PetscLogEventEnd(PC_Apply,pc,x,y,0);
    365:   return(0);
    366: }

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

    373:    Collective on PC and Vec

    375:    Input Parameters:
    376: +  pc - the preconditioner context
    377: -  x - input vector

    379:    Output Parameter:
    380: .  y - output vector

    382:    Notes:
    383:    Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.

    385:    Level: developer

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

    389: .seealso: PCApply(), PCApplySymmetricRight()
    390: @*/
    391: PetscErrorCode PETSCKSP_DLLEXPORT PCApplySymmetricLeft(PC pc,Vec x,Vec y)
    392: {

    399:   if (!pc->ops->applysymmetricleft) SETERRQ(PETSC_ERR_SUP,"PC does not have left symmetric apply");

    401:   if (pc->setupcalled < 2) {
    402:     PCSetUp(pc);
    403:   }

    405:   PetscLogEventBegin(PC_ApplySymmetricLeft,pc,x,y,0);
    406:   (*pc->ops->applysymmetricleft)(pc,x,y);
    407:   PetscLogEventEnd(PC_ApplySymmetricLeft,pc,x,y,0);
    408:   return(0);
    409: }

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

    416:    Collective on PC and Vec

    418:    Input Parameters:
    419: +  pc - the preconditioner context
    420: -  x - input vector

    422:    Output Parameter:
    423: .  y - output vector

    425:    Level: developer

    427:    Notes:
    428:    Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.

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

    432: .seealso: PCApply(), PCApplySymmetricLeft()
    433: @*/
    434: PetscErrorCode PETSCKSP_DLLEXPORT PCApplySymmetricRight(PC pc,Vec x,Vec y)
    435: {

    442:   if (!pc->ops->applysymmetricright) SETERRQ(PETSC_ERR_SUP,"PC does not have left symmetric apply");

    444:   if (pc->setupcalled < 2) {
    445:     PCSetUp(pc);
    446:   }

    448:   PetscLogEventBegin(PC_ApplySymmetricRight,pc,x,y,0);
    449:   (*pc->ops->applysymmetricright)(pc,x,y);
    450:   PetscLogEventEnd(PC_ApplySymmetricRight,pc,x,y,0);
    451:   return(0);
    452: }

    456: /*@
    457:    PCApplyTranspose - Applies the transpose of preconditioner to a vector.

    459:    Collective on PC and Vec

    461:    Input Parameters:
    462: +  pc - the preconditioner context
    463: -  x - input vector

    465:    Output Parameter:
    466: .  y - output vector

    468:    Level: developer

    470: .keywords: PC, apply, transpose

    472: .seealso: PCApply(), PCApplyBAorAB(), PCApplyBAorABTranspose(), PCHasApplyTranspose()
    473: @*/
    474: PetscErrorCode PETSCKSP_DLLEXPORT PCApplyTranspose(PC pc,Vec x,Vec y)
    475: {

    482:   if (x == y) SETERRQ(PETSC_ERR_ARG_IDN,"x and y must be different vectors");
    483:   if (!pc->ops->applytranspose) SETERRQ(PETSC_ERR_SUP," ");

    485:   if (pc->setupcalled < 2) {
    486:     PCSetUp(pc);
    487:   }

    489:   PetscLogEventBegin(PC_Apply,pc,x,y,0);
    490:   (*pc->ops->applytranspose)(pc,x,y);
    491:   PetscLogEventEnd(PC_Apply,pc,x,y,0);
    492:   return(0);
    493: }

    497: /*@
    498:    PCHasApplyTranspose - Test whether the preconditioner has a transpose apply operation

    500:    Collective on PC and Vec

    502:    Input Parameters:
    503: .  pc - the preconditioner context

    505:    Output Parameter:
    506: .  flg - PETSC_TRUE if a transpose operation is defined

    508:    Level: developer

    510: .keywords: PC, apply, transpose

    512: .seealso: PCApplyTranspose()
    513: @*/
    514: PetscErrorCode PETSCKSP_DLLEXPORT PCHasApplyTranspose(PC pc,PetscTruth *flg)
    515: {
    519:   *flg = (PetscTruth) (pc->ops->applytranspose == 0);
    520:   return(0);
    521: }

    525: /*@
    526:    PCApplyBAorAB - Applies the preconditioner and operator to a vector.

    528:    Collective on PC and Vec

    530:    Input Parameters:
    531: +  pc - the preconditioner context
    532: .  side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
    533: .  x - input vector
    534: -  work - work vector

    536:    Output Parameter:
    537: .  y - output vector

    539:    Level: developer

    541: .keywords: PC, apply, operator

    543: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorABTranspose()
    544: @*/
    545: PetscErrorCode PETSCKSP_DLLEXPORT PCApplyBAorAB(PC pc,PCSide side,Vec x,Vec y,Vec work)
    546: {

    554:   if (x == y) SETERRQ(PETSC_ERR_ARG_IDN,"x and y must be different vectors");
    555:   if (side != PC_LEFT && side != PC_SYMMETRIC && side != PC_RIGHT) {
    556:     SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Side must be right, left, or symmetric");
    557:   }
    558:   if (pc->diagonalscale && side == PC_SYMMETRIC) {
    559:     SETERRQ(PETSC_ERR_SUP,"Cannot include diagonal scaling with symmetric preconditioner application");
    560:   }

    562:   if (pc->setupcalled < 2) {
    563:     PCSetUp(pc);
    564:   }

    566:   if (pc->diagonalscale) {
    567:     if (pc->ops->applyBA) {
    568:       Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
    569:       VecDuplicate(x,&work2);
    570:       PCDiagonalScaleRight(pc,x,work2);
    571:       (*pc->ops->applyBA)(pc,side,work2,y,work);
    572:       PCDiagonalScaleLeft(pc,y,y);
    573:       VecDestroy(work2);
    574:     } else if (side == PC_RIGHT) {
    575:       PCDiagonalScaleRight(pc,x,y);
    576:       PCApply(pc,y,work);
    577:       MatMult(pc->mat,work,y);
    578:       PCDiagonalScaleLeft(pc,y,y);
    579:     } else if (side == PC_LEFT) {
    580:       PCDiagonalScaleRight(pc,x,y);
    581:       MatMult(pc->mat,y,work);
    582:       PCApply(pc,work,y);
    583:       PCDiagonalScaleLeft(pc,y,y);
    584:     } else if (side == PC_SYMMETRIC) {
    585:       SETERRQ(PETSC_ERR_SUP,"Cannot provide diagonal scaling with symmetric application of preconditioner");
    586:     }
    587:   } else {
    588:     if (pc->ops->applyBA) {
    589:       (*pc->ops->applyBA)(pc,side,x,y,work);
    590:     } else if (side == PC_RIGHT) {
    591:       PCApply(pc,x,work);
    592:       MatMult(pc->mat,work,y);
    593:     } else if (side == PC_LEFT) {
    594:       MatMult(pc->mat,x,work);
    595:       PCApply(pc,work,y);
    596:     } else if (side == PC_SYMMETRIC) {
    597:       /* There's an extra copy here; maybe should provide 2 work vectors instead? */
    598:       PCApplySymmetricRight(pc,x,work);
    599:       MatMult(pc->mat,work,y);
    600:       VecCopy(y,work);
    601:       PCApplySymmetricLeft(pc,work,y);
    602:     }
    603:   }
    604:   return(0);
    605: }

    609: /*@
    610:    PCApplyBAorABTranspose - Applies the transpose of the preconditioner
    611:    and operator to a vector. That is, applies tr(B) * tr(A) with left preconditioning,
    612:    not tr(B*A) = tr(A)*tr(B).

    614:    Collective on PC and Vec

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

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

    625:    Level: developer

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

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

    640:   if (x == y) SETERRQ(PETSC_ERR_ARG_IDN,"x and y must be different vectors");
    641:   if (pc->ops->applyBAtranspose) {
    642:     (*pc->ops->applyBAtranspose)(pc,side,x,y,work);
    643:     return(0);
    644:   }
    645:   if (side != PC_LEFT && side != PC_RIGHT) {
    646:     SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Side must be right or left");
    647:   }

    649:   if (pc->setupcalled < 2) {
    650:     PCSetUp(pc);
    651:   }

    653:   if (side == PC_RIGHT) {
    654:     PCApplyTranspose(pc,x,work);
    655:     MatMultTranspose(pc->mat,work,y);
    656:   } else if (side == PC_LEFT) {
    657:     MatMultTranspose(pc->mat,x,work);
    658:     PCApplyTranspose(pc,work,y);
    659:   }
    660:   /* add support for PC_SYMMETRIC */
    661:   return(0); /* actually will never get here */
    662: }

    664: /* -------------------------------------------------------------------------------*/

    668: /*@
    669:    PCApplyRichardsonExists - Determines whether a particular preconditioner has a
    670:    built-in fast application of Richardson's method.

    672:    Not Collective

    674:    Input Parameter:
    675: .  pc - the preconditioner

    677:    Output Parameter:
    678: .  exists - PETSC_TRUE or PETSC_FALSE

    680:    Level: developer

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

    684: .seealso: PCApplyRichardson()
    685: @*/
    686: PetscErrorCode PETSCKSP_DLLEXPORT PCApplyRichardsonExists(PC pc,PetscTruth *exists)
    687: {
    691:   if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
    692:   else                    *exists = PETSC_FALSE;
    693:   return(0);
    694: }

    698: /*@
    699:    PCApplyRichardson - Applies several steps of Richardson iteration with
    700:    the particular preconditioner. This routine is usually used by the
    701:    Krylov solvers and not the application code directly.

    703:    Collective on PC

    705:    Input Parameters:
    706: +  pc  - the preconditioner context
    707: .  x   - the initial guess
    708: .  w   - one work vector
    709: .  rtol - relative decrease in residual norm convergence criteria
    710: .  abstol - absolute residual norm convergence criteria
    711: .  dtol - divergence residual norm increase criteria
    712: -  its - the number of iterations to apply.

    714:    Output Parameter:
    715: .  y - the solution

    717:    Notes:
    718:    Most preconditioners do not support this function. Use the command
    719:    PCApplyRichardsonExists() to determine if one does.

    721:    Except for the multigrid PC this routine ignores the convergence tolerances
    722:    and always runs for the number of iterations
    723:
    724:    Level: developer

    726: .keywords: PC, apply, Richardson

    728: .seealso: PCApplyRichardsonExists()
    729: @*/
    730: PetscErrorCode PETSCKSP_DLLEXPORT PCApplyRichardson(PC pc,Vec x,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its)
    731: {

    739:   if (!pc->ops->applyrichardson) SETERRQ(PETSC_ERR_SUP," ");

    741:   if (pc->setupcalled < 2) {
    742:     PCSetUp(pc);
    743:   }

    745:   (*pc->ops->applyrichardson)(pc,x,y,w,rtol,abstol,dtol,its);
    746:   return(0);
    747: }

    749: /*
    750:       a setupcall of 0 indicates never setup,
    751:                      1 needs to be resetup,
    752:                      2 does not need any changes.
    753: */
    756: /*@
    757:    PCSetUp - Prepares for the use of a preconditioner.

    759:    Collective on PC

    761:    Input Parameter:
    762: .  pc - the preconditioner context

    764:    Level: developer

    766: .keywords: PC, setup

    768: .seealso: PCCreate(), PCApply(), PCDestroy()
    769: @*/
    770: PetscErrorCode PETSCKSP_DLLEXPORT PCSetUp(PC pc)
    771: {
    773:   const char     *def;

    778:   if (pc->setupcalled > 1) {
    779:     PetscLogInfo((pc,"PCSetUp:Setting PC with identical preconditioner\n"));
    780:     return(0);
    781:   } else if (!pc->setupcalled) {
    782:     PetscLogInfo((pc,"PCSetUp:Setting up new PC\n"));
    783:   } else if (pc->flag == SAME_NONZERO_PATTERN) {
    784:     PetscLogInfo((pc,"PCSetUp:Setting up PC with same nonzero pattern\n"));
    785:   } else {
    786:     PetscLogInfo((pc,"PCSetUp:Setting up PC with different nonzero pattern\n"));
    787:   }

    789:   PetscLogEventBegin(PC_SetUp,pc,0,0,0);
    790:   if (!pc->mat) {SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Matrix must be set first");}

    792:   if (!pc->type_name) {
    793:     PCGetDefaultType_Private(pc,&def);
    794:     PCSetType(pc,def);
    795:   }

    797:   if (pc->ops->setup) {
    798:     (*pc->ops->setup)(pc);
    799:   }
    800:   pc->setupcalled = 2;
    801:   PetscLogEventEnd(PC_SetUp,pc,0,0,0);
    802:   return(0);
    803: }

    807: /*@
    808:    PCSetUpOnBlocks - Sets up the preconditioner for each block in
    809:    the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
    810:    methods.

    812:    Collective on PC

    814:    Input Parameters:
    815: .  pc - the preconditioner context

    817:    Level: developer

    819: .keywords: PC, setup, blocks

    821: .seealso: PCCreate(), PCApply(), PCDestroy(), PCSetUp()
    822: @*/
    823: PetscErrorCode PETSCKSP_DLLEXPORT PCSetUpOnBlocks(PC pc)
    824: {

    829:   if (!pc->ops->setuponblocks) return(0);
    830:   PetscLogEventBegin(PC_SetUpOnBlocks,pc,0,0,0);
    831:   (*pc->ops->setuponblocks)(pc);
    832:   PetscLogEventEnd(PC_SetUpOnBlocks,pc,0,0,0);
    833:   return(0);
    834: }

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

    845:    Collective on PC

    847:    Input Parameters:
    848: +  pc - the preconditioner context
    849: .  func - routine for modifying the submatrices
    850: -  ctx - optional user-defined context (may be null)

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

    855: .  row - an array of index sets that contain the global row numbers
    856:          that comprise each local submatrix
    857: .  col - an array of index sets that contain the global column numbers
    858:          that comprise each local submatrix
    859: .  submat - array of local submatrices
    860: -  ctx - optional user-defined context for private data for the
    861:          user-defined func routine (may be null)

    863:    Notes:
    864:    PCSetModifySubMatrices() MUST be called before KSPSetUp() and
    865:    KSPSolve().

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

    871:    Level: advanced

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

    875: .seealso: PCModifySubMatrices()
    876: @*/
    877: PetscErrorCode PETSCKSP_DLLEXPORT PCSetModifySubMatrices(PC pc,PetscErrorCode (*func)(PC,PetscInt,const IS[],const IS[],Mat[],void*),void *ctx)
    878: {
    881:   pc->modifysubmatrices  = func;
    882:   pc->modifysubmatricesP = ctx;
    883:   return(0);
    884: }

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

    892:    Collective on PC

    894:    Input Parameters:
    895: +  pc - the preconditioner context
    896: .  nsub - the number of local submatrices
    897: .  row - an array of index sets that contain the global row numbers
    898:          that comprise each local submatrix
    899: .  col - an array of index sets that contain the global column numbers
    900:          that comprise each local submatrix
    901: .  submat - array of local submatrices
    902: -  ctx - optional user-defined context for private data for the
    903:          user-defined routine (may be null)

    905:    Output Parameter:
    906: .  submat - array of local submatrices (the entries of which may
    907:             have been modified)

    909:    Notes:
    910:    The user should NOT generally call this routine, as it will
    911:    automatically be called within certain preconditioners (currently
    912:    block Jacobi, additive Schwarz) if set.

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

    919:    Level: developer

    921: .keywords: PC, modify, submatrices

    923: .seealso: PCSetModifySubMatrices()
    924: @*/
    925: PetscErrorCode PETSCKSP_DLLEXPORT PCModifySubMatrices(PC pc,PetscInt nsub,const IS row[],const IS col[],Mat submat[],void *ctx)
    926: {

    930:   if (!pc->modifysubmatrices) return(0);
    931:   PetscLogEventBegin(PC_ModifySubMatrices,pc,0,0,0);
    932:   (*pc->modifysubmatrices)(pc,nsub,row,col,submat,ctx);
    933:   PetscLogEventEnd(PC_ModifySubMatrices,pc,0,0,0);
    934:   return(0);
    935: }

    939: /*@
    940:    PCSetOperators - Sets the matrix associated with the linear system and
    941:    a (possibly) different one associated with the preconditioner.

    943:    Collective on PC and Mat

    945:    Input Parameters:
    946: +  pc - the preconditioner context
    947: .  Amat - the matrix associated with the linear system
    948: .  Pmat - the matrix to be used in constructing the preconditioner, usually
    949:           the same as Amat.
    950: -  flag - flag indicating information about the preconditioner matrix structure
    951:    during successive linear solves.  This flag is ignored the first time a
    952:    linear system is solved, and thus is irrelevant when solving just one linear
    953:    system.

    955:    Notes:
    956:    The flag can be used to eliminate unnecessary work in the preconditioner
    957:    during the repeated solution of linear systems of the same size.  The
    958:    available options are
    959: +    SAME_PRECONDITIONER -
    960:        Pmat is identical during successive linear solves.
    961:        This option is intended for folks who are using
    962:        different Amat and Pmat matrices and wish to reuse the
    963:        same preconditioner matrix.  For example, this option
    964:        saves work by not recomputing incomplete factorization
    965:        for ILU/ICC preconditioners.
    966: .     SAME_NONZERO_PATTERN -
    967:        Pmat has the same nonzero structure during
    968:        successive linear solves.
    969: -     DIFFERENT_NONZERO_PATTERN -
    970:        Pmat does not have the same nonzero structure.

    972:    Caution:
    973:    If you specify SAME_NONZERO_PATTERN, PETSc believes your assertion
    974:    and does not check the structure of the matrix.  If you erroneously
    975:    claim that the structure is the same when it actually is not, the new
    976:    preconditioner will not function correctly.  Thus, use this optimization
    977:    feature carefully!

    979:    If in doubt about whether your preconditioner matrix has changed
    980:    structure or not, use the flag DIFFERENT_NONZERO_PATTERN.

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

    988:    Level: intermediate

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

    992: .seealso: PCGetOperators(), MatZeroEntries()
    993:  @*/
    994: PetscErrorCode PETSCKSP_DLLEXPORT PCSetOperators(PC pc,Mat Amat,Mat Pmat,MatStructure flag)
    995: {
    997:   PetscTruth     isbjacobi,isrowbs;

    1004:   /*
    1005:       BlockSolve95 cannot use default BJacobi preconditioning
    1006:   */
    1007:   if (Amat) {
    1008:     PetscTypeCompare((PetscObject)Amat,MATMPIROWBS,&isrowbs);
    1009:     if (isrowbs) {
    1010:       PetscTypeCompare((PetscObject)pc,PCBJACOBI,&isbjacobi);
    1011:       if (isbjacobi) {
    1012:         PCSetType(pc,PCILU);
    1013:         PetscLogInfo((pc,"PCSetOperators:Switching default PC to PCILU since BS95 doesn't support PCBJACOBI\n"));
    1014:       }
    1015:     }
    1016:   }

    1018:   pc->mat  = Amat;
    1019:   pc->pmat = Pmat;
    1020:   if (pc->setupcalled == 2 && flag != SAME_PRECONDITIONER) {
    1021:     pc->setupcalled = 1;
    1022:   }
    1023:   pc->flag = flag;
    1024:   return(0);
    1025: }

    1029: /*@C
    1030:    PCGetOperators - Gets the matrix associated with the linear system and
    1031:    possibly a different one associated with the preconditioner.

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

    1035:    Input Parameter:
    1036: .  pc - the preconditioner context

    1038:    Output Parameters:
    1039: +  mat - the matrix associated with the linear system
    1040: .  pmat - matrix associated with the preconditioner, usually the same
    1041:           as mat.
    1042: -  flag - flag indicating information about the preconditioner
    1043:           matrix structure.  See PCSetOperators() for details.

    1045:    Level: intermediate

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

    1049: .seealso: PCSetOperators()
    1050: @*/
    1051: PetscErrorCode PETSCKSP_DLLEXPORT PCGetOperators(PC pc,Mat *mat,Mat *pmat,MatStructure *flag)
    1052: {
    1055:   if (mat)  *mat  = pc->mat;
    1056:   if (pmat) *pmat = pc->pmat;
    1057:   if (flag) *flag = pc->flag;
    1058:   return(0);
    1059: }

    1063: /*@C
    1064:    PCGetFactoredMatrix - Gets the factored matrix from the
    1065:    preconditioner context.  This routine is valid only for the LU,
    1066:    incomplete LU, Cholesky, and incomplete Cholesky methods.

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

    1070:    Input Parameters:
    1071: .  pc - the preconditioner context

    1073:    Output parameters:
    1074: .  mat - the factored matrix

    1076:    Level: advanced

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

    1080: .keywords: PC, get, factored, matrix
    1081: @*/
    1082: PetscErrorCode PETSCKSP_DLLEXPORT PCGetFactoredMatrix(PC pc,Mat *mat)
    1083: {

    1089:   if (pc->ops->getfactoredmatrix) {
    1090:     (*pc->ops->getfactoredmatrix)(pc,mat);
    1091:   }
    1092:   return(0);
    1093: }

    1097: /*@C
    1098:    PCSetOptionsPrefix - Sets the prefix used for searching for all
    1099:    PC options in the database.

    1101:    Collective on PC

    1103:    Input Parameters:
    1104: +  pc - the preconditioner context
    1105: -  prefix - the prefix string to prepend to all PC option requests

    1107:    Notes:
    1108:    A hyphen (-) must NOT be given at the beginning of the prefix name.
    1109:    The first character of all runtime options is AUTOMATICALLY the
    1110:    hyphen.

    1112:    Level: advanced

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

    1116: .seealso: PCAppendOptionsPrefix(), PCGetOptionsPrefix()
    1117: @*/
    1118: PetscErrorCode PETSCKSP_DLLEXPORT PCSetOptionsPrefix(PC pc,const char prefix[])
    1119: {

    1124:   PetscObjectSetOptionsPrefix((PetscObject)pc,prefix);
    1125:   return(0);
    1126: }

    1130: /*@C
    1131:    PCAppendOptionsPrefix - Appends to the prefix used for searching for all
    1132:    PC options in the database.

    1134:    Collective on PC

    1136:    Input Parameters:
    1137: +  pc - the preconditioner context
    1138: -  prefix - the prefix string to prepend to all PC option requests

    1140:    Notes:
    1141:    A hyphen (-) must NOT be given at the beginning of the prefix name.
    1142:    The first character of all runtime options is AUTOMATICALLY the
    1143:    hyphen.

    1145:    Level: advanced

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

    1149: .seealso: PCSetOptionsPrefix(), PCGetOptionsPrefix()
    1150: @*/
    1151: PetscErrorCode PETSCKSP_DLLEXPORT PCAppendOptionsPrefix(PC pc,const char prefix[])
    1152: {

    1157:   PetscObjectAppendOptionsPrefix((PetscObject)pc,prefix);
    1158:   return(0);
    1159: }

    1163: /*@C
    1164:    PCGetOptionsPrefix - Gets the prefix used for searching for all
    1165:    PC options in the database.

    1167:    Not Collective

    1169:    Input Parameters:
    1170: .  pc - the preconditioner context

    1172:    Output Parameters:
    1173: .  prefix - pointer to the prefix string used, is returned

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

    1178:    Level: advanced

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

    1182: .seealso: PCSetOptionsPrefix(), PCAppendOptionsPrefix()
    1183: @*/
    1184: PetscErrorCode PETSCKSP_DLLEXPORT PCGetOptionsPrefix(PC pc,const char *prefix[])
    1185: {

    1191:   PetscObjectGetOptionsPrefix((PetscObject)pc,prefix);
    1192:   return(0);
    1193: }

    1197: /*@
    1198:    PCPreSolve - Optional pre-solve phase, intended for any
    1199:    preconditioner-specific actions that must be performed before
    1200:    the iterative solve itself.

    1202:    Collective on PC

    1204:    Input Parameters:
    1205: +  pc - the preconditioner context
    1206: -  ksp - the Krylov subspace context

    1208:    Level: developer

    1210:    Sample of Usage:
    1211: .vb
    1212:     PCPreSolve(pc,ksp);
    1213:     KSPSolve(ksp,b,x);
    1214:     PCPostSolve(pc,ksp);
    1215: .ve

    1217:    Notes:
    1218:    The pre-solve phase is distinct from the PCSetUp() phase.

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

    1222: .keywords: PC, pre-solve

    1224: .seealso: PCPostSolve()
    1225: @*/
    1226: PetscErrorCode PETSCKSP_DLLEXPORT PCPreSolve(PC pc,KSP ksp)
    1227: {
    1229:   Vec            x,rhs;
    1230:   Mat            A,B;

    1235:   KSPGetSolution(ksp,&x);
    1236:   KSPGetRhs(ksp,&rhs);
    1237:   /*
    1238:       Scale the system and have the matrices use the scaled form
    1239:     only if the two matrices are actually the same (and hence
    1240:     have the same scaling
    1241:   */
    1242:   PCGetOperators(pc,&A,&B,PETSC_NULL);
    1243:   if (A == B) {
    1244:     MatScaleSystem(pc->mat,x,rhs);
    1245:     MatUseScaledForm(pc->mat,PETSC_TRUE);
    1246:   }

    1248:   if (pc->ops->presolve) {
    1249:     (*pc->ops->presolve)(pc,ksp,rhs,x);
    1250:   }
    1251:   return(0);
    1252: }

    1256: /*@
    1257:    PCPostSolve - Optional post-solve phase, intended for any
    1258:    preconditioner-specific actions that must be performed after
    1259:    the iterative solve itself.

    1261:    Collective on PC

    1263:    Input Parameters:
    1264: +  pc - the preconditioner context
    1265: -  ksp - the Krylov subspace context

    1267:    Sample of Usage:
    1268: .vb
    1269:     PCPreSolve(pc,ksp);
    1270:     KSPSolve(ksp,b,x);
    1271:     PCPostSolve(pc,ksp);
    1272: .ve

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

    1277:    Level: developer

    1279: .keywords: PC, post-solve

    1281: .seealso: PCPreSolve(), KSPSolve()
    1282: @*/
    1283: PetscErrorCode PETSCKSP_DLLEXPORT PCPostSolve(PC pc,KSP ksp)
    1284: {
    1286:   Vec            x,rhs;
    1287:   Mat            A,B;

    1292:   KSPGetSolution(ksp,&x);
    1293:   KSPGetRhs(ksp,&rhs);
    1294:   if (pc->ops->postsolve) {
    1295:      (*pc->ops->postsolve)(pc,ksp,rhs,x);
    1296:   }

    1298:   /*
    1299:       Scale the system and have the matrices use the scaled form
    1300:     only if the two matrices are actually the same (and hence
    1301:     have the same scaling
    1302:   */
    1303:   PCGetOperators(pc,&A,&B,PETSC_NULL);
    1304:   if (A == B) {
    1305:     MatUnScaleSystem(pc->mat,x,rhs);
    1306:     MatUseScaledForm(pc->mat,PETSC_FALSE);
    1307:   }
    1308:   return(0);
    1309: }

    1313: /*@C
    1314:    PCView - Prints the PC data structure.

    1316:    Collective on PC

    1318:    Input Parameters:
    1319: +  PC - the PC context
    1320: -  viewer - optional visualization context

    1322:    Note:
    1323:    The available visualization contexts include
    1324: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
    1325: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
    1326:          output where only the first processor opens
    1327:          the file.  All other processors send their
    1328:          data to the first processor to print.

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

    1333:    Level: developer

    1335: .keywords: PC, view

    1337: .seealso: KSPView(), PetscViewerASCIIOpen()
    1338: @*/
    1339: PetscErrorCode PETSCKSP_DLLEXPORT PCView(PC pc,PetscViewer viewer)
    1340: {
    1341:   PCType            cstr;
    1342:   PetscErrorCode    ierr;
    1343:   PetscTruth        mat_exists,iascii,isstring;
    1344:   PetscViewerFormat format;

    1348:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(pc->comm);

    1352:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
    1353:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
    1354:   if (iascii) {
    1355:     PetscViewerGetFormat(viewer,&format);
    1356:     if (pc->prefix) {
    1357:       PetscViewerASCIIPrintf(viewer,"PC Object:(%s)\n",pc->prefix);
    1358:     } else {
    1359:       PetscViewerASCIIPrintf(viewer,"PC Object:\n");
    1360:     }
    1361:     PCGetType(pc,&cstr);
    1362:     if (cstr) {
    1363:       PetscViewerASCIIPrintf(viewer,"  type: %s\n",cstr);
    1364:     } else {
    1365:       PetscViewerASCIIPrintf(viewer,"  type: not yet set\n");
    1366:     }
    1367:     if (pc->ops->view) {
    1368:       PetscViewerASCIIPushTab(viewer);
    1369:       (*pc->ops->view)(pc,viewer);
    1370:       PetscViewerASCIIPopTab(viewer);
    1371:     }
    1372:     PetscObjectExists((PetscObject)pc->mat,&mat_exists);
    1373:     if (mat_exists) {
    1374:       PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
    1375:       if (pc->pmat == pc->mat) {
    1376:         PetscViewerASCIIPrintf(viewer,"  linear system matrix = precond matrix:\n");
    1377:         PetscViewerASCIIPushTab(viewer);
    1378:         MatView(pc->mat,viewer);
    1379:         PetscViewerASCIIPopTab(viewer);
    1380:       } else {
    1381:         PetscObjectExists((PetscObject)pc->pmat,&mat_exists);
    1382:         if (mat_exists) {
    1383:           PetscViewerASCIIPrintf(viewer,"  linear system matrix followed by preconditioner matrix:\n");
    1384:         } else {
    1385:           PetscViewerASCIIPrintf(viewer,"  linear system matrix:\n");
    1386:         }
    1387:         PetscViewerASCIIPushTab(viewer);
    1388:         MatView(pc->mat,viewer);
    1389:         if (mat_exists) {MatView(pc->pmat,viewer);}
    1390:         PetscViewerASCIIPopTab(viewer);
    1391:       }
    1392:       PetscViewerPopFormat(viewer);
    1393:     }
    1394:   } else if (isstring) {
    1395:     PCGetType(pc,&cstr);
    1396:     PetscViewerStringSPrintf(viewer," %-7.7s",cstr);
    1397:     if (pc->ops->view) {(*pc->ops->view)(pc,viewer);}
    1398:   } else {
    1399:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by PC",((PetscObject)viewer)->type_name);
    1400:   }
    1401:   return(0);
    1402: }

    1406: /*@C
    1407:   PCRegister - See PCRegisterDynamic()

    1409:   Level: advanced
    1410: @*/
    1411: PetscErrorCode PETSCKSP_DLLEXPORT PCRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(PC))
    1412: {
    1414:   char           fullname[PETSC_MAX_PATH_LEN];

    1418:   PetscFListConcat(path,name,fullname);
    1419:   PetscFListAdd(&PCList,sname,fullname,(void (*)(void))function);
    1420:   return(0);
    1421: }

    1425: /*@
    1426:     PCComputeExplicitOperator - Computes the explicit preconditioned operator.

    1428:     Collective on PC

    1430:     Input Parameter:
    1431: .   pc - the preconditioner object

    1433:     Output Parameter:
    1434: .   mat - the explict preconditioned operator

    1436:     Notes:
    1437:     This computation is done by applying the operators to columns of the
    1438:     identity matrix.

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

    1444:     Level: advanced
    1445:
    1446: .keywords: PC, compute, explicit, operator

    1448: @*/
    1449: PetscErrorCode PETSCKSP_DLLEXPORT PCComputeExplicitOperator(PC pc,Mat *mat)
    1450: {
    1451:   Vec            in,out;
    1453:   PetscInt       i,M,m,*rows,start,end;
    1454:   PetscMPIInt    size;
    1455:   MPI_Comm       comm;
    1456:   PetscScalar    *array,zero = 0.0,one = 1.0;

    1462:   comm = pc->comm;
    1463:   MPI_Comm_size(comm,&size);

    1465:   if (!pc->pmat) SETERRQ(PETSC_ERR_ORDER,"You must call KSPSetOperators() or PCSetOperators() before this call");
    1466:   MatGetVecs(pc->pmat,&in,0);
    1467:   VecDuplicate(in,&out);
    1468:   VecGetOwnershipRange(in,&start,&end);
    1469:   VecGetSize(in,&M);
    1470:   VecGetLocalSize(in,&m);
    1471:   PetscMalloc((m+1)*sizeof(PetscInt),&rows);
    1472:   for (i=0; i<m; i++) {rows[i] = start + i;}

    1474:   MatCreate(comm,mat);
    1475:   MatSetSizes(*mat,m,m,M,M);
    1476:   if (size == 1) {
    1477:     MatSetType(*mat,MATSEQDENSE);
    1478:     MatSeqDenseSetPreallocation(*mat,PETSC_NULL);
    1479:   } else {
    1480:     MatSetType(*mat,MATMPIAIJ);
    1481:     MatMPIAIJSetPreallocation(*mat,0,PETSC_NULL,0,PETSC_NULL);
    1482:   }

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

    1486:     VecSet(in,zero);
    1487:     VecSetValues(in,1,&i,&one,INSERT_VALUES);
    1488:     VecAssemblyBegin(in);
    1489:     VecAssemblyEnd(in);

    1491:     /* should fix, allowing user to choose side */
    1492:     PCApply(pc,in,out);
    1493:
    1494:     VecGetArray(out,&array);
    1495:     MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
    1496:     VecRestoreArray(out,&array);

    1498:   }
    1499:   PetscFree(rows);
    1500:   VecDestroy(out);
    1501:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
    1502:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
    1503:   return(0);
    1504: }