Actual source code: telescope.c

petsc-3.10.4 2019-02-26
Report Typos and Errors

  2:  #include <petsc/private/matimpl.h>
  3:  #include <petsc/private/pcimpl.h>
  4:  #include <petscksp.h>
  5:  #include <petscdm.h>
  6: #include "../src/ksp/pc/impls/telescope/telescope.h"

  8: static PetscBool  cited = PETSC_FALSE;
  9: static const char citation[] =
 10: "@inproceedings{MaySananRuppKnepleySmith2016,\n"
 11: "  title     = {Extreme-Scale Multigrid Components within PETSc},\n"
 12: "  author    = {Dave A. May and Patrick Sanan and Karl Rupp and Matthew G. Knepley and Barry F. Smith},\n"
 13: "  booktitle = {Proceedings of the Platform for Advanced Scientific Computing Conference},\n"
 14: "  series    = {PASC '16},\n"
 15: "  isbn      = {978-1-4503-4126-4},\n"
 16: "  location  = {Lausanne, Switzerland},\n"
 17: "  pages     = {5:1--5:12},\n"
 18: "  articleno = {5},\n"
 19: "  numpages  = {12},\n"
 20: "  url       = {http://doi.acm.org/10.1145/2929908.2929913},\n"
 21: "  doi       = {10.1145/2929908.2929913},\n"
 22: "  acmid     = {2929913},\n"
 23: "  publisher = {ACM},\n"
 24: "  address   = {New York, NY, USA},\n"
 25: "  keywords  = {GPU, HPC, agglomeration, coarse-level solver, multigrid, parallel computing, preconditioning},\n"
 26: "  year      = {2016}\n"
 27: "}\n";

 29: /*
 30:  PCTelescopeSetUp_default()
 31:  PCTelescopeMatCreate_default()

 33:  default

 35:  // scatter in
 36:  x(comm) -> xtmp(comm)

 38:  xred(subcomm) <- xtmp
 39:  yred(subcomm)

 41:  yred(subcomm) --> xtmp

 43:  // scatter out
 44:  xtmp(comm) -> y(comm)
 45: */

 47: PetscBool isActiveRank(PetscSubcomm scomm)
 48: {
 49:   if (scomm->color == 0) { return PETSC_TRUE; }
 50:   else { return PETSC_FALSE; }
 51: }

 53: DM private_PCTelescopeGetSubDM(PC_Telescope sred)
 54: {
 55:   DM subdm = NULL;

 57:   if (!isActiveRank(sred->psubcomm)) { subdm = NULL; }
 58:   else {
 59:     switch (sred->sr_type) {
 60:     case TELESCOPE_DEFAULT: subdm = NULL;
 61:       break;
 62:     case TELESCOPE_DMDA:    subdm = ((PC_Telescope_DMDACtx*)sred->dm_ctx)->dmrepart;
 63:       break;
 64:     case TELESCOPE_DMPLEX:  subdm = NULL;
 65:       break;
 66:     }
 67:   }
 68:   return(subdm);
 69: }

 71: PetscErrorCode PCTelescopeSetUp_default(PC pc,PC_Telescope sred)
 72: {
 74:   PetscInt       m,M,bs,st,ed;
 75:   Vec            x,xred,yred,xtmp;
 76:   Mat            B;
 77:   MPI_Comm       comm,subcomm;
 78:   VecScatter     scatter;
 79:   IS             isin;

 82:   PetscInfo(pc,"PCTelescope: setup (default)\n");
 83:   comm = PetscSubcommParent(sred->psubcomm);
 84:   subcomm = PetscSubcommChild(sred->psubcomm);

 86:   PCGetOperators(pc,NULL,&B);
 87:   MatGetSize(B,&M,NULL);
 88:   MatGetBlockSize(B,&bs);
 89:   MatCreateVecs(B,&x,NULL);

 91:   xred = NULL;
 92:   m    = 0;
 93:   if (isActiveRank(sred->psubcomm)) {
 94:     VecCreate(subcomm,&xred);
 95:     VecSetSizes(xred,PETSC_DECIDE,M);
 96:     VecSetBlockSize(xred,bs);
 97:     VecSetFromOptions(xred);
 98:     VecGetLocalSize(xred,&m);
 99:   }

101:   yred = NULL;
102:   if (isActiveRank(sred->psubcomm)) {
103:     VecDuplicate(xred,&yred);
104:   }

106:   VecCreate(comm,&xtmp);
107:   VecSetSizes(xtmp,m,PETSC_DECIDE);
108:   VecSetBlockSize(xtmp,bs);
109:   VecSetType(xtmp,((PetscObject)x)->type_name);

111:   if (isActiveRank(sred->psubcomm)) {
112:     VecGetOwnershipRange(xred,&st,&ed);
113:     ISCreateStride(comm,(ed-st),st,1,&isin);
114:   } else {
115:     VecGetOwnershipRange(x,&st,&ed);
116:     ISCreateStride(comm,0,st,1,&isin);
117:   }
118:   ISSetBlockSize(isin,bs);

120:   VecScatterCreate(x,isin,xtmp,NULL,&scatter);

122:   sred->isin    = isin;
123:   sred->scatter = scatter;
124:   sred->xred    = xred;
125:   sred->yred    = yred;
126:   sred->xtmp    = xtmp;
127:   VecDestroy(&x);
128:   return(0);
129: }

