xref: /petsc/src/ksp/pc/impls/telescope/telescope_coarsedm.c (revision a8d69d7b0124b1e6ce75950a93e6ff079980e86f)
18d9f7141SDave May 
28d9f7141SDave May #include <petsc/private/matimpl.h>
38d9f7141SDave May #include <petsc/private/pcimpl.h>
48d9f7141SDave May #include <petsc/private/dmimpl.h>
58d9f7141SDave May #include <petscksp.h>           /*I "petscksp.h" I*/
68d9f7141SDave May #include <petscdm.h>
78d9f7141SDave May #include <petscdmda.h>
88d9f7141SDave May #include <petscdmshell.h>
98d9f7141SDave May 
108d9f7141SDave May #include "../src/ksp/pc/impls/telescope/telescope.h"
118d9f7141SDave May 
128d9f7141SDave May static PetscBool  cited = PETSC_FALSE;
138d9f7141SDave May static const char citation[] =
148d9f7141SDave May "@inproceedings{MaySananRuppKnepleySmith2016,\n"
158d9f7141SDave May "  title     = {Extreme-Scale Multigrid Components within PETSc},\n"
168d9f7141SDave May "  author    = {Dave A. May and Patrick Sanan and Karl Rupp and Matthew G. Knepley and Barry F. Smith},\n"
178d9f7141SDave May "  booktitle = {Proceedings of the Platform for Advanced Scientific Computing Conference},\n"
188d9f7141SDave May "  series    = {PASC '16},\n"
198d9f7141SDave May "  isbn      = {978-1-4503-4126-4},\n"
208d9f7141SDave May "  location  = {Lausanne, Switzerland},\n"
218d9f7141SDave May "  pages     = {5:1--5:12},\n"
228d9f7141SDave May "  articleno = {5},\n"
238d9f7141SDave May "  numpages  = {12},\n"
24*a8d69d7bSBarry Smith "  url       = {https://doi.acm.org/10.1145/2929908.2929913},\n"
258d9f7141SDave May "  doi       = {10.1145/2929908.2929913},\n"
268d9f7141SDave May "  acmid     = {2929913},\n"
278d9f7141SDave May "  publisher = {ACM},\n"
288d9f7141SDave May "  address   = {New York, NY, USA},\n"
298d9f7141SDave May "  keywords  = {GPU, HPC, agglomeration, coarse-level solver, multigrid, parallel computing, preconditioning},\n"
308d9f7141SDave May "  year      = {2016}\n"
318d9f7141SDave May "}\n";
328d9f7141SDave May 
338d9f7141SDave May typedef struct {
348d9f7141SDave May   DM              dm_fine,dm_coarse; /* these DM's should be topologically identical but use different communicators */
358d9f7141SDave May   Mat             permutation;
368d9f7141SDave May   Vec             xp;
378d9f7141SDave May   PetscErrorCode  (*fp_dm_field_scatter)(DM,Vec,ScatterMode,DM,Vec);
388d9f7141SDave May   PetscErrorCode  (*fp_dm_state_scatter)(DM,ScatterMode,DM);
398d9f7141SDave May   void            *dmksp_context_determined;
408d9f7141SDave May   void            *dmksp_context_user;
418d9f7141SDave May } PC_Telescope_CoarseDMCtx;
428d9f7141SDave May 
438d9f7141SDave May 
448d9f7141SDave May PetscErrorCode PCTelescopeSetUp_scatters_CoarseDM(PC pc,PC_Telescope sred,PC_Telescope_CoarseDMCtx *ctx)
458d9f7141SDave May {
468d9f7141SDave May   PetscErrorCode ierr;
478d9f7141SDave May   Vec            xred,yred,xtmp,x,xp;
488d9f7141SDave May   VecScatter     scatter;
498d9f7141SDave May   IS             isin;
508d9f7141SDave May   Mat            B;
518d9f7141SDave May   PetscInt       m,bs,st,ed;
528d9f7141SDave May   MPI_Comm       comm;
538d9f7141SDave May 
548d9f7141SDave May   PetscFunctionBegin;
558d9f7141SDave May   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
568d9f7141SDave May   ierr = PCGetOperators(pc,NULL,&B);CHKERRQ(ierr);
578d9f7141SDave May   ierr = MatCreateVecs(B,&x,NULL);CHKERRQ(ierr);
588d9f7141SDave May   ierr = MatGetBlockSize(B,&bs);CHKERRQ(ierr);
598d9f7141SDave May   ierr = VecDuplicate(x,&xp);CHKERRQ(ierr);
608d9f7141SDave May   m = 0;
618d9f7141SDave May   xred = NULL;
628d9f7141SDave May   yred = NULL;
6357f12427SDave May   if (PCTelescope_isActiveRank(sred)) {
648d9f7141SDave May     ierr = DMCreateGlobalVector(ctx->dm_coarse,&xred);CHKERRQ(ierr);
658d9f7141SDave May     ierr = VecDuplicate(xred,&yred);CHKERRQ(ierr);
668d9f7141SDave May     ierr = VecGetOwnershipRange(xred,&st,&ed);CHKERRQ(ierr);
678d9f7141SDave May     ierr = ISCreateStride(comm,ed-st,st,1,&isin);CHKERRQ(ierr);
688d9f7141SDave May     ierr = VecGetLocalSize(xred,&m);CHKERRQ(ierr);
698d9f7141SDave May   } else {
708d9f7141SDave May     ierr = VecGetOwnershipRange(x,&st,&ed);CHKERRQ(ierr);
718d9f7141SDave May     ierr = ISCreateStride(comm,0,st,1,&isin);CHKERRQ(ierr);
728d9f7141SDave May   }
738d9f7141SDave May   ierr = ISSetBlockSize(isin,bs);CHKERRQ(ierr);
748d9f7141SDave May   ierr = VecCreate(comm,&xtmp);CHKERRQ(ierr);
758d9f7141SDave May   ierr = VecSetSizes(xtmp,m,PETSC_DECIDE);CHKERRQ(ierr);
768d9f7141SDave May   ierr = VecSetBlockSize(xtmp,bs);CHKERRQ(ierr);
778d9f7141SDave May   ierr = VecSetType(xtmp,((PetscObject)x)->type_name);CHKERRQ(ierr);
789448b7f1SJunchao Zhang   ierr = VecScatterCreate(x,isin,xtmp,NULL,&scatter);CHKERRQ(ierr);
798d9f7141SDave May   sred->xred    = xred;
808d9f7141SDave May   sred->yred    = yred;
818d9f7141SDave May   sred->isin    = isin;
828d9f7141SDave May   sred->scatter = scatter;
838d9f7141SDave May   sred->xtmp    = xtmp;
848d9f7141SDave May   ctx->xp       = xp;
858d9f7141SDave May   ierr = VecDestroy(&x);CHKERRQ(ierr);
868d9f7141SDave May   PetscFunctionReturn(0);
878d9f7141SDave May }
888d9f7141SDave May 
898d9f7141SDave May PetscErrorCode PCTelescopeSetUp_CoarseDM(PC pc,PC_Telescope sred)
908d9f7141SDave May {
918d9f7141SDave May   PC_Telescope_CoarseDMCtx *ctx;
928d9f7141SDave May   DM                       dm,dm_coarse = NULL;
938d9f7141SDave May   MPI_Comm                 comm;
948d9f7141SDave May   PetscBool                has_perm,has_kspcomputeoperators,using_kspcomputeoperators;
958d9f7141SDave May   PetscErrorCode           ierr;
968d9f7141SDave May 
978d9f7141SDave May   PetscFunctionBegin;
988d9f7141SDave May   ierr = PetscInfo(pc,"PCTelescope: setup (CoarseDM)\n");CHKERRQ(ierr);
998d9f7141SDave May   ierr = PetscMalloc1(1,&ctx);CHKERRQ(ierr);
1008d9f7141SDave May   ierr = PetscMemzero(ctx,sizeof(PC_Telescope_CoarseDMCtx));CHKERRQ(ierr);
1018d9f7141SDave May   sred->dm_ctx = (void*)ctx;
1028d9f7141SDave May 
1038d9f7141SDave May   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
1048d9f7141SDave May   ierr = PCGetDM(pc,&dm);CHKERRQ(ierr);
1058d9f7141SDave May   ierr = DMGetCoarseDM(dm,&dm_coarse);CHKERRQ(ierr);
1068d9f7141SDave May   ctx->dm_fine   = dm;
1078d9f7141SDave May   ctx->dm_coarse = dm_coarse;
1088d9f7141SDave May 
1098d9f7141SDave May   /* attach coarse dm to ksp on sub communicator */
11057f12427SDave May   if (PCTelescope_isActiveRank(sred)) {
1118d9f7141SDave May     ierr = KSPSetDM(sred->ksp,ctx->dm_coarse);CHKERRQ(ierr);
1128d9f7141SDave May     if (sred->ignore_kspcomputeoperators) {
1138d9f7141SDave May       ierr = KSPSetDMActive(sred->ksp,PETSC_FALSE);CHKERRQ(ierr);
1148d9f7141SDave May     }
1158d9f7141SDave May   }
1168d9f7141SDave May 
1178d9f7141SDave May   /* check if there is a method to provide a permutation */
1188d9f7141SDave May   has_perm = PETSC_FALSE;
1198d9f7141SDave May   has_kspcomputeoperators = PETSC_FALSE;
1208d9f7141SDave May   using_kspcomputeoperators = PETSC_FALSE;
1218d9f7141SDave May 
1228d9f7141SDave May   /* if no permutation is provided, we must rely on KSPSetComputeOperators */
1238d9f7141SDave May   {
12457f12427SDave May     PetscErrorCode (*dmfine_kspfunc)(KSP,Mat,Mat,void*) = NULL;
1258d9f7141SDave May     void           *dmfine_kspctx = NULL,*dmcoarse_kspctx = NULL;
1268d9f7141SDave May     void           *dmfine_appctx = NULL,*dmcoarse_appctx = NULL;
1278d9f7141SDave May     void           *dmfine_shellctx = NULL,*dmcoarse_shellctx = NULL;
1288d9f7141SDave May 
1298d9f7141SDave May     ierr = DMKSPGetComputeOperators(dm,&dmfine_kspfunc,&dmfine_kspctx);CHKERRQ(ierr);
1308d9f7141SDave May     if (dmfine_kspfunc) { has_kspcomputeoperators = PETSC_TRUE; }
1318d9f7141SDave May 
1328d9f7141SDave May     ierr = DMGetApplicationContext(ctx->dm_fine,&dmfine_appctx);CHKERRQ(ierr);
1338d9f7141SDave May     ierr = DMShellGetContext(ctx->dm_fine,&dmfine_shellctx);CHKERRQ(ierr);
1348d9f7141SDave May 
1358d9f7141SDave May     /* need to define dmcoarse_kspctx */
1368d9f7141SDave May     if (dmfine_kspfunc && !sred->ignore_kspcomputeoperators) {
1378d9f7141SDave May 
1388d9f7141SDave May       ierr = PetscInfo(pc,"PCTelescope: KSPSetComputeOperators fetched from parent DM\n");CHKERRQ(ierr);
13957f12427SDave May       if (PCTelescope_isActiveRank(sred)) {
1408d9f7141SDave May         ierr = DMGetApplicationContext(ctx->dm_coarse,&dmcoarse_appctx);CHKERRQ(ierr);
1418d9f7141SDave May         ierr = DMShellGetContext(ctx->dm_coarse,&dmcoarse_shellctx);CHKERRQ(ierr);
1428d9f7141SDave May       }
1438d9f7141SDave May 
1448d9f7141SDave May       /* Assume that if the fine operator didn't require any context, neither will the coarse */
1458d9f7141SDave May       if (!dmfine_kspctx) {
1468d9f7141SDave May         dmcoarse_kspctx = NULL;
1478d9f7141SDave May         ierr = PetscInfo(pc,"PCTelescope: KSPSetComputeOperators using NULL context\n");CHKERRQ(ierr);
1488d9f7141SDave May       } else {
1498d9f7141SDave May 
1508d9f7141SDave May         ierr = PetscInfo(pc,"PCTelescope: KSPSetComputeOperators detected non-NULL context from parent DM \n");CHKERRQ(ierr);
15157f12427SDave May         if (PCTelescope_isActiveRank(sred)) {
1528d9f7141SDave May 
1538d9f7141SDave May           if (dmfine_kspctx == dmfine_appctx) {
1548d9f7141SDave May             dmcoarse_kspctx = dmcoarse_appctx;
1558d9f7141SDave May             ierr = PetscInfo(pc,"PCTelescope: KSPSetComputeOperators using context from DM->ApplicationContext\n");CHKERRQ(ierr);
1568d9f7141SDave May             if (!dmcoarse_kspctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Non NULL dmfine->kspctx == dmfine->appctx. NULL dmcoarse->appctx found. Likely this is an error");
1578d9f7141SDave May           } else if (dmfine_kspctx == dmfine_shellctx) {
1588d9f7141SDave May             dmcoarse_kspctx = dmcoarse_shellctx;
1598d9f7141SDave May             ierr = PetscInfo(pc,"PCTelescope: KSPSetComputeOperators using context from DMShell->Context\n");CHKERRQ(ierr);
1608d9f7141SDave May             if (!dmcoarse_kspctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Non NULL dmfine->kspctx == dmfine.shell->ctx. NULL dmcoarse.shell->ctx found. Likely this is an error");
1618d9f7141SDave May           }
1628d9f7141SDave May           ctx->dmksp_context_determined = dmcoarse_kspctx;
1638d9f7141SDave May 
1648d9f7141SDave May           /* look for user provided method to fetch the context */
1658d9f7141SDave May           {
1668d9f7141SDave May             PetscErrorCode (*fp_get_coarsedm_context)(DM,void**) = NULL;
1678d9f7141SDave May             void *dmcoarse_context_user = NULL;
1688d9f7141SDave May             char dmcoarse_method[PETSC_MAX_PATH_LEN];
1698d9f7141SDave May 
1708d9f7141SDave May             ierr = PetscSNPrintf(dmcoarse_method,sizeof(dmcoarse_method),"PCTelescopeGetCoarseDMKSPContext");CHKERRQ(ierr);
1718d9f7141SDave May             ierr = PetscObjectQueryFunction((PetscObject)ctx->dm_coarse,dmcoarse_method,&fp_get_coarsedm_context);CHKERRQ(ierr);
1728d9f7141SDave May             if (fp_get_coarsedm_context) {
1738d9f7141SDave May               ierr = PetscInfo(pc,"PCTelescope: Found composed method PCTelescopeGetCoarseDMKSPContext from coarse DM\n");CHKERRQ(ierr);
1748d9f7141SDave May               ierr = fp_get_coarsedm_context(ctx->dm_coarse,&dmcoarse_context_user);CHKERRQ(ierr);
1758d9f7141SDave May               ctx->dmksp_context_user = dmcoarse_context_user;
1768d9f7141SDave May               dmcoarse_kspctx = dmcoarse_context_user;
1778d9f7141SDave May             } else {
1788d9f7141SDave May               ierr = PetscInfo(pc,"PCTelescope: Failed to find composed method PCTelescopeGetCoarseDMKSPContext from coarse DM\n");CHKERRQ(ierr);
1798d9f7141SDave May             }
1808d9f7141SDave May           }
1818d9f7141SDave May 
1828d9f7141SDave May           if (!dmcoarse_kspctx) {
1838d9f7141SDave May             ierr = PetscInfo(pc,"PCTelescope: KSPSetComputeOperators failed to determine the context to use on sub-communicator\n");CHKERRQ(ierr);
1848d9f7141SDave May             SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot determine which context with use for KSPSetComputeOperators() on sub-communicator");
1858d9f7141SDave May           }
1868d9f7141SDave May         }
1878d9f7141SDave May       }
1888d9f7141SDave May     }
1898d9f7141SDave May 
1908d9f7141SDave May     if (dmfine_kspfunc && !sred->ignore_kspcomputeoperators) {
1918d9f7141SDave May       using_kspcomputeoperators = PETSC_TRUE;
1928d9f7141SDave May 
19357f12427SDave May       if (PCTelescope_isActiveRank(sred)) {
1948d9f7141SDave May         /* sub ksp inherits dmksp_func and context provided by user */
1958d9f7141SDave May         ierr = KSPSetComputeOperators(sred->ksp,dmfine_kspfunc,dmcoarse_kspctx);CHKERRQ(ierr);
196c5083d92SDave May         /*ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)ctx->dmrepart);CHKERRQ(ierr);*/
1978d9f7141SDave May         ierr = KSPSetDMActive(sred->ksp,PETSC_TRUE);CHKERRQ(ierr);
1988d9f7141SDave May       }
1998d9f7141SDave May     }
2008d9f7141SDave May   }
2018d9f7141SDave May 
2028d9f7141SDave May   if (!has_perm && has_kspcomputeoperators && !using_kspcomputeoperators) SETERRQ(comm,PETSC_ERR_SUP,"No method to permute an operator was found on the parent DM. A method for KSPSetComputeOperators() was provided but it was requested to be ignored. Telescope setup cannot proceed");
2038d9f7141SDave May   if (!has_perm && !has_kspcomputeoperators) SETERRQ(comm,PETSC_ERR_SUP,"No method to permute an operator was found on the parent DM. No method for KSPSetComputeOperators() was provided. Telescope setup cannot proceed");
2048d9f7141SDave May 
2058d9f7141SDave May   {
2068d9f7141SDave May     char dmfine_method[PETSC_MAX_PATH_LEN];
2078d9f7141SDave May 
2088d9f7141SDave May     ierr = PetscSNPrintf(dmfine_method,sizeof(dmfine_method),"PCTelescopeFieldScatter");CHKERRQ(ierr);
2098d9f7141SDave May     ierr = PetscObjectQueryFunction((PetscObject)ctx->dm_fine,dmfine_method,&ctx->fp_dm_field_scatter);CHKERRQ(ierr);
2108d9f7141SDave May 
2118d9f7141SDave May     ierr = PetscSNPrintf(dmfine_method,sizeof(dmfine_method),"PCTelescopeStateScatter");CHKERRQ(ierr);
2128d9f7141SDave May     ierr = PetscObjectQueryFunction((PetscObject)ctx->dm_fine,dmfine_method,&ctx->fp_dm_state_scatter);CHKERRQ(ierr);
2138d9f7141SDave May   }
2148d9f7141SDave May 
2158d9f7141SDave May   if (ctx->fp_dm_state_scatter) {
2168d9f7141SDave May     ierr = PetscInfo(pc,"PCTelescope: Found composed method PCTelescopeStateScatter from parent DM\n");CHKERRQ(ierr);
2178d9f7141SDave May   } else {
2188d9f7141SDave May     ierr = PetscInfo(pc,"PCTelescope: Failed to find composed method PCTelescopeStateScatter from parent DM\n");CHKERRQ(ierr);
2198d9f7141SDave May   }
2208d9f7141SDave May 
2218d9f7141SDave May   if (ctx->fp_dm_field_scatter) {
2228d9f7141SDave May     ierr = PetscInfo(pc,"PCTelescope: Found composed method PCTelescopeFieldScatter from parent DM\n");CHKERRQ(ierr);
2238d9f7141SDave May   } else {
2248d9f7141SDave May     ierr = PetscInfo(pc,"PCTelescope: Failed to find composed method PCTelescopeFieldScatter from parent DM\n");CHKERRQ(ierr);
2258d9f7141SDave May     SETERRQ(comm,PETSC_ERR_SUP,"No method to scatter fields between the parent DM and coarse DM was found. Must call PetscObjectComposeFunction() with the parent DM. Telescope setup cannot proceed");
2268d9f7141SDave May   }
2278d9f7141SDave May 
2288d9f7141SDave May   /*ierr = PCTelescopeSetUp_permutation_CoarseDM(pc,sred,ctx);CHKERRQ(ierr);*/
2298d9f7141SDave May   ierr = PCTelescopeSetUp_scatters_CoarseDM(pc,sred,ctx);CHKERRQ(ierr);
2308d9f7141SDave May   PetscFunctionReturn(0);
2318d9f7141SDave May }
2328d9f7141SDave May 
2338d9f7141SDave May PetscErrorCode PCApply_Telescope_CoarseDM(PC pc,Vec x,Vec y)
2348d9f7141SDave May {
2358d9f7141SDave May   PC_Telescope             sred = (PC_Telescope)pc->data;
2368d9f7141SDave May   PetscErrorCode           ierr;
2378d9f7141SDave May   Vec                      xred,yred;
2388d9f7141SDave May   PC_Telescope_CoarseDMCtx *ctx;
2398d9f7141SDave May 
2408d9f7141SDave May   PetscFunctionBegin;
2418d9f7141SDave May   ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx;
2428d9f7141SDave May   xred = sred->xred;
2438d9f7141SDave May   yred = sred->yred;
2448d9f7141SDave May 
2458d9f7141SDave May   ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr);
2468d9f7141SDave May 
2478d9f7141SDave May   if (ctx->fp_dm_state_scatter) {
2488d9f7141SDave May     ierr = ctx->fp_dm_state_scatter(ctx->dm_fine,SCATTER_FORWARD,ctx->dm_coarse);CHKERRQ(ierr);
2498d9f7141SDave May   }
2508d9f7141SDave May 
2518d9f7141SDave May   ierr = ctx->fp_dm_field_scatter(ctx->dm_fine,x,SCATTER_FORWARD,ctx->dm_coarse,xred);CHKERRQ(ierr);
2528d9f7141SDave May 
2538d9f7141SDave May   /* solve */
25457f12427SDave May   if (PCTelescope_isActiveRank(sred)) {
2558d9f7141SDave May     ierr = KSPSolve(sred->ksp,xred,yred);CHKERRQ(ierr);
2568d9f7141SDave May   }
2578d9f7141SDave May 
2588d9f7141SDave May   ierr = ctx->fp_dm_field_scatter(ctx->dm_fine,y,SCATTER_REVERSE,ctx->dm_coarse,yred);CHKERRQ(ierr);
2598d9f7141SDave May   PetscFunctionReturn(0);
2608d9f7141SDave May }
2618d9f7141SDave May 
2628d9f7141SDave May PetscErrorCode PCTelescopeSubNullSpaceCreate_CoarseDM(PC pc,PC_Telescope sred,MatNullSpace nullspace,MatNullSpace *sub_nullspace)
2638d9f7141SDave May {
2648d9f7141SDave May   PetscErrorCode           ierr;
2658d9f7141SDave May   PetscBool                has_const;
2668d9f7141SDave May   PetscInt                 k,n = 0;
2678d9f7141SDave May   const Vec                *vecs;
2688d9f7141SDave May   Vec                      *sub_vecs = NULL;
2698d9f7141SDave May   MPI_Comm                 subcomm;
2708d9f7141SDave May   PC_Telescope_CoarseDMCtx *ctx;
2718d9f7141SDave May 
2728d9f7141SDave May   PetscFunctionBegin;
2738d9f7141SDave May   ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx;
2748d9f7141SDave May   subcomm = sred->subcomm;
2758d9f7141SDave May   ierr = MatNullSpaceGetVecs(nullspace,&has_const,&n,&vecs);CHKERRQ(ierr);
2768d9f7141SDave May 
27757f12427SDave May   if (PCTelescope_isActiveRank(sred)) {
2788d9f7141SDave May     /* create new vectors */
2798d9f7141SDave May     if (n) {
2808d9f7141SDave May       ierr = VecDuplicateVecs(sred->xred,n,&sub_vecs);CHKERRQ(ierr);
2818d9f7141SDave May     }
2828d9f7141SDave May   }
2838d9f7141SDave May 
2848d9f7141SDave May   /* copy entries */
2858d9f7141SDave May   for (k=0; k<n; k++) {
2868d9f7141SDave May     ierr = ctx->fp_dm_field_scatter(ctx->dm_fine,vecs[k],SCATTER_FORWARD,ctx->dm_coarse,sub_vecs[k]);CHKERRQ(ierr);
2878d9f7141SDave May   }
2888d9f7141SDave May 
28957f12427SDave May   if (PCTelescope_isActiveRank(sred)) {
2908d9f7141SDave May     /* create new (near) nullspace for redundant object */
2918d9f7141SDave May     ierr = MatNullSpaceCreate(subcomm,has_const,n,sub_vecs,sub_nullspace);CHKERRQ(ierr);
2928d9f7141SDave May     ierr = VecDestroyVecs(n,&sub_vecs);CHKERRQ(ierr);
2938d9f7141SDave May   }
2948d9f7141SDave May   PetscFunctionReturn(0);
2958d9f7141SDave May }
2968d9f7141SDave May 
2978d9f7141SDave May PetscErrorCode PCTelescopeMatNullSpaceCreate_CoarseDM(PC pc,PC_Telescope sred,Mat sub_mat)
2988d9f7141SDave May {
2998d9f7141SDave May   PetscErrorCode           ierr;
3008d9f7141SDave May   Mat                      B;
3018d9f7141SDave May   PC_Telescope_CoarseDMCtx *ctx;
3028d9f7141SDave May 
3038d9f7141SDave May   PetscFunctionBegin;
3048d9f7141SDave May   ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx;
3058d9f7141SDave May   ierr = PCGetOperators(pc,NULL,&B);CHKERRQ(ierr);
3068d9f7141SDave May   {
3078d9f7141SDave May     MatNullSpace nullspace,sub_nullspace;
3088d9f7141SDave May     ierr = MatGetNullSpace(B,&nullspace);CHKERRQ(ierr);
3098d9f7141SDave May     if (nullspace) {
3108d9f7141SDave May       ierr = PetscInfo(pc,"PCTelescope: generating nullspace (CoarseDM)\n");CHKERRQ(ierr);
3118d9f7141SDave May       ierr = PCTelescopeSubNullSpaceCreate_CoarseDM(pc,sred,nullspace,&sub_nullspace);CHKERRQ(ierr);
3128d9f7141SDave May 
3138d9f7141SDave May       /* attach any user nullspace removal methods and contexts */
31457f12427SDave May       if (PCTelescope_isActiveRank(sred)) {
3158d9f7141SDave May         void *context = NULL;
3168d9f7141SDave May         if (nullspace->remove && !nullspace->rmctx){
3178d9f7141SDave May           ierr = MatNullSpaceSetFunction(sub_nullspace,nullspace->remove,context);CHKERRQ(ierr);
3188d9f7141SDave May         } else if (nullspace->remove && nullspace->rmctx) {
3198d9f7141SDave May           char           dmcoarse_method[PETSC_MAX_PATH_LEN];
3208d9f7141SDave May           PetscErrorCode (*fp_get_coarsedm_context)(DM,void**) = NULL;
3218d9f7141SDave May 
3228d9f7141SDave May           ierr = PetscSNPrintf(dmcoarse_method,sizeof(dmcoarse_method),"PCTelescopeGetCoarseDMNullSpaceUserContext");CHKERRQ(ierr);
3238d9f7141SDave May           ierr = PetscObjectQueryFunction((PetscObject)ctx->dm_coarse,dmcoarse_method,&fp_get_coarsedm_context);CHKERRQ(ierr);
3248d9f7141SDave May           if (!context) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Propagation of user null-space removal method with non-NULL context requires the coarse DM be composed with a function named \"%s\"",dmcoarse_method);
3258d9f7141SDave May           ierr = MatNullSpaceSetFunction(sub_nullspace,nullspace->remove,context);CHKERRQ(ierr);
3268d9f7141SDave May         }
3278d9f7141SDave May       }
3288d9f7141SDave May 
32957f12427SDave May       if (PCTelescope_isActiveRank(sred)) {
3308d9f7141SDave May         ierr = MatSetNullSpace(sub_mat,sub_nullspace);CHKERRQ(ierr);
3318d9f7141SDave May         ierr = MatNullSpaceDestroy(&sub_nullspace);CHKERRQ(ierr);
3328d9f7141SDave May       }
3338d9f7141SDave May     }
3348d9f7141SDave May   }
3358d9f7141SDave May   {
3368d9f7141SDave May     MatNullSpace nearnullspace,sub_nearnullspace;
3378d9f7141SDave May     ierr = MatGetNearNullSpace(B,&nearnullspace);CHKERRQ(ierr);
3388d9f7141SDave May     if (nearnullspace) {
3398d9f7141SDave May       ierr = PetscInfo(pc,"PCTelescope: generating near nullspace (CoarseDM)\n");CHKERRQ(ierr);
3408d9f7141SDave May       ierr = PCTelescopeSubNullSpaceCreate_CoarseDM(pc,sred,nearnullspace,&sub_nearnullspace);CHKERRQ(ierr);
3418d9f7141SDave May 
3428d9f7141SDave May       /* attach any user nullspace removal methods and contexts */
34357f12427SDave May       if (PCTelescope_isActiveRank(sred)) {
3448d9f7141SDave May         void *context = NULL;
3458d9f7141SDave May         if (nearnullspace->remove && !nearnullspace->rmctx){
3468d9f7141SDave May           ierr = MatNullSpaceSetFunction(sub_nearnullspace,nearnullspace->remove,context);CHKERRQ(ierr);
3478d9f7141SDave May         } else if (nearnullspace->remove && nearnullspace->rmctx) {
3488d9f7141SDave May           char           dmcoarse_method[PETSC_MAX_PATH_LEN];
3498d9f7141SDave May           PetscErrorCode (*fp_get_coarsedm_context)(DM,void**) = NULL;
3508d9f7141SDave May 
3518d9f7141SDave May           ierr = PetscSNPrintf(dmcoarse_method,sizeof(dmcoarse_method),"PCTelescopeGetCoarseDMNearNullSpaceUserContext");CHKERRQ(ierr);
3528d9f7141SDave May           ierr = PetscObjectQueryFunction((PetscObject)ctx->dm_coarse,dmcoarse_method,&fp_get_coarsedm_context);CHKERRQ(ierr);
3538d9f7141SDave May           if (!context) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Propagation of user near null-space removal method with non-NULL context requires the coarse DM be composed with a function named \"%s\"",dmcoarse_method);
3548d9f7141SDave May           ierr = MatNullSpaceSetFunction(sub_nearnullspace,nearnullspace->remove,context);CHKERRQ(ierr);
3558d9f7141SDave May         }
3568d9f7141SDave May       }
3578d9f7141SDave May 
35857f12427SDave May       if (PCTelescope_isActiveRank(sred)) {
3598d9f7141SDave May         ierr = MatSetNearNullSpace(sub_mat,sub_nearnullspace);CHKERRQ(ierr);
3608d9f7141SDave May         ierr = MatNullSpaceDestroy(&sub_nearnullspace);CHKERRQ(ierr);
3618d9f7141SDave May       }
3628d9f7141SDave May     }
3638d9f7141SDave May   }
3648d9f7141SDave May   PetscFunctionReturn(0);
3658d9f7141SDave May }
3668d9f7141SDave May 
3678d9f7141SDave May PetscErrorCode PCReset_Telescope_CoarseDM(PC pc)
3688d9f7141SDave May {
3698d9f7141SDave May   PetscErrorCode           ierr;
3708d9f7141SDave May   PC_Telescope             sred = (PC_Telescope)pc->data;
3718d9f7141SDave May   PC_Telescope_CoarseDMCtx *ctx;
3728d9f7141SDave May 
3738d9f7141SDave May   PetscFunctionBegin;
3748d9f7141SDave May   ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx;
3758d9f7141SDave May   ctx->dm_fine = NULL; /* since I did not increment the ref counter we set these to NULL */
3768d9f7141SDave May   ctx->dm_coarse = NULL; /* since I did not increment the ref counter we set these to NULL */
3778d9f7141SDave May   ctx->permutation = NULL; /* this will be fetched from the dm so no need to call destroy */
3788d9f7141SDave May   ierr = VecDestroy(&ctx->xp);CHKERRQ(ierr);
3798d9f7141SDave May   ctx->fp_dm_field_scatter = NULL;
3808d9f7141SDave May   ctx->fp_dm_state_scatter = NULL;
3818d9f7141SDave May   ctx->dmksp_context_determined = NULL;
3828d9f7141SDave May   ctx->dmksp_context_user = NULL;
3838d9f7141SDave May   PetscFunctionReturn(0);
3848d9f7141SDave May }
3858d9f7141SDave May 
3868d9f7141SDave May PetscErrorCode PCApplyRichardson_Telescope_CoarseDM(PC pc,Vec x,Vec y,Vec w,PetscReal rtol,PetscReal abstol,PetscReal dtol,PetscInt its,PetscBool zeroguess,PetscInt *outits,PCRichardsonConvergedReason *reason)
3878d9f7141SDave May {
3888d9f7141SDave May   PC_Telescope             sred = (PC_Telescope)pc->data;
3898d9f7141SDave May   PetscErrorCode           ierr;
3908d9f7141SDave May   Vec                      yred = NULL;
3918d9f7141SDave May   PetscBool                default_init_guess_value = PETSC_FALSE;
3928d9f7141SDave May   PC_Telescope_CoarseDMCtx *ctx;
3938d9f7141SDave May 
39457f12427SDave May   PetscFunctionBegin;
3958d9f7141SDave May   ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx;
3968d9f7141SDave May   yred = sred->yred;
3978d9f7141SDave May 
3988d9f7141SDave May   if (its > 1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PCApplyRichardson_Telescope_CoarseDM only supports max_it = 1");
3998d9f7141SDave May   *reason = (PCRichardsonConvergedReason)0;
4008d9f7141SDave May 
4018d9f7141SDave May   if (!zeroguess) {
4028d9f7141SDave May     ierr = PetscInfo(pc,"PCTelescopeCoarseDM: Scattering y for non-zero-initial guess\n");CHKERRQ(ierr);
4038d9f7141SDave May 
4048d9f7141SDave May     ierr = ctx->fp_dm_field_scatter(ctx->dm_fine,y,SCATTER_FORWARD,ctx->dm_coarse,yred);CHKERRQ(ierr);
4058d9f7141SDave May   }
4068d9f7141SDave May 
40757f12427SDave May   if (PCTelescope_isActiveRank(sred)) {
4088d9f7141SDave May     ierr = KSPGetInitialGuessNonzero(sred->ksp,&default_init_guess_value);CHKERRQ(ierr);
4098d9f7141SDave May     if (!zeroguess) ierr = KSPSetInitialGuessNonzero(sred->ksp,PETSC_TRUE);CHKERRQ(ierr);
4108d9f7141SDave May   }
4118d9f7141SDave May 
4128d9f7141SDave May   ierr = PCApply_Telescope_CoarseDM(pc,x,y);CHKERRQ(ierr);
4138d9f7141SDave May 
41457f12427SDave May   if (PCTelescope_isActiveRank(sred)) {
4158d9f7141SDave May     ierr = KSPSetInitialGuessNonzero(sred->ksp,default_init_guess_value);CHKERRQ(ierr);
4168d9f7141SDave May   }
4178d9f7141SDave May 
4188d9f7141SDave May   if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS;
4198d9f7141SDave May   *outits = 1;
4208d9f7141SDave May   PetscFunctionReturn(0);
4218d9f7141SDave May }
422