Actual source code: mgfunc.c

petsc-master 2017-05-26
Report Typos and Errors

  2:  #include <petsc/private/pcmgimpl.h>

  4: /* ---------------------------------------------------------------------------*/
  5: /*@C
  6:    PCMGResidualDefault - Default routine to calculate the residual.

  8:    Collective on Mat and Vec

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

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

 18:    Level: developer

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

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

 29:   MatResidual(mat,b,x,r);
 30:   return(0);
 31: }

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

 36:    Not Collective

 38:    Input Parameter:
 39: .  pc - the multigrid context

 41:    Output Parameter:
 42: .  ksp - the coarse grid solver context

 44:    Level: advanced

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

 55:   *ksp =  mglevels[0]->smoothd;
 56:   return(0);
 57: }

 59: /*@C
 60:    PCMGSetResidual - Sets the function to be used to calculate the residual
 61:    on the lth level.

 63:    Logically Collective on PC and Mat

 65:    Input Parameters:
 66: +  pc       - the multigrid context
 67: .  l        - the level (0 is coarsest) to supply
 68: .  residual - function used to form residual, if none is provided the previously provide one is used, if no
 69:               previous one were provided then a default is used
 70: -  mat      - matrix associated with residual

 72:    Level: advanced

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

 76: .seealso: PCMGResidualDefault()
 77: @*/
 78: PetscErrorCode  PCMGSetResidual(PC pc,PetscInt l,PetscErrorCode (*residual)(Mat,Vec,Vec,Vec),Mat mat)
 79: {
 80:   PC_MG          *mg        = (PC_MG*)pc->data;
 81:   PC_MG_Levels   **mglevels = mg->levels;

 86:   if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
 87:   if (residual) mglevels[l]->residual = residual;
 88:   if (!mglevels[l]->residual) mglevels[l]->residual = PCMGResidualDefault;
 89:   if (mat) {PetscObjectReference((PetscObject)mat);}
 90:   MatDestroy(&mglevels[l]->A);
 91:   mglevels[l]->A = mat;
 92:   return(0);
 93: }

 95: /*@
 96:    PCMGSetInterpolation - Sets the function to be used to calculate the
 97:    interpolation from l-1 to the lth level

 99:    Logically Collective on PC and Mat

101:    Input Parameters:
102: +  pc  - the multigrid context
103: .  mat - the interpolation operator
104: -  l   - the level (0 is coarsest) to supply [do not supply 0]

106:    Level: advanced

108:    Notes:
109:           Usually this is the same matrix used also to set the restriction
110:     for the same level.

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

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

117: .seealso: PCMGSetRestriction()
118: @*/
119: PetscErrorCode  PCMGSetInterpolation(PC pc,PetscInt l,Mat mat)
120: {
121:   PC_MG          *mg        = (PC_MG*)pc->data;
122:   PC_MG_Levels   **mglevels = mg->levels;

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

132:   mglevels[l]->interpolate = mat;
133:   return(0);
134: }

136: /*@
137:    PCMGGetInterpolation - Gets the function to be used to calculate the
138:    interpolation from l-1 to the lth level

140:    Logically Collective on PC

142:    Input Parameters:
143: +  pc - the multigrid context
144: -  l - the level (0 is coarsest) to supply [Do not supply 0]

146:    Output Parameter:
147: .  mat - the interpolation matrix, can be NULL

149:    Level: advanced

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

153: .seealso: PCMGGetRestriction(), PCMGSetInterpolation(), PCMGGetRScale()
154: @*/
155: PetscErrorCode  PCMGGetInterpolation(PC pc,PetscInt l,Mat *mat)
156: {
157:   PC_MG          *mg        = (PC_MG*)pc->data;
158:   PC_MG_Levels   **mglevels = mg->levels;

164:   if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
165:   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);
166:   if (!mglevels[l]->interpolate) {
167:     if (!mglevels[l]->restrct) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call PCMGSetInterpolation() or PCMGSetRestriction()");
168:     PCMGSetInterpolation(pc,l,mglevels[l]->restrct);
169:   }
170:   if (mat) *mat = mglevels[l]->interpolate;
171:   return(0);
172: }

