Actual source code: precon.c

  1: #define PETSCKSP_DLL

  3: /*
  4:     The PC (preconditioner) interface routines, callable by users.
  5: */
 6:  #include private/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: /*@
 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  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:   if (pc->pmat) {MatDestroy(pc->pmat);}
 88:   if (pc->mat) {MatDestroy(pc->mat);}

 90:   PetscHeaderDestroy(pc);
 91:   return(0);
 92: }

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

100:    Collective on PC

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

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

108:    Level: developer

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

114: .keywords: PC

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

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

133:    Collective on PC

135:    Input Parameters:
136: +  pc - the preconditioner context
137: -  s - scaling vector

139:    Level: intermediate

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

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

147: .keywords: PC

149: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCDiagonalScale()
150: @*/
151: PetscErrorCode  PCDiagonalScaleSet(PC pc,Vec s)
152: {

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

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

178:    Collective on PC

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

185:    Level: intermediate

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

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

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

195: .keywords: PC

197: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleSet(), PCDiagonalScaleRight(), PCDiagonalScale()
198: @*/
199: PetscErrorCode  PCDiagonalScaleLeft(PC pc,Vec in,Vec out)
200: {

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

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

220:    Collective on PC

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

227:    Level: intermediate

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

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

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

237: .keywords: PC

239: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleSet(), PCDiagonalScale()
240: @*/
241: PetscErrorCode  PCDiagonalScaleRight(PC pc,Vec in,Vec out)
242: {

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

259: static PetscErrorCode PCPublish_Petsc(PetscObject obj)
260: {
262:   return(0);
263: }

267: /*@
268:    PCCreate - Creates a preconditioner context.

270:    Collective on MPI_Comm

272:    Input Parameter:
273: .  comm - MPI communicator 

275:    Output Parameter:
276: .  pc - location to put the preconditioner context

278:    Notes:
279:    The default preconditioner on one processor is PCILU with 0 fill on more 
280:    then one it is PCBJACOBI with ILU() on each processor.

282:    Level: developer

284: .keywords: PC, create, context

286: .seealso: PCSetUp(), PCApply(), PCDestroy()
287: @*/
288: PetscErrorCode  PCCreate(MPI_Comm comm,PC *newpc)
289: {
290:   PC             pc;

295:   *newpc = 0;
296: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
297:   PCInitializePackage(PETSC_NULL);
298: #endif

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

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

322:   pc->modifysubmatrices   = 0;
323:   pc->modifysubmatricesP  = 0;
324:   *newpc                  = pc;
325:   PetscPublishAll(pc);
326:   return(0);

328: }

330: /* -------------------------------------------------------------------------------*/

334: /*@
335:    PCApply - Applies the preconditioner to a vector.

337:    Collective on PC and Vec

339:    Input Parameters:
340: +  pc - the preconditioner context
341: -  x - input vector

343:    Output Parameter:
344: .  y - output vector

346:    Level: developer

348: .keywords: PC, apply

350: .seealso: PCApplyTranspose(), PCApplyBAorAB()
351: @*/
352: PetscErrorCode  PCApply(PC pc,Vec x,Vec y)
353: {

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

362:   if (pc->setupcalled < 2) {
363:     PCSetUp(pc);
364:   }
366:   (*pc->ops->apply)(pc,x,y);
368:   return(0);
369: }

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

376:    Collective on PC and Vec

378:    Input Parameters:
379: +  pc - the preconditioner context
380: -  x - input vector

382:    Output Parameter:
383: .  y - output vector

385:    Notes:
386:    Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.

388:    Level: developer

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

392: .seealso: PCApply(), PCApplySymmetricRight()
393: @*/
394: PetscErrorCode  PCApplySymmetricLeft(PC pc,Vec x,Vec y)
395: {

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

404:   if (pc->setupcalled < 2) {
405:     PCSetUp(pc);
406:   }

409:   (*pc->ops->applysymmetricleft)(pc,x,y);
411:   return(0);
412: }

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

419:    Collective on PC and Vec

421:    Input Parameters:
422: +  pc - the preconditioner context
423: -  x - input vector

425:    Output Parameter:
426: .  y - output vector

428:    Level: developer

430:    Notes:
431:    Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.

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

435: .seealso: PCApply(), PCApplySymmetricLeft()
436: @*/
437: PetscErrorCode  PCApplySymmetricRight(PC pc,Vec x,Vec y)
438: {

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

447:   if (pc->setupcalled < 2) {
448:     PCSetUp(pc);
449:   }

452:   (*pc->ops->applysymmetricright)(pc,x,y);
454:   return(0);
455: }

459: /*@
460:    PCApplyTranspose - Applies the transpose of preconditioner to a vector.

462:    Collective on PC and Vec

464:    Input Parameters:
465: +  pc - the preconditioner context
466: -  x - input vector

468:    Output Parameter:
469: .  y - output vector

471:    Level: developer

473: .keywords: PC, apply, transpose

475: .seealso: PCApply(), PCApplyBAorAB(), PCApplyBAorABTranspose(), PCHasApplyTranspose()
476: @*/
477: PetscErrorCode  PCApplyTranspose(PC pc,Vec x,Vec y)
478: {

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

488:   if (pc->setupcalled < 2) {
489:     PCSetUp(pc);
490:   }

493:   (*pc->ops->applytranspose)(pc,x,y);
495:   return(0);
496: }

500: /*@
501:    PCHasApplyTranspose - Test whether the preconditioner has a transpose apply operation

503:    Collective on PC and Vec

505:    Input Parameters:
506: .  pc - the preconditioner context

508:    Output Parameter:
509: .  flg - PETSC_TRUE if a transpose operation is defined

511:    Level: developer

513: .keywords: PC, apply, transpose

515: .seealso: PCApplyTranspose()
516: @*/
517: PetscErrorCode  PCHasApplyTranspose(PC pc,PetscTruth *flg)
518: {
522:   *flg = (PetscTruth) (pc->ops->applytranspose == 0);
523:   return(0);
524: }

528: /*@
529:    PCApplyBAorAB - Applies the preconditioner and operator to a vector. 

531:    Collective on PC and Vec

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

539:    Output Parameter:
540: .  y - output vector

542:    Level: developer

544: .keywords: PC, apply, operator

546: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorABTranspose()
547: @*/
548: PetscErrorCode  PCApplyBAorAB(PC pc,PCSide side,Vec x,Vec y,Vec work)
549: {

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

565:   if (pc->setupcalled < 2) {
566:     PCSetUp(pc);
567:   }

569:   if (pc->diagonalscale) {
570:     if (pc->ops->applyBA) {
571:       Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
572:       VecDuplicate(x,&work2);
573:       PCDiagonalScaleRight(pc,x,work2);
574:       (*pc->ops->applyBA)(pc,side,work2,y,work);
575:       PCDiagonalScaleLeft(pc,y,y);
576:       VecDestroy(work2);
577:     } else if (side == PC_RIGHT) {
578:       PCDiagonalScaleRight(pc,x,y);
579:       PCApply(pc,y,work);
580:       MatMult(pc->mat,work,y);
581:       PCDiagonalScaleLeft(pc,y,y);
582:     } else if (side == PC_LEFT) {
583:       PCDiagonalScaleRight(pc,x,y);
584:       MatMult(pc->mat,y,work);
585:       PCApply(pc,work,y);
586:       PCDiagonalScaleLeft(pc,y,y);
587:     } else if (side == PC_SYMMETRIC) {
588:       SETERRQ(PETSC_ERR_SUP,"Cannot provide diagonal scaling with symmetric application of preconditioner");
589:     }
590:   } else {
591:     if (pc->ops->applyBA) {
592:       (*pc->ops->applyBA)(pc,side,x,y,work);
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: }

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: PetscErrorCode  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: /* -------------------------------------------------------------------------------*/

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: PetscErrorCode  PCApplyRichardsonExists(PC pc,PetscTruth *exists)
690: {
694:   if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
695:   else                    *exists = PETSC_FALSE;
696:   return(0);
697: }

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: .  abstol - 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: PetscErrorCode  PCApplyRichardson(PC pc,Vec x,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt 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,abstol,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: */
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: PetscErrorCode  PCSetUp(PC pc)
774: {
776:   const char     *def;


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

793:   if (!pc->mat) {SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Matrix must be set first");}

795:   if (!pc->type_name) {
796:     PCGetDefaultType_Private(pc,&def);
797:     PCSetType(pc,def);
798:   }

800:   if (pc->ops->setup) {
801:     (*pc->ops->setup)(pc);
802:   }
803:   pc->setupcalled = 2;
805:   return(0);
806: }

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

815:    Collective on PC

817:    Input Parameters:
818: .  pc - the preconditioner context

820:    Level: developer

822: .keywords: PC, setup, blocks

824: .seealso: PCCreate(), PCApply(), PCDestroy(), PCSetUp()
825: @*/
826: PetscErrorCode  PCSetUpOnBlocks(PC pc)
827: {

832:   if (!pc->ops->setuponblocks) return(0);
834:   (*pc->ops->setuponblocks)(pc);
836:   return(0);
837: }

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

848:    Collective on PC

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

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

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

866:    Notes:
867:    PCSetModifySubMatrices() MUST be called before KSPSetUp() and
868:    KSPSolve().

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

874:    Level: advanced

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

878: .seealso: PCModifySubMatrices()
879: @*/
880: PetscErrorCode  PCSetModifySubMatrices(PC pc,PetscErrorCode (*func)(PC,PetscInt,const IS[],const IS[],Mat[],void*),void *ctx)
881: {
884:   pc->modifysubmatrices  = func;
885:   pc->modifysubmatricesP = ctx;
886:   return(0);
887: }

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

895:    Collective on PC

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

908:    Output Parameter:
909: .  submat - array of local submatrices (the entries of which may
910:             have been modified)

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

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

922:    Level: developer

924: .keywords: PC, modify, submatrices

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

933:   if (!pc->modifysubmatrices) return(0);
935:   (*pc->modifysubmatrices)(pc,nsub,row,col,submat,ctx);
937:   return(0);
938: }

942: /*@
943:    PCSetOperators - Sets the matrix associated with the linear system and 
944:    a (possibly) different one associated with the preconditioner.

946:    Collective on PC and Mat

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

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

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

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

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

991:    Level: intermediate

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

995: .seealso: PCGetOperators(), MatZeroEntries()
996:  @*/
997: PetscErrorCode  PCSetOperators(PC pc,Mat Amat,Mat Pmat,MatStructure flag)
998: {
1000:   PetscTruth     isbjacobi,isrowbs;


1009:   /*
1010:       BlockSolve95 cannot use default BJacobi preconditioning
1011:   */
1012:   if (Amat) {
1013:     PetscTypeCompare((PetscObject)Amat,MATMPIROWBS,&isrowbs);
1014:     if (isrowbs) {
1015:       PetscTypeCompare((PetscObject)pc,PCBJACOBI,&isbjacobi);
1016:       if (isbjacobi) {
1017:         PCSetType(pc,PCILU);
1018:         PetscInfo(pc,"Switching default PC to PCILU since BS95 doesn't support PCBJACOBI\n");
1019:       }
1020:     }
1021:   }

1023:   /* reference first in case the matrices are the same */
1024:   if (Amat) {PetscObjectReference((PetscObject)Amat);}
1025:   if (pc->mat) {MatDestroy(pc->mat);}
1026:   if (Pmat) {PetscObjectReference((PetscObject)Pmat);}
1027:   if (pc->pmat) {MatDestroy(pc->pmat);}
1028:   pc->mat  = Amat;
1029:   pc->pmat = Pmat;

1031:   if (pc->setupcalled == 2 && flag != SAME_PRECONDITIONER) {
1032:     pc->setupcalled = 1;
1033:   }
1034:   pc->flag = flag;
1035:   return(0);
1036: }

1040: /*@C
1041:    PCGetOperators - Gets the matrix associated with the linear system and
1042:    possibly a different one associated with the preconditioner.

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

1046:    Input Parameter:
1047: .  pc - the preconditioner context

1049:    Output Parameters:
1050: +  mat - the matrix associated with the linear system
1051: .  pmat - matrix associated with the preconditioner, usually the same
1052:           as mat. 
1053: -  flag - flag indicating information about the preconditioner
1054:           matrix structure.  See PCSetOperators() for details.

1056:    Level: intermediate

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

1066: $         KSP/PCGetOperators(ksp/pc,&mat,PETSC_NULL,PETSC_NULL); is equivalent to
1067: $           set size, type, etc of mat

1069: $         MatCreate(comm,&mat);
1070: $         KSP/PCSetOperators(ksp/pc,mat,mat,SAME_NONZERO_PATTERN);
1071: $         PetscObjectDereference((PetscObject)mat);
1072: $           set size, type, etc of mat

1074:      and

1076: $         KSP/PCGetOperators(ksp/pc,&mat,&pmat,PETSC_NULL); is equivalent to
1077: $           set size, type, etc of mat and pmat

1079: $         MatCreate(comm,&mat);
1080: $         MatCreate(comm,&pmat);
1081: $         KSP/PCSetOperators(ksp/pc,mat,pmat,SAME_NONZERO_PATTERN);
1082: $         PetscObjectDereference((PetscObject)mat);
1083: $         PetscObjectDereference((PetscObject)pmat);
1084: $           set size, type, etc of mat and pmat

1086:     The rational for this support is so that when creating a TS, SNES, or KSP the hierarchy
1087:     of underlying objects (i.e. SNES, KSP, PC, Mat) and their livespans can be completely 
1088:     managed by the top most level object (i.e. the TS, SNES, or KSP). Another way to look
1089:     at this is when you create a SNES you do not NEED to create a KSP and attach it to 
1090:     the SNES object (the SNES object manages it for you). Similarly when you create a KSP
1091:     you do not need to attach a PC to it (the KSP object manages the PC object for you).
1092:     Thus, why should YOU have to create the Mat and attach it to the SNES/KSP/PC, when
1093:     it can be created for you?
1094:      

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

1098: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperatorsSet()
1099: @*/
1100: PetscErrorCode  PCGetOperators(PC pc,Mat *mat,Mat *pmat,MatStructure *flag)
1101: {

1106:   if (mat) {
1107:     if (!pc->mat) {
1108:       MatCreate(pc->comm,&pc->mat);
1109:       if (!pc->pmat && !pmat) { /* user did NOT request pmat, so make same as mat */
1110:         pc->pmat = pc->mat;
1111:         PetscObjectReference((PetscObject)pc->pmat);
1112:       }
1113:     }
1114:     *mat  = pc->mat;
1115:   }
1116:   if (pmat) {
1117:     if (!pc->pmat) {
1118:       MatCreate(pc->comm,&pc->mat);
1119:       if (!pc->mat && !mat) { /* user did NOT request mat, so make same as pmat */
1120:         pc->mat = pc->pmat;
1121:         PetscObjectReference((PetscObject)pc->mat);
1122:       }
1123:     }
1124:     *pmat = pc->pmat;
1125:   }
1126:   if (flag) *flag = pc->flag;
1127:   return(0);
1128: }

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

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

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

1141:    Output Parameters:
1142: +  mat - the matrix associated with the linear system was set
1143: -  pmat - matrix associated with the preconditioner was set, usually the same

1145:    Level: intermediate

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

1149: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperators()
1150: @*/
1151: PetscErrorCode  PCGetOperatorsSet(PC pc,PetscTruth *mat,PetscTruth *pmat)
1152: {
1155:   if (mat)  *mat  = (pc->mat)  ? PETSC_TRUE : PETSC_FALSE;
1156:   if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1157:   return(0);
1158: }

1162: /*@
1163:    PCGetFactoredMatrix - Gets the factored matrix from the
1164:    preconditioner context.  This routine is valid only for the LU, 
1165:    incomplete LU, Cholesky, and incomplete Cholesky methods.

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

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

1172:    Output parameters:
1173: .  mat - the factored matrix

1175:    Level: advanced

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

1179: .keywords: PC, get, factored, matrix
1180: @*/
1181: PetscErrorCode  PCGetFactoredMatrix(PC pc,Mat *mat)
1182: {

1188:   if (pc->ops->getfactoredmatrix) {
1189:     (*pc->ops->getfactoredmatrix)(pc,mat);
1190:   }
1191:   return(0);
1192: }

1196: /*@C
1197:    PCSetOptionsPrefix - Sets the prefix used for searching for all 
1198:    PC options in the database.

1200:    Collective on PC

1202:    Input Parameters:
1203: +  pc - the preconditioner context
1204: -  prefix - the prefix string to prepend to all PC option requests

1206:    Notes:
1207:    A hyphen (-) must NOT be given at the beginning of the prefix name.
1208:    The first character of all runtime options is AUTOMATICALLY the
1209:    hyphen.

1211:    Level: advanced

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

1215: .seealso: PCAppendOptionsPrefix(), PCGetOptionsPrefix()
1216: @*/
1217: PetscErrorCode  PCSetOptionsPrefix(PC pc,const char prefix[])
1218: {

1223:   PetscObjectSetOptionsPrefix((PetscObject)pc,prefix);
1224:   return(0);
1225: }

1229: /*@C
1230:    PCAppendOptionsPrefix - Appends to the prefix used for searching for all 
1231:    PC options in the database.

1233:    Collective on PC

1235:    Input Parameters:
1236: +  pc - the preconditioner context
1237: -  prefix - the prefix string to prepend to all PC option requests

1239:    Notes:
1240:    A hyphen (-) must NOT be given at the beginning of the prefix name.
1241:    The first character of all runtime options is AUTOMATICALLY the
1242:    hyphen.

1244:    Level: advanced

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

1248: .seealso: PCSetOptionsPrefix(), PCGetOptionsPrefix()
1249: @*/
1250: PetscErrorCode  PCAppendOptionsPrefix(PC pc,const char prefix[])
1251: {

1256:   PetscObjectAppendOptionsPrefix((PetscObject)pc,prefix);
1257:   return(0);
1258: }

1262: /*@C
1263:    PCGetOptionsPrefix - Gets the prefix used for searching for all 
1264:    PC options in the database.

1266:    Not Collective

1268:    Input Parameters:
1269: .  pc - the preconditioner context

1271:    Output Parameters:
1272: .  prefix - pointer to the prefix string used, is returned

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

1277:    Level: advanced

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

1281: .seealso: PCSetOptionsPrefix(), PCAppendOptionsPrefix()
1282: @*/
1283: PetscErrorCode  PCGetOptionsPrefix(PC pc,const char *prefix[])
1284: {

1290:   PetscObjectGetOptionsPrefix((PetscObject)pc,prefix);
1291:   return(0);
1292: }

1296: /*@
1297:    PCPreSolve - Optional pre-solve phase, intended for any
1298:    preconditioner-specific actions that must be performed before 
1299:    the iterative solve itself.

1301:    Collective on PC

1303:    Input Parameters:
1304: +  pc - the preconditioner context
1305: -  ksp - the Krylov subspace context

1307:    Level: developer

1309:    Sample of Usage:
1310: .vb
1311:     PCPreSolve(pc,ksp);
1312:     KSPSolve(ksp,b,x);
1313:     PCPostSolve(pc,ksp);
1314: .ve

1316:    Notes:
1317:    The pre-solve phase is distinct from the PCSetUp() phase.

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

1321: .keywords: PC, pre-solve

1323: .seealso: PCPostSolve()
1324: @*/
1325: PetscErrorCode  PCPreSolve(PC pc,KSP ksp)
1326: {
1328:   Vec            x,rhs;
1329:   Mat            A,B;

1334:   KSPGetSolution(ksp,&x);
1335:   KSPGetRhs(ksp,&rhs);
1336:   /*
1337:       Scale the system and have the matrices use the scaled form
1338:     only if the two matrices are actually the same (and hence
1339:     have the same scaling
1340:   */
1341:   PCGetOperators(pc,&A,&B,PETSC_NULL);
1342:   if (A == B) {
1343:     MatScaleSystem(pc->mat,rhs,x);
1344:     MatUseScaledForm(pc->mat,PETSC_TRUE);
1345:   }

1347:   if (pc->ops->presolve) {
1348:     (*pc->ops->presolve)(pc,ksp,rhs,x);
1349:   }
1350:   return(0);
1351: }

1355: /*@
1356:    PCPostSolve - Optional post-solve phase, intended for any
1357:    preconditioner-specific actions that must be performed after
1358:    the iterative solve itself.

1360:    Collective on PC

1362:    Input Parameters:
1363: +  pc - the preconditioner context
1364: -  ksp - the Krylov subspace context

1366:    Sample of Usage:
1367: .vb
1368:     PCPreSolve(pc,ksp);
1369:     KSPSolve(ksp,b,x);
1370:     PCPostSolve(pc,ksp);
1371: .ve

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

1376:    Level: developer

1378: .keywords: PC, post-solve

1380: .seealso: PCPreSolve(), KSPSolve()
1381: @*/
1382: PetscErrorCode  PCPostSolve(PC pc,KSP ksp)
1383: {
1385:   Vec            x,rhs;
1386:   Mat            A,B;

1391:   KSPGetSolution(ksp,&x);
1392:   KSPGetRhs(ksp,&rhs);
1393:   if (pc->ops->postsolve) {
1394:      (*pc->ops->postsolve)(pc,ksp,rhs,x);
1395:   }

1397:   /*
1398:       Scale the system and have the matrices use the scaled form
1399:     only if the two matrices are actually the same (and hence
1400:     have the same scaling
1401:   */
1402:   PCGetOperators(pc,&A,&B,PETSC_NULL);
1403:   if (A == B) {
1404:     MatUnScaleSystem(pc->mat,rhs,x);
1405:     MatUseScaledForm(pc->mat,PETSC_FALSE);
1406:   }
1407:   return(0);
1408: }

1412: /*@C
1413:    PCView - Prints the PC data structure.

1415:    Collective on PC

1417:    Input Parameters:
1418: +  PC - the PC context
1419: -  viewer - optional visualization context

1421:    Note:
1422:    The available visualization contexts include
1423: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
1424: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1425:          output where only the first processor opens
1426:          the file.  All other processors send their 
1427:          data to the first processor to print. 

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

1432:    Level: developer

1434: .keywords: PC, view

1436: .seealso: KSPView(), PetscViewerASCIIOpen()
1437: @*/
1438: PetscErrorCode  PCView(PC pc,PetscViewer viewer)
1439: {
1440:   PCType            cstr;
1441:   PetscErrorCode    ierr;
1442:   PetscTruth        mat_exists,iascii,isstring;
1443:   PetscViewerFormat format;

1447:   if (!viewer) {
1448:     PetscViewerASCIIGetStdout(pc->comm,&viewer);
1449:   }

1453:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
1454:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
1455:   if (iascii) {
1456:     PetscViewerGetFormat(viewer,&format);
1457:     if (pc->prefix) {
1458:       PetscViewerASCIIPrintf(viewer,"PC Object:(%s)\n",pc->prefix);
1459:     } else {
1460:       PetscViewerASCIIPrintf(viewer,"PC Object:\n");
1461:     }
1462:     PCGetType(pc,&cstr);
1463:     if (cstr) {
1464:       PetscViewerASCIIPrintf(viewer,"  type: %s\n",cstr);
1465:     } else {
1466:       PetscViewerASCIIPrintf(viewer,"  type: not yet set\n");
1467:     }
1468:     if (pc->ops->view) {
1469:       PetscViewerASCIIPushTab(viewer);
1470:       (*pc->ops->view)(pc,viewer);
1471:       PetscViewerASCIIPopTab(viewer);
1472:     }
1473:     PetscObjectExists((PetscObject)pc->mat,&mat_exists);
1474:     if (mat_exists) {
1475:       PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
1476:       if (pc->pmat == pc->mat) {
1477:         PetscViewerASCIIPrintf(viewer,"  linear system matrix = precond matrix:\n");
1478:         PetscViewerASCIIPushTab(viewer);
1479:         MatView(pc->mat,viewer);
1480:         PetscViewerASCIIPopTab(viewer);
1481:       } else {
1482:         PetscObjectExists((PetscObject)pc->pmat,&mat_exists);
1483:         if (mat_exists) {
1484:           PetscViewerASCIIPrintf(viewer,"  linear system matrix followed by preconditioner matrix:\n");
1485:         } else {
1486:           PetscViewerASCIIPrintf(viewer,"  linear system matrix:\n");
1487:         }
1488:         PetscViewerASCIIPushTab(viewer);
1489:         MatView(pc->mat,viewer);
1490:         if (mat_exists) {MatView(pc->pmat,viewer);}
1491:         PetscViewerASCIIPopTab(viewer);
1492:       }
1493:       PetscViewerPopFormat(viewer);
1494:     }
1495:   } else if (isstring) {
1496:     PCGetType(pc,&cstr);
1497:     PetscViewerStringSPrintf(viewer," %-7.7s",cstr);
1498:     if (pc->ops->view) {(*pc->ops->view)(pc,viewer);}
1499:   } else {
1500:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by PC",((PetscObject)viewer)->type_name);
1501:   }
1502:   return(0);
1503: }

1507: /*@C
1508:   PCRegister - See PCRegisterDynamic()

1510:   Level: advanced
1511: @*/
1512: PetscErrorCode  PCRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(PC))
1513: {
1515:   char           fullname[PETSC_MAX_PATH_LEN];


1519:   PetscFListConcat(path,name,fullname);
1520:   PetscFListAdd(&PCList,sname,fullname,(void (*)(void))function);
1521:   return(0);
1522: }

1526: /*@
1527:     PCComputeExplicitOperator - Computes the explicit preconditioned operator.  

1529:     Collective on PC

1531:     Input Parameter:
1532: .   pc - the preconditioner object

1534:     Output Parameter:
1535: .   mat - the explict preconditioned operator

1537:     Notes:
1538:     This computation is done by applying the operators to columns of the 
1539:     identity matrix.

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

1545:     Level: advanced
1546:    
1547: .keywords: PC, compute, explicit, operator

1549: @*/
1550: PetscErrorCode  PCComputeExplicitOperator(PC pc,Mat *mat)
1551: {
1552:   Vec            in,out;
1554:   PetscInt       i,M,m,*rows,start,end;
1555:   PetscMPIInt    size;
1556:   MPI_Comm       comm;
1557:   PetscScalar    *array,one = 1.0;
1558: 

1563:   comm = pc->comm;
1564:   MPI_Comm_size(comm,&size);

1566:   if (!pc->pmat) SETERRQ(PETSC_ERR_ORDER,"You must call KSPSetOperators() or PCSetOperators() before this call");
1567:   MatGetVecs(pc->pmat,&in,0);
1568:   VecDuplicate(in,&out);
1569:   VecGetOwnershipRange(in,&start,&end);
1570:   VecGetSize(in,&M);
1571:   VecGetLocalSize(in,&m);
1572:   PetscMalloc((m+1)*sizeof(PetscInt),&rows);
1573:   for (i=0; i<m; i++) {rows[i] = start + i;}

1575:   MatCreate(comm,mat);
1576:   MatSetSizes(*mat,m,m,M,M);
1577:   if (size == 1) {
1578:     MatSetType(*mat,MATSEQDENSE);
1579:     MatSeqDenseSetPreallocation(*mat,PETSC_NULL);
1580:   } else {
1581:     MatSetType(*mat,MATMPIAIJ);
1582:     MatMPIAIJSetPreallocation(*mat,0,PETSC_NULL,0,PETSC_NULL);
1583:   }

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

1587:     VecSet(in,0.0);
1588:     VecSetValues(in,1,&i,&one,INSERT_VALUES);
1589:     VecAssemblyBegin(in);
1590:     VecAssemblyEnd(in);

1592:     /* should fix, allowing user to choose side */
1593:     PCApply(pc,in,out);
1594: 
1595:     VecGetArray(out,&array);
1596:     MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
1597:     VecRestoreArray(out,&array);

1599:   }
1600:   PetscFree(rows);
1601:   VecDestroy(out);
1602:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
1603:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
1604:   return(0);
1605: }