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) {
 39:       MatDuplicate(mglevels->B,MAT_DO_NOT_COPY_VALUES,&mglevels->R);
 40:     }
 41:     if (!transpose) {
 42:       if (matapp) (*mglevels->matresidual)(mglevels->A,mglevels->B,mglevels->X,mglevels->R);
 43:       else (*mglevels->residual)(mglevels->A,mglevels->b,mglevels->x,mglevels->r);
 44:     } else {
 45:       if (matapp) (*mglevels->matresidualtranspose)(mglevels->A,mglevels->B,mglevels->X,mglevels->R);
 46:       else (*mglevels->residualtranspose)(mglevels->A,mglevels->b,mglevels->x,mglevels->r);
 47:     }
 48:     if (mglevels->eventresidual) PetscLogEventEnd(mglevels->eventresidual,0,0,0,0);

 50:     /* if on finest level and have convergence criteria set */
 51:     if (mglevels->level == mglevels->levels-1 && mg->ttol && reason) {
 52:       PetscReal rnorm;
 53:       VecNorm(mglevels->r,NORM_2,&rnorm);
 54:       if (rnorm <= mg->ttol) {
 55:         if (rnorm < mg->abstol) {
 56:           *reason = PCRICHARDSON_CONVERGED_ATOL;
 57:           PetscInfo(pc,"Linear solver has converged. Residual norm %g is less than absolute tolerance %g\n",(double)rnorm,(double)mg->abstol);
 58:         } else {
 59:           *reason = PCRICHARDSON_CONVERGED_RTOL;
 60:           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);
 61:         }
 62:         return 0;
 63:       }
 64:     }

 66:     mgc = *(mglevelsin - 1);
 67:     if (mglevels->eventinterprestrict) PetscLogEventBegin(mglevels->eventinterprestrict,0,0,0,0);
 68:     if (!transpose) {
 69:       if (matapp) MatMatRestrict(mglevels->restrct,mglevels->R,&mgc->B);
 70:       else MatRestrict(mglevels->restrct,mglevels->r,mgc->b);
 71:     } else {
 72:       if (matapp) MatMatRestrict(mglevels->interpolate,mglevels->R,&mgc->B);
 73:       else MatRestrict(mglevels->interpolate,mglevels->r,mgc->b);
 74:     }
 75:     if (mglevels->eventinterprestrict) PetscLogEventEnd(mglevels->eventinterprestrict,0,0,0,0);
 76:     if (matapp) {
 77:       if (!mgc->X) {
 78:         MatDuplicate(mgc->B,MAT_DO_NOT_COPY_VALUES,&mgc->X);
 79:       } else {
 80:         MatZeroEntries(mgc->X);
 81:       }
 82:     } else {
 83:       VecZeroEntries(mgc->x);
 84:     }
 85:     while (cycles--) {
 86:       PCMGMCycle_Private(pc,mglevelsin-1,transpose,matapp,reason);
 87:     }
 88:     if (mglevels->eventinterprestrict) PetscLogEventBegin(mglevels->eventinterprestrict,0,0,0,0);
 89:     if (!transpose) {
 90:       if (matapp) MatMatInterpolateAdd(mglevels->interpolate,mgc->X,mglevels->X,&mglevels->X);
 91:       else MatInterpolateAdd(mglevels->interpolate,mgc->x,mglevels->x,mglevels->x);
 92:     } else {
 93:       MatInterpolateAdd(mglevels->restrct,mgc->x,mglevels->x,mglevels->x);
 94:     }
 95:     if (mglevels->eventinterprestrict) PetscLogEventEnd(mglevels->eventinterprestrict,0,0,0,0);
 96:     if (mglevels->eventsmoothsolve) PetscLogEventBegin(mglevels->eventsmoothsolve,0,0,0,0);
 97:     if (!transpose) {
 98:       if (matapp) {
 99:         KSPMatSolve(mglevels->smoothu,mglevels->B,mglevels->X);    /* post smooth */
100:         KSPCheckSolve(mglevels->smoothu,pc,NULL);
101:       } else {
102:         KSPSolve(mglevels->smoothu,mglevels->b,mglevels->x);    /* post smooth */
103:         KSPCheckSolve(mglevels->smoothu,pc,mglevels->x);
104:       }
105:     } else {
107:       KSPSolveTranspose(mglevels->smoothd,mglevels->b,mglevels->x);    /* post smooth */
108:       KSPCheckSolve(mglevels->smoothd,pc,mglevels->x);
109:     }
110:     if (mglevels->cr) {
112:       /* TODO Turn on copy and turn off noisy if we have an exact solution
113:       VecCopy(mglevels->x, mglevels->crx);
114:       VecCopy(mglevels->b, mglevels->crb); */
115:       KSPSetNoisy_Private(mglevels->crx);
116:       KSPSolve(mglevels->cr,mglevels->crb,mglevels->crx);    /* compatible relaxation */
117:       KSPCheckSolve(mglevels->cr,pc,mglevels->crx);
118:     }
119:     if (mglevels->eventsmoothsolve) PetscLogEventEnd(mglevels->eventsmoothsolve,0,0,0,0);
120:   }
121:   return 0;
122: }