174: /*@
175:    PCMGSetRestriction - Sets the function to be used to restrict vector
176:    from level l to l-1.

178:    Logically Collective on PC and Mat

180:    Input Parameters:
181: +  pc - the multigrid context
182: .  l - the level (0 is coarsest) to supply [Do not supply 0]
183: -  mat - the restriction matrix

185:    Level: advanced

187:    Notes:
188:           Usually this is the same matrix used also to set the interpolation
189:     for the same level.

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

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

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

199: .seealso: PCMGSetInterpolation()
200: @*/
201: PetscErrorCode  PCMGSetRestriction(PC pc,PetscInt l,Mat mat)
202: {
204:   PC_MG          *mg        = (PC_MG*)pc->data;
205:   PC_MG_Levels   **mglevels = mg->levels;

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

215:   mglevels[l]->restrct = mat;
216:   return(0);
217: }

219: /*@
220:    PCMGGetRestriction - Gets the function to be used to restrict vector
221:    from level l to l-1.

223:    Logically Collective on PC and Mat

225:    Input Parameters:
226: +  pc - the multigrid context
227: -  l - the level (0 is coarsest) to supply [Do not supply 0]

229:    Output Parameter:
230: .  mat - the restriction matrix

232:    Level: advanced

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

236: .seealso: PCMGGetInterpolation(), PCMGSetRestriction(), PCMGGetRScale()
237: @*/
238: PetscErrorCode  PCMGGetRestriction(PC pc,PetscInt l,Mat *mat)
239: {
240:   PC_MG          *mg        = (PC_MG*)pc->data;
241:   PC_MG_Levels   **mglevels = mg->levels;

247:   if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
248:   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);
249:   if (!mglevels[l]->restrct) {
250:     if (!mglevels[l]->interpolate) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call PCMGSetRestriction() or PCMGSetInterpolation()");
251:     PCMGSetRestriction(pc,l,mglevels[l]->interpolate);
252:   }
253:   if (mat) *mat = mglevels[l]->restrct;
254:   return(0);
255: }

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

260:    Logically Collective on PC and Vec

262:    Input Parameters:
263: +  pc - the multigrid context
264: -  l - the level (0 is coarsest) to supply [Do not supply 0]
265: .  rscale - the scaling

267:    Level: advanced

269:    Notes:
270:        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.

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

274: .seealso: PCMGSetInterpolation(), PCMGSetRestriction(), PCMGGetRScale()
275: @*/
276: PetscErrorCode  PCMGSetRScale(PC pc,PetscInt l,Vec rscale)
277: {
279:   PC_MG          *mg        = (PC_MG*)pc->data;
280:   PC_MG_Levels   **mglevels = mg->levels;

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

289:   mglevels[l]->rscale = rscale;
290:   return(0);
291: }

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

296:    Collective on PC

298:    Input Parameters:
299: +  pc - the multigrid context
300: .  rscale - the scaling
301: -  l - the level (0 is coarsest) to supply [Do not supply 0]

303:    Level: advanced

305:    Notes:
306:        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.

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

310: .seealso: PCMGSetInterpolation(), PCMGGetRestriction()
311: @*/
312: PetscErrorCode PCMGGetRScale(PC pc,PetscInt l,Vec *rscale)
313: {
315:   PC_MG          *mg        = (PC_MG*)pc->data;
316:   PC_MG_Levels   **mglevels = mg->levels;

320:   if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
321:   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);
322:   if (!mglevels[l]->rscale) {
323:     Mat      R;
324:     Vec      X,Y,coarse,fine;
325:     PetscInt M,N;
326:     PCMGGetRestriction(pc,l,&R);
327:     MatCreateVecs(R,&X,&Y);
328:     MatGetSize(R,&M,&N);
329:     if (M < N) {
330:       fine = X;
331:       coarse = Y;
332:     } else if (N < M) {
333:       fine = Y; coarse = X;
334:     } else SETERRQ(PetscObjectComm((PetscObject)R),PETSC_ERR_SUP,"Restriction matrix is square, cannot determine which Vec is coarser");
335:     VecSet(fine,1.);
336:     MatRestrict(R,fine,coarse);
337:     VecDestroy(&fine);
338:     VecReciprocal(coarse);
339:     mglevels[l]->rscale = coarse;
340:   }
341:   *rscale = mglevels[l]->rscale;
342:   return(0);
343: }