131: PetscErrorCode PCTelescopeMatCreate_default(PC pc,PC_Telescope sred,MatReuse reuse,Mat *A)
132: {
134:   MPI_Comm       comm,subcomm;
135:   Mat            Bred,B;
136:   PetscInt       nr,nc;
137:   IS             isrow,iscol;
138:   Mat            Blocal,*_Blocal;

141:   PetscInfo(pc,"PCTelescope: updating the redundant preconditioned operator (default)\n");
142:   PetscObjectGetComm((PetscObject)pc,&comm);
143:   subcomm = PetscSubcommChild(sred->psubcomm);
144:   PCGetOperators(pc,NULL,&B);
145:   MatGetSize(B,&nr,&nc);
146:   isrow = sred->isin;
147:   ISCreateStride(comm,nc,0,1,&iscol);
148:   MatCreateSubMatrices(B,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&_Blocal);
149:   Blocal = *_Blocal;
150:   PetscFree(_Blocal);
151:   Bred = NULL;
152:   if (isActiveRank(sred->psubcomm)) {
153:     PetscInt mm;

155:     if (reuse != MAT_INITIAL_MATRIX) { Bred = *A; }

157:     MatGetSize(Blocal,&mm,NULL);
158:     MatCreateMPIMatConcatenateSeqMat(subcomm,Blocal,mm,reuse,&Bred);
159:   }
160:   *A = Bred;
161:   ISDestroy(&iscol);
162:   MatDestroy(&Blocal);
163:   return(0);
164: }

166: static PetscErrorCode PCTelescopeSubNullSpaceCreate_Telescope(PC pc,PC_Telescope sred,MatNullSpace nullspace,MatNullSpace *sub_nullspace)
167: {
168:   PetscErrorCode   ierr;
169:   PetscBool        has_const;
170:   const Vec        *vecs;
171:   Vec              *sub_vecs = NULL;
172:   PetscInt         i,k,n = 0;
173:   MPI_Comm         subcomm;

176:   subcomm = PetscSubcommChild(sred->psubcomm);
177:   MatNullSpaceGetVecs(nullspace,&has_const,&n,&vecs);

179:   if (isActiveRank(sred->psubcomm)) {
180:     if (n) {
181:       VecDuplicateVecs(sred->xred,n,&sub_vecs);
182:     }
183:   }

185:   /* copy entries */
186:   for (k=0; k<n; k++) {
187:     const PetscScalar *x_array;
188:     PetscScalar       *LA_sub_vec;
189:     PetscInt          st,ed;

191:     /* pull in vector x->xtmp */
192:     VecScatterBegin(sred->scatter,vecs[k],sred->xtmp,INSERT_VALUES,SCATTER_FORWARD);
193:     VecScatterEnd(sred->scatter,vecs[k],sred->xtmp,INSERT_VALUES,SCATTER_FORWARD);
194:     if (sub_vecs) {
195:       /* copy vector entries into xred */
196:       VecGetArrayRead(sred->xtmp,&x_array);
197:       if (sub_vecs[k]) {
198:         VecGetOwnershipRange(sub_vecs[k],&st,&ed);
199:         VecGetArray(sub_vecs[k],&LA_sub_vec);
200:         for (i=0; i<ed-st; i++) {
201:           LA_sub_vec[i] = x_array[i];
202:         }
203:         VecRestoreArray(sub_vecs[k],&LA_sub_vec);
204:       }
205:       VecRestoreArrayRead(sred->xtmp,&x_array);
206:     }
207:   }

209:   if (isActiveRank(sred->psubcomm)) {
210:     /* create new (near) nullspace for redundant object */
211:     MatNullSpaceCreate(subcomm,has_const,n,sub_vecs,sub_nullspace);
212:     VecDestroyVecs(n,&sub_vecs);
213:     if (nullspace->remove) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Propagation of custom remove callbacks not supported when propagating (near) nullspaces with PCTelescope");
214:     if (nullspace->rmctx) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Propagation of custom remove callback context not supported when propagating (near) nullspaces with PCTelescope");
215:   }
216:   return(0);
217: }

219: static PetscErrorCode PCTelescopeMatNullSpaceCreate_default(PC pc,PC_Telescope sred,Mat sub_mat)
220: {
221:   PetscErrorCode   ierr;
222:   Mat              B;

225:   PCGetOperators(pc,NULL,&B);

227:   /* Propagate the nullspace if it exists */
228:   {
229:     MatNullSpace nullspace,sub_nullspace;
230:     MatGetNullSpace(B,&nullspace);
231:     if (nullspace) {
232:       PetscInfo(pc,"PCTelescope: generating nullspace (default)\n");
233:       PCTelescopeSubNullSpaceCreate_Telescope(pc,sred,nullspace,&sub_nullspace);
234:       if (isActiveRank(sred->psubcomm)) {
235:         MatSetNullSpace(sub_mat,sub_nullspace);
236:         MatNullSpaceDestroy(&sub_nullspace);
237:       }
238:     }
239:   }

241:   /* Propagate the near nullspace if it exists */
242:   {
243:     MatNullSpace nearnullspace,sub_nearnullspace;
244:     MatGetNearNullSpace(B,&nearnullspace);
245:     if (nearnullspace) {
246:       PetscInfo(pc,"PCTelescope: generating near nullspace (default)\n");
247:       PCTelescopeSubNullSpaceCreate_Telescope(pc,sred,nearnullspace,&sub_nearnullspace);
248:       if (isActiveRank(sred->psubcomm)) {
249:         MatSetNearNullSpace(sub_mat,sub_nearnullspace);
250:         MatNullSpaceDestroy(&sub_nearnullspace);
251:       }
252:     }
253:   }
254:   return(0);
255: }

