Actual source code: mg.c
2: /*
3: Defines the multigrid preconditioner interface.
4: */
5: #include <petsc/private/pcmgimpl.h>
6: #include <petsc/private/kspimpl.h>
7: #include <petscdm.h>
8: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC, PetscBool *);
10: /*
11: Contains the list of registered coarse space construction routines
12: */
13: PetscFunctionList PCMGCoarseList = NULL;
15: PetscErrorCode PCMGMCycle_Private(PC pc, PC_MG_Levels **mglevelsin, PetscBool transpose, PetscBool matapp, PCRichardsonConvergedReason *reason)
16: {
17: PC_MG *mg = (PC_MG *)pc->data;
18: PC_MG_Levels *mgc, *mglevels = *mglevelsin;
19: PetscInt cycles = (mglevels->level == 1) ? 1 : (PetscInt)mglevels->cycles;
21: if (mglevels->eventsmoothsolve) PetscLogEventBegin(mglevels->eventsmoothsolve, 0, 0, 0, 0);
22: if (!transpose) {
23: if (matapp) {
24: KSPMatSolve(mglevels->smoothd, mglevels->B, mglevels->X); /* pre-smooth */
25: KSPCheckSolve(mglevels->smoothd, pc, NULL);
26: } else {
27: KSPSolve(mglevels->smoothd, mglevels->b, mglevels->x); /* pre-smooth */
28: KSPCheckSolve(mglevels->smoothd, pc, mglevels->x);
29: }
30: } else {
32: KSPSolveTranspose(mglevels->smoothu, mglevels->b, mglevels->x); /* transpose of post-smooth */
33: KSPCheckSolve(mglevels->smoothu, pc, mglevels->x);
34: }
35: if (mglevels->eventsmoothsolve) PetscLogEventEnd(mglevels->eventsmoothsolve, 0, 0, 0, 0);
36: if (mglevels->level) { /* not the coarsest grid */
37: if (mglevels->eventresidual) PetscLogEventBegin(mglevels->eventresidual, 0, 0, 0, 0);
38: if (matapp && !mglevels->R) MatDuplicate(mglevels->B, MAT_DO_NOT_COPY_VALUES, &mglevels->R);
39: if (!transpose) {
40: if (matapp) (*mglevels->matresidual)(mglevels->A, mglevels->B, mglevels->X, mglevels->R);
41: else (*mglevels->residual)(mglevels->A, mglevels->b, mglevels->x, mglevels->r);
42: } else {
43: if (matapp) (*mglevels->matresidualtranspose)(mglevels->A, mglevels->B, mglevels->X, mglevels->R);
44: else (*mglevels->residualtranspose)(mglevels->A, mglevels->b, mglevels->x, mglevels->r);
45: }
46: if (mglevels->eventresidual) PetscLogEventEnd(mglevels->eventresidual, 0, 0, 0, 0);
48: /* if on finest level and have convergence criteria set */
49: if (mglevels->level == mglevels->levels - 1 && mg->ttol && reason) {
50: PetscReal rnorm;
51: VecNorm(mglevels->r, NORM_2, &rnorm);
52: if (rnorm <= mg->ttol) {
53: if (rnorm < mg->abstol) {
54: *reason = PCRICHARDSON_CONVERGED_ATOL;
55: PetscInfo(pc, "Linear solver has converged. Residual norm %g is less than absolute tolerance %g\n", (double)rnorm, (double)mg->abstol);
56: } else {
57: *reason = PCRICHARDSON_CONVERGED_RTOL;
58: PetscInfo(pc, "Linear solver has converged. Residual norm %g is less than relative tolerance times initial residual norm %g\n", (double)rnorm, (double)mg->ttol);
59: }
60: return 0;
61: }
62: }
64: mgc = *(mglevelsin - 1);
65: if (mglevels->eventinterprestrict) PetscLogEventBegin(mglevels->eventinterprestrict, 0, 0, 0, 0);
66: if (!transpose) {
67: if (matapp) MatMatRestrict(mglevels->restrct, mglevels->R, &mgc->B);
68: else MatRestrict(mglevels->restrct, mglevels->r, mgc->b);
69: } else {
70: if (matapp) MatMatRestrict(mglevels->interpolate, mglevels->R, &mgc->B);
71: else MatRestrict(mglevels->interpolate, mglevels->r, mgc->b);
72: }
73: if (mglevels->eventinterprestrict) PetscLogEventEnd(mglevels->eventinterprestrict, 0, 0, 0, 0);
74: if (matapp) {
75: if (!mgc->X) {
76: MatDuplicate(mgc->B, MAT_DO_NOT_COPY_VALUES, &mgc->X);
77: } else {
78: MatZeroEntries(mgc->X);
79: }
80: } else {
81: VecZeroEntries(mgc->x);
82: }
83: while (cycles--) PCMGMCycle_Private(pc, mglevelsin - 1, transpose, matapp, reason);
84: if (mglevels->eventinterprestrict) PetscLogEventBegin(mglevels->eventinterprestrict, 0, 0, 0, 0);
85: if (!transpose) {
86: if (matapp) MatMatInterpolateAdd(mglevels->interpolate, mgc->X, mglevels->X, &mglevels->X);
87: else MatInterpolateAdd(mglevels->interpolate, mgc->x, mglevels->x, mglevels->x);
88: } else {
89: MatInterpolateAdd(mglevels->restrct, mgc->x, mglevels->x, mglevels->x);
90: }
91: if (mglevels->eventinterprestrict) PetscLogEventEnd(mglevels->eventinterprestrict, 0, 0, 0, 0);
92: if (mglevels->eventsmoothsolve) PetscLogEventBegin(mglevels->eventsmoothsolve, 0, 0, 0, 0);
93: if (!transpose) {
94: if (matapp) {
95: KSPMatSolve(mglevels->smoothu, mglevels->B, mglevels->X); /* post smooth */
96: KSPCheckSolve(mglevels->smoothu, pc, NULL);
97: } else {
98: KSPSolve(mglevels->smoothu, mglevels->b, mglevels->x); /* post smooth */
99: KSPCheckSolve(mglevels->smoothu, pc, mglevels->x);
100: }
101: } else {
103: KSPSolveTranspose(mglevels->smoothd, mglevels->b, mglevels->x); /* post smooth */
104: KSPCheckSolve(mglevels->smoothd, pc, mglevels->x);
105: }
106: if (mglevels->cr) {
108: /* TODO Turn on copy and turn off noisy if we have an exact solution
109: VecCopy(mglevels->x, mglevels->crx);
110: VecCopy(mglevels->b, mglevels->crb); */
111: KSPSetNoisy_Private(mglevels->crx);
112: KSPSolve(mglevels->cr, mglevels->crb, mglevels->crx); /* compatible relaxation */
113: KSPCheckSolve(mglevels->cr, pc, mglevels->crx);
114: }
115: if (mglevels->eventsmoothsolve) PetscLogEventEnd(mglevels->eventsmoothsolve, 0, 0, 0, 0);
116: }
117: return 0;
118: }
120: static PetscErrorCode PCApplyRichardson_MG(PC pc, Vec b, Vec x, Vec w, PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt its, PetscBool zeroguess, PetscInt *outits, PCRichardsonConvergedReason *reason)
121: {
122: PC_MG *mg = (PC_MG *)pc->data;
123: PC_MG_Levels **mglevels = mg->levels;
124: PC tpc;
125: PetscBool changeu, changed;
126: PetscInt levels = mglevels[0]->levels, i;
128: /* When the DM is supplying the matrix then it will not exist until here */
129: for (i = 0; i < levels; i++) {
130: if (!mglevels[i]->A) {
131: KSPGetOperators(mglevels[i]->smoothu, &mglevels[i]->A, NULL);
132: PetscObjectReference((PetscObject)mglevels[i]->A);
133: }
134: }
136: KSPGetPC(mglevels[levels - 1]->smoothd, &tpc);
137: PCPreSolveChangeRHS(tpc, &changed);
138: KSPGetPC(mglevels[levels - 1]->smoothu, &tpc);
139: PCPreSolveChangeRHS(tpc, &changeu);
140: if (!changed && !changeu) {
141: VecDestroy(&mglevels[levels - 1]->b);
142: mglevels[levels - 1]->b = b;
143: } else { /* if the smoother changes the rhs during PreSolve, we cannot use the input vector */
144: if (!mglevels[levels - 1]->b) {
145: Vec *vec;
147: KSPCreateVecs(mglevels[levels - 1]->smoothd, 1, &vec, 0, NULL);
148: mglevels[levels - 1]->b = *vec;
149: PetscFree(vec);
150: }
151: VecCopy(b, mglevels[levels - 1]->b);
152: }
153: mglevels[levels - 1]->x = x;
155: mg->rtol = rtol;
156: mg->abstol = abstol;
157: mg->dtol = dtol;
158: if (rtol) {
159: /* compute initial residual norm for relative convergence test */
160: PetscReal rnorm;
161: if (zeroguess) {
162: VecNorm(b, NORM_2, &rnorm);
163: } else {
164: (*mglevels[levels - 1]->residual)(mglevels[levels - 1]->A, b, x, w);
165: VecNorm(w, NORM_2, &rnorm);
166: }
167: mg->ttol = PetscMax(rtol * rnorm, abstol);
168: } else if (abstol) mg->ttol = abstol;
169: else mg->ttol = 0.0;
171: /* since smoother is applied to full system, not just residual we need to make sure that smoothers don't
172: stop prematurely due to small residual */
173: for (i = 1; i < levels; i++) {
174: KSPSetTolerances(mglevels[i]->smoothu, 0, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
175: if (mglevels[i]->smoothu != mglevels[i]->smoothd) {
176: /* For Richardson the initial guess is nonzero since it is solving in each cycle the original system not just applying as a preconditioner */
177: KSPSetInitialGuessNonzero(mglevels[i]->smoothd, PETSC_TRUE);
178: KSPSetTolerances(mglevels[i]->smoothd, 0, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
179: }
180: }
182: *reason = (PCRichardsonConvergedReason)0;
183: for (i = 0; i < its; i++) {
184: PCMGMCycle_Private(pc, mglevels + levels - 1, PETSC_FALSE, PETSC_FALSE, reason);
185: if (*reason) break;
186: }
187: if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS;
188: *outits = i;
189: if (!changed && !changeu) mglevels[levels - 1]->b = NULL;
190: return 0;
191: }
193: PetscErrorCode PCReset_MG(PC pc)
194: {
195: PC_MG *mg = (PC_MG *)pc->data;
196: PC_MG_Levels **mglevels = mg->levels;
197: PetscInt i, n;
199: if (mglevels) {
200: n = mglevels[0]->levels;
201: for (i = 0; i < n - 1; i++) {
202: VecDestroy(&mglevels[i + 1]->r);
203: VecDestroy(&mglevels[i]->b);
204: VecDestroy(&mglevels[i]->x);
205: MatDestroy(&mglevels[i + 1]->R);
206: MatDestroy(&mglevels[i]->B);
207: MatDestroy(&mglevels[i]->X);
208: VecDestroy(&mglevels[i]->crx);
209: VecDestroy(&mglevels[i]->crb);
210: MatDestroy(&mglevels[i + 1]->restrct);
211: MatDestroy(&mglevels[i + 1]->interpolate);
212: MatDestroy(&mglevels[i + 1]->inject);
213: VecDestroy(&mglevels[i + 1]->rscale);
214: }
215: VecDestroy(&mglevels[n - 1]->crx);
216: VecDestroy(&mglevels[n - 1]->crb);
217: /* this is not null only if the smoother on the finest level
218: changes the rhs during PreSolve */
219: VecDestroy(&mglevels[n - 1]->b);
220: MatDestroy(&mglevels[n - 1]->B);
222: for (i = 0; i < n; i++) {
223: MatDestroy(&mglevels[i]->coarseSpace);
224: MatDestroy(&mglevels[i]->A);
225: if (mglevels[i]->smoothd != mglevels[i]->smoothu) KSPReset(mglevels[i]->smoothd);
226: KSPReset(mglevels[i]->smoothu);
227: if (mglevels[i]->cr) KSPReset(mglevels[i]->cr);
228: }
229: mg->Nc = 0;
230: }
231: return 0;
232: }
234: /* Implementing CR
236: We only want to make corrections that ``do not change'' the coarse solution. What we mean by not changing is that if I prolong my coarse solution to the fine grid and then inject that fine solution back to the coarse grid, I get the same answer. Injection is what Brannick calls R. We want the complementary projector to Inj, which we will call S, after Brannick, so that Inj S = 0. Now the orthogonal projector onto the range of Inj^T is
238: Inj^T (Inj Inj^T)^{-1} Inj
240: and if Inj is a VecScatter, as it is now in PETSc, we have
242: Inj^T Inj
244: and
246: S = I - Inj^T Inj
248: since
250: Inj S = Inj - (Inj Inj^T) Inj = 0.
252: Brannick suggests
254: A \to S^T A S \qquad\mathrm{and}\qquad M \to S^T M S
256: but I do not think his :math:`S^T S = I` is correct. Our S is an orthogonal projector, so :math:`S^T S = S^2 = S`. We will use
258: M^{-1} A \to S M^{-1} A S
260: In fact, since it is somewhat hard in PETSc to do the symmetric application, we will just apply S on the left.
262: Check: || Inj P - I ||_F < tol
263: Check: In general, Inj Inj^T = I
264: */
266: typedef struct {
267: PC mg; /* The PCMG object */
268: PetscInt l; /* The multigrid level for this solver */
269: Mat Inj; /* The injection matrix */
270: Mat S; /* I - Inj^T Inj */
271: } CRContext;
273: static PetscErrorCode CRSetup_Private(PC pc)
274: {
275: CRContext *ctx;
276: Mat It;
279: PCShellGetContext(pc, &ctx);
280: PCMGGetInjection(ctx->mg, ctx->l, &It);
282: MatCreateTranspose(It, &ctx->Inj);
283: MatCreateNormal(ctx->Inj, &ctx->S);
284: MatScale(ctx->S, -1.0);
285: MatShift(ctx->S, 1.0);
286: return 0;
287: }
289: static PetscErrorCode CRApply_Private(PC pc, Vec x, Vec y)
290: {
291: CRContext *ctx;
294: PCShellGetContext(pc, &ctx);
295: MatMult(ctx->S, x, y);
296: return 0;
297: }
299: static PetscErrorCode CRDestroy_Private(PC pc)
300: {
301: CRContext *ctx;
304: PCShellGetContext(pc, &ctx);
305: MatDestroy(&ctx->Inj);
306: MatDestroy(&ctx->S);
307: PetscFree(ctx);
308: PCShellSetContext(pc, NULL);
309: return 0;
310: }
312: static PetscErrorCode CreateCR_Private(PC pc, PetscInt l, PC *cr)
313: {
314: CRContext *ctx;
317: PCCreate(PetscObjectComm((PetscObject)pc), cr);
318: PetscObjectSetName((PetscObject)*cr, "S (complementary projector to injection)");
319: PetscCalloc1(1, &ctx);
320: ctx->mg = pc;
321: ctx->l = l;
322: PCSetType(*cr, PCSHELL);
323: PCShellSetContext(*cr, ctx);
324: PCShellSetApply(*cr, CRApply_Private);
325: PCShellSetSetUp(*cr, CRSetup_Private);
326: PCShellSetDestroy(*cr, CRDestroy_Private);
327: return 0;
328: }
330: PetscErrorCode PCMGSetLevels_MG(PC pc, PetscInt levels, MPI_Comm *comms)
331: {
332: PC_MG *mg = (PC_MG *)pc->data;
333: MPI_Comm comm;
334: PC_MG_Levels **mglevels = mg->levels;
335: PCMGType mgtype = mg->am;
336: PetscInt mgctype = (PetscInt)PC_MG_CYCLE_V;
337: PetscInt i;
338: PetscMPIInt size;
339: const char *prefix;
340: PC ipc;
341: PetscInt n;
345: if (mg->nlevels == levels) return 0;
346: PetscObjectGetComm((PetscObject)pc, &comm);
347: if (mglevels) {
348: mgctype = mglevels[0]->cycles;
349: /* changing the number of levels so free up the previous stuff */
350: PCReset_MG(pc);
351: n = mglevels[0]->levels;
352: for (i = 0; i < n; i++) {
353: if (mglevels[i]->smoothd != mglevels[i]->smoothu) KSPDestroy(&mglevels[i]->smoothd);
354: KSPDestroy(&mglevels[i]->smoothu);
355: KSPDestroy(&mglevels[i]->cr);
356: PetscFree(mglevels[i]);
357: }
358: PetscFree(mg->levels);
359: }
361: mg->nlevels = levels;
363: PetscMalloc1(levels, &mglevels);
365: PCGetOptionsPrefix(pc, &prefix);
367: mg->stageApply = 0;
368: for (i = 0; i < levels; i++) {
369: PetscNew(&mglevels[i]);
371: mglevels[i]->level = i;
372: mglevels[i]->levels = levels;
373: mglevels[i]->cycles = mgctype;
374: mg->default_smoothu = 2;
375: mg->default_smoothd = 2;
376: mglevels[i]->eventsmoothsetup = 0;
377: mglevels[i]->eventsmoothsolve = 0;
378: mglevels[i]->eventresidual = 0;
379: mglevels[i]->eventinterprestrict = 0;
381: if (comms) comm = comms[i];
382: if (comm != MPI_COMM_NULL) {
383: KSPCreate(comm, &mglevels[i]->smoothd);
384: KSPSetErrorIfNotConverged(mglevels[i]->smoothd, pc->erroriffailure);
385: PetscObjectIncrementTabLevel((PetscObject)mglevels[i]->smoothd, (PetscObject)pc, levels - i);
386: KSPSetOptionsPrefix(mglevels[i]->smoothd, prefix);
387: PetscObjectComposedDataSetInt((PetscObject)mglevels[i]->smoothd, PetscMGLevelId, mglevels[i]->level);
388: if (i || levels == 1) {
389: char tprefix[128];
391: KSPSetType(mglevels[i]->smoothd, KSPCHEBYSHEV);
392: KSPSetConvergenceTest(mglevels[i]->smoothd, KSPConvergedSkip, NULL, NULL);
393: KSPSetNormType(mglevels[i]->smoothd, KSP_NORM_NONE);
394: KSPGetPC(mglevels[i]->smoothd, &ipc);
395: PCSetType(ipc, PCSOR);
396: KSPSetTolerances(mglevels[i]->smoothd, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, mg->default_smoothd);
398: PetscSNPrintf(tprefix, 128, "mg_levels_%d_", (int)i);
399: KSPAppendOptionsPrefix(mglevels[i]->smoothd, tprefix);
400: } else {
401: KSPAppendOptionsPrefix(mglevels[0]->smoothd, "mg_coarse_");
403: /* coarse solve is (redundant) LU by default; set shifttype NONZERO to avoid annoying zero-pivot in LU preconditioner */
404: KSPSetType(mglevels[0]->smoothd, KSPPREONLY);
405: KSPGetPC(mglevels[0]->smoothd, &ipc);
406: MPI_Comm_size(comm, &size);
407: if (size > 1) {
408: PCSetType(ipc, PCREDUNDANT);
409: } else {
410: PCSetType(ipc, PCLU);
411: }
412: PCFactorSetShiftType(ipc, MAT_SHIFT_INBLOCKS);
413: }
414: }
415: mglevels[i]->smoothu = mglevels[i]->smoothd;
416: mg->rtol = 0.0;
417: mg->abstol = 0.0;
418: mg->dtol = 0.0;
419: mg->ttol = 0.0;
420: mg->cyclesperpcapply = 1;
421: }
422: mg->levels = mglevels;
423: PCMGSetType(pc, mgtype);
424: return 0;
425: }
427: /*@C
428: PCMGSetLevels - Sets the number of levels to use with `PCMG`.
429: Must be called before any other `PCMG` routine.
431: Logically Collective
433: Input Parameters:
434: + pc - the preconditioner context
435: . levels - the number of levels
436: - comms - optional communicators for each level; this is to allow solving the coarser problems
437: on smaller sets of processes. For processes that are not included in the computation
438: you must pass `MPI_COMM_NULL`. Use comms = NULL to specify that all processes
439: should participate in each level of problem.
441: Level: intermediate
443: Notes:
444: If the number of levels is one then the multigrid uses the -mg_levels prefix
445: for setting the level options rather than the -mg_coarse prefix.
447: You can free the information in comms after this routine is called.
449: The array of MPI communicators must contain `MPI_COMM_NULL` for those ranks that at each level
450: are not participating in the coarser solve. For example, with 2 levels and 1 and 2 ranks on
451: the two levels, rank 0 in the original communicator will pass in an array of 2 communicators
452: of size 2 and 1, while rank 1 in the original communicator will pass in array of 2 communicators
453: the first of size 2 and the second of value `MPI_COMM_NULL` since the rank 1 does not participate
454: in the coarse grid solve.
456: Since each coarser level may have a new `MPI_Comm` with fewer ranks than the previous, one
457: must take special care in providing the restriction and interpolation operation. We recommend
458: providing these as two step operations; first perform a standard restriction or interpolation on
459: the full number of ranks for that level and then use an MPI call to copy the resulting vector
460: array entries (after calls to VecGetArray()) to the smaller or larger number of ranks, note in both
461: cases the MPI calls must be made on the larger of the two communicators. Traditional MPI send and
462: receives or MPI_AlltoAllv() could be used to do the reshuffling of the vector entries.
464: Fortran Note:
465: Use comms = `PETSC_NULL_MPI_COMM` as the equivalent of NULL in the C interface. Note `PETSC_NULL_MPI_COMM`
466: is not `MPI_COMM_NULL`. It is more like `PETSC_NULL_INTEGER`, `PETSC_NULL_REAL` etc.
468: .seealso: `PCMGSetType()`, `PCMGGetLevels()`
469: @*/
470: PetscErrorCode PCMGSetLevels(PC pc, PetscInt levels, MPI_Comm *comms)
471: {
474: PetscTryMethod(pc, "PCMGSetLevels_C", (PC, PetscInt, MPI_Comm *), (pc, levels, comms));
475: return 0;
476: }
478: PetscErrorCode PCDestroy_MG(PC pc)
479: {
480: PC_MG *mg = (PC_MG *)pc->data;
481: PC_MG_Levels **mglevels = mg->levels;
482: PetscInt i, n;
484: PCReset_MG(pc);
485: if (mglevels) {
486: n = mglevels[0]->levels;
487: for (i = 0; i < n; i++) {
488: if (mglevels[i]->smoothd != mglevels[i]->smoothu) KSPDestroy(&mglevels[i]->smoothd);
489: KSPDestroy(&mglevels[i]->smoothu);
490: KSPDestroy(&mglevels[i]->cr);
491: PetscFree(mglevels[i]);
492: }
493: PetscFree(mg->levels);
494: }
495: PetscFree(pc->data);
496: PetscObjectComposeFunction((PetscObject)pc, "PCGetInterpolations_C", NULL);
497: PetscObjectComposeFunction((PetscObject)pc, "PCGetCoarseOperators_C", NULL);
498: PetscObjectComposeFunction((PetscObject)pc, "PCMGSetGalerkin_C", NULL);
499: PetscObjectComposeFunction((PetscObject)pc, "PCMGGetLevels_C", NULL);
500: PetscObjectComposeFunction((PetscObject)pc, "PCMGSetLevels_C", NULL);
501: PetscObjectComposeFunction((PetscObject)pc, "PCGetInterpolations_C", NULL);
502: PetscObjectComposeFunction((PetscObject)pc, "PCGetCoarseOperators_C", NULL);
503: PetscObjectComposeFunction((PetscObject)pc, "PCMGSetAdaptInterpolation_C", NULL);
504: PetscObjectComposeFunction((PetscObject)pc, "PCMGGetAdaptInterpolation_C", NULL);
505: PetscObjectComposeFunction((PetscObject)pc, "PCMGSetAdaptCR_C", NULL);
506: PetscObjectComposeFunction((PetscObject)pc, "PCMGGetAdaptCR_C", NULL);
507: PetscObjectComposeFunction((PetscObject)pc, "PCMGSetAdaptCoarseSpaceType_C", NULL);
508: PetscObjectComposeFunction((PetscObject)pc, "PCMGGetAdaptCoarseSpaceType_C", NULL);
509: return 0;
510: }
512: /*
513: PCApply_MG - Runs either an additive, multiplicative, Kaskadic
514: or full cycle of multigrid.
516: Note:
517: A simple wrapper which calls PCMGMCycle(),PCMGACycle(), or PCMGFCycle().
518: */
519: static PetscErrorCode PCApply_MG_Internal(PC pc, Vec b, Vec x, Mat B, Mat X, PetscBool transpose)
520: {
521: PC_MG *mg = (PC_MG *)pc->data;
522: PC_MG_Levels **mglevels = mg->levels;
523: PC tpc;
524: PetscInt levels = mglevels[0]->levels, i;
525: PetscBool changeu, changed, matapp;
527: matapp = (PetscBool)(B && X);
528: if (mg->stageApply) PetscLogStagePush(mg->stageApply);
529: /* When the DM is supplying the matrix then it will not exist until here */
530: for (i = 0; i < levels; i++) {
531: if (!mglevels[i]->A) {
532: KSPGetOperators(mglevels[i]->smoothu, &mglevels[i]->A, NULL);
533: PetscObjectReference((PetscObject)mglevels[i]->A);
534: }
535: }
537: KSPGetPC(mglevels[levels - 1]->smoothd, &tpc);
538: PCPreSolveChangeRHS(tpc, &changed);
539: KSPGetPC(mglevels[levels - 1]->smoothu, &tpc);
540: PCPreSolveChangeRHS(tpc, &changeu);
541: if (!changeu && !changed) {
542: if (matapp) {
543: MatDestroy(&mglevels[levels - 1]->B);
544: mglevels[levels - 1]->B = B;
545: } else {
546: VecDestroy(&mglevels[levels - 1]->b);
547: mglevels[levels - 1]->b = b;
548: }
549: } else { /* if the smoother changes the rhs during PreSolve, we cannot use the input vector */
550: if (matapp) {
551: if (mglevels[levels - 1]->B) {
552: PetscInt N1, N2;
553: PetscBool flg;
555: MatGetSize(mglevels[levels - 1]->B, NULL, &N1);
556: MatGetSize(B, NULL, &N2);
557: PetscObjectTypeCompare((PetscObject)mglevels[levels - 1]->B, ((PetscObject)B)->type_name, &flg);
558: if (N1 != N2 || !flg) MatDestroy(&mglevels[levels - 1]->B);
559: }
560: if (!mglevels[levels - 1]->B) {
561: MatDuplicate(B, MAT_COPY_VALUES, &mglevels[levels - 1]->B);
562: } else {
563: MatCopy(B, mglevels[levels - 1]->B, SAME_NONZERO_PATTERN);
564: }
565: } else {
566: if (!mglevels[levels - 1]->b) {
567: Vec *vec;
569: KSPCreateVecs(mglevels[levels - 1]->smoothd, 1, &vec, 0, NULL);
570: mglevels[levels - 1]->b = *vec;
571: PetscFree(vec);
572: }
573: VecCopy(b, mglevels[levels - 1]->b);
574: }
575: }
576: if (matapp) {
577: mglevels[levels - 1]->X = X;
578: } else {
579: mglevels[levels - 1]->x = x;
580: }
582: /* If coarser Xs are present, it means we have already block applied the PC at least once
583: Reset operators if sizes/type do no match */
584: if (matapp && levels > 1 && mglevels[levels - 2]->X) {
585: PetscInt Xc, Bc;
586: PetscBool flg;
588: MatGetSize(mglevels[levels - 2]->X, NULL, &Xc);
589: MatGetSize(mglevels[levels - 1]->B, NULL, &Bc);
590: PetscObjectTypeCompare((PetscObject)mglevels[levels - 2]->X, ((PetscObject)mglevels[levels - 1]->X)->type_name, &flg);
591: if (Xc != Bc || !flg) {
592: MatDestroy(&mglevels[levels - 1]->R);
593: for (i = 0; i < levels - 1; i++) {
594: MatDestroy(&mglevels[i]->R);
595: MatDestroy(&mglevels[i]->B);
596: MatDestroy(&mglevels[i]->X);
597: }
598: }
599: }
601: if (mg->am == PC_MG_MULTIPLICATIVE) {
602: if (matapp) MatZeroEntries(X);
603: else VecZeroEntries(x);
604: for (i = 0; i < mg->cyclesperpcapply; i++) PCMGMCycle_Private(pc, mglevels + levels - 1, transpose, matapp, NULL);
605: } else if (mg->am == PC_MG_ADDITIVE) {
606: PCMGACycle_Private(pc, mglevels, transpose, matapp);
607: } else if (mg->am == PC_MG_KASKADE) {
608: PCMGKCycle_Private(pc, mglevels, transpose, matapp);
609: } else {
610: PCMGFCycle_Private(pc, mglevels, transpose, matapp);
611: }
612: if (mg->stageApply) PetscLogStagePop();
613: if (!changeu && !changed) {
614: if (matapp) {
615: mglevels[levels - 1]->B = NULL;
616: } else {
617: mglevels[levels - 1]->b = NULL;
618: }
619: }
620: return 0;
621: }
623: static PetscErrorCode PCApply_MG(PC pc, Vec b, Vec x)
624: {
625: PCApply_MG_Internal(pc, b, x, NULL, NULL, PETSC_FALSE);
626: return 0;
627: }
629: static PetscErrorCode PCApplyTranspose_MG(PC pc, Vec b, Vec x)
630: {
631: PCApply_MG_Internal(pc, b, x, NULL, NULL, PETSC_TRUE);
632: return 0;
633: }
635: static PetscErrorCode PCMatApply_MG(PC pc, Mat b, Mat x)
636: {
637: PCApply_MG_Internal(pc, NULL, NULL, b, x, PETSC_FALSE);
638: return 0;
639: }
641: PetscErrorCode PCSetFromOptions_MG(PC pc, PetscOptionItems *PetscOptionsObject)
642: {
643: PetscInt levels, cycles;
644: PetscBool flg, flg2;
645: PC_MG *mg = (PC_MG *)pc->data;
646: PC_MG_Levels **mglevels;
647: PCMGType mgtype;
648: PCMGCycleType mgctype;
649: PCMGGalerkinType gtype;
650: PCMGCoarseSpaceType coarseSpaceType;
652: levels = PetscMax(mg->nlevels, 1);
653: PetscOptionsHeadBegin(PetscOptionsObject, "Multigrid options");
654: PetscOptionsInt("-pc_mg_levels", "Number of Levels", "PCMGSetLevels", levels, &levels, &flg);
655: if (!flg && !mg->levels && pc->dm) {
656: DMGetRefineLevel(pc->dm, &levels);
657: levels++;
658: mg->usedmfornumberoflevels = PETSC_TRUE;
659: }
660: PCMGSetLevels(pc, levels, NULL);
661: mglevels = mg->levels;
663: mgctype = (PCMGCycleType)mglevels[0]->cycles;
664: PetscOptionsEnum("-pc_mg_cycle_type", "V cycle or for W-cycle", "PCMGSetCycleType", PCMGCycleTypes, (PetscEnum)mgctype, (PetscEnum *)&mgctype, &flg);
665: if (flg) PCMGSetCycleType(pc, mgctype);
666: gtype = mg->galerkin;
667: PetscOptionsEnum("-pc_mg_galerkin", "Use Galerkin process to compute coarser operators", "PCMGSetGalerkin", PCMGGalerkinTypes, (PetscEnum)gtype, (PetscEnum *)>ype, &flg);
668: if (flg) PCMGSetGalerkin(pc, gtype);
669: coarseSpaceType = mg->coarseSpaceType;
670: PetscOptionsEnum("-pc_mg_adapt_interp_coarse_space", "Type of adaptive coarse space: none, polynomial, harmonic, eigenvector, generalized_eigenvector, gdsw", "PCMGSetAdaptCoarseSpaceType", PCMGCoarseSpaceTypes, (PetscEnum)coarseSpaceType, (PetscEnum *)&coarseSpaceType, &flg);
671: if (flg) PCMGSetAdaptCoarseSpaceType(pc, coarseSpaceType);
672: PetscOptionsInt("-pc_mg_adapt_interp_n", "Size of the coarse space for adaptive interpolation", "PCMGSetCoarseSpace", mg->Nc, &mg->Nc, &flg);
673: PetscOptionsBool("-pc_mg_mesp_monitor", "Monitor the multilevel eigensolver", "PCMGSetAdaptInterpolation", PETSC_FALSE, &mg->mespMonitor, &flg);
674: flg2 = PETSC_FALSE;
675: PetscOptionsBool("-pc_mg_adapt_cr", "Monitor coarse space quality using Compatible Relaxation (CR)", "PCMGSetAdaptCR", PETSC_FALSE, &flg2, &flg);
676: if (flg) PCMGSetAdaptCR(pc, flg2);
677: flg = PETSC_FALSE;
678: PetscOptionsBool("-pc_mg_distinct_smoothup", "Create separate smoothup KSP and append the prefix _up", "PCMGSetDistinctSmoothUp", PETSC_FALSE, &flg, NULL);
679: if (flg) PCMGSetDistinctSmoothUp(pc);
680: mgtype = mg->am;
681: PetscOptionsEnum("-pc_mg_type", "Multigrid type", "PCMGSetType", PCMGTypes, (PetscEnum)mgtype, (PetscEnum *)&mgtype, &flg);
682: if (flg) PCMGSetType(pc, mgtype);
683: if (mg->am == PC_MG_MULTIPLICATIVE) {
684: PetscOptionsInt("-pc_mg_multiplicative_cycles", "Number of cycles for each preconditioner step", "PCMGMultiplicativeSetCycles", mg->cyclesperpcapply, &cycles, &flg);
685: if (flg) PCMGMultiplicativeSetCycles(pc, cycles);
686: }
687: flg = PETSC_FALSE;
688: PetscOptionsBool("-pc_mg_log", "Log times for each multigrid level", "None", flg, &flg, NULL);
689: if (flg) {
690: PetscInt i;
691: char eventname[128];
693: levels = mglevels[0]->levels;
694: for (i = 0; i < levels; i++) {
695: sprintf(eventname, "MGSetup Level %d", (int)i);
696: PetscLogEventRegister(eventname, ((PetscObject)pc)->classid, &mglevels[i]->eventsmoothsetup);
697: sprintf(eventname, "MGSmooth Level %d", (int)i);
698: PetscLogEventRegister(eventname, ((PetscObject)pc)->classid, &mglevels[i]->eventsmoothsolve);
699: if (i) {
700: sprintf(eventname, "MGResid Level %d", (int)i);
701: PetscLogEventRegister(eventname, ((PetscObject)pc)->classid, &mglevels[i]->eventresidual);
702: sprintf(eventname, "MGInterp Level %d", (int)i);
703: PetscLogEventRegister(eventname, ((PetscObject)pc)->classid, &mglevels[i]->eventinterprestrict);
704: }
705: }
707: #if defined(PETSC_USE_LOG)
708: {
709: const char *sname = "MG Apply";
710: PetscStageLog stageLog;
711: PetscInt st;
713: PetscLogGetStageLog(&stageLog);
714: for (st = 0; st < stageLog->numStages; ++st) {
715: PetscBool same;
717: PetscStrcmp(stageLog->stageInfo[st].name, sname, &same);
718: if (same) mg->stageApply = st;
719: }
720: if (!mg->stageApply) PetscLogStageRegister(sname, &mg->stageApply);
721: }
722: #endif
723: }
724: PetscOptionsHeadEnd();
725: /* Check option consistency */
726: PCMGGetGalerkin(pc, >ype);
727: PCMGGetAdaptInterpolation(pc, &flg);
729: return 0;
730: }
732: const char *const PCMGTypes[] = {"MULTIPLICATIVE", "ADDITIVE", "FULL", "KASKADE", "PCMGType", "PC_MG", NULL};
733: const char *const PCMGCycleTypes[] = {"invalid", "v", "w", "PCMGCycleType", "PC_MG_CYCLE", NULL};
734: const char *const PCMGGalerkinTypes[] = {"both", "pmat", "mat", "none", "external", "PCMGGalerkinType", "PC_MG_GALERKIN", NULL};
735: const char *const PCMGCoarseSpaceTypes[] = {"none", "polynomial", "harmonic", "eigenvector", "generalized_eigenvector", "gdsw", "PCMGCoarseSpaceType", "PCMG_ADAPT_NONE", NULL};
737: #include <petscdraw.h>
738: PetscErrorCode PCView_MG(PC pc, PetscViewer viewer)
739: {
740: PC_MG *mg = (PC_MG *)pc->data;
741: PC_MG_Levels **mglevels = mg->levels;
742: PetscInt levels = mglevels ? mglevels[0]->levels : 0, i;
743: PetscBool iascii, isbinary, isdraw;
745: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
746: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary);
747: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
748: if (iascii) {
749: const char *cyclename = levels ? (mglevels[0]->cycles == PC_MG_CYCLE_V ? "v" : "w") : "unknown";
750: PetscViewerASCIIPrintf(viewer, " type is %s, levels=%" PetscInt_FMT " cycles=%s\n", PCMGTypes[mg->am], levels, cyclename);
751: if (mg->am == PC_MG_MULTIPLICATIVE) PetscViewerASCIIPrintf(viewer, " Cycles per PCApply=%" PetscInt_FMT "\n", mg->cyclesperpcapply);
752: if (mg->galerkin == PC_MG_GALERKIN_BOTH) {
753: PetscViewerASCIIPrintf(viewer, " Using Galerkin computed coarse grid matrices\n");
754: } else if (mg->galerkin == PC_MG_GALERKIN_PMAT) {
755: PetscViewerASCIIPrintf(viewer, " Using Galerkin computed coarse grid matrices for pmat\n");
756: } else if (mg->galerkin == PC_MG_GALERKIN_MAT) {
757: PetscViewerASCIIPrintf(viewer, " Using Galerkin computed coarse grid matrices for mat\n");
758: } else if (mg->galerkin == PC_MG_GALERKIN_EXTERNAL) {
759: PetscViewerASCIIPrintf(viewer, " Using externally compute Galerkin coarse grid matrices\n");
760: } else {
761: PetscViewerASCIIPrintf(viewer, " Not using Galerkin computed coarse grid matrices\n");
762: }
763: if (mg->view) (*mg->view)(pc, viewer);
764: for (i = 0; i < levels; i++) {
765: if (i) {
766: PetscViewerASCIIPrintf(viewer, "Down solver (pre-smoother) on level %" PetscInt_FMT " -------------------------------\n", i);
767: } else {
768: PetscViewerASCIIPrintf(viewer, "Coarse grid solver -- level %" PetscInt_FMT " -------------------------------\n", i);
769: }
770: PetscViewerASCIIPushTab(viewer);
771: KSPView(mglevels[i]->smoothd, viewer);
772: PetscViewerASCIIPopTab(viewer);
773: if (i && mglevels[i]->smoothd == mglevels[i]->smoothu) {
774: PetscViewerASCIIPrintf(viewer, "Up solver (post-smoother) same as down solver (pre-smoother)\n");
775: } else if (i) {
776: PetscViewerASCIIPrintf(viewer, "Up solver (post-smoother) on level %" PetscInt_FMT " -------------------------------\n", i);
777: PetscViewerASCIIPushTab(viewer);
778: KSPView(mglevels[i]->smoothu, viewer);
779: PetscViewerASCIIPopTab(viewer);
780: }
781: if (i && mglevels[i]->cr) {
782: PetscViewerASCIIPrintf(viewer, "CR solver on level %" PetscInt_FMT " -------------------------------\n", i);
783: PetscViewerASCIIPushTab(viewer);
784: KSPView(mglevels[i]->cr, viewer);
785: PetscViewerASCIIPopTab(viewer);
786: }
787: }
788: } else if (isbinary) {
789: for (i = levels - 1; i >= 0; i--) {
790: KSPView(mglevels[i]->smoothd, viewer);
791: if (i && mglevels[i]->smoothd != mglevels[i]->smoothu) KSPView(mglevels[i]->smoothu, viewer);
792: }
793: } else if (isdraw) {
794: PetscDraw draw;
795: PetscReal x, w, y, bottom, th;
796: PetscViewerDrawGetDraw(viewer, 0, &draw);
797: PetscDrawGetCurrentPoint(draw, &x, &y);
798: PetscDrawStringGetSize(draw, NULL, &th);
799: bottom = y - th;
800: for (i = levels - 1; i >= 0; i--) {
801: if (!mglevels[i]->smoothu || (mglevels[i]->smoothu == mglevels[i]->smoothd)) {
802: PetscDrawPushCurrentPoint(draw, x, bottom);
803: KSPView(mglevels[i]->smoothd, viewer);
804: PetscDrawPopCurrentPoint(draw);
805: } else {
806: w = 0.5 * PetscMin(1.0 - x, x);
807: PetscDrawPushCurrentPoint(draw, x + w, bottom);
808: KSPView(mglevels[i]->smoothd, viewer);
809: PetscDrawPopCurrentPoint(draw);
810: PetscDrawPushCurrentPoint(draw, x - w, bottom);
811: KSPView(mglevels[i]->smoothu, viewer);
812: PetscDrawPopCurrentPoint(draw);
813: }
814: PetscDrawGetBoundingBox(draw, NULL, &bottom, NULL, NULL);
815: bottom -= th;
816: }
817: }
818: return 0;
819: }
821: #include <petsc/private/kspimpl.h>
823: /*
824: Calls setup for the KSP on each level
825: */
826: PetscErrorCode PCSetUp_MG(PC pc)
827: {
828: PC_MG *mg = (PC_MG *)pc->data;
829: PC_MG_Levels **mglevels = mg->levels;
830: PetscInt i, n;
831: PC cpc;
832: PetscBool dump = PETSC_FALSE, opsset, use_amat, missinginterpolate = PETSC_FALSE;
833: Mat dA, dB;
834: Vec tvec;
835: DM *dms;
836: PetscViewer viewer = NULL;
837: PetscBool dAeqdB = PETSC_FALSE, needRestricts = PETSC_FALSE, doCR = PETSC_FALSE;
838: PetscBool adaptInterpolation = mg->adaptInterpolation;
841: n = mglevels[0]->levels;
842: /* FIX: Move this to PCSetFromOptions_MG? */
843: if (mg->usedmfornumberoflevels) {
844: PetscInt levels;
845: DMGetRefineLevel(pc->dm, &levels);
846: levels++;
847: if (levels > n) { /* the problem is now being solved on a finer grid */
848: PCMGSetLevels(pc, levels, NULL);
849: n = levels;
850: PCSetFromOptions(pc); /* it is bad to call this here, but otherwise will never be called for the new hierarchy */
851: mglevels = mg->levels;
852: }
853: }
854: KSPGetPC(mglevels[0]->smoothd, &cpc);
856: /* If user did not provide fine grid operators OR operator was not updated since last global KSPSetOperators() */
857: /* so use those from global PC */
858: /* Is this what we always want? What if user wants to keep old one? */
859: KSPGetOperatorsSet(mglevels[n - 1]->smoothd, NULL, &opsset);
860: if (opsset) {
861: Mat mmat;
862: KSPGetOperators(mglevels[n - 1]->smoothd, NULL, &mmat);
863: if (mmat == pc->pmat) opsset = PETSC_FALSE;
864: }
866: /* Create CR solvers */
867: PCMGGetAdaptCR(pc, &doCR);
868: if (doCR) {
869: const char *prefix;
871: PCGetOptionsPrefix(pc, &prefix);
872: for (i = 1; i < n; ++i) {
873: PC ipc, cr;
874: char crprefix[128];
876: KSPCreate(PetscObjectComm((PetscObject)pc), &mglevels[i]->cr);
877: KSPSetErrorIfNotConverged(mglevels[i]->cr, PETSC_FALSE);
878: PetscObjectIncrementTabLevel((PetscObject)mglevels[i]->cr, (PetscObject)pc, n - i);
879: KSPSetOptionsPrefix(mglevels[i]->cr, prefix);
880: PetscObjectComposedDataSetInt((PetscObject)mglevels[i]->cr, PetscMGLevelId, mglevels[i]->level);
881: KSPSetType(mglevels[i]->cr, KSPCHEBYSHEV);
882: KSPSetConvergenceTest(mglevels[i]->cr, KSPConvergedSkip, NULL, NULL);
883: KSPSetNormType(mglevels[i]->cr, KSP_NORM_PRECONDITIONED);
884: KSPGetPC(mglevels[i]->cr, &ipc);
886: PCSetType(ipc, PCCOMPOSITE);
887: PCCompositeSetType(ipc, PC_COMPOSITE_MULTIPLICATIVE);
888: PCCompositeAddPCType(ipc, PCSOR);
889: CreateCR_Private(pc, i, &cr);
890: PCCompositeAddPC(ipc, cr);
891: PCDestroy(&cr);
893: KSPSetTolerances(mglevels[i]->cr, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, mg->default_smoothd);
894: KSPSetInitialGuessNonzero(mglevels[i]->cr, PETSC_TRUE);
895: PetscSNPrintf(crprefix, 128, "mg_levels_%d_cr_", (int)i);
896: KSPAppendOptionsPrefix(mglevels[i]->cr, crprefix);
897: }
898: }
900: if (!opsset) {
901: PCGetUseAmat(pc, &use_amat);
902: if (use_amat) {
903: PetscInfo(pc, "Using outer operators to define finest grid operator \n because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");
904: KSPSetOperators(mglevels[n - 1]->smoothd, pc->mat, pc->pmat);
905: } else {
906: PetscInfo(pc, "Using matrix (pmat) operators to define finest grid operator \n because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");
907: KSPSetOperators(mglevels[n - 1]->smoothd, pc->pmat, pc->pmat);
908: }
909: }
911: for (i = n - 1; i > 0; i--) {
912: if (!(mglevels[i]->interpolate || mglevels[i]->restrct)) {
913: missinginterpolate = PETSC_TRUE;
914: break;
915: }
916: }
918: KSPGetOperators(mglevels[n - 1]->smoothd, &dA, &dB);
919: if (dA == dB) dAeqdB = PETSC_TRUE;
920: if (mg->galerkin == PC_MG_GALERKIN_NONE || ((mg->galerkin == PC_MG_GALERKIN_PMAT || mg->galerkin == PC_MG_GALERKIN_MAT) && !dAeqdB)) {
921: needRestricts = PETSC_TRUE; /* user must compute either mat, pmat, or both so must restrict x to coarser levels */
922: }
924: if (pc->dm && !pc->setupcalled) {
925: /* finest smoother also gets DM but it is not active, independent of whether galerkin==PC_MG_GALERKIN_EXTERNAL */
926: KSPSetDM(mglevels[n - 1]->smoothd, pc->dm);
927: KSPSetDMActive(mglevels[n - 1]->smoothd, PETSC_FALSE);
928: if (mglevels[n - 1]->smoothd != mglevels[n - 1]->smoothu) {
929: KSPSetDM(mglevels[n - 1]->smoothu, pc->dm);
930: KSPSetDMActive(mglevels[n - 1]->smoothu, PETSC_FALSE);
931: }
932: if (mglevels[n - 1]->cr) {
933: KSPSetDM(mglevels[n - 1]->cr, pc->dm);
934: KSPSetDMActive(mglevels[n - 1]->cr, PETSC_FALSE);
935: }
936: }
938: /*
939: Skipping if user has provided all interpolation/restriction needed (since DM might not be able to produce them (when coming from SNES/TS)
940: Skipping for externally managed hierarchy (such as ML and GAMG). Cleaner logic here would be great. Wrap ML/GAMG as DMs?
941: */
942: if (missinginterpolate && mg->galerkin != PC_MG_GALERKIN_EXTERNAL && !pc->setupcalled) {
943: /* first see if we can compute a coarse space */
944: if (mg->coarseSpaceType == PCMG_ADAPT_GDSW) {
945: for (i = n - 2; i > -1; i--) {
946: if (!mglevels[i + 1]->restrct && !mglevels[i + 1]->interpolate) {
947: PCMGComputeCoarseSpace_Internal(pc, i + 1, mg->coarseSpaceType, mg->Nc, NULL, &mglevels[i + 1]->coarseSpace);
948: PCMGSetInterpolation(pc, i + 1, mglevels[i + 1]->coarseSpace);
949: }
950: }
951: } else { /* construct the interpolation from the DMs */
952: Mat p;
953: Vec rscale;
954: PetscMalloc1(n, &dms);
955: dms[n - 1] = pc->dm;
956: /* Separately create them so we do not get DMKSP interference between levels */
957: for (i = n - 2; i > -1; i--) DMCoarsen(dms[i + 1], MPI_COMM_NULL, &dms[i]);
958: for (i = n - 2; i > -1; i--) {
959: DMKSP kdm;
960: PetscBool dmhasrestrict, dmhasinject;
961: KSPSetDM(mglevels[i]->smoothd, dms[i]);
962: if (!needRestricts) KSPSetDMActive(mglevels[i]->smoothd, PETSC_FALSE);
963: if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
964: KSPSetDM(mglevels[i]->smoothu, dms[i]);
965: if (!needRestricts) KSPSetDMActive(mglevels[i]->smoothu, PETSC_FALSE);
966: }
967: if (mglevels[i]->cr) {
968: KSPSetDM(mglevels[i]->cr, dms[i]);
969: if (!needRestricts) KSPSetDMActive(mglevels[i]->cr, PETSC_FALSE);
970: }
971: DMGetDMKSPWrite(dms[i], &kdm);
972: /* Ugly hack so that the next KSPSetUp() will use the RHS that we set. A better fix is to change dmActive to take
973: * a bitwise OR of computing the matrix, RHS, and initial iterate. */
974: kdm->ops->computerhs = NULL;
975: kdm->rhsctx = NULL;
976: if (!mglevels[i + 1]->interpolate) {
977: DMCreateInterpolation(dms[i], dms[i + 1], &p, &rscale);
978: PCMGSetInterpolation(pc, i + 1, p);
979: if (rscale) PCMGSetRScale(pc, i + 1, rscale);
980: VecDestroy(&rscale);
981: MatDestroy(&p);
982: }
983: DMHasCreateRestriction(dms[i], &dmhasrestrict);
984: if (dmhasrestrict && !mglevels[i + 1]->restrct) {
985: DMCreateRestriction(dms[i], dms[i + 1], &p);
986: PCMGSetRestriction(pc, i + 1, p);
987: MatDestroy(&p);
988: }
989: DMHasCreateInjection(dms[i], &dmhasinject);
990: if (dmhasinject && !mglevels[i + 1]->inject) {
991: DMCreateInjection(dms[i], dms[i + 1], &p);
992: PCMGSetInjection(pc, i + 1, p);
993: MatDestroy(&p);
994: }
995: }
997: for (i = n - 2; i > -1; i--) DMDestroy(&dms[i]);
998: PetscFree(dms);
999: }
1000: }
1002: if (mg->galerkin < PC_MG_GALERKIN_NONE) {
1003: Mat A, B;
1004: PetscBool doA = PETSC_FALSE, doB = PETSC_FALSE;
1005: MatReuse reuse = MAT_INITIAL_MATRIX;
1007: if (mg->galerkin == PC_MG_GALERKIN_PMAT || mg->galerkin == PC_MG_GALERKIN_BOTH) doB = PETSC_TRUE;
1008: if (mg->galerkin == PC_MG_GALERKIN_MAT || (mg->galerkin == PC_MG_GALERKIN_BOTH && dA != dB)) doA = PETSC_TRUE;
1009: if (pc->setupcalled) reuse = MAT_REUSE_MATRIX;
1010: for (i = n - 2; i > -1; i--) {
1012: if (!mglevels[i + 1]->interpolate) PCMGSetInterpolation(pc, i + 1, mglevels[i + 1]->restrct);
1013: if (!mglevels[i + 1]->restrct) PCMGSetRestriction(pc, i + 1, mglevels[i + 1]->interpolate);
1014: if (reuse == MAT_REUSE_MATRIX) KSPGetOperators(mglevels[i]->smoothd, &A, &B);
1015: if (doA) MatGalerkin(mglevels[i + 1]->restrct, dA, mglevels[i + 1]->interpolate, reuse, 1.0, &A);
1016: if (doB) MatGalerkin(mglevels[i + 1]->restrct, dB, mglevels[i + 1]->interpolate, reuse, 1.0, &B);
1017: /* the management of the PetscObjectReference() and PetscObjecDereference() below is rather delicate */
1018: if (!doA && dAeqdB) {
1019: if (reuse == MAT_INITIAL_MATRIX) PetscObjectReference((PetscObject)B);
1020: A = B;
1021: } else if (!doA && reuse == MAT_INITIAL_MATRIX) {
1022: KSPGetOperators(mglevels[i]->smoothd, &A, NULL);
1023: PetscObjectReference((PetscObject)A);
1024: }
1025: if (!doB && dAeqdB) {
1026: if (reuse == MAT_INITIAL_MATRIX) PetscObjectReference((PetscObject)A);
1027: B = A;
1028: } else if (!doB && reuse == MAT_INITIAL_MATRIX) {
1029: KSPGetOperators(mglevels[i]->smoothd, NULL, &B);
1030: PetscObjectReference((PetscObject)B);
1031: }
1032: if (reuse == MAT_INITIAL_MATRIX) {
1033: KSPSetOperators(mglevels[i]->smoothd, A, B);
1034: PetscObjectDereference((PetscObject)A);
1035: PetscObjectDereference((PetscObject)B);
1036: }
1037: dA = A;
1038: dB = B;
1039: }
1040: }
1042: /* Adapt interpolation matrices */
1043: if (adaptInterpolation) {
1044: for (i = 0; i < n; ++i) {
1045: if (!mglevels[i]->coarseSpace) PCMGComputeCoarseSpace_Internal(pc, i, mg->coarseSpaceType, mg->Nc, !i ? NULL : mglevels[i - 1]->coarseSpace, &mglevels[i]->coarseSpace);
1046: if (i) PCMGAdaptInterpolator_Internal(pc, i, mglevels[i - 1]->smoothu, mglevels[i]->smoothu, mglevels[i - 1]->coarseSpace, mglevels[i]->coarseSpace);
1047: }
1048: for (i = n - 2; i > -1; --i) PCMGRecomputeLevelOperators_Internal(pc, i);
1049: }
1051: if (needRestricts && pc->dm) {
1052: for (i = n - 2; i >= 0; i--) {
1053: DM dmfine, dmcoarse;
1054: Mat Restrict, Inject;
1055: Vec rscale;
1056: KSPGetDM(mglevels[i + 1]->smoothd, &dmfine);
1057: KSPGetDM(mglevels[i]->smoothd, &dmcoarse);
1058: PCMGGetRestriction(pc, i + 1, &Restrict);
1059: PCMGGetRScale(pc, i + 1, &rscale);
1060: PCMGGetInjection(pc, i + 1, &Inject);
1061: DMRestrict(dmfine, Restrict, rscale, Inject, dmcoarse);
1062: }
1063: }
1065: if (!pc->setupcalled) {
1066: for (i = 0; i < n; i++) KSPSetFromOptions(mglevels[i]->smoothd);
1067: for (i = 1; i < n; i++) {
1068: if (mglevels[i]->smoothu && (mglevels[i]->smoothu != mglevels[i]->smoothd)) KSPSetFromOptions(mglevels[i]->smoothu);
1069: if (mglevels[i]->cr) KSPSetFromOptions(mglevels[i]->cr);
1070: }
1071: /* insure that if either interpolation or restriction is set the other other one is set */
1072: for (i = 1; i < n; i++) {
1073: PCMGGetInterpolation(pc, i, NULL);
1074: PCMGGetRestriction(pc, i, NULL);
1075: }
1076: for (i = 0; i < n - 1; i++) {
1077: if (!mglevels[i]->b) {
1078: Vec *vec;
1079: KSPCreateVecs(mglevels[i]->smoothd, 1, &vec, 0, NULL);
1080: PCMGSetRhs(pc, i, *vec);
1081: VecDestroy(vec);
1082: PetscFree(vec);
1083: }
1084: if (!mglevels[i]->r && i) {
1085: VecDuplicate(mglevels[i]->b, &tvec);
1086: PCMGSetR(pc, i, tvec);
1087: VecDestroy(&tvec);
1088: }
1089: if (!mglevels[i]->x) {
1090: VecDuplicate(mglevels[i]->b, &tvec);
1091: PCMGSetX(pc, i, tvec);
1092: VecDestroy(&tvec);
1093: }
1094: if (doCR) {
1095: VecDuplicate(mglevels[i]->b, &mglevels[i]->crx);
1096: VecDuplicate(mglevels[i]->b, &mglevels[i]->crb);
1097: VecZeroEntries(mglevels[i]->crb);
1098: }
1099: }
1100: if (n != 1 && !mglevels[n - 1]->r) {
1101: /* PCMGSetR() on the finest level if user did not supply it */
1102: Vec *vec;
1103: KSPCreateVecs(mglevels[n - 1]->smoothd, 1, &vec, 0, NULL);
1104: PCMGSetR(pc, n - 1, *vec);
1105: VecDestroy(vec);
1106: PetscFree(vec);
1107: }
1108: if (doCR) {
1109: VecDuplicate(mglevels[n - 1]->r, &mglevels[n - 1]->crx);
1110: VecDuplicate(mglevels[n - 1]->r, &mglevels[n - 1]->crb);
1111: VecZeroEntries(mglevels[n - 1]->crb);
1112: }
1113: }
1115: if (pc->dm) {
1116: /* need to tell all the coarser levels to rebuild the matrix using the DM for that level */
1117: for (i = 0; i < n - 1; i++) {
1118: if (mglevels[i]->smoothd->setupstage != KSP_SETUP_NEW) mglevels[i]->smoothd->setupstage = KSP_SETUP_NEWMATRIX;
1119: }
1120: }
1121: // We got here (PCSetUp_MG) because the matrix has changed, which means the smoother needs to be set up again (e.g.,
1122: // new diagonal for Jacobi). Setting it here allows it to be logged under PCSetUp rather than deep inside a PCApply.
1123: if (mglevels[n - 1]->smoothd->setupstage != KSP_SETUP_NEW) mglevels[n - 1]->smoothd->setupstage = KSP_SETUP_NEWMATRIX;
1125: for (i = 1; i < n; i++) {
1126: if (mglevels[i]->smoothu == mglevels[i]->smoothd || mg->am == PC_MG_FULL || mg->am == PC_MG_KASKADE || mg->cyclesperpcapply > 1) {
1127: /* if doing only down then initial guess is zero */
1128: KSPSetInitialGuessNonzero(mglevels[i]->smoothd, PETSC_TRUE);
1129: }
1130: if (mglevels[i]->cr) KSPSetInitialGuessNonzero(mglevels[i]->cr, PETSC_TRUE);
1131: if (mglevels[i]->eventsmoothsetup) PetscLogEventBegin(mglevels[i]->eventsmoothsetup, 0, 0, 0, 0);
1132: KSPSetUp(mglevels[i]->smoothd);
1133: if (mglevels[i]->smoothd->reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
1134: if (mglevels[i]->eventsmoothsetup) PetscLogEventEnd(mglevels[i]->eventsmoothsetup, 0, 0, 0, 0);
1135: if (!mglevels[i]->residual) {
1136: Mat mat;
1137: KSPGetOperators(mglevels[i]->smoothd, &mat, NULL);
1138: PCMGSetResidual(pc, i, PCMGResidualDefault, mat);
1139: }
1140: if (!mglevels[i]->residualtranspose) {
1141: Mat mat;
1142: KSPGetOperators(mglevels[i]->smoothd, &mat, NULL);
1143: PCMGSetResidualTranspose(pc, i, PCMGResidualTransposeDefault, mat);
1144: }
1145: }
1146: for (i = 1; i < n; i++) {
1147: if (mglevels[i]->smoothu && mglevels[i]->smoothu != mglevels[i]->smoothd) {
1148: Mat downmat, downpmat;
1150: /* check if operators have been set for up, if not use down operators to set them */
1151: KSPGetOperatorsSet(mglevels[i]->smoothu, &opsset, NULL);
1152: if (!opsset) {
1153: KSPGetOperators(mglevels[i]->smoothd, &downmat, &downpmat);
1154: KSPSetOperators(mglevels[i]->smoothu, downmat, downpmat);
1155: }
1157: KSPSetInitialGuessNonzero(mglevels[i]->smoothu, PETSC_TRUE);
1158: if (mglevels[i]->eventsmoothsetup) PetscLogEventBegin(mglevels[i]->eventsmoothsetup, 0, 0, 0, 0);
1159: KSPSetUp(mglevels[i]->smoothu);
1160: if (mglevels[i]->smoothu->reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
1161: if (mglevels[i]->eventsmoothsetup) PetscLogEventEnd(mglevels[i]->eventsmoothsetup, 0, 0, 0, 0);
1162: }
1163: if (mglevels[i]->cr) {
1164: Mat downmat, downpmat;
1166: /* check if operators have been set for up, if not use down operators to set them */
1167: KSPGetOperatorsSet(mglevels[i]->cr, &opsset, NULL);
1168: if (!opsset) {
1169: KSPGetOperators(mglevels[i]->smoothd, &downmat, &downpmat);
1170: KSPSetOperators(mglevels[i]->cr, downmat, downpmat);
1171: }
1173: KSPSetInitialGuessNonzero(mglevels[i]->cr, PETSC_TRUE);
1174: if (mglevels[i]->eventsmoothsetup) PetscLogEventBegin(mglevels[i]->eventsmoothsetup, 0, 0, 0, 0);
1175: KSPSetUp(mglevels[i]->cr);
1176: if (mglevels[i]->cr->reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
1177: if (mglevels[i]->eventsmoothsetup) PetscLogEventEnd(mglevels[i]->eventsmoothsetup, 0, 0, 0, 0);
1178: }
1179: }
1181: if (mglevels[0]->eventsmoothsetup) PetscLogEventBegin(mglevels[0]->eventsmoothsetup, 0, 0, 0, 0);
1182: KSPSetUp(mglevels[0]->smoothd);
1183: if (mglevels[0]->smoothd->reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
1184: if (mglevels[0]->eventsmoothsetup) PetscLogEventEnd(mglevels[0]->eventsmoothsetup, 0, 0, 0, 0);
1186: /*
1187: Dump the interpolation/restriction matrices plus the
1188: Jacobian/stiffness on each level. This allows MATLAB users to
1189: easily check if the Galerkin condition A_c = R A_f R^T is satisfied.
1191: Only support one or the other at the same time.
1192: */
1193: #if defined(PETSC_USE_SOCKET_VIEWER)
1194: PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_mg_dump_matlab", &dump, NULL);
1195: if (dump) viewer = PETSC_VIEWER_SOCKET_(PetscObjectComm((PetscObject)pc));
1196: dump = PETSC_FALSE;
1197: #endif
1198: PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_mg_dump_binary", &dump, NULL);
1199: if (dump) viewer = PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)pc));
1201: if (viewer) {
1202: for (i = 1; i < n; i++) MatView(mglevels[i]->restrct, viewer);
1203: for (i = 0; i < n; i++) {
1204: KSPGetPC(mglevels[i]->smoothd, &pc);
1205: MatView(pc->mat, viewer);
1206: }
1207: }
1208: return 0;
1209: }
1211: /* -------------------------------------------------------------------------------------*/
1213: PetscErrorCode PCMGGetLevels_MG(PC pc, PetscInt *levels)
1214: {
1215: PC_MG *mg = (PC_MG *)pc->data;
1217: *levels = mg->nlevels;
1218: return 0;
1219: }
1221: /*@
1222: PCMGGetLevels - Gets the number of levels to use with `PCMG`.
1224: Not Collective
1226: Input Parameter:
1227: . pc - the preconditioner context
1229: Output parameter:
1230: . levels - the number of levels
1232: Level: advanced
1234: .seealso: `PCMG`, `PCMGSetLevels()`
1235: @*/
1236: PetscErrorCode PCMGGetLevels(PC pc, PetscInt *levels)
1237: {
1240: *levels = 0;
1241: PetscTryMethod(pc, "PCMGGetLevels_C", (PC, PetscInt *), (pc, levels));
1242: return 0;
1243: }
1245: /*@
1246: PCMGGetGridComplexity - compute operator and grid complexity of the `PCMG` hierarchy
1248: Input Parameter:
1249: . pc - the preconditioner context
1251: Output Parameters:
1252: + gc - grid complexity = sum_i(n_i) / n_0
1253: - oc - operator complexity = sum_i(nnz_i) / nnz_0
1255: Level: advanced
1257: Note:
1258: This is often call the operator complexity in multigrid literature
1260: .seealso: `PCMG`, `PCMGGetLevels()`, `PCMGSetLevels()`
1261: @*/
1262: PetscErrorCode PCMGGetGridComplexity(PC pc, PetscReal *gc, PetscReal *oc)
1263: {
1264: PC_MG *mg = (PC_MG *)pc->data;
1265: PC_MG_Levels **mglevels = mg->levels;
1266: PetscInt lev, N;
1267: PetscLogDouble nnz0 = 0, sgc = 0, soc = 0, n0 = 0;
1268: MatInfo info;
1273: if (!pc->setupcalled) {
1274: if (gc) *gc = 0;
1275: if (oc) *oc = 0;
1276: return 0;
1277: }
1279: for (lev = 0; lev < mg->nlevels; lev++) {
1280: Mat dB;
1281: KSPGetOperators(mglevels[lev]->smoothd, NULL, &dB);
1282: MatGetInfo(dB, MAT_GLOBAL_SUM, &info); /* global reduction */
1283: MatGetSize(dB, &N, NULL);
1284: sgc += N;
1285: soc += info.nz_used;
1286: if (lev == mg->nlevels - 1) {
1287: nnz0 = info.nz_used;
1288: n0 = N;
1289: }
1290: }
1291: if (n0 > 0 && gc) *gc = (PetscReal)(sgc / n0);
1292: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number for grid points on finest level is not available");
1293: if (nnz0 > 0 && oc) *oc = (PetscReal)(soc / nnz0);
1294: return 0;
1295: }
1297: /*@
1298: PCMGSetType - Determines the form of multigrid to use:
1299: multiplicative, additive, full, or the Kaskade algorithm.
1301: Logically Collective
1303: Input Parameters:
1304: + pc - the preconditioner context
1305: - form - multigrid form, one of `PC_MG_MULTIPLICATIVE`, `PC_MG_ADDITIVE`, `PC_MG_FULL`, `PC_MG_KASKADE`
1307: Options Database Key:
1308: . -pc_mg_type <form> - Sets <form>, one of multiplicative,
1309: additive, full, kaskade
1311: Level: advanced
1313: .seealso: `PCMGType`, `PCMG`, `PCMGGetLevels()`, `PCMGSetLevels()`, `PCMGGetType()`, `PCMGCycleType`
1314: @*/
1315: PetscErrorCode PCMGSetType(PC pc, PCMGType form)
1316: {
1317: PC_MG *mg = (PC_MG *)pc->data;
1321: mg->am = form;
1322: if (form == PC_MG_MULTIPLICATIVE) pc->ops->applyrichardson = PCApplyRichardson_MG;
1323: else pc->ops->applyrichardson = NULL;
1324: return 0;
1325: }
1327: /*@
1328: PCMGGetType - Finds the form of multigrid the `PCMG` is using multiplicative, additive, full, or the Kaskade algorithm.
1330: Logically Collective
1332: Input Parameter:
1333: . pc - the preconditioner context
1335: Output Parameter:
1336: . type - one of `PC_MG_MULTIPLICATIVE`, `PC_MG_ADDITIVE`, `PC_MG_FULL`, `PC_MG_KASKADE`, `PCMGCycleType`
1338: Level: advanced
1340: .seealso: `PCMGType`, `PCMG`, `PCMGGetLevels()`, `PCMGSetLevels()`, `PCMGSetType()`
1341: @*/
1342: PetscErrorCode PCMGGetType(PC pc, PCMGType *type)
1343: {
1344: PC_MG *mg = (PC_MG *)pc->data;
1347: *type = mg->am;
1348: return 0;
1349: }
1351: /*@
1352: PCMGSetCycleType - Sets the type cycles to use. Use `PCMGSetCycleTypeOnLevel()` for more
1353: complicated cycling.
1355: Logically Collective
1357: Input Parameters:
1358: + pc - the multigrid context
1359: - n - either `PC_MG_CYCLE_V` or `PC_MG_CYCLE_W`
1361: Options Database Key:
1362: . -pc_mg_cycle_type <v,w> - provide the cycle desired
1364: Level: advanced
1366: .seealso: `PCMG`, `PCMGSetCycleTypeOnLevel()`, `PCMGType`, `PCMGCycleType`
1367: @*/
1368: PetscErrorCode PCMGSetCycleType(PC pc, PCMGCycleType n)
1369: {
1370: PC_MG *mg = (PC_MG *)pc->data;
1371: PC_MG_Levels **mglevels = mg->levels;
1372: PetscInt i, levels;
1377: levels = mglevels[0]->levels;
1378: for (i = 0; i < levels; i++) mglevels[i]->cycles = n;
1379: return 0;
1380: }
1382: /*@
1383: PCMGMultiplicativeSetCycles - Sets the number of cycles to use for each preconditioner step
1384: of multigrid when `PCMGType` is `PC_MG_MULTIPLICATIVE`
1386: Logically Collective
1388: Input Parameters:
1389: + pc - the multigrid context
1390: - n - number of cycles (default is 1)
1392: Options Database Key:
1393: . -pc_mg_multiplicative_cycles n - set the number of cycles
1395: Level: advanced
1397: Note:
1398: This is not associated with setting a v or w cycle, that is set with `PCMGSetCycleType()`
1400: .seealso: `PCMGSetCycleTypeOnLevel()`, `PCMGSetCycleType()`, `PCMGCycleType`, `PCMGType`
1401: @*/
1402: PetscErrorCode PCMGMultiplicativeSetCycles(PC pc, PetscInt n)
1403: {
1404: PC_MG *mg = (PC_MG *)pc->data;
1408: mg->cyclesperpcapply = n;
1409: return 0;
1410: }
1412: PetscErrorCode PCMGSetGalerkin_MG(PC pc, PCMGGalerkinType use)
1413: {
1414: PC_MG *mg = (PC_MG *)pc->data;
1416: mg->galerkin = use;
1417: return 0;
1418: }
1420: /*@
1421: PCMGSetGalerkin - Causes the coarser grid matrices to be computed from the
1422: finest grid via the Galerkin process: A_i-1 = r_i * A_i * p_i
1424: Logically Collective
1426: Input Parameters:
1427: + pc - the multigrid context
1428: - use - one of `PC_MG_GALERKIN_BOTH`, `PC_MG_GALERKIN_PMAT`, `PC_MG_GALERKIN_MAT`, or `PC_MG_GALERKIN_NONE`
1430: Options Database Key:
1431: . -pc_mg_galerkin <both,pmat,mat,none> - set the matrices to form via the Galerkin process
1433: Level: intermediate
1435: Note:
1436: Some codes that use `PCMG` such as `PCGAMG` use Galerkin internally while constructing the hierarchy and thus do not
1437: use the `PCMG` construction of the coarser grids.
1439: .seealso: `PCMG`, `PCMGGetGalerkin()`, `PCMGGalerkinType`
1440: @*/
1441: PetscErrorCode PCMGSetGalerkin(PC pc, PCMGGalerkinType use)
1442: {
1444: PetscTryMethod(pc, "PCMGSetGalerkin_C", (PC, PCMGGalerkinType), (pc, use));
1445: return 0;
1446: }
1448: /*@
1449: PCMGGetGalerkin - Checks if Galerkin multigrid is being used, i.e. A_i-1 = r_i * A_i * p_i
1451: Not Collective
1453: Input Parameter:
1454: . pc - the multigrid context
1456: Output Parameter:
1457: . galerkin - one of `PC_MG_GALERKIN_BOTH`,`PC_MG_GALERKIN_PMAT`,`PC_MG_GALERKIN_MAT`, `PC_MG_GALERKIN_NONE`, or `PC_MG_GALERKIN_EXTERNAL`
1459: Level: intermediate
1461: .seealso: `PCMG`, `PCMGSetGalerkin()`, `PCMGGalerkinType`
1462: @*/
1463: PetscErrorCode PCMGGetGalerkin(PC pc, PCMGGalerkinType *galerkin)
1464: {
1465: PC_MG *mg = (PC_MG *)pc->data;
1468: *galerkin = mg->galerkin;
1469: return 0;
1470: }
1472: PetscErrorCode PCMGSetAdaptInterpolation_MG(PC pc, PetscBool adapt)
1473: {
1474: PC_MG *mg = (PC_MG *)pc->data;
1476: mg->adaptInterpolation = adapt;
1477: return 0;
1478: }
1480: PetscErrorCode PCMGGetAdaptInterpolation_MG(PC pc, PetscBool *adapt)
1481: {
1482: PC_MG *mg = (PC_MG *)pc->data;
1484: *adapt = mg->adaptInterpolation;
1485: return 0;
1486: }
1488: PetscErrorCode PCMGSetAdaptCoarseSpaceType_MG(PC pc, PCMGCoarseSpaceType ctype)
1489: {
1490: PC_MG *mg = (PC_MG *)pc->data;
1492: mg->adaptInterpolation = ctype != PCMG_ADAPT_NONE ? PETSC_TRUE : PETSC_FALSE;
1493: mg->coarseSpaceType = ctype;
1494: return 0;
1495: }
1497: PetscErrorCode PCMGGetAdaptCoarseSpaceType_MG(PC pc, PCMGCoarseSpaceType *ctype)
1498: {
1499: PC_MG *mg = (PC_MG *)pc->data;
1501: *ctype = mg->coarseSpaceType;
1502: return 0;
1503: }
1505: PetscErrorCode PCMGSetAdaptCR_MG(PC pc, PetscBool cr)
1506: {
1507: PC_MG *mg = (PC_MG *)pc->data;
1509: mg->compatibleRelaxation = cr;
1510: return 0;
1511: }
1513: PetscErrorCode PCMGGetAdaptCR_MG(PC pc, PetscBool *cr)
1514: {
1515: PC_MG *mg = (PC_MG *)pc->data;
1517: *cr = mg->compatibleRelaxation;
1518: return 0;
1519: }
1521: /*@C
1522: PCMGSetAdaptCoarseSpaceType - Set the type of adaptive coarse space.
1524: Adapts or creates the interpolator based upon a vector space which should be accurately captured by the next coarser mesh, and thus accurately interpolated.
1526: Logically Collective
1528: Input Parameters:
1529: + pc - the multigrid context
1530: - ctype - the type of coarse space
1532: Options Database Keys:
1533: + -pc_mg_adapt_interp_n <int> - The number of modes to use
1534: - -pc_mg_adapt_interp_coarse_space <type> - The type of coarse space: none, polynomial, harmonic, eigenvector, generalized_eigenvector, gdsw
1536: Level: intermediate
1538: .seealso: `PCMG`, `PCMGCoarseSpaceType`, `PCMGGetAdaptCoarseSpaceType()`, `PCMGSetGalerkin()`, `PCMGSetAdaptInterpolation()`
1539: @*/
1540: PetscErrorCode PCMGSetAdaptCoarseSpaceType(PC pc, PCMGCoarseSpaceType ctype)
1541: {
1544: PetscTryMethod(pc, "PCMGSetAdaptCoarseSpaceType_C", (PC, PCMGCoarseSpaceType), (pc, ctype));
1545: return 0;
1546: }
1548: /*@C
1549: PCMGGetAdaptCoarseSpaceType - Get the type of adaptive coarse space.
1551: Not Collective
1553: Input Parameter:
1554: . pc - the multigrid context
1556: Output Parameter:
1557: . ctype - the type of coarse space
1559: Level: intermediate
1561: .seealso: `PCMG`, `PCMGCoarseSpaceType`, `PCMGSetAdaptCoarseSpaceType()`, `PCMGSetGalerkin()`, `PCMGSetAdaptInterpolation()`
1562: @*/
1563: PetscErrorCode PCMGGetAdaptCoarseSpaceType(PC pc, PCMGCoarseSpaceType *ctype)
1564: {
1567: PetscUseMethod(pc, "PCMGGetAdaptCoarseSpaceType_C", (PC, PCMGCoarseSpaceType *), (pc, ctype));
1568: return 0;
1569: }
1571: /* MATT: REMOVE? */
1572: /*@
1573: PCMGSetAdaptInterpolation - Adapt the interpolator based upon a vector space which should be accurately captured by the next coarser mesh, and thus accurately interpolated.
1575: Logically Collective
1577: Input Parameters:
1578: + pc - the multigrid context
1579: - adapt - flag for adaptation of the interpolator
1581: Options Database Keys:
1582: + -pc_mg_adapt_interp - Turn on adaptation
1583: . -pc_mg_adapt_interp_n <int> - The number of modes to use, should be divisible by dimension
1584: - -pc_mg_adapt_interp_coarse_space <type> - The type of coarse space: polynomial, harmonic, eigenvector, generalized_eigenvector
1586: Level: intermediate
1588: .seealso: `PCMG`, `PCMGGetAdaptInterpolation()`, `PCMGSetGalerkin()`, `PCMGGetAdaptCoarseSpaceType()`, `PCMGSetAdaptCoarseSpaceType()`
1589: @*/
1590: PetscErrorCode PCMGSetAdaptInterpolation(PC pc, PetscBool adapt)
1591: {
1593: PetscTryMethod(pc, "PCMGSetAdaptInterpolation_C", (PC, PetscBool), (pc, adapt));
1594: return 0;
1595: }
1597: /*@
1598: PCMGGetAdaptInterpolation - Get the flag to adapt the interpolator based upon a vector space which should be accurately captured by the next coarser mesh,
1599: and thus accurately interpolated.
1601: Not collective
1603: Input Parameter:
1604: . pc - the multigrid context
1606: Output Parameter:
1607: . adapt - flag for adaptation of the interpolator
1609: Level: intermediate
1611: .seealso: `PCMG`, `PCMGSetAdaptInterpolation()`, `PCMGSetGalerkin()`, `PCMGGetAdaptCoarseSpaceType()`, `PCMGSetAdaptCoarseSpaceType()`
1612: @*/
1613: PetscErrorCode PCMGGetAdaptInterpolation(PC pc, PetscBool *adapt)
1614: {
1617: PetscUseMethod(pc, "PCMGGetAdaptInterpolation_C", (PC, PetscBool *), (pc, adapt));
1618: return 0;
1619: }
1621: /*@
1622: PCMGSetAdaptCR - Monitor the coarse space quality using an auxiliary solve with compatible relaxation.
1624: Logically Collective
1626: Input Parameters:
1627: + pc - the multigrid context
1628: - cr - flag for compatible relaxation
1630: Options Database Key:
1631: . -pc_mg_adapt_cr - Turn on compatible relaxation
1633: Level: intermediate
1635: .seealso: `PCMG`, `PCMGGetAdaptCR()`, `PCMGSetAdaptInterpolation()`, `PCMGSetGalerkin()`, `PCMGGetAdaptCoarseSpaceType()`, `PCMGSetAdaptCoarseSpaceType()`
1636: @*/
1637: PetscErrorCode PCMGSetAdaptCR(PC pc, PetscBool cr)
1638: {
1640: PetscTryMethod(pc, "PCMGSetAdaptCR_C", (PC, PetscBool), (pc, cr));
1641: return 0;
1642: }
1644: /*@
1645: PCMGGetAdaptCR - Get the flag to monitor coarse space quality using an auxiliary solve with compatible relaxation.
1647: Not collective
1649: Input Parameter:
1650: . pc - the multigrid context
1652: Output Parameter:
1653: . cr - flag for compatible relaxaion
1655: Level: intermediate
1657: .seealso: `PCMGSetAdaptCR()`, `PCMGGetAdaptInterpolation()`, `PCMGSetGalerkin()`, `PCMGGetAdaptCoarseSpaceType()`, `PCMGSetAdaptCoarseSpaceType()`
1658: @*/
1659: PetscErrorCode PCMGGetAdaptCR(PC pc, PetscBool *cr)
1660: {
1663: PetscUseMethod(pc, "PCMGGetAdaptCR_C", (PC, PetscBool *), (pc, cr));
1664: return 0;
1665: }
1667: /*@
1668: PCMGSetNumberSmooth - Sets the number of pre and post-smoothing steps to use
1669: on all levels. Use `PCMGDistinctSmoothUp()` to create separate up and down smoothers if you want different numbers of
1670: pre- and post-smoothing steps.
1672: Logically Collective
1674: Input Parameters:
1675: + mg - the multigrid context
1676: - n - the number of smoothing steps
1678: Options Database Key:
1679: . -mg_levels_ksp_max_it <n> - Sets number of pre and post-smoothing steps
1681: Level: advanced
1683: Note:
1684: This does not set a value on the coarsest grid, since we assume that there is no separate smooth up on the coarsest grid.
1686: .seealso: `PCMG`, `PCMGSetDistinctSmoothUp()`
1687: @*/
1688: PetscErrorCode PCMGSetNumberSmooth(PC pc, PetscInt n)
1689: {
1690: PC_MG *mg = (PC_MG *)pc->data;
1691: PC_MG_Levels **mglevels = mg->levels;
1692: PetscInt i, levels;
1697: levels = mglevels[0]->levels;
1699: for (i = 1; i < levels; i++) {
1700: KSPSetTolerances(mglevels[i]->smoothu, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, n);
1701: KSPSetTolerances(mglevels[i]->smoothd, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, n);
1702: mg->default_smoothu = n;
1703: mg->default_smoothd = n;
1704: }
1705: return 0;
1706: }
1708: /*@
1709: PCMGSetDistinctSmoothUp - sets the up (post) smoother to be a separate `KSP` from the down (pre) smoother on all levels
1710: and adds the suffix _up to the options name
1712: Logically Collective
1714: Input Parameter:
1715: . pc - the preconditioner context
1717: Options Database Key:
1718: . -pc_mg_distinct_smoothup <bool> - use distinct smoothing objects
1720: Level: advanced
1722: Note:
1723: This does not set a value on the coarsest grid, since we assume that there is no separate smooth up on the coarsest grid.
1725: .seealso: `PCMG`, `PCMGSetNumberSmooth()`
1726: @*/
1727: PetscErrorCode PCMGSetDistinctSmoothUp(PC pc)
1728: {
1729: PC_MG *mg = (PC_MG *)pc->data;
1730: PC_MG_Levels **mglevels = mg->levels;
1731: PetscInt i, levels;
1732: KSP subksp;
1736: levels = mglevels[0]->levels;
1738: for (i = 1; i < levels; i++) {
1739: const char *prefix = NULL;
1740: /* make sure smoother up and down are different */
1741: PCMGGetSmootherUp(pc, i, &subksp);
1742: KSPGetOptionsPrefix(mglevels[i]->smoothd, &prefix);
1743: KSPSetOptionsPrefix(subksp, prefix);
1744: KSPAppendOptionsPrefix(subksp, "up_");
1745: }
1746: return 0;
1747: }
1749: /* No new matrices are created, and the coarse operator matrices are the references to the original ones */
1750: PetscErrorCode PCGetInterpolations_MG(PC pc, PetscInt *num_levels, Mat *interpolations[])
1751: {
1752: PC_MG *mg = (PC_MG *)pc->data;
1753: PC_MG_Levels **mglevels = mg->levels;
1754: Mat *mat;
1755: PetscInt l;
1758: PetscMalloc1(mg->nlevels, &mat);
1759: for (l = 1; l < mg->nlevels; l++) {
1760: mat[l - 1] = mglevels[l]->interpolate;
1761: PetscObjectReference((PetscObject)mat[l - 1]);
1762: }
1763: *num_levels = mg->nlevels;
1764: *interpolations = mat;
1765: return 0;
1766: }
1768: /* No new matrices are created, and the coarse operator matrices are the references to the original ones */
1769: PetscErrorCode PCGetCoarseOperators_MG(PC pc, PetscInt *num_levels, Mat *coarseOperators[])
1770: {
1771: PC_MG *mg = (PC_MG *)pc->data;
1772: PC_MG_Levels **mglevels = mg->levels;
1773: PetscInt l;
1774: Mat *mat;
1777: PetscMalloc1(mg->nlevels, &mat);
1778: for (l = 0; l < mg->nlevels - 1; l++) {
1779: KSPGetOperators(mglevels[l]->smoothd, NULL, &(mat[l]));
1780: PetscObjectReference((PetscObject)mat[l]);
1781: }
1782: *num_levels = mg->nlevels;
1783: *coarseOperators = mat;
1784: return 0;
1785: }
1787: /*@C
1788: PCMGRegisterCoarseSpaceConstructor - Adds a method to the `PCMG` package for coarse space construction.
1790: Not collective
1792: Input Parameters:
1793: + name - name of the constructor
1794: - function - constructor routine
1796: Calling sequence for the routine:
1797: $ my_csp(PC pc, PetscInt l, DM dm, KSP smooth, PetscInt Nc, Mat initGuess, Mat *coarseSp)
1798: + pc - The PC object
1799: . l - The multigrid level, 0 is the coarse level
1800: . dm - The DM for this level
1801: . smooth - The level smoother
1802: . Nc - The size of the coarse space
1803: . initGuess - Basis for an initial guess for the space
1804: - coarseSp - A basis for the computed coarse space
1806: Level: advanced
1808: Developer Note:
1809: How come this is not used by `PCGAMG`?
1811: .seealso: `PCMG`, `PCMGGetCoarseSpaceConstructor()`, `PCRegister()`
1812: @*/
1813: PetscErrorCode PCMGRegisterCoarseSpaceConstructor(const char name[], PetscErrorCode (*function)(PC, PetscInt, DM, KSP, PetscInt, Mat, Mat *))
1814: {
1815: PCInitializePackage();
1816: PetscFunctionListAdd(&PCMGCoarseList, name, function);
1817: return 0;
1818: }
1820: /*@C
1821: PCMGGetCoarseSpaceConstructor - Returns the given coarse space construction method.
1823: Not collective
1825: Input Parameter:
1826: . name - name of the constructor
1828: Output Parameter:
1829: . function - constructor routine
1831: Level: advanced
1833: .seealso: `PCMG`, `PCMGRegisterCoarseSpaceConstructor()`, `PCRegister()`
1834: @*/
1835: PetscErrorCode PCMGGetCoarseSpaceConstructor(const char name[], PetscErrorCode (**function)(PC, PetscInt, DM, KSP, PetscInt, Mat, Mat *))
1836: {
1837: PetscFunctionListFind(PCMGCoarseList, name, function);
1838: return 0;
1839: }
1841: /*MC
1842: PCMG - Use multigrid preconditioning. This preconditioner requires you provide additional
1843: information about the coarser grid matrices and restriction/interpolation operators.
1845: Options Database Keys:
1846: + -pc_mg_levels <nlevels> - number of levels including finest
1847: . -pc_mg_cycle_type <v,w> - provide the cycle desired
1848: . -pc_mg_type <additive,multiplicative,full,kaskade> - multiplicative is the default
1849: . -pc_mg_log - log information about time spent on each level of the solver
1850: . -pc_mg_distinct_smoothup - configure up (after interpolation) and down (before restriction) smoothers separately (with different options prefixes)
1851: . -pc_mg_galerkin <both,pmat,mat,none> - use Galerkin process to compute coarser operators, i.e. Acoarse = R A R'
1852: . -pc_mg_multiplicative_cycles - number of cycles to use as the preconditioner (defaults to 1)
1853: . -pc_mg_dump_matlab - dumps the matrices for each level and the restriction/interpolation matrices
1854: to the Socket viewer for reading from MATLAB.
1855: - -pc_mg_dump_binary - dumps the matrices for each level and the restriction/interpolation matrices
1856: to the binary output file called binaryoutput
1858: Notes:
1859: If one uses a Krylov method such `KSPGMRES` or `KSPCG` as the smoother then one must use `KSPFGMRES`, `KSPGCR`, or `KSPRICHARDSON` as the outer Krylov method
1861: When run with a single level the smoother options are used on that level NOT the coarse grid solver options
1863: When run with `KSPRICHARDSON` the convergence test changes slightly if monitor is turned on. The iteration count may change slightly. This
1864: is because without monitoring the residual norm is computed WITHIN each multigrid cycle on the finest level after the pre-smoothing
1865: (because the residual has just been computed for the multigrid algorithm and is hence available for free) while with monitoring the
1866: residual is computed at the end of each cycle.
1868: Level: intermediate
1870: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCMGType`, `PCEXOTIC`, `PCGAMG`, `PCML`, `PCHYPRE`
1871: `PCMGSetLevels()`, `PCMGGetLevels()`, `PCMGSetType()`, `PCMGSetCycleType()`,
1872: `PCMGSetDistinctSmoothUp()`, `PCMGGetCoarseSolve()`, `PCMGSetResidual()`, `PCMGSetInterpolation()`,
1873: `PCMGSetRestriction()`, `PCMGGetSmoother()`, `PCMGGetSmootherUp()`, `PCMGGetSmootherDown()`,
1874: `PCMGSetCycleTypeOnLevel()`, `PCMGSetRhs()`, `PCMGSetX()`, `PCMGSetR()`,
1875: `PCMGSetAdaptCR()`, `PCMGGetAdaptInterpolation()`, `PCMGSetGalerkin()`, `PCMGGetAdaptCoarseSpaceType()`, `PCMGSetAdaptCoarseSpaceType()`
1876: M*/
1878: PETSC_EXTERN PetscErrorCode PCCreate_MG(PC pc)
1879: {
1880: PC_MG *mg;
1882: PetscNew(&mg);
1883: pc->data = mg;
1884: mg->nlevels = -1;
1885: mg->am = PC_MG_MULTIPLICATIVE;
1886: mg->galerkin = PC_MG_GALERKIN_NONE;
1887: mg->adaptInterpolation = PETSC_FALSE;
1888: mg->Nc = -1;
1889: mg->eigenvalue = -1;
1891: pc->useAmat = PETSC_TRUE;
1893: pc->ops->apply = PCApply_MG;
1894: pc->ops->applytranspose = PCApplyTranspose_MG;
1895: pc->ops->matapply = PCMatApply_MG;
1896: pc->ops->setup = PCSetUp_MG;
1897: pc->ops->reset = PCReset_MG;
1898: pc->ops->destroy = PCDestroy_MG;
1899: pc->ops->setfromoptions = PCSetFromOptions_MG;
1900: pc->ops->view = PCView_MG;
1902: PetscObjectComposedDataRegister(&mg->eigenvalue);
1903: PetscObjectComposeFunction((PetscObject)pc, "PCMGSetGalerkin_C", PCMGSetGalerkin_MG);
1904: PetscObjectComposeFunction((PetscObject)pc, "PCMGGetLevels_C", PCMGGetLevels_MG);
1905: PetscObjectComposeFunction((PetscObject)pc, "PCMGSetLevels_C", PCMGSetLevels_MG);
1906: PetscObjectComposeFunction((PetscObject)pc, "PCGetInterpolations_C", PCGetInterpolations_MG);
1907: PetscObjectComposeFunction((PetscObject)pc, "PCGetCoarseOperators_C", PCGetCoarseOperators_MG);
1908: PetscObjectComposeFunction((PetscObject)pc, "PCMGSetAdaptInterpolation_C", PCMGSetAdaptInterpolation_MG);
1909: PetscObjectComposeFunction((PetscObject)pc, "PCMGGetAdaptInterpolation_C", PCMGGetAdaptInterpolation_MG);
1910: PetscObjectComposeFunction((PetscObject)pc, "PCMGSetAdaptCR_C", PCMGSetAdaptCR_MG);
1911: PetscObjectComposeFunction((PetscObject)pc, "PCMGGetAdaptCR_C", PCMGGetAdaptCR_MG);
1912: PetscObjectComposeFunction((PetscObject)pc, "PCMGSetAdaptCoarseSpaceType_C", PCMGSetAdaptCoarseSpaceType_MG);
1913: PetscObjectComposeFunction((PetscObject)pc, "PCMGGetAdaptCoarseSpaceType_C", PCMGGetAdaptCoarseSpaceType_MG);
1914: return 0;
1915: }