1
2 /*
3 Defines the multigrid preconditioner interface.
4 */
5 #include <petsc/private/pcmgimpl.h> /*I "petscksp.h" I*/
6 #include <petscdm.h>
7 PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC,PetscBool*);
8
9 /*
10 Contains the list of registered coarse space construction routines
11 */
12 PetscFunctionList PCMGCoarseList = NULL;
13
PCMGMCycle_Private(PC pc,PC_MG_Levels ** mglevelsin,PCRichardsonConvergedReason * reason)14 PetscErrorCode PCMGMCycle_Private(PC pc,PC_MG_Levels **mglevelsin,PCRichardsonConvergedReason *reason)
15 {
16 PC_MG *mg = (PC_MG*)pc->data;
17 PC_MG_Levels *mgc,*mglevels = *mglevelsin;
18 PetscErrorCode ierr;
19 PetscInt cycles = (mglevels->level == 1) ? 1 : (PetscInt) mglevels->cycles;
20
21 PetscFunctionBegin;
22 if (mglevels->eventsmoothsolve) {ierr = PetscLogEventBegin(mglevels->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
23 ierr = KSPSolve(mglevels->smoothd,mglevels->b,mglevels->x);CHKERRQ(ierr); /* pre-smooth */
24 ierr = KSPCheckSolve(mglevels->smoothd,pc,mglevels->x);CHKERRQ(ierr);
25 if (mglevels->eventsmoothsolve) {ierr = PetscLogEventEnd(mglevels->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
26 if (mglevels->level) { /* not the coarsest grid */
27 if (mglevels->eventresidual) {ierr = PetscLogEventBegin(mglevels->eventresidual,0,0,0,0);CHKERRQ(ierr);}
28 ierr = (*mglevels->residual)(mglevels->A,mglevels->b,mglevels->x,mglevels->r);CHKERRQ(ierr);
29 if (mglevels->eventresidual) {ierr = PetscLogEventEnd(mglevels->eventresidual,0,0,0,0);CHKERRQ(ierr);}
30
31 /* if on finest level and have convergence criteria set */
32 if (mglevels->level == mglevels->levels-1 && mg->ttol && reason) {
33 PetscReal rnorm;
34 ierr = VecNorm(mglevels->r,NORM_2,&rnorm);CHKERRQ(ierr);
35 if (rnorm <= mg->ttol) {
36 if (rnorm < mg->abstol) {
37 *reason = PCRICHARDSON_CONVERGED_ATOL;
38 ierr = PetscInfo2(pc,"Linear solver has converged. Residual norm %g is less than absolute tolerance %g\n",(double)rnorm,(double)mg->abstol);CHKERRQ(ierr);
39 } else {
40 *reason = PCRICHARDSON_CONVERGED_RTOL;
41 ierr = PetscInfo2(pc,"Linear solver has converged. Residual norm %g is less than relative tolerance times initial residual norm %g\n",(double)rnorm,(double)mg->ttol);CHKERRQ(ierr);
42 }
43 PetscFunctionReturn(0);
44 }
45 }
46
47 mgc = *(mglevelsin - 1);
48 if (mglevels->eventinterprestrict) {ierr = PetscLogEventBegin(mglevels->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
49 ierr = MatRestrict(mglevels->restrct,mglevels->r,mgc->b);CHKERRQ(ierr);
50 if (mglevels->eventinterprestrict) {ierr = PetscLogEventEnd(mglevels->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
51 ierr = VecSet(mgc->x,0.0);CHKERRQ(ierr);
52 while (cycles--) {
53 ierr = PCMGMCycle_Private(pc,mglevelsin-1,reason);CHKERRQ(ierr);
54 }
55 if (mglevels->eventinterprestrict) {ierr = PetscLogEventBegin(mglevels->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
56 ierr = MatInterpolateAdd(mglevels->interpolate,mgc->x,mglevels->x,mglevels->x);CHKERRQ(ierr);
57 if (mglevels->eventinterprestrict) {ierr = PetscLogEventEnd(mglevels->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
58 if (mglevels->eventsmoothsolve) {ierr = PetscLogEventBegin(mglevels->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
59 ierr = KSPSolve(mglevels->smoothu,mglevels->b,mglevels->x);CHKERRQ(ierr); /* post smooth */
60 ierr = KSPCheckSolve(mglevels->smoothu,pc,mglevels->x);CHKERRQ(ierr);
61 if (mglevels->eventsmoothsolve) {ierr = PetscLogEventEnd(mglevels->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
62 }
63 PetscFunctionReturn(0);
64 }
65
PCApplyRichardson_MG(PC pc,Vec b,Vec x,Vec w,PetscReal rtol,PetscReal abstol,PetscReal dtol,PetscInt its,PetscBool zeroguess,PetscInt * outits,PCRichardsonConvergedReason * reason)66 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)
67 {
68 PC_MG *mg = (PC_MG*)pc->data;
69 PC_MG_Levels **mglevels = mg->levels;
70 PetscErrorCode ierr;
71 PC tpc;
72 PetscBool changeu,changed;
73 PetscInt levels = mglevels[0]->levels,i;
74
75 PetscFunctionBegin;
76 /* When the DM is supplying the matrix then it will not exist until here */
77 for (i=0; i<levels; i++) {
78 if (!mglevels[i]->A) {
79 ierr = KSPGetOperators(mglevels[i]->smoothu,&mglevels[i]->A,NULL);CHKERRQ(ierr);
80 ierr = PetscObjectReference((PetscObject)mglevels[i]->A);CHKERRQ(ierr);
81 }
82 }
83
84 ierr = KSPGetPC(mglevels[levels-1]->smoothd,&tpc);CHKERRQ(ierr);
85 ierr = PCPreSolveChangeRHS(tpc,&changed);CHKERRQ(ierr);
86 ierr = KSPGetPC(mglevels[levels-1]->smoothu,&tpc);CHKERRQ(ierr);
87 ierr = PCPreSolveChangeRHS(tpc,&changeu);CHKERRQ(ierr);
88 if (!changed && !changeu) {
89 ierr = VecDestroy(&mglevels[levels-1]->b);CHKERRQ(ierr);
90 mglevels[levels-1]->b = b;
91 } else { /* if the smoother changes the rhs during PreSolve, we cannot use the input vector */
92 if (!mglevels[levels-1]->b) {
93 Vec *vec;
94
95 ierr = KSPCreateVecs(mglevels[levels-1]->smoothd,1,&vec,0,NULL);CHKERRQ(ierr);
96 mglevels[levels-1]->b = *vec;
97 ierr = PetscFree(vec);CHKERRQ(ierr);
98 }
99 ierr = VecCopy(b,mglevels[levels-1]->b);CHKERRQ(ierr);
100 }
101 mglevels[levels-1]->x = x;
102
103 mg->rtol = rtol;
104 mg->abstol = abstol;
105 mg->dtol = dtol;
106 if (rtol) {
107 /* compute initial residual norm for relative convergence test */
108 PetscReal rnorm;
109 if (zeroguess) {
110 ierr = VecNorm(b,NORM_2,&rnorm);CHKERRQ(ierr);
111 } else {
112 ierr = (*mglevels[levels-1]->residual)(mglevels[levels-1]->A,b,x,w);CHKERRQ(ierr);
113 ierr = VecNorm(w,NORM_2,&rnorm);CHKERRQ(ierr);
114 }
115 mg->ttol = PetscMax(rtol*rnorm,abstol);
116 } else if (abstol) mg->ttol = abstol;
117 else mg->ttol = 0.0;
118
119 /* since smoother is applied to full system, not just residual we need to make sure that smoothers don't
120 stop prematurely due to small residual */
121 for (i=1; i<levels; i++) {
122 ierr = KSPSetTolerances(mglevels[i]->smoothu,0,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
123 if (mglevels[i]->smoothu != mglevels[i]->smoothd) {
124 /* For Richardson the initial guess is nonzero since it is solving in each cycle the original system not just applying as a preconditioner */
125 ierr = KSPSetInitialGuessNonzero(mglevels[i]->smoothd,PETSC_TRUE);CHKERRQ(ierr);
126 ierr = KSPSetTolerances(mglevels[i]->smoothd,0,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
127 }
128 }
129
130 *reason = (PCRichardsonConvergedReason)0;
131 for (i=0; i<its; i++) {
132 ierr = PCMGMCycle_Private(pc,mglevels+levels-1,reason);CHKERRQ(ierr);
133 if (*reason) break;
134 }
135 if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS;
136 *outits = i;
137 if (!changed && !changeu) mglevels[levels-1]->b = NULL;
138 PetscFunctionReturn(0);
139 }
140
PCReset_MG(PC pc)141 PetscErrorCode PCReset_MG(PC pc)
142 {
143 PC_MG *mg = (PC_MG*)pc->data;
144 PC_MG_Levels **mglevels = mg->levels;
145 PetscErrorCode ierr;
146 PetscInt i,c,n;
147
148 PetscFunctionBegin;
149 if (mglevels) {
150 n = mglevels[0]->levels;
151 for (i=0; i<n-1; i++) {
152 ierr = VecDestroy(&mglevels[i+1]->r);CHKERRQ(ierr);
153 ierr = VecDestroy(&mglevels[i]->b);CHKERRQ(ierr);
154 ierr = VecDestroy(&mglevels[i]->x);CHKERRQ(ierr);
155 ierr = MatDestroy(&mglevels[i+1]->restrct);CHKERRQ(ierr);
156 ierr = MatDestroy(&mglevels[i+1]->interpolate);CHKERRQ(ierr);
157 ierr = MatDestroy(&mglevels[i+1]->inject);CHKERRQ(ierr);
158 ierr = VecDestroy(&mglevels[i+1]->rscale);CHKERRQ(ierr);
159 }
160 /* this is not null only if the smoother on the finest level
161 changes the rhs during PreSolve */
162 ierr = VecDestroy(&mglevels[n-1]->b);CHKERRQ(ierr);
163
164 for (i=0; i<n; i++) {
165 if (mglevels[i]->coarseSpace) for (c = 0; c < mg->Nc; ++c) {ierr = VecDestroy(&mglevels[i]->coarseSpace[c]);CHKERRQ(ierr);}
166 ierr = PetscFree(mglevels[i]->coarseSpace);CHKERRQ(ierr);
167 mglevels[i]->coarseSpace = NULL;
168 ierr = MatDestroy(&mglevels[i]->A);CHKERRQ(ierr);
169 if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
170 ierr = KSPReset(mglevels[i]->smoothd);CHKERRQ(ierr);
171 }
172 ierr = KSPReset(mglevels[i]->smoothu);CHKERRQ(ierr);
173 }
174 mg->Nc = 0;
175 }
176 PetscFunctionReturn(0);
177 }
178
PCMGSetLevels_MG(PC pc,PetscInt levels,MPI_Comm * comms)179 PetscErrorCode PCMGSetLevels_MG(PC pc,PetscInt levels,MPI_Comm *comms)
180 {
181 PetscErrorCode ierr;
182 PC_MG *mg = (PC_MG*)pc->data;
183 MPI_Comm comm;
184 PC_MG_Levels **mglevels = mg->levels;
185 PCMGType mgtype = mg->am;
186 PetscInt mgctype = (PetscInt) PC_MG_CYCLE_V;
187 PetscInt i;
188 PetscMPIInt size;
189 const char *prefix;
190 PC ipc;
191 PetscInt n;
192
193 PetscFunctionBegin;
194 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
195 PetscValidLogicalCollectiveInt(pc,levels,2);
196 if (mg->nlevels == levels) PetscFunctionReturn(0);
197 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
198 if (mglevels) {
199 mgctype = mglevels[0]->cycles;
200 /* changing the number of levels so free up the previous stuff */
201 ierr = PCReset_MG(pc);CHKERRQ(ierr);
202 n = mglevels[0]->levels;
203 for (i=0; i<n; i++) {
204 if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
205 ierr = KSPDestroy(&mglevels[i]->smoothd);CHKERRQ(ierr);
206 }
207 ierr = KSPDestroy(&mglevels[i]->smoothu);CHKERRQ(ierr);
208 ierr = PetscFree(mglevels[i]);CHKERRQ(ierr);
209 }
210 ierr = PetscFree(mg->levels);CHKERRQ(ierr);
211 }
212
213 mg->nlevels = levels;
214
215 ierr = PetscMalloc1(levels,&mglevels);CHKERRQ(ierr);
216 ierr = PetscLogObjectMemory((PetscObject)pc,levels*(sizeof(PC_MG*)));CHKERRQ(ierr);
217
218 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr);
219
220 mg->stageApply = 0;
221 for (i=0; i<levels; i++) {
222 ierr = PetscNewLog(pc,&mglevels[i]);CHKERRQ(ierr);
223
224 mglevels[i]->level = i;
225 mglevels[i]->levels = levels;
226 mglevels[i]->cycles = mgctype;
227 mg->default_smoothu = 2;
228 mg->default_smoothd = 2;
229 mglevels[i]->eventsmoothsetup = 0;
230 mglevels[i]->eventsmoothsolve = 0;
231 mglevels[i]->eventresidual = 0;
232 mglevels[i]->eventinterprestrict = 0;
233
234 if (comms) comm = comms[i];
235 ierr = KSPCreate(comm,&mglevels[i]->smoothd);CHKERRQ(ierr);
236 ierr = KSPSetErrorIfNotConverged(mglevels[i]->smoothd,pc->erroriffailure);CHKERRQ(ierr);
237 ierr = PetscObjectIncrementTabLevel((PetscObject)mglevels[i]->smoothd,(PetscObject)pc,levels-i);CHKERRQ(ierr);
238 ierr = KSPSetOptionsPrefix(mglevels[i]->smoothd,prefix);CHKERRQ(ierr);
239 ierr = PetscObjectComposedDataSetInt((PetscObject) mglevels[i]->smoothd, PetscMGLevelId, mglevels[i]->level);CHKERRQ(ierr);
240 if (i || levels == 1) {
241 char tprefix[128];
242
243 ierr = KSPSetType(mglevels[i]->smoothd,KSPCHEBYSHEV);CHKERRQ(ierr);
244 ierr = KSPSetConvergenceTest(mglevels[i]->smoothd,KSPConvergedSkip,NULL,NULL);CHKERRQ(ierr);
245 ierr = KSPSetNormType(mglevels[i]->smoothd,KSP_NORM_NONE);CHKERRQ(ierr);
246 ierr = KSPGetPC(mglevels[i]->smoothd,&ipc);CHKERRQ(ierr);
247 ierr = PCSetType(ipc,PCSOR);CHKERRQ(ierr);
248 ierr = KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT, mg->default_smoothd);CHKERRQ(ierr);
249
250 sprintf(tprefix,"mg_levels_%d_",(int)i);
251 ierr = KSPAppendOptionsPrefix(mglevels[i]->smoothd,tprefix);CHKERRQ(ierr);
252 } else {
253 ierr = KSPAppendOptionsPrefix(mglevels[0]->smoothd,"mg_coarse_");CHKERRQ(ierr);
254
255 /* coarse solve is (redundant) LU by default; set shifttype NONZERO to avoid annoying zero-pivot in LU preconditioner */
256 ierr = KSPSetType(mglevels[0]->smoothd,KSPPREONLY);CHKERRQ(ierr);
257 ierr = KSPGetPC(mglevels[0]->smoothd,&ipc);CHKERRQ(ierr);
258 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
259 if (size > 1) {
260 ierr = PCSetType(ipc,PCREDUNDANT);CHKERRQ(ierr);
261 } else {
262 ierr = PCSetType(ipc,PCLU);CHKERRQ(ierr);
263 }
264 ierr = PCFactorSetShiftType(ipc,MAT_SHIFT_INBLOCKS);CHKERRQ(ierr);
265 }
266 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)mglevels[i]->smoothd);CHKERRQ(ierr);
267
268 mglevels[i]->smoothu = mglevels[i]->smoothd;
269 mg->rtol = 0.0;
270 mg->abstol = 0.0;
271 mg->dtol = 0.0;
272 mg->ttol = 0.0;
273 mg->cyclesperpcapply = 1;
274 }
275 mg->levels = mglevels;
276 ierr = PCMGSetType(pc,mgtype);CHKERRQ(ierr);
277 PetscFunctionReturn(0);
278 }
279
280 /*@C
281 PCMGSetLevels - Sets the number of levels to use with MG.
282 Must be called before any other MG routine.
283
284 Logically Collective on PC
285
286 Input Parameters:
287 + pc - the preconditioner context
288 . levels - the number of levels
289 - comms - optional communicators for each level; this is to allow solving the coarser problems
290 on smaller sets of processors.
291
292 Level: intermediate
293
294 Notes:
295 If the number of levels is one then the multigrid uses the -mg_levels prefix
296 for setting the level options rather than the -mg_coarse prefix.
297
298 .seealso: PCMGSetType(), PCMGGetLevels()
299 @*/
PCMGSetLevels(PC pc,PetscInt levels,MPI_Comm * comms)300 PetscErrorCode PCMGSetLevels(PC pc,PetscInt levels,MPI_Comm *comms)
301 {
302 PetscErrorCode ierr;
303
304 PetscFunctionBegin;
305 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
306 if (comms) PetscValidPointer(comms,3);
307 ierr = PetscTryMethod(pc,"PCMGSetLevels_C",(PC,PetscInt,MPI_Comm*),(pc,levels,comms));CHKERRQ(ierr);
308 PetscFunctionReturn(0);
309 }
310
311
PCDestroy_MG(PC pc)312 PetscErrorCode PCDestroy_MG(PC pc)
313 {
314 PetscErrorCode ierr;
315 PC_MG *mg = (PC_MG*)pc->data;
316 PC_MG_Levels **mglevels = mg->levels;
317 PetscInt i,n;
318
319 PetscFunctionBegin;
320 ierr = PCReset_MG(pc);CHKERRQ(ierr);
321 if (mglevels) {
322 n = mglevels[0]->levels;
323 for (i=0; i<n; i++) {
324 if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
325 ierr = KSPDestroy(&mglevels[i]->smoothd);CHKERRQ(ierr);
326 }
327 ierr = KSPDestroy(&mglevels[i]->smoothu);CHKERRQ(ierr);
328 ierr = PetscFree(mglevels[i]);CHKERRQ(ierr);
329 }
330 ierr = PetscFree(mg->levels);CHKERRQ(ierr);
331 }
332 ierr = PetscFree(pc->data);CHKERRQ(ierr);
333 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGetInterpolations_C",NULL);CHKERRQ(ierr);
334 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGetCoarseOperators_C",NULL);CHKERRQ(ierr);
335 PetscFunctionReturn(0);
336 }
337
338
339
340 extern PetscErrorCode PCMGACycle_Private(PC,PC_MG_Levels**);
341 extern PetscErrorCode PCMGFCycle_Private(PC,PC_MG_Levels**);
342 extern PetscErrorCode PCMGKCycle_Private(PC,PC_MG_Levels**);
343
344 /*
345 PCApply_MG - Runs either an additive, multiplicative, Kaskadic
346 or full cycle of multigrid.
347
348 Note:
349 A simple wrapper which calls PCMGMCycle(),PCMGACycle(), or PCMGFCycle().
350 */
PCApply_MG(PC pc,Vec b,Vec x)351 static PetscErrorCode PCApply_MG(PC pc,Vec b,Vec x)
352 {
353 PC_MG *mg = (PC_MG*)pc->data;
354 PC_MG_Levels **mglevels = mg->levels;
355 PetscErrorCode ierr;
356 PC tpc;
357 PetscInt levels = mglevels[0]->levels,i;
358 PetscBool changeu,changed;
359
360 PetscFunctionBegin;
361 if (mg->stageApply) {ierr = PetscLogStagePush(mg->stageApply);CHKERRQ(ierr);}
362 /* When the DM is supplying the matrix then it will not exist until here */
363 for (i=0; i<levels; i++) {
364 if (!mglevels[i]->A) {
365 ierr = KSPGetOperators(mglevels[i]->smoothu,&mglevels[i]->A,NULL);CHKERRQ(ierr);
366 ierr = PetscObjectReference((PetscObject)mglevels[i]->A);CHKERRQ(ierr);
367 }
368 }
369
370 ierr = KSPGetPC(mglevels[levels-1]->smoothd,&tpc);CHKERRQ(ierr);
371 ierr = PCPreSolveChangeRHS(tpc,&changed);CHKERRQ(ierr);
372 ierr = KSPGetPC(mglevels[levels-1]->smoothu,&tpc);CHKERRQ(ierr);
373 ierr = PCPreSolveChangeRHS(tpc,&changeu);CHKERRQ(ierr);
374 if (!changeu && !changed) {
375 ierr = VecDestroy(&mglevels[levels-1]->b);CHKERRQ(ierr);
376 mglevels[levels-1]->b = b;
377 } else { /* if the smoother changes the rhs during PreSolve, we cannot use the input vector */
378 if (!mglevels[levels-1]->b) {
379 Vec *vec;
380
381 ierr = KSPCreateVecs(mglevels[levels-1]->smoothd,1,&vec,0,NULL);CHKERRQ(ierr);
382 mglevels[levels-1]->b = *vec;
383 ierr = PetscFree(vec);CHKERRQ(ierr);
384 }
385 ierr = VecCopy(b,mglevels[levels-1]->b);CHKERRQ(ierr);
386 }
387 mglevels[levels-1]->x = x;
388
389 if (mg->am == PC_MG_MULTIPLICATIVE) {
390 ierr = VecSet(x,0.0);CHKERRQ(ierr);
391 for (i=0; i<mg->cyclesperpcapply; i++) {
392 ierr = PCMGMCycle_Private(pc,mglevels+levels-1,NULL);CHKERRQ(ierr);
393 }
394 } else if (mg->am == PC_MG_ADDITIVE) {
395 ierr = PCMGACycle_Private(pc,mglevels);CHKERRQ(ierr);
396 } else if (mg->am == PC_MG_KASKADE) {
397 ierr = PCMGKCycle_Private(pc,mglevels);CHKERRQ(ierr);
398 } else {
399 ierr = PCMGFCycle_Private(pc,mglevels);CHKERRQ(ierr);
400 }
401 if (mg->stageApply) {ierr = PetscLogStagePop();CHKERRQ(ierr);}
402 if (!changeu && !changed) mglevels[levels-1]->b = NULL;
403 PetscFunctionReturn(0);
404 }
405
406
PCSetFromOptions_MG(PetscOptionItems * PetscOptionsObject,PC pc)407 PetscErrorCode PCSetFromOptions_MG(PetscOptionItems *PetscOptionsObject,PC pc)
408 {
409 PetscErrorCode ierr;
410 PetscInt levels,cycles;
411 PetscBool flg, flg2;
412 PC_MG *mg = (PC_MG*)pc->data;
413 PC_MG_Levels **mglevels;
414 PCMGType mgtype;
415 PCMGCycleType mgctype;
416 PCMGGalerkinType gtype;
417
418 PetscFunctionBegin;
419 levels = PetscMax(mg->nlevels,1);
420 ierr = PetscOptionsHead(PetscOptionsObject,"Multigrid options");CHKERRQ(ierr);
421 ierr = PetscOptionsInt("-pc_mg_levels","Number of Levels","PCMGSetLevels",levels,&levels,&flg);CHKERRQ(ierr);
422 if (!flg && !mg->levels && pc->dm) {
423 ierr = DMGetRefineLevel(pc->dm,&levels);CHKERRQ(ierr);
424 levels++;
425 mg->usedmfornumberoflevels = PETSC_TRUE;
426 }
427 ierr = PCMGSetLevels(pc,levels,NULL);CHKERRQ(ierr);
428 mglevels = mg->levels;
429
430 mgctype = (PCMGCycleType) mglevels[0]->cycles;
431 ierr = PetscOptionsEnum("-pc_mg_cycle_type","V cycle or for W-cycle","PCMGSetCycleType",PCMGCycleTypes,(PetscEnum)mgctype,(PetscEnum*)&mgctype,&flg);CHKERRQ(ierr);
432 if (flg) {
433 ierr = PCMGSetCycleType(pc,mgctype);CHKERRQ(ierr);
434 }
435 gtype = mg->galerkin;
436 ierr = PetscOptionsEnum("-pc_mg_galerkin","Use Galerkin process to compute coarser operators","PCMGSetGalerkin",PCMGGalerkinTypes,(PetscEnum)gtype,(PetscEnum*)>ype,&flg);CHKERRQ(ierr);
437 if (flg) {
438 ierr = PCMGSetGalerkin(pc,gtype);CHKERRQ(ierr);
439 }
440 flg2 = PETSC_FALSE;
441 ierr = PetscOptionsBool("-pc_mg_adapt_interp","Adapt interpolation using some coarse space","PCMGSetAdaptInterpolation",PETSC_FALSE,&flg2,&flg);CHKERRQ(ierr);
442 if (flg) {ierr = PCMGSetAdaptInterpolation(pc, flg2);CHKERRQ(ierr);}
443 ierr = PetscOptionsInt("-pc_mg_adapt_interp_n","Size of the coarse space for adaptive interpolation","PCMGSetCoarseSpace",mg->Nc,&mg->Nc,&flg);CHKERRQ(ierr);
444 ierr = 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);CHKERRQ(ierr);
445 ierr = PetscOptionsBool("-pc_mg_mesp_monitor","Monitor the multilevel eigensolver","PCMGSetAdaptInterpolation",PETSC_FALSE,&mg->mespMonitor,&flg);CHKERRQ(ierr);
446 flg = PETSC_FALSE;
447 ierr = PetscOptionsBool("-pc_mg_distinct_smoothup","Create separate smoothup KSP and append the prefix _up","PCMGSetDistinctSmoothUp",PETSC_FALSE,&flg,NULL);CHKERRQ(ierr);
448 if (flg) {
449 ierr = PCMGSetDistinctSmoothUp(pc);CHKERRQ(ierr);
450 }
451 mgtype = mg->am;
452 ierr = PetscOptionsEnum("-pc_mg_type","Multigrid type","PCMGSetType",PCMGTypes,(PetscEnum)mgtype,(PetscEnum*)&mgtype,&flg);CHKERRQ(ierr);
453 if (flg) {
454 ierr = PCMGSetType(pc,mgtype);CHKERRQ(ierr);
455 }
456 if (mg->am == PC_MG_MULTIPLICATIVE) {
457 ierr = PetscOptionsInt("-pc_mg_multiplicative_cycles","Number of cycles for each preconditioner step","PCMGMultiplicativeSetCycles",mg->cyclesperpcapply,&cycles,&flg);CHKERRQ(ierr);
458 if (flg) {
459 ierr = PCMGMultiplicativeSetCycles(pc,cycles);CHKERRQ(ierr);
460 }
461 }
462 flg = PETSC_FALSE;
463 ierr = PetscOptionsBool("-pc_mg_log","Log times for each multigrid level","None",flg,&flg,NULL);CHKERRQ(ierr);
464 if (flg) {
465 PetscInt i;
466 char eventname[128];
467
468 levels = mglevels[0]->levels;
469 for (i=0; i<levels; i++) {
470 sprintf(eventname,"MGSetup Level %d",(int)i);
471 ierr = PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventsmoothsetup);CHKERRQ(ierr);
472 sprintf(eventname,"MGSmooth Level %d",(int)i);
473 ierr = PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventsmoothsolve);CHKERRQ(ierr);
474 if (i) {
475 sprintf(eventname,"MGResid Level %d",(int)i);
476 ierr = PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventresidual);CHKERRQ(ierr);
477 sprintf(eventname,"MGInterp Level %d",(int)i);
478 ierr = PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventinterprestrict);CHKERRQ(ierr);
479 }
480 }
481
482 #if defined(PETSC_USE_LOG)
483 {
484 const char *sname = "MG Apply";
485 PetscStageLog stageLog;
486 PetscInt st;
487
488 ierr = PetscLogGetStageLog(&stageLog);CHKERRQ(ierr);
489 for (st = 0; st < stageLog->numStages; ++st) {
490 PetscBool same;
491
492 ierr = PetscStrcmp(stageLog->stageInfo[st].name, sname, &same);CHKERRQ(ierr);
493 if (same) mg->stageApply = st;
494 }
495 if (!mg->stageApply) {
496 ierr = PetscLogStageRegister(sname, &mg->stageApply);CHKERRQ(ierr);
497 }
498 }
499 #endif
500 }
501 ierr = PetscOptionsTail();CHKERRQ(ierr);
502 /* Check option consistency */
503 ierr = PCMGGetGalerkin(pc, >ype);CHKERRQ(ierr);
504 ierr = PCMGGetAdaptInterpolation(pc, &flg);CHKERRQ(ierr);
505 if (flg && (gtype >= PC_MG_GALERKIN_NONE)) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_INCOMP, "Must use Galerkin coarse operators when adapting the interpolator");
506 PetscFunctionReturn(0);
507 }
508
509 const char *const PCMGTypes[] = {"MULTIPLICATIVE","ADDITIVE","FULL","KASKADE","PCMGType","PC_MG",NULL};
510 const char *const PCMGCycleTypes[] = {"invalid","v","w","PCMGCycleType","PC_MG_CYCLE",NULL};
511 const char *const PCMGGalerkinTypes[] = {"both","pmat","mat","none","external","PCMGGalerkinType","PC_MG_GALERKIN",NULL};
512 const char *const PCMGCoarseSpaceTypes[] = {"polynomial","harmonic","eigenvector","generalized_eigenvector","PCMGCoarseSpaceType","PCMG_POLYNOMIAL",NULL};
513
514 #include <petscdraw.h>
PCView_MG(PC pc,PetscViewer viewer)515 PetscErrorCode PCView_MG(PC pc,PetscViewer viewer)
516 {
517 PC_MG *mg = (PC_MG*)pc->data;
518 PC_MG_Levels **mglevels = mg->levels;
519 PetscErrorCode ierr;
520 PetscInt levels = mglevels ? mglevels[0]->levels : 0,i;
521 PetscBool iascii,isbinary,isdraw;
522
523 PetscFunctionBegin;
524 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
525 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr);
526 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr);
527 if (iascii) {
528 const char *cyclename = levels ? (mglevels[0]->cycles == PC_MG_CYCLE_V ? "v" : "w") : "unknown";
529 ierr = PetscViewerASCIIPrintf(viewer," type is %s, levels=%D cycles=%s\n", PCMGTypes[mg->am],levels,cyclename);CHKERRQ(ierr);
530 if (mg->am == PC_MG_MULTIPLICATIVE) {
531 ierr = PetscViewerASCIIPrintf(viewer," Cycles per PCApply=%d\n",mg->cyclesperpcapply);CHKERRQ(ierr);
532 }
533 if (mg->galerkin == PC_MG_GALERKIN_BOTH) {
534 ierr = PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid matrices\n");CHKERRQ(ierr);
535 } else if (mg->galerkin == PC_MG_GALERKIN_PMAT) {
536 ierr = PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid matrices for pmat\n");CHKERRQ(ierr);
537 } else if (mg->galerkin == PC_MG_GALERKIN_MAT) {
538 ierr = PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid matrices for mat\n");CHKERRQ(ierr);
539 } else if (mg->galerkin == PC_MG_GALERKIN_EXTERNAL) {
540 ierr = PetscViewerASCIIPrintf(viewer," Using externally compute Galerkin coarse grid matrices\n");CHKERRQ(ierr);
541 } else {
542 ierr = PetscViewerASCIIPrintf(viewer," Not using Galerkin computed coarse grid matrices\n");CHKERRQ(ierr);
543 }
544 if (mg->view){
545 ierr = (*mg->view)(pc,viewer);CHKERRQ(ierr);
546 }
547 for (i=0; i<levels; i++) {
548 if (!i) {
549 ierr = PetscViewerASCIIPrintf(viewer,"Coarse grid solver -- level -------------------------------\n",i);CHKERRQ(ierr);
550 } else {
551 ierr = PetscViewerASCIIPrintf(viewer,"Down solver (pre-smoother) on level %D -------------------------------\n",i);CHKERRQ(ierr);
552 }
553 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
554 ierr = KSPView(mglevels[i]->smoothd,viewer);CHKERRQ(ierr);
555 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
556 if (i && mglevels[i]->smoothd == mglevels[i]->smoothu) {
557 ierr = PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) same as down solver (pre-smoother)\n");CHKERRQ(ierr);
558 } else if (i) {
559 ierr = PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) on level %D -------------------------------\n",i);CHKERRQ(ierr);
560 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
561 ierr = KSPView(mglevels[i]->smoothu,viewer);CHKERRQ(ierr);
562 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
563 }
564 }
565 } else if (isbinary) {
566 for (i=levels-1; i>=0; i--) {
567 ierr = KSPView(mglevels[i]->smoothd,viewer);CHKERRQ(ierr);
568 if (i && mglevels[i]->smoothd != mglevels[i]->smoothu) {
569 ierr = KSPView(mglevels[i]->smoothu,viewer);CHKERRQ(ierr);
570 }
571 }
572 } else if (isdraw) {
573 PetscDraw draw;
574 PetscReal x,w,y,bottom,th;
575 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr);
576 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr);
577 ierr = PetscDrawStringGetSize(draw,NULL,&th);CHKERRQ(ierr);
578 bottom = y - th;
579 for (i=levels-1; i>=0; i--) {
580 if (!mglevels[i]->smoothu || (mglevels[i]->smoothu == mglevels[i]->smoothd)) {
581 ierr = PetscDrawPushCurrentPoint(draw,x,bottom);CHKERRQ(ierr);
582 ierr = KSPView(mglevels[i]->smoothd,viewer);CHKERRQ(ierr);
583 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr);
584 } else {
585 w = 0.5*PetscMin(1.0-x,x);
586 ierr = PetscDrawPushCurrentPoint(draw,x+w,bottom);CHKERRQ(ierr);
587 ierr = KSPView(mglevels[i]->smoothd,viewer);CHKERRQ(ierr);
588 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr);
589 ierr = PetscDrawPushCurrentPoint(draw,x-w,bottom);CHKERRQ(ierr);
590 ierr = KSPView(mglevels[i]->smoothu,viewer);CHKERRQ(ierr);
591 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr);
592 }
593 ierr = PetscDrawGetBoundingBox(draw,NULL,&bottom,NULL,NULL);CHKERRQ(ierr);
594 bottom -= th;
595 }
596 }
597 PetscFunctionReturn(0);
598 }
599
600 #include <petsc/private/kspimpl.h>
601
602 /*
603 Calls setup for the KSP on each level
604 */
PCSetUp_MG(PC pc)605 PetscErrorCode PCSetUp_MG(PC pc)
606 {
607 PC_MG *mg = (PC_MG*)pc->data;
608 PC_MG_Levels **mglevels = mg->levels;
609 PetscErrorCode ierr;
610 PetscInt i,n;
611 PC cpc;
612 PetscBool dump = PETSC_FALSE,opsset,use_amat,missinginterpolate = PETSC_FALSE;
613 Mat dA,dB;
614 Vec tvec;
615 DM *dms;
616 PetscViewer viewer = NULL;
617 PetscBool dAeqdB = PETSC_FALSE, needRestricts = PETSC_FALSE;
618
619 PetscFunctionBegin;
620 if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels with PCMGSetLevels() before setting up");
621 n = mglevels[0]->levels;
622 /* FIX: Move this to PCSetFromOptions_MG? */
623 if (mg->usedmfornumberoflevels) {
624 PetscInt levels;
625 ierr = DMGetRefineLevel(pc->dm,&levels);CHKERRQ(ierr);
626 levels++;
627 if (levels > n) { /* the problem is now being solved on a finer grid */
628 ierr = PCMGSetLevels(pc,levels,NULL);CHKERRQ(ierr);
629 n = levels;
630 ierr = PCSetFromOptions(pc);CHKERRQ(ierr); /* it is bad to call this here, but otherwise will never be called for the new hierarchy */
631 mglevels = mg->levels;
632 }
633 }
634 ierr = KSPGetPC(mglevels[0]->smoothd,&cpc);CHKERRQ(ierr);
635
636
637 /* If user did not provide fine grid operators OR operator was not updated since last global KSPSetOperators() */
638 /* so use those from global PC */
639 /* Is this what we always want? What if user wants to keep old one? */
640 ierr = KSPGetOperatorsSet(mglevels[n-1]->smoothd,NULL,&opsset);CHKERRQ(ierr);
641 if (opsset) {
642 Mat mmat;
643 ierr = KSPGetOperators(mglevels[n-1]->smoothd,NULL,&mmat);CHKERRQ(ierr);
644 if (mmat == pc->pmat) opsset = PETSC_FALSE;
645 }
646
647 if (!opsset) {
648 ierr = PCGetUseAmat(pc,&use_amat);CHKERRQ(ierr);
649 if (use_amat) {
650 ierr = PetscInfo(pc,"Using outer operators to define finest grid operator \n because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");CHKERRQ(ierr);
651 ierr = KSPSetOperators(mglevels[n-1]->smoothd,pc->mat,pc->pmat);CHKERRQ(ierr);
652 } else {
653 ierr = PetscInfo(pc,"Using matrix (pmat) operators to define finest grid operator \n because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");CHKERRQ(ierr);
654 ierr = KSPSetOperators(mglevels[n-1]->smoothd,pc->pmat,pc->pmat);CHKERRQ(ierr);
655 }
656 }
657
658 for (i=n-1; i>0; i--) {
659 if (!(mglevels[i]->interpolate || mglevels[i]->restrct)) {
660 missinginterpolate = PETSC_TRUE;
661 continue;
662 }
663 }
664
665 ierr = KSPGetOperators(mglevels[n-1]->smoothd,&dA,&dB);CHKERRQ(ierr);
666 if (dA == dB) dAeqdB = PETSC_TRUE;
667 if ((mg->galerkin == PC_MG_GALERKIN_NONE) || (((mg->galerkin == PC_MG_GALERKIN_PMAT) || (mg->galerkin == PC_MG_GALERKIN_MAT)) && !dAeqdB)) {
668 needRestricts = PETSC_TRUE; /* user must compute either mat, pmat, or both so must restrict x to coarser levels */
669 }
670
671
672 /*
673 Skipping if user has provided all interpolation/restriction needed (since DM might not be able to produce them (when coming from SNES/TS)
674 Skipping for galerkin==2 (externally managed hierarchy such as ML and GAMG). Cleaner logic here would be great. Wrap ML/GAMG as DMs?
675 */
676 if (missinginterpolate && pc->dm && mg->galerkin != PC_MG_GALERKIN_EXTERNAL && !pc->setupcalled) {
677 /* construct the interpolation from the DMs */
678 Mat p;
679 Vec rscale;
680 ierr = PetscMalloc1(n,&dms);CHKERRQ(ierr);
681 dms[n-1] = pc->dm;
682 /* Separately create them so we do not get DMKSP interference between levels */
683 for (i=n-2; i>-1; i--) {ierr = DMCoarsen(dms[i+1],MPI_COMM_NULL,&dms[i]);CHKERRQ(ierr);}
684 /*
685 Force the mat type of coarse level operator to be AIJ because usually we want to use LU for coarse level.
686 Notice that it can be overwritten by -mat_type because KSPSetUp() reads command line options.
687 But it is safe to use -dm_mat_type.
688
689 The mat type should not be hardcoded like this, we need to find a better way.
690 ierr = DMSetMatType(dms[0],MATAIJ);CHKERRQ(ierr);
691 */
692 for (i=n-2; i>-1; i--) {
693 DMKSP kdm;
694 PetscBool dmhasrestrict, dmhasinject;
695 ierr = KSPSetDM(mglevels[i]->smoothd,dms[i]);CHKERRQ(ierr);
696 if (!needRestricts) {ierr = KSPSetDMActive(mglevels[i]->smoothd,PETSC_FALSE);CHKERRQ(ierr);}
697 if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
698 ierr = KSPSetDM(mglevels[i]->smoothu,dms[i]);CHKERRQ(ierr);
699 if (!needRestricts) {ierr = KSPSetDMActive(mglevels[i]->smoothu,PETSC_FALSE);CHKERRQ(ierr);}
700 }
701 ierr = DMGetDMKSPWrite(dms[i],&kdm);CHKERRQ(ierr);
702 /* Ugly hack so that the next KSPSetUp() will use the RHS that we set. A better fix is to change dmActive to take
703 * a bitwise OR of computing the matrix, RHS, and initial iterate. */
704 kdm->ops->computerhs = NULL;
705 kdm->rhsctx = NULL;
706 if (!mglevels[i+1]->interpolate) {
707 ierr = DMCreateInterpolation(dms[i],dms[i+1],&p,&rscale);CHKERRQ(ierr);
708 ierr = PCMGSetInterpolation(pc,i+1,p);CHKERRQ(ierr);
709 if (rscale) {ierr = PCMGSetRScale(pc,i+1,rscale);CHKERRQ(ierr);}
710 ierr = VecDestroy(&rscale);CHKERRQ(ierr);
711 ierr = MatDestroy(&p);CHKERRQ(ierr);
712 }
713 ierr = DMHasCreateRestriction(dms[i],&dmhasrestrict);CHKERRQ(ierr);
714 if (dmhasrestrict && !mglevels[i+1]->restrct){
715 ierr = DMCreateRestriction(dms[i],dms[i+1],&p);CHKERRQ(ierr);
716 ierr = PCMGSetRestriction(pc,i+1,p);CHKERRQ(ierr);
717 ierr = MatDestroy(&p);CHKERRQ(ierr);
718 }
719 ierr = DMHasCreateInjection(dms[i],&dmhasinject);CHKERRQ(ierr);
720 if (dmhasinject && !mglevels[i+1]->inject){
721 ierr = DMCreateInjection(dms[i],dms[i+1],&p);CHKERRQ(ierr);
722 ierr = PCMGSetInjection(pc,i+1,p);CHKERRQ(ierr);
723 ierr = MatDestroy(&p);CHKERRQ(ierr);
724 }
725 }
726
727 for (i=n-2; i>-1; i--) {ierr = DMDestroy(&dms[i]);CHKERRQ(ierr);}
728 ierr = PetscFree(dms);CHKERRQ(ierr);
729 }
730
731 if (pc->dm && !pc->setupcalled) {
732 /* finest smoother also gets DM but it is not active, independent of whether galerkin==PC_MG_GALERKIN_EXTERNAL */
733 ierr = KSPSetDM(mglevels[n-1]->smoothd,pc->dm);CHKERRQ(ierr);
734 ierr = KSPSetDMActive(mglevels[n-1]->smoothd,PETSC_FALSE);CHKERRQ(ierr);
735 if (mglevels[n-1]->smoothd != mglevels[n-1]->smoothu) {
736 ierr = KSPSetDM(mglevels[n-1]->smoothu,pc->dm);CHKERRQ(ierr);
737 ierr = KSPSetDMActive(mglevels[n-1]->smoothu,PETSC_FALSE);CHKERRQ(ierr);
738 }
739 }
740
741 if (mg->galerkin < PC_MG_GALERKIN_NONE) {
742 Mat A,B;
743 PetscBool doA = PETSC_FALSE,doB = PETSC_FALSE;
744 MatReuse reuse = MAT_INITIAL_MATRIX;
745
746 if ((mg->galerkin == PC_MG_GALERKIN_PMAT) || (mg->galerkin == PC_MG_GALERKIN_BOTH)) doB = PETSC_TRUE;
747 if ((mg->galerkin == PC_MG_GALERKIN_MAT) || ((mg->galerkin == PC_MG_GALERKIN_BOTH) && (dA != dB))) doA = PETSC_TRUE;
748 if (pc->setupcalled) reuse = MAT_REUSE_MATRIX;
749 for (i=n-2; i>-1; i--) {
750 if (!mglevels[i+1]->restrct && !mglevels[i+1]->interpolate) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must provide interpolation or restriction for each MG level except level 0");
751 if (!mglevels[i+1]->interpolate) {
752 ierr = PCMGSetInterpolation(pc,i+1,mglevels[i+1]->restrct);CHKERRQ(ierr);
753 }
754 if (!mglevels[i+1]->restrct) {
755 ierr = PCMGSetRestriction(pc,i+1,mglevels[i+1]->interpolate);CHKERRQ(ierr);
756 }
757 if (reuse == MAT_REUSE_MATRIX) {
758 ierr = KSPGetOperators(mglevels[i]->smoothd,&A,&B);CHKERRQ(ierr);
759 }
760 if (doA) {
761 ierr = MatGalerkin(mglevels[i+1]->restrct,dA,mglevels[i+1]->interpolate,reuse,1.0,&A);CHKERRQ(ierr);
762 }
763 if (doB) {
764 ierr = MatGalerkin(mglevels[i+1]->restrct,dB,mglevels[i+1]->interpolate,reuse,1.0,&B);CHKERRQ(ierr);
765 }
766 /* the management of the PetscObjectReference() and PetscObjecDereference() below is rather delicate */
767 if (!doA && dAeqdB) {
768 if (reuse == MAT_INITIAL_MATRIX) {ierr = PetscObjectReference((PetscObject)B);CHKERRQ(ierr);}
769 A = B;
770 } else if (!doA && reuse == MAT_INITIAL_MATRIX) {
771 ierr = KSPGetOperators(mglevels[i]->smoothd,&A,NULL);CHKERRQ(ierr);
772 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
773 }
774 if (!doB && dAeqdB) {
775 if (reuse == MAT_INITIAL_MATRIX) {ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);}
776 B = A;
777 } else if (!doB && reuse == MAT_INITIAL_MATRIX) {
778 ierr = KSPGetOperators(mglevels[i]->smoothd,NULL,&B);CHKERRQ(ierr);
779 ierr = PetscObjectReference((PetscObject)B);CHKERRQ(ierr);
780 }
781 if (reuse == MAT_INITIAL_MATRIX) {
782 ierr = KSPSetOperators(mglevels[i]->smoothd,A,B);CHKERRQ(ierr);
783 ierr = PetscObjectDereference((PetscObject)A);CHKERRQ(ierr);
784 ierr = PetscObjectDereference((PetscObject)B);CHKERRQ(ierr);
785 }
786 dA = A;
787 dB = B;
788 }
789 }
790
791
792 /* Adapt interpolation matrices */
793 if (mg->adaptInterpolation) {
794 mg->Nc = mg->Nc < 0 ? 6 : mg->Nc; /* Default to 6 modes */
795 for (i = 0; i < n; ++i) {
796 ierr = PCMGComputeCoarseSpace_Internal(pc, i, mg->coarseSpaceType, mg->Nc, !i ? NULL : mglevels[i-1]->coarseSpace, &mglevels[i]->coarseSpace);CHKERRQ(ierr);
797 if (i) {ierr = PCMGAdaptInterpolator_Internal(pc, i, mglevels[i-1]->smoothu, mglevels[i]->smoothu, mg->Nc, mglevels[i-1]->coarseSpace, mglevels[i]->coarseSpace);CHKERRQ(ierr);}
798 }
799 for (i = n-2; i > -1; --i) {
800 ierr = PCMGRecomputeLevelOperators_Internal(pc, i);CHKERRQ(ierr);
801 }
802 }
803
804 if (needRestricts && pc->dm) {
805 for (i=n-2; i>=0; i--) {
806 DM dmfine,dmcoarse;
807 Mat Restrict,Inject;
808 Vec rscale;
809 ierr = KSPGetDM(mglevels[i+1]->smoothd,&dmfine);CHKERRQ(ierr);
810 ierr = KSPGetDM(mglevels[i]->smoothd,&dmcoarse);CHKERRQ(ierr);
811 ierr = PCMGGetRestriction(pc,i+1,&Restrict);CHKERRQ(ierr);
812 ierr = PCMGGetRScale(pc,i+1,&rscale);CHKERRQ(ierr);
813 ierr = PCMGGetInjection(pc,i+1,&Inject);CHKERRQ(ierr);
814 ierr = DMRestrict(dmfine,Restrict,rscale,Inject,dmcoarse);CHKERRQ(ierr);
815 }
816 }
817
818 if (!pc->setupcalled) {
819 for (i=0; i<n; i++) {
820 ierr = KSPSetFromOptions(mglevels[i]->smoothd);CHKERRQ(ierr);
821 }
822 for (i=1; i<n; i++) {
823 if (mglevels[i]->smoothu && (mglevels[i]->smoothu != mglevels[i]->smoothd)) {
824 ierr = KSPSetFromOptions(mglevels[i]->smoothu);CHKERRQ(ierr);
825 }
826 }
827 /* insure that if either interpolation or restriction is set the other other one is set */
828 for (i=1; i<n; i++) {
829 ierr = PCMGGetInterpolation(pc,i,NULL);CHKERRQ(ierr);
830 ierr = PCMGGetRestriction(pc,i,NULL);CHKERRQ(ierr);
831 }
832 for (i=0; i<n-1; i++) {
833 if (!mglevels[i]->b) {
834 Vec *vec;
835 ierr = KSPCreateVecs(mglevels[i]->smoothd,1,&vec,0,NULL);CHKERRQ(ierr);
836 ierr = PCMGSetRhs(pc,i,*vec);CHKERRQ(ierr);
837 ierr = VecDestroy(vec);CHKERRQ(ierr);
838 ierr = PetscFree(vec);CHKERRQ(ierr);
839 }
840 if (!mglevels[i]->r && i) {
841 ierr = VecDuplicate(mglevels[i]->b,&tvec);CHKERRQ(ierr);
842 ierr = PCMGSetR(pc,i,tvec);CHKERRQ(ierr);
843 ierr = VecDestroy(&tvec);CHKERRQ(ierr);
844 }
845 if (!mglevels[i]->x) {
846 ierr = VecDuplicate(mglevels[i]->b,&tvec);CHKERRQ(ierr);
847 ierr = PCMGSetX(pc,i,tvec);CHKERRQ(ierr);
848 ierr = VecDestroy(&tvec);CHKERRQ(ierr);
849 }
850 }
851 if (n != 1 && !mglevels[n-1]->r) {
852 /* PCMGSetR() on the finest level if user did not supply it */
853 Vec *vec;
854 ierr = KSPCreateVecs(mglevels[n-1]->smoothd,1,&vec,0,NULL);CHKERRQ(ierr);
855 ierr = PCMGSetR(pc,n-1,*vec);CHKERRQ(ierr);
856 ierr = VecDestroy(vec);CHKERRQ(ierr);
857 ierr = PetscFree(vec);CHKERRQ(ierr);
858 }
859 }
860
861 if (pc->dm) {
862 /* need to tell all the coarser levels to rebuild the matrix using the DM for that level */
863 for (i=0; i<n-1; i++) {
864 if (mglevels[i]->smoothd->setupstage != KSP_SETUP_NEW) mglevels[i]->smoothd->setupstage = KSP_SETUP_NEWMATRIX;
865 }
866 }
867
868 for (i=1; i<n; i++) {
869 if (mglevels[i]->smoothu == mglevels[i]->smoothd || mg->am == PC_MG_FULL || mg->am == PC_MG_KASKADE || mg->cyclesperpcapply > 1){
870 /* if doing only down then initial guess is zero */
871 ierr = KSPSetInitialGuessNonzero(mglevels[i]->smoothd,PETSC_TRUE);CHKERRQ(ierr);
872 }
873 if (mglevels[i]->eventsmoothsetup) {ierr = PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);}
874 ierr = KSPSetUp(mglevels[i]->smoothd);CHKERRQ(ierr);
875 if (mglevels[i]->smoothd->reason == KSP_DIVERGED_PC_FAILED) {
876 pc->failedreason = PC_SUBPC_ERROR;
877 }
878 if (mglevels[i]->eventsmoothsetup) {ierr = PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);}
879 if (!mglevels[i]->residual) {
880 Mat mat;
881 ierr = KSPGetOperators(mglevels[i]->smoothd,&mat,NULL);CHKERRQ(ierr);
882 ierr = PCMGSetResidual(pc,i,PCMGResidualDefault,mat);CHKERRQ(ierr);
883 }
884 }
885 for (i=1; i<n; i++) {
886 if (mglevels[i]->smoothu && mglevels[i]->smoothu != mglevels[i]->smoothd) {
887 Mat downmat,downpmat;
888
889 /* check if operators have been set for up, if not use down operators to set them */
890 ierr = KSPGetOperatorsSet(mglevels[i]->smoothu,&opsset,NULL);CHKERRQ(ierr);
891 if (!opsset) {
892 ierr = KSPGetOperators(mglevels[i]->smoothd,&downmat,&downpmat);CHKERRQ(ierr);
893 ierr = KSPSetOperators(mglevels[i]->smoothu,downmat,downpmat);CHKERRQ(ierr);
894 }
895
896 ierr = KSPSetInitialGuessNonzero(mglevels[i]->smoothu,PETSC_TRUE);CHKERRQ(ierr);
897 if (mglevels[i]->eventsmoothsetup) {ierr = PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);}
898 ierr = KSPSetUp(mglevels[i]->smoothu);CHKERRQ(ierr);
899 if (mglevels[i]->smoothu->reason == KSP_DIVERGED_PC_FAILED) {
900 pc->failedreason = PC_SUBPC_ERROR;
901 }
902 if (mglevels[i]->eventsmoothsetup) {ierr = PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);}
903 }
904 }
905
906 if (mglevels[0]->eventsmoothsetup) {ierr = PetscLogEventBegin(mglevels[0]->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);}
907 ierr = KSPSetUp(mglevels[0]->smoothd);CHKERRQ(ierr);
908 if (mglevels[0]->smoothd->reason == KSP_DIVERGED_PC_FAILED) {
909 pc->failedreason = PC_SUBPC_ERROR;
910 }
911 if (mglevels[0]->eventsmoothsetup) {ierr = PetscLogEventEnd(mglevels[0]->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);}
912
913 /*
914 Dump the interpolation/restriction matrices plus the
915 Jacobian/stiffness on each level. This allows MATLAB users to
916 easily check if the Galerkin condition A_c = R A_f R^T is satisfied.
917
918 Only support one or the other at the same time.
919 */
920 #if defined(PETSC_USE_SOCKET_VIEWER)
921 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_mg_dump_matlab",&dump,NULL);CHKERRQ(ierr);
922 if (dump) viewer = PETSC_VIEWER_SOCKET_(PetscObjectComm((PetscObject)pc));
923 dump = PETSC_FALSE;
924 #endif
925 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_mg_dump_binary",&dump,NULL);CHKERRQ(ierr);
926 if (dump) viewer = PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)pc));
927
928 if (viewer) {
929 for (i=1; i<n; i++) {
930 ierr = MatView(mglevels[i]->restrct,viewer);CHKERRQ(ierr);
931 }
932 for (i=0; i<n; i++) {
933 ierr = KSPGetPC(mglevels[i]->smoothd,&pc);CHKERRQ(ierr);
934 ierr = MatView(pc->mat,viewer);CHKERRQ(ierr);
935 }
936 }
937 PetscFunctionReturn(0);
938 }
939
940 /* -------------------------------------------------------------------------------------*/
941
PCMGGetLevels_MG(PC pc,PetscInt * levels)942 PetscErrorCode PCMGGetLevels_MG(PC pc, PetscInt *levels)
943 {
944 PC_MG *mg = (PC_MG *) pc->data;
945
946 PetscFunctionBegin;
947 *levels = mg->nlevels;
948 PetscFunctionReturn(0);
949 }
950
951 /*@
952 PCMGGetLevels - Gets the number of levels to use with MG.
953
954 Not Collective
955
956 Input Parameter:
957 . pc - the preconditioner context
958
959 Output parameter:
960 . levels - the number of levels
961
962 Level: advanced
963
964 .seealso: PCMGSetLevels()
965 @*/
PCMGGetLevels(PC pc,PetscInt * levels)966 PetscErrorCode PCMGGetLevels(PC pc,PetscInt *levels)
967 {
968 PetscErrorCode ierr;
969
970 PetscFunctionBegin;
971 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
972 PetscValidIntPointer(levels,2);
973 *levels = 0;
974 ierr = PetscTryMethod(pc,"PCMGGetLevels_C",(PC,PetscInt*),(pc,levels));CHKERRQ(ierr);
975 PetscFunctionReturn(0);
976 }
977
978 /*@
979 PCMGSetType - Determines the form of multigrid to use:
980 multiplicative, additive, full, or the Kaskade algorithm.
981
982 Logically Collective on PC
983
984 Input Parameters:
985 + pc - the preconditioner context
986 - form - multigrid form, one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,
987 PC_MG_FULL, PC_MG_KASKADE
988
989 Options Database Key:
990 . -pc_mg_type <form> - Sets <form>, one of multiplicative,
991 additive, full, kaskade
992
993 Level: advanced
994
995 .seealso: PCMGSetLevels()
996 @*/
PCMGSetType(PC pc,PCMGType form)997 PetscErrorCode PCMGSetType(PC pc,PCMGType form)
998 {
999 PC_MG *mg = (PC_MG*)pc->data;
1000
1001 PetscFunctionBegin;
1002 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1003 PetscValidLogicalCollectiveEnum(pc,form,2);
1004 mg->am = form;
1005 if (form == PC_MG_MULTIPLICATIVE) pc->ops->applyrichardson = PCApplyRichardson_MG;
1006 else pc->ops->applyrichardson = NULL;
1007 PetscFunctionReturn(0);
1008 }
1009
1010 /*@
1011 PCMGGetType - Determines the form of multigrid to use:
1012 multiplicative, additive, full, or the Kaskade algorithm.
1013
1014 Logically Collective on PC
1015
1016 Input Parameter:
1017 . pc - the preconditioner context
1018
1019 Output Parameter:
1020 . type - one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,PC_MG_FULL, PC_MG_KASKADE
1021
1022
1023 Level: advanced
1024
1025 .seealso: PCMGSetLevels()
1026 @*/
PCMGGetType(PC pc,PCMGType * type)1027 PetscErrorCode PCMGGetType(PC pc,PCMGType *type)
1028 {
1029 PC_MG *mg = (PC_MG*)pc->data;
1030
1031 PetscFunctionBegin;
1032 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1033 *type = mg->am;
1034 PetscFunctionReturn(0);
1035 }
1036
1037 /*@
1038 PCMGSetCycleType - Sets the type cycles to use. Use PCMGSetCycleTypeOnLevel() for more
1039 complicated cycling.
1040
1041 Logically Collective on PC
1042
1043 Input Parameters:
1044 + pc - the multigrid context
1045 - n - either PC_MG_CYCLE_V or PC_MG_CYCLE_W
1046
1047 Options Database Key:
1048 . -pc_mg_cycle_type <v,w> - provide the cycle desired
1049
1050 Level: advanced
1051
1052 .seealso: PCMGSetCycleTypeOnLevel()
1053 @*/
PCMGSetCycleType(PC pc,PCMGCycleType n)1054 PetscErrorCode PCMGSetCycleType(PC pc,PCMGCycleType n)
1055 {
1056 PC_MG *mg = (PC_MG*)pc->data;
1057 PC_MG_Levels **mglevels = mg->levels;
1058 PetscInt i,levels;
1059
1060 PetscFunctionBegin;
1061 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1062 PetscValidLogicalCollectiveEnum(pc,n,2);
1063 if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Must set MG levels with PCMGSetLevels() before calling");
1064 levels = mglevels[0]->levels;
1065 for (i=0; i<levels; i++) mglevels[i]->cycles = n;
1066 PetscFunctionReturn(0);
1067 }
1068
1069 /*@
1070 PCMGMultiplicativeSetCycles - Sets the number of cycles to use for each preconditioner step
1071 of multigrid when PCMGType of PC_MG_MULTIPLICATIVE is used
1072
1073 Logically Collective on PC
1074
1075 Input Parameters:
1076 + pc - the multigrid context
1077 - n - number of cycles (default is 1)
1078
1079 Options Database Key:
1080 . -pc_mg_multiplicative_cycles n
1081
1082 Level: advanced
1083
1084 Notes:
1085 This is not associated with setting a v or w cycle, that is set with PCMGSetCycleType()
1086
1087 .seealso: PCMGSetCycleTypeOnLevel(), PCMGSetCycleType()
1088 @*/
PCMGMultiplicativeSetCycles(PC pc,PetscInt n)1089 PetscErrorCode PCMGMultiplicativeSetCycles(PC pc,PetscInt n)
1090 {
1091 PC_MG *mg = (PC_MG*)pc->data;
1092
1093 PetscFunctionBegin;
1094 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1095 PetscValidLogicalCollectiveInt(pc,n,2);
1096 mg->cyclesperpcapply = n;
1097 PetscFunctionReturn(0);
1098 }
1099
PCMGSetGalerkin_MG(PC pc,PCMGGalerkinType use)1100 PetscErrorCode PCMGSetGalerkin_MG(PC pc,PCMGGalerkinType use)
1101 {
1102 PC_MG *mg = (PC_MG*)pc->data;
1103
1104 PetscFunctionBegin;
1105 mg->galerkin = use;
1106 PetscFunctionReturn(0);
1107 }
1108
1109 /*@
1110 PCMGSetGalerkin - Causes the coarser grid matrices to be computed from the
1111 finest grid via the Galerkin process: A_i-1 = r_i * A_i * p_i
1112
1113 Logically Collective on PC
1114
1115 Input Parameters:
1116 + pc - the multigrid context
1117 - use - one of PC_MG_GALERKIN_BOTH,PC_MG_GALERKIN_PMAT,PC_MG_GALERKIN_MAT, or PC_MG_GALERKIN_NONE
1118
1119 Options Database Key:
1120 . -pc_mg_galerkin <both,pmat,mat,none>
1121
1122 Level: intermediate
1123
1124 Notes:
1125 Some codes that use PCMG such as PCGAMG use Galerkin internally while constructing the hierarchy and thus do not
1126 use the PCMG construction of the coarser grids.
1127
1128 .seealso: PCMGGetGalerkin(), PCMGGalerkinType
1129
1130 @*/
PCMGSetGalerkin(PC pc,PCMGGalerkinType use)1131 PetscErrorCode PCMGSetGalerkin(PC pc,PCMGGalerkinType use)
1132 {
1133 PetscErrorCode ierr;
1134
1135 PetscFunctionBegin;
1136 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1137 ierr = PetscTryMethod(pc,"PCMGSetGalerkin_C",(PC,PCMGGalerkinType),(pc,use));CHKERRQ(ierr);
1138 PetscFunctionReturn(0);
1139 }
1140
1141 /*@
1142 PCMGGetGalerkin - Checks if Galerkin multigrid is being used, i.e.
1143 A_i-1 = r_i * A_i * p_i
1144
1145 Not Collective
1146
1147 Input Parameter:
1148 . pc - the multigrid context
1149
1150 Output Parameter:
1151 . galerkin - one of PC_MG_GALERKIN_BOTH,PC_MG_GALERKIN_PMAT,PC_MG_GALERKIN_MAT, PC_MG_GALERKIN_NONE, or PC_MG_GALERKIN_EXTERNAL
1152
1153 Level: intermediate
1154
1155 .seealso: PCMGSetGalerkin(), PCMGGalerkinType
1156
1157 @*/
PCMGGetGalerkin(PC pc,PCMGGalerkinType * galerkin)1158 PetscErrorCode PCMGGetGalerkin(PC pc,PCMGGalerkinType *galerkin)
1159 {
1160 PC_MG *mg = (PC_MG*)pc->data;
1161
1162 PetscFunctionBegin;
1163 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1164 *galerkin = mg->galerkin;
1165 PetscFunctionReturn(0);
1166 }
1167
PCMGSetAdaptInterpolation_MG(PC pc,PetscBool adapt)1168 PetscErrorCode PCMGSetAdaptInterpolation_MG(PC pc, PetscBool adapt)
1169 {
1170 PC_MG *mg = (PC_MG *) pc->data;
1171
1172 PetscFunctionBegin;
1173 mg->adaptInterpolation = adapt;
1174 PetscFunctionReturn(0);
1175 }
1176
PCMGGetAdaptInterpolation_MG(PC pc,PetscBool * adapt)1177 PetscErrorCode PCMGGetAdaptInterpolation_MG(PC pc, PetscBool *adapt)
1178 {
1179 PC_MG *mg = (PC_MG *) pc->data;
1180
1181 PetscFunctionBegin;
1182 *adapt = mg->adaptInterpolation;
1183 PetscFunctionReturn(0);
1184 }
1185
1186 /*@
1187 PCMGSetAdaptInterpolation - Adapt the interpolator based upon a vector space which should be accurately captured by the next coarser mesh, and thus accurately interpolated.
1188
1189 Logically Collective on PC
1190
1191 Input Parameters:
1192 + pc - the multigrid context
1193 - adapt - flag for adaptation of the interpolator
1194
1195 Options Database Keys:
1196 + -pc_mg_adapt_interp - Turn on adaptation
1197 . -pc_mg_adapt_interp_n <int> - The number of modes to use, should be divisible by dimension
1198 - -pc_mg_adapt_interp_coarse_space <type> - The type of coarse space: polynomial, harmonic, eigenvector, generalized_eigenvector
1199
1200 Level: intermediate
1201
1202 .keywords: MG, set, Galerkin
1203 .seealso: PCMGGetAdaptInterpolation(), PCMGSetGalerkin()
1204 @*/
PCMGSetAdaptInterpolation(PC pc,PetscBool adapt)1205 PetscErrorCode PCMGSetAdaptInterpolation(PC pc, PetscBool adapt)
1206 {
1207 PetscErrorCode ierr;
1208
1209 PetscFunctionBegin;
1210 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
1211 ierr = PetscTryMethod(pc,"PCMGSetAdaptInterpolation_C",(PC,PetscBool),(pc,adapt));CHKERRQ(ierr);
1212 PetscFunctionReturn(0);
1213 }
1214
1215 /*@
1216 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.
1217
1218 Not collective
1219
1220 Input Parameter:
1221 . pc - the multigrid context
1222
1223 Output Parameter:
1224 . adapt - flag for adaptation of the interpolator
1225
1226 Level: intermediate
1227
1228 .keywords: MG, set, Galerkin
1229 .seealso: PCMGSetAdaptInterpolation(), PCMGSetGalerkin()
1230 @*/
PCMGGetAdaptInterpolation(PC pc,PetscBool * adapt)1231 PetscErrorCode PCMGGetAdaptInterpolation(PC pc, PetscBool *adapt)
1232 {
1233 PetscErrorCode ierr;
1234
1235 PetscFunctionBegin;
1236 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
1237 PetscValidPointer(adapt, 2);
1238 ierr = PetscUseMethod(pc,"PCMGGetAdaptInterpolation_C",(PC,PetscBool*),(pc,adapt));CHKERRQ(ierr);
1239 PetscFunctionReturn(0);
1240 }
1241
1242 /*@
1243 PCMGSetNumberSmooth - Sets the number of pre and post-smoothing steps to use
1244 on all levels. Use PCMGDistinctSmoothUp() to create separate up and down smoothers if you want different numbers of
1245 pre- and post-smoothing steps.
1246
1247 Logically Collective on PC
1248
1249 Input Parameters:
1250 + mg - the multigrid context
1251 - n - the number of smoothing steps
1252
1253 Options Database Key:
1254 . -mg_levels_ksp_max_it <n> - Sets number of pre and post-smoothing steps
1255
1256 Level: advanced
1257
1258 Notes:
1259 this does not set a value on the coarsest grid, since we assume that
1260 there is no separate smooth up on the coarsest grid.
1261
1262 .seealso: PCMGSetDistinctSmoothUp()
1263 @*/
PCMGSetNumberSmooth(PC pc,PetscInt n)1264 PetscErrorCode PCMGSetNumberSmooth(PC pc,PetscInt n)
1265 {
1266 PC_MG *mg = (PC_MG*)pc->data;
1267 PC_MG_Levels **mglevels = mg->levels;
1268 PetscErrorCode ierr;
1269 PetscInt i,levels;
1270
1271 PetscFunctionBegin;
1272 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1273 PetscValidLogicalCollectiveInt(pc,n,2);
1274 if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Must set MG levels with PCMGSetLevels() before calling");
1275 levels = mglevels[0]->levels;
1276
1277 for (i=1; i<levels; i++) {
1278 ierr = KSPSetTolerances(mglevels[i]->smoothu,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);CHKERRQ(ierr);
1279 ierr = KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);CHKERRQ(ierr);
1280 mg->default_smoothu = n;
1281 mg->default_smoothd = n;
1282 }
1283 PetscFunctionReturn(0);
1284 }
1285
1286 /*@
1287 PCMGSetDistinctSmoothUp - sets the up (post) smoother to be a separate KSP from the down (pre) smoother on all levels
1288 and adds the suffix _up to the options name
1289
1290 Logically Collective on PC
1291
1292 Input Parameters:
1293 . pc - the preconditioner context
1294
1295 Options Database Key:
1296 . -pc_mg_distinct_smoothup
1297
1298 Level: advanced
1299
1300 Notes:
1301 this does not set a value on the coarsest grid, since we assume that
1302 there is no separate smooth up on the coarsest grid.
1303
1304 .seealso: PCMGSetNumberSmooth()
1305 @*/
PCMGSetDistinctSmoothUp(PC pc)1306 PetscErrorCode PCMGSetDistinctSmoothUp(PC pc)
1307 {
1308 PC_MG *mg = (PC_MG*)pc->data;
1309 PC_MG_Levels **mglevels = mg->levels;
1310 PetscErrorCode ierr;
1311 PetscInt i,levels;
1312 KSP subksp;
1313
1314 PetscFunctionBegin;
1315 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1316 if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Must set MG levels with PCMGSetLevels() before calling");
1317 levels = mglevels[0]->levels;
1318
1319 for (i=1; i<levels; i++) {
1320 const char *prefix = NULL;
1321 /* make sure smoother up and down are different */
1322 ierr = PCMGGetSmootherUp(pc,i,&subksp);CHKERRQ(ierr);
1323 ierr = KSPGetOptionsPrefix(mglevels[i]->smoothd,&prefix);CHKERRQ(ierr);
1324 ierr = KSPSetOptionsPrefix(subksp,prefix);CHKERRQ(ierr);
1325 ierr = KSPAppendOptionsPrefix(subksp,"up_");CHKERRQ(ierr);
1326 }
1327 PetscFunctionReturn(0);
1328 }
1329
1330 /* No new matrices are created, and the coarse operator matrices are the references to the original ones */
PCGetInterpolations_MG(PC pc,PetscInt * num_levels,Mat * interpolations[])1331 PetscErrorCode PCGetInterpolations_MG(PC pc,PetscInt *num_levels,Mat *interpolations[])
1332 {
1333 PC_MG *mg = (PC_MG*)pc->data;
1334 PC_MG_Levels **mglevels = mg->levels;
1335 Mat *mat;
1336 PetscInt l;
1337 PetscErrorCode ierr;
1338
1339 PetscFunctionBegin;
1340 if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
1341 ierr = PetscMalloc1(mg->nlevels,&mat);CHKERRQ(ierr);
1342 for (l=1; l< mg->nlevels; l++) {
1343 mat[l-1] = mglevels[l]->interpolate;
1344 ierr = PetscObjectReference((PetscObject)mat[l-1]);CHKERRQ(ierr);
1345 }
1346 *num_levels = mg->nlevels;
1347 *interpolations = mat;
1348 PetscFunctionReturn(0);
1349 }
1350
1351 /* No new matrices are created, and the coarse operator matrices are the references to the original ones */
PCGetCoarseOperators_MG(PC pc,PetscInt * num_levels,Mat * coarseOperators[])1352 PetscErrorCode PCGetCoarseOperators_MG(PC pc,PetscInt *num_levels,Mat *coarseOperators[])
1353 {
1354 PC_MG *mg = (PC_MG*)pc->data;
1355 PC_MG_Levels **mglevels = mg->levels;
1356 PetscInt l;
1357 Mat *mat;
1358 PetscErrorCode ierr;
1359
1360 PetscFunctionBegin;
1361 if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
1362 ierr = PetscMalloc1(mg->nlevels,&mat);CHKERRQ(ierr);
1363 for (l=0; l<mg->nlevels-1; l++) {
1364 ierr = KSPGetOperators(mglevels[l]->smoothd,NULL,&(mat[l]));CHKERRQ(ierr);
1365 ierr = PetscObjectReference((PetscObject)mat[l]);CHKERRQ(ierr);
1366 }
1367 *num_levels = mg->nlevels;
1368 *coarseOperators = mat;
1369 PetscFunctionReturn(0);
1370 }
1371
1372 /*@C
1373 PCMGRegisterCoarseSpaceConstructor - Adds a method to the PCMG package for coarse space construction.
1374
1375 Not collective
1376
1377 Input Parameters:
1378 + name - name of the constructor
1379 - function - constructor routine
1380
1381 Notes:
1382 Calling sequence for the routine:
1383 $ my_csp(PC pc, PetscInt l, DM dm, KSP smooth, PetscInt Nc, const Vec initGuess[], Vec **coarseSp)
1384 $ pc - The PC object
1385 $ l - The multigrid level, 0 is the coarse level
1386 $ dm - The DM for this level
1387 $ smooth - The level smoother
1388 $ Nc - The size of the coarse space
1389 $ initGuess - Basis for an initial guess for the space
1390 $ coarseSp - A basis for the computed coarse space
1391
1392 Level: advanced
1393
1394 .seealso: PCMGGetCoarseSpaceConstructor(), PCRegister()
1395 @*/
PCMGRegisterCoarseSpaceConstructor(const char name[],PetscErrorCode (* function)(PC,PetscInt,DM,KSP,PetscInt,const Vec[],Vec **))1396 PetscErrorCode PCMGRegisterCoarseSpaceConstructor(const char name[], PetscErrorCode (*function)(PC, PetscInt, DM, KSP, PetscInt, const Vec[], Vec **))
1397 {
1398 PetscErrorCode ierr;
1399
1400 PetscFunctionBegin;
1401 ierr = PCInitializePackage();CHKERRQ(ierr);
1402 ierr = PetscFunctionListAdd(&PCMGCoarseList,name,function);CHKERRQ(ierr);
1403 PetscFunctionReturn(0);
1404 }
1405
1406 /*@C
1407 PCMGGetCoarseSpaceConstructor - Returns the given coarse space construction method.
1408
1409 Not collective
1410
1411 Input Parameter:
1412 . name - name of the constructor
1413
1414 Output Parameter:
1415 . function - constructor routine
1416
1417 Notes:
1418 Calling sequence for the routine:
1419 $ my_csp(PC pc, PetscInt l, DM dm, KSP smooth, PetscInt Nc, const Vec initGuess[], Vec **coarseSp)
1420 $ pc - The PC object
1421 $ l - The multigrid level, 0 is the coarse level
1422 $ dm - The DM for this level
1423 $ smooth - The level smoother
1424 $ Nc - The size of the coarse space
1425 $ initGuess - Basis for an initial guess for the space
1426 $ coarseSp - A basis for the computed coarse space
1427
1428 Level: advanced
1429
1430 .seealso: PCMGRegisterCoarseSpaceConstructor(), PCRegister()
1431 @*/
PCMGGetCoarseSpaceConstructor(const char name[],PetscErrorCode (** function)(PC,PetscInt,DM,KSP,PetscInt,const Vec[],Vec **))1432 PetscErrorCode PCMGGetCoarseSpaceConstructor(const char name[], PetscErrorCode (**function)(PC, PetscInt, DM, KSP, PetscInt, const Vec[], Vec **))
1433 {
1434 PetscErrorCode ierr;
1435
1436 PetscFunctionBegin;
1437 ierr = PetscFunctionListFind(PCMGCoarseList,name,function);CHKERRQ(ierr);
1438 PetscFunctionReturn(0);
1439 }
1440
1441 /* ----------------------------------------------------------------------------------------*/
1442
1443 /*MC
1444 PCMG - Use multigrid preconditioning. This preconditioner requires you provide additional
1445 information about the coarser grid matrices and restriction/interpolation operators.
1446
1447 Options Database Keys:
1448 + -pc_mg_levels <nlevels> - number of levels including finest
1449 . -pc_mg_cycle_type <v,w> - provide the cycle desired
1450 . -pc_mg_type <additive,multiplicative,full,kaskade> - multiplicative is the default
1451 . -pc_mg_log - log information about time spent on each level of the solver
1452 . -pc_mg_distinct_smoothup - configure up (after interpolation) and down (before restriction) smoothers separately (with different options prefixes)
1453 . -pc_mg_galerkin <both,pmat,mat,none> - use Galerkin process to compute coarser operators, i.e. Acoarse = R A R'
1454 . -pc_mg_multiplicative_cycles - number of cycles to use as the preconditioner (defaults to 1)
1455 . -pc_mg_dump_matlab - dumps the matrices for each level and the restriction/interpolation matrices
1456 to the Socket viewer for reading from MATLAB.
1457 - -pc_mg_dump_binary - dumps the matrices for each level and the restriction/interpolation matrices
1458 to the binary output file called binaryoutput
1459
1460 Notes:
1461 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
1462
1463 When run with a single level the smoother options are used on that level NOT the coarse grid solver options
1464
1465 When run with KSPRICHARDSON the convergence test changes slightly if monitor is turned on. The iteration count may change slightly. This
1466 is because without monitoring the residual norm is computed WITHIN each multigrid cycle on the finest level after the pre-smoothing
1467 (because the residual has just been computed for the multigrid algorithm and is hence available for free) while with monitoring the
1468 residual is computed at the end of each cycle.
1469
1470 Level: intermediate
1471
1472 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType, PCEXOTIC, PCGAMG, PCML, PCHYPRE
1473 PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), PCMGSetCycleType(),
1474 PCMGSetDistinctSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(),
1475 PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(),
1476 PCMGSetCycleTypeOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR()
1477 M*/
1478
PCCreate_MG(PC pc)1479 PETSC_EXTERN PetscErrorCode PCCreate_MG(PC pc)
1480 {
1481 PC_MG *mg;
1482 PetscErrorCode ierr;
1483
1484 PetscFunctionBegin;
1485 ierr = PetscNewLog(pc,&mg);CHKERRQ(ierr);
1486 pc->data = (void*)mg;
1487 mg->nlevels = -1;
1488 mg->am = PC_MG_MULTIPLICATIVE;
1489 mg->galerkin = PC_MG_GALERKIN_NONE;
1490 mg->adaptInterpolation = PETSC_FALSE;
1491 mg->Nc = -1;
1492 mg->eigenvalue = -1;
1493
1494 pc->useAmat = PETSC_TRUE;
1495
1496 pc->ops->apply = PCApply_MG;
1497 pc->ops->setup = PCSetUp_MG;
1498 pc->ops->reset = PCReset_MG;
1499 pc->ops->destroy = PCDestroy_MG;
1500 pc->ops->setfromoptions = PCSetFromOptions_MG;
1501 pc->ops->view = PCView_MG;
1502
1503 ierr = PetscObjectComposedDataRegister(&mg->eigenvalue);CHKERRQ(ierr);
1504 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGSetGalerkin_C",PCMGSetGalerkin_MG);CHKERRQ(ierr);
1505 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGGetLevels_C",PCMGGetLevels_MG);CHKERRQ(ierr);
1506 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGSetLevels_C",PCMGSetLevels_MG);CHKERRQ(ierr);
1507 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGetInterpolations_C",PCGetInterpolations_MG);CHKERRQ(ierr);
1508 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGetCoarseOperators_C",PCGetCoarseOperators_MG);CHKERRQ(ierr);
1509 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGSetAdaptInterpolation_C",PCMGSetAdaptInterpolation_MG);CHKERRQ(ierr);
1510 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGGetAdaptInterpolation_C",PCMGGetAdaptInterpolation_MG);CHKERRQ(ierr);
1511 PetscFunctionReturn(0);
1512 }
1513