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