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 *)&gtype, &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, &gtype);
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: }