124: 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)
125: {
126:   PC_MG          *mg        = (PC_MG*)pc->data;
127:   PC_MG_Levels   **mglevels = mg->levels;
128:   PC             tpc;
129:   PetscBool      changeu,changed;
130:   PetscInt       levels = mglevels[0]->levels,i;

132:   /* When the DM is supplying the matrix then it will not exist until here */
133:   for (i=0; i<levels; i++) {
134:     if (!mglevels[i]->A) {
135:       KSPGetOperators(mglevels[i]->smoothu,&mglevels[i]->A,NULL);
136:       PetscObjectReference((PetscObject)mglevels[i]->A);
137:     }
138:   }

140:   KSPGetPC(mglevels[levels-1]->smoothd,&tpc);
141:   PCPreSolveChangeRHS(tpc,&changed);
142:   KSPGetPC(mglevels[levels-1]->smoothu,&tpc);
143:   PCPreSolveChangeRHS(tpc,&changeu);
144:   if (!changed && !changeu) {
145:     VecDestroy(&mglevels[levels-1]->b);
146:     mglevels[levels-1]->b = b;
147:   } else { /* if the smoother changes the rhs during PreSolve, we cannot use the input vector */
148:     if (!mglevels[levels-1]->b) {
149:       Vec *vec;

151:       KSPCreateVecs(mglevels[levels-1]->smoothd,1,&vec,0,NULL);
152:       mglevels[levels-1]->b = *vec;
153:       PetscFree(vec);
154:     }
155:     VecCopy(b,mglevels[levels-1]->b);
156:   }
157:   mglevels[levels-1]->x = x;

159:   mg->rtol   = rtol;
160:   mg->abstol = abstol;
161:   mg->dtol   = dtol;
162:   if (rtol) {
163:     /* compute initial residual norm for relative convergence test */
164:     PetscReal rnorm;
165:     if (zeroguess) {
166:       VecNorm(b,NORM_2,&rnorm);
167:     } else {
168:       (*mglevels[levels-1]->residual)(mglevels[levels-1]->A,b,x,w);
169:       VecNorm(w,NORM_2,&rnorm);
170:     }
171:     mg->ttol = PetscMax(rtol*rnorm,abstol);
172:   } else if (abstol) mg->ttol = abstol;
173:   else mg->ttol = 0.0;

175:   /* since smoother is applied to full system, not just residual we need to make sure that smoothers don't
176:      stop prematurely due to small residual */
177:   for (i=1; i<levels; i++) {
178:     KSPSetTolerances(mglevels[i]->smoothu,0,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
179:     if (mglevels[i]->smoothu != mglevels[i]->smoothd) {
180:       /* For Richardson the initial guess is nonzero since it is solving in each cycle the original system not just applying as a preconditioner */
181:       KSPSetInitialGuessNonzero(mglevels[i]->smoothd,PETSC_TRUE);
182:       KSPSetTolerances(mglevels[i]->smoothd,0,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
183:     }
184:   }

186:   *reason = (PCRichardsonConvergedReason)0;
187:   for (i=0; i<its; i++) {
188:     PCMGMCycle_Private(pc,mglevels+levels-1,PETSC_FALSE,PETSC_FALSE,reason);
189:     if (*reason) break;
190:   }
191:   if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS;
192:   *outits = i;
193:   if (!changed && !changeu) mglevels[levels-1]->b = NULL;
194:   return 0;
195: }

197: PetscErrorCode PCReset_MG(PC pc)
198: {
199:   PC_MG          *mg        = (PC_MG*)pc->data;
200:   PC_MG_Levels   **mglevels = mg->levels;
201:   PetscInt       i,c,n;

203:   if (mglevels) {
204:     n = mglevels[0]->levels;
205:     for (i=0; i<n-1; i++) {
206:       VecDestroy(&mglevels[i+1]->r);
207:       VecDestroy(&mglevels[i]->b);
208:       VecDestroy(&mglevels[i]->x);
209:       MatDestroy(&mglevels[i+1]->R);
210:       MatDestroy(&mglevels[i]->B);
211:       MatDestroy(&mglevels[i]->X);
212:       VecDestroy(&mglevels[i]->crx);
213:       VecDestroy(&mglevels[i]->crb);
214:       MatDestroy(&mglevels[i+1]->restrct);
215:       MatDestroy(&mglevels[i+1]->interpolate);
216:       MatDestroy(&mglevels[i+1]->inject);
217:       VecDestroy(&mglevels[i+1]->rscale);
218:     }
219:     VecDestroy(&mglevels[n-1]->crx);
220:     VecDestroy(&mglevels[n-1]->crb);
221:     /* this is not null only if the smoother on the finest level
222:        changes the rhs during PreSolve */
223:     VecDestroy(&mglevels[n-1]->b);
224:     MatDestroy(&mglevels[n-1]->B);

226:     for (i=0; i<n; i++) {
227:       if (mglevels[i]->coarseSpace) for (c = 0; c < mg->Nc; ++c) VecDestroy(&mglevels[i]->coarseSpace[c]);
228:       PetscFree(mglevels[i]->coarseSpace);
229:       mglevels[i]->coarseSpace = NULL;
230:       MatDestroy(&mglevels[i]->A);
231:       if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
232:         KSPReset(mglevels[i]->smoothd);
233:       }
234:       KSPReset(mglevels[i]->smoothu);
235:       if (mglevels[i]->cr) KSPReset(mglevels[i]->cr);
236:     }
237:     mg->Nc = 0;
238:   }
239:   return 0;
240: }

242: /* Implementing CR

244: 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

246:   Inj^T (Inj Inj^T)^{-1} Inj

248: and if Inj is a VecScatter, as it is now in PETSc, we have

250:   Inj^T Inj

252: and

254:   S = I - Inj^T Inj

256: since

258:   Inj S = Inj - (Inj Inj^T) Inj = 0.

260: Brannick suggests

262:   A \to S^T A S  \qquad\mathrm{and}\qquad M \to S^T M S

264: 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

266:   M^{-1} A \to S M^{-1} A S

268: In fact, since it is somewhat hard in PETSc to do the symmetric application, we will just apply S on the left.

270:   Check: || Inj P - I ||_F < tol
271:   Check: In general, Inj Inj^T = I
272: */

274: typedef struct {
275:   PC       mg;  /* The PCMG object */
276:   PetscInt l;   /* The multigrid level for this solver */
277:   Mat      Inj; /* The injection matrix */
278:   Mat      S;   /* I - Inj^T Inj */
279: } CRContext;

281: static PetscErrorCode CRSetup_Private(PC pc)
282: {
283:   CRContext     *ctx;
284:   Mat            It;

287:   PCShellGetContext(pc, &ctx);
288:   PCMGGetInjection(ctx->mg, ctx->l, &It);
290:   MatCreateTranspose(It, &ctx->Inj);
291:   MatCreateNormal(ctx->Inj, &ctx->S);
292:   MatScale(ctx->S, -1.0);
293:   MatShift(ctx->S,  1.0);
294:   return 0;
295: }

297: static PetscErrorCode CRApply_Private(PC pc, Vec x, Vec y)
298: {
299:   CRContext     *ctx;

302:   PCShellGetContext(pc, &ctx);
303:   MatMult(ctx->S, x, y);
304:   return 0;
305: }

307: static PetscErrorCode CRDestroy_Private(PC pc)
308: {
309:   CRContext     *ctx;

312:   PCShellGetContext(pc, &ctx);
313:   MatDestroy(&ctx->Inj);
314:   MatDestroy(&ctx->S);
315:   PetscFree(ctx);
316:   PCShellSetContext(pc, NULL);
317:   return 0;
318: }

320: static PetscErrorCode CreateCR_Private(PC pc, PetscInt l, PC *cr)
321: {
322:   CRContext     *ctx;

325:   PCCreate(PetscObjectComm((PetscObject) pc), cr);
326:   PetscObjectSetName((PetscObject) *cr, "S (complementary projector to injection)");
327:   PetscCalloc1(1, &ctx);
328:   ctx->mg = pc;
329:   ctx->l  = l;
330:   PCSetType(*cr, PCSHELL);
331:   PCShellSetContext(*cr, ctx);
332:   PCShellSetApply(*cr, CRApply_Private);
333:   PCShellSetSetUp(*cr, CRSetup_Private);
334:   PCShellSetDestroy(*cr, CRDestroy_Private);
335:   return 0;
336: }

338: PetscErrorCode PCMGSetLevels_MG(PC pc,PetscInt levels,MPI_Comm *comms)
339: {
340:   PC_MG          *mg        = (PC_MG*)pc->data;
341:   MPI_Comm       comm;
342:   PC_MG_Levels   **mglevels = mg->levels;
343:   PCMGType       mgtype     = mg->am;
344:   PetscInt       mgctype    = (PetscInt) PC_MG_CYCLE_V;
345:   PetscInt       i;
346:   PetscMPIInt    size;
347:   const char     *prefix;
348:   PC             ipc;
349:   PetscInt       n;

353:   if (mg->nlevels == levels) return 0;
354:   PetscObjectGetComm((PetscObject)pc,&comm);
355:   if (mglevels) {
356:     mgctype = mglevels[0]->cycles;
357:     /* changing the number of levels so free up the previous stuff */
358:     PCReset_MG(pc);
359:     n    = mglevels[0]->levels;
360:     for (i=0; i<n; i++) {
361:       if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
362:         KSPDestroy(&mglevels[i]->smoothd);
363:       }
364:       KSPDestroy(&mglevels[i]->smoothu);
365:       KSPDestroy(&mglevels[i]->cr);
366:       PetscFree(mglevels[i]);
367:     }
368:     PetscFree(mg->levels);
369:   }

371:   mg->nlevels = levels;

373:   PetscMalloc1(levels,&mglevels);
374:   PetscLogObjectMemory((PetscObject)pc,levels*(sizeof(PC_MG*)));

376:   PCGetOptionsPrefix(pc,&prefix);

378:   mg->stageApply = 0;
379:   for (i=0; i<levels; i++) {
380:     PetscNewLog(pc,&mglevels[i]);

382:     mglevels[i]->level               = i;
383:     mglevels[i]->levels              = levels;
384:     mglevels[i]->cycles              = mgctype;
385:     mg->default_smoothu              = 2;
386:     mg->default_smoothd              = 2;
387:     mglevels[i]->eventsmoothsetup    = 0;
388:     mglevels[i]->eventsmoothsolve    = 0;
389:     mglevels[i]->eventresidual       = 0;
390:     mglevels[i]->eventinterprestrict = 0;

392:     if (comms) comm = comms[i];
393:     if (comm != MPI_COMM_NULL) {
394:       KSPCreate(comm,&mglevels[i]->smoothd);
395:       KSPSetErrorIfNotConverged(mglevels[i]->smoothd,pc->erroriffailure);
396:       PetscObjectIncrementTabLevel((PetscObject)mglevels[i]->smoothd,(PetscObject)pc,levels-i);
397:       KSPSetOptionsPrefix(mglevels[i]->smoothd,prefix);
398:       PetscObjectComposedDataSetInt((PetscObject) mglevels[i]->smoothd, PetscMGLevelId, mglevels[i]->level);
399:       if (i || levels == 1) {
400:         char tprefix[128];

402:         KSPSetType(mglevels[i]->smoothd,KSPCHEBYSHEV);
403:         KSPSetConvergenceTest(mglevels[i]->smoothd,KSPConvergedSkip,NULL,NULL);
404:         KSPSetNormType(mglevels[i]->smoothd,KSP_NORM_NONE);
405:         KSPGetPC(mglevels[i]->smoothd,&ipc);
406:         PCSetType(ipc,PCSOR);
407:         KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT, mg->default_smoothd);

409:         PetscSNPrintf(tprefix,128,"mg_levels_%d_",(int)i);
410:         KSPAppendOptionsPrefix(mglevels[i]->smoothd,tprefix);
411:       } else {
412:         KSPAppendOptionsPrefix(mglevels[0]->smoothd,"mg_coarse_");

414:         /* coarse solve is (redundant) LU by default; set shifttype NONZERO to avoid annoying zero-pivot in LU preconditioner */
415:         KSPSetType(mglevels[0]->smoothd,KSPPREONLY);
416:         KSPGetPC(mglevels[0]->smoothd,&ipc);
417:         MPI_Comm_size(comm,&size);
418:         if (size > 1) {
419:           PCSetType(ipc,PCREDUNDANT);
420:         } else {
421:           PCSetType(ipc,PCLU);
422:         }
423:         PCFactorSetShiftType(ipc,MAT_SHIFT_INBLOCKS);
424:       }
425:       PetscLogObjectParent((PetscObject)pc,(PetscObject)mglevels[i]->smoothd);
426:     }
427:     mglevels[i]->smoothu = mglevels[i]->smoothd;
428:     mg->rtol             = 0.0;
429:     mg->abstol           = 0.0;
430:     mg->dtol             = 0.0;
431:     mg->ttol             = 0.0;
432:     mg->cyclesperpcapply = 1;
433:   }
434:   mg->levels = mglevels;
435:   PCMGSetType(pc,mgtype);
436:   return 0;
437: }

439: /*@C
440:    PCMGSetLevels - Sets the number of levels to use with MG.
441:    Must be called before any other MG routine.

443:    Logically Collective on PC

445:    Input Parameters:
446: +  pc - the preconditioner context
447: .  levels - the number of levels
448: -  comms - optional communicators for each level; this is to allow solving the coarser problems
449:            on smaller sets of processes. For processes that are not included in the computation
450:            you must pass MPI_COMM_NULL. Use comms = NULL to specify that all processes
451:            should participate in each level of problem.

453:    Level: intermediate

455:    Notes:
456:      If the number of levels is one then the multigrid uses the -mg_levels prefix
457:      for setting the level options rather than the -mg_coarse prefix.

459:      You can free the information in comms after this routine is called.

461:      The array of MPI communicators must contain MPI_COMM_NULL for those ranks that at each level
462:      are not participating in the coarser solve. For example, with 2 levels and 1 and 2 ranks on
463:      the two levels, rank 0 in the original communicator will pass in an array of 2 communicators
464:      of size 2 and 1, while rank 1 in the original communicator will pass in array of 2 communicators
465:      the first of size 2 and the second of value MPI_COMM_NULL since the rank 1 does not participate
466:      in the coarse grid solve.

468:      Since each coarser level may have a new MPI_Comm with fewer ranks than the previous, one
469:      must take special care in providing the restriction and interpolation operation. We recommend
470:      providing these as two step operations; first perform a standard restriction or interpolation on
471:      the full number of ranks for that level and then use an MPI call to copy the resulting vector
472:      array entries (after calls to VecGetArray()) to the smaller or larger number of ranks, note in both
473:      cases the MPI calls must be made on the larger of the two communicators. Traditional MPI send and
474:      recieves or MPI_AlltoAllv() could be used to do the reshuffling of the vector entries.

476:    Fortran Notes:
477:      Use comms = PETSC_NULL_MPI_COMM as the equivalent of NULL in the C interface. Note PETSC_NULL_MPI_COMM
478:      is not MPI_COMM_NULL. It is more like PETSC_NULL_INTEGER, PETSC_NULL_REAL etc.

480: .seealso: PCMGSetType(), PCMGGetLevels()
481: @*/
482: PetscErrorCode PCMGSetLevels(PC pc,PetscInt levels,MPI_Comm *comms)
483: {
486:   PetscTryMethod(pc,"PCMGSetLevels_C",(PC,PetscInt,MPI_Comm*),(pc,levels,comms));
487:   return 0;
488: }

490: PetscErrorCode PCDestroy_MG(PC pc)
491: {
492:   PC_MG          *mg        = (PC_MG*)pc->data;
493:   PC_MG_Levels   **mglevels = mg->levels;
494:   PetscInt       i,n;

496:   PCReset_MG(pc);
497:   if (mglevels) {
498:     n = mglevels[0]->levels;
499:     for (i=0; i<n; i++) {
500:       if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
501:         KSPDestroy(&mglevels[i]->smoothd);
502:       }
503:       KSPDestroy(&mglevels[i]->smoothu);
504:       KSPDestroy(&mglevels[i]->cr);
505:       PetscFree(mglevels[i]);
506:     }
507:     PetscFree(mg->levels);
508:   }
509:   PetscFree(pc->data);
510:   PetscObjectComposeFunction((PetscObject)pc,"PCGetInterpolations_C",NULL);
511:   PetscObjectComposeFunction((PetscObject)pc,"PCGetCoarseOperators_C",NULL);
512:   return 0;
513: }