257: static PetscErrorCode PCView_Telescope(PC pc,PetscViewer viewer)
258: {
259:   PC_Telescope     sred = (PC_Telescope)pc->data;
260:   PetscErrorCode   ierr;
261:   PetscBool        iascii,isstring;
262:   PetscViewer      subviewer;

265:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
266:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
267:   if (iascii) {
268:     if (!sred->psubcomm) {
269:       PetscViewerASCIIPrintf(viewer,"  preconditioner not yet setup\n");
270:     } else {
271:       MPI_Comm    comm,subcomm;
272:       PetscMPIInt comm_size,subcomm_size;
273:       DM          dm,subdm;

275:       PCGetDM(pc,&dm);
276:       subdm = private_PCTelescopeGetSubDM(sred);
277:       comm = PetscSubcommParent(sred->psubcomm);
278:       subcomm = PetscSubcommChild(sred->psubcomm);
279:       MPI_Comm_size(comm,&comm_size);
280:       MPI_Comm_size(subcomm,&subcomm_size);

282:       PetscViewerASCIIPrintf(viewer,"  parent comm size reduction factor = %D\n",sred->redfactor);
283:       PetscViewerASCIIPrintf(viewer,"  comm_size = %d , subcomm_size = %d\n",(int)comm_size,(int)subcomm_size);
284:       switch (sred->subcommtype) {
285:         case PETSC_SUBCOMM_INTERLACED :
286:           PetscViewerASCIIPrintf(viewer,"  subcomm type: interlaced\n",sred->subcommtype);
287:           break;
288:         case PETSC_SUBCOMM_CONTIGUOUS :
289:           PetscViewerASCIIPrintf(viewer,"  subcomm type: contiguous\n",sred->subcommtype);
290:           break;
291:         default :
292:           SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"General subcomm type not supported by PCTelescope");
293:       }
294:       PetscViewerGetSubViewer(viewer,subcomm,&subviewer);
295:       if (isActiveRank(sred->psubcomm)) {
296:         PetscViewerASCIIPushTab(subviewer);

298:         if (dm && sred->ignore_dm) {
299:           PetscViewerASCIIPrintf(subviewer,"  ignoring DM\n");
300:         }
301:         if (sred->ignore_kspcomputeoperators) {
302:           PetscViewerASCIIPrintf(subviewer,"  ignoring KSPComputeOperators\n");
303:         }
304:         switch (sred->sr_type) {
305:         case TELESCOPE_DEFAULT:
306:           PetscViewerASCIIPrintf(subviewer,"  using default setup\n");
307:           break;
308:         case TELESCOPE_DMDA:
309:           PetscViewerASCIIPrintf(subviewer,"  DMDA detected\n");
310:           DMView_DA_Short(subdm,subviewer);
311:           break;
312:         case TELESCOPE_DMPLEX:
313:           PetscViewerASCIIPrintf(subviewer,"  DMPLEX detected\n");
314:           break;
315:         }

317:         KSPView(sred->ksp,subviewer);
318:         PetscViewerASCIIPopTab(subviewer);
319:       }
320:       PetscViewerRestoreSubViewer(viewer,subcomm,&subviewer);
321:     }
322:   }
323:   return(0);
324: }

326: static PetscErrorCode PCSetUp_Telescope(PC pc)
327: {
328:   PC_Telescope      sred = (PC_Telescope)pc->data;
329:   PetscErrorCode    ierr;
330:   MPI_Comm          comm,subcomm=0;
331:   PCTelescopeType   sr_type;

334:   PetscObjectGetComm((PetscObject)pc,&comm);

336:   /* subcomm definition */
337:   if (!pc->setupcalled) {
338:     if (!sred->psubcomm) {
339:       PetscSubcommCreate(comm,&sred->psubcomm);
340:       PetscSubcommSetNumber(sred->psubcomm,sred->redfactor);
341:       PetscSubcommSetType(sred->psubcomm,sred->subcommtype);
342:       PetscLogObjectMemory((PetscObject)pc,sizeof(PetscSubcomm));
343:     }
344:   }
345:   subcomm = PetscSubcommChild(sred->psubcomm);

347:   /* internal KSP */
348:   if (!pc->setupcalled) {
349:     const char *prefix;

351:     if (isActiveRank(sred->psubcomm)) {
352:       KSPCreate(subcomm,&sred->ksp);
353:       KSPSetErrorIfNotConverged(sred->ksp,pc->erroriffailure);
354:       PetscObjectIncrementTabLevel((PetscObject)sred->ksp,(PetscObject)pc,1);
355:       PetscLogObjectParent((PetscObject)pc,(PetscObject)sred->ksp);
356:       KSPSetType(sred->ksp,KSPPREONLY);
357:       PCGetOptionsPrefix(pc,&prefix);
358:       KSPSetOptionsPrefix(sred->ksp,prefix);
359:       KSPAppendOptionsPrefix(sred->ksp,"telescope_");
360:     }
361:   }
362:   /* Determine type of setup/update */
363:   if (!pc->setupcalled) {
364:     PetscBool has_dm,same;
365:     DM        dm;

367:     sr_type = TELESCOPE_DEFAULT;
368:     has_dm = PETSC_FALSE;
369:     PCGetDM(pc,&dm);
370:     if (dm) { has_dm = PETSC_TRUE; }
371:     if (has_dm) {
372:       /* check for dmda */
373:       PetscObjectTypeCompare((PetscObject)dm,DMDA,&same);
374:       if (same) {
375:         PetscInfo(pc,"PCTelescope: found DMDA\n");
376:         sr_type = TELESCOPE_DMDA;
377:       }
378:       /* check for dmplex */
379:       PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&same);
380:       if (same) {
381:         PetscInfo(pc,"PCTelescope: found DMPLEX\n");
382:         sr_type = TELESCOPE_DMPLEX;
383:       }
384:     }

386:     if (sred->ignore_dm) {
387:       PetscInfo(pc,"PCTelescope: ignore DM\n");
388:       sr_type = TELESCOPE_DEFAULT;
389:     }
390:     sred->sr_type = sr_type;
391:   } else {
392:     sr_type = sred->sr_type;
393:   }

