Actual source code: precon.c

  1: /*$Id: precon.c,v 1.216 2001/08/21 21:03:13 bsmith Exp $*/
  2: /*
  3:     The PC (preconditioner) interface routines, callable by users.
  4: */
 5:  #include src/sles/pc/pcimpl.h

  7: /* Logging support */
  8: int PC_COOKIE;
  9: int PC_SetUp, PC_SetUpOnBlocks, PC_Apply, PC_ApplyCoarse, PC_ApplyMultiple, PC_ApplySymmetricLeft;
 10: int PC_ApplySymmetricRight, PC_ModifySubMatrices;

 12: EXTERN int SLESInitializePackage(char *);

 14: #undef __FUNCT__  
 16: /*@C
 17:    PCNullSpaceAttach - attaches a null space to a preconditioner object.
 18:         This null space will be removed from the resulting vector whenever
 19:         the preconditioner is applied.

 21:    Collective on PC

 23:    Input Parameters:
 24: +  pc - the preconditioner context
 25: -  nullsp - the null space object

 27:    Level: developer

 29:    Notes:
 30:       Overwrites any previous null space that may have been attached

 32:   Users manual sections:
 33: .   Section 4.15 Solving Singular Systems

 35: .keywords: PC, destroy, null space

 37: .seealso: PCCreate(), PCSetUp(), MatNullSpaceCreate(), MatNullSpace

 39: @*/
 40: int PCNullSpaceAttach(PC pc,MatNullSpace nullsp)
 41: {


 48:   if (pc->nullsp) {
 49:     MatNullSpaceDestroy(pc->nullsp);
 50:   }
 51:   pc->nullsp = nullsp;
 52:   PetscObjectReference((PetscObject)nullsp);
 53:   return(0);
 54: }

 56: #undef __FUNCT__  
 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: int 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->nullsp)             {MatNullSpaceDestroy(pc->nullsp);}
 85:   if (pc->diagonalscaleright) {VecDestroy(pc->diagonalscaleright);}
 86:   if (pc->diagonalscaleleft)  {VecDestroy(pc->diagonalscaleleft);}

 88:   PetscLogObjectDestroy(pc);
 89:   PetscHeaderDestroy(pc);
 90:   return(0);
 91: }

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

 99:    Collective on PC

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

104:    Output Parameter:
105: .  flag - PETSC_TRUE if it applies the scaling

107:    Level: developer

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

113: .keywords: PC

115: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCDiagonalScaleSet()
116: @*/
117: int PCDiagonalScale(PC pc,PetscTruth *flag)
118: {
121:   *flag = pc->diagonalscale;
122:   return(0);
123: }

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

131:    Collective on PC

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

137:    Level: intermediate

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

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

145: .keywords: PC

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

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

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

176:    Collective on PC

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

183:    Level: intermediate

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

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

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

193: .keywords: PC

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

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

213: #undef __FUNCT__  
215: /*@C
216:    PCDiagonalScaleRight - Scales a vector by the right scaling as needed by certain time-stepping codes.

218:    Collective on PC

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

225:    Level: intermediate

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

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

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

235: .keywords: PC

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

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

255: #undef __FUNCT__  
257: static int PCPublish_Petsc(PetscObject obj)
258: {
259: #if defined(PETSC_HAVE_AMS)
260:   PC          v = (PC) obj;
261:   int         ierr;
262: #endif


266: #if defined(PETSC_HAVE_AMS)
267:   /* if it is already published then return */
268:   if (v->amem >=0) return(0);

270:   PetscObjectPublishBaseBegin(obj);
271:   PetscObjectPublishBaseEnd(obj);
272: #endif

274:   return(0);
275: }

277: #undef __FUNCT__  
279: /*@C
280:    PCCreate - Creates a preconditioner context.

282:    Collective on MPI_Comm

284:    Input Parameter:
285: .  comm - MPI communicator 

287:    Output Parameter:
288: .  pc - location to put the preconditioner context

290:    Notes:
291:    The default preconditioner on one processor is PCILU with 0 fill on more 
292:    then one it is PCBJACOBI with ILU() on each processor.

294:    Level: developer

296: .keywords: PC, create, context

298: .seealso: PCSetUp(), PCApply(), PCDestroy()
299: @*/
300: int PCCreate(MPI_Comm comm,PC *newpc)
301: {
302:   PC  pc;

307:   *newpc = 0;
308: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
309:   SLESInitializePackage(PETSC_NULL);
310: #endif

312:   PetscHeaderCreate(pc,_p_PC,struct _PCOps,PC_COOKIE,-1,"PC",comm,PCDestroy,PCView);
313:   PetscLogObjectCreate(pc);
314:   pc->bops->publish      = PCPublish_Petsc;
315:   pc->vec                = 0;
316:   pc->mat                = 0;
317:   pc->setupcalled        = 0;
318:   pc->nullsp             = 0;
319:   pc->data               = 0;
320:   pc->diagonalscale      = PETSC_FALSE;
321:   pc->diagonalscaleleft  = 0;
322:   pc->diagonalscaleright = 0;

324:   pc->ops->destroy             = 0;
325:   pc->ops->apply               = 0;
326:   pc->ops->applytranspose      = 0;
327:   pc->ops->applyBA             = 0;
328:   pc->ops->applyBAtranspose    = 0;
329:   pc->ops->applyrichardson     = 0;
330:   pc->ops->view                = 0;
331:   pc->ops->getfactoredmatrix   = 0;
332:   pc->ops->applysymmetricright = 0;
333:   pc->ops->applysymmetricleft  = 0;
334:   pc->ops->setuponblocks       = 0;

336:   pc->modifysubmatrices   = 0;
337:   pc->modifysubmatricesP  = 0;
338:   *newpc                  = pc;
339:   PetscPublishAll(pc);
340:   return(0);

342: }