515: /*
516:    PCApply_MG - Runs either an additive, multiplicative, Kaskadic
517:              or full cycle of multigrid.

519:   Note:
520:   A simple wrapper which calls PCMGMCycle(),PCMGACycle(), or PCMGFCycle().
521: */
522: static PetscErrorCode PCApply_MG_Internal(PC pc,Vec b,Vec x,Mat B,Mat X,PetscBool transpose)
523: {
524:   PC_MG          *mg        = (PC_MG*)pc->data;
525:   PC_MG_Levels   **mglevels = mg->levels;
526:   PC             tpc;
527:   PetscInt       levels = mglevels[0]->levels,i;
528:   PetscBool      changeu,changed,matapp;

530:   matapp = (PetscBool)(B && X);
531:   if (mg->stageApply) PetscLogStagePush(mg->stageApply);
532:   /* When the DM is supplying the matrix then it will not exist until here */
533:   for (i=0; i<levels; i++) {
534:     if (!mglevels[i]->A) {
535:       KSPGetOperators(mglevels[i]->smoothu,&mglevels[i]->A,NULL);
536:       PetscObjectReference((PetscObject)mglevels[i]->A);
537:     }
538:   }

540:   KSPGetPC(mglevels[levels-1]->smoothd,&tpc);
541:   PCPreSolveChangeRHS(tpc,&changed);
542:   KSPGetPC(mglevels[levels-1]->smoothu,&tpc);
543:   PCPreSolveChangeRHS(tpc,&changeu);
544:   if (!changeu && !changed) {
545:     if (matapp) {
546:       MatDestroy(&mglevels[levels-1]->B);
547:       mglevels[levels-1]->B = B;
548:     } else {
549:       VecDestroy(&mglevels[levels-1]->b);
550:       mglevels[levels-1]->b = b;
551:     }
552:   } else { /* if the smoother changes the rhs during PreSolve, we cannot use the input vector */
553:     if (matapp) {
554:       if (mglevels[levels-1]->B) {
555:         PetscInt  N1,N2;
556:         PetscBool flg;

558:         MatGetSize(mglevels[levels-1]->B,NULL,&N1);
559:         MatGetSize(B,NULL,&N2);
560:         PetscObjectTypeCompare((PetscObject)mglevels[levels-1]->B,((PetscObject)B)->type_name,&flg);
561:         if (N1 != N2 || !flg) {
562:           MatDestroy(&mglevels[levels-1]->B);
563:         }
564:       }
565:       if (!mglevels[levels-1]->B) {
566:         MatDuplicate(B,MAT_COPY_VALUES,&mglevels[levels-1]->B);
567:       } else {
568:         MatCopy(B,mglevels[levels-1]->B,SAME_NONZERO_PATTERN);
569:       }
570:     } else {
571:       if (!mglevels[levels-1]->b) {
572:         Vec *vec;

574:         KSPCreateVecs(mglevels[levels-1]->smoothd,1,&vec,0,NULL);
575:         mglevels[levels-1]->b = *vec;
576:         PetscFree(vec);
577:       }
578:       VecCopy(b,mglevels[levels-1]->b);
579:     }
580:   }
581:   if (matapp) { mglevels[levels-1]->X = X; }
582:   else { mglevels[levels-1]->x = x; }

584:   /* If coarser Xs are present, it means we have already block applied the PC at least once
585:      Reset operators if sizes/type do no match */
586:   if (matapp && levels > 1 && mglevels[levels-2]->X) {
587:     PetscInt  Xc,Bc;
588:     PetscBool flg;

590:     MatGetSize(mglevels[levels-2]->X,NULL,&Xc);
591:     MatGetSize(mglevels[levels-1]->B,NULL,&Bc);
592:     PetscObjectTypeCompare((PetscObject)mglevels[levels-2]->X,((PetscObject)mglevels[levels-1]->X)->type_name,&flg);
593:     if (Xc != Bc || !flg) {
594:       MatDestroy(&mglevels[levels-1]->R);
595:       for (i=0;i<levels-1;i++) {
596:         MatDestroy(&mglevels[i]->R);
597:         MatDestroy(&mglevels[i]->B);
598:         MatDestroy(&mglevels[i]->X);
599:       }
600:     }
601:   }

603:   if (mg->am == PC_MG_MULTIPLICATIVE) {
604:     if (matapp) MatZeroEntries(X);
605:     else VecZeroEntries(x);
606:     for (i=0; i<mg->cyclesperpcapply; i++) {
607:       PCMGMCycle_Private(pc,mglevels+levels-1,transpose,matapp,NULL);
608:     }
609:   } else if (mg->am == PC_MG_ADDITIVE) {
610:     PCMGACycle_Private(pc,mglevels,transpose,matapp);
611:   } else if (mg->am == PC_MG_KASKADE) {
612:     PCMGKCycle_Private(pc,mglevels,transpose,matapp);
613:   } else {
614:     PCMGFCycle_Private(pc,mglevels,transpose,matapp);
615:   }
616:   if (mg->stageApply) PetscLogStagePop();
617:   if (!changeu && !changed) {
618:     if (matapp) { mglevels[levels-1]->B = NULL; }
619:     else { mglevels[levels-1]->b = NULL; }
620:   }
621:   return 0;
622: }

624: static PetscErrorCode PCApply_MG(PC pc,Vec b,Vec x)
625: {
626:   PCApply_MG_Internal(pc,b,x,NULL,NULL,PETSC_FALSE);
627:   return 0;
628: }

630: static PetscErrorCode PCApplyTranspose_MG(PC pc,Vec b,Vec x)
631: {
632:   PCApply_MG_Internal(pc,b,x,NULL,NULL,PETSC_TRUE);
633:   return 0;
634: }

636: static PetscErrorCode PCMatApply_MG(PC pc,Mat b,Mat x)
637: {
638:   PCApply_MG_Internal(pc,NULL,NULL,b,x,PETSC_FALSE);
639:   return 0;
640: }

642: PetscErrorCode PCSetFromOptions_MG(PetscOptionItems *PetscOptionsObject,PC pc)
643: {
644:   PetscInt         levels,cycles;
645:   PetscBool        flg, flg2;
646:   PC_MG            *mg = (PC_MG*)pc->data;
647:   PC_MG_Levels     **mglevels;
648:   PCMGType         mgtype;
649:   PCMGCycleType    mgctype;
650:   PCMGGalerkinType gtype;

652:   levels = PetscMax(mg->nlevels,1);
653:   PetscOptionsHead(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) {
666:     PCMGSetCycleType(pc,mgctype);
667:   }
668:   gtype = mg->galerkin;
669:   PetscOptionsEnum("-pc_mg_galerkin","Use Galerkin process to compute coarser operators","PCMGSetGalerkin",PCMGGalerkinTypes,(PetscEnum)gtype,(PetscEnum*)&gtype,&flg);
670:   if (flg) {
671:     PCMGSetGalerkin(pc,gtype);
672:   }
673:   flg2 = PETSC_FALSE;
674:   PetscOptionsBool("-pc_mg_adapt_interp","Adapt interpolation using some coarse space","PCMGSetAdaptInterpolation",PETSC_FALSE,&flg2,&flg);
675:   if (flg) PCMGSetAdaptInterpolation(pc, flg2);
676:   PetscOptionsInt("-pc_mg_adapt_interp_n","Size of the coarse space for adaptive interpolation","PCMGSetCoarseSpace",mg->Nc,&mg->Nc,&flg);
677:   PetscOptionsEnum("-pc_mg_adapt_interp_coarse_space","Type of coarse space: polynomial, harmonic, eigenvector, generalized_eigenvector","PCMGSetAdaptCoarseSpaceType",PCMGCoarseSpaceTypes,(PetscEnum)mg->coarseSpaceType,(PetscEnum*)&mg->coarseSpaceType,&flg);
678:   PetscOptionsBool("-pc_mg_mesp_monitor","Monitor the multilevel eigensolver","PCMGSetAdaptInterpolation",PETSC_FALSE,&mg->mespMonitor,&flg);
679:   flg2 = PETSC_FALSE;
680:   PetscOptionsBool("-pc_mg_adapt_cr","Monitor coarse space quality using Compatible Relaxation (CR)","PCMGSetAdaptCR",PETSC_FALSE,&flg2,&flg);
681:   if (flg) PCMGSetAdaptCR(pc, flg2);
682:   flg = PETSC_FALSE;
683:   PetscOptionsBool("-pc_mg_distinct_smoothup","Create separate smoothup KSP and append the prefix _up","PCMGSetDistinctSmoothUp",PETSC_FALSE,&flg,NULL);
684:   if (flg) {
685:     PCMGSetDistinctSmoothUp(pc);
686:   }
687:   mgtype = mg->am;
688:   PetscOptionsEnum("-pc_mg_type","Multigrid type","PCMGSetType",PCMGTypes,(PetscEnum)mgtype,(PetscEnum*)&mgtype,&flg);
689:   if (flg) {
690:     PCMGSetType(pc,mgtype);
691:   }
692:   if (mg->am == PC_MG_MULTIPLICATIVE) {
693:     PetscOptionsInt("-pc_mg_multiplicative_cycles","Number of cycles for each preconditioner step","PCMGMultiplicativeSetCycles",mg->cyclesperpcapply,&cycles,&flg);
694:     if (flg) {
695:       PCMGMultiplicativeSetCycles(pc,cycles);
696:     }
697:   }
698:   flg  = PETSC_FALSE;
699:   PetscOptionsBool("-pc_mg_log","Log times for each multigrid level","None",flg,&flg,NULL);
700:   if (flg) {
701:     PetscInt i;
702:     char     eventname[128];

704:     levels = mglevels[0]->levels;
705:     for (i=0; i<levels; i++) {
706:       sprintf(eventname,"MGSetup Level %d",(int)i);
707:       PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventsmoothsetup);
708:       sprintf(eventname,"MGSmooth Level %d",(int)i);
709:       PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventsmoothsolve);
710:       if (i) {
711:         sprintf(eventname,"MGResid Level %d",(int)i);
712:         PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventresidual);
713:         sprintf(eventname,"MGInterp Level %d",(int)i);
714:         PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventinterprestrict);
715:       }
716:     }

718: #if defined(PETSC_USE_LOG)
719:     {
720:       const char    *sname = "MG Apply";
721:       PetscStageLog stageLog;
722:       PetscInt      st;

724:       PetscLogGetStageLog(&stageLog);
725:       for (st = 0; st < stageLog->numStages; ++st) {
726:         PetscBool same;

728:         PetscStrcmp(stageLog->stageInfo[st].name, sname, &same);
729:         if (same) mg->stageApply = st;
730:       }
731:       if (!mg->stageApply) {
732:         PetscLogStageRegister(sname, &mg->stageApply);
733:       }
734:     }
735: #endif
736:   }
737:   PetscOptionsTail();
738:   /* Check option consistency */
739:   PCMGGetGalerkin(pc, &gtype);
740:   PCMGGetAdaptInterpolation(pc, &flg);
742:   return 0;
743: }