395:   /* set function pointers for repartition setup, matrix creation/update, matrix (near) nullspace, and reset functionality */
396:   switch (sr_type) {
397:   case TELESCOPE_DEFAULT:
398:     sred->pctelescope_setup_type              = PCTelescopeSetUp_default;
399:     sred->pctelescope_matcreate_type          = PCTelescopeMatCreate_default;
400:     sred->pctelescope_matnullspacecreate_type = PCTelescopeMatNullSpaceCreate_default;
401:     sred->pctelescope_reset_type              = NULL;
402:     break;
403:   case TELESCOPE_DMDA:
404:     pc->ops->apply                            = PCApply_Telescope_dmda;
405:     pc->ops->applyrichardson                  = PCApplyRichardson_Telescope_dmda;
406:     sred->pctelescope_setup_type              = PCTelescopeSetUp_dmda;
407:     sred->pctelescope_matcreate_type          = PCTelescopeMatCreate_dmda;
408:     sred->pctelescope_matnullspacecreate_type = PCTelescopeMatNullSpaceCreate_dmda;
409:     sred->pctelescope_reset_type              = PCReset_Telescope_dmda;
410:     break;
411:   case TELESCOPE_DMPLEX: SETERRQ(comm,PETSC_ERR_SUP,"Support for DMPLEX is currently not available");
412:     break;
413:   default: SETERRQ(comm,PETSC_ERR_SUP,"Only support for repartitioning DMDA is provided");
414:     break;
415:   }

417:   /* setup */
418:   if (sred->pctelescope_setup_type) {
419:     sred->pctelescope_setup_type(pc,sred);
420:   }
421:   /* update */
422:   if (!pc->setupcalled) {
423:     if (sred->pctelescope_matcreate_type) {
424:       sred->pctelescope_matcreate_type(pc,sred,MAT_INITIAL_MATRIX,&sred->Bred);
425:     }
426:     if (sred->pctelescope_matnullspacecreate_type) {
427:       sred->pctelescope_matnullspacecreate_type(pc,sred,sred->Bred);
428:     }
429:   } else {
430:     if (sred->pctelescope_matcreate_type) {
431:       sred->pctelescope_matcreate_type(pc,sred,MAT_REUSE_MATRIX,&sred->Bred);
432:     }
433:   }

435:   /* common - no construction */
436:   if (isActiveRank(sred->psubcomm)) {
437:     KSPSetOperators(sred->ksp,sred->Bred,sred->Bred);
438:     if (pc->setfromoptionscalled && !pc->setupcalled){
439:       KSPSetFromOptions(sred->ksp);
440:     }
441:   }
442:   return(0);
443: }

445: static PetscErrorCode PCApply_Telescope(PC pc,Vec x,Vec y)
446: {
447:   PC_Telescope      sred = (PC_Telescope)pc->data;
448:   PetscErrorCode    ierr;
449:   Vec               xtmp,xred,yred;
450:   PetscInt          i,st,ed;
451:   VecScatter        scatter;
452:   PetscScalar       *array;
453:   const PetscScalar *x_array;

456:   PetscCitationsRegister(citation,&cited);

458:   xtmp    = sred->xtmp;
459:   scatter = sred->scatter;
460:   xred    = sred->xred;
461:   yred    = sred->yred;

463:   /* pull in vector x->xtmp */
464:   VecScatterBegin(scatter,x,xtmp,INSERT_VALUES,SCATTER_FORWARD);
465:   VecScatterEnd(scatter,x,xtmp,INSERT_VALUES,SCATTER_FORWARD);

467:   /* copy vector entries into xred */
468:   VecGetArrayRead(xtmp,&x_array);
469:   if (xred) {
470:     PetscScalar *LA_xred;
471:     VecGetOwnershipRange(xred,&st,&ed);
472:     VecGetArray(xred,&LA_xred);
473:     for (i=0; i<ed-st; i++) {
474:       LA_xred[i] = x_array[i];
475:     }
476:     VecRestoreArray(xred,&LA_xred);
477:   }
478:   VecRestoreArrayRead(xtmp,&x_array);
479:   /* solve */
480:   if (isActiveRank(sred->psubcomm)) {
481:     KSPSolve(sred->ksp,xred,yred);
482:   }
483:   /* return vector */
484:   VecGetArray(xtmp,&array);
485:   if (yred) {
486:     const PetscScalar *LA_yred;
487:     VecGetOwnershipRange(yred,&st,&ed);
488:     VecGetArrayRead(yred,&LA_yred);
489:     for (i=0; i<ed-st; i++) {
490:       array[i] = LA_yred[i];
491:     }
492:     VecRestoreArrayRead(yred,&LA_yred);
493:   }
494:   VecRestoreArray(xtmp,&array);
495:   VecScatterBegin(scatter,xtmp,y,INSERT_VALUES,SCATTER_REVERSE);
496:   VecScatterEnd(scatter,xtmp,y,INSERT_VALUES,SCATTER_REVERSE);
497:   return(0);
498: }