344: /* -------------------------------------------------------------------------------*/

346: #undef __FUNCT__  
348: /*@
349:    PCApply - Applies the preconditioner to a vector.

351:    Collective on PC and Vec

353:    Input Parameters:
354: +  pc - the preconditioner context
355: -  x - input vector

357:    Output Parameter:
358: .  y - output vector

360:    Level: developer

362: .keywords: PC, apply

364: .seealso: PCApplyTranspose(), PCApplyBAorAB()
365: @*/
366: int PCApply(PC pc,Vec x,Vec y)
367: {
368:   int        ierr;

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

376:   if (pc->setupcalled < 2) {
377:     PCSetUp(pc);
378:   }

380:   PetscLogEventBegin(PC_Apply,pc,x,y,0);
381:   (*pc->ops->apply)(pc,x,y);
382:   PetscLogEventEnd(PC_Apply,pc,x,y,0);

384:   /* Remove null space from preconditioned vector y */
385:   if (pc->nullsp) {
386:     MatNullSpaceRemove(pc->nullsp,y,PETSC_NULL);
387:   }
388:   return(0);
389: }

391: #undef __FUNCT__  
393: /*@
394:    PCApplySymmetricLeft - Applies the left part of a symmetric preconditioner to a vector.

396:    Collective on PC and Vec

398:    Input Parameters:
399: +  pc - the preconditioner context
400: -  x - input vector

402:    Output Parameter:
403: .  y - output vector

405:    Notes:
406:    Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.

408:    Level: developer

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

412: .seealso: PCApply(), PCApplySymmetricRight()
413: @*/
414: int PCApplySymmetricLeft(PC pc,Vec x,Vec y)
415: {

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

424:   if (pc->setupcalled < 2) {
425:     PCSetUp(pc);
426:   }

428:   PetscLogEventBegin(PC_ApplySymmetricLeft,pc,x,y,0);
429:   (*pc->ops->applysymmetricleft)(pc,x,y);
430:   PetscLogEventEnd(PC_ApplySymmetricLeft,pc,x,y,0);
431:   return(0);
432: }

434: #undef __FUNCT__  
436: /*@
437:    PCApplySymmetricRight - Applies the right part of a symmetric preconditioner to a vector.

439:    Collective on PC and Vec

441:    Input Parameters:
442: +  pc - the preconditioner context
443: -  x - input vector

445:    Output Parameter:
446: .  y - output vector

448:    Level: developer

450:    Notes:
451:    Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.

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

455: .seealso: PCApply(), PCApplySymmetricLeft()
456: @*/
457: int PCApplySymmetricRight(PC pc,Vec x,Vec y)
458: {

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

467:   if (pc->setupcalled < 2) {
468:     PCSetUp(pc);
469:   }

471:   PetscLogEventBegin(PC_ApplySymmetricRight,pc,x,y,0);
472:   (*pc->ops->applysymmetricright)(pc,x,y);
473:   PetscLogEventEnd(PC_ApplySymmetricRight,pc,x,y,0);
474:   return(0);
475: }

477: #undef __FUNCT__  
479: /*@
480:    PCApplyTranspose - Applies the transpose of preconditioner to a vector.

482:    Collective on PC and Vec

484:    Input Parameters:
485: +  pc - the preconditioner context
486: -  x - input vector

488:    Output Parameter:
489: .  y - output vector

491:    Level: developer

493: .keywords: PC, apply, transpose

495: .seealso: PCApplyTranspose(), PCApplyBAorAB(), PCApplyBAorABTranspose()
496: @*/
497: int PCApplyTranspose(PC pc,Vec x,Vec y)
498: {

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

508:   if (pc->setupcalled < 2) {
509:     PCSetUp(pc);
510:   }

512:   PetscLogEventBegin(PC_Apply,pc,x,y,0);
513:   (*pc->ops->applytranspose)(pc,x,y);
514:   PetscLogEventEnd(PC_Apply,pc,x,y,0);
515:   return(0);
516: }

