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