500: static PetscErrorCode PCApplyRichardson_Telescope(PC pc,Vec x,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool zeroguess,PetscInt *outits,PCRichardsonConvergedReason *reason)
501: {
502:   PC_Telescope      sred = (PC_Telescope)pc->data;
503:   PetscErrorCode    ierr;
504:   Vec               xtmp,yred;
505:   PetscInt          i,st,ed;
506:   VecScatter        scatter;
507:   const PetscScalar *x_array;
508:   PetscBool         default_init_guess_value;

511:   xtmp    = sred->xtmp;
512:   scatter = sred->scatter;
513:   yred    = sred->yred;

515:   if (its > 1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PCApplyRichardson_Telescope only supports max_it = 1");
516:   *reason = (PCRichardsonConvergedReason)0;

518:   if (!zeroguess) {
519:     PetscInfo(pc,"PCTelescope: Scattering y for non-zero initial guess\n");
520:     /* pull in vector y->xtmp */
521:     VecScatterBegin(scatter,y,xtmp,INSERT_VALUES,SCATTER_FORWARD);
522:     VecScatterEnd(scatter,y,xtmp,INSERT_VALUES,SCATTER_FORWARD);

524:     /* copy vector entries into xred */
525:     VecGetArrayRead(xtmp,&x_array);
526:     if (yred) {
527:       PetscScalar *LA_yred;
528:       VecGetOwnershipRange(yred,&st,&ed);
529:       VecGetArray(yred,&LA_yred);
530:       for (i=0; i<ed-st; i++) {
531:         LA_yred[i] = x_array[i];
532:       }
533:       VecRestoreArray(yred,&LA_yred);
534:     }
535:     VecRestoreArrayRead(xtmp,&x_array);
536:   }

538:   if (isActiveRank(sred->psubcomm)) {
539:     KSPGetInitialGuessNonzero(sred->ksp,&default_init_guess_value);
540:     if (!zeroguess) KSPSetInitialGuessNonzero(sred->ksp,PETSC_TRUE);
541:   }

543:   PCApply_Telescope(pc,x,y);

545:   if (isActiveRank(sred->psubcomm)) {
546:     KSPSetInitialGuessNonzero(sred->ksp,default_init_guess_value);
547:   }

549:   if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS;
550:   *outits = 1;
551:   return(0);
552: }

554: static PetscErrorCode PCReset_Telescope(PC pc)
555: {
556:   PC_Telescope   sred = (PC_Telescope)pc->data;

559:   ISDestroy(&sred->isin);
560:   VecScatterDestroy(&sred->scatter);
561:   VecDestroy(&sred->xred);
562:   VecDestroy(&sred->yred);
563:   VecDestroy(&sred->xtmp);
564:   MatDestroy(&sred->Bred);
565:   KSPReset(sred->ksp);
566:   if (sred->pctelescope_reset_type) {
567:     sred->pctelescope_reset_type(pc);
568:   }
569:   return(0);
570: }

572: static PetscErrorCode PCDestroy_Telescope(PC pc)
573: {
574:   PC_Telescope     sred = (PC_Telescope)pc->data;
575:   PetscErrorCode   ierr;

578:   PCReset_Telescope(pc);
579:   KSPDestroy(&sred->ksp);
580:   PetscSubcommDestroy(&sred->psubcomm);
581:   PetscFree(sred->dm_ctx);
582:   PetscFree(pc->data);
583:   return(0);
584: }

586: static PetscErrorCode PCSetFromOptions_Telescope(PetscOptionItems *PetscOptionsObject,PC pc)
587: {
588:   PC_Telescope     sred = (PC_Telescope)pc->data;
589:   PetscErrorCode   ierr;
590:   MPI_Comm         comm;
591:   PetscMPIInt      size;
592:   PetscBool        flg;
593:   PetscSubcommType subcommtype;

596:   PetscObjectGetComm((PetscObject)pc,&comm);
597:   MPI_Comm_size(comm,&size);
598:   PetscOptionsHead(PetscOptionsObject,"Telescope options");
599:   PetscOptionsEnum("-pc_telescope_subcomm_type","Subcomm type (interlaced or contiguous)","PCTelescopeSetSubcommType",PetscSubcommTypes,(PetscEnum)sred->subcommtype,(PetscEnum*)&subcommtype,&flg);
600:   if (flg) {
601:     PCTelescopeSetSubcommType(pc,subcommtype);
602:   }
603:   PetscOptionsInt("-pc_telescope_reduction_factor","Factor to reduce comm size by","PCTelescopeSetReductionFactor",sred->redfactor,&sred->redfactor,0);
604:   if (sred->redfactor > size) SETERRQ(comm,PETSC_ERR_ARG_WRONG,"-pc_telescope_reduction_factor <= comm size");
605:   PetscOptionsBool("-pc_telescope_ignore_dm","Ignore any DM attached to the PC","PCTelescopeSetIgnoreDM",sred->ignore_dm,&sred->ignore_dm,0);
606:   PetscOptionsBool("-pc_telescope_ignore_kspcomputeoperators","Ignore method used to compute A","PCTelescopeSetIgnoreKSPComputeOperators",sred->ignore_kspcomputeoperators,&sred->ignore_kspcomputeoperators,0);
607:   PetscOptionsTail();
608:   return(0);
609: }

611: /* PC simplementation specific API's */

613: static PetscErrorCode PCTelescopeGetKSP_Telescope(PC pc,KSP *ksp)
614: {
615:   PC_Telescope red = (PC_Telescope)pc->data;
617:   if (ksp) *ksp = red->ksp;
618:   return(0);
619: }

621: static PetscErrorCode PCTelescopeGetSubcommType_Telescope(PC pc,PetscSubcommType *subcommtype)
622: {
623:   PC_Telescope red = (PC_Telescope)pc->data;
625:   if (subcommtype) *subcommtype = red->subcommtype;
626:   return(0);
627: }

629: static PetscErrorCode PCTelescopeSetSubcommType_Telescope(PC pc,PetscSubcommType subcommtype)
630: {
631:   PC_Telescope     red = (PC_Telescope)pc->data;

634:   if (pc->setupcalled) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"You cannot change the subcommunicator type for PCTelescope after it has been set up.");
635:   red->subcommtype = subcommtype;
636:   return(0);
637: }

639: static PetscErrorCode PCTelescopeGetReductionFactor_Telescope(PC pc,PetscInt *fact)
640: {
641:   PC_Telescope red = (PC_Telescope)pc->data;
643:   if (fact) *fact = red->redfactor;
644:   return(0);
645: }

