Actual source code: hmg.c
1: #include <petscdm.h>
2: #include <petscctable.h>
3: #include <petsc/private/matimpl.h>
4: #include <petsc/private/pcmgimpl.h>
5: #include <petsc/private/pcimpl.h>
7: typedef struct {
8: PC innerpc; /* A MG inner PC (Hypre or PCGAMG) to setup interpolations and coarse operators */
9: char* innerpctype; /* PCGAMG or PCHYPRE */
10: PetscBool reuseinterp; /* A flag indicates if or not to reuse the interpolations */
11: PetscBool subcoarsening; /* If or not to use a subspace-based coarsening algorithm */
12: PetscBool usematmaij; /* If or not to use MatMAIJ for saving memory */
13: PetscInt component; /* Which subspace is used for the subspace-based coarsening algorithm? */
14: } PC_HMG;
16: PetscErrorCode PCSetFromOptions_HMG(PetscOptionItems*,PC);
17: PetscErrorCode PCReset_MG(PC);
19: static PetscErrorCode PCHMGExtractSubMatrix_Private(Mat pmat,Mat *submat,MatReuse reuse,PetscInt component,PetscInt blocksize)
20: {
21: IS isrow;
22: PetscInt rstart,rend;
23: MPI_Comm comm;
25: PetscObjectGetComm((PetscObject)pmat,&comm);
27: MatGetOwnershipRange(pmat,&rstart,&rend);
29: ISCreateStride(comm,(rend-rstart)/blocksize,rstart+component,blocksize,&isrow);
30: MatCreateSubMatrix(pmat,isrow,isrow,reuse,submat);
31: ISDestroy(&isrow);
32: return 0;
33: }
35: static PetscErrorCode PCHMGExpandInterpolation_Private(Mat subinterp, Mat *interp, PetscInt blocksize)
36: {
37: PetscInt subrstart,subrend,subrowsize,subcolsize,subcstart,subcend,rowsize,colsize;
38: PetscInt subrow,row,nz,*d_nnz,*o_nnz,i,j,dnz,onz,max_nz,*indices;
39: const PetscInt *idx;
40: const PetscScalar *values;
41: MPI_Comm comm;
43: PetscObjectGetComm((PetscObject)subinterp,&comm);
44: MatGetOwnershipRange(subinterp,&subrstart,&subrend);
45: subrowsize = subrend-subrstart;
46: rowsize = subrowsize*blocksize;
47: PetscCalloc2(rowsize,&d_nnz,rowsize,&o_nnz);
48: MatGetOwnershipRangeColumn(subinterp,&subcstart,&subcend);
49: subcolsize = subcend - subcstart;
50: colsize = subcolsize*blocksize;
51: max_nz = 0;
52: for (subrow=subrstart;subrow<subrend;subrow++) {
53: MatGetRow(subinterp,subrow,&nz,&idx,NULL);
54: if (max_nz<nz) max_nz = nz;
55: dnz = 0; onz = 0;
56: for (i=0;i<nz;i++) {
57: if (idx[i]>=subcstart && idx[i]<subcend) dnz++;
58: else onz++;
59: }
60: for (i=0;i<blocksize;i++) {
61: d_nnz[(subrow-subrstart)*blocksize+i] = dnz;
62: o_nnz[(subrow-subrstart)*blocksize+i] = onz;
63: }
64: MatRestoreRow(subinterp,subrow,&nz,&idx,NULL);
65: }
66: MatCreateAIJ(comm,rowsize,colsize,PETSC_DETERMINE,PETSC_DETERMINE,0,d_nnz,0,o_nnz,interp);
67: MatSetOption(*interp,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
68: MatSetOption(*interp,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
69: MatSetOption(*interp,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
70: MatSetFromOptions(*interp);
72: MatSetUp(*interp);
73: PetscFree2(d_nnz,o_nnz);
74: PetscMalloc1(max_nz,&indices);
75: for (subrow=subrstart; subrow<subrend; subrow++) {
76: MatGetRow(subinterp,subrow,&nz,&idx,&values);
77: for (i=0;i<blocksize;i++) {
78: row = subrow*blocksize+i;
79: for (j=0;j<nz;j++) {
80: indices[j] = idx[j]*blocksize+i;
81: }
82: MatSetValues(*interp,1,&row,nz,indices,values,INSERT_VALUES);
83: }
84: MatRestoreRow(subinterp,subrow,&nz,&idx,&values);
85: }
86: PetscFree(indices);
87: MatAssemblyBegin(*interp,MAT_FINAL_ASSEMBLY);
88: MatAssemblyEnd(*interp,MAT_FINAL_ASSEMBLY);
89: return 0;
90: }
92: PetscErrorCode PCSetUp_HMG(PC pc)
93: {
94: PetscErrorCode ierr;
95: Mat PA, submat;
96: PC_MG *mg = (PC_MG*)pc->data;
97: PC_HMG *hmg = (PC_HMG*) mg->innerctx;
98: MPI_Comm comm;
99: PetscInt level;
100: PetscInt num_levels;
101: Mat *operators,*interpolations;
102: PetscInt blocksize;
103: const char *prefix;
104: PCMGGalerkinType galerkin;
106: PetscObjectGetComm((PetscObject)pc,&comm);
107: if (pc->setupcalled) {
108: if (hmg->reuseinterp) {
109: /* If we did not use Galerkin in the last call or we have a different sparsity pattern now,
110: * we have to build from scratch
111: * */
112: PCMGGetGalerkin(pc,&galerkin);
113: if (galerkin == PC_MG_GALERKIN_NONE || pc->flag != SAME_NONZERO_PATTERN) pc->setupcalled = PETSC_FALSE;
114: PCMGSetGalerkin(pc,PC_MG_GALERKIN_PMAT);
115: PCSetUp_MG(pc);
116: return 0;
117: } else {
118: PCReset_MG(pc);
119: pc->setupcalled = PETSC_FALSE;
120: }
121: }
123: /* Create an inner PC (GAMG or HYPRE) */
124: if (!hmg->innerpc) {
125: PCCreate(comm,&hmg->innerpc);
126: /* If users do not set an inner pc type, we need to set a default value */
127: if (!hmg->innerpctype) {
128: /* If hypre is available, use hypre, otherwise, use gamg */
129: #if PETSC_HAVE_HYPRE
130: PetscStrallocpy(PCHYPRE,&(hmg->innerpctype));
131: #else
132: PetscStrallocpy(PCGAMG,&(hmg->innerpctype));
133: #endif
134: }
135: PCSetType(hmg->innerpc,hmg->innerpctype);
136: }
137: PCGetOperators(pc,NULL,&PA);
138: /* Users need to correctly set a block size of matrix in order to use subspace coarsening */
139: MatGetBlockSize(PA,&blocksize);
140: if (blocksize<=1) hmg->subcoarsening = PETSC_FALSE;
141: /* Extract a submatrix for constructing subinterpolations */
142: if (hmg->subcoarsening) {
143: PCHMGExtractSubMatrix_Private(PA,&submat,MAT_INITIAL_MATRIX,hmg->component,blocksize);
144: PA = submat;
145: }
146: PCSetOperators(hmg->innerpc,PA,PA);
147: if (hmg->subcoarsening) {
148: MatDestroy(&PA);
149: }
150: /* Setup inner PC correctly. During this step, matrix will be coarsened */
151: PCSetUseAmat(hmg->innerpc,PETSC_FALSE);
152: PetscObjectGetOptionsPrefix((PetscObject)pc,&prefix);
153: PetscObjectSetOptionsPrefix((PetscObject)hmg->innerpc,prefix);
154: PetscObjectAppendOptionsPrefix((PetscObject)hmg->innerpc,"hmg_inner_");
155: PCSetFromOptions(hmg->innerpc);
156: PCSetUp(hmg->innerpc);
158: /* Obtain interpolations IN PLACE. For BoomerAMG, (I,J,data) is reused to avoid memory overhead */
159: PCGetInterpolations(hmg->innerpc,&num_levels,&interpolations);
160: /* We can reuse the coarse operators when we do the full space coarsening */
161: if (!hmg->subcoarsening) {
162: PCGetCoarseOperators(hmg->innerpc,&num_levels,&operators);
163: }
165: PCDestroy(&hmg->innerpc);
166: hmg->innerpc = NULL;
167: PCMGSetLevels_MG(pc,num_levels,NULL);
168: /* Set coarse matrices and interpolations to PCMG */
169: for (level=num_levels-1; level>0; level--) {
170: Mat P=NULL, pmat=NULL;
171: Vec b, x,r;
172: if (hmg->subcoarsening) {
173: if (hmg->usematmaij) {
174: MatCreateMAIJ(interpolations[level-1],blocksize,&P);
175: MatDestroy(&interpolations[level-1]);
176: } else {
177: /* Grow interpolation. In the future, we should use MAIJ */
178: PCHMGExpandInterpolation_Private(interpolations[level-1],&P,blocksize);
179: MatDestroy(&interpolations[level-1]);
180: }
181: } else {
182: P = interpolations[level-1];
183: }
184: MatCreateVecs(P,&b,&r);
185: PCMGSetInterpolation(pc,level,P);
186: PCMGSetRestriction(pc,level,P);
187: MatDestroy(&P);
188: /* We reuse the matrices when we do not do subspace coarsening */
189: if ((level-1)>=0 && !hmg->subcoarsening) {
190: pmat = operators[level-1];
191: PCMGSetOperators(pc,level-1,pmat,pmat);
192: MatDestroy(&pmat);
193: }
194: PCMGSetRhs(pc,level-1,b);
196: PCMGSetR(pc,level,r);
197: VecDestroy(&r);
199: VecDuplicate(b,&x);
200: PCMGSetX(pc,level-1,x);
201: VecDestroy(&x);
202: VecDestroy(&b);
203: }
204: PetscFree(interpolations);
205: if (!hmg->subcoarsening) {
206: PetscFree(operators);
207: }
208: /* Turn Galerkin off when we already have coarse operators */
209: PCMGSetGalerkin(pc,hmg->subcoarsening ? PC_MG_GALERKIN_PMAT:PC_MG_GALERKIN_NONE);
210: PCSetDM(pc,NULL);
211: PCSetUseAmat(pc,PETSC_FALSE);
212: PetscObjectOptionsBegin((PetscObject)pc);
213: PCSetFromOptions_MG(PetscOptionsObject,pc); /* should be called in PCSetFromOptions_HMG(), but cannot be called prior to PCMGSetLevels() */
214: PetscOptionsEnd();
215: PCSetUp_MG(pc);
216: return 0;
217: }
219: PetscErrorCode PCDestroy_HMG(PC pc)
220: {
221: PC_MG *mg = (PC_MG*)pc->data;
222: PC_HMG *hmg = (PC_HMG*) mg->innerctx;
224: PCDestroy(&hmg->innerpc);
225: PetscFree(hmg->innerpctype);
226: PetscFree(hmg);
227: PCDestroy_MG(pc);
229: PetscObjectComposeFunction((PetscObject)pc,"PCHMGSetReuseInterpolation_C",NULL);
230: PetscObjectComposeFunction((PetscObject)pc,"PCHMGSetUseSubspaceCoarsening_C",NULL);
231: PetscObjectComposeFunction((PetscObject)pc,"PCHMGSetInnerPCType_C",NULL);
232: PetscObjectComposeFunction((PetscObject)pc,"PCHMGSetCoarseningComponent_C",NULL);
233: return 0;
234: }
236: PetscErrorCode PCView_HMG(PC pc,PetscViewer viewer)
237: {
238: PC_MG *mg = (PC_MG*)pc->data;
239: PC_HMG *hmg = (PC_HMG*) mg->innerctx;
240: PetscBool iascii;
242: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
243: if (iascii) {
244: PetscViewerASCIIPrintf(viewer," Reuse interpolation: %s\n",hmg->reuseinterp? "true":"false");
245: PetscViewerASCIIPrintf(viewer," Use subspace coarsening: %s\n",hmg->subcoarsening? "true":"false");
246: PetscViewerASCIIPrintf(viewer," Coarsening component: %D \n",hmg->component);
247: PetscViewerASCIIPrintf(viewer," Use MatMAIJ: %s \n",hmg->usematmaij? "true":"false");
248: PetscViewerASCIIPrintf(viewer," Inner PC type: %s \n",hmg->innerpctype);
249: }
250: PCView_MG(pc,viewer);
251: return 0;
252: }
254: PetscErrorCode PCSetFromOptions_HMG(PetscOptionItems *PetscOptionsObject,PC pc)
255: {
256: PC_MG *mg = (PC_MG*)pc->data;
257: PC_HMG *hmg = (PC_HMG*) mg->innerctx;
259: PetscOptionsHead(PetscOptionsObject,"HMG");
260: PetscOptionsBool("-pc_hmg_reuse_interpolation","Reuse the interpolation operators when possible (cheaper, weaker when matrix entries change a lot)","PCHMGSetReuseInterpolation",hmg->reuseinterp,&hmg->reuseinterp,NULL);
261: PetscOptionsBool("-pc_hmg_use_subspace_coarsening","Use the subspace coarsening to compute the interpolations","PCHMGSetUseSubspaceCoarsening",hmg->subcoarsening,&hmg->subcoarsening,NULL);
262: PetscOptionsBool("-pc_hmg_use_matmaij","Use MatMAIJ store interpolation for saving memory","PCHMGSetInnerPCType",hmg->usematmaij,&hmg->usematmaij,NULL);
263: PetscOptionsInt("-pc_hmg_coarsening_component","Which component is chosen for the subspace-based coarsening algorithm","PCHMGSetCoarseningComponent",hmg->component,&hmg->component,NULL);
264: PetscOptionsTail();
265: return 0;
266: }
268: static PetscErrorCode PCHMGSetReuseInterpolation_HMG(PC pc, PetscBool reuse)
269: {
270: PC_MG *mg = (PC_MG*)pc->data;
271: PC_HMG *hmg = (PC_HMG*) mg->innerctx;
273: hmg->reuseinterp = reuse;
274: return 0;
275: }
277: /*@
278: PCHMGSetReuseInterpolation - Reuse interpolation matrices in HMG
280: Logically Collective on PC
282: Input Parameters:
283: + pc - the HMG context
284: - reuse - True indicates that HMG will reuse the interpolations
286: Options Database Keys:
287: . -pc_hmg_reuse_interpolation <true | false> - Whether or not to reuse the interpolations. If true, it potentially save the compute time.
289: Level: beginner
291: .keywords: HMG, multigrid, interpolation, reuse, set
293: .seealso: PCHMG
294: @*/
295: PetscErrorCode PCHMGSetReuseInterpolation(PC pc, PetscBool reuse)
296: {
298: PetscUseMethod(pc,"PCHMGSetReuseInterpolation_C",(PC,PetscBool),(pc,reuse));
299: return 0;
300: }
302: static PetscErrorCode PCHMGSetUseSubspaceCoarsening_HMG(PC pc, PetscBool subspace)
303: {
304: PC_MG *mg = (PC_MG*)pc->data;
305: PC_HMG *hmg = (PC_HMG*) mg->innerctx;
307: hmg->subcoarsening = subspace;
308: return 0;
309: }
311: /*@
312: PCHMGSetUseSubspaceCoarsening - Use subspace coarsening in HMG
314: Logically Collective on PC
316: Input Parameters:
317: + pc - the HMG context
318: - reuse - True indicates that HMG will use the subspace coarsening
320: Options Database Keys:
321: . -pc_hmg_use_subspace_coarsening <true | false> - Whether or not to use subspace coarsening (that is, coarsen a submatrix).
323: Level: beginner
325: .keywords: HMG, multigrid, interpolation, subspace, coarsening
327: .seealso: PCHMG
328: @*/
329: PetscErrorCode PCHMGSetUseSubspaceCoarsening(PC pc, PetscBool subspace)
330: {
332: PetscUseMethod(pc,"PCHMGSetUseSubspaceCoarsening_C",(PC,PetscBool),(pc,subspace));
333: return 0;
334: }
336: static PetscErrorCode PCHMGSetInnerPCType_HMG(PC pc, PCType type)
337: {
338: PC_MG *mg = (PC_MG*)pc->data;
339: PC_HMG *hmg = (PC_HMG*) mg->innerctx;
341: PetscStrallocpy(type,&(hmg->innerpctype));
342: return 0;
343: }
345: /*@C
346: PCHMGSetInnerPCType - Set an inner PC type
348: Logically Collective on PC
350: Input Parameters:
351: + pc - the HMG context
352: - type - <hypre, gamg> coarsening algorithm
354: Options Database Keys:
355: . -hmg_inner_pc_type <hypre, gamg> - What method is used to coarsen matrix
357: Level: beginner
359: .keywords: HMG, multigrid, interpolation, coarsening
361: .seealso: PCHMG, PCType
362: @*/
363: PetscErrorCode PCHMGSetInnerPCType(PC pc, PCType type)
364: {
366: PetscUseMethod(pc,"PCHMGSetInnerPCType_C",(PC,PCType),(pc,type));
367: return 0;
368: }
370: static PetscErrorCode PCHMGSetCoarseningComponent_HMG(PC pc, PetscInt component)
371: {
372: PC_MG *mg = (PC_MG*)pc->data;
373: PC_HMG *hmg = (PC_HMG*) mg->innerctx;
375: hmg->component = component;
376: return 0;
377: }
379: /*@
380: PCHMGSetCoarseningComponent - Set which component is used for the subspace-based coarsening algorithm
382: Logically Collective on PC
384: Input Parameters:
385: + pc - the HMG context
386: - component - which component PC will coarsen
388: Options Database Keys:
389: . -pc_hmg_coarsening_component - Which component is chosen for the subspace-based coarsening algorithm
391: Level: beginner
393: .keywords: HMG, multigrid, interpolation, coarsening, component
395: .seealso: PCHMG, PCType
396: @*/
397: PetscErrorCode PCHMGSetCoarseningComponent(PC pc, PetscInt component)
398: {
400: PetscUseMethod(pc,"PCHMGSetCoarseningComponent_C",(PC,PetscInt),(pc,component));
401: return 0;
402: }
404: static PetscErrorCode PCHMGUseMatMAIJ_HMG(PC pc, PetscBool usematmaij)
405: {
406: PC_MG *mg = (PC_MG*)pc->data;
407: PC_HMG *hmg = (PC_HMG*) mg->innerctx;
409: hmg->usematmaij = usematmaij;
410: return 0;
411: }
413: /*@
414: PCHMGUseMatMAIJ - Set a flag that indicates if or not to use MatMAIJ for interpolations for saving memory
416: Logically Collective on PC
418: Input Parameters:
419: + pc - the HMG context
420: - usematmaij - if or not to use MatMAIJ for interpolations. By default, it is true for saving memory
422: Options Database Keys:
423: . -pc_hmg_use_matmaij - <true | false >
425: Level: beginner
427: .keywords: HMG, multigrid, interpolation, coarsening, MatMAIJ
429: .seealso: PCHMG, PCType
430: @*/
431: PetscErrorCode PCHMGUseMatMAIJ(PC pc, PetscBool usematmaij)
432: {
434: PetscUseMethod(pc,"PCHMGUseMatMAIJ_C",(PC,PetscBool),(pc,usematmaij));
435: return 0;
436: }
438: /*MC
439: PCHMG - Hybrid of PETSc preconditioners (such as ASM, BJacobi, SOR, etc.) and Hypre BoomerAMG, GAMG or other multilevel methods. BoomerAMG, GAMG
440: or other multilevel methods is used to coarsen matrix and generate a sequence of coarse matrices and interpolations. The matrices and
441: interpolations are employed to construct PCMG, and then any available PETSc preconditioners can be chosen as smoothers and the coarse solver.
443: Options Database Keys:
444: + -pc_hmg_reuse_interpolation <true | false> - Whether or not to reuse the interpolations. If true, it potentially save the compute time.
445: . -pc_hmg_use_subspace_coarsening <true | false> - Whether or not to use subspace coarsening (that is, coarsen a submatrix).
446: . -hmg_inner_pc_type <hypre, gamg, ...> - What method is used to coarsen matrix
447: - -pc_hmg_use_matmaij <true | false> - Whether or not to use MatMAIJ for multicomponent problems for saving memory
449: Notes:
450: For multicomponent problems, we can just coarsen one submatrix associated with one particular component. In this way, the preconditioner setup
451: time is significantly reduced. One typical use case is neutron transport equations. There are many variables on each mesh vertex due to the
452: of angle and energy. Each variable, in fact, corresponds to the same PDEs but with different material properties.
454: Level: beginner
456: References:
457: . * - Fande Kong, Yaqi Wang, Derek R Gaston, Cody J Permann, Andrew E Slaughter, Alexander D Lindsay, Richard C Martineau, A highly parallel multilevel
458: Newton-Krylov-Schwarz method with subspace-based coarsening and partition-based balancing for the multigroup neutron transport equations on
459: 3D unstructured meshes, arXiv preprint arXiv:1903.03659, 2019
461: .seealso: PCCreate(), PCSetType(), PCType, PC, PCMG, PCHYPRE, PCHMG, PCGetCoarseOperators(), PCGetInterpolations(), PCHMGSetReuseInterpolation(), PCHMGSetUseSubspaceCoarsening(),
462: PCHMGSetInnerPCType()
464: M*/
465: PETSC_EXTERN PetscErrorCode PCCreate_HMG(PC pc)
466: {
467: PC_HMG *hmg;
468: PC_MG *mg;
470: /* if type was previously mg; must manually destroy it because call to PCSetType(pc,PCMG) will not destroy it */
471: if (pc->ops->destroy) {
472: (*pc->ops->destroy)(pc);
473: pc->data = NULL;
474: }
475: PetscFree(((PetscObject)pc)->type_name);
477: PCSetType(pc,PCMG);
478: PetscObjectChangeTypeName((PetscObject)pc, PCHMG);
479: PetscNew(&hmg);
481: mg = (PC_MG*) pc->data;
482: mg->innerctx = hmg;
483: hmg->reuseinterp = PETSC_FALSE;
484: hmg->subcoarsening = PETSC_FALSE;
485: hmg->usematmaij = PETSC_TRUE;
486: hmg->component = 0;
487: hmg->innerpc = NULL;
489: pc->ops->setfromoptions = PCSetFromOptions_HMG;
490: pc->ops->view = PCView_HMG;
491: pc->ops->destroy = PCDestroy_HMG;
492: pc->ops->setup = PCSetUp_HMG;
494: PetscObjectComposeFunction((PetscObject)pc,"PCHMGSetReuseInterpolation_C",PCHMGSetReuseInterpolation_HMG);
495: PetscObjectComposeFunction((PetscObject)pc,"PCHMGSetUseSubspaceCoarsening_C",PCHMGSetUseSubspaceCoarsening_HMG);
496: PetscObjectComposeFunction((PetscObject)pc,"PCHMGSetInnerPCType_C",PCHMGSetInnerPCType_HMG);
497: PetscObjectComposeFunction((PetscObject)pc,"PCHMGSetCoarseningComponent_C",PCHMGSetCoarseningComponent_HMG);
498: PetscObjectComposeFunction((PetscObject)pc,"PCHMGUseMatMAIJ_C",PCHMGUseMatMAIJ_HMG);
499: return 0;
500: }