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) { 109 PetscCall(KSPSetDMActive(sred->ksp,PETSC_FALSE)); 110 } 111 } 112 113 /* check if there is a method to provide a permutation */ 114 has_perm = PETSC_FALSE; 115 has_kspcomputeoperators = PETSC_FALSE; 116 using_kspcomputeoperators = PETSC_FALSE; 117 118 /* if no permutation is provided, we must rely on KSPSetComputeOperators */ 119 { 120 PetscErrorCode (*dmfine_kspfunc)(KSP,Mat,Mat,void*) = NULL; 121 void *dmfine_kspctx = NULL,*dmcoarse_kspctx = NULL; 122 void *dmfine_appctx = NULL,*dmcoarse_appctx = NULL; 123 void *dmfine_shellctx = NULL,*dmcoarse_shellctx = NULL; 124 125 PetscCall(DMKSPGetComputeOperators(dm,&dmfine_kspfunc,&dmfine_kspctx)); 126 if (dmfine_kspfunc) { has_kspcomputeoperators = PETSC_TRUE; } 127 128 PetscCall(DMGetApplicationContext(ctx->dm_fine,&dmfine_appctx)); 129 PetscCall(DMShellGetContext(ctx->dm_fine,&dmfine_shellctx)); 130 131 /* need to define dmcoarse_kspctx */ 132 if (dmfine_kspfunc && !sred->ignore_kspcomputeoperators) { 133 134 PetscCall(PetscInfo(pc,"PCTelescope: KSPSetComputeOperators fetched from parent DM\n")); 135 if (PCTelescope_isActiveRank(sred)) { 136 PetscCall(DMGetApplicationContext(ctx->dm_coarse,&dmcoarse_appctx)); 137 PetscCall(DMShellGetContext(ctx->dm_coarse,&dmcoarse_shellctx)); 138 } 139 140 /* Assume that if the fine operator didn't require any context, neither will the coarse */ 141 if (!dmfine_kspctx) { 142 dmcoarse_kspctx = NULL; 143 PetscCall(PetscInfo(pc,"PCTelescope: KSPSetComputeOperators using NULL context\n")); 144 } else { 145 146 PetscCall(PetscInfo(pc,"PCTelescope: KSPSetComputeOperators detected non-NULL context from parent DM \n")); 147 if (PCTelescope_isActiveRank(sred)) { 148 149 if (dmfine_kspctx == dmfine_appctx) { 150 dmcoarse_kspctx = dmcoarse_appctx; 151 PetscCall(PetscInfo(pc,"PCTelescope: KSPSetComputeOperators using context from DM->ApplicationContext\n")); 152 PetscCheck(dmcoarse_kspctx,PETSC_COMM_SELF,PETSC_ERR_USER,"Non NULL dmfine->kspctx == dmfine->appctx. NULL dmcoarse->appctx found. Likely this is an error"); 153 } else if (dmfine_kspctx == dmfine_shellctx) { 154 dmcoarse_kspctx = dmcoarse_shellctx; 155 PetscCall(PetscInfo(pc,"PCTelescope: KSPSetComputeOperators using context from DMShell->Context\n")); 156 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"); 157 } 158 ctx->dmksp_context_determined = dmcoarse_kspctx; 159 160 /* look for user provided method to fetch the context */ 161 { 162 PetscErrorCode (*fp_get_coarsedm_context)(DM,void**) = NULL; 163 void *dmcoarse_context_user = NULL; 164 char dmcoarse_method[PETSC_MAX_PATH_LEN]; 165 166 PetscCall(PetscSNPrintf(dmcoarse_method,sizeof(dmcoarse_method),"PCTelescopeGetCoarseDMKSPContext")); 167 PetscCall(PetscObjectQueryFunction((PetscObject)ctx->dm_coarse,dmcoarse_method,&fp_get_coarsedm_context)); 168 if (fp_get_coarsedm_context) { 169 PetscCall(PetscInfo(pc,"PCTelescope: Found composed method PCTelescopeGetCoarseDMKSPContext from coarse DM\n")); 170 PetscCall(fp_get_coarsedm_context(ctx->dm_coarse,&dmcoarse_context_user)); 171 ctx->dmksp_context_user = dmcoarse_context_user; 172 dmcoarse_kspctx = dmcoarse_context_user; 173 } else { 174 PetscCall(PetscInfo(pc,"PCTelescope: Failed to find composed method PCTelescopeGetCoarseDMKSPContext from coarse DM\n")); 175 } 176 } 177 178 if (!dmcoarse_kspctx) { 179 PetscCall(PetscInfo(pc,"PCTelescope: KSPSetComputeOperators failed to determine the context to use on sub-communicator\n")); 180 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot determine which context with use for KSPSetComputeOperators() on sub-communicator"); 181 } 182 } 183 } 184 } 185 186 if (dmfine_kspfunc && !sred->ignore_kspcomputeoperators) { 187 using_kspcomputeoperators = PETSC_TRUE; 188 189 if (PCTelescope_isActiveRank(sred)) { 190 /* sub ksp inherits dmksp_func and context provided by user */ 191 PetscCall(KSPSetComputeOperators(sred->ksp,dmfine_kspfunc,dmcoarse_kspctx)); 192 /* PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)ctx->dmrepart)); */ 193 PetscCall(KSPSetDMActive(sred->ksp,PETSC_TRUE)); 194 } 195 } 196 } 197 198 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"); 199 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"); 200 201 { 202 char dmfine_method[PETSC_MAX_PATH_LEN]; 203 204 PetscCall(PetscSNPrintf(dmfine_method,sizeof(dmfine_method),"PCTelescopeFieldScatter")); 205 PetscCall(PetscObjectQueryFunction((PetscObject)ctx->dm_fine,dmfine_method,&ctx->fp_dm_field_scatter)); 206 207 PetscCall(PetscSNPrintf(dmfine_method,sizeof(dmfine_method),"PCTelescopeStateScatter")); 208 PetscCall(PetscObjectQueryFunction((PetscObject)ctx->dm_fine,dmfine_method,&ctx->fp_dm_state_scatter)); 209 } 210 211 if (ctx->fp_dm_state_scatter) { 212 PetscCall(PetscInfo(pc,"PCTelescope: Found composed method PCTelescopeStateScatter from parent DM\n")); 213 } else { 214 PetscCall(PetscInfo(pc,"PCTelescope: Failed to find composed method PCTelescopeStateScatter from parent DM\n")); 215 } 216 217 if (ctx->fp_dm_field_scatter) { 218 PetscCall(PetscInfo(pc,"PCTelescope: Found composed method PCTelescopeFieldScatter from parent DM\n")); 219 } else { 220 PetscCall(PetscInfo(pc,"PCTelescope: Failed to find composed method PCTelescopeFieldScatter from parent DM\n")); 221 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"); 222 } 223 224 /* PetscCall(PCTelescopeSetUp_permutation_CoarseDM(pc,sred,ctx)); */ 225 PetscCall(PCTelescopeSetUp_scatters_CoarseDM(pc,sred,ctx)); 226 PetscFunctionReturn(0); 227 } 228 229 PetscErrorCode PCApply_Telescope_CoarseDM(PC pc,Vec x,Vec y) 230 { 231 PC_Telescope sred = (PC_Telescope)pc->data; 232 Vec xred,yred; 233 PC_Telescope_CoarseDMCtx *ctx; 234 235 PetscFunctionBegin; 236 ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx; 237 xred = sred->xred; 238 yred = sred->yred; 239 240 PetscCall(PetscCitationsRegister(citation,&cited)); 241 242 if (ctx->fp_dm_state_scatter) { 243 PetscCall(ctx->fp_dm_state_scatter(ctx->dm_fine,SCATTER_FORWARD,ctx->dm_coarse)); 244 } 245 246 PetscCall(ctx->fp_dm_field_scatter(ctx->dm_fine,x,SCATTER_FORWARD,ctx->dm_coarse,xred)); 247 248 /* solve */ 249 if (PCTelescope_isActiveRank(sred)) { 250 PetscCall(KSPSolve(sred->ksp,xred,yred)); 251 } 252 253 PetscCall(ctx->fp_dm_field_scatter(ctx->dm_fine,y,SCATTER_REVERSE,ctx->dm_coarse,yred)); 254 PetscFunctionReturn(0); 255 } 256 257 PetscErrorCode PCTelescopeSubNullSpaceCreate_CoarseDM(PC pc,PC_Telescope sred,MatNullSpace nullspace,MatNullSpace *sub_nullspace) 258 { 259 PetscBool has_const; 260 PetscInt k,n = 0; 261 const Vec *vecs; 262 Vec *sub_vecs = NULL; 263 MPI_Comm subcomm; 264 PC_Telescope_CoarseDMCtx *ctx; 265 266 PetscFunctionBegin; 267 ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx; 268 subcomm = sred->subcomm; 269 PetscCall(MatNullSpaceGetVecs(nullspace,&has_const,&n,&vecs)); 270 271 if (PCTelescope_isActiveRank(sred)) { 272 /* create new vectors */ 273 if (n) { 274 PetscCall(VecDuplicateVecs(sred->xred,n,&sub_vecs)); 275 } 276 } 277 278 /* copy entries */ 279 for (k=0; k<n; k++) { 280 PetscCall(ctx->fp_dm_field_scatter(ctx->dm_fine,vecs[k],SCATTER_FORWARD,ctx->dm_coarse,sub_vecs[k])); 281 } 282 283 if (PCTelescope_isActiveRank(sred)) { 284 /* create new (near) nullspace for redundant object */ 285 PetscCall(MatNullSpaceCreate(subcomm,has_const,n,sub_vecs,sub_nullspace)); 286 PetscCall(VecDestroyVecs(n,&sub_vecs)); 287 } 288 PetscFunctionReturn(0); 289 } 290 291 PetscErrorCode PCTelescopeMatNullSpaceCreate_CoarseDM(PC pc,PC_Telescope sred,Mat sub_mat) 292 { 293 Mat B; 294 PC_Telescope_CoarseDMCtx *ctx; 295 296 PetscFunctionBegin; 297 ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx; 298 PetscCall(PCGetOperators(pc,NULL,&B)); 299 { 300 MatNullSpace nullspace,sub_nullspace; 301 PetscCall(MatGetNullSpace(B,&nullspace)); 302 if (nullspace) { 303 PetscCall(PetscInfo(pc,"PCTelescope: generating nullspace (CoarseDM)\n")); 304 PetscCall(PCTelescopeSubNullSpaceCreate_CoarseDM(pc,sred,nullspace,&sub_nullspace)); 305 306 /* attach any user nullspace removal methods and contexts */ 307 if (PCTelescope_isActiveRank(sred)) { 308 void *context = NULL; 309 if (nullspace->remove && !nullspace->rmctx) { 310 PetscCall(MatNullSpaceSetFunction(sub_nullspace,nullspace->remove,context)); 311 } else if (nullspace->remove && nullspace->rmctx) { 312 char dmcoarse_method[PETSC_MAX_PATH_LEN]; 313 PetscErrorCode (*fp_get_coarsedm_context)(DM,void**) = NULL; 314 315 PetscCall(PetscSNPrintf(dmcoarse_method,sizeof(dmcoarse_method),"PCTelescopeGetCoarseDMNullSpaceUserContext")); 316 PetscCall(PetscObjectQueryFunction((PetscObject)ctx->dm_coarse,dmcoarse_method,&fp_get_coarsedm_context)); 317 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); 318 PetscCall(MatNullSpaceSetFunction(sub_nullspace,nullspace->remove,context)); 319 } 320 } 321 322 if (PCTelescope_isActiveRank(sred)) { 323 PetscCall(MatSetNullSpace(sub_mat,sub_nullspace)); 324 PetscCall(MatNullSpaceDestroy(&sub_nullspace)); 325 } 326 } 327 } 328 { 329 MatNullSpace nearnullspace,sub_nearnullspace; 330 PetscCall(MatGetNearNullSpace(B,&nearnullspace)); 331 if (nearnullspace) { 332 PetscCall(PetscInfo(pc,"PCTelescope: generating near nullspace (CoarseDM)\n")); 333 PetscCall(PCTelescopeSubNullSpaceCreate_CoarseDM(pc,sred,nearnullspace,&sub_nearnullspace)); 334 335 /* attach any user nullspace removal methods and contexts */ 336 if (PCTelescope_isActiveRank(sred)) { 337 void *context = NULL; 338 if (nearnullspace->remove && !nearnullspace->rmctx) { 339 PetscCall(MatNullSpaceSetFunction(sub_nearnullspace,nearnullspace->remove,context)); 340 } else if (nearnullspace->remove && nearnullspace->rmctx) { 341 char dmcoarse_method[PETSC_MAX_PATH_LEN]; 342 PetscErrorCode (*fp_get_coarsedm_context)(DM,void**) = NULL; 343 344 PetscCall(PetscSNPrintf(dmcoarse_method,sizeof(dmcoarse_method),"PCTelescopeGetCoarseDMNearNullSpaceUserContext")); 345 PetscCall(PetscObjectQueryFunction((PetscObject)ctx->dm_coarse,dmcoarse_method,&fp_get_coarsedm_context)); 346 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); 347 PetscCall(MatNullSpaceSetFunction(sub_nearnullspace,nearnullspace->remove,context)); 348 } 349 } 350 351 if (PCTelescope_isActiveRank(sred)) { 352 PetscCall(MatSetNearNullSpace(sub_mat,sub_nearnullspace)); 353 PetscCall(MatNullSpaceDestroy(&sub_nearnullspace)); 354 } 355 } 356 } 357 PetscFunctionReturn(0); 358 } 359 360 PetscErrorCode PCReset_Telescope_CoarseDM(PC pc) 361 { 362 PC_Telescope sred = (PC_Telescope)pc->data; 363 PC_Telescope_CoarseDMCtx *ctx; 364 365 PetscFunctionBegin; 366 ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx; 367 ctx->dm_fine = NULL; /* since I did not increment the ref counter we set these to NULL */ 368 ctx->dm_coarse = NULL; /* since I did not increment the ref counter we set these to NULL */ 369 ctx->permutation = NULL; /* this will be fetched from the dm so no need to call destroy */ 370 PetscCall(VecDestroy(&ctx->xp)); 371 ctx->fp_dm_field_scatter = NULL; 372 ctx->fp_dm_state_scatter = NULL; 373 ctx->dmksp_context_determined = NULL; 374 ctx->dmksp_context_user = NULL; 375 PetscFunctionReturn(0); 376 } 377 378 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) 379 { 380 PC_Telescope sred = (PC_Telescope)pc->data; 381 Vec yred = NULL; 382 PetscBool default_init_guess_value = PETSC_FALSE; 383 PC_Telescope_CoarseDMCtx *ctx; 384 385 PetscFunctionBegin; 386 ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx; 387 yred = sred->yred; 388 389 PetscCheckFalse(its > 1,PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PCApplyRichardson_Telescope_CoarseDM only supports max_it = 1"); 390 *reason = (PCRichardsonConvergedReason)0; 391 392 if (!zeroguess) { 393 PetscCall(PetscInfo(pc,"PCTelescopeCoarseDM: Scattering y for non-zero-initial guess\n")); 394 395 PetscCall(ctx->fp_dm_field_scatter(ctx->dm_fine,y,SCATTER_FORWARD,ctx->dm_coarse,yred)); 396 } 397 398 if (PCTelescope_isActiveRank(sred)) { 399 PetscCall(KSPGetInitialGuessNonzero(sred->ksp,&default_init_guess_value)); 400 if (!zeroguess) PetscCall(KSPSetInitialGuessNonzero(sred->ksp,PETSC_TRUE)); 401 } 402 403 PetscCall(PCApply_Telescope_CoarseDM(pc,x,y)); 404 405 if (PCTelescope_isActiveRank(sred)) { 406 PetscCall(KSPSetInitialGuessNonzero(sred->ksp,default_init_guess_value)); 407 } 408 409 if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS; 410 *outits = 1; 411 PetscFunctionReturn(0); 412 } 413