1 2 #include <petsc/private/matimpl.h> 3 #include <petsc/private/pcimpl.h> 4 #include <petsc/private/dmimpl.h> 5 #include <petscksp.h> /*I "petscksp.h" I*/ 6 #include <petscdm.h> 7 #include <petscdmda.h> 8 #include <petscdmshell.h> 9 10 #include "../src/ksp/pc/impls/telescope/telescope.h" 11 12 static PetscBool cited = PETSC_FALSE; 13 static const char citation[] = 14 "@inproceedings{MaySananRuppKnepleySmith2016,\n" 15 " title = {Extreme-Scale Multigrid Components within PETSc},\n" 16 " author = {Dave A. May and Patrick Sanan and Karl Rupp and Matthew G. Knepley and Barry F. Smith},\n" 17 " booktitle = {Proceedings of the Platform for Advanced Scientific Computing Conference},\n" 18 " series = {PASC '16},\n" 19 " isbn = {978-1-4503-4126-4},\n" 20 " location = {Lausanne, Switzerland},\n" 21 " pages = {5:1--5:12},\n" 22 " articleno = {5},\n" 23 " numpages = {12},\n" 24 " url = {https://doi.acm.org/10.1145/2929908.2929913},\n" 25 " doi = {10.1145/2929908.2929913},\n" 26 " acmid = {2929913},\n" 27 " publisher = {ACM},\n" 28 " address = {New York, NY, USA},\n" 29 " keywords = {GPU, HPC, agglomeration, coarse-level solver, multigrid, parallel computing, preconditioning},\n" 30 " year = {2016}\n" 31 "}\n"; 32 33 typedef struct { 34 DM dm_fine,dm_coarse; /* these DM's should be topologically identical but use different communicators */ 35 Mat permutation; 36 Vec xp; 37 PetscErrorCode (*fp_dm_field_scatter)(DM,Vec,ScatterMode,DM,Vec); 38 PetscErrorCode (*fp_dm_state_scatter)(DM,ScatterMode,DM); 39 void *dmksp_context_determined; 40 void *dmksp_context_user; 41 } PC_Telescope_CoarseDMCtx; 42 43 PetscErrorCode PCTelescopeSetUp_scatters_CoarseDM(PC pc,PC_Telescope sred,PC_Telescope_CoarseDMCtx *ctx) 44 { 45 Vec xred,yred,xtmp,x,xp; 46 VecScatter scatter; 47 IS isin; 48 Mat B; 49 PetscInt m,bs,st,ed; 50 MPI_Comm comm; 51 52 PetscFunctionBegin; 53 PetscCall(PetscObjectGetComm((PetscObject)pc,&comm)); 54 PetscCall(PCGetOperators(pc,NULL,&B)); 55 PetscCall(MatCreateVecs(B,&x,NULL)); 56 PetscCall(MatGetBlockSize(B,&bs)); 57 PetscCall(VecDuplicate(x,&xp)); 58 m = 0; 59 xred = NULL; 60 yred = NULL; 61 if (PCTelescope_isActiveRank(sred)) { 62 PetscCall(DMCreateGlobalVector(ctx->dm_coarse,&xred)); 63 PetscCall(VecDuplicate(xred,&yred)); 64 PetscCall(VecGetOwnershipRange(xred,&st,&ed)); 65 PetscCall(ISCreateStride(comm,ed-st,st,1,&isin)); 66 PetscCall(VecGetLocalSize(xred,&m)); 67 } else { 68 PetscCall(VecGetOwnershipRange(x,&st,&ed)); 69 PetscCall(ISCreateStride(comm,0,st,1,&isin)); 70 } 71 PetscCall(ISSetBlockSize(isin,bs)); 72 PetscCall(VecCreate(comm,&xtmp)); 73 PetscCall(VecSetSizes(xtmp,m,PETSC_DECIDE)); 74 PetscCall(VecSetBlockSize(xtmp,bs)); 75 PetscCall(VecSetType(xtmp,((PetscObject)x)->type_name)); 76 PetscCall(VecScatterCreate(x,isin,xtmp,NULL,&scatter)); 77 sred->xred = xred; 78 sred->yred = yred; 79 sred->isin = isin; 80 sred->scatter = scatter; 81 sred->xtmp = xtmp; 82 ctx->xp = xp; 83 PetscCall(VecDestroy(&x)); 84 PetscFunctionReturn(0); 85 } 86 87 PetscErrorCode PCTelescopeSetUp_CoarseDM(PC pc,PC_Telescope sred) 88 { 89 PC_Telescope_CoarseDMCtx *ctx; 90 DM dm,dm_coarse = NULL; 91 MPI_Comm comm; 92 PetscBool has_perm,has_kspcomputeoperators,using_kspcomputeoperators; 93 94 PetscFunctionBegin; 95 PetscCall(PetscInfo(pc,"PCTelescope: setup (CoarseDM)\n")); 96 PetscCall(PetscNew(&ctx)); 97 sred->dm_ctx = (void*)ctx; 98 99 PetscCall(PetscObjectGetComm((PetscObject)pc,&comm)); 100 PetscCall(PCGetDM(pc,&dm)); 101 PetscCall(DMGetCoarseDM(dm,&dm_coarse)); 102 ctx->dm_fine = dm; 103 ctx->dm_coarse = dm_coarse; 104 105 /* attach coarse dm to ksp on sub communicator */ 106 if (PCTelescope_isActiveRank(sred)) { 107 PetscCall(KSPSetDM(sred->ksp,ctx->dm_coarse)); 108 if (sred->ignore_kspcomputeoperators) PetscCall(KSPSetDMActive(sred->ksp,PETSC_FALSE)); 109 } 110 111 /* check if there is a method to provide a permutation */ 112 has_perm = PETSC_FALSE; 113 has_kspcomputeoperators = PETSC_FALSE; 114 using_kspcomputeoperators = PETSC_FALSE; 115 116 /* if no permutation is provided, we must rely on KSPSetComputeOperators */ 117 { 118 PetscErrorCode (*dmfine_kspfunc)(KSP,Mat,Mat,void*) = NULL; 119 void *dmfine_kspctx = NULL,*dmcoarse_kspctx = NULL; 120 void *dmfine_appctx = NULL,*dmcoarse_appctx = NULL; 121 void *dmfine_shellctx = NULL,*dmcoarse_shellctx = NULL; 122 123 PetscCall(DMKSPGetComputeOperators(dm,&dmfine_kspfunc,&dmfine_kspctx)); 124 if (dmfine_kspfunc) { has_kspcomputeoperators = PETSC_TRUE; } 125 126 PetscCall(DMGetApplicationContext(ctx->dm_fine,&dmfine_appctx)); 127 PetscCall(DMShellGetContext(ctx->dm_fine,&dmfine_shellctx)); 128 129 /* need to define dmcoarse_kspctx */ 130 if (dmfine_kspfunc && !sred->ignore_kspcomputeoperators) { 131 132 PetscCall(PetscInfo(pc,"PCTelescope: KSPSetComputeOperators fetched from parent DM\n")); 133 if (PCTelescope_isActiveRank(sred)) { 134 PetscCall(DMGetApplicationContext(ctx->dm_coarse,&dmcoarse_appctx)); 135 PetscCall(DMShellGetContext(ctx->dm_coarse,&dmcoarse_shellctx)); 136 } 137 138 /* Assume that if the fine operator didn't require any context, neither will the coarse */ 139 if (!dmfine_kspctx) { 140 dmcoarse_kspctx = NULL; 141 PetscCall(PetscInfo(pc,"PCTelescope: KSPSetComputeOperators using NULL context\n")); 142 } else { 143 144 PetscCall(PetscInfo(pc,"PCTelescope: KSPSetComputeOperators detected non-NULL context from parent DM \n")); 145 if (PCTelescope_isActiveRank(sred)) { 146 147 if (dmfine_kspctx == dmfine_appctx) { 148 dmcoarse_kspctx = dmcoarse_appctx; 149 PetscCall(PetscInfo(pc,"PCTelescope: KSPSetComputeOperators using context from DM->ApplicationContext\n")); 150 PetscCheck(dmcoarse_kspctx,PETSC_COMM_SELF,PETSC_ERR_USER,"Non NULL dmfine->kspctx == dmfine->appctx. NULL dmcoarse->appctx found. Likely this is an error"); 151 } else if (dmfine_kspctx == dmfine_shellctx) { 152 dmcoarse_kspctx = dmcoarse_shellctx; 153 PetscCall(PetscInfo(pc,"PCTelescope: KSPSetComputeOperators using context from DMShell->Context\n")); 154 PetscCheck(dmcoarse_kspctx,PETSC_COMM_SELF,PETSC_ERR_USER,"Non NULL dmfine->kspctx == dmfine.shell->ctx. NULL dmcoarse.shell->ctx found. Likely this is an error"); 155 } 156 ctx->dmksp_context_determined = dmcoarse_kspctx; 157 158 /* look for user provided method to fetch the context */ 159 { 160 PetscErrorCode (*fp_get_coarsedm_context)(DM,void**) = NULL; 161 void *dmcoarse_context_user = NULL; 162 char dmcoarse_method[PETSC_MAX_PATH_LEN]; 163 164 PetscCall(PetscSNPrintf(dmcoarse_method,sizeof(dmcoarse_method),"PCTelescopeGetCoarseDMKSPContext")); 165 PetscCall(PetscObjectQueryFunction((PetscObject)ctx->dm_coarse,dmcoarse_method,&fp_get_coarsedm_context)); 166 if (fp_get_coarsedm_context) { 167 PetscCall(PetscInfo(pc,"PCTelescope: Found composed method PCTelescopeGetCoarseDMKSPContext from coarse DM\n")); 168 PetscCall(fp_get_coarsedm_context(ctx->dm_coarse,&dmcoarse_context_user)); 169 ctx->dmksp_context_user = dmcoarse_context_user; 170 dmcoarse_kspctx = dmcoarse_context_user; 171 } else { 172 PetscCall(PetscInfo(pc,"PCTelescope: Failed to find composed method PCTelescopeGetCoarseDMKSPContext from coarse DM\n")); 173 } 174 } 175 176 if (!dmcoarse_kspctx) { 177 PetscCall(PetscInfo(pc,"PCTelescope: KSPSetComputeOperators failed to determine the context to use on sub-communicator\n")); 178 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot determine which context with use for KSPSetComputeOperators() on sub-communicator"); 179 } 180 } 181 } 182 } 183 184 if (dmfine_kspfunc && !sred->ignore_kspcomputeoperators) { 185 using_kspcomputeoperators = PETSC_TRUE; 186 187 if (PCTelescope_isActiveRank(sred)) { 188 /* sub ksp inherits dmksp_func and context provided by user */ 189 PetscCall(KSPSetComputeOperators(sred->ksp,dmfine_kspfunc,dmcoarse_kspctx)); 190 /* PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)ctx->dmrepart)); */ 191 PetscCall(KSPSetDMActive(sred->ksp,PETSC_TRUE)); 192 } 193 } 194 } 195 196 PetscCheck(has_perm || !has_kspcomputeoperators || using_kspcomputeoperators,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"); 197 PetscCheck(has_perm || has_kspcomputeoperators,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"); 198 199 { 200 char dmfine_method[PETSC_MAX_PATH_LEN]; 201 202 PetscCall(PetscSNPrintf(dmfine_method,sizeof(dmfine_method),"PCTelescopeFieldScatter")); 203 PetscCall(PetscObjectQueryFunction((PetscObject)ctx->dm_fine,dmfine_method,&ctx->fp_dm_field_scatter)); 204 205 PetscCall(PetscSNPrintf(dmfine_method,sizeof(dmfine_method),"PCTelescopeStateScatter")); 206 PetscCall(PetscObjectQueryFunction((PetscObject)ctx->dm_fine,dmfine_method,&ctx->fp_dm_state_scatter)); 207 } 208 209 if (ctx->fp_dm_state_scatter) { 210 PetscCall(PetscInfo(pc,"PCTelescope: Found composed method PCTelescopeStateScatter from parent DM\n")); 211 } else { 212 PetscCall(PetscInfo(pc,"PCTelescope: Failed to find composed method PCTelescopeStateScatter from parent DM\n")); 213 } 214 215 if (ctx->fp_dm_field_scatter) { 216 PetscCall(PetscInfo(pc,"PCTelescope: Found composed method PCTelescopeFieldScatter from parent DM\n")); 217 } else { 218 PetscCall(PetscInfo(pc,"PCTelescope: Failed to find composed method PCTelescopeFieldScatter from parent DM\n")); 219 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"); 220 } 221 222 /* PetscCall(PCTelescopeSetUp_permutation_CoarseDM(pc,sred,ctx)); */ 223 PetscCall(PCTelescopeSetUp_scatters_CoarseDM(pc,sred,ctx)); 224 PetscFunctionReturn(0); 225 } 226 227 PetscErrorCode PCApply_Telescope_CoarseDM(PC pc,Vec x,Vec y) 228 { 229 PC_Telescope sred = (PC_Telescope)pc->data; 230 Vec xred,yred; 231 PC_Telescope_CoarseDMCtx *ctx; 232 233 PetscFunctionBegin; 234 ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx; 235 xred = sred->xred; 236 yred = sred->yred; 237 238 PetscCall(PetscCitationsRegister(citation,&cited)); 239 240 if (ctx->fp_dm_state_scatter) PetscCall(ctx->fp_dm_state_scatter(ctx->dm_fine,SCATTER_FORWARD,ctx->dm_coarse)); 241 242 PetscCall(ctx->fp_dm_field_scatter(ctx->dm_fine,x,SCATTER_FORWARD,ctx->dm_coarse,xred)); 243 244 /* solve */ 245 if (PCTelescope_isActiveRank(sred)) { 246 PetscCall(KSPSolve(sred->ksp,xred,yred)); 247 } 248 249 PetscCall(ctx->fp_dm_field_scatter(ctx->dm_fine,y,SCATTER_REVERSE,ctx->dm_coarse,yred)); 250 PetscFunctionReturn(0); 251 } 252 253 PetscErrorCode PCTelescopeSubNullSpaceCreate_CoarseDM(PC pc,PC_Telescope sred,MatNullSpace nullspace,MatNullSpace *sub_nullspace) 254 { 255 PetscBool has_const; 256 PetscInt k,n = 0; 257 const Vec *vecs; 258 Vec *sub_vecs = NULL; 259 MPI_Comm subcomm; 260 PC_Telescope_CoarseDMCtx *ctx; 261 262 PetscFunctionBegin; 263 ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx; 264 subcomm = sred->subcomm; 265 PetscCall(MatNullSpaceGetVecs(nullspace,&has_const,&n,&vecs)); 266 267 if (PCTelescope_isActiveRank(sred)) { 268 /* create new vectors */ 269 if (n) { 270 PetscCall(VecDuplicateVecs(sred->xred,n,&sub_vecs)); 271 } 272 } 273 274 /* copy entries */ 275 for (k=0; k<n; k++) { 276 PetscCall(ctx->fp_dm_field_scatter(ctx->dm_fine,vecs[k],SCATTER_FORWARD,ctx->dm_coarse,sub_vecs[k])); 277 } 278 279 if (PCTelescope_isActiveRank(sred)) { 280 /* create new (near) nullspace for redundant object */ 281 PetscCall(MatNullSpaceCreate(subcomm,has_const,n,sub_vecs,sub_nullspace)); 282 PetscCall(VecDestroyVecs(n,&sub_vecs)); 283 } 284 PetscFunctionReturn(0); 285 } 286 287 PetscErrorCode PCTelescopeMatNullSpaceCreate_CoarseDM(PC pc,PC_Telescope sred,Mat sub_mat) 288 { 289 Mat B; 290 PC_Telescope_CoarseDMCtx *ctx; 291 292 PetscFunctionBegin; 293 ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx; 294 PetscCall(PCGetOperators(pc,NULL,&B)); 295 { 296 MatNullSpace nullspace,sub_nullspace; 297 PetscCall(MatGetNullSpace(B,&nullspace)); 298 if (nullspace) { 299 PetscCall(PetscInfo(pc,"PCTelescope: generating nullspace (CoarseDM)\n")); 300 PetscCall(PCTelescopeSubNullSpaceCreate_CoarseDM(pc,sred,nullspace,&sub_nullspace)); 301 302 /* attach any user nullspace removal methods and contexts */ 303 if (PCTelescope_isActiveRank(sred)) { 304 void *context = NULL; 305 if (nullspace->remove && !nullspace->rmctx) { 306 PetscCall(MatNullSpaceSetFunction(sub_nullspace,nullspace->remove,context)); 307 } else if (nullspace->remove && nullspace->rmctx) { 308 char dmcoarse_method[PETSC_MAX_PATH_LEN]; 309 PetscErrorCode (*fp_get_coarsedm_context)(DM,void**) = NULL; 310 311 PetscCall(PetscSNPrintf(dmcoarse_method,sizeof(dmcoarse_method),"PCTelescopeGetCoarseDMNullSpaceUserContext")); 312 PetscCall(PetscObjectQueryFunction((PetscObject)ctx->dm_coarse,dmcoarse_method,&fp_get_coarsedm_context)); 313 PetscCheck(context,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); 314 PetscCall(MatNullSpaceSetFunction(sub_nullspace,nullspace->remove,context)); 315 } 316 } 317 318 if (PCTelescope_isActiveRank(sred)) { 319 PetscCall(MatSetNullSpace(sub_mat,sub_nullspace)); 320 PetscCall(MatNullSpaceDestroy(&sub_nullspace)); 321 } 322 } 323 } 324 { 325 MatNullSpace nearnullspace,sub_nearnullspace; 326 PetscCall(MatGetNearNullSpace(B,&nearnullspace)); 327 if (nearnullspace) { 328 PetscCall(PetscInfo(pc,"PCTelescope: generating near nullspace (CoarseDM)\n")); 329 PetscCall(PCTelescopeSubNullSpaceCreate_CoarseDM(pc,sred,nearnullspace,&sub_nearnullspace)); 330 331 /* attach any user nullspace removal methods and contexts */ 332 if (PCTelescope_isActiveRank(sred)) { 333 void *context = NULL; 334 if (nearnullspace->remove && !nearnullspace->rmctx) { 335 PetscCall(MatNullSpaceSetFunction(sub_nearnullspace,nearnullspace->remove,context)); 336 } else if (nearnullspace->remove && nearnullspace->rmctx) { 337 char dmcoarse_method[PETSC_MAX_PATH_LEN]; 338 PetscErrorCode (*fp_get_coarsedm_context)(DM,void**) = NULL; 339 340 PetscCall(PetscSNPrintf(dmcoarse_method,sizeof(dmcoarse_method),"PCTelescopeGetCoarseDMNearNullSpaceUserContext")); 341 PetscCall(PetscObjectQueryFunction((PetscObject)ctx->dm_coarse,dmcoarse_method,&fp_get_coarsedm_context)); 342 PetscCheck(context,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); 343 PetscCall(MatNullSpaceSetFunction(sub_nearnullspace,nearnullspace->remove,context)); 344 } 345 } 346 347 if (PCTelescope_isActiveRank(sred)) { 348 PetscCall(MatSetNearNullSpace(sub_mat,sub_nearnullspace)); 349 PetscCall(MatNullSpaceDestroy(&sub_nearnullspace)); 350 } 351 } 352 } 353 PetscFunctionReturn(0); 354 } 355 356 PetscErrorCode PCReset_Telescope_CoarseDM(PC pc) 357 { 358 PC_Telescope sred = (PC_Telescope)pc->data; 359 PC_Telescope_CoarseDMCtx *ctx; 360 361 PetscFunctionBegin; 362 ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx; 363 ctx->dm_fine = NULL; /* since I did not increment the ref counter we set these to NULL */ 364 ctx->dm_coarse = NULL; /* since I did not increment the ref counter we set these to NULL */ 365 ctx->permutation = NULL; /* this will be fetched from the dm so no need to call destroy */ 366 PetscCall(VecDestroy(&ctx->xp)); 367 ctx->fp_dm_field_scatter = NULL; 368 ctx->fp_dm_state_scatter = NULL; 369 ctx->dmksp_context_determined = NULL; 370 ctx->dmksp_context_user = NULL; 371 PetscFunctionReturn(0); 372 } 373 374 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) 375 { 376 PC_Telescope sred = (PC_Telescope)pc->data; 377 Vec yred = NULL; 378 PetscBool default_init_guess_value = PETSC_FALSE; 379 PC_Telescope_CoarseDMCtx *ctx; 380 381 PetscFunctionBegin; 382 ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx; 383 yred = sred->yred; 384 385 PetscCheck(its <= 1,PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PCApplyRichardson_Telescope_CoarseDM only supports max_it = 1"); 386 *reason = (PCRichardsonConvergedReason)0; 387 388 if (!zeroguess) { 389 PetscCall(PetscInfo(pc,"PCTelescopeCoarseDM: Scattering y for non-zero-initial guess\n")); 390 391 PetscCall(ctx->fp_dm_field_scatter(ctx->dm_fine,y,SCATTER_FORWARD,ctx->dm_coarse,yred)); 392 } 393 394 if (PCTelescope_isActiveRank(sred)) { 395 PetscCall(KSPGetInitialGuessNonzero(sred->ksp,&default_init_guess_value)); 396 if (!zeroguess) PetscCall(KSPSetInitialGuessNonzero(sred->ksp,PETSC_TRUE)); 397 } 398 399 PetscCall(PCApply_Telescope_CoarseDM(pc,x,y)); 400 401 if (PCTelescope_isActiveRank(sred)) { 402 PetscCall(KSPSetInitialGuessNonzero(sred->ksp,default_init_guess_value)); 403 } 404 405 if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS; 406 *outits = 1; 407 PetscFunctionReturn(0); 408 } 409