647: static PetscErrorCode PCTelescopeSetReductionFactor_Telescope(PC pc,PetscInt fact)
648: {
649:   PC_Telescope     red = (PC_Telescope)pc->data;
650:   PetscMPIInt      size;
651:   PetscErrorCode   ierr;

654:   MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
655:   if (fact <= 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Reduction factor of telescoping PC %D must be positive",fact);
656:   if (fact > size) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Reduction factor of telescoping PC %D must be <= comm.size",fact);
657:   red->redfactor = fact;
658:   return(0);
659: }

661: static PetscErrorCode PCTelescopeGetIgnoreDM_Telescope(PC pc,PetscBool *v)
662: {
663:   PC_Telescope red = (PC_Telescope)pc->data;
665:   if (v) *v = red->ignore_dm;
666:   return(0);
667: }

669: static PetscErrorCode PCTelescopeSetIgnoreDM_Telescope(PC pc,PetscBool v)
670: {
671:   PC_Telescope red = (PC_Telescope)pc->data;
673:   red->ignore_dm = v;
674:   return(0);
675: }

677: static PetscErrorCode PCTelescopeGetIgnoreKSPComputeOperators_Telescope(PC pc,PetscBool *v)
678: {
679:   PC_Telescope red = (PC_Telescope)pc->data;
681:   if (v) *v = red->ignore_kspcomputeoperators;
682:   return(0);
683: }

685: static PetscErrorCode PCTelescopeSetIgnoreKSPComputeOperators_Telescope(PC pc,PetscBool v)
686: {
687:   PC_Telescope red = (PC_Telescope)pc->data;
689:   red->ignore_kspcomputeoperators = v;
690:   return(0);
691: }

693: static PetscErrorCode PCTelescopeGetDM_Telescope(PC pc,DM *dm)
694: {
695:   PC_Telescope red = (PC_Telescope)pc->data;
697:   *dm = private_PCTelescopeGetSubDM(red);
698:   return(0);
699: }

701: /*@
702:  PCTelescopeGetKSP - Gets the KSP created by the telescoping PC.

704:  Not Collective

706:  Input Parameter:
707: .  pc - the preconditioner context

709:  Output Parameter:
710: .  subksp - the KSP defined the smaller set of processes

712:  Level: advanced

714: .keywords: PC, telescopting solve
715: @*/
716: PetscErrorCode PCTelescopeGetKSP(PC pc,KSP *subksp)
717: {
720:   PetscUseMethod(pc,"PCTelescopeGetKSP_C",(PC,KSP*),(pc,subksp));
721:   return(0);
722: }

724: /*@
725:  PCTelescopeGetReductionFactor - Gets the factor by which the original number of processes has been reduced by.

727:  Not Collective

729:  Input Parameter:
730: .  pc - the preconditioner context

732:  Output Parameter:
733: .  fact - the reduction factor

735:  Level: advanced

737: .keywords: PC, telescoping solve
738: @*/
739: PetscErrorCode PCTelescopeGetReductionFactor(PC pc,PetscInt *fact)
740: {
743:   PetscUseMethod(pc,"PCTelescopeGetReductionFactor_C",(PC,PetscInt*),(pc,fact));
744:   return(0);
745: }

747: /*@
748:  PCTelescopeSetReductionFactor - Sets the factor by which the original number of processes has been reduced by.

750:  Not Collective

752:  Input Parameter:
753: .  pc - the preconditioner context

755:  Output Parameter:
756: .  fact - the reduction factor

758:  Level: advanced

760: .keywords: PC, telescoping solve
761: @*/
762: PetscErrorCode PCTelescopeSetReductionFactor(PC pc,PetscInt fact)
763: {
766:   PetscTryMethod(pc,"PCTelescopeSetReductionFactor_C",(PC,PetscInt),(pc,fact));
767:   return(0);
768: }

770: /*@
771:  PCTelescopeGetIgnoreDM - Get the flag indicating if any DM attached to the PC will be used.

773:  Not Collective

775:  Input Parameter:
776: .  pc - the preconditioner context

778:  Output Parameter:
779: .  v - the flag

781:  Level: advanced

783: .keywords: PC, telescoping solve
784: @*/
785: PetscErrorCode PCTelescopeGetIgnoreDM(PC pc,PetscBool *v)
786: {
789:   PetscUseMethod(pc,"PCTelescopeGetIgnoreDM_C",(PC,PetscBool*),(pc,v));
790:   return(0);
791: }

793: /*@
794:  PCTelescopeSetIgnoreDM - Set a flag to ignore any DM attached to the PC.

796:  Not Collective

798:  Input Parameter:
799: .  pc - the preconditioner context

801:  Output Parameter:
802: .  v - Use PETSC_TRUE to ignore any DM

804:  Level: advanced

806: .keywords: PC, telescoping solve
807: @*/
808: PetscErrorCode PCTelescopeSetIgnoreDM(PC pc,PetscBool v)
809: {
812:   PetscTryMethod(pc,"PCTelescopeSetIgnoreDM_C",(PC,PetscBool),(pc,v));
813:   return(0);
814: }

816: /*@
817:  PCTelescopeGetIgnoreKSPComputeOperators - Get the flag indicating if KSPComputeOperators will be used.

819:  Not Collective

821:  Input Parameter:
822: .  pc - the preconditioner context

824:  Output Parameter:
825: .  v - the flag

827:  Level: advanced

829: .keywords: PC, telescoping solve
830: @*/
831: PetscErrorCode PCTelescopeGetIgnoreKSPComputeOperators(PC pc,PetscBool *v)
832: {
835:   PetscUseMethod(pc,"PCTelescopeGetIgnoreKSPComputeOperators_C",(PC,PetscBool*),(pc,v));
836:   return(0);
837: }

839: /*@
840:  PCTelescopeSetIgnoreKSPComputeOperators - Set a flag to ignore KSPComputeOperators.

842:  Not Collective

844:  Input Parameter:
845: .  pc - the preconditioner context

847:  Output Parameter:
848: .  v - Use PETSC_TRUE to ignore the method (if defined) set via KSPSetComputeOperators on pc

850:  Level: advanced

852: .keywords: PC, telescoping solve
853: @*/
854: PetscErrorCode PCTelescopeSetIgnoreKSPComputeOperators(PC pc,PetscBool v)
855: {
858:   PetscTryMethod(pc,"PCTelescopeSetIgnoreKSPComputeOperators_C",(PC,PetscBool),(pc,v));
859:   return(0);
860: }

862: /*@
863:  PCTelescopeGetDM - Get the re-partitioned DM attached to the sub KSP.

865:  Not Collective

867:  Input Parameter:
868: .  pc - the preconditioner context

870:  Output Parameter:
871: .  subdm - The re-partitioned DM

873:  Level: advanced

875: .keywords: PC, telescoping solve
876: @*/
877: PetscErrorCode PCTelescopeGetDM(PC pc,DM *subdm)
878: {
881:   PetscUseMethod(pc,"PCTelescopeGetDM_C",(PC,DM*),(pc,subdm));
882:   return(0);
883: }

885: /*@
886:  PCTelescopeSetSubcommType - set subcommunicator type (interlaced or contiguous)

888:  Logically Collective

890:  Input Parameter:
891: +  pc - the preconditioner context
892: -  subcommtype - the subcommunicator type (see PetscSubcommType)

894:  Level: advanced

896: .keywords: PC, telescoping solve

898: .seealso: PetscSubcommType, PetscSubcomm, PCTELESCOPE
899: @*/
900: PetscErrorCode PCTelescopeSetSubcommType(PC pc, PetscSubcommType subcommtype)
901: {
904:   PetscTryMethod(pc,"PCTelescopeSetSubcommType_C",(PC,PetscSubcommType),(pc,subcommtype));
905:   return(0);
906: }

908: /*@
909:  PCTelescopeGetSubcommType - Get the subcommunicator type (interlaced or contiguous)

911:  Not Collective

913:  Input Parameter:
914: .  pc - the preconditioner context

916:  Output Parameter:
917: .  subcommtype - the subcommunicator type (see PetscSubcommType)

919:  Level: advanced

921: .keywords: PC, telescoping solve

923: .seealso: PetscSubcomm, PetscSubcommType, PCTELESCOPE
924: @*/
925: PetscErrorCode PCTelescopeGetSubcommType(PC pc, PetscSubcommType *subcommtype)
926: {
929:   PetscUseMethod(pc,"PCTelescopeGetSubcommType_C",(PC,PetscSubcommType*),(pc,subcommtype));
930:   return(0);
931: }

933: /* -------------------------------------------------------------------------------------*/
934: /*MC
935:    PCTELESCOPE - Runs a KSP solver on a sub-group of processors. MPI processes not in the sub-communicator are idle during the solve.

937:    Options Database:
938: +  -pc_telescope_reduction_factor <r> - factor to use communicator size by. e.g. with 64 MPI processes and r=4, the new sub-communicator will have 64/4 = 16 ranks.
939: -  -pc_telescope_ignore_dm  - flag to indicate whether an attached DM should be ignored
940: -  -pc_telescope_subcomm_type <interlaced,contiguous> - how to define the reduced communicator. see PetscSubcomm for more.

942:    Level: advanced

944:    Notes:
945:    The preconditioner is deemed telescopic as it only calls KSPSolve() on a single
946:    sub-communicator, in contrast with PCREDUNDANT which calls KSPSolve() on N sub-communicators.
947:    This means there will be MPI processes which will be idle during the application of this preconditioner.

949:    The default KSP is PREONLY. If a DM is attached to the PC, it is re-partitioned on the sub-communicator.
950:    Both the Bmat operator and the right hand side vector are permuted into the new DOF ordering defined by the re-partitioned DM.
951:    Currently only support for re-partitioning a DMDA is provided.
952:    Any nullspace attached to the original Bmat operator is extracted, re-partitioned and set on the repartitioned Bmat operator.
953:    KSPSetComputeOperators() is not propagated to the sub KSP.
954:    Currently there is no support for the flag -pc_use_amat

956:    Assuming that the parent preconditioner (PC) is defined on a communicator c, this implementation
957:    creates a child sub-communicator (c') containing fewer MPI processes than the original parent preconditioner (PC).

959:   Developer Notes:
960:    During PCSetup, the B operator is scattered onto c'.
961:    Within PCApply, the RHS vector (x) is scattered into a redundant vector, xred (defined on c').
962:    Then, KSPSolve() is executed on the c' communicator.

964:    The communicator used within the telescoping preconditioner is defined by a PetscSubcomm using the INTERLACED 
965:    creation routine by default (this can be changed with -pc_telescope_subcomm_type). We run the sub KSP on only the ranks within the communicator which have a color equal to zero.

967:    The telescoping preconditioner is aware of nullspaces and near nullspaces which are attached to the B operator.
968:    In the case where B has a (near) nullspace attached, the (near) nullspace vectors are extracted from B and mapped into
969:    a new (near) nullspace, defined on the sub-communicator, which is attached to B' (the B operator which was scattered to c')

971:    The telescoping preconditioner is aware of an attached DM. In the event that the DM is of type DMDA (2D or 3D - 
972:    1D support for 1D DMDAs is not provided), a new DMDA is created on c' (e.g. it is re-partitioned), and this new DM 
973:    is attached the sub KSPSolve(). The design of telescope is such that it should be possible to extend support 
974:    for re-partitioning other to DM's (e.g. DMPLEX). The user can supply a flag to ignore attached DMs.

976:    By default, B' is defined by simply fusing rows from different MPI processes

978:    When a DMDA is attached to the parent preconditioner, B' is defined by: (i) performing a symmetric permutation of B
979:    into the ordering defined by the DMDA on c', (ii) extracting the local chunks via MatCreateSubMatrices(), (iii) fusing the
980:    locally (sequential) matrices defined on the ranks common to c and c' into B' using MatCreateMPIMatConcatenateSeqMat()

982:    Limitations/improvements include the following.
983:    VecPlaceArray() could be used within PCApply() to improve efficiency and reduce memory usage.

985:    The symmetric permutation used when a DMDA is encountered is performed via explicitly assmbleming a permutation matrix P,
986:    and performing P^T.A.P. Possibly it might be more efficient to use MatPermute(). We opted to use P^T.A.P as it appears
987:    VecPermute() does not supported for the use case required here. By computing P, one can permute both the operator and RHS in a 
988:    consistent manner.

990:    Mapping of vectors is performed in the following way.
991:    Suppose the parent comm size was 4, and we set a reduction factor of 2; this would give a comm size on c' of 2.
992:    Using the interlaced creation routine, the ranks in c with color = 0 will be rank 0 and 2.
993:    We perform the scatter to the sub-comm in the following way.
994:    [1] Given a vector x defined on comm c

996:    rank(c) : _________ 0 ______  ________ 1 _______  ________ 2 _____________ ___________ 3 __________
997:          x : [0, 1, 2, 3, 4, 5] [6, 7, 8, 9, 10, 11] [12, 13, 14, 15, 16, 17] [18, 19, 20, 21, 22, 23]

999:    scatter to xtmp defined also on comm c so that we have the following values

1001:    rank(c) : ___________________ 0 ________________  _1_  ______________________ 2 _______________________  __3_
1002:       xtmp : [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] [  ] [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23] [  ]

1004:    The entries on rank 1 and 3 (ranks which do not have a color = 0 in c') have no values


1007:    [2] Copy the value from rank 0, 2 (indices with respect to comm c) into the vector xred which is defined on communicator c'.
1008:    Ranks 0 and 2 are the only ranks in the subcomm which have a color = 0.

1010:     rank(c') : ___________________ 0 _______________  ______________________ 1 _____________________
1011:       xred : [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]


1014:   Contributed by Dave May

1016:   Reference:
1017:   Dave A. May, Patrick Sanan, Karl Rupp, Matthew G. Knepley, and Barry F. Smith, "Extreme-Scale Multigrid Components within PETSc". 2016. In Proceedings of the Platform for Advanced Scientific Computing Conference (PASC '16). DOI: 10.1145/2929908.2929913

1019: .seealso:  PCTelescopeGetKSP(), PCTelescopeGetDM(), PCTelescopeGetReductionFactor(), PCTelescopeSetReductionFactor(), PCTelescopeGetIgnoreDM(), PCTelescopeSetIgnoreDM(), PCREDUNDANT
1020: M*/
1021: PETSC_EXTERN PetscErrorCode PCCreate_Telescope(PC pc)
1022: {
1023:   PetscErrorCode       ierr;
1024:   struct _PC_Telescope *sred;

1027:   PetscNewLog(pc,&sred);
1028:   sred->subcommtype    = PETSC_SUBCOMM_INTERLACED;
1029:   sred->redfactor      = 1;
1030:   sred->ignore_dm      = PETSC_FALSE;
1031:   sred->ignore_kspcomputeoperators = PETSC_FALSE;
1032:   pc->data             = (void*)sred;

1034:   pc->ops->apply           = PCApply_Telescope;
1035:   pc->ops->applytranspose  = NULL;
1036:   pc->ops->applyrichardson = PCApplyRichardson_Telescope;
1037:   pc->ops->setup           = PCSetUp_Telescope;
1038:   pc->ops->destroy         = PCDestroy_Telescope;
1039:   pc->ops->reset           = PCReset_Telescope;
1040:   pc->ops->setfromoptions  = PCSetFromOptions_Telescope;
1041:   pc->ops->view            = PCView_Telescope;

1043:   sred->pctelescope_setup_type              = PCTelescopeSetUp_default;
1044:   sred->pctelescope_matcreate_type          = PCTelescopeMatCreate_default;
1045:   sred->pctelescope_matnullspacecreate_type = PCTelescopeMatNullSpaceCreate_default;
1046:   sred->pctelescope_reset_type              = NULL;

1048:   PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeGetKSP_C",PCTelescopeGetKSP_Telescope);
1049:   PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeGetSubcommType_C",PCTelescopeGetSubcommType_Telescope);
1050:   PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeSetSubcommType_C",PCTelescopeSetSubcommType_Telescope);
1051:   PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeGetReductionFactor_C",PCTelescopeGetReductionFactor_Telescope);
1052:   PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeSetReductionFactor_C",PCTelescopeSetReductionFactor_Telescope);
1053:   PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeGetIgnoreDM_C",PCTelescopeGetIgnoreDM_Telescope);
1054:   PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeSetIgnoreDM_C",PCTelescopeSetIgnoreDM_Telescope);
1055:   PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeGetIgnoreKSPComputeOperators_C",PCTelescopeGetIgnoreKSPComputeOperators_Telescope);
1056:   PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeSetIgnoreKSPComputeOperators_C",PCTelescopeSetIgnoreKSPComputeOperators_Telescope);
1057:   PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeGetDM_C",PCTelescopeGetDM_Telescope);
1058:   return(0);
1059: }