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