345: /*@
346:    PCMGGetSmoother - Gets the KSP context to be used as smoother for
347:    both pre- and post-smoothing.  Call both PCMGGetSmootherUp() and
348:    PCMGGetSmootherDown() to use different functions for pre- and
349:    post-smoothing.

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

353:    Input Parameters:
354: +  pc - the multigrid context
355: -  l - the level (0 is coarsest) to supply

357:    Ouput Parameters:
358: .  ksp - the smoother

360:    Notes:
361:    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.
362:    You can also modify smoother options by calling the various KSPSetXXX() options on this ksp. In addition you can call KSPGetPC(ksp,&pc)
363:    and modify PC options for the smoother; for example PCSetType(pc,PCSOR); to use SOR smoothing.

365:    Level: advanced

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

369: .seealso: PCMGGetSmootherUp(), PCMGGetSmootherDown()
370: @*/
371: PetscErrorCode  PCMGGetSmoother(PC pc,PetscInt l,KSP *ksp)
372: {
373:   PC_MG        *mg        = (PC_MG*)pc->data;
374:   PC_MG_Levels **mglevels = mg->levels;

378:   *ksp = mglevels[l]->smoothd;
379:   return(0);
380: }

382: /*@
383:    PCMGGetSmootherUp - Gets the KSP context to be used as smoother after
384:    coarse grid correction (post-smoother).

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

388:    Input Parameters:
389: +  pc - the multigrid context
390: -  l  - the level (0 is coarsest) to supply

392:    Ouput Parameters:
393: .  ksp - the smoother

395:    Level: advanced

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

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

402: .seealso: PCMGGetSmootherUp(), PCMGGetSmootherDown()
403: @*/
404: PetscErrorCode  PCMGGetSmootherUp(PC pc,PetscInt l,KSP *ksp)
405: {
406:   PC_MG          *mg        = (PC_MG*)pc->data;
407:   PC_MG_Levels   **mglevels = mg->levels;
409:   const char     *prefix;
410:   MPI_Comm       comm;

414:   /*
415:      This is called only if user wants a different pre-smoother from post.
416:      Thus we check if a different one has already been allocated,
417:      if not we allocate it.
418:   */
419:   if (!l) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"There is no such thing as a up smoother on the coarse grid");
420:   if (mglevels[l]->smoothu == mglevels[l]->smoothd) {
421:     KSPType     ksptype;
422:     PCType      pctype;
423:     PC          ipc;
424:     PetscReal   rtol,abstol,dtol;
425:     PetscInt    maxits;
426:     KSPNormType normtype;
427:     PetscObjectGetComm((PetscObject)mglevels[l]->smoothd,&comm);
428:     KSPGetOptionsPrefix(mglevels[l]->smoothd,&prefix);
429:     KSPGetTolerances(mglevels[l]->smoothd,&rtol,&abstol,&dtol,&maxits);
430:     KSPGetType(mglevels[l]->smoothd,&ksptype);
431:     KSPGetNormType(mglevels[l]->smoothd,&normtype);
432:     KSPGetPC(mglevels[l]->smoothd,&ipc);
433:     PCGetType(ipc,&pctype);

435:     KSPCreate(comm,&mglevels[l]->smoothu);
436:     KSPSetErrorIfNotConverged(mglevels[l]->smoothu,pc->erroriffailure);
437:     PetscObjectIncrementTabLevel((PetscObject)mglevels[l]->smoothu,(PetscObject)pc,mglevels[0]->levels-l);
438:     KSPSetOptionsPrefix(mglevels[l]->smoothu,prefix);
439:     KSPSetTolerances(mglevels[l]->smoothu,rtol,abstol,dtol,maxits);
440:     KSPSetType(mglevels[l]->smoothu,ksptype);
441:     KSPSetNormType(mglevels[l]->smoothu,normtype);
442:     KSPSetConvergenceTest(mglevels[l]->smoothu,KSPConvergedSkip,NULL,NULL);
443:     KSPGetPC(mglevels[l]->smoothu,&ipc);
444:     PCSetType(ipc,pctype);
445:     PetscLogObjectParent((PetscObject)pc,(PetscObject)mglevels[l]->smoothu);
446:     PetscObjectComposedDataSetInt((PetscObject) mglevels[l]->smoothu, PetscMGLevelId, mglevels[l]->level);
447:   }
448:   if (ksp) *ksp = mglevels[l]->smoothu;
449:   return(0);
450: }

