Actual source code: mgfunc.c
2: #include <petsc/private/pcmgimpl.h>
4: /* ---------------------------------------------------------------------------*/
5: /*@C
6: PCMGResidualDefault - Default routine to calculate the residual.
8: Collective on Mat
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: .seealso: PCMGSetResidual()
21: @*/
22: PetscErrorCode PCMGResidualDefault(Mat mat,Vec b,Vec x,Vec r)
23: {
24: MatResidual(mat,b,x,r);
25: return 0;
26: }
28: /*@C
29: PCMGResidualTransposeDefault - Default routine to calculate the residual of the transposed linear system
31: Collective on Mat
33: Input Parameters:
34: + mat - the matrix
35: . b - the right-hand-side
36: - x - the approximate solution
38: Output Parameter:
39: . r - location to store the residual
41: Level: developer
43: .seealso: PCMGSetResidualTranspose()
44: @*/
45: PetscErrorCode PCMGResidualTransposeDefault(Mat mat,Vec b,Vec x,Vec r)
46: {
47: MatMultTranspose(mat,x,r);
48: VecAYPX(r,-1.0,b);
49: return 0;
50: }
52: /*@C
53: PCMGMatResidualDefault - Default routine to calculate the residual.
55: Collective on Mat
57: Input Parameters:
58: + mat - the matrix
59: . b - the right-hand-side
60: - x - the approximate solution
62: Output Parameter:
63: . r - location to store the residual
65: Level: developer
67: .seealso: PCMGSetMatResidual()
68: @*/
69: PetscErrorCode PCMGMatResidualDefault(Mat mat,Mat b,Mat x,Mat r)
70: {
71: MatMatMult(mat,x,MAT_REUSE_MATRIX,PETSC_DEFAULT,&r);
72: MatAYPX(r,-1.0,b,UNKNOWN_NONZERO_PATTERN);
73: return 0;
74: }
76: /*@C
77: PCMGMatResidualTransposeDefault - Default routine to calculate the residual of the transposed linear system
79: Collective on Mat
81: Input Parameters:
82: + mat - the matrix
83: . b - the right-hand-side
84: - x - the approximate solution
86: Output Parameter:
87: . r - location to store the residual
89: Level: developer
91: .seealso: PCMGSetMatResidualTranspose()
92: @*/
93: PetscErrorCode PCMGMatResidualTransposeDefault(Mat mat,Mat b,Mat x,Mat r)
94: {
95: MatTransposeMatMult(mat,x,MAT_REUSE_MATRIX,PETSC_DEFAULT,&r);
96: MatAYPX(r,-1.0,b,UNKNOWN_NONZERO_PATTERN);
97: return 0;
98: }
99: /*@
100: PCMGGetCoarseSolve - Gets the solver context to be used on the coarse grid.
102: Not Collective
104: Input Parameter:
105: . pc - the multigrid context
107: Output Parameter:
108: . ksp - the coarse grid solver context
110: Level: advanced
112: .seealso: PCMGGetSmootherUp(), PCMGGetSmootherDown(), PCMGGetSmoother()
113: @*/
114: PetscErrorCode PCMGGetCoarseSolve(PC pc,KSP *ksp)
115: {
116: PC_MG *mg = (PC_MG*)pc->data;
117: PC_MG_Levels **mglevels = mg->levels;
120: *ksp = mglevels[0]->smoothd;
121: return 0;
122: }
124: /*@C
125: PCMGSetResidual - Sets the function to be used to calculate the residual
126: on the lth level.
128: Logically Collective on PC
130: Input Parameters:
131: + pc - the multigrid context
132: . l - the level (0 is coarsest) to supply
133: . residual - function used to form residual, if none is provided the previously provide one is used, if no
134: previous one were provided then a default is used
135: - mat - matrix associated with residual
137: Level: advanced
139: .seealso: PCMGResidualDefault()
140: @*/
141: PetscErrorCode PCMGSetResidual(PC pc,PetscInt l,PetscErrorCode (*residual)(Mat,Vec,Vec,Vec),Mat mat)
142: {
143: PC_MG *mg = (PC_MG*)pc->data;
144: PC_MG_Levels **mglevels = mg->levels;
148: if (residual) mglevels[l]->residual = residual;
149: if (!mglevels[l]->residual) mglevels[l]->residual = PCMGResidualDefault;
150: mglevels[l]->matresidual = PCMGMatResidualDefault;
151: if (mat) PetscObjectReference((PetscObject)mat);
152: MatDestroy(&mglevels[l]->A);
153: mglevels[l]->A = mat;
154: return 0;
155: }
157: /*@C
158: PCMGSetResidualTranspose - Sets the function to be used to calculate the residual of the transposed linear system
159: on the lth level.
161: Logically Collective on PC
163: Input Parameters:
164: + pc - the multigrid context
165: . l - the level (0 is coarsest) to supply
166: . residualt - function used to form transpose of residual, if none is provided the previously provide one is used, if no
167: previous one were provided then a default is used
168: - mat - matrix associated with residual
170: Level: advanced
172: .seealso: PCMGResidualTransposeDefault()
173: @*/
174: PetscErrorCode PCMGSetResidualTranspose(PC pc,PetscInt l,PetscErrorCode (*residualt)(Mat,Vec,Vec,Vec),Mat mat)
175: {
176: PC_MG *mg = (PC_MG*)pc->data;
177: PC_MG_Levels **mglevels = mg->levels;
181: if (residualt) mglevels[l]->residualtranspose = residualt;
182: if (!mglevels[l]->residualtranspose) mglevels[l]->residualtranspose = PCMGResidualTransposeDefault;
183: mglevels[l]->matresidualtranspose = PCMGMatResidualTransposeDefault;
184: if (mat) PetscObjectReference((PetscObject)mat);
185: MatDestroy(&mglevels[l]->A);
186: mglevels[l]->A = mat;
187: return 0;
188: }
190: /*@
191: PCMGSetInterpolation - Sets the function to be used to calculate the
192: interpolation from l-1 to the lth level
194: Logically Collective on PC
196: Input Parameters:
197: + pc - the multigrid context
198: . mat - the interpolation operator
199: - l - the level (0 is coarsest) to supply [do not supply 0]
201: Level: advanced
203: Notes:
204: Usually this is the same matrix used also to set the restriction
205: for the same level.
207: One can pass in the interpolation matrix or its transpose; PETSc figures
208: out from the matrix size which one it is.
210: .seealso: PCMGSetRestriction()
211: @*/
212: PetscErrorCode PCMGSetInterpolation(PC pc,PetscInt l,Mat mat)
213: {
214: PC_MG *mg = (PC_MG*)pc->data;
215: PC_MG_Levels **mglevels = mg->levels;
220: PetscObjectReference((PetscObject)mat);
221: MatDestroy(&mglevels[l]->interpolate);
223: mglevels[l]->interpolate = mat;
224: return 0;
225: }
227: /*@
228: PCMGSetOperators - Sets operator and preconditioning matrix for lth level
230: Logically Collective on PC
232: Input Parameters:
233: + pc - the multigrid context
234: . Amat - the operator
235: . pmat - the preconditioning operator
236: - l - the level (0 is the coarsest) to supply
238: Level: advanced
240: .keywords: multigrid, set, interpolate, level
242: .seealso: PCMGSetRestriction(), PCMGSetInterpolation()
243: @*/
244: PetscErrorCode PCMGSetOperators(PC pc,PetscInt l,Mat Amat,Mat Pmat)
245: {
246: PC_MG *mg = (PC_MG*)pc->data;
247: PC_MG_Levels **mglevels = mg->levels;
253: KSPSetOperators(mglevels[l]->smoothd,Amat,Pmat);
254: return 0;
255: }
257: /*@
258: PCMGGetInterpolation - Gets the function to be used to calculate the
259: interpolation from l-1 to the lth level
261: Logically Collective on PC
263: Input Parameters:
264: + pc - the multigrid context
265: - l - the level (0 is coarsest) to supply [Do not supply 0]
267: Output Parameter:
268: . mat - the interpolation matrix, can be NULL
270: Level: advanced
272: .seealso: PCMGGetRestriction(), PCMGSetInterpolation(), PCMGGetRScale()
273: @*/
274: PetscErrorCode PCMGGetInterpolation(PC pc,PetscInt l,Mat *mat)
275: {
276: PC_MG *mg = (PC_MG*)pc->data;
277: PC_MG_Levels **mglevels = mg->levels;
283: if (!mglevels[l]->interpolate) {
285: PCMGSetInterpolation(pc,l,mglevels[l]->restrct);
286: }
287: if (mat) *mat = mglevels[l]->interpolate;
288: return 0;
289: }
291: /*@
292: PCMGSetRestriction - Sets the function to be used to restrict dual vectors
293: from level l to l-1.
295: Logically Collective on PC
297: Input Parameters:
298: + pc - the multigrid context
299: . l - the level (0 is coarsest) to supply [Do not supply 0]
300: - mat - the restriction matrix
302: Level: advanced
304: Notes:
305: Usually this is the same matrix used also to set the interpolation
306: for the same level.
308: One can pass in the interpolation matrix or its transpose; PETSc figures
309: out from the matrix size which one it is.
311: If you do not set this, the transpose of the Mat set with PCMGSetInterpolation()
312: is used.
314: .seealso: PCMGSetInterpolation()
315: @*/
316: PetscErrorCode PCMGSetRestriction(PC pc,PetscInt l,Mat mat)
317: {
318: PC_MG *mg = (PC_MG*)pc->data;
319: PC_MG_Levels **mglevels = mg->levels;
325: PetscObjectReference((PetscObject)mat);
326: MatDestroy(&mglevels[l]->restrct);
328: mglevels[l]->restrct = mat;
329: return 0;
330: }
332: /*@
333: PCMGGetRestriction - Gets the function to be used to restrict dual vectors
334: from level l to l-1.
336: Logically Collective on PC
338: Input Parameters:
339: + pc - the multigrid context
340: - l - the level (0 is coarsest) to supply [Do not supply 0]
342: Output Parameter:
343: . mat - the restriction matrix
345: Level: advanced
347: .seealso: PCMGGetInterpolation(), PCMGSetRestriction(), PCMGGetRScale(), PCMGGetInjection()
348: @*/
349: PetscErrorCode PCMGGetRestriction(PC pc,PetscInt l,Mat *mat)
350: {
351: PC_MG *mg = (PC_MG*)pc->data;
352: PC_MG_Levels **mglevels = mg->levels;
358: if (!mglevels[l]->restrct) {
360: PCMGSetRestriction(pc,l,mglevels[l]->interpolate);
361: }
362: if (mat) *mat = mglevels[l]->restrct;
363: return 0;
364: }
366: /*@
367: PCMGSetRScale - Sets the pointwise scaling for the restriction operator from level l to l-1.
369: Logically Collective on PC
371: Input Parameters:
372: + pc - the multigrid context
373: - l - the level (0 is coarsest) to supply [Do not supply 0]
374: . rscale - the scaling
376: Level: advanced
378: Notes:
379: 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. It is preferable to use PCMGSetInjection() to control moving primal vectors.
381: .seealso: PCMGSetInterpolation(), PCMGSetRestriction(), PCMGGetRScale(), PCMGSetInjection()
382: @*/
383: PetscErrorCode PCMGSetRScale(PC pc,PetscInt l,Vec rscale)
384: {
385: PC_MG *mg = (PC_MG*)pc->data;
386: PC_MG_Levels **mglevels = mg->levels;
391: PetscObjectReference((PetscObject)rscale);
392: VecDestroy(&mglevels[l]->rscale);
394: mglevels[l]->rscale = rscale;
395: return 0;
396: }
398: /*@
399: PCMGGetRScale - Gets the pointwise scaling for the restriction operator from level l to l-1.
401: Collective on PC
403: Input Parameters:
404: + pc - the multigrid context
405: . rscale - the scaling
406: - l - the level (0 is coarsest) to supply [Do not supply 0]
408: Level: advanced
410: Notes:
411: 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. It is preferable to use PCMGGetInjection() to control moving primal vectors.
413: .seealso: PCMGSetInterpolation(), PCMGGetRestriction(), PCMGGetInjection()
414: @*/
415: PetscErrorCode PCMGGetRScale(PC pc,PetscInt l,Vec *rscale)
416: {
417: PC_MG *mg = (PC_MG*)pc->data;
418: PC_MG_Levels **mglevels = mg->levels;
423: if (!mglevels[l]->rscale) {
424: Mat R;
425: Vec X,Y,coarse,fine;
426: PetscInt M,N;
427: PCMGGetRestriction(pc,l,&R);
428: MatCreateVecs(R,&X,&Y);
429: MatGetSize(R,&M,&N);
430: if (M < N) {
431: fine = X;
432: coarse = Y;
433: } else if (N < M) {
434: fine = Y; coarse = X;
435: } else SETERRQ(PetscObjectComm((PetscObject)R),PETSC_ERR_SUP,"Restriction matrix is square, cannot determine which Vec is coarser");
436: VecSet(fine,1.);
437: MatRestrict(R,fine,coarse);
438: VecDestroy(&fine);
439: VecReciprocal(coarse);
440: mglevels[l]->rscale = coarse;
441: }
442: *rscale = mglevels[l]->rscale;
443: return 0;
444: }
446: /*@
447: PCMGSetInjection - Sets the function to be used to inject primal vectors
448: from level l to l-1.
450: Logically Collective on PC
452: Input Parameters:
453: + pc - the multigrid context
454: . l - the level (0 is coarsest) to supply [Do not supply 0]
455: - mat - the injection matrix
457: Level: advanced
459: .seealso: PCMGSetRestriction()
460: @*/
461: PetscErrorCode PCMGSetInjection(PC pc,PetscInt l,Mat mat)
462: {
463: PC_MG *mg = (PC_MG*)pc->data;
464: PC_MG_Levels **mglevels = mg->levels;
470: PetscObjectReference((PetscObject)mat);
471: MatDestroy(&mglevels[l]->inject);
473: mglevels[l]->inject = mat;
474: return 0;
475: }
477: /*@
478: PCMGGetInjection - Gets the function to be used to inject primal vectors
479: from level l to l-1.
481: Logically Collective on PC
483: Input Parameters:
484: + pc - the multigrid context
485: - l - the level (0 is coarsest) to supply [Do not supply 0]
487: Output Parameter:
488: . mat - the restriction matrix (may be NULL if no injection is available).
490: Level: advanced
492: .seealso: PCMGSetInjection(), PCMGetGetRestriction()
493: @*/
494: PetscErrorCode PCMGGetInjection(PC pc,PetscInt l,Mat *mat)
495: {
496: PC_MG *mg = (PC_MG*)pc->data;
497: PC_MG_Levels **mglevels = mg->levels;
503: if (mat) *mat = mglevels[l]->inject;
504: return 0;
505: }
507: /*@
508: PCMGGetSmoother - Gets the KSP context to be used as smoother for
509: both pre- and post-smoothing. Call both PCMGGetSmootherUp() and
510: PCMGGetSmootherDown() to use different functions for pre- and
511: post-smoothing.
513: Not Collective, KSP returned is parallel if PC is
515: Input Parameters:
516: + pc - the multigrid context
517: - l - the level (0 is coarsest) to supply
519: Output Parameter:
520: . ksp - the smoother
522: Notes:
523: 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.
524: You can also modify smoother options by calling the various KSPSetXXX() options on this ksp. In addition you can call KSPGetPC(ksp,&pc)
525: and modify PC options for the smoother; for example PCSetType(pc,PCSOR); to use SOR smoothing.
527: Level: advanced
529: .seealso: PCMGGetSmootherUp(), PCMGGetSmootherDown(), PCMGGetCoarseSolve()
530: @*/
531: PetscErrorCode PCMGGetSmoother(PC pc,PetscInt l,KSP *ksp)
532: {
533: PC_MG *mg = (PC_MG*)pc->data;
534: PC_MG_Levels **mglevels = mg->levels;
537: *ksp = mglevels[l]->smoothd;
538: return 0;
539: }
541: /*@
542: PCMGGetSmootherUp - Gets the KSP context to be used as smoother after
543: coarse grid correction (post-smoother).
545: Not Collective, KSP returned is parallel if PC is
547: Input Parameters:
548: + pc - the multigrid context
549: - l - the level (0 is coarsest) to supply
551: Output Parameter:
552: . ksp - the smoother
554: Level: advanced
556: Notes:
557: calling this will result in a different pre and post smoother so you may need to
558: set options on the pre smoother also
560: .seealso: PCMGGetSmootherUp(), PCMGGetSmootherDown()
561: @*/
562: PetscErrorCode PCMGGetSmootherUp(PC pc,PetscInt l,KSP *ksp)
563: {
564: PC_MG *mg = (PC_MG*)pc->data;
565: PC_MG_Levels **mglevels = mg->levels;
566: const char *prefix;
567: MPI_Comm comm;
570: /*
571: This is called only if user wants a different pre-smoother from post.
572: Thus we check if a different one has already been allocated,
573: if not we allocate it.
574: */
576: if (mglevels[l]->smoothu == mglevels[l]->smoothd) {
577: KSPType ksptype;
578: PCType pctype;
579: PC ipc;
580: PetscReal rtol,abstol,dtol;
581: PetscInt maxits;
582: KSPNormType normtype;
583: PetscObjectGetComm((PetscObject)mglevels[l]->smoothd,&comm);
584: KSPGetOptionsPrefix(mglevels[l]->smoothd,&prefix);
585: KSPGetTolerances(mglevels[l]->smoothd,&rtol,&abstol,&dtol,&maxits);
586: KSPGetType(mglevels[l]->smoothd,&ksptype);
587: KSPGetNormType(mglevels[l]->smoothd,&normtype);
588: KSPGetPC(mglevels[l]->smoothd,&ipc);
589: PCGetType(ipc,&pctype);
591: KSPCreate(comm,&mglevels[l]->smoothu);
592: KSPSetErrorIfNotConverged(mglevels[l]->smoothu,pc->erroriffailure);
593: PetscObjectIncrementTabLevel((PetscObject)mglevels[l]->smoothu,(PetscObject)pc,mglevels[0]->levels-l);
594: KSPSetOptionsPrefix(mglevels[l]->smoothu,prefix);
595: KSPSetTolerances(mglevels[l]->smoothu,rtol,abstol,dtol,maxits);
596: KSPSetType(mglevels[l]->smoothu,ksptype);
597: KSPSetNormType(mglevels[l]->smoothu,normtype);
598: KSPSetConvergenceTest(mglevels[l]->smoothu,KSPConvergedSkip,NULL,NULL);
599: KSPGetPC(mglevels[l]->smoothu,&ipc);
600: PCSetType(ipc,pctype);
601: PetscLogObjectParent((PetscObject)pc,(PetscObject)mglevels[l]->smoothu);
602: PetscObjectComposedDataSetInt((PetscObject) mglevels[l]->smoothu, PetscMGLevelId, mglevels[l]->level);
603: }
604: if (ksp) *ksp = mglevels[l]->smoothu;
605: return 0;
606: }
608: /*@
609: PCMGGetSmootherDown - Gets the KSP context to be used as smoother before
610: coarse grid correction (pre-smoother).
612: Not Collective, KSP returned is parallel if PC is
614: Input Parameters:
615: + pc - the multigrid context
616: - l - the level (0 is coarsest) to supply
618: Output Parameter:
619: . ksp - the smoother
621: Level: advanced
623: Notes:
624: calling this will result in a different pre and post smoother so you may need to
625: set options on the post smoother also
627: .seealso: PCMGGetSmootherUp(), PCMGGetSmoother()
628: @*/
629: PetscErrorCode PCMGGetSmootherDown(PC pc,PetscInt l,KSP *ksp)
630: {
631: PC_MG *mg = (PC_MG*)pc->data;
632: PC_MG_Levels **mglevels = mg->levels;
635: /* make sure smoother up and down are different */
636: if (l) {
637: PCMGGetSmootherUp(pc,l,NULL);
638: }
639: *ksp = mglevels[l]->smoothd;
640: return 0;
641: }
643: /*@
644: PCMGSetCycleTypeOnLevel - Sets the type of cycle (aka cycle index) to run on the specified level.
646: Logically Collective on PC
648: Input Parameters:
649: + pc - the multigrid context
650: . l - the level (0 is coarsest)
651: - c - either PC_MG_CYCLE_V or PC_MG_CYCLE_W
653: Level: advanced
655: .seealso: PCMGSetCycleType()
656: @*/
657: PetscErrorCode PCMGSetCycleTypeOnLevel(PC pc,PetscInt l,PCMGCycleType c)
658: {
659: PC_MG *mg = (PC_MG*)pc->data;
660: PC_MG_Levels **mglevels = mg->levels;
666: mglevels[l]->cycles = c;
667: return 0;
668: }
670: /*@
671: PCMGSetRhs - Sets the vector to be used to store the right-hand side on a particular level.
673: Logically Collective on PC
675: Input Parameters:
676: + pc - the multigrid context
677: . l - the level (0 is coarsest) this is to be used for
678: - c - the Vec
680: Level: advanced
682: Notes:
683: If this is not provided PETSc will automatically generate one. You do not need to keep a reference to this vector if you do not need it. PCDestroy() will properly free it.
685: .keywords: MG, multigrid, set, right-hand-side, rhs, level
686: .seealso: PCMGSetX(), PCMGSetR()
687: @*/
688: PetscErrorCode PCMGSetRhs(PC pc,PetscInt l,Vec c)
689: {
690: PC_MG *mg = (PC_MG*)pc->data;
691: PC_MG_Levels **mglevels = mg->levels;
696: PetscObjectReference((PetscObject)c);
697: VecDestroy(&mglevels[l]->b);
699: mglevels[l]->b = c;
700: return 0;
701: }
703: /*@
704: PCMGSetX - Sets the vector to be used to store the solution on a particular level.
706: Logically Collective on PC
708: Input Parameters:
709: + pc - the multigrid context
710: . l - the level (0 is coarsest) this is to be used for (do not supply the finest level)
711: - c - the Vec
713: Level: advanced
715: Notes:
716: If this is not provided PETSc will automatically generate one. You do not need to keep a reference to this vector if you do not need it. PCDestroy() will properly free it.
718: .keywords: MG, multigrid, set, solution, level
719: .seealso: PCMGSetRhs(), PCMGSetR()
720: @*/
721: PetscErrorCode PCMGSetX(PC pc,PetscInt l,Vec c)
722: {
723: PC_MG *mg = (PC_MG*)pc->data;
724: PC_MG_Levels **mglevels = mg->levels;
729: PetscObjectReference((PetscObject)c);
730: VecDestroy(&mglevels[l]->x);
732: mglevels[l]->x = c;
733: return 0;
734: }
736: /*@
737: PCMGSetR - Sets the vector to be used to store the residual on a particular level.
739: Logically Collective on PC
741: Input Parameters:
742: + pc - the multigrid context
743: . l - the level (0 is coarsest) this is to be used for
744: - c - the Vec
746: Level: advanced
748: Notes:
749: If this is not provided PETSc will automatically generate one. You do not need to keep a reference to this vector if you do not need it. PCDestroy() will properly free it.
751: .keywords: MG, multigrid, set, residual, level
752: .seealso: PCMGSetRhs(), PCMGSetX()
753: @*/
754: PetscErrorCode PCMGSetR(PC pc,PetscInt l,Vec c)
755: {
756: PC_MG *mg = (PC_MG*)pc->data;
757: PC_MG_Levels **mglevels = mg->levels;
762: PetscObjectReference((PetscObject)c);
763: VecDestroy(&mglevels[l]->r);
765: mglevels[l]->r = c;
766: return 0;
767: }