745: const char *const PCMGTypes[] = {"MULTIPLICATIVE","ADDITIVE","FULL","KASKADE","PCMGType","PC_MG",NULL};
746: const char *const PCMGCycleTypes[] = {"invalid","v","w","PCMGCycleType","PC_MG_CYCLE",NULL};
747: const char *const PCMGGalerkinTypes[] = {"both","pmat","mat","none","external","PCMGGalerkinType","PC_MG_GALERKIN",NULL};
748: const char *const PCMGCoarseSpaceTypes[] = {"polynomial","harmonic","eigenvector","generalized_eigenvector","PCMGCoarseSpaceType","PCMG_POLYNOMIAL",NULL};

750: #include <petscdraw.h>
751: PetscErrorCode PCView_MG(PC pc,PetscViewer viewer)
752: {
753:   PC_MG          *mg        = (PC_MG*)pc->data;
754:   PC_MG_Levels   **mglevels = mg->levels;
755:   PetscInt       levels = mglevels ? mglevels[0]->levels : 0,i;
756:   PetscBool      iascii,isbinary,isdraw;

758:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
759:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
760:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
761:   if (iascii) {
762:     const char *cyclename = levels ? (mglevels[0]->cycles == PC_MG_CYCLE_V ? "v" : "w") : "unknown";
763:     PetscViewerASCIIPrintf(viewer,"  type is %s, levels=%D cycles=%s\n", PCMGTypes[mg->am],levels,cyclename);
764:     if (mg->am == PC_MG_MULTIPLICATIVE) {
765:       PetscViewerASCIIPrintf(viewer,"    Cycles per PCApply=%d\n",mg->cyclesperpcapply);
766:     }
767:     if (mg->galerkin == PC_MG_GALERKIN_BOTH) {
768:       PetscViewerASCIIPrintf(viewer,"    Using Galerkin computed coarse grid matrices\n");
769:     } else if (mg->galerkin == PC_MG_GALERKIN_PMAT) {
770:       PetscViewerASCIIPrintf(viewer,"    Using Galerkin computed coarse grid matrices for pmat\n");
771:     } else if (mg->galerkin == PC_MG_GALERKIN_MAT) {
772:       PetscViewerASCIIPrintf(viewer,"    Using Galerkin computed coarse grid matrices for mat\n");
773:     } else if (mg->galerkin == PC_MG_GALERKIN_EXTERNAL) {
774:       PetscViewerASCIIPrintf(viewer,"    Using externally compute Galerkin coarse grid matrices\n");
775:     } else {
776:       PetscViewerASCIIPrintf(viewer,"    Not using Galerkin computed coarse grid matrices\n");
777:     }
778:     if (mg->view) {
779:       (*mg->view)(pc,viewer);
780:     }
781:     for (i=0; i<levels; i++) {
782:       if (!i) {
783:         PetscViewerASCIIPrintf(viewer,"Coarse grid solver -- level -------------------------------\n",i);
784:       } else {
785:         PetscViewerASCIIPrintf(viewer,"Down solver (pre-smoother) on level %D -------------------------------\n",i);
786:       }
787:       PetscViewerASCIIPushTab(viewer);
788:       KSPView(mglevels[i]->smoothd,viewer);
789:       PetscViewerASCIIPopTab(viewer);
790:       if (i && mglevels[i]->smoothd == mglevels[i]->smoothu) {
791:         PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) same as down solver (pre-smoother)\n");
792:       } else if (i) {
793:         PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) on level %D -------------------------------\n",i);
794:         PetscViewerASCIIPushTab(viewer);
795:         KSPView(mglevels[i]->smoothu,viewer);
796:         PetscViewerASCIIPopTab(viewer);
797:       }
798:       if (i && mglevels[i]->cr) {
799:         PetscViewerASCIIPrintf(viewer,"CR solver on level %D -------------------------------\n",i);
800:         PetscViewerASCIIPushTab(viewer);
801:         KSPView(mglevels[i]->cr,viewer);
802:         PetscViewerASCIIPopTab(viewer);
803:       }
804:     }
805:   } else if (isbinary) {
806:     for (i=levels-1; i>=0; i--) {
807:       KSPView(mglevels[i]->smoothd,viewer);
808:       if (i && mglevels[i]->smoothd != mglevels[i]->smoothu) {
809:         KSPView(mglevels[i]->smoothu,viewer);
810:       }
811:     }
812:   } else if (isdraw) {
813:     PetscDraw draw;
814:     PetscReal x,w,y,bottom,th;
815:     PetscViewerDrawGetDraw(viewer,0,&draw);
816:     PetscDrawGetCurrentPoint(draw,&x,&y);
817:     PetscDrawStringGetSize(draw,NULL,&th);
818:     bottom = y - th;
819:     for (i=levels-1; i>=0; i--) {
820:       if (!mglevels[i]->smoothu || (mglevels[i]->smoothu == mglevels[i]->smoothd)) {
821:         PetscDrawPushCurrentPoint(draw,x,bottom);
822:         KSPView(mglevels[i]->smoothd,viewer);
823:         PetscDrawPopCurrentPoint(draw);
824:       } else {
825:         w    = 0.5*PetscMin(1.0-x,x);
826:         PetscDrawPushCurrentPoint(draw,x+w,bottom);
827:         KSPView(mglevels[i]->smoothd,viewer);
828:         PetscDrawPopCurrentPoint(draw);
829:         PetscDrawPushCurrentPoint(draw,x-w,bottom);
830:         KSPView(mglevels[i]->smoothu,viewer);
831:         PetscDrawPopCurrentPoint(draw);
832:       }
833:       PetscDrawGetBoundingBox(draw,NULL,&bottom,NULL,NULL);
834:       bottom -= th;
835:     }
836:   }
837:   return 0;
838: }

840: #include <petsc/private/kspimpl.h>

842: /*
843:     Calls setup for the KSP on each level
844: */
845: PetscErrorCode PCSetUp_MG(PC pc)
846: {
847:   PC_MG          *mg        = (PC_MG*)pc->data;
848:   PC_MG_Levels   **mglevels = mg->levels;
849:   PetscInt       i,n;
850:   PC             cpc;
851:   PetscBool      dump = PETSC_FALSE,opsset,use_amat,missinginterpolate = PETSC_FALSE;
852:   Mat            dA,dB;
853:   Vec            tvec;
854:   DM             *dms;
855:   PetscViewer    viewer = NULL;
856:   PetscBool      dAeqdB = PETSC_FALSE, needRestricts = PETSC_FALSE, doCR = PETSC_FALSE;

859:   n = mglevels[0]->levels;
860:   /* FIX: Move this to PCSetFromOptions_MG? */
861:   if (mg->usedmfornumberoflevels) {
862:     PetscInt levels;
863:     DMGetRefineLevel(pc->dm,&levels);
864:     levels++;
865:     if (levels > n) { /* the problem is now being solved on a finer grid */
866:       PCMGSetLevels(pc,levels,NULL);
867:       n        = levels;
868:       PCSetFromOptions(pc); /* it is bad to call this here, but otherwise will never be called for the new hierarchy */
869:       mglevels = mg->levels;
870:     }
871:   }
872:   KSPGetPC(mglevels[0]->smoothd,&cpc);

874:   /* If user did not provide fine grid operators OR operator was not updated since last global KSPSetOperators() */
875:   /* so use those from global PC */
876:   /* Is this what we always want? What if user wants to keep old one? */
877:   KSPGetOperatorsSet(mglevels[n-1]->smoothd,NULL,&opsset);
878:   if (opsset) {
879:     Mat mmat;
880:     KSPGetOperators(mglevels[n-1]->smoothd,NULL,&mmat);
881:     if (mmat == pc->pmat) opsset = PETSC_FALSE;
882:   }

884:   /* Create CR solvers */
885:   PCMGGetAdaptCR(pc, &doCR);
886:   if (doCR) {
887:     const char *prefix;

889:     PCGetOptionsPrefix(pc, &prefix);
890:     for (i = 1; i < n; ++i) {
891:       PC   ipc, cr;
892:       char crprefix[128];

894:       KSPCreate(PetscObjectComm((PetscObject) pc), &mglevels[i]->cr);
895:       KSPSetErrorIfNotConverged(mglevels[i]->cr, PETSC_FALSE);
896:       PetscObjectIncrementTabLevel((PetscObject) mglevels[i]->cr, (PetscObject) pc, n-i);
897:       KSPSetOptionsPrefix(mglevels[i]->cr, prefix);
898:       PetscObjectComposedDataSetInt((PetscObject) mglevels[i]->cr, PetscMGLevelId, mglevels[i]->level);
899:       KSPSetType(mglevels[i]->cr, KSPCHEBYSHEV);
900:       KSPSetConvergenceTest(mglevels[i]->cr, KSPConvergedSkip, NULL, NULL);
901:       KSPSetNormType(mglevels[i]->cr, KSP_NORM_PRECONDITIONED);
902:       KSPGetPC(mglevels[i]->cr, &ipc);

904:       PCSetType(ipc, PCCOMPOSITE);
905:       PCCompositeSetType(ipc, PC_COMPOSITE_MULTIPLICATIVE);
906:       PCCompositeAddPCType(ipc, PCSOR);
907:       CreateCR_Private(pc, i, &cr);
908:       PCCompositeAddPC(ipc, cr);
909:       PCDestroy(&cr);

911:       KSPSetTolerances(mglevels[i]->cr, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, mg->default_smoothd);
912:       KSPSetInitialGuessNonzero(mglevels[i]->cr, PETSC_TRUE);
913:       PetscSNPrintf(crprefix, 128, "mg_levels_%d_cr_", (int) i);
914:       KSPAppendOptionsPrefix(mglevels[i]->cr, crprefix);
915:       PetscLogObjectParent((PetscObject) pc, (PetscObject) mglevels[i]->cr);
916:     }
917:   }

919:   if (!opsset) {
920:     PCGetUseAmat(pc,&use_amat);
921:     if (use_amat) {
922:       PetscInfo(pc,"Using outer operators to define finest grid operator \n  because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");
923:       KSPSetOperators(mglevels[n-1]->smoothd,pc->mat,pc->pmat);
924:     } else {
925:       PetscInfo(pc,"Using matrix (pmat) operators to define finest grid operator \n  because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");
926:       KSPSetOperators(mglevels[n-1]->smoothd,pc->pmat,pc->pmat);
927:     }
928:   }

930:   for (i=n-1; i>0; i--) {
931:     if (!(mglevels[i]->interpolate || mglevels[i]->restrct)) {
932:       missinginterpolate = PETSC_TRUE;
933:       continue;
934:     }
935:   }

937:   KSPGetOperators(mglevels[n-1]->smoothd,&dA,&dB);
938:   if (dA == dB) dAeqdB = PETSC_TRUE;
939:   if ((mg->galerkin == PC_MG_GALERKIN_NONE) || (((mg->galerkin == PC_MG_GALERKIN_PMAT) || (mg->galerkin == PC_MG_GALERKIN_MAT)) && !dAeqdB)) {
940:     needRestricts = PETSC_TRUE;  /* user must compute either mat, pmat, or both so must restrict x to coarser levels */
941:   }

943:   /*
944:    Skipping if user has provided all interpolation/restriction needed (since DM might not be able to produce them (when coming from SNES/TS)
945:    Skipping for galerkin==2 (externally managed hierarchy such as ML and GAMG). Cleaner logic here would be great. Wrap ML/GAMG as DMs?
946:   */
947:   if (missinginterpolate && pc->dm && mg->galerkin != PC_MG_GALERKIN_EXTERNAL && !pc->setupcalled) {
948:         /* construct the interpolation from the DMs */
949:     Mat p;
950:     Vec rscale;
951:     PetscMalloc1(n,&dms);
952:     dms[n-1] = pc->dm;
953:     /* Separately create them so we do not get DMKSP interference between levels */
954:     for (i=n-2; i>-1; i--) DMCoarsen(dms[i+1],MPI_COMM_NULL,&dms[i]);
955:         /*
956:            Force the mat type of coarse level operator to be AIJ because usually we want to use LU for coarse level.
957:            Notice that it can be overwritten by -mat_type because KSPSetUp() reads command line options.
958:            But it is safe to use -dm_mat_type.

960:            The mat type should not be hardcoded like this, we need to find a better way.
961:     DMSetMatType(dms[0],MATAIJ);
962:     */
963:     for (i=n-2; i>-1; i--) {
964:       DMKSP     kdm;
965:       PetscBool dmhasrestrict, dmhasinject;
966:       KSPSetDM(mglevels[i]->smoothd,dms[i]);
967:       if (!needRestricts) KSPSetDMActive(mglevels[i]->smoothd,PETSC_FALSE);
968:       if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
969:         KSPSetDM(mglevels[i]->smoothu,dms[i]);
970:         if (!needRestricts) KSPSetDMActive(mglevels[i]->smoothu,PETSC_FALSE);
971:       }
972:       if (mglevels[i]->cr) {
973:         KSPSetDM(mglevels[i]->cr,dms[i]);
974:         if (!needRestricts) KSPSetDMActive(mglevels[i]->cr,PETSC_FALSE);
975:       }
976:       DMGetDMKSPWrite(dms[i],&kdm);
977:       /* Ugly hack so that the next KSPSetUp() will use the RHS that we set. A better fix is to change dmActive to take
978:        * a bitwise OR of computing the matrix, RHS, and initial iterate. */
979:       kdm->ops->computerhs = NULL;
980:       kdm->rhsctx          = NULL;
981:       if (!mglevels[i+1]->interpolate) {
982:         DMCreateInterpolation(dms[i],dms[i+1],&p,&rscale);
983:         PCMGSetInterpolation(pc,i+1,p);
984:         if (rscale) PCMGSetRScale(pc,i+1,rscale);
985:         VecDestroy(&rscale);
986:         MatDestroy(&p);
987:       }
988:       DMHasCreateRestriction(dms[i],&dmhasrestrict);
989:       if (dmhasrestrict && !mglevels[i+1]->restrct) {
990:         DMCreateRestriction(dms[i],dms[i+1],&p);
991:         PCMGSetRestriction(pc,i+1,p);
992:         MatDestroy(&p);
993:       }
994:       DMHasCreateInjection(dms[i],&dmhasinject);
995:       if (dmhasinject && !mglevels[i+1]->inject) {
996:         DMCreateInjection(dms[i],dms[i+1],&p);
997:         PCMGSetInjection(pc,i+1,p);
998:         MatDestroy(&p);
999:       }
1000:     }

1002:     for (i=n-2; i>-1; i--) DMDestroy(&dms[i]);
1003:     PetscFree(dms);
1004:   }

