Actual source code: mgfunc.c

petsc-dev 2014-08-28
Report Typos and Errors
  2: #include <../src/ksp/pc/impls/mg/mgimpl.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:     MatGetVecs(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:    Level: advanced

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

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

393:   *ksp = mglevels[l]->smoothd;
394:   return(0);
395: }

399: /*@
400:    PCMGGetSmootherUp - Gets the KSP context to be used as smoother after
401:    coarse grid correction (post-smoother).

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

405:    Input Parameters:
406: +  pc - the multigrid context
407: -  l  - the level (0 is coarsest) to supply

409:    Ouput Parameters:
410: .  ksp - the smoother

412:    Level: advanced

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

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

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

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

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

469: /*@
470:    PCMGGetSmootherDown - Gets the KSP context to be used as smoother before
471:    coarse grid correction (pre-smoother).

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

475:    Input Parameters:
476: +  pc - the multigrid context
477: -  l  - the level (0 is coarsest) to supply

479:    Ouput Parameters:
480: .  ksp - the smoother

482:    Level: advanced

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

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

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

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

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

512:    Logically Collective on PC

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

519:    Level: advanced

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

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

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

541: /*@
542:    PCMGSetRhs - Sets the vector space to be used to store the right-hand side
543:    on a particular level.

545:    Logically Collective on PC and Vec

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

552:    Level: advanced

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

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

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

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

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

576:   mglevels[l]->b = c;
577:   return(0);
578: }

582: /*@
583:    PCMGSetX - Sets the vector space to be used to store the solution on a
584:    particular level.

586:    Logically Collective on PC and Vec

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

593:    Level: advanced

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

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

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

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

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

617:   mglevels[l]->x = c;
618:   return(0);
619: }

623: /*@
624:    PCMGSetR - Sets the vector space to be used to store the residual on a
625:    particular level.

627:    Logically Collective on PC and Vec

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

634:    Level: advanced

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

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

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

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

656:   mglevels[l]->r = c;
657:   return(0);
658: }