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(PC, PetscOptionItems *);
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;
56: onz = 0;
57: for (i = 0; i < nz; i++) {
58: if (idx[i] >= subcstart && idx[i] < subcend) dnz++;
59: else onz++;
60: }
61: for (i = 0; i < blocksize; i++) {
62: d_nnz[(subrow - subrstart) * blocksize + i] = dnz;
63: o_nnz[(subrow - subrstart) * blocksize + i] = onz;
64: }
65: MatRestoreRow(subinterp, subrow, &nz, &idx, NULL);
66: }
67: MatCreateAIJ(comm, rowsize, colsize, PETSC_DETERMINE, PETSC_DETERMINE, 0, d_nnz, 0, o_nnz, interp);
68: MatSetOption(*interp, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE);
69: MatSetOption(*interp, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);
70: MatSetOption(*interp, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
71: MatSetFromOptions(*interp);
73: MatSetUp(*interp);
74: PetscFree2(d_nnz, o_nnz);
75: PetscMalloc1(max_nz, &indices);
76: for (subrow = subrstart; subrow < subrend; subrow++) {
77: MatGetRow(subinterp, subrow, &nz, &idx, &values);
78: for (i = 0; i < blocksize; i++) {
79: row = subrow * blocksize + i;
80: for (j = 0; j < nz; j++) indices[j] = idx[j] * blocksize + i;
81: MatSetValues(*interp, 1, &row, nz, indices, values, INSERT_VALUES);
82: }
83: MatRestoreRow(subinterp, subrow, &nz, &idx, &values);
84: }
85: PetscFree(indices);
86: MatAssemblyBegin(*interp, MAT_FINAL_ASSEMBLY);
87: MatAssemblyEnd(*interp, MAT_FINAL_ASSEMBLY);
88: return 0;
89: }
91: PetscErrorCode PCSetUp_HMG(PC pc)
92: {
93: Mat PA, submat;
94: PC_MG *mg = (PC_MG *)pc->data;
95: PC_HMG *hmg = (PC_HMG *)mg->innerctx;
96: MPI_Comm comm;
97: PetscInt level;
98: PetscInt num_levels;
99: Mat *operators, *interpolations;
100: PetscInt blocksize;
101: const char *prefix;
102: PCMGGalerkinType galerkin;
104: PetscObjectGetComm((PetscObject)pc, &comm);
105: if (pc->setupcalled) {
106: if (hmg->reuseinterp) {
107: /* If we did not use Galerkin in the last call or we have a different sparsity pattern now,
108: * we have to build from scratch
109: * */
110: PCMGGetGalerkin(pc, &galerkin);
111: if (galerkin == PC_MG_GALERKIN_NONE || pc->flag != SAME_NONZERO_PATTERN) pc->setupcalled = PETSC_FALSE;
112: PCMGSetGalerkin(pc, PC_MG_GALERKIN_PMAT);
113: PCSetUp_MG(pc);
114: return 0;
115: } else {
116: PCReset_MG(pc);
117: pc->setupcalled = PETSC_FALSE;
118: }
119: }
121: /* Create an inner PC (GAMG or HYPRE) */
122: if (!hmg->innerpc) {
123: PCCreate(comm, &hmg->innerpc);
124: /* If users do not set an inner pc type, we need to set a default value */
125: if (!hmg->innerpctype) {
126: /* If hypre is available, use hypre, otherwise, use gamg */
127: #if PETSC_HAVE_HYPRE
128: PetscStrallocpy(PCHYPRE, &(hmg->innerpctype));
129: #else
130: PetscStrallocpy(PCGAMG, &(hmg->innerpctype));
131: #endif
132: }
133: PCSetType(hmg->innerpc, hmg->innerpctype);
134: }
135: PCGetOperators(pc, NULL, &PA);
136: /* Users need to correctly set a block size of matrix in order to use subspace coarsening */
137: MatGetBlockSize(PA, &blocksize);
138: if (blocksize <= 1) hmg->subcoarsening = PETSC_FALSE;
139: /* Extract a submatrix for constructing subinterpolations */
140: if (hmg->subcoarsening) {
141: PCHMGExtractSubMatrix_Private(PA, &submat, MAT_INITIAL_MATRIX, hmg->component, blocksize);
142: PA = submat;
143: }
144: PCSetOperators(hmg->innerpc, PA, PA);
145: if (hmg->subcoarsening) MatDestroy(&PA);
146: /* Setup inner PC correctly. During this step, matrix will be coarsened */
147: PCSetUseAmat(hmg->innerpc, PETSC_FALSE);
148: PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix);
149: PetscObjectSetOptionsPrefix((PetscObject)hmg->innerpc, prefix);
150: PetscObjectAppendOptionsPrefix((PetscObject)hmg->innerpc, "hmg_inner_");
151: PCSetFromOptions(hmg->innerpc);
152: PCSetUp(hmg->innerpc);
154: /* Obtain interpolations IN PLACE. For BoomerAMG, (I,J,data) is reused to avoid memory overhead */
155: PCGetInterpolations(hmg->innerpc, &num_levels, &interpolations);
156: /* We can reuse the coarse operators when we do the full space coarsening */
157: if (!hmg->subcoarsening) PCGetCoarseOperators(hmg->innerpc, &num_levels, &operators);
159: PCDestroy(&hmg->innerpc);
160: hmg->innerpc = NULL;
161: PCMGSetLevels_MG(pc, num_levels, NULL);
162: /* Set coarse matrices and interpolations to PCMG */
163: for (level = num_levels - 1; level > 0; level--) {
164: Mat P = NULL, pmat = NULL;
165: Vec b, x, r;
166: if (hmg->subcoarsening) {
167: if (hmg->usematmaij) {
168: MatCreateMAIJ(interpolations[level - 1], blocksize, &P);
169: MatDestroy(&interpolations[level - 1]);
170: } else {
171: /* Grow interpolation. In the future, we should use MAIJ */
172: PCHMGExpandInterpolation_Private(interpolations[level - 1], &P, blocksize);
173: MatDestroy(&interpolations[level - 1]);
174: }
175: } else {
176: P = interpolations[level - 1];
177: }
178: MatCreateVecs(P, &b, &r);
179: PCMGSetInterpolation(pc, level, P);
180: PCMGSetRestriction(pc, level, P);
181: MatDestroy(&P);
182: /* We reuse the matrices when we do not do subspace coarsening */
183: if ((level - 1) >= 0 && !hmg->subcoarsening) {
184: pmat = operators[level - 1];
185: PCMGSetOperators(pc, level - 1, pmat, pmat);
186: MatDestroy(&pmat);
187: }
188: PCMGSetRhs(pc, level - 1, b);
190: PCMGSetR(pc, level, r);
191: VecDestroy(&r);
193: VecDuplicate(b, &x);
194: PCMGSetX(pc, level - 1, x);
195: VecDestroy(&x);
196: VecDestroy(&b);
197: }
198: PetscFree(interpolations);
199: if (!hmg->subcoarsening) PetscFree(operators);
200: /* Turn Galerkin off when we already have coarse operators */
201: PCMGSetGalerkin(pc, hmg->subcoarsening ? PC_MG_GALERKIN_PMAT : PC_MG_GALERKIN_NONE);
202: PCSetDM(pc, NULL);
203: PCSetUseAmat(pc, PETSC_FALSE);
204: PetscObjectOptionsBegin((PetscObject)pc);
205: PCSetFromOptions_MG(pc, PetscOptionsObject); /* should be called in PCSetFromOptions_HMG(), but cannot be called prior to PCMGSetLevels() */
206: PetscOptionsEnd();
207: PCSetUp_MG(pc);
208: return 0;
209: }
211: PetscErrorCode PCDestroy_HMG(PC pc)
212: {
213: PC_MG *mg = (PC_MG *)pc->data;
214: PC_HMG *hmg = (PC_HMG *)mg->innerctx;
216: PCDestroy(&hmg->innerpc);
217: PetscFree(hmg->innerpctype);
218: PetscFree(hmg);
219: PCDestroy_MG(pc);
221: PetscObjectComposeFunction((PetscObject)pc, "PCHMGSetReuseInterpolation_C", NULL);
222: PetscObjectComposeFunction((PetscObject)pc, "PCHMGSetUseSubspaceCoarsening_C", NULL);
223: PetscObjectComposeFunction((PetscObject)pc, "PCHMGSetInnerPCType_C", NULL);
224: PetscObjectComposeFunction((PetscObject)pc, "PCHMGSetCoarseningComponent_C", NULL);
225: PetscObjectComposeFunction((PetscObject)pc, "PCHMGUseMatMAIJ_C", NULL);
226: return 0;
227: }
229: PetscErrorCode PCView_HMG(PC pc, PetscViewer viewer)
230: {
231: PC_MG *mg = (PC_MG *)pc->data;
232: PC_HMG *hmg = (PC_HMG *)mg->innerctx;
233: PetscBool iascii;
235: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
236: if (iascii) {
237: PetscViewerASCIIPrintf(viewer, " Reuse interpolation: %s\n", hmg->reuseinterp ? "true" : "false");
238: PetscViewerASCIIPrintf(viewer, " Use subspace coarsening: %s\n", hmg->subcoarsening ? "true" : "false");
239: PetscViewerASCIIPrintf(viewer, " Coarsening component: %" PetscInt_FMT " \n", hmg->component);
240: PetscViewerASCIIPrintf(viewer, " Use MatMAIJ: %s \n", hmg->usematmaij ? "true" : "false");
241: PetscViewerASCIIPrintf(viewer, " Inner PC type: %s \n", hmg->innerpctype);
242: }
243: PCView_MG(pc, viewer);
244: return 0;
245: }
247: PetscErrorCode PCSetFromOptions_HMG(PC pc, PetscOptionItems *PetscOptionsObject)
248: {
249: PC_MG *mg = (PC_MG *)pc->data;
250: PC_HMG *hmg = (PC_HMG *)mg->innerctx;
252: PetscOptionsHeadBegin(PetscOptionsObject, "HMG");
253: 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);
254: PetscOptionsBool("-pc_hmg_use_subspace_coarsening", "Use the subspace coarsening to compute the interpolations", "PCHMGSetUseSubspaceCoarsening", hmg->subcoarsening, &hmg->subcoarsening, NULL);
255: PetscOptionsBool("-pc_hmg_use_matmaij", "Use MatMAIJ store interpolation for saving memory", "PCHMGSetInnerPCType", hmg->usematmaij, &hmg->usematmaij, NULL);
256: PetscOptionsInt("-pc_hmg_coarsening_component", "Which component is chosen for the subspace-based coarsening algorithm", "PCHMGSetCoarseningComponent", hmg->component, &hmg->component, NULL);
257: PetscOptionsHeadEnd();
258: return 0;
259: }
261: static PetscErrorCode PCHMGSetReuseInterpolation_HMG(PC pc, PetscBool reuse)
262: {
263: PC_MG *mg = (PC_MG *)pc->data;
264: PC_HMG *hmg = (PC_HMG *)mg->innerctx;
266: hmg->reuseinterp = reuse;
267: return 0;
268: }
270: /*@
271: PCHMGSetReuseInterpolation - Reuse the interpolation matrices in `PCHMG` after changing the matrices numerical values
273: Logically Collective
275: Input Parameters:
276: + pc - the `PCHMG` context
277: - reuse - `PETSC_TRUE` indicates that `PCHMG` will reuse the interpolations
279: Options Database Key:
280: . -pc_hmg_reuse_interpolation <true | false> - Whether or not to reuse the interpolations. If true, it potentially save the compute time.
282: Level: beginner
284: .seealso: `PCHMG`, `PCGAMG`, `PCHMGSetUseSubspaceCoarsening()`, `PCHMGSetCoarseningComponent()`, `PCHMGSetInnerPCType()`
285: @*/
286: PetscErrorCode PCHMGSetReuseInterpolation(PC pc, PetscBool reuse)
287: {
289: PetscUseMethod(pc, "PCHMGSetReuseInterpolation_C", (PC, PetscBool), (pc, reuse));
290: return 0;
291: }
293: static PetscErrorCode PCHMGSetUseSubspaceCoarsening_HMG(PC pc, PetscBool subspace)
294: {
295: PC_MG *mg = (PC_MG *)pc->data;
296: PC_HMG *hmg = (PC_HMG *)mg->innerctx;
298: hmg->subcoarsening = subspace;
299: return 0;
300: }
302: /*@
303: PCHMGSetUseSubspaceCoarsening - Use subspace coarsening in `PCHMG`
305: Logically Collective
307: Input Parameters:
308: + pc - the `PCHMG` context
309: - reuse - `PETSC_TRUE` indicates that `PCHMG` will use the subspace coarsening
311: Options Database Key:
312: . -pc_hmg_use_subspace_coarsening <true | false> - Whether or not to use subspace coarsening (that is, coarsen a submatrix).
314: Level: beginner
316: .seealso: `PCHMG`, `PCHMGSetReuseInterpolation()`, `PCHMGSetCoarseningComponent()`, `PCHMGSetInnerPCType()`
317: @*/
318: PetscErrorCode PCHMGSetUseSubspaceCoarsening(PC pc, PetscBool subspace)
319: {
321: PetscUseMethod(pc, "PCHMGSetUseSubspaceCoarsening_C", (PC, PetscBool), (pc, subspace));
322: return 0;
323: }
325: static PetscErrorCode PCHMGSetInnerPCType_HMG(PC pc, PCType type)
326: {
327: PC_MG *mg = (PC_MG *)pc->data;
328: PC_HMG *hmg = (PC_HMG *)mg->innerctx;
330: PetscStrallocpy(type, &(hmg->innerpctype));
331: return 0;
332: }
334: /*@C
335: PCHMGSetInnerPCType - Set an inner `PC` type
337: Logically Collective
339: Input Parameters:
340: + pc - the `PCHMG` context
341: - type - `PCHYPRE` or `PCGAMG` coarsening algorithm
343: Options Database Key:
344: . -hmg_inner_pc_type <hypre, gamg> - What method is used to coarsen matrix
346: Level: beginner
348: .seealso: `PCHMG`, `PCType`, `PCHMGSetReuseInterpolation()`, `PCHMGSetUseSubspaceCoarsening()`, `PCHMGSetCoarseningComponent()`
349: @*/
350: PetscErrorCode PCHMGSetInnerPCType(PC pc, PCType type)
351: {
353: PetscUseMethod(pc, "PCHMGSetInnerPCType_C", (PC, PCType), (pc, type));
354: return 0;
355: }
357: static PetscErrorCode PCHMGSetCoarseningComponent_HMG(PC pc, PetscInt component)
358: {
359: PC_MG *mg = (PC_MG *)pc->data;
360: PC_HMG *hmg = (PC_HMG *)mg->innerctx;
362: hmg->component = component;
363: return 0;
364: }
366: /*@
367: PCHMGSetCoarseningComponent - Set which component of the PDE is used for the subspace-based coarsening algorithm
369: Logically Collective
371: Input Parameters:
372: + pc - the `PCHMG` context
373: - component - which component `PC` will coarsen
375: Options Database Key:
376: . -pc_hmg_coarsening_component <i> - Which component is chosen for the subspace-based coarsening algorithm
378: Level: beginner
380: .seealso: `PCHMG`, `PCType`, `PCGAMG`, `PCHMGSetReuseInterpolation()`, `PCHMGSetUseSubspaceCoarsening()`, `PCHMGSetInnerPCType()`
381: @*/
382: PetscErrorCode PCHMGSetCoarseningComponent(PC pc, PetscInt component)
383: {
385: PetscUseMethod(pc, "PCHMGSetCoarseningComponent_C", (PC, PetscInt), (pc, component));
386: return 0;
387: }
389: static PetscErrorCode PCHMGUseMatMAIJ_HMG(PC pc, PetscBool usematmaij)
390: {
391: PC_MG *mg = (PC_MG *)pc->data;
392: PC_HMG *hmg = (PC_HMG *)mg->innerctx;
394: hmg->usematmaij = usematmaij;
395: return 0;
396: }
398: /*@
399: PCHMGUseMatMAIJ - Set a flag that indicates if or not to use `MATMAIJ` for the interpolation matrices for saving memory
401: Logically Collective
403: Input Parameters:
404: + pc - the `PCHMG` context
405: - usematmaij - `PETSC_TRUE` (default) to use `MATMAIJ` for interpolations.
407: Options Database Key:
408: . -pc_hmg_use_matmaij - <true | false >
410: Level: beginner
412: .seealso: `PCHMG`, `PCType`, `PCGAMG`
413: @*/
414: PetscErrorCode PCHMGUseMatMAIJ(PC pc, PetscBool usematmaij)
415: {
417: PetscUseMethod(pc, "PCHMGUseMatMAIJ_C", (PC, PetscBool), (pc, usematmaij));
418: return 0;
419: }
421: /*MC
422: PCHMG - For multiple component PDE problems constructs a hierarchy of restriction operators to coarse grid problems using the submatrix of
423: a single component with either `PCHYPRE` or `PCGAMG`. The same restriction operators are used for each of the components of the PDE with `PCMG`
424: resulting in a much more efficient to build and apply preconditioner than using `PCGAMG` on the entire system.
426: Options Database Keys:
427: + -pc_hmg_reuse_interpolation <true | false> - Whether or not to reuse the interpolations for new matrix values. It can potentially save compute time.
428: . -pc_hmg_use_subspace_coarsening <true | false> - Whether or not to use subspace coarsening (that is, coarsen a submatrix).
429: . -hmg_inner_pc_type <hypre, gamg, ...> - What method is used to coarsen matrix
430: - -pc_hmg_use_matmaij <true | false> - Whether or not to use `MATMAIJ` for multicomponent problems for saving memory
432: Level: intermediate
434: Note:
435: `MatSetBlockSize()` must be called on the linear system matrix to set the number of components of the PDE.
437: References:
438: . * - Fande Kong, Yaqi Wang, Derek R Gaston, Cody J Permann, Andrew E Slaughter, Alexander D Lindsay, Richard C Martineau, A highly parallel multilevel
439: Newton-Krylov-Schwarz method with subspace-based coarsening and partition-based balancing for the multigroup neutron transport equations on
440: 3D unstructured meshes, arXiv preprint arXiv:1903.03659, 2019
442: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCMG`, `PCHYPRE`, `PCHMG`, `PCGetCoarseOperators()`, `PCGetInterpolations()`,
443: `PCHMGSetReuseInterpolation()`, `PCHMGSetUseSubspaceCoarsening()`, `PCHMGSetInnerPCType()`
444: M*/
445: PETSC_EXTERN PetscErrorCode PCCreate_HMG(PC pc)
446: {
447: PC_HMG *hmg;
448: PC_MG *mg;
450: /* if type was previously mg; must manually destroy it because call to PCSetType(pc,PCMG) will not destroy it */
451: PetscTryTypeMethod(pc, destroy);
452: pc->data = NULL;
453: PetscFree(((PetscObject)pc)->type_name);
455: PCSetType(pc, PCMG);
456: PetscObjectChangeTypeName((PetscObject)pc, PCHMG);
457: PetscNew(&hmg);
459: mg = (PC_MG *)pc->data;
460: mg->innerctx = hmg;
461: hmg->reuseinterp = PETSC_FALSE;
462: hmg->subcoarsening = PETSC_FALSE;
463: hmg->usematmaij = PETSC_TRUE;
464: hmg->component = 0;
465: hmg->innerpc = NULL;
467: pc->ops->setfromoptions = PCSetFromOptions_HMG;
468: pc->ops->view = PCView_HMG;
469: pc->ops->destroy = PCDestroy_HMG;
470: pc->ops->setup = PCSetUp_HMG;
472: PetscObjectComposeFunction((PetscObject)pc, "PCHMGSetReuseInterpolation_C", PCHMGSetReuseInterpolation_HMG);
473: PetscObjectComposeFunction((PetscObject)pc, "PCHMGSetUseSubspaceCoarsening_C", PCHMGSetUseSubspaceCoarsening_HMG);
474: PetscObjectComposeFunction((PetscObject)pc, "PCHMGSetInnerPCType_C", PCHMGSetInnerPCType_HMG);
475: PetscObjectComposeFunction((PetscObject)pc, "PCHMGSetCoarseningComponent_C", PCHMGSetCoarseningComponent_HMG);
476: PetscObjectComposeFunction((PetscObject)pc, "PCHMGUseMatMAIJ_C", PCHMGUseMatMAIJ_HMG);
477: return 0;
478: }