452: /*@
453:    PCMGGetSmootherDown - Gets the KSP context to be used as smoother before
454:    coarse grid correction (pre-smoother).

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

458:    Input Parameters:
459: +  pc - the multigrid context
460: -  l  - the level (0 is coarsest) to supply

462:    Ouput Parameters:
463: .  ksp - the smoother

465:    Level: advanced

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

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

472: .seealso: PCMGGetSmootherUp(), PCMGGetSmoother()
473: @*/
474: PetscErrorCode  PCMGGetSmootherDown(PC pc,PetscInt l,KSP *ksp)
475: {
477:   PC_MG          *mg        = (PC_MG*)pc->data;
478:   PC_MG_Levels   **mglevels = mg->levels;

482:   /* make sure smoother up and down are different */
483:   if (l) {
484:     PCMGGetSmootherUp(pc,l,NULL);
485:   }
486:   *ksp = mglevels[l]->smoothd;
487:   return(0);
488: }

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

493:    Logically Collective on PC

495:    Input Parameters:
496: +  pc - the multigrid context
497: .  l  - the level (0 is coarsest) this is to be used for
498: -  n  - the number of cycles

500:    Level: advanced

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

504: .seealso: PCMGSetCycles()
505: @*/
506: PetscErrorCode  PCMGSetCyclesOnLevel(PC pc,PetscInt l,PetscInt c)
507: {
508:   PC_MG        *mg        = (PC_MG*)pc->data;
509:   PC_MG_Levels **mglevels = mg->levels;

513:   if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
516:   mglevels[l]->cycles = c;
517:   return(0);
518: }

520: /*@
521:    PCMGSetRhs - Sets the vector space to be used to store the right-hand side
522:    on a particular level.

524:    Logically Collective on PC and Vec

526:    Input Parameters:
527: +  pc - the multigrid context
528: .  l  - the level (0 is coarsest) this is to be used for
529: -  c  - the space

531:    Level: advanced

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

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

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

540: .seealso: PCMGSetX(), PCMGSetR()
541: @*/
542: PetscErrorCode  PCMGSetRhs(PC pc,PetscInt l,Vec c)
543: {
545:   PC_MG          *mg        = (PC_MG*)pc->data;
546:   PC_MG_Levels   **mglevels = mg->levels;

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

555:   mglevels[l]->b = c;
556:   return(0);
557: }

559: /*@
560:    PCMGSetX - Sets the vector space to be used to store the solution on a
561:    particular level.

563:    Logically Collective on PC and Vec

565:    Input Parameters:
566: +  pc - the multigrid context
567: .  l - the level (0 is coarsest) this is to be used for (do not supply the finest level)
568: -  c - the space

570:    Level: advanced

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

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

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

579: .seealso: PCMGSetRhs(), PCMGSetR()
580: @*/
581: PetscErrorCode  PCMGSetX(PC pc,PetscInt l,Vec c)
582: {
584:   PC_MG          *mg        = (PC_MG*)pc->data;
585:   PC_MG_Levels   **mglevels = mg->levels;

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

594:   mglevels[l]->x = c;
595:   return(0);
596: }

598: /*@
599:    PCMGSetR - Sets the vector space to be used to store the residual on a
600:    particular level.

602:    Logically Collective on PC and Vec

604:    Input Parameters:
605: +  pc - the multigrid context
606: .  l - the level (0 is coarsest) this is to be used for
607: -  c - the space

609:    Level: advanced

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

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

616: .keywords: MG, multigrid, set, residual, level
617: @*/
618: PetscErrorCode  PCMGSetR(PC pc,PetscInt l,Vec c)
619: {
621:   PC_MG          *mg        = (PC_MG*)pc->data;
622:   PC_MG_Levels   **mglevels = mg->levels;

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

631:   mglevels[l]->r = c;
632:   return(0);
633: }