1006:   if (pc->dm && !pc->setupcalled) {
1007:     /* finest smoother also gets DM but it is not active, independent of whether galerkin==PC_MG_GALERKIN_EXTERNAL */
1008:     KSPSetDM(mglevels[n-1]->smoothd,pc->dm);
1009:     KSPSetDMActive(mglevels[n-1]->smoothd,PETSC_FALSE);
1010:     if (mglevels[n-1]->smoothd != mglevels[n-1]->smoothu) {
1011:       KSPSetDM(mglevels[n-1]->smoothu,pc->dm);
1012:       KSPSetDMActive(mglevels[n-1]->smoothu,PETSC_FALSE);
1013:     }
1014:     if (mglevels[n-1]->cr) {
1015:       KSPSetDM(mglevels[n-1]->cr,pc->dm);
1016:       KSPSetDMActive(mglevels[n-1]->cr,PETSC_FALSE);
1017:     }
1018:   }

1020:   if (mg->galerkin < PC_MG_GALERKIN_NONE) {
1021:     Mat       A,B;
1022:     PetscBool doA = PETSC_FALSE,doB = PETSC_FALSE;
1023:     MatReuse  reuse = MAT_INITIAL_MATRIX;

1025:     if ((mg->galerkin == PC_MG_GALERKIN_PMAT) || (mg->galerkin == PC_MG_GALERKIN_BOTH)) doB = PETSC_TRUE;
1026:     if ((mg->galerkin == PC_MG_GALERKIN_MAT) || ((mg->galerkin == PC_MG_GALERKIN_BOTH) && (dA != dB))) doA = PETSC_TRUE;
1027:     if (pc->setupcalled) reuse = MAT_REUSE_MATRIX;
1028:     for (i=n-2; i>-1; i--) {
1030:       if (!mglevels[i+1]->interpolate) {
1031:         PCMGSetInterpolation(pc,i+1,mglevels[i+1]->restrct);
1032:       }
1033:       if (!mglevels[i+1]->restrct) {
1034:         PCMGSetRestriction(pc,i+1,mglevels[i+1]->interpolate);
1035:       }
1036:       if (reuse == MAT_REUSE_MATRIX) {
1037:         KSPGetOperators(mglevels[i]->smoothd,&A,&B);
1038:       }
1039:       if (doA) {
1040:         MatGalerkin(mglevels[i+1]->restrct,dA,mglevels[i+1]->interpolate,reuse,1.0,&A);
1041:       }
1042:       if (doB) {
1043:         MatGalerkin(mglevels[i+1]->restrct,dB,mglevels[i+1]->interpolate,reuse,1.0,&B);
1044:       }
1045:       /* the management of the PetscObjectReference() and PetscObjecDereference() below is rather delicate */
1046:       if (!doA && dAeqdB) {
1047:         if (reuse == MAT_INITIAL_MATRIX) PetscObjectReference((PetscObject)B);
1048:         A = B;
1049:       } else if (!doA && reuse == MAT_INITIAL_MATRIX) {
1050:         KSPGetOperators(mglevels[i]->smoothd,&A,NULL);
1051:         PetscObjectReference((PetscObject)A);
1052:       }
1053:       if (!doB && dAeqdB) {
1054:         if (reuse == MAT_INITIAL_MATRIX) PetscObjectReference((PetscObject)A);
1055:         B = A;
1056:       } else if (!doB && reuse == MAT_INITIAL_MATRIX) {
1057:         KSPGetOperators(mglevels[i]->smoothd,NULL,&B);
1058:         PetscObjectReference((PetscObject)B);
1059:       }
1060:       if (reuse == MAT_INITIAL_MATRIX) {
1061:         KSPSetOperators(mglevels[i]->smoothd,A,B);
1062:         PetscObjectDereference((PetscObject)A);
1063:         PetscObjectDereference((PetscObject)B);
1064:       }
1065:       dA = A;
1066:       dB = B;
1067:     }
1068:   }

1070:   /* Adapt interpolation matrices */
1071:   if (mg->adaptInterpolation) {
1072:     mg->Nc = mg->Nc < 0 ? 6 : mg->Nc; /* Default to 6 modes */
1073:     for (i = 0; i < n; ++i) {
1074:       PCMGComputeCoarseSpace_Internal(pc, i, mg->coarseSpaceType, mg->Nc, !i ? NULL : mglevels[i-1]->coarseSpace, &mglevels[i]->coarseSpace);
1075:       if (i) PCMGAdaptInterpolator_Internal(pc, i, mglevels[i-1]->smoothu, mglevels[i]->smoothu, mg->Nc, mglevels[i-1]->coarseSpace, mglevels[i]->coarseSpace);
1076:     }
1077:     for (i = n-2; i > -1; --i) {
1078:       PCMGRecomputeLevelOperators_Internal(pc, i);
1079:     }
1080:   }

1082:   if (needRestricts && pc->dm) {
1083:     for (i=n-2; i>=0; i--) {
1084:       DM  dmfine,dmcoarse;
1085:       Mat Restrict,Inject;
1086:       Vec rscale;
1087:       KSPGetDM(mglevels[i+1]->smoothd,&dmfine);
1088:       KSPGetDM(mglevels[i]->smoothd,&dmcoarse);
1089:       PCMGGetRestriction(pc,i+1,&Restrict);
1090:       PCMGGetRScale(pc,i+1,&rscale);
1091:       PCMGGetInjection(pc,i+1,&Inject);
1092:       DMRestrict(dmfine,Restrict,rscale,Inject,dmcoarse);
1093:     }
1094:   }

1096:   if (!pc->setupcalled) {
1097:     for (i=0; i<n; i++) {
1098:       KSPSetFromOptions(mglevels[i]->smoothd);
1099:     }
1100:     for (i=1; i<n; i++) {
1101:       if (mglevels[i]->smoothu && (mglevels[i]->smoothu != mglevels[i]->smoothd)) {
1102:         KSPSetFromOptions(mglevels[i]->smoothu);
1103:       }
1104:       if (mglevels[i]->cr) {
1105:         KSPSetFromOptions(mglevels[i]->cr);
1106:       }
1107:     }
1108:     /* insure that if either interpolation or restriction is set the other other one is set */
1109:     for (i=1; i<n; i++) {
1110:       PCMGGetInterpolation(pc,i,NULL);
1111:       PCMGGetRestriction(pc,i,NULL);
1112:     }
1113:     for (i=0; i<n-1; i++) {
1114:       if (!mglevels[i]->b) {
1115:         Vec *vec;
1116:         KSPCreateVecs(mglevels[i]->smoothd,1,&vec,0,NULL);
1117:         PCMGSetRhs(pc,i,*vec);
1118:         VecDestroy(vec);
1119:         PetscFree(vec);
1120:       }
1121:       if (!mglevels[i]->r && i) {
1122:         VecDuplicate(mglevels[i]->b,&tvec);
1123:         PCMGSetR(pc,i,tvec);
1124:         VecDestroy(&tvec);
1125:       }
1126:       if (!mglevels[i]->x) {
1127:         VecDuplicate(mglevels[i]->b,&tvec);
1128:         PCMGSetX(pc,i,tvec);
1129:         VecDestroy(&tvec);
1130:       }
1131:       if (doCR) {
1132:         VecDuplicate(mglevels[i]->b,&mglevels[i]->crx);
1133:         VecDuplicate(mglevels[i]->b,&mglevels[i]->crb);
1134:         VecZeroEntries(mglevels[i]->crb);
1135:       }
1136:     }
1137:     if (n != 1 && !mglevels[n-1]->r) {
1138:       /* PCMGSetR() on the finest level if user did not supply it */
1139:       Vec *vec;
1140:       KSPCreateVecs(mglevels[n-1]->smoothd,1,&vec,0,NULL);
1141:       PCMGSetR(pc,n-1,*vec);
1142:       VecDestroy(vec);
1143:       PetscFree(vec);
1144:     }
1145:     if (doCR) {
1146:       VecDuplicate(mglevels[n-1]->r, &mglevels[n-1]->crx);
1147:       VecDuplicate(mglevels[n-1]->r, &mglevels[n-1]->crb);
1148:       VecZeroEntries(mglevels[n-1]->crb);
1149:     }
1150:   }

1152:   if (pc->dm) {
1153:     /* need to tell all the coarser levels to rebuild the matrix using the DM for that level */
1154:     for (i=0; i<n-1; i++) {
1155:       if (mglevels[i]->smoothd->setupstage != KSP_SETUP_NEW) mglevels[i]->smoothd->setupstage = KSP_SETUP_NEWMATRIX;
1156:     }
1157:   }
1158:   // We got here (PCSetUp_MG) because the matrix has changed, which means the smoother needs to be set up again (e.g.,
1159:   // new diagonal for Jacobi). Setting it here allows it to be logged under PCSetUp rather than deep inside a PCApply.
1160:   if (mglevels[n-1]->smoothd->setupstage != KSP_SETUP_NEW) mglevels[n-1]->smoothd->setupstage = KSP_SETUP_NEWMATRIX;

1162:   for (i=1; i<n; i++) {
1163:     if (mglevels[i]->smoothu == mglevels[i]->smoothd || mg->am == PC_MG_FULL || mg->am == PC_MG_KASKADE || mg->cyclesperpcapply > 1) {
1164:       /* if doing only down then initial guess is zero */
1165:       KSPSetInitialGuessNonzero(mglevels[i]->smoothd,PETSC_TRUE);
1166:     }
1167:     if (mglevels[i]->cr) KSPSetInitialGuessNonzero(mglevels[i]->cr,PETSC_TRUE);
1168:     if (mglevels[i]->eventsmoothsetup) PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);
1169:     KSPSetUp(mglevels[i]->smoothd);
1170:     if (mglevels[i]->smoothd->reason == KSP_DIVERGED_PC_FAILED) {
1171:       pc->failedreason = PC_SUBPC_ERROR;
1172:     }
1173:     if (mglevels[i]->eventsmoothsetup) PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);
1174:     if (!mglevels[i]->residual) {
1175:       Mat mat;
1176:       KSPGetOperators(mglevels[i]->smoothd,&mat,NULL);
1177:       PCMGSetResidual(pc,i,PCMGResidualDefault,mat);
1178:     }
1179:     if (!mglevels[i]->residualtranspose) {
1180:       Mat mat;
1181:       KSPGetOperators(mglevels[i]->smoothd,&mat,NULL);
1182:       PCMGSetResidualTranspose(pc,i,PCMGResidualTransposeDefault,mat);
1183:     }
1184:   }
1185:   for (i=1; i<n; i++) {
1186:     if (mglevels[i]->smoothu && mglevels[i]->smoothu != mglevels[i]->smoothd) {
1187:       Mat downmat,downpmat;

1189:       /* check if operators have been set for up, if not use down operators to set them */
1190:       KSPGetOperatorsSet(mglevels[i]->smoothu,&opsset,NULL);
1191:       if (!opsset) {
1192:         KSPGetOperators(mglevels[i]->smoothd,&downmat,&downpmat);
1193:         KSPSetOperators(mglevels[i]->smoothu,downmat,downpmat);
1194:       }

1196:       KSPSetInitialGuessNonzero(mglevels[i]->smoothu,PETSC_TRUE);
1197:       if (mglevels[i]->eventsmoothsetup) PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);
1198:       KSPSetUp(mglevels[i]->smoothu);
1199:       if (mglevels[i]->smoothu->reason == KSP_DIVERGED_PC_FAILED) {
1200:         pc->failedreason = PC_SUBPC_ERROR;
1201:       }
1202:       if (mglevels[i]->eventsmoothsetup) PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);
1203:     }
1204:     if (mglevels[i]->cr) {
1205:       Mat downmat,downpmat;

1207:       /* check if operators have been set for up, if not use down operators to set them */
1208:       KSPGetOperatorsSet(mglevels[i]->cr,&opsset,NULL);
1209:       if (!opsset) {
1210:         KSPGetOperators(mglevels[i]->smoothd,&downmat,&downpmat);
1211:         KSPSetOperators(mglevels[i]->cr,downmat,downpmat);
1212:       }

1214:       KSPSetInitialGuessNonzero(mglevels[i]->cr,PETSC_TRUE);
1215:       if (mglevels[i]->eventsmoothsetup) PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);
1216:       KSPSetUp(mglevels[i]->cr);
1217:       if (mglevels[i]->cr->reason == KSP_DIVERGED_PC_FAILED) {
1218:         pc->failedreason = PC_SUBPC_ERROR;
1219:       }
1220:       if (mglevels[i]->eventsmoothsetup) PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);
1221:     }
1222:   }

1224:   if (mglevels[0]->eventsmoothsetup) PetscLogEventBegin(mglevels[0]->eventsmoothsetup,0,0,0,0);
1225:   KSPSetUp(mglevels[0]->smoothd);
1226:   if (mglevels[0]->smoothd->reason == KSP_DIVERGED_PC_FAILED) {
1227:     pc->failedreason = PC_SUBPC_ERROR;
1228:   }
1229:   if (mglevels[0]->eventsmoothsetup) PetscLogEventEnd(mglevels[0]->eventsmoothsetup,0,0,0,0);

1231:   /*
1232:      Dump the interpolation/restriction matrices plus the
1233:    Jacobian/stiffness on each level. This allows MATLAB users to
1234:    easily check if the Galerkin condition A_c = R A_f R^T is satisfied.

1236:    Only support one or the other at the same time.
1237:   */
1238: #if defined(PETSC_USE_SOCKET_VIEWER)
1239:   PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_mg_dump_matlab",&dump,NULL);
1240:   if (dump) viewer = PETSC_VIEWER_SOCKET_(PetscObjectComm((PetscObject)pc));
1241:   dump = PETSC_FALSE;
1242: #endif
1243:   PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_mg_dump_binary",&dump,NULL);
1244:   if (dump) viewer = PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)pc));

1246:   if (viewer) {
1247:     for (i=1; i<n; i++) {
1248:       MatView(mglevels[i]->restrct,viewer);
1249:     }
1250:     for (i=0; i<n; i++) {
1251:       KSPGetPC(mglevels[i]->smoothd,&pc);
1252:       MatView(pc->mat,viewer);
1253:     }
1254:   }
1255:   return 0;
1256: }

1258: /* -------------------------------------------------------------------------------------*/

1260: PetscErrorCode PCMGGetLevels_MG(PC pc, PetscInt *levels)
1261: {
1262:   PC_MG *mg = (PC_MG *) pc->data;

1264:   *levels = mg->nlevels;
1265:   return 0;
1266: }

1268: /*@
1269:    PCMGGetLevels - Gets the number of levels to use with MG.

1271:    Not Collective

1273:    Input Parameter:
1274: .  pc - the preconditioner context

1276:    Output parameter:
1277: .  levels - the number of levels

1279:    Level: advanced

1281: .seealso: PCMGSetLevels()
1282: @*/
1283: PetscErrorCode PCMGGetLevels(PC pc,PetscInt *levels)
1284: {
1287:   *levels = 0;
1288:   PetscTryMethod(pc,"PCMGGetLevels_C",(PC,PetscInt*),(pc,levels));
1289:   return 0;
1290: }

1292: /*@
1293:    PCMGGetGridComplexity - compute operator and grid complexity of MG hierarchy

1295:    Input Parameter:
1296: .  pc - the preconditioner context

1298:    Output Parameters:
1299: +  gc - grid complexity = sum_i(n_i) / n_0
1300: -  oc - operator complexity = sum_i(nnz_i) / nnz_0

1302:    Level: advanced

1304: .seealso: PCMGGetLevels()
1305: @*/
1306: PetscErrorCode PCMGGetGridComplexity(PC pc, PetscReal *gc, PetscReal *oc)
1307: {
1308:   PC_MG          *mg      = (PC_MG*)pc->data;
1309:   PC_MG_Levels   **mglevels = mg->levels;
1310:   PetscInt       lev,N;
1311:   PetscLogDouble nnz0 = 0, sgc = 0, soc = 0, n0 = 0;
1312:   MatInfo        info;

1317:   if (!pc->setupcalled) {
1318:     if (gc) *gc = 0;
1319:     if (oc) *oc = 0;
1320:     return 0;
1321:   }
1323:   for (lev=0; lev<mg->nlevels; lev++) {
1324:     Mat dB;
1325:     KSPGetOperators(mglevels[lev]->smoothd,NULL,&dB);
1326:     MatGetInfo(dB,MAT_GLOBAL_SUM,&info); /* global reduction */
1327:     MatGetSize(dB,&N,NULL);
1328:     sgc += N;
1329:     soc += info.nz_used;
1330:     if (lev==mg->nlevels-1) {nnz0 = info.nz_used; n0 = N;}
1331:   }
1332:   if (n0 > 0 && gc) *gc = (PetscReal)(sgc/n0);
1333:   else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Number for grid points on finest level is not available");
1334:   if (nnz0 > 0 && oc) *oc = (PetscReal)(soc/nnz0);
1335:   return 0;
1336: }

1338: /*@
1339:    PCMGSetType - Determines the form of multigrid to use:
1340:    multiplicative, additive, full, or the Kaskade algorithm.

1342:    Logically Collective on PC

1344:    Input Parameters:
1345: +  pc - the preconditioner context
1346: -  form - multigrid form, one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,
1347:    PC_MG_FULL, PC_MG_KASKADE

1349:    Options Database Key:
1350: .  -pc_mg_type <form> - Sets <form>, one of multiplicative,
1351:    additive, full, kaskade

1353:    Level: advanced

1355: .seealso: PCMGSetLevels()
1356: @*/
1357: PetscErrorCode  PCMGSetType(PC pc,PCMGType form)
1358: {
1359:   PC_MG *mg = (PC_MG*)pc->data;

1363:   mg->am = form;
1364:   if (form == PC_MG_MULTIPLICATIVE) pc->ops->applyrichardson = PCApplyRichardson_MG;
1365:   else pc->ops->applyrichardson = NULL;
1366:   return 0;
1367: }

1369: /*@
1370:    PCMGGetType - Determines the form of multigrid to use:
1371:    multiplicative, additive, full, or the Kaskade algorithm.

1373:    Logically Collective on PC

1375:    Input Parameter:
1376: .  pc - the preconditioner context

1378:    Output Parameter:
1379: .  type - one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,PC_MG_FULL, PC_MG_KASKADE

1381:    Level: advanced

1383: .seealso: PCMGSetLevels()
1384: @*/
1385: PetscErrorCode  PCMGGetType(PC pc,PCMGType *type)
1386: {
1387:   PC_MG *mg = (PC_MG*)pc->data;

1390:   *type = mg->am;
1391:   return 0;
1392: }

1394: /*@
1395:    PCMGSetCycleType - Sets the type cycles to use.  Use PCMGSetCycleTypeOnLevel() for more
1396:    complicated cycling.

1398:    Logically Collective on PC

1400:    Input Parameters:
1401: +  pc - the multigrid context
1402: -  n - either PC_MG_CYCLE_V or PC_MG_CYCLE_W

1404:    Options Database Key:
1405: .  -pc_mg_cycle_type <v,w> - provide the cycle desired

1407:    Level: advanced

1409: .seealso: PCMGSetCycleTypeOnLevel()
1410: @*/
1411: PetscErrorCode  PCMGSetCycleType(PC pc,PCMGCycleType n)
1412: {
1413:   PC_MG        *mg        = (PC_MG*)pc->data;
1414:   PC_MG_Levels **mglevels = mg->levels;
1415:   PetscInt     i,levels;

1420:   levels = mglevels[0]->levels;
1421:   for (i=0; i<levels; i++) mglevels[i]->cycles = n;
1422:   return 0;
1423: }

1425: /*@
1426:    PCMGMultiplicativeSetCycles - Sets the number of cycles to use for each preconditioner step
1427:          of multigrid when PCMGType of PC_MG_MULTIPLICATIVE is used

1429:    Logically Collective on PC

1431:    Input Parameters:
1432: +  pc - the multigrid context
1433: -  n - number of cycles (default is 1)

1435:    Options Database Key:
1436: .  -pc_mg_multiplicative_cycles n - set the number of cycles

1438:    Level: advanced

1440:    Notes:
1441:     This is not associated with setting a v or w cycle, that is set with PCMGSetCycleType()

1443: .seealso: PCMGSetCycleTypeOnLevel(), PCMGSetCycleType()
1444: @*/
1445: PetscErrorCode  PCMGMultiplicativeSetCycles(PC pc,PetscInt n)
1446: {
1447:   PC_MG        *mg        = (PC_MG*)pc->data;

1451:   mg->cyclesperpcapply = n;
1452:   return 0;
1453: }

1455: PetscErrorCode PCMGSetGalerkin_MG(PC pc,PCMGGalerkinType use)
1456: {
1457:   PC_MG *mg = (PC_MG*)pc->data;

1459:   mg->galerkin = use;
1460:   return 0;
1461: }

1463: /*@
1464:    PCMGSetGalerkin - Causes the coarser grid matrices to be computed from the
1465:       finest grid via the Galerkin process: A_i-1 = r_i * A_i * p_i

1467:    Logically Collective on PC

1469:    Input Parameters:
1470: +  pc - the multigrid context
1471: -  use - one of PC_MG_GALERKIN_BOTH,PC_MG_GALERKIN_PMAT,PC_MG_GALERKIN_MAT, or PC_MG_GALERKIN_NONE

1473:    Options Database Key:
1474: .  -pc_mg_galerkin <both,pmat,mat,none> - set the matrices to form via the Galerkin process

1476:    Level: intermediate

1478:    Notes:
1479:     Some codes that use PCMG such as PCGAMG use Galerkin internally while constructing the hierarchy and thus do not
1480:      use the PCMG construction of the coarser grids.

1482: .seealso: PCMGGetGalerkin(), PCMGGalerkinType

1484: @*/
1485: PetscErrorCode PCMGSetGalerkin(PC pc,PCMGGalerkinType use)
1486: {
1488:   PetscTryMethod(pc,"PCMGSetGalerkin_C",(PC,PCMGGalerkinType),(pc,use));
1489:   return 0;
1490: }

1492: /*@
1493:    PCMGGetGalerkin - Checks if Galerkin multigrid is being used, i.e.
1494:       A_i-1 = r_i * A_i * p_i

1496:    Not Collective

1498:    Input Parameter:
1499: .  pc - the multigrid context

1501:    Output Parameter:
1502: .  galerkin - one of PC_MG_GALERKIN_BOTH,PC_MG_GALERKIN_PMAT,PC_MG_GALERKIN_MAT, PC_MG_GALERKIN_NONE, or PC_MG_GALERKIN_EXTERNAL

1504:    Level: intermediate

1506: .seealso: PCMGSetGalerkin(), PCMGGalerkinType

1508: @*/
1509: PetscErrorCode  PCMGGetGalerkin(PC pc,PCMGGalerkinType  *galerkin)
1510: {
1511:   PC_MG *mg = (PC_MG*)pc->data;

1514:   *galerkin = mg->galerkin;
1515:   return 0;
1516: }

1518: PetscErrorCode PCMGSetAdaptInterpolation_MG(PC pc, PetscBool adapt)
1519: {
1520:   PC_MG *mg = (PC_MG *) pc->data;

1522:   mg->adaptInterpolation = adapt;
1523:   return 0;
1524: }

1526: PetscErrorCode PCMGGetAdaptInterpolation_MG(PC pc, PetscBool *adapt)
1527: {
1528:   PC_MG *mg = (PC_MG *) pc->data;

1530:   *adapt = mg->adaptInterpolation;
1531:   return 0;
1532: }

1534: PetscErrorCode PCMGSetAdaptCR_MG(PC pc, PetscBool cr)
1535: {
1536:   PC_MG *mg = (PC_MG *) pc->data;

1538:   mg->compatibleRelaxation = cr;
1539:   return 0;
1540: }

1542: PetscErrorCode PCMGGetAdaptCR_MG(PC pc, PetscBool *cr)
1543: {
1544:   PC_MG *mg = (PC_MG *) pc->data;

1546:   *cr = mg->compatibleRelaxation;
1547:   return 0;
1548: }

1550: /*@
1551:   PCMGSetAdaptInterpolation - Adapt the interpolator based upon a vector space which should be accurately captured by the next coarser mesh, and thus accurately interpolated.

1553:   Logically Collective on PC

1555:   Input Parameters:
1556: + pc    - the multigrid context
1557: - adapt - flag for adaptation of the interpolator

1559:   Options Database Keys:
1560: + -pc_mg_adapt_interp                     - Turn on adaptation
1561: . -pc_mg_adapt_interp_n <int>             - The number of modes to use, should be divisible by dimension
1562: - -pc_mg_adapt_interp_coarse_space <type> - The type of coarse space: polynomial, harmonic, eigenvector, generalized_eigenvector

1564:   Level: intermediate

1566: .keywords: MG, set, Galerkin
1567: .seealso: PCMGGetAdaptInterpolation(), PCMGSetGalerkin()
1568: @*/
1569: PetscErrorCode PCMGSetAdaptInterpolation(PC pc, PetscBool adapt)
1570: {
1572:   PetscTryMethod(pc,"PCMGSetAdaptInterpolation_C",(PC,PetscBool),(pc,adapt));
1573:   return 0;
1574: }

1576: /*@
1577:   PCMGGetAdaptInterpolation - Get the flag to adapt the interpolator based upon a vector space which should be accurately captured by the next coarser mesh, and thus accurately interpolated.

1579:   Not collective

1581:   Input Parameter:
1582: . pc    - the multigrid context

1584:   Output Parameter:
1585: . adapt - flag for adaptation of the interpolator

1587:   Level: intermediate

1589: .keywords: MG, set, Galerkin
1590: .seealso: PCMGSetAdaptInterpolation(), PCMGSetGalerkin()
1591: @*/
1592: PetscErrorCode PCMGGetAdaptInterpolation(PC pc, PetscBool *adapt)
1593: {
1596:   PetscUseMethod(pc,"PCMGGetAdaptInterpolation_C",(PC,PetscBool*),(pc,adapt));
1597:   return 0;
1598: }

1600: /*@
1601:   PCMGSetAdaptCR - Monitor the coarse space quality using an auxiliary solve with compatible relaxation.

1603:   Logically Collective on PC

1605:   Input Parameters:
1606: + pc - the multigrid context
1607: - cr - flag for compatible relaxation

1609:   Options Database Keys:
1610: . -pc_mg_adapt_cr - Turn on compatible relaxation

1612:   Level: intermediate

1614: .keywords: MG, set, Galerkin
1615: .seealso: PCMGGetAdaptCR(), PCMGSetAdaptInterpolation(), PCMGSetGalerkin()
1616: @*/
1617: PetscErrorCode PCMGSetAdaptCR(PC pc, PetscBool cr)
1618: {
1620:   PetscTryMethod(pc,"PCMGSetAdaptCR_C",(PC,PetscBool),(pc,cr));
1621:   return 0;
1622: }

1624: /*@
1625:   PCMGGetAdaptCR - Get the flag to monitor coarse space quality using an auxiliary solve with compatible relaxation.

1627:   Not collective

1629:   Input Parameter:
1630: . pc    - the multigrid context

1632:   Output Parameter:
1633: . cr - flag for compatible relaxaion

1635:   Level: intermediate

1637: .keywords: MG, set, Galerkin
1638: .seealso: PCMGSetAdaptCR(), PCMGGetAdaptInterpolation(), PCMGSetGalerkin()
1639: @*/
1640: PetscErrorCode PCMGGetAdaptCR(PC pc, PetscBool *cr)
1641: {
1644:   PetscUseMethod(pc,"PCMGGetAdaptCR_C",(PC,PetscBool*),(pc,cr));
1645:   return 0;
1646: }

1648: /*@
1649:    PCMGSetNumberSmooth - Sets the number of pre and post-smoothing steps to use
1650:    on all levels.  Use PCMGDistinctSmoothUp() to create separate up and down smoothers if you want different numbers of
1651:    pre- and post-smoothing steps.

1653:    Logically Collective on PC

1655:    Input Parameters:
1656: +  mg - the multigrid context
1657: -  n - the number of smoothing steps

1659:    Options Database Key:
1660: .  -mg_levels_ksp_max_it <n> - Sets number of pre and post-smoothing steps

1662:    Level: advanced

1664:    Notes:
1665:     this does not set a value on the coarsest grid, since we assume that
1666:     there is no separate smooth up on the coarsest grid.

1668: .seealso: PCMGSetDistinctSmoothUp()
1669: @*/
1670: PetscErrorCode  PCMGSetNumberSmooth(PC pc,PetscInt n)
1671: {
1672:   PC_MG          *mg        = (PC_MG*)pc->data;
1673:   PC_MG_Levels   **mglevels = mg->levels;
1674:   PetscInt       i,levels;

1679:   levels = mglevels[0]->levels;

1681:   for (i=1; i<levels; i++) {
1682:     KSPSetTolerances(mglevels[i]->smoothu,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);
1683:     KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);
1684:     mg->default_smoothu = n;
1685:     mg->default_smoothd = n;
1686:   }
1687:   return 0;
1688: }

1690: /*@
1691:    PCMGSetDistinctSmoothUp - sets the up (post) smoother to be a separate KSP from the down (pre) smoother on all levels
1692:        and adds the suffix _up to the options name

1694:    Logically Collective on PC

1696:    Input Parameters:
1697: .  pc - the preconditioner context

1699:    Options Database Key:
1700: .  -pc_mg_distinct_smoothup <bool> - use distinct smoothing objects

1702:    Level: advanced

1704:    Notes:
1705:     this does not set a value on the coarsest grid, since we assume that
1706:     there is no separate smooth up on the coarsest grid.

1708: .seealso: PCMGSetNumberSmooth()
1709: @*/
1710: PetscErrorCode  PCMGSetDistinctSmoothUp(PC pc)
1711: {
1712:   PC_MG          *mg        = (PC_MG*)pc->data;
1713:   PC_MG_Levels   **mglevels = mg->levels;
1714:   PetscInt       i,levels;
1715:   KSP            subksp;

1719:   levels = mglevels[0]->levels;

1721:   for (i=1; i<levels; i++) {
1722:     const char *prefix = NULL;
1723:     /* make sure smoother up and down are different */
1724:     PCMGGetSmootherUp(pc,i,&subksp);
1725:     KSPGetOptionsPrefix(mglevels[i]->smoothd,&prefix);
1726:     KSPSetOptionsPrefix(subksp,prefix);
1727:     KSPAppendOptionsPrefix(subksp,"up_");
1728:   }
1729:   return 0;
1730: }

1732: /* No new matrices are created, and the coarse operator matrices are the references to the original ones */
1733: PetscErrorCode  PCGetInterpolations_MG(PC pc,PetscInt *num_levels,Mat *interpolations[])
1734: {
1735:   PC_MG          *mg        = (PC_MG*)pc->data;
1736:   PC_MG_Levels   **mglevels = mg->levels;
1737:   Mat            *mat;
1738:   PetscInt       l;

1741:   PetscMalloc1(mg->nlevels,&mat);
1742:   for (l=1; l< mg->nlevels; l++) {
1743:     mat[l-1] = mglevels[l]->interpolate;
1744:     PetscObjectReference((PetscObject)mat[l-1]);
1745:   }
1746:   *num_levels = mg->nlevels;
1747:   *interpolations = mat;
1748:   return 0;
1749: }

1751: /* No new matrices are created, and the coarse operator matrices are the references to the original ones */
1752: PetscErrorCode  PCGetCoarseOperators_MG(PC pc,PetscInt *num_levels,Mat *coarseOperators[])
1753: {
1754:   PC_MG          *mg        = (PC_MG*)pc->data;
1755:   PC_MG_Levels   **mglevels = mg->levels;
1756:   PetscInt       l;
1757:   Mat            *mat;

1760:   PetscMalloc1(mg->nlevels,&mat);
1761:   for (l=0; l<mg->nlevels-1; l++) {
1762:     KSPGetOperators(mglevels[l]->smoothd,NULL,&(mat[l]));
1763:     PetscObjectReference((PetscObject)mat[l]);
1764:   }
1765:   *num_levels = mg->nlevels;
1766:   *coarseOperators = mat;
1767:   return 0;
1768: }

1770: /*@C
1771:   PCMGRegisterCoarseSpaceConstructor -  Adds a method to the PCMG package for coarse space construction.

1773:   Not collective

1775:   Input Parameters:
1776: + name     - name of the constructor
1777: - function - constructor routine

1779:   Notes:
1780:   Calling sequence for the routine:
1781: $ my_csp(PC pc, PetscInt l, DM dm, KSP smooth, PetscInt Nc, const Vec initGuess[], Vec **coarseSp)
1782: $   pc        - The PC object
1783: $   l         - The multigrid level, 0 is the coarse level
1784: $   dm        - The DM for this level
1785: $   smooth    - The level smoother
1786: $   Nc        - The size of the coarse space
1787: $   initGuess - Basis for an initial guess for the space
1788: $   coarseSp  - A basis for the computed coarse space

1790:   Level: advanced

1792: .seealso: PCMGGetCoarseSpaceConstructor(), PCRegister()
1793: @*/
1794: PetscErrorCode PCMGRegisterCoarseSpaceConstructor(const char name[], PetscErrorCode (*function)(PC, PetscInt, DM, KSP, PetscInt, const Vec[], Vec **))
1795: {
1796:   PCInitializePackage();
1797:   PetscFunctionListAdd(&PCMGCoarseList,name,function);
1798:   return 0;
1799: }

1801: /*@C
1802:   PCMGGetCoarseSpaceConstructor -  Returns the given coarse space construction method.

1804:   Not collective

1806:   Input Parameter:
1807: . name     - name of the constructor

1809:   Output Parameter:
1810: . function - constructor routine

1812:   Notes:
1813:   Calling sequence for the routine:
1814: $ my_csp(PC pc, PetscInt l, DM dm, KSP smooth, PetscInt Nc, const Vec initGuess[], Vec **coarseSp)
1815: $   pc        - The PC object
1816: $   l         - The multigrid level, 0 is the coarse level
1817: $   dm        - The DM for this level
1818: $   smooth    - The level smoother
1819: $   Nc        - The size of the coarse space
1820: $   initGuess - Basis for an initial guess for the space
1821: $   coarseSp  - A basis for the computed coarse space

1823:   Level: advanced

1825: .seealso: PCMGRegisterCoarseSpaceConstructor(), PCRegister()
1826: @*/
1827: PetscErrorCode PCMGGetCoarseSpaceConstructor(const char name[], PetscErrorCode (**function)(PC, PetscInt, DM, KSP, PetscInt, const Vec[], Vec **))
1828: {
1829:   PetscFunctionListFind(PCMGCoarseList,name,function);
1830:   return 0;
1831: }

1833: /* ----------------------------------------------------------------------------------------*/

1835: /*MC
1836:    PCMG - Use multigrid preconditioning. This preconditioner requires you provide additional
1837:     information about the coarser grid matrices and restriction/interpolation operators.

1839:    Options Database Keys:
1840: +  -pc_mg_levels <nlevels> - number of levels including finest
1841: .  -pc_mg_cycle_type <v,w> - provide the cycle desired
1842: .  -pc_mg_type <additive,multiplicative,full,kaskade> - multiplicative is the default
1843: .  -pc_mg_log - log information about time spent on each level of the solver
1844: .  -pc_mg_distinct_smoothup - configure up (after interpolation) and down (before restriction) smoothers separately (with different options prefixes)
1845: .  -pc_mg_galerkin <both,pmat,mat,none> - use Galerkin process to compute coarser operators, i.e. Acoarse = R A R'
1846: .  -pc_mg_multiplicative_cycles - number of cycles to use as the preconditioner (defaults to 1)
1847: .  -pc_mg_dump_matlab - dumps the matrices for each level and the restriction/interpolation matrices
1848:                         to the Socket viewer for reading from MATLAB.
1849: -  -pc_mg_dump_binary - dumps the matrices for each level and the restriction/interpolation matrices
1850:                         to the binary output file called binaryoutput

1852:    Notes:
1853:     If one uses a Krylov method such GMRES or CG as the smoother then one must use KSPFGMRES, KSPGCR, or KSPRICHARDSON as the outer Krylov method

1855:        When run with a single level the smoother options are used on that level NOT the coarse grid solver options

1857:        When run with KSPRICHARDSON the convergence test changes slightly if monitor is turned on. The iteration count may change slightly. This
1858:        is because without monitoring the residual norm is computed WITHIN each multigrid cycle on the finest level after the pre-smoothing
1859:        (because the residual has just been computed for the multigrid algorithm and is hence available for free) while with monitoring the
1860:        residual is computed at the end of each cycle.

1862:    Level: intermediate

1864: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType, PCEXOTIC, PCGAMG, PCML, PCHYPRE
1865:            PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), PCMGSetCycleType(),
1866:            PCMGSetDistinctSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(),
1867:            PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(),
1868:            PCMGSetCycleTypeOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR()
1869: M*/

1871: PETSC_EXTERN PetscErrorCode PCCreate_MG(PC pc)
1872: {
1873:   PC_MG          *mg;

1875:   PetscNewLog(pc,&mg);
1876:   pc->data     = mg;
1877:   mg->nlevels  = -1;
1878:   mg->am       = PC_MG_MULTIPLICATIVE;
1879:   mg->galerkin = PC_MG_GALERKIN_NONE;
1880:   mg->adaptInterpolation = PETSC_FALSE;
1881:   mg->Nc                 = -1;
1882:   mg->eigenvalue         = -1;

1884:   pc->useAmat = PETSC_TRUE;

1886:   pc->ops->apply          = PCApply_MG;
1887:   pc->ops->applytranspose = PCApplyTranspose_MG;
1888:   pc->ops->matapply       = PCMatApply_MG;
1889:   pc->ops->setup          = PCSetUp_MG;
1890:   pc->ops->reset          = PCReset_MG;
1891:   pc->ops->destroy        = PCDestroy_MG;
1892:   pc->ops->setfromoptions = PCSetFromOptions_MG;
1893:   pc->ops->view           = PCView_MG;

1895:   PetscObjectComposedDataRegister(&mg->eigenvalue);
1896:   PetscObjectComposeFunction((PetscObject)pc,"PCMGSetGalerkin_C",PCMGSetGalerkin_MG);
1897:   PetscObjectComposeFunction((PetscObject)pc,"PCMGGetLevels_C",PCMGGetLevels_MG);
1898:   PetscObjectComposeFunction((PetscObject)pc,"PCMGSetLevels_C",PCMGSetLevels_MG);
1899:   PetscObjectComposeFunction((PetscObject)pc,"PCGetInterpolations_C",PCGetInterpolations_MG);
1900:   PetscObjectComposeFunction((PetscObject)pc,"PCGetCoarseOperators_C",PCGetCoarseOperators_MG);
1901:   PetscObjectComposeFunction((PetscObject)pc,"PCMGSetAdaptInterpolation_C",PCMGSetAdaptInterpolation_MG);
1902:   PetscObjectComposeFunction((PetscObject)pc,"PCMGGetAdaptInterpolation_C",PCMGGetAdaptInterpolation_MG);
1903:   PetscObjectComposeFunction((PetscObject)pc,"PCMGSetAdaptCR_C",PCMGSetAdaptCR_MG);
1904:   PetscObjectComposeFunction((PetscObject)pc,"PCMGGetAdaptCR_C",PCMGGetAdaptCR_MG);
1905:   return 0;
1906: }