Actual source code: hmg.c

petsc-master 2019-08-22
Report Typos and Errors
  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;
 12:   PetscBool        usematmaij;
 13: } PC_HMG;

 15: PetscErrorCode PCSetFromOptions_HMG(PetscOptionItems*,PC);
 16: PetscErrorCode PCReset_MG(PC);

 18: static PetscErrorCode PCHMGExtractSubMatrix_Private(Mat pmat,Mat *submat,MatReuse reuse,PetscInt blocksize)
 19: {
 20:   IS             isrow;
 22:   PetscInt       rstart,rend;
 23:   MPI_Comm       comm;

 26:   PetscObjectGetComm((PetscObject)pmat,&comm);
 27:   MatGetOwnershipRange(pmat,&rstart,&rend);
 28:   if ((rend-rstart)%blocksize != 0) SETERRQ3(comm,PETSC_ERR_ARG_INCOMP,"Block size %d is inconsisent for [%d, %d) \n",blocksize,rstart,rend);
 29:   ISCreateStride(comm,(rend-rstart)/blocksize,rstart,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:   PetscErrorCode        ierr;
 42:   MPI_Comm              comm;

 45:   PetscObjectGetComm((PetscObject)subinterp,&comm);
 46:   MatGetOwnershipRange(subinterp,&subrstart,&subrend);
 47:   subrowsize = subrend-subrstart;
 48:   rowsize = subrowsize*blocksize;
 49:   PetscCalloc2(rowsize,&d_nnz,rowsize,&o_nnz);
 50:   MatGetOwnershipRangeColumn(subinterp,&subcstart,&subcend);
 51:   subcolsize = subcend - subcstart;
 52:   colsize    = subcolsize*blocksize;
 53:   max_nz = 0;
 54:   for (subrow=subrstart;subrow<subrend;subrow++) {
 55:     MatGetRow(subinterp,subrow,&nz,&idx,NULL);
 56:     if (max_nz<nz) max_nz = nz;
 57:     dnz = 0; onz = 0;
 58:     for (i=0;i<nz;i++) {
 59:       if(idx[i]>=subcstart && idx[i]<subcend) dnz++;
 60:       else onz++;
 61:     }
 62:     for (i=0;i<blocksize;i++) {
 63:       d_nnz[(subrow-subrstart)*blocksize+i] = dnz;
 64:       o_nnz[(subrow-subrstart)*blocksize+i] = onz;
 65:     }
 66:     MatRestoreRow(subinterp,subrow,&nz,&idx,NULL);
 67:   }
 68:   MatCreateAIJ(comm,rowsize,colsize,PETSC_DETERMINE,PETSC_DETERMINE,0,d_nnz,0,o_nnz,interp);
 69:   MatSetOption(*interp,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
 70:   MatSetOption(*interp,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
 71:   MatSetOption(*interp,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
 72:   MatSetFromOptions(*interp);

 74:   MatSetUp(*interp);
 75:   PetscFree2(d_nnz,o_nnz);
 76:   PetscMalloc1(max_nz,&indices);
 77:   for (subrow=subrstart; subrow<subrend; subrow++) {
 78:     MatGetRow(subinterp,subrow,&nz,&idx,&values);
 79:     for (i=0;i<blocksize;i++) {
 80:       row = subrow*blocksize+i;
 81:       for (j=0;j<nz;j++) {
 82:         indices[j] = idx[j]*blocksize+i;
 83:       }
 84:       MatSetValues(*interp,1,&row,nz,indices,values,INSERT_VALUES);
 85:     }
 86:     MatRestoreRow(subinterp,subrow,&nz,&idx,&values);
 87:   }
 88:   PetscFree(indices);
 89:   MatAssemblyBegin(*interp,MAT_FINAL_ASSEMBLY);
 90:   MatAssemblyEnd(*interp,MAT_FINAL_ASSEMBLY);
 91:   return(0);
 92: }

 94: PetscErrorCode PCSetUp_HMG(PC pc)
 95: {
 96:   PetscErrorCode     ierr;
 97:   Mat                PA, submat;
 98:   PC_MG              *mg   = (PC_MG*)pc->data;
 99:   PC_HMG             *hmg   = (PC_HMG*) mg->innerctx;
100:   MPI_Comm           comm;
101:   PetscInt           level;
102:   PetscInt           num_levels;
103:   Mat                *operators,*interpolations;
104:   PetscInt           blocksize;
105:   const char         *prefix;
106:   PCMGGalerkinType   galerkin;

109:   PetscObjectGetComm((PetscObject)pc,&comm);
110:   if (pc->setupcalled) {
111:    if (hmg->reuseinterp) {
112:      /* If we did not use Galerkin in the last call or we have a different sparsity pattern now,
113:       * we have to build from scratch
114:       * */
115:      PCMGGetGalerkin(pc,&galerkin);
116:      if (galerkin == PC_MG_GALERKIN_NONE || pc->flag != SAME_NONZERO_PATTERN) pc->setupcalled = PETSC_FALSE;
117:      PCMGSetGalerkin(pc,PC_MG_GALERKIN_PMAT);
118:      PCSetUp_MG(pc);
119:      return(0);
120:     } else {
121:      PCReset_MG(pc);
122:      pc->setupcalled = PETSC_FALSE;
123:     }
124:   }

126:   /* Create an inner PC (GAMG or HYPRE) */
127:   if (!hmg->innerpc) {
128:     PCCreate(comm,&hmg->innerpc);
129:     /* If users do not set an inner pc type, we need to set a default value */
130:     if (!hmg->innerpctype) {
131:     /* If hypre is available, use hypre, otherwise, use gamg */
132: #if PETSC_HAVE_HYPRE
133:       PetscStrallocpy(PCHYPRE,&(hmg->innerpctype));
134: #else
135:       PetscStrallocpy(PCGAMG,&(hmg->innerpctype));
136: #endif
137:     }
138:     PCSetType(hmg->innerpc,hmg->innerpctype);
139:   }
140:   PCGetOperators(pc,NULL,&PA);
141:   /* Users need to correctly set a block size of matrix in order to use subspace coarsening */
142:   MatGetBlockSize(PA,&blocksize);
143:   if (blocksize<=1) hmg->subcoarsening = PETSC_FALSE;
144:   /* Extract a submatrix for constructing subinterpolations */
145:   if (hmg->subcoarsening) {
146:     PCHMGExtractSubMatrix_Private(PA,&submat,MAT_INITIAL_MATRIX,blocksize);
147:     PA = submat;
148:   }
149:   PCSetOperators(hmg->innerpc,PA,PA);
150:   if (hmg->subcoarsening) {
151:    MatDestroy(&PA);
152:   }
153:   /* Setup inner PC correctly. During this step, matrix will be coarsened */
154:   PCSetUseAmat(hmg->innerpc,PETSC_FALSE);
155:   PetscObjectGetOptionsPrefix((PetscObject)pc,&prefix);
156:   PetscObjectSetOptionsPrefix((PetscObject)hmg->innerpc,prefix);
157:   PetscObjectAppendOptionsPrefix((PetscObject)hmg->innerpc,"hmg_inner_");
158:   PCSetFromOptions(hmg->innerpc);
159:   PCSetUp(hmg->innerpc);

161:   /* Obtain interpolations IN PLACE. For BoomerAMG, (I,J,data) is reused to avoid memory overhead */
162:   PCGetInterpolations(hmg->innerpc,&num_levels,&interpolations);
163:   /* We can reuse the coarse operators when we do the full space coarsening */
164:   if (!hmg->subcoarsening) {
165:     PCGetCoarseOperators(hmg->innerpc,&num_levels,&operators);
166:   }

168:   PCDestroy(&hmg->innerpc);
169:   hmg->innerpc = 0;
170:   PCMGSetLevels_MG(pc,num_levels,NULL);
171:   /* Set coarse matrices and interpolations to PCMG */
172:   for (level=num_levels-1; level>0; level--) {
173:     Mat P=0, pmat=0;
174:     Vec b, x,r;
175:     if (hmg->subcoarsening) {
176:       if (hmg->usematmaij) {
177:         MatCreateMAIJ(interpolations[level-1],blocksize,&P);
178:         MatDestroy(&interpolations[level-1]);
179:       } else {
180:         /* Grow interpolation. In the future, we should use MAIJ */
181:         PCHMGExpandInterpolation_Private(interpolations[level-1],&P,blocksize);
182:         MatDestroy(&interpolations[level-1]);
183:       }
184:     } else {
185:       P = interpolations[level-1];
186:     }
187:     MatCreateVecs(P,&b,&r);
188:     PCMGSetInterpolation(pc,level,P);
189:     PCMGSetRestriction(pc,level,P);
190:     MatDestroy(&P);
191:     /* We reuse the matrices when we do not do subspace coarsening */
192:     if ((level-1)>=0 && !hmg->subcoarsening) {
193:       pmat = operators[level-1];
194:       PCMGSetOperators(pc,level-1,pmat,pmat);
195:       MatDestroy(&pmat);
196:     }
197:     PCMGSetRhs(pc,level-1,b);

199:     PCMGSetR(pc,level,r);
200:     VecDestroy(&r);

202:     VecDuplicate(b,&x);
203:     PCMGSetX(pc,level-1,x);
204:     VecDestroy(&x);
205:     VecDestroy(&b);
206:   }
207:   PetscFree(interpolations);
208:   if (!hmg->subcoarsening) {
209:     PetscFree(operators);
210:   }
211:   /* Turn Galerkin off when we already have coarse operators */
212:   PCMGSetGalerkin(pc,hmg->subcoarsening ? PC_MG_GALERKIN_PMAT:PC_MG_GALERKIN_NONE);
213:   PCSetDM(pc,NULL);
214:   PCSetUseAmat(pc,PETSC_FALSE);
215:   PetscObjectOptionsBegin((PetscObject)pc);
216:   PCSetFromOptions_MG(PetscOptionsObject,pc); /* should be called in PCSetFromOptions_HMG(), but cannot be called prior to PCMGSetLevels() */
217:   PetscOptionsEnd();
218:   PCSetUp_MG(pc);
219:   return(0);
220: }

222: PetscErrorCode PCDestroy_HMG(PC pc)
223: {
225:   PC_MG          *mg  = (PC_MG*)pc->data;
226:   PC_HMG         *hmg = (PC_HMG*) mg->innerctx;

229:   PCDestroy(&hmg->innerpc);
230:   PetscFree(hmg->innerpctype);
231:   PetscFree(hmg);
232:   PCDestroy_MG(pc);

234:   PetscObjectComposeFunction((PetscObject)pc,"PCHMGSetReuseInterpolation_C",NULL);
235:   PetscObjectComposeFunction((PetscObject)pc,"PCHMGSetUseSubspaceCoarsening_C",NULL);
236:   PetscObjectComposeFunction((PetscObject)pc,"PCHMGSetInnerPCType_C",NULL);
237:   return(0);
238: }

240: PetscErrorCode PCView_HMG(PC pc,PetscViewer viewer)
241: {
242:   PC_MG          *mg = (PC_MG*)pc->data;
243:   PC_HMG         *hmg = (PC_HMG*) mg->innerctx;
245:   PetscBool      iascii;

248:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
249:   if (iascii) {
250:     PetscViewerASCIIPrintf(viewer," Reuse interpolation: %s\n",hmg->reuseinterp? "true":"false");
251:     PetscViewerASCIIPrintf(viewer," Use subspace coarsening: %s\n",hmg->subcoarsening? "true":"false");
252:     PetscViewerASCIIPrintf(viewer," Inner PC type: %s \n",hmg->innerpctype);
253:   }
254:   PCView_MG(pc,viewer);
255:   return(0);
256: }

258: PetscErrorCode PCSetFromOptions_HMG(PetscOptionItems *PetscOptionsObject,PC pc)
259: {
261:   PC_MG          *mg = (PC_MG*)pc->data;
262:   PC_HMG         *hmg = (PC_HMG*) mg->innerctx;

265:   PetscOptionsHead(PetscOptionsObject,"HMG");
266:   PetscOptionsBool("-pc_hmg_reuse_interpolation","Reuse the interpolation operators when possible (cheaper, weaker when matrix entries change a lot)","None",hmg->reuseinterp,&hmg->reuseinterp,NULL);
267:   PetscOptionsBool("-pc_hmg_use_subspace_coarsening","Use the subspace coarsening to compute the interpolations","None",hmg->subcoarsening,&hmg->subcoarsening,NULL);
268:   PetscOptionsBool("-pc_hmg_use_matmaij","Use MatMAIJ store interpolation for saving memory","None",hmg->usematmaij,&hmg->usematmaij,NULL);
269:   PetscOptionsTail();
270:   return(0);
271: }

273: static PetscErrorCode PCHMGSetReuseInterpolation_HMG(PC pc, PetscBool reuse)
274: {
275:   PC_MG          *mg  = (PC_MG*)pc->data;
276:   PC_HMG         *hmg = (PC_HMG*) mg->innerctx;

279:   hmg->reuseinterp = reuse;
280:   return(0);
281: }

283: /*MC
284:    PCHMGSetReuseInterpolation - Reuse interpolation matrices in HMG

286:    Logically Collective on PC

288:    Input Parameters:
289: +  pc - the HMG context
290: -  reuse - True indicates that HMG will reuse the interpolations

292:    Options Database Keys:
293: +  -pc_hmg_reuse_interpolation <true | false> - Whether or not to reuse the interpolations. If true, it potentially save the compute time.

295:    Level: beginner

297: .keywords: HMG, multigrid, interpolation, reuse, set

299: .seealso: PCHMG
300: M*/
301: PetscErrorCode PCHMGSetReuseInterpolation(PC pc, PetscBool reuse)
302: {

307:   PetscUseMethod(pc,"PCHMGSetReuseInterpolation_C",(PC,PetscBool),(pc,reuse));
308:   return(0);
309: }

311: static PetscErrorCode PCHMGSetUseSubspaceCoarsening_HMG(PC pc, PetscBool subspace)
312: {
313:   PC_MG          *mg  = (PC_MG*)pc->data;
314:   PC_HMG         *hmg = (PC_HMG*) mg->innerctx;

317:   hmg->subcoarsening = subspace;
318:   return(0);
319: }

321: /*MC
322:    PCHMGSetUseSubspaceCoarsening - Use subspace coarsening in HMG

324:    Logically Collective on PC

326:    Input Parameters:
327: +  pc - the HMG context
328: -  reuse - True indicates that HMG will use the subspace coarsening

330:    Options Database Keys:
331: +  -pc_hmg_use_subspace_coarsening  <true | false> - Whether or not to use subspace coarsening (that is, coarsen a submatrix).

333:    Level: beginner

335: .keywords: HMG, multigrid, interpolation, subspace, coarsening

337: .seealso: PCHMG
338: M*/
339: PetscErrorCode PCHMGSetUseSubspaceCoarsening(PC pc, PetscBool subspace)
340: {

345:   PetscUseMethod(pc,"PCHMGSetUseSubspaceCoarsening_C",(PC,PetscBool),(pc,subspace));
346:   return(0);
347: }

349: static PetscErrorCode PCHMGSetInnerPCType_HMG(PC pc, PCType type)
350: {
351:   PC_MG           *mg  = (PC_MG*)pc->data;
352:   PC_HMG          *hmg = (PC_HMG*) mg->innerctx;
353:   PetscErrorCode  ierr;

356:   PetscStrallocpy(type,&(hmg->innerpctype));
357:   return(0);
358: }

360: /*MC
361:    PCHMGSetInnerPCType - Set an inner PC type

363:    Logically Collective on PC

365:    Input Parameters:
366: +  pc - the HMG context
367: -  type - <hypre, gamg> coarsening algorithm

369:    Options Database Keys:
370: +  -hmg_inner_pc_type <hypre, gamg> - What method is used to coarsen matrix

372:    Level: beginner

374: .keywords: HMG, multigrid, interpolation, coarsening

376: .seealso: PCHMG, PCType
377: M*/
378: PetscErrorCode PCHMGSetInnerPCType(PC pc, PCType type)
379: {
380:   PetscErrorCode  ierr;

384:   PetscUseMethod(pc,"PCHMGSetInnerPCType_C",(PC,PCType),(pc,type));
385:   return(0);
386: }

388: /*MC
389:    PCHMG - Hybrid of PETSc preconditioners (such as ASM, BJacobi, SOR, etc.) and Hypre BoomerAMG, GAMG or other multilevel methods. BoomerAMG, GAMG
390:            or other multilevel methods is used to coarsen matrix and generate a sequence of coarse matrices and interpolations. The matrices and
391:            interpolations are employed to construct PCMG, and then any available PETSc preconditioners can be chosen as smoothers and the coarse solver.

393:    Options Database Keys:
394: +  -pc_hmg_reuse_interpolation <true | false> - Whether or not to reuse the interpolations. If true, it potentially save the compute time.
395: .  -pc_hmg_use_subspace_coarsening  <true | false> - Whether or not to use subspace coarsening (that is, coarsen a submatrix).
396: .  -hmg_inner_pc_type <hypre, gamg, ...> - What method is used to coarsen matrix
397: -  -pc_hmg_use_matmaij <true | false> - Whether or not to use MatMAIJ for multicomponent problems for saving memory


400:    Notes:
401:     For multicomponent problems, we can just coarsen one submatrix associated with one particular component. In this way, the preconditioner setup
402:     time is significantly reduced. One typical use case is neutron transport equations. There are many variables on each mesh vertex due to the
403:     of angle and energy. Each variable, in fact, corresponds to the same PDEs but with different material properties.

405:    Level: beginner

407:    Concepts: Hybrid of ASM and MG, Subspace Coarsening

409:     References:
410: +   1. - Fande Kong, Yaqi Wang, Derek R Gaston, Cody J Permann, Andrew E Slaughter, Alexander D Lindsay, Richard C Martineau, A highly parallel multilevel
411:     Newton-Krylov-Schwarz method with subspace-based coarsening and partition-based balancing for the multigroup neutron transport equations on
412:     3D unstructured meshes, arXiv preprint arXiv:1903.03659, 2019

414: .seealso:  PCCreate(), PCSetType(), PCType, PC, PCMG, PCHYPRE, PCHMG, PCGetCoarseOperators(), PCGetInterpolations(), PCHMGSetReuseInterpolation(), PCHMGSetUseSubspaceCoarsening(),
415:            PCHMGSetInnerPCType()

417: M*/
418: PETSC_EXTERN PetscErrorCode PCCreate_HMG(PC pc)
419: {
421:   PC_HMG         *hmg;
422:   PC_MG          *mg;

425:   /* if type was previously mg; must manually destroy it because call to PCSetType(pc,PCMG) will not destroy it */
426:   if (pc->ops->destroy) {
427:      (*pc->ops->destroy)(pc);
428:     pc->data = 0;
429:   }
430:   PetscFree(((PetscObject)pc)->type_name);

432:   PCSetType(pc,PCMG);
433:   PetscObjectChangeTypeName((PetscObject)pc, PCHMG);
434:   PetscNew(&hmg);

436:   mg                      = (PC_MG*) pc->data;
437:   mg->innerctx            = hmg;
438:   hmg->reuseinterp        = PETSC_FALSE;
439:   hmg->subcoarsening      = PETSC_FALSE;
440:   hmg->usematmaij         = PETSC_TRUE;
441:   hmg->innerpc            = NULL;

443:   pc->ops->setfromoptions = PCSetFromOptions_HMG;
444:   pc->ops->view           = PCView_HMG;
445:   pc->ops->destroy        = PCDestroy_HMG;
446:   pc->ops->setup          = PCSetUp_HMG;

448:   PetscObjectComposeFunction((PetscObject)pc,"PCHMGSetReuseInterpolation_C",PCHMGSetReuseInterpolation_HMG);
449:   PetscObjectComposeFunction((PetscObject)pc,"PCHMGSetUseSubspaceCoarsening_C",PCHMGSetUseSubspaceCoarsening_HMG);
450:   PetscObjectComposeFunction((PetscObject)pc,"PCHMGSetInnerPCType_C",PCHMGSetInnerPCType_HMG);
451:   return(0);
452: }