Actual source code: mgfunc.c

petsc-master 2015-07-30
Report Typos and Errors
  2: #include <petsc/private/pcmgimpl.h>       /*I "petscksp.h" I*/

  4: /* ---------------------------------------------------------------------------*/
  7: /*@C
  8:    PCMGResidualDefault - Default routine to calculate the residual.

 10:    Collective on Mat and Vec

 12:    Input Parameters:
 13: +  mat - the matrix
 14: .  b   - the right-hand-side
 15: -  x   - the approximate solution

 17:    Output Parameter:
 18: .  r - location to store the residual

 20:    Level: developer

 22: .keywords: MG, default, multigrid, residual

 24: .seealso: PCMGSetResidual()
 25: @*/
 26: PetscErrorCode  PCMGResidualDefault(Mat mat,Vec b,Vec x,Vec r)
 27: {

 31:   MatResidual(mat,b,x,r);
 32:   return(0);
 33: }

 37: /*@
 38:    PCMGGetCoarseSolve - Gets the solver context to be used on the coarse grid.

 40:    Not Collective

 42:    Input Parameter:
 43: .  pc - the multigrid context

 45:    Output Parameter:
 46: .  ksp - the coarse grid solver context

 48:    Level: advanced

 50: .keywords: MG, multigrid, get, coarse grid
 51: @*/
 52: PetscErrorCode  PCMGGetCoarseSolve(PC pc,KSP *ksp)
 53: {
 54:   PC_MG        *mg        = (PC_MG*)pc->data;
 55:   PC_MG_Levels **mglevels = mg->levels;

 59:   *ksp =  mglevels[0]->smoothd;
 60:   return(0);
 61: }

 65: /*@C
 66:    PCMGSetResidual - Sets the function to be used to calculate the residual
 67:    on the lth level.

 69:    Logically Collective on PC and Mat

 71:    Input Parameters:
 72: +  pc       - the multigrid context
 73: .  l        - the level (0 is coarsest) to supply
 74: .  residual - function used to form residual, if none is provided the previously provide one is used, if no
 75:               previous one were provided then a default is used
 76: -  mat      - matrix associated with residual

 78:    Level: advanced

 80: .keywords:  MG, set, multigrid, residual, level

 82: .seealso: PCMGResidualDefault()
 83: @*/
 84: PetscErrorCode  PCMGSetResidual(PC pc,PetscInt l,PetscErrorCode (*residual)(Mat,Vec,Vec,Vec),Mat mat)
 85: {
 86:   PC_MG          *mg        = (PC_MG*)pc->data;
 87:   PC_MG_Levels   **mglevels = mg->levels;

 92:   if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
 93:   if (residual) mglevels[l]->residual = residual;
 94:   if (!mglevels[l]->residual) mglevels[l]->residual = PCMGResidualDefault;
 95:   if (mat) {PetscObjectReference((PetscObject)mat);}
 96:   MatDestroy(&mglevels[l]->A);
 97:   mglevels[l]->A = mat;
 98:   return(0);
 99: }

103: /*@
104:    PCMGSetInterpolation - Sets the function to be used to calculate the
105:    interpolation from l-1 to the lth level

107:    Logically Collective on PC and Mat

109:    Input Parameters:
110: +  pc  - the multigrid context
111: .  mat - the interpolation operator
112: -  l   - the level (0 is coarsest) to supply [do not supply 0]

114:    Level: advanced

116:    Notes:
117:           Usually this is the same matrix used also to set the restriction
118:     for the same level.

120:           One can pass in the interpolation matrix or its transpose; PETSc figures
121:     out from the matrix size which one it is.

123: .keywords:  multigrid, set, interpolate, level

125: .seealso: PCMGSetRestriction()
126: @*/
127: PetscErrorCode  PCMGSetInterpolation(PC pc,PetscInt l,Mat mat)
128: {
129:   PC_MG          *mg        = (PC_MG*)pc->data;
130:   PC_MG_Levels   **mglevels = mg->levels;

135:   if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
136:   if (!l) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Do not set interpolation routine for coarsest level");
137:   PetscObjectReference((PetscObject)mat);
138:   MatDestroy(&mglevels[l]->interpolate);

140:   mglevels[l]->interpolate = mat;
141:   return(0);
142: }

146: /*@
147:    PCMGGetInterpolation - Gets the function to be used to calculate the
148:    interpolation from l-1 to the lth level

150:    Logically Collective on PC

152:    Input Parameters:
153: +  pc - the multigrid context
154: -  l - the level (0 is coarsest) to supply [Do not supply 0]

156:    Output Parameter:
157: .  mat - the interpolation matrix

159:    Level: advanced

161: .keywords: MG, get, multigrid, interpolation, level

163: .seealso: PCMGGetRestriction(), PCMGSetInterpolation(), PCMGGetRScale()
164: @*/
165: PetscErrorCode  PCMGGetInterpolation(PC pc,PetscInt l,Mat *mat)
166: {
167:   PC_MG          *mg        = (PC_MG*)pc->data;
168:   PC_MG_Levels   **mglevels = mg->levels;

174:   if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
175:   if (l <= 0 || mg->nlevels <= l) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Level %D must be in range {1,...,%D}",l,mg->nlevels-1);
176:   if (!mglevels[l]->interpolate) {
177:     if (!mglevels[l]->restrct) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call PCMGSetInterpolation() or PCMGSetInterpolation()");
178:     PCMGSetInterpolation(pc,l,mglevels[l]->restrct);
179:   }
180:   *mat = mglevels[l]->interpolate;
181:   return(0);
182: }

186: /*@
187:    PCMGSetRestriction - Sets the function to be used to restrict vector
188:    from level l to l-1.

190:    Logically Collective on PC and Mat

192:    Input Parameters:
193: +  pc - the multigrid context
194: .  l - the level (0 is coarsest) to supply [Do not supply 0]
195: -  mat - the restriction matrix

197:    Level: advanced

199:    Notes:
200:           Usually this is the same matrix used also to set the interpolation
201:     for the same level.

203:           One can pass in the interpolation matrix or its transpose; PETSc figures
204:     out from the matrix size which one it is.

206:          If you do not set this, the transpose of the Mat set with PCMGSetInterpolation()
207:     is used.

209: .keywords: MG, set, multigrid, restriction, level

211: .seealso: PCMGSetInterpolation()
212: @*/
213: PetscErrorCode  PCMGSetRestriction(PC pc,PetscInt l,Mat mat)
214: {
216:   PC_MG          *mg        = (PC_MG*)pc->data;
217:   PC_MG_Levels   **mglevels = mg->levels;

222:   if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
223:   if (!l) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Do not set restriction routine for coarsest level");
224:   PetscObjectReference((PetscObject)mat);
225:   MatDestroy(&mglevels[l]->restrct);

227:   mglevels[l]->restrct = mat;
228:   return(0);
229: }

233: /*@
234:    PCMGGetRestriction - Gets the function to be used to restrict vector
235:    from level l to l-1.

237:    Logically Collective on PC and Mat

239:    Input Parameters:
240: +  pc - the multigrid context
241: -  l - the level (0 is coarsest) to supply [Do not supply 0]

243:    Output Parameter:
244: .  mat - the restriction matrix

246:    Level: advanced

248: .keywords: MG, get, multigrid, restriction, level

250: .seealso: PCMGGetInterpolation(), PCMGSetRestriction(), PCMGGetRScale()
251: @*/
252: PetscErrorCode  PCMGGetRestriction(PC pc,PetscInt l,Mat *mat)
253: {
254:   PC_MG          *mg        = (PC_MG*)pc->data;
255:   PC_MG_Levels   **mglevels = mg->levels;

261:   if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
262:   if (l <= 0 || mg->nlevels <= l) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Level %D must be in range {1,...,%D}",l,mg->nlevels-1);
263:   if (!mglevels[l]->restrct) {
264:     if (!mglevels[l]->interpolate) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call PCMGSetRestriction() or PCMGSetInterpolation()");
265:     PCMGSetRestriction(pc,l,mglevels[l]->interpolate);
266:   }
267:   *mat = mglevels[l]->restrct;
268:   return(0);
269: }

273: /*@
274:    PCMGSetRScale - Sets the pointwise scaling for the restriction operator from level l to l-1.

276:    Logically Collective on PC and Vec

278:    Input Parameters:
279: +  pc - the multigrid context
280: -  l - the level (0 is coarsest) to supply [Do not supply 0]
281: .  rscale - the scaling

283:    Level: advanced

285:    Notes:
286:        When evaluating a function on a coarse level one does not want to do F(R * x) one does F(rscale * R * x) where rscale is 1 over the row sums of R.

288: .keywords: MG, set, multigrid, restriction, level

290: .seealso: PCMGSetInterpolation(), PCMGSetRestriction(), PCMGGetRScale()
291: @*/
292: PetscErrorCode  PCMGSetRScale(PC pc,PetscInt l,Vec rscale)
293: {
295:   PC_MG          *mg        = (PC_MG*)pc->data;
296:   PC_MG_Levels   **mglevels = mg->levels;

300:   if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
301:   if (l <= 0 || mg->nlevels <= l) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Level %D must be in range {1,...,%D}",l,mg->nlevels-1);
302:   PetscObjectReference((PetscObject)rscale);
303:   VecDestroy(&mglevels[l]->rscale);

305:   mglevels[l]->rscale = rscale;
306:   return(0);
307: }

311: /*@
312:    PCMGGetRScale - Gets the pointwise scaling for the restriction operator from level l to l-1.

314:    Collective on PC

316:    Input Parameters:
317: +  pc - the multigrid context
318: .  rscale - the scaling
319: -  l - the level (0 is coarsest) to supply [Do not supply 0]

321:    Level: advanced

323:    Notes:
324:        When evaluating a function on a coarse level one does not want to do F(R * x) one does F(rscale * R * x) where rscale is 1 over the row sums of R.

326: .keywords: MG, set, multigrid, restriction, level

328: .seealso: PCMGSetInterpolation(), PCMGGetRestriction()
329: @*/
330: PetscErrorCode PCMGGetRScale(PC pc,PetscInt l,Vec *rscale)
331: {
333:   PC_MG          *mg        = (PC_MG*)pc->data;
334:   PC_MG_Levels   **mglevels = mg->levels;

338:   if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
339:   if (l <= 0 || mg->nlevels <= l) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Level %D must be in range {1,...,%D}",l,mg->nlevels-1);
340:   if (!mglevels[l]->rscale) {
341:     Mat      R;
342:     Vec      X,Y,coarse,fine;
343:     PetscInt M,N;
344:     PCMGGetRestriction(pc,l,&R);
345:     MatCreateVecs(R,&X,&Y);
346:     MatGetSize(R,&M,&N);
347:     if (M < N) {
348:       fine = X;
349:       coarse = Y;
350:     } else if (N < M) {
351:       fine = Y; coarse = X;
352:     } else SETERRQ(PetscObjectComm((PetscObject)R),PETSC_ERR_SUP,"Restriction matrix is square, cannot determine which Vec is coarser");
353:     VecSet(fine,1.);
354:     MatRestrict(R,fine,coarse);
355:     VecDestroy(&fine);
356:     VecReciprocal(coarse);
357:     mglevels[l]->rscale = coarse;
358:   }
359:   *rscale = mglevels[l]->rscale;
360:   return(0);
361: }

365: /*@
366:    PCMGGetSmoother - Gets the KSP context to be used as smoother for
367:    both pre- and post-smoothing.  Call both PCMGGetSmootherUp() and
368:    PCMGGetSmootherDown() to use different functions for pre- and
369:    post-smoothing.

371:    Not Collective, KSP returned is parallel if PC is

373:    Input Parameters:
374: +  pc - the multigrid context
375: -  l - the level (0 is coarsest) to supply

377:    Ouput Parameters:
378: .  ksp - the smoother

380:    Notes:
381:    Once you have called this routine, you can call KSPSetOperators(ksp,...) on the resulting ksp to provide the operators for the smoother for this level.
382:    You can also modify smoother options by calling the various KSPSetXXX() options on this ksp. In addition you can call KSPGetPC(ksp,&pc)
383:    and modify PC options for the smoother; for example PCSetType(pc,PCSOR); to use SOR smoothing.

385:    Level: advanced

387: .keywords: MG, get, multigrid, level, smoother, pre-smoother, post-smoother

389: .seealso: PCMGGetSmootherUp(), PCMGGetSmootherDown()
390: @*/
391: PetscErrorCode  PCMGGetSmoother(PC pc,PetscInt l,KSP *ksp)
392: {
393:   PC_MG        *mg        = (PC_MG*)pc->data;
394:   PC_MG_Levels **mglevels = mg->levels;

398:   *ksp = mglevels[l]->smoothd;
399:   return(0);
400: }

404: /*@
405:    PCMGGetSmootherUp - Gets the KSP context to be used as smoother after
406:    coarse grid correction (post-smoother).

408:    Not Collective, KSP returned is parallel if PC is

410:    Input Parameters:
411: +  pc - the multigrid context
412: -  l  - the level (0 is coarsest) to supply

414:    Ouput Parameters:
415: .  ksp - the smoother

417:    Level: advanced

419:    Notes: calling this will result in a different pre and post smoother so you may need to
420:          set options on the pre smoother also

422: .keywords: MG, multigrid, get, smoother, up, post-smoother, level

424: .seealso: PCMGGetSmootherUp(), PCMGGetSmootherDown()
425: @*/
426: PetscErrorCode  PCMGGetSmootherUp(PC pc,PetscInt l,KSP *ksp)
427: {
428:   PC_MG          *mg        = (PC_MG*)pc->data;
429:   PC_MG_Levels   **mglevels = mg->levels;
431:   const char     *prefix;
432:   MPI_Comm       comm;

436:   /*
437:      This is called only if user wants a different pre-smoother from post.
438:      Thus we check if a different one has already been allocated,
439:      if not we allocate it.
440:   */
441:   if (!l) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"There is no such thing as a up smoother on the coarse grid");
442:   if (mglevels[l]->smoothu == mglevels[l]->smoothd) {
443:     KSPType     ksptype;
444:     PCType      pctype;
445:     PC          ipc;
446:     PetscReal   rtol,abstol,dtol;
447:     PetscInt    maxits;
448:     KSPNormType normtype;
449:     PetscObjectGetComm((PetscObject)mglevels[l]->smoothd,&comm);
450:     KSPGetOptionsPrefix(mglevels[l]->smoothd,&prefix);
451:     KSPGetTolerances(mglevels[l]->smoothd,&rtol,&abstol,&dtol,&maxits);
452:     KSPGetType(mglevels[l]->smoothd,&ksptype);
453:     KSPGetNormType(mglevels[l]->smoothd,&normtype);
454:     KSPGetPC(mglevels[l]->smoothd,&ipc);
455:     PCGetType(ipc,&pctype);

457:     KSPCreate(comm,&mglevels[l]->smoothu);
458:     KSPSetErrorIfNotConverged(mglevels[l]->smoothu,pc->erroriffailure);
459:     PetscObjectIncrementTabLevel((PetscObject)mglevels[l]->smoothu,(PetscObject)pc,mglevels[0]->levels-l);
460:     KSPSetOptionsPrefix(mglevels[l]->smoothu,prefix);
461:     KSPSetTolerances(mglevels[l]->smoothu,rtol,abstol,dtol,maxits);
462:     KSPSetType(mglevels[l]->smoothu,ksptype);
463:     KSPSetNormType(mglevels[l]->smoothu,normtype);
464:     KSPSetConvergenceTest(mglevels[l]->smoothu,KSPConvergedSkip,NULL,NULL);
465:     KSPGetPC(mglevels[l]->smoothu,&ipc);
466:     PCSetType(ipc,pctype);
467:     PetscLogObjectParent((PetscObject)pc,(PetscObject)mglevels[l]->smoothu);
468:   }
469:   if (ksp) *ksp = mglevels[l]->smoothu;
470:   return(0);
471: }

475: /*@
476:    PCMGGetSmootherDown - Gets the KSP context to be used as smoother before
477:    coarse grid correction (pre-smoother).

479:    Not Collective, KSP returned is parallel if PC is

481:    Input Parameters:
482: +  pc - the multigrid context
483: -  l  - the level (0 is coarsest) to supply

485:    Ouput Parameters:
486: .  ksp - the smoother

488:    Level: advanced

490:    Notes: calling this will result in a different pre and post smoother so you may need to
491:          set options on the post smoother also

493: .keywords: MG, multigrid, get, smoother, down, pre-smoother, level

495: .seealso: PCMGGetSmootherUp(), PCMGGetSmoother()
496: @*/
497: PetscErrorCode  PCMGGetSmootherDown(PC pc,PetscInt l,KSP *ksp)
498: {
500:   PC_MG          *mg        = (PC_MG*)pc->data;
501:   PC_MG_Levels   **mglevels = mg->levels;

505:   /* make sure smoother up and down are different */
506:   if (l) {
507:     PCMGGetSmootherUp(pc,l,NULL);
508:   }
509:   *ksp = mglevels[l]->smoothd;
510:   return(0);
511: }

515: /*@
516:    PCMGSetCyclesOnLevel - Sets the number of cycles to run on this level.

518:    Logically Collective on PC

520:    Input Parameters:
521: +  pc - the multigrid context
522: .  l  - the level (0 is coarsest) this is to be used for
523: -  n  - the number of cycles

525:    Level: advanced

527: .keywords: MG, multigrid, set, cycles, V-cycle, W-cycle, level

529: .seealso: PCMGSetCycles()
530: @*/
531: PetscErrorCode  PCMGSetCyclesOnLevel(PC pc,PetscInt l,PetscInt c)
532: {
533:   PC_MG        *mg        = (PC_MG*)pc->data;
534:   PC_MG_Levels **mglevels = mg->levels;

538:   if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
541:   mglevels[l]->cycles = c;
542:   return(0);
543: }

547: /*@
548:    PCMGSetRhs - Sets the vector space to be used to store the right-hand side
549:    on a particular level.

551:    Logically Collective on PC and Vec

553:    Input Parameters:
554: +  pc - the multigrid context
555: .  l  - the level (0 is coarsest) this is to be used for
556: -  c  - the space

558:    Level: advanced

560:    Notes: If this is not provided PETSc will automatically generate one.

562:           You do not need to keep a reference to this vector if you do
563:           not need it PCDestroy() will properly free it.

565: .keywords: MG, multigrid, set, right-hand-side, rhs, level

567: .seealso: PCMGSetX(), PCMGSetR()
568: @*/
569: PetscErrorCode  PCMGSetRhs(PC pc,PetscInt l,Vec c)
570: {
572:   PC_MG          *mg        = (PC_MG*)pc->data;
573:   PC_MG_Levels   **mglevels = mg->levels;

577:   if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
578:   if (l == mglevels[0]->levels-1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Do not set rhs for finest level");
579:   PetscObjectReference((PetscObject)c);
580:   VecDestroy(&mglevels[l]->b);

582:   mglevels[l]->b = c;
583:   return(0);
584: }

588: /*@
589:    PCMGSetX - Sets the vector space to be used to store the solution on a
590:    particular level.

592:    Logically Collective on PC and Vec

594:    Input Parameters:
595: +  pc - the multigrid context
596: .  l - the level (0 is coarsest) this is to be used for (do not supply the finest level)
597: -  c - the space

599:    Level: advanced

601:    Notes: If this is not provided PETSc will automatically generate one.

603:           You do not need to keep a reference to this vector if you do
604:           not need it PCDestroy() will properly free it.

606: .keywords: MG, multigrid, set, solution, level

608: .seealso: PCMGSetRhs(), PCMGSetR()
609: @*/
610: PetscErrorCode  PCMGSetX(PC pc,PetscInt l,Vec c)
611: {
613:   PC_MG          *mg        = (PC_MG*)pc->data;
614:   PC_MG_Levels   **mglevels = mg->levels;

618:   if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
619:   if (l == mglevels[0]->levels-1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Do not set x for finest level");
620:   PetscObjectReference((PetscObject)c);
621:   VecDestroy(&mglevels[l]->x);

623:   mglevels[l]->x = c;
624:   return(0);
625: }

629: /*@
630:    PCMGSetR - Sets the vector space to be used to store the residual on a
631:    particular level.

633:    Logically Collective on PC and Vec

635:    Input Parameters:
636: +  pc - the multigrid context
637: .  l - the level (0 is coarsest) this is to be used for
638: -  c - the space

640:    Level: advanced

642:    Notes: If this is not provided PETSc will automatically generate one.

644:           You do not need to keep a reference to this vector if you do
645:           not need it PCDestroy() will properly free it.

647: .keywords: MG, multigrid, set, residual, level
648: @*/
649: PetscErrorCode  PCMGSetR(PC pc,PetscInt l,Vec c)
650: {
652:   PC_MG          *mg        = (PC_MG*)pc->data;
653:   PC_MG_Levels   **mglevels = mg->levels;

657:   if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
658:   if (!l) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Need not set residual vector for coarse grid");
659:   PetscObjectReference((PetscObject)c);
660:   VecDestroy(&mglevels[l]->r);

662:   mglevels[l]->r = c;
663:   return(0);
664: }