Actual source code: mgfunc.c

petsc-3.4.4 2014-03-13
  2: #include <../src/ksp/pc/impls/mg/mgimpl.h>       /*I "petscksp.h" I*/

  6: /*@C
  7:    PCMGResidual_Default - Default routine to calculate the residual.

  9:    Collective on Mat and Vec

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

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

 19:    Level: developer

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

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

 30:   MatMult(mat,x,r);
 31:   VecAYPX(r,-1.0,b);
 32:   return(0);
 33: }

 35: /* ---------------------------------------------------------------------------*/

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

 42:    Not Collective

 44:    Input Parameter:
 45: .  pc - the multigrid context

 47:    Output Parameter:
 48: .  ksp - the coarse grid solver context

 50:    Level: advanced

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

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

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

 71:    Logically Collective on PC and Mat

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

 80:    Level: advanced

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

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

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

100:   mglevels[l]->A = mat;
101:   return(0);
102: }

106: /*@
107:    PCMGSetInterpolation - Sets the function to be used to calculate the
108:    interpolation from l-1 to the lth level

110:    Logically Collective on PC and Mat

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

117:    Level: advanced

119:    Notes:
120:           Usually this is the same matrix used also to set the restriction
121:     for the same level.

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

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

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

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

143:   mglevels[l]->interpolate = mat;
144:   return(0);
145: }

149: /*@
150:    PCMGGetInterpolation - Gets the function to be used to calculate the
151:    interpolation from l-1 to the lth level

153:    Logically Collective on PC

155:    Input Parameters:
156: +  pc - the multigrid context
157: -  l - the level (0 is coarsest) to supply [Do not supply 0]

159:    Output Parameter:
160: .  mat - the interpolation matrix

162:    Level: advanced

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

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

177:   if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
178:   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);
179:   if (!mglevels[l]->interpolate) {
180:     if (!mglevels[l]->restrct) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call PCMGSetInterpolation() or PCMGSetInterpolation()");
181:     PCMGSetInterpolation(pc,l,mglevels[l]->restrct);
182:   }
183:   *mat = mglevels[l]->interpolate;
184:   return(0);
185: }

189: /*@
190:    PCMGSetRestriction - Sets the function to be used to restrict vector
191:    from level l to l-1.

193:    Logically Collective on PC and Mat

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

200:    Level: advanced

202:    Notes:
203:           Usually this is the same matrix used also to set the interpolation
204:     for the same level.

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

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

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

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

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

230:   mglevels[l]->restrct = mat;
231:   return(0);
232: }

236: /*@
237:    PCMGGetRestriction - Gets the function to be used to restrict vector
238:    from level l to l-1.

240:    Logically Collective on PC and Mat

242:    Input Parameters:
243: +  pc - the multigrid context
244: -  l - the level (0 is coarsest) to supply [Do not supply 0]

246:    Output Parameter:
247: .  mat - the restriction matrix

249:    Level: advanced

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

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

264:   if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
265:   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);
266:   if (!mglevels[l]->restrct) {
267:     if (!mglevels[l]->interpolate) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call PCMGSetRestriction() or PCMGSetInterpolation()");
268:     PCMGSetRestriction(pc,l,mglevels[l]->interpolate);
269:   }
270:   *mat = mglevels[l]->restrct;
271:   return(0);
272: }

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

279:    Logically Collective on PC and Vec

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

286:    Level: advanced

288:    Notes:
289:        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.

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

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

303:   if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
304:   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);
305:   PetscObjectReference((PetscObject)rscale);
306:   VecDestroy(&mglevels[l]->rscale);

308:   mglevels[l]->rscale = rscale;
309:   return(0);
310: }

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

317:    Collective on PC

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

324:    Level: advanced

326:    Notes:
327:        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.

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

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

341:   if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
342:   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);
343:   if (!mglevels[l]->rscale) {
344:     Mat      R;
345:     Vec      X,Y,coarse,fine;
346:     PetscInt M,N;
347:     PCMGGetRestriction(pc,l,&R);
348:     MatGetVecs(R,&X,&Y);
349:     MatGetSize(R,&M,&N);
350:     if (M < N) {
351:       fine = X;
352:       coarse = Y;
353:     } else if (N < M) {
354:       fine = Y; coarse = X;
355:     } else SETERRQ(PetscObjectComm((PetscObject)R),PETSC_ERR_SUP,"Restriction matrix is square, cannot determine which Vec is coarser");
356:     VecSet(fine,1.);
357:     MatRestrict(R,fine,coarse);
358:     VecDestroy(&fine);
359:     VecReciprocal(coarse);
360:     mglevels[l]->rscale = coarse;
361:   }
362:   *rscale = mglevels[l]->rscale;
363:   return(0);
364: }

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

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

376:    Input Parameters:
377: +  pc - the multigrid context
378: -  l - the level (0 is coarsest) to supply

380:    Ouput Parameters:
381: .  ksp - the smoother

383:    Level: advanced

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

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

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

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

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

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

412:    Ouput Parameters:
413: .  ksp - the smoother

415:    Level: advanced

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

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

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

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

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

472: /*@
473:    PCMGGetSmootherDown - Gets the KSP context to be used as smoother before
474:    coarse grid correction (pre-smoother).

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

478:    Input Parameters:
479: +  pc - the multigrid context
480: -  l  - the level (0 is coarsest) to supply

482:    Ouput Parameters:
483: .  ksp - the smoother

485:    Level: advanced

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

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

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

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

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

515:    Logically Collective on PC

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

522:    Level: advanced

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

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

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

544: /*@
545:    PCMGSetRhs - Sets the vector space to be used to store the right-hand side
546:    on a particular level.

548:    Logically Collective on PC and Vec

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

555:    Level: advanced

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

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

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

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

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

579:   mglevels[l]->b = c;
580:   return(0);
581: }

585: /*@
586:    PCMGSetX - Sets the vector space to be used to store the solution on a
587:    particular level.

589:    Logically Collective on PC and Vec

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

596:    Level: advanced

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

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

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

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

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

620:   mglevels[l]->x = c;
621:   return(0);
622: }

626: /*@
627:    PCMGSetR - Sets the vector space to be used to store the residual on a
628:    particular level.

630:    Logically Collective on PC and Vec

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

637:    Level: advanced

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

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

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

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

659:   mglevels[l]->r = c;
660:   return(0);
661: }