518: #undef __FUNCT__  
520: /*@
521:    PCApplyBAorAB - Applies the preconditioner and operator to a vector. 

523:    Collective on PC and Vec

525:    Input Parameters:
526: +  pc - the preconditioner context
527: .  side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
528: .  x - input vector
529: -  work - work vector

531:    Output Parameter:
532: .  y - output vector

534:    Level: developer

536: .keywords: PC, apply, operator

538: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorABTranspose()
539: @*/
540: int PCApplyBAorAB(PC pc,PCSide side,Vec x,Vec y,Vec work)
541: {
542:   int        ierr;

549:   if (x == y) SETERRQ(PETSC_ERR_ARG_IDN,"x and y must be different vectors");
550:   if (side != PC_LEFT && side != PC_SYMMETRIC && side != PC_RIGHT) {
551:     SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Side must be right, left, or symmetric");
552:   }
553:   if (pc->diagonalscale && side == PC_SYMMETRIC) {
554:     SETERRQ(1,"Cannot include diagonal scaling with symmetric preconditioner application");
555:   }

557:   if (pc->setupcalled < 2) {
558:     PCSetUp(pc);
559:   }

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

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

617:    Collective on PC and Vec

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

625:    Output Parameter:
626: .  y - output vector

628:    Level: developer

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

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

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

652:   if (pc->setupcalled < 2) {
653:     PCSetUp(pc);
654:   }

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

667: /* -------------------------------------------------------------------------------*/

669: #undef __FUNCT__  
671: /*@
672:    PCApplyRichardsonExists - Determines whether a particular preconditioner has a 
673:    built-in fast application of Richardson's method.

675:    Not Collective

677:    Input Parameter:
678: .  pc - the preconditioner

680:    Output Parameter:
681: .  exists - PETSC_TRUE or PETSC_FALSE

683:    Level: developer

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

687: .seealso: PCApplyRichardson()
688: @*/
689: int PCApplyRichardsonExists(PC pc,PetscTruth *exists)
690: {
694:   if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
695:   else                    *exists = PETSC_FALSE;
696:   return(0);
697: }

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

706:    Collective on PC

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

717:    Output Parameter:
718: .  y - the solution

720:    Notes: 
721:    Most preconditioners do not support this function. Use the command
722:    PCApplyRichardsonExists() to determine if one does.

724:    Except for the multigrid PC this routine ignores the convergence tolerances
725:    and always runs for the number of iterations
726:  
727:    Level: developer

729: .keywords: PC, apply, Richardson

731: .seealso: PCApplyRichardsonExists()
732: @*/
733: int PCApplyRichardson(PC pc,Vec x,Vec y,Vec w,PetscReal rtol,PetscReal atol, PetscReal dtol,int its)
734: {

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

744:   if (pc->setupcalled < 2) {
745:     PCSetUp(pc);
746:   }

748:   (*pc->ops->applyrichardson)(pc,x,y,w,rtol,atol,dtol,its);
749:   return(0);
750: }

752: /* 
753:       a setupcall of 0 indicates never setup, 
754:                      1 needs to be resetup,
755:                      2 does not need any changes.
756: */
757: #undef __FUNCT__  
759: /*@
760:    PCSetUp - Prepares for the use of a preconditioner.

762:    Collective on PC

764:    Input Parameter:
765: .  pc - the preconditioner context

767:    Level: developer

769: .keywords: PC, setup

771: .seealso: PCCreate(), PCApply(), PCDestroy()
772: @*/
773: int PCSetUp(PC pc)
774: {
776:   PetscTruth flg;


781:   if (pc->setupcalled > 1) {
782:     PetscLogInfo(pc,"PCSetUp:Setting PC with identical preconditionern");
783:     return(0);
784:   } else if (pc->setupcalled == 0) {
785:     PetscLogInfo(pc,"PCSetUp:Setting up new PCn");
786:   } else if (pc->flag == SAME_NONZERO_PATTERN) {
787:     PetscLogInfo(pc,"PCSetUp:Setting up PC with same nonzero patternn");
788:   } else {
789:     PetscLogInfo(pc,"PCSetUp:Setting up PC with different nonzero patternn");
790:   }

792:   PetscLogEventBegin(PC_SetUp,pc,0,0,0);
793:   if (!pc->vec) {SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Vector must be set first");}
794:   if (!pc->mat) {SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Matrix must be set first");}

796:   if (!pc->type_name) {
797:     int size;

799:     MPI_Comm_size(pc->comm,&size);
800:     if (size == 1) {
801:       PetscTypeCompare((PetscObject)pc->pmat,MATSEQSBAIJ,&flg);
802:       if (flg) { /* for sbaij mat */
803:         PCSetType(pc,PCICC);
804:       } else {
805:         PCSetType(pc,PCILU);
806:       }
807:     } else {
808:       PCSetType(pc,PCBJACOBI);
809:     }
810:   }
811:   if (pc->ops->setup) {
812:     (*pc->ops->setup)(pc);
813:   }
814:   pc->setupcalled = 2;
815:   if (pc->nullsp) {
816:     PetscTruth test;
817:     PetscOptionsHasName(pc->prefix,"-pc_test_null_space",&test);
818:     if (test) {
819:       MatNullSpaceTest(pc->nullsp,pc->mat);
820:     }
821:   }
822:   PetscLogEventEnd(PC_SetUp,pc,0,0,0);
823:   return(0);
824: }

826: #undef __FUNCT__  
828: /*@
829:    PCSetUpOnBlocks - Sets up the preconditioner for each block in
830:    the block Jacobi, block Gauss-Seidel, and overlapping Schwarz 
831:    methods.

833:    Collective on PC

835:    Input Parameters:
836: .  pc - the preconditioner context

838:    Level: developer

840: .keywords: PC, setup, blocks

842: .seealso: PCCreate(), PCApply(), PCDestroy(), PCSetUp()
843: @*/
844: int PCSetUpOnBlocks(PC pc)
845: {

850:   if (!pc->ops->setuponblocks) return(0);
851:   PetscLogEventBegin(PC_SetUpOnBlocks,pc,0,0,0);
852:   (*pc->ops->setuponblocks)(pc);
853:   PetscLogEventEnd(PC_SetUpOnBlocks,pc,0,0,0);
854:   return(0);
855: }

857: #undef __FUNCT__  
859: /*@
860:    PCSetModifySubMatrices - Sets a user-defined routine for modifying the
861:    submatrices that arise within certain subdomain-based preconditioners.
862:    The basic submatrices are extracted from the preconditioner matrix as
863:    usual; the user can then alter these (for example, to set different boundary
864:    conditions for each submatrix) before they are used for the local solves.

866:    Collective on PC

868:    Input Parameters:
869: +  pc - the preconditioner context
870: .  func - routine for modifying the submatrices
871: -  ctx - optional user-defined context (may be null)

873:    Calling sequence of func:
874: $     func (PC pc,int nsub,IS *row,IS *col,Mat *submat,void *ctx);

876: .  row - an array of index sets that contain the global row numbers
877:          that comprise each local submatrix
878: .  col - an array of index sets that contain the global column numbers
879:          that comprise each local submatrix
880: .  submat - array of local submatrices
881: -  ctx - optional user-defined context for private data for the 
882:          user-defined func routine (may be null)

884:    Notes:
885:    PCSetModifySubMatrices() MUST be called before SLESSetUp() and
886:    SLESSolve().

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

892:    Level: advanced

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

896: .seealso: PCModifySubMatrices()
897: @*/
898: int PCSetModifySubMatrices(PC pc,int(*func)(PC,int,IS*,IS*,Mat*,void*),void *ctx)
899: {
902:   pc->modifysubmatrices  = func;
903:   pc->modifysubmatricesP = ctx;
904:   return(0);
905: }

907: #undef __FUNCT__  
909: /*@
910:    PCModifySubMatrices - Calls an optional user-defined routine within 
911:    certain preconditioners if one has been set with PCSetModifySubMarices().

913:    Collective on PC

915:    Input Parameters:
916: +  pc - the preconditioner context
917: .  nsub - the number of local submatrices
918: .  row - an array of index sets that contain the global row numbers
919:          that comprise each local submatrix
920: .  col - an array of index sets that contain the global column numbers
921:          that comprise each local submatrix
922: .  submat - array of local submatrices
923: -  ctx - optional user-defined context for private data for the 
924:          user-defined routine (may be null)

926:    Output Parameter:
927: .  submat - array of local submatrices (the entries of which may
928:             have been modified)

930:    Notes:
931:    The user should NOT generally call this routine, as it will
932:    automatically be called within certain preconditioners (currently
933:    block Jacobi, additive Schwarz) if set.

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

940:    Level: developer

942: .keywords: PC, modify, submatrices

944: .seealso: PCSetModifySubMatrices()
945: @*/
946: int PCModifySubMatrices(PC pc,int nsub,IS *row,IS *col,Mat *submat,void *ctx)
947: {

951:   if (!pc->modifysubmatrices) return(0);
952:   PetscLogEventBegin(PC_ModifySubMatrices,pc,0,0,0);
953:   (*pc->modifysubmatrices)(pc,nsub,row,col,submat,ctx);
954:   PetscLogEventEnd(PC_ModifySubMatrices,pc,0,0,0);
955:   return(0);
956: }

958: #undef __FUNCT__  
960: /*@
961:    PCSetOperators - Sets the matrix associated with the linear system and 
962:    a (possibly) different one associated with the preconditioner.

964:    Collective on PC and Mat

966:    Input Parameters:
967: +  pc - the preconditioner context
968: .  Amat - the matrix associated with the linear system
969: .  Pmat - the matrix to be used in constructing the preconditioner, usually
970:           the same as Amat. 
971: -  flag - flag indicating information about the preconditioner matrix structure
972:    during successive linear solves.  This flag is ignored the first time a
973:    linear system is solved, and thus is irrelevant when solving just one linear
974:    system.

976:    Notes: 
977:    The flag can be used to eliminate unnecessary work in the preconditioner 
978:    during the repeated solution of linear systems of the same size.  The 
979:    available options are
980: +    SAME_PRECONDITIONER -
981:        Pmat is identical during successive linear solves.
982:        This option is intended for folks who are using
983:        different Amat and Pmat matrices and wish to reuse the
984:        same preconditioner matrix.  For example, this option
985:        saves work by not recomputing incomplete factorization
986:        for ILU/ICC preconditioners.
987: .     SAME_NONZERO_PATTERN -
988:        Pmat has the same nonzero structure during
989:        successive linear solves. 
990: -     DIFFERENT_NONZERO_PATTERN -
991:        Pmat does not have the same nonzero structure.

993:    Caution:
994:    If you specify SAME_NONZERO_PATTERN, PETSc believes your assertion
995:    and does not check the structure of the matrix.  If you erroneously
996:    claim that the structure is the same when it actually is not, the new
997:    preconditioner will not function correctly.  Thus, use this optimization
998:    feature carefully!

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

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

1009:    Level: intermediate

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

1013: .seealso: PCGetOperators(), MatZeroEntries()
1014:  @*/
1015: int PCSetOperators(PC pc,Mat Amat,Mat Pmat,MatStructure flag)
1016: {
1017:   int        ierr;
1018:   PetscTruth isbjacobi,isrowbs;


1025:   /*
1026:       BlockSolve95 cannot use default BJacobi preconditioning
1027:   */
1028:   PetscTypeCompare((PetscObject)Amat,MATMPIROWBS,&isrowbs);
1029:   if (isrowbs) {
1030:     PetscTypeCompare((PetscObject)pc,PCBJACOBI,&isbjacobi);
1031:     if (isbjacobi) {
1032:       PCSetType(pc,PCILU);
1033:       PetscLogInfo(pc,"PCSetOperators:Switching default PC to PCILU since BS95 doesn't support PCBJACOBIn");
1034:     }
1035:   }

1037:   pc->mat  = Amat;
1038:   pc->pmat = Pmat;
1039:   if (pc->setupcalled == 2 && flag != SAME_PRECONDITIONER) {
1040:     pc->setupcalled = 1;
1041:   }
1042:   pc->flag = flag;
1043:   return(0);
1044: }

1046: #undef __FUNCT__  
1048: /*@C
1049:    PCGetOperators - Gets the matrix associated with the linear system and
1050:    possibly a different one associated with the preconditioner.

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

1054:    Input Parameter:
1055: .  pc - the preconditioner context

1057:    Output Parameters:
1058: +  mat - the matrix associated with the linear system
1059: .  pmat - matrix associated with the preconditioner, usually the same
1060:           as mat. 
1061: -  flag - flag indicating information about the preconditioner
1062:           matrix structure.  See PCSetOperators() for details.

1064:    Level: intermediate

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

1068: .seealso: PCSetOperators()
1069: @*/
1070: int PCGetOperators(PC pc,Mat *mat,Mat *pmat,MatStructure *flag)
1071: {
1074:   if (mat)  *mat  = pc->mat;
1075:   if (pmat) *pmat = pc->pmat;
1076:   if (flag) *flag = pc->flag;
1077:   return(0);
1078: }

1080: #undef __FUNCT__  
1082: /*@
1083:    PCSetVector - Sets a vector associated with the preconditioner.

1085:    Collective on PC and Vec

1087:    Input Parameters:
1088: +  pc - the preconditioner context
1089: -  vec - the vector

1091:    Notes:
1092:    The vector must be set so that the preconditioner knows what type
1093:    of vector to allocate if necessary.

1095:    Level: intermediate

1097: .keywords: PC, set, vector

1099: .seealso: PCGetVector()

1101: @*/
1102: int PCSetVector(PC pc,Vec vec)
1103: {
1108:   pc->vec = vec;
1109:   return(0);
1110: }

1112: #undef __FUNCT__  
1114: /*@
1115:    PCGetVector - Gets a vector associated with the preconditioner; if the 
1116:    vector was not get set it will return a 0 pointer.

1118:    Not collective, but vector is shared by all processors that share the PC

1120:    Input Parameter:
1121: .  pc - the preconditioner context

1123:    Output Parameter:
1124: .  vec - the vector

1126:    Level: intermediate

1128: .keywords: PC, get, vector

1130: .seealso: PCSetVector()

1132: @*/
1133: int PCGetVector(PC pc,Vec *vec)
1134: {
1137:   *vec = pc->vec;
1138:   return(0);
1139: }

1141: #undef __FUNCT__  
1143: /*@C 
1144:    PCGetFactoredMatrix - Gets the factored matrix from the
1145:    preconditioner context.  This routine is valid only for the LU, 
1146:    incomplete LU, Cholesky, and incomplete Cholesky methods.

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

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

1153:    Output parameters:
1154: .  mat - the factored matrix

1156:    Level: advanced

1158: .keywords: PC, get, factored, matrix
1159: @*/
1160: int PCGetFactoredMatrix(PC pc,Mat *mat)
1161: {

1166:   if (pc->ops->getfactoredmatrix) {
1167:     (*pc->ops->getfactoredmatrix)(pc,mat);
1168:   }
1169:   return(0);
1170: }

1172: #undef __FUNCT__  
1174: /*@C
1175:    PCSetOptionsPrefix - Sets the prefix used for searching for all 
1176:    PC options in the database.

1178:    Collective on PC

1180:    Input Parameters:
1181: +  pc - the preconditioner context
1182: -  prefix - the prefix string to prepend to all PC option requests

1184:    Notes:
1185:    A hyphen (-) must NOT be given at the beginning of the prefix name.
1186:    The first character of all runtime options is AUTOMATICALLY the
1187:    hyphen.

1189:    Level: advanced

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

1193: .seealso: PCAppendOptionsPrefix(), PCGetOptionsPrefix()
1194: @*/
1195: int PCSetOptionsPrefix(PC pc,char *prefix)
1196: {

1201:   PetscObjectSetOptionsPrefix((PetscObject)pc,prefix);
1202:   return(0);
1203: }

1205: #undef __FUNCT__  
1207: /*@C
1208:    PCAppendOptionsPrefix - Appends to the prefix used for searching for all 
1209:    PC options in the database.

1211:    Collective on PC

1213:    Input Parameters:
1214: +  pc - the preconditioner context
1215: -  prefix - the prefix string to prepend to all PC option requests

1217:    Notes:
1218:    A hyphen (-) must NOT be given at the beginning of the prefix name.
1219:    The first character of all runtime options is AUTOMATICALLY the
1220:    hyphen.

1222:    Level: advanced

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

1226: .seealso: PCSetOptionsPrefix(), PCGetOptionsPrefix()
1227: @*/
1228: int PCAppendOptionsPrefix(PC pc,char *prefix)
1229: {

1234:   PetscObjectAppendOptionsPrefix((PetscObject)pc,prefix);
1235:   return(0);
1236: }

1238: #undef __FUNCT__  
1240: /*@C
1241:    PCGetOptionsPrefix - Gets the prefix used for searching for all 
1242:    PC options in the database.

1244:    Not Collective

1246:    Input Parameters:
1247: .  pc - the preconditioner context

1249:    Output Parameters:
1250: .  prefix - pointer to the prefix string used, is returned

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

1255:    Level: advanced

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

1259: .seealso: PCSetOptionsPrefix(), PCAppendOptionsPrefix()
1260: @*/
1261: int PCGetOptionsPrefix(PC pc,char **prefix)
1262: {

1267:   PetscObjectGetOptionsPrefix((PetscObject)pc,prefix);
1268:   return(0);
1269: }

1271: #undef __FUNCT__  
1273: /*@
1274:    PCPreSolve - Optional pre-solve phase, intended for any
1275:    preconditioner-specific actions that must be performed before 
1276:    the iterative solve itself.

1278:    Collective on PC

1280:    Input Parameters:
1281: +  pc - the preconditioner context
1282: -  ksp - the Krylov subspace context

1284:    Level: developer

1286:    Sample of Usage:
1287: .vb
1288:     PCPreSolve(pc,ksp);
1289:     KSPSolve(ksp,its);
1290:     PCPostSolve(pc,ksp);
1291: .ve

1293:    Notes:
1294:    The pre-solve phase is distinct from the PCSetUp() phase.

1296:    SLESSolve() calls this directly, so is rarely called by the user.

1298: .keywords: PC, pre-solve

1300: .seealso: PCPostSolve()
1301: @*/
1302: int PCPreSolve(PC pc,KSP ksp)
1303: {
1305:   Vec x,rhs;
1306:   Mat A,B;


1311:   KSPGetSolution(ksp,&x);
1312:   KSPGetRhs(ksp,&rhs);
1313:   /*
1314:       Scale the system and have the matrices use the scaled form
1315:     only if the two matrices are actually the same (and hence
1316:     have the same scaling
1317:   */
1318:   PCGetOperators(pc,&A,&B,PETSC_NULL);
1319:   if (A == B) {
1320:     MatScaleSystem(pc->mat,x,rhs);
1321:     MatUseScaledForm(pc->mat,PETSC_TRUE);
1322:   }

1324:   if (pc->ops->presolve) {
1325:     (*pc->ops->presolve)(pc,ksp,x,rhs);
1326:   }
1327:   return(0);
1328: }

1330: #undef __FUNCT__  
1332: /*@
1333:    PCPostSolve - Optional post-solve phase, intended for any
1334:    preconditioner-specific actions that must be performed after
1335:    the iterative solve itself.

1337:    Collective on PC

1339:    Input Parameters:
1340: +  pc - the preconditioner context
1341: -  ksp - the Krylov subspace context

1343:    Sample of Usage:
1344: .vb
1345:     PCPreSolve(pc,ksp);
1346:     KSPSolve(ksp,its);
1347:     PCPostSolve(pc,ksp);
1348: .ve

1350:    Note:
1351:    SLESSolve() calls this routine directly, so it is rarely called by the user.

1353:    Level: developer

1355: .keywords: PC, post-solve

1357: .seealso: PCPreSolve(), SLESSolve()
1358: @*/
1359: int PCPostSolve(PC pc,KSP ksp)
1360: {
1362:   Vec x,rhs;
1363:   Mat A,B;

1367:   KSPGetSolution(ksp,&x);
1368:   KSPGetRhs(ksp,&rhs);
1369:   if (pc->ops->postsolve) {
1370:      (*pc->ops->postsolve)(pc,ksp,x,rhs);
1371:   }

1373:   /*
1374:       Scale the system and have the matrices use the scaled form
1375:     only if the two matrices are actually the same (and hence
1376:     have the same scaling
1377:   */
1378:   PCGetOperators(pc,&A,&B,PETSC_NULL);
1379:   if (A == B) {
1380:     MatUnScaleSystem(pc->mat,x,rhs);
1381:     MatUseScaledForm(pc->mat,PETSC_FALSE);
1382:   }
1383:   return(0);
1384: }

1386: #undef __FUNCT__  
1388: /*@C
1389:    PCView - Prints the PC data structure.

1391:    Collective on PC

1393:    Input Parameters:
1394: +  PC - the PC context
1395: -  viewer - optional visualization context

1397:    Note:
1398:    The available visualization contexts include
1399: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
1400: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1401:          output where only the first processor opens
1402:          the file.  All other processors send their 
1403:          data to the first processor to print. 

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

1408:    Level: developer

1410: .keywords: PC, view

1412: .seealso: KSPView(), PetscViewerASCIIOpen()
1413: @*/
1414: int PCView(PC pc,PetscViewer viewer)
1415: {
1416:   PCType            cstr;
1417:   int               ierr;
1418:   PetscTruth        mat_exists,isascii,isstring;
1419:   PetscViewerFormat format;

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

1427:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
1428:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
1429:   if (isascii) {
1430:     PetscViewerGetFormat(viewer,&format);
1431:     if (pc->prefix) {
1432:       PetscViewerASCIIPrintf(viewer,"PC Object:(%s)n",pc->prefix);
1433:     } else {
1434:       PetscViewerASCIIPrintf(viewer,"PC Object:n");
1435:     }
1436:     PCGetType(pc,&cstr);
1437:     if (cstr) {
1438:       PetscViewerASCIIPrintf(viewer,"  type: %sn",cstr);
1439:     } else {
1440:       PetscViewerASCIIPrintf(viewer,"  type: not yet setn");
1441:     }
1442:     if (pc->ops->view) {
1443:       PetscViewerASCIIPushTab(viewer);
1444:       (*pc->ops->view)(pc,viewer);
1445:       PetscViewerASCIIPopTab(viewer);
1446:     }
1447:     PetscObjectExists((PetscObject)pc->mat,&mat_exists);
1448:     if (mat_exists) {
1449:       PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
1450:       if (pc->pmat == pc->mat) {
1451:         PetscViewerASCIIPrintf(viewer,"  linear system matrix = precond matrix:n");
1452:         PetscViewerASCIIPushTab(viewer);
1453:         MatView(pc->mat,viewer);
1454:         PetscViewerASCIIPopTab(viewer);
1455:       } else {
1456:         PetscObjectExists((PetscObject)pc->pmat,&mat_exists);
1457:         if (mat_exists) {
1458:           PetscViewerASCIIPrintf(viewer,"  linear system matrix followed by preconditioner matrix:n");
1459:         } else {
1460:           PetscViewerASCIIPrintf(viewer,"  linear system matrix:n");
1461:         }
1462:         PetscViewerASCIIPushTab(viewer);
1463:         MatView(pc->mat,viewer);
1464:         if (mat_exists) {MatView(pc->pmat,viewer);}
1465:         PetscViewerASCIIPopTab(viewer);
1466:       }
1467:       PetscViewerPopFormat(viewer);
1468:     }
1469:   } else if (isstring) {
1470:     PCGetType(pc,&cstr);
1471:     PetscViewerStringSPrintf(viewer," %-7.7s",cstr);
1472:     if (pc->ops->view) {(*pc->ops->view)(pc,viewer);}
1473:   } else {
1474:     SETERRQ1(1,"Viewer type %s not supported by PC",((PetscObject)viewer)->type_name);
1475:   }
1476:   return(0);
1477: }

1479: /*MC
1480:    PCRegisterDynamic - Adds a method to the preconditioner package.

1482:    Synopsis:
1483:    int PCRegisterDynamic(char *name_solver,char *path,char *name_create,int (*routine_create)(PC))

1485:    Not collective

1487:    Input Parameters:
1488: +  name_solver - name of a new user-defined solver
1489: .  path - path (either absolute or relative) the library containing this solver
1490: .  name_create - name of routine to create method context
1491: -  routine_create - routine to create method context

1493:    Notes:
1494:    PCRegisterDynamic() may be called multiple times to add several user-defined preconditioners.

1496:    If dynamic libraries are used, then the fourth input argument (routine_create)
1497:    is ignored.

1499:    Sample usage:
1500: .vb
1501:    PCRegisterDynamic("my_solver","/home/username/my_lib/lib/libO/solaris/mylib",
1502:               "MySolverCreate",MySolverCreate);
1503: .ve

1505:    Then, your solver can be chosen with the procedural interface via
1506: $     PCSetType(pc,"my_solver")
1507:    or at runtime via the option
1508: $     -pc_type my_solver

1510:    Level: advanced

1512:    Notes: ${PETSC_ARCH}, ${PETSC_DIR}, ${PETSC_LIB_DIR}, ${BOPT}, or ${any environmental variable}
1513:            occuring in pathname will be replaced with appropriate values.
1514:          If your function is not being put into a shared library then use PCRegister() instead

1516: .keywords: PC, register

1518: .seealso: PCRegisterAll(), PCRegisterDestroy()
1519: M*/

1521: #undef __FUNCT__  
1523: /*@C
1524:   PCRegister - See PCRegisterDynamic()

1526:   Level: advanced
1527: @*/
1528: int PCRegister(char *sname,char *path,char *name,int (*function)(PC))
1529: {
1530:   int  ierr;
1531:   char fullname[256];


1535:   PetscFListConcat(path,name,fullname);
1536:   PetscFListAdd(&PCList,sname,fullname,(void (*)(void))function);
1537:   return(0);
1538: }

1540: #undef __FUNCT__  
1542: /*@
1543:     PCComputeExplicitOperator - Computes the explicit preconditioned operator.  

1545:     Collective on PC

1547:     Input Parameter:
1548: .   pc - the preconditioner object

1550:     Output Parameter:
1551: .   mat - the explict preconditioned operator

1553:     Notes:
1554:     This computation is done by applying the operators to columns of the 
1555:     identity matrix.

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

1561:     Level: advanced
1562:    
1563: .keywords: PC, compute, explicit, operator

1565: @*/
1566: int PCComputeExplicitOperator(PC pc,Mat *mat)
1567: {
1568:   Vec      in,out;
1569:   int      ierr,i,M,m,size,*rows,start,end;
1570:   MPI_Comm comm;
1571:   PetscScalar   *array,zero = 0.0,one = 1.0;


1577:   comm = pc->comm;
1578:   MPI_Comm_size(comm,&size);

1580:   PCGetVector(pc,&in);
1581:   VecDuplicate(in,&out);
1582:   VecGetOwnershipRange(in,&start,&end);
1583:   VecGetSize(in,&M);
1584:   VecGetLocalSize(in,&m);
1585:   PetscMalloc((m+1)*sizeof(int),&rows);
1586:   for (i=0; i<m; i++) {rows[i] = start + i;}

1588:   if (size == 1) {
1589:     MatCreateSeqDense(comm,M,M,PETSC_NULL,mat);
1590:   } else {
1591:     MatCreateMPIAIJ(comm,m,m,M,M,0,0,0,0,mat);
1592:   }

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

1596:     VecSet(&zero,in);
1597:     VecSetValues(in,1,&i,&one,INSERT_VALUES);
1598:     VecAssemblyBegin(in);
1599:     VecAssemblyEnd(in);

1601:     PCApply(pc,in,out);
1602: 
1603:     VecGetArray(out,&array);
1604:     MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
1605:     VecRestoreArray(out,&array);

1607:   }
1608:   PetscFree(rows);
1609:   VecDestroy(out);
1610:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
1611:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
1612:   return(0);
1613: }