1 2 #include <petsc/private/matimpl.h> 3 #include <petsc/private/pcimpl.h> 4 #include <petscksp.h> /*I "petscksp.h" I*/ 5 #include <petscdm.h> /*I "petscdm.h" I*/ 6 #include "../src/ksp/pc/impls/telescope/telescope.h" 7 8 static PetscBool cited = PETSC_FALSE; 9 static const char citation[] = 10 "@inproceedings{MaySananRuppKnepleySmith2016,\n" 11 " title = {Extreme-Scale Multigrid Components within PETSc},\n" 12 " author = {Dave A. May and Patrick Sanan and Karl Rupp and Matthew G. Knepley and Barry F. Smith},\n" 13 " booktitle = {Proceedings of the Platform for Advanced Scientific Computing Conference},\n" 14 " series = {PASC '16},\n" 15 " isbn = {978-1-4503-4126-4},\n" 16 " location = {Lausanne, Switzerland},\n" 17 " pages = {5:1--5:12},\n" 18 " articleno = {5},\n" 19 " numpages = {12},\n" 20 " url = {http://doi.acm.org/10.1145/2929908.2929913},\n" 21 " doi = {10.1145/2929908.2929913},\n" 22 " acmid = {2929913},\n" 23 " publisher = {ACM},\n" 24 " address = {New York, NY, USA},\n" 25 " keywords = {GPU, HPC, agglomeration, coarse-level solver, multigrid, parallel computing, preconditioning},\n" 26 " year = {2016}\n" 27 "}\n"; 28 29 /* 30 PCTelescopeSetUp_default() 31 PCTelescopeMatCreate_default() 32 33 default 34 35 // scatter in 36 x(comm) -> xtmp(comm) 37 38 xred(subcomm) <- xtmp 39 yred(subcomm) 40 41 yred(subcomm) --> xtmp 42 43 // scatter out 44 xtmp(comm) -> y(comm) 45 */ 46 47 PetscBool PetscSubComm_isActiveRank(PetscSubcomm scomm) 48 { 49 if (scomm->color == 0) { return PETSC_TRUE; } 50 else { return PETSC_FALSE; } 51 } 52 53 PetscBool isActiveRank(PC_Telescope sred) 54 { 55 if (sred->psubcomm) { 56 return(PetscSubComm_isActiveRank(sred->psubcomm)); 57 } else { 58 if (sred->subcomm != MPI_COMM_NULL) { return PETSC_TRUE; } 59 else { return PETSC_FALSE; } 60 } 61 } 62 63 DM private_PCTelescopeGetSubDM(PC_Telescope sred) 64 { 65 DM subdm = NULL; 66 67 if (!isActiveRank(sred)) { subdm = NULL; } 68 else { 69 switch (sred->sr_type) { 70 case TELESCOPE_DEFAULT: subdm = NULL; 71 break; 72 case TELESCOPE_DMDA: subdm = ((PC_Telescope_DMDACtx*)sred->dm_ctx)->dmrepart; 73 break; 74 case TELESCOPE_DMPLEX: subdm = NULL; 75 break; 76 } 77 } 78 return(subdm); 79 } 80 81 PetscErrorCode PCTelescopeSetUp_default(PC pc,PC_Telescope sred) 82 { 83 PetscErrorCode ierr; 84 PetscInt m,M,bs,st,ed; 85 Vec x,xred,yred,xtmp; 86 Mat B; 87 MPI_Comm comm,subcomm; 88 VecScatter scatter; 89 IS isin; 90 91 PetscFunctionBegin; 92 ierr = PetscInfo(pc,"PCTelescope: setup (default)\n");CHKERRQ(ierr); 93 comm = PetscSubcommParent(sred->psubcomm); 94 subcomm = PetscSubcommChild(sred->psubcomm); 95 96 ierr = PCGetOperators(pc,NULL,&B);CHKERRQ(ierr); 97 ierr = MatGetSize(B,&M,NULL);CHKERRQ(ierr); 98 ierr = MatGetBlockSize(B,&bs);CHKERRQ(ierr); 99 ierr = MatCreateVecs(B,&x,NULL);CHKERRQ(ierr); 100 101 xred = NULL; 102 m = 0; 103 if (isActiveRank(sred)) { 104 ierr = VecCreate(subcomm,&xred);CHKERRQ(ierr); 105 ierr = VecSetSizes(xred,PETSC_DECIDE,M);CHKERRQ(ierr); 106 ierr = VecSetBlockSize(xred,bs);CHKERRQ(ierr); 107 ierr = VecSetFromOptions(xred);CHKERRQ(ierr); 108 ierr = VecGetLocalSize(xred,&m);CHKERRQ(ierr); 109 } 110 111 yred = NULL; 112 if (isActiveRank(sred)) { 113 ierr = VecDuplicate(xred,&yred);CHKERRQ(ierr); 114 } 115 116 ierr = VecCreate(comm,&xtmp);CHKERRQ(ierr); 117 ierr = VecSetSizes(xtmp,m,PETSC_DECIDE);CHKERRQ(ierr); 118 ierr = VecSetBlockSize(xtmp,bs);CHKERRQ(ierr); 119 ierr = VecSetType(xtmp,((PetscObject)x)->type_name);CHKERRQ(ierr); 120 121 if (isActiveRank(sred)) { 122 ierr = VecGetOwnershipRange(xred,&st,&ed);CHKERRQ(ierr); 123 ierr = ISCreateStride(comm,(ed-st),st,1,&isin);CHKERRQ(ierr); 124 } else { 125 ierr = VecGetOwnershipRange(x,&st,&ed);CHKERRQ(ierr); 126 ierr = ISCreateStride(comm,0,st,1,&isin);CHKERRQ(ierr); 127 } 128 ierr = ISSetBlockSize(isin,bs);CHKERRQ(ierr); 129 130 ierr = VecScatterCreateWithData(x,isin,xtmp,NULL,&scatter);CHKERRQ(ierr); 131 132 sred->isin = isin; 133 sred->scatter = scatter; 134 sred->xred = xred; 135 sred->yred = yred; 136 sred->xtmp = xtmp; 137 ierr = VecDestroy(&x);CHKERRQ(ierr); 138 PetscFunctionReturn(0); 139 } 140 141 PetscErrorCode PCTelescopeMatCreate_default(PC pc,PC_Telescope sred,MatReuse reuse,Mat *A) 142 { 143 PetscErrorCode ierr; 144 MPI_Comm comm,subcomm; 145 Mat Bred,B; 146 PetscInt nr,nc; 147 IS isrow,iscol; 148 Mat Blocal,*_Blocal; 149 150 PetscFunctionBegin; 151 ierr = PetscInfo(pc,"PCTelescope: updating the redundant preconditioned operator (default)\n");CHKERRQ(ierr); 152 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 153 subcomm = PetscSubcommChild(sred->psubcomm); 154 ierr = PCGetOperators(pc,NULL,&B);CHKERRQ(ierr); 155 ierr = MatGetSize(B,&nr,&nc);CHKERRQ(ierr); 156 isrow = sred->isin; 157 ierr = ISCreateStride(comm,nc,0,1,&iscol);CHKERRQ(ierr); 158 ierr = MatCreateSubMatrices(B,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&_Blocal);CHKERRQ(ierr); 159 Blocal = *_Blocal; 160 ierr = PetscFree(_Blocal);CHKERRQ(ierr); 161 Bred = NULL; 162 if (isActiveRank(sred)) { 163 PetscInt mm; 164 165 if (reuse != MAT_INITIAL_MATRIX) { Bred = *A; } 166 167 ierr = MatGetSize(Blocal,&mm,NULL);CHKERRQ(ierr); 168 ierr = MatCreateMPIMatConcatenateSeqMat(subcomm,Blocal,mm,reuse,&Bred);CHKERRQ(ierr); 169 } 170 *A = Bred; 171 ierr = ISDestroy(&iscol);CHKERRQ(ierr); 172 ierr = MatDestroy(&Blocal);CHKERRQ(ierr); 173 PetscFunctionReturn(0); 174 } 175 176 static PetscErrorCode PCTelescopeSubNullSpaceCreate_Telescope(PC pc,PC_Telescope sred,MatNullSpace nullspace,MatNullSpace *sub_nullspace) 177 { 178 PetscErrorCode ierr; 179 PetscBool has_const; 180 const Vec *vecs; 181 Vec *sub_vecs = NULL; 182 PetscInt i,k,n = 0; 183 MPI_Comm subcomm; 184 185 PetscFunctionBegin; 186 subcomm = PetscSubcommChild(sred->psubcomm); 187 ierr = MatNullSpaceGetVecs(nullspace,&has_const,&n,&vecs);CHKERRQ(ierr); 188 189 if (isActiveRank(sred)) { 190 if (n) { 191 ierr = VecDuplicateVecs(sred->xred,n,&sub_vecs);CHKERRQ(ierr); 192 } 193 } 194 195 /* copy entries */ 196 for (k=0; k<n; k++) { 197 const PetscScalar *x_array; 198 PetscScalar *LA_sub_vec; 199 PetscInt st,ed; 200 201 /* pull in vector x->xtmp */ 202 ierr = VecScatterBegin(sred->scatter,vecs[k],sred->xtmp,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 203 ierr = VecScatterEnd(sred->scatter,vecs[k],sred->xtmp,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 204 if (sub_vecs) { 205 /* copy vector entries into xred */ 206 ierr = VecGetArrayRead(sred->xtmp,&x_array);CHKERRQ(ierr); 207 if (sub_vecs[k]) { 208 ierr = VecGetOwnershipRange(sub_vecs[k],&st,&ed);CHKERRQ(ierr); 209 ierr = VecGetArray(sub_vecs[k],&LA_sub_vec);CHKERRQ(ierr); 210 for (i=0; i<ed-st; i++) { 211 LA_sub_vec[i] = x_array[i]; 212 } 213 ierr = VecRestoreArray(sub_vecs[k],&LA_sub_vec);CHKERRQ(ierr); 214 } 215 ierr = VecRestoreArrayRead(sred->xtmp,&x_array);CHKERRQ(ierr); 216 } 217 } 218 219 if (isActiveRank(sred)) { 220 /* create new (near) nullspace for redundant object */ 221 ierr = MatNullSpaceCreate(subcomm,has_const,n,sub_vecs,sub_nullspace);CHKERRQ(ierr); 222 ierr = VecDestroyVecs(n,&sub_vecs);CHKERRQ(ierr); 223 if (nullspace->remove) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Propagation of custom remove callbacks not supported when propagating (near) nullspaces with PCTelescope"); 224 if (nullspace->rmctx) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Propagation of custom remove callback context not supported when propagating (near) nullspaces with PCTelescope"); 225 } 226 PetscFunctionReturn(0); 227 } 228 229 static PetscErrorCode PCTelescopeMatNullSpaceCreate_default(PC pc,PC_Telescope sred,Mat sub_mat) 230 { 231 PetscErrorCode ierr; 232 Mat B; 233 234 PetscFunctionBegin; 235 ierr = PCGetOperators(pc,NULL,&B);CHKERRQ(ierr); 236 237 /* Propagate the nullspace if it exists */ 238 { 239 MatNullSpace nullspace,sub_nullspace; 240 ierr = MatGetNullSpace(B,&nullspace);CHKERRQ(ierr); 241 if (nullspace) { 242 ierr = PetscInfo(pc,"PCTelescope: generating nullspace (default)\n");CHKERRQ(ierr); 243 ierr = PCTelescopeSubNullSpaceCreate_Telescope(pc,sred,nullspace,&sub_nullspace);CHKERRQ(ierr); 244 if (isActiveRank(sred)) { 245 ierr = MatSetNullSpace(sub_mat,sub_nullspace);CHKERRQ(ierr); 246 ierr = MatNullSpaceDestroy(&sub_nullspace);CHKERRQ(ierr); 247 } 248 } 249 } 250 251 /* Propagate the near nullspace if it exists */ 252 { 253 MatNullSpace nearnullspace,sub_nearnullspace; 254 ierr = MatGetNearNullSpace(B,&nearnullspace);CHKERRQ(ierr); 255 if (nearnullspace) { 256 ierr = PetscInfo(pc,"PCTelescope: generating near nullspace (default)\n");CHKERRQ(ierr); 257 ierr = PCTelescopeSubNullSpaceCreate_Telescope(pc,sred,nearnullspace,&sub_nearnullspace);CHKERRQ(ierr); 258 if (isActiveRank(sred)) { 259 ierr = MatSetNearNullSpace(sub_mat,sub_nearnullspace);CHKERRQ(ierr); 260 ierr = MatNullSpaceDestroy(&sub_nearnullspace);CHKERRQ(ierr); 261 } 262 } 263 } 264 PetscFunctionReturn(0); 265 } 266 267 static PetscErrorCode PCView_Telescope(PC pc,PetscViewer viewer) 268 { 269 PC_Telescope sred = (PC_Telescope)pc->data; 270 PetscErrorCode ierr; 271 PetscBool iascii,isstring; 272 PetscViewer subviewer; 273 274 PetscFunctionBegin; 275 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 276 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr); 277 if (iascii) { 278 if (!sred->psubcomm) { 279 ierr = PetscViewerASCIIPrintf(viewer," preconditioner not yet setup\n");CHKERRQ(ierr); 280 } else { 281 MPI_Comm comm,subcomm; 282 PetscMPIInt comm_size,subcomm_size; 283 DM dm,subdm; 284 285 ierr = PCGetDM(pc,&dm);CHKERRQ(ierr); 286 subdm = private_PCTelescopeGetSubDM(sred); 287 comm = PetscSubcommParent(sred->psubcomm); 288 subcomm = PetscSubcommChild(sred->psubcomm); 289 ierr = MPI_Comm_size(comm,&comm_size);CHKERRQ(ierr); 290 ierr = MPI_Comm_size(subcomm,&subcomm_size);CHKERRQ(ierr); 291 292 ierr = PetscViewerASCIIPrintf(viewer," parent comm size reduction factor = %D\n",sred->redfactor);CHKERRQ(ierr); 293 ierr = PetscViewerASCIIPrintf(viewer," comm_size = %d , subcomm_size = %d\n",(int)comm_size,(int)subcomm_size);CHKERRQ(ierr); 294 switch (sred->subcommtype) { 295 case PETSC_SUBCOMM_INTERLACED : 296 ierr = PetscViewerASCIIPrintf(viewer," subcomm type: interlaced\n",sred->subcommtype);CHKERRQ(ierr); 297 break; 298 case PETSC_SUBCOMM_CONTIGUOUS : 299 ierr = PetscViewerASCIIPrintf(viewer," subcomm type: contiguous\n",sred->subcommtype);CHKERRQ(ierr); 300 break; 301 default : 302 SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"General subcomm type not supported by PCTelescope"); 303 } 304 ierr = PetscViewerGetSubViewer(viewer,subcomm,&subviewer);CHKERRQ(ierr); 305 if (isActiveRank(sred)) { 306 ierr = PetscViewerASCIIPushTab(subviewer);CHKERRQ(ierr); 307 308 if (dm && sred->ignore_dm) { 309 ierr = PetscViewerASCIIPrintf(subviewer," ignoring DM\n");CHKERRQ(ierr); 310 } 311 if (sred->ignore_kspcomputeoperators) { 312 ierr = PetscViewerASCIIPrintf(subviewer," ignoring KSPComputeOperators\n");CHKERRQ(ierr); 313 } 314 switch (sred->sr_type) { 315 case TELESCOPE_DEFAULT: 316 ierr = PetscViewerASCIIPrintf(subviewer," using default setup\n");CHKERRQ(ierr); 317 break; 318 case TELESCOPE_DMDA: 319 ierr = PetscViewerASCIIPrintf(subviewer," DMDA detected\n");CHKERRQ(ierr); 320 ierr = DMView_DA_Short(subdm,subviewer);CHKERRQ(ierr); 321 break; 322 case TELESCOPE_DMPLEX: 323 ierr = PetscViewerASCIIPrintf(subviewer," DMPLEX detected\n");CHKERRQ(ierr); 324 break; 325 } 326 327 ierr = KSPView(sred->ksp,subviewer);CHKERRQ(ierr); 328 ierr = PetscViewerASCIIPopTab(subviewer);CHKERRQ(ierr); 329 } 330 ierr = PetscViewerRestoreSubViewer(viewer,subcomm,&subviewer);CHKERRQ(ierr); 331 } 332 } 333 PetscFunctionReturn(0); 334 } 335 336 static PetscErrorCode PCSetUp_Telescope(PC pc) 337 { 338 PC_Telescope sred = (PC_Telescope)pc->data; 339 PetscErrorCode ierr; 340 MPI_Comm comm,subcomm=0; 341 PCTelescopeType sr_type; 342 343 PetscFunctionBegin; 344 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 345 346 /* subcomm definition */ 347 if (!pc->setupcalled) { 348 if (!sred->psubcomm) { 349 ierr = PetscSubcommCreate(comm,&sred->psubcomm);CHKERRQ(ierr); 350 ierr = PetscSubcommSetNumber(sred->psubcomm,sred->redfactor);CHKERRQ(ierr); 351 ierr = PetscSubcommSetType(sred->psubcomm,sred->subcommtype);CHKERRQ(ierr); 352 ierr = PetscLogObjectMemory((PetscObject)pc,sizeof(PetscSubcomm));CHKERRQ(ierr); 353 sred->subcomm = PetscSubcommChild(sred->psubcomm); 354 } 355 } 356 subcomm = sred->subcomm; 357 358 /* internal KSP */ 359 if (!pc->setupcalled) { 360 const char *prefix; 361 362 if (isActiveRank(sred)) { 363 ierr = KSPCreate(subcomm,&sred->ksp);CHKERRQ(ierr); 364 ierr = KSPSetErrorIfNotConverged(sred->ksp,pc->erroriffailure);CHKERRQ(ierr); 365 ierr = PetscObjectIncrementTabLevel((PetscObject)sred->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 366 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)sred->ksp);CHKERRQ(ierr); 367 ierr = KSPSetType(sred->ksp,KSPPREONLY);CHKERRQ(ierr); 368 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr); 369 ierr = KSPSetOptionsPrefix(sred->ksp,prefix);CHKERRQ(ierr); 370 ierr = KSPAppendOptionsPrefix(sred->ksp,"telescope_");CHKERRQ(ierr); 371 } 372 } 373 /* Determine type of setup/update */ 374 if (!pc->setupcalled) { 375 PetscBool has_dm,same; 376 DM dm; 377 378 sr_type = TELESCOPE_DEFAULT; 379 has_dm = PETSC_FALSE; 380 ierr = PCGetDM(pc,&dm);CHKERRQ(ierr); 381 if (dm) { has_dm = PETSC_TRUE; } 382 if (has_dm) { 383 /* check for dmda */ 384 ierr = PetscObjectTypeCompare((PetscObject)dm,DMDA,&same);CHKERRQ(ierr); 385 if (same) { 386 ierr = PetscInfo(pc,"PCTelescope: found DMDA\n");CHKERRQ(ierr); 387 sr_type = TELESCOPE_DMDA; 388 } 389 /* check for dmplex */ 390 ierr = PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&same);CHKERRQ(ierr); 391 if (same) { 392 ierr = PetscInfo(pc,"PCTelescope: found DMPLEX\n");CHKERRQ(ierr); 393 sr_type = TELESCOPE_DMPLEX; 394 } 395 } 396 397 if (sred->ignore_dm) { 398 ierr = PetscInfo(pc,"PCTelescope: ignore DM\n");CHKERRQ(ierr); 399 sr_type = TELESCOPE_DEFAULT; 400 } 401 sred->sr_type = sr_type; 402 } else { 403 sr_type = sred->sr_type; 404 } 405 406 /* set function pointers for repartition setup, matrix creation/update, matrix (near) nullspace, and reset functionality */ 407 switch (sr_type) { 408 case TELESCOPE_DEFAULT: 409 sred->pctelescope_setup_type = PCTelescopeSetUp_default; 410 sred->pctelescope_matcreate_type = PCTelescopeMatCreate_default; 411 sred->pctelescope_matnullspacecreate_type = PCTelescopeMatNullSpaceCreate_default; 412 sred->pctelescope_reset_type = NULL; 413 break; 414 case TELESCOPE_DMDA: 415 pc->ops->apply = PCApply_Telescope_dmda; 416 pc->ops->applyrichardson = PCApplyRichardson_Telescope_dmda; 417 sred->pctelescope_setup_type = PCTelescopeSetUp_dmda; 418 sred->pctelescope_matcreate_type = PCTelescopeMatCreate_dmda; 419 sred->pctelescope_matnullspacecreate_type = PCTelescopeMatNullSpaceCreate_dmda; 420 sred->pctelescope_reset_type = PCReset_Telescope_dmda; 421 break; 422 case TELESCOPE_DMPLEX: SETERRQ(comm,PETSC_ERR_SUP,"Support for DMPLEX is currently not available"); 423 break; 424 default: SETERRQ(comm,PETSC_ERR_SUP,"Only support for repartitioning DMDA is provided"); 425 break; 426 } 427 428 /* setup */ 429 if (sred->pctelescope_setup_type) { 430 ierr = sred->pctelescope_setup_type(pc,sred);CHKERRQ(ierr); 431 } 432 /* update */ 433 if (!pc->setupcalled) { 434 if (sred->pctelescope_matcreate_type) { 435 ierr = sred->pctelescope_matcreate_type(pc,sred,MAT_INITIAL_MATRIX,&sred->Bred);CHKERRQ(ierr); 436 } 437 if (sred->pctelescope_matnullspacecreate_type) { 438 ierr = sred->pctelescope_matnullspacecreate_type(pc,sred,sred->Bred);CHKERRQ(ierr); 439 } 440 } else { 441 if (sred->pctelescope_matcreate_type) { 442 ierr = sred->pctelescope_matcreate_type(pc,sred,MAT_REUSE_MATRIX,&sred->Bred);CHKERRQ(ierr); 443 } 444 } 445 446 /* common - no construction */ 447 if (isActiveRank(sred)) { 448 ierr = KSPSetOperators(sred->ksp,sred->Bred,sred->Bred);CHKERRQ(ierr); 449 if (pc->setfromoptionscalled && !pc->setupcalled){ 450 ierr = KSPSetFromOptions(sred->ksp);CHKERRQ(ierr); 451 } 452 } 453 PetscFunctionReturn(0); 454 } 455 456 static PetscErrorCode PCApply_Telescope(PC pc,Vec x,Vec y) 457 { 458 PC_Telescope sred = (PC_Telescope)pc->data; 459 PetscErrorCode ierr; 460 Vec xtmp,xred,yred; 461 PetscInt i,st,ed; 462 VecScatter scatter; 463 PetscScalar *array; 464 const PetscScalar *x_array; 465 466 PetscFunctionBegin; 467 ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr); 468 469 xtmp = sred->xtmp; 470 scatter = sred->scatter; 471 xred = sred->xred; 472 yred = sred->yred; 473 474 /* pull in vector x->xtmp */ 475 ierr = VecScatterBegin(scatter,x,xtmp,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 476 ierr = VecScatterEnd(scatter,x,xtmp,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 477 478 /* copy vector entries into xred */ 479 ierr = VecGetArrayRead(xtmp,&x_array);CHKERRQ(ierr); 480 if (xred) { 481 PetscScalar *LA_xred; 482 ierr = VecGetOwnershipRange(xred,&st,&ed);CHKERRQ(ierr); 483 ierr = VecGetArray(xred,&LA_xred);CHKERRQ(ierr); 484 for (i=0; i<ed-st; i++) { 485 LA_xred[i] = x_array[i]; 486 } 487 ierr = VecRestoreArray(xred,&LA_xred);CHKERRQ(ierr); 488 } 489 ierr = VecRestoreArrayRead(xtmp,&x_array);CHKERRQ(ierr); 490 /* solve */ 491 if (isActiveRank(sred)) { 492 ierr = KSPSolve(sred->ksp,xred,yred);CHKERRQ(ierr); 493 ierr = KSPCheckSolve(sred->ksp,pc,yred);CHKERRQ(ierr); 494 } 495 /* return vector */ 496 ierr = VecGetArray(xtmp,&array);CHKERRQ(ierr); 497 if (yred) { 498 const PetscScalar *LA_yred; 499 ierr = VecGetOwnershipRange(yred,&st,&ed);CHKERRQ(ierr); 500 ierr = VecGetArrayRead(yred,&LA_yred);CHKERRQ(ierr); 501 for (i=0; i<ed-st; i++) { 502 array[i] = LA_yred[i]; 503 } 504 ierr = VecRestoreArrayRead(yred,&LA_yred);CHKERRQ(ierr); 505 } 506 ierr = VecRestoreArray(xtmp,&array);CHKERRQ(ierr); 507 ierr = VecScatterBegin(scatter,xtmp,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 508 ierr = VecScatterEnd(scatter,xtmp,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 509 PetscFunctionReturn(0); 510 } 511 512 static PetscErrorCode PCApplyRichardson_Telescope(PC pc,Vec x,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool zeroguess,PetscInt *outits,PCRichardsonConvergedReason *reason) 513 { 514 PC_Telescope sred = (PC_Telescope)pc->data; 515 PetscErrorCode ierr; 516 Vec xtmp,yred; 517 PetscInt i,st,ed; 518 VecScatter scatter; 519 const PetscScalar *x_array; 520 PetscBool default_init_guess_value; 521 522 PetscFunctionBegin; 523 xtmp = sred->xtmp; 524 scatter = sred->scatter; 525 yred = sred->yred; 526 527 if (its > 1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PCApplyRichardson_Telescope only supports max_it = 1"); 528 *reason = (PCRichardsonConvergedReason)0; 529 530 if (!zeroguess) { 531 ierr = PetscInfo(pc,"PCTelescope: Scattering y for non-zero initial guess\n");CHKERRQ(ierr); 532 /* pull in vector y->xtmp */ 533 ierr = VecScatterBegin(scatter,y,xtmp,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 534 ierr = VecScatterEnd(scatter,y,xtmp,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 535 536 /* copy vector entries into xred */ 537 ierr = VecGetArrayRead(xtmp,&x_array);CHKERRQ(ierr); 538 if (yred) { 539 PetscScalar *LA_yred; 540 ierr = VecGetOwnershipRange(yred,&st,&ed);CHKERRQ(ierr); 541 ierr = VecGetArray(yred,&LA_yred);CHKERRQ(ierr); 542 for (i=0; i<ed-st; i++) { 543 LA_yred[i] = x_array[i]; 544 } 545 ierr = VecRestoreArray(yred,&LA_yred);CHKERRQ(ierr); 546 } 547 ierr = VecRestoreArrayRead(xtmp,&x_array);CHKERRQ(ierr); 548 } 549 550 if (isActiveRank(sred)) { 551 ierr = KSPGetInitialGuessNonzero(sred->ksp,&default_init_guess_value);CHKERRQ(ierr); 552 if (!zeroguess) ierr = KSPSetInitialGuessNonzero(sred->ksp,PETSC_TRUE);CHKERRQ(ierr); 553 } 554 555 ierr = PCApply_Telescope(pc,x,y);CHKERRQ(ierr); 556 557 if (isActiveRank(sred)) { 558 ierr = KSPSetInitialGuessNonzero(sred->ksp,default_init_guess_value);CHKERRQ(ierr); 559 } 560 561 if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS; 562 *outits = 1; 563 PetscFunctionReturn(0); 564 } 565 566 static PetscErrorCode PCReset_Telescope(PC pc) 567 { 568 PC_Telescope sred = (PC_Telescope)pc->data; 569 PetscErrorCode ierr; 570 571 ierr = ISDestroy(&sred->isin);CHKERRQ(ierr); 572 ierr = VecScatterDestroy(&sred->scatter);CHKERRQ(ierr); 573 ierr = VecDestroy(&sred->xred);CHKERRQ(ierr); 574 ierr = VecDestroy(&sred->yred);CHKERRQ(ierr); 575 ierr = VecDestroy(&sred->xtmp);CHKERRQ(ierr); 576 ierr = MatDestroy(&sred->Bred);CHKERRQ(ierr); 577 ierr = KSPReset(sred->ksp);CHKERRQ(ierr); 578 if (sred->pctelescope_reset_type) { 579 ierr = sred->pctelescope_reset_type(pc);CHKERRQ(ierr); 580 } 581 PetscFunctionReturn(0); 582 } 583 584 static PetscErrorCode PCDestroy_Telescope(PC pc) 585 { 586 PC_Telescope sred = (PC_Telescope)pc->data; 587 PetscErrorCode ierr; 588 589 PetscFunctionBegin; 590 ierr = PCReset_Telescope(pc);CHKERRQ(ierr); 591 ierr = KSPDestroy(&sred->ksp);CHKERRQ(ierr); 592 ierr = PetscSubcommDestroy(&sred->psubcomm);CHKERRQ(ierr); 593 ierr = PetscFree(sred->dm_ctx);CHKERRQ(ierr); 594 ierr = PetscFree(pc->data);CHKERRQ(ierr); 595 PetscFunctionReturn(0); 596 } 597 598 static PetscErrorCode PCSetFromOptions_Telescope(PetscOptionItems *PetscOptionsObject,PC pc) 599 { 600 PC_Telescope sred = (PC_Telescope)pc->data; 601 PetscErrorCode ierr; 602 MPI_Comm comm; 603 PetscMPIInt size; 604 PetscBool flg; 605 PetscSubcommType subcommtype; 606 607 PetscFunctionBegin; 608 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 609 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 610 ierr = PetscOptionsHead(PetscOptionsObject,"Telescope options");CHKERRQ(ierr); 611 ierr = PetscOptionsEnum("-pc_telescope_subcomm_type","Subcomm type (interlaced or contiguous)","PCTelescopeSetSubcommType",PetscSubcommTypes,(PetscEnum)sred->subcommtype,(PetscEnum*)&subcommtype,&flg);CHKERRQ(ierr); 612 if (flg) { 613 ierr = PCTelescopeSetSubcommType(pc,subcommtype);CHKERRQ(ierr); 614 } 615 ierr = PetscOptionsInt("-pc_telescope_reduction_factor","Factor to reduce comm size by","PCTelescopeSetReductionFactor",sred->redfactor,&sred->redfactor,0);CHKERRQ(ierr); 616 if (sred->redfactor > size) SETERRQ(comm,PETSC_ERR_ARG_WRONG,"-pc_telescope_reduction_factor <= comm size"); 617 ierr = PetscOptionsBool("-pc_telescope_ignore_dm","Ignore any DM attached to the PC","PCTelescopeSetIgnoreDM",sred->ignore_dm,&sred->ignore_dm,0);CHKERRQ(ierr); 618 ierr = PetscOptionsBool("-pc_telescope_ignore_kspcomputeoperators","Ignore method used to compute A","PCTelescopeSetIgnoreKSPComputeOperators",sred->ignore_kspcomputeoperators,&sred->ignore_kspcomputeoperators,0);CHKERRQ(ierr); 619 ierr = PetscOptionsTail();CHKERRQ(ierr); 620 PetscFunctionReturn(0); 621 } 622 623 /* PC simplementation specific API's */ 624 625 static PetscErrorCode PCTelescopeGetKSP_Telescope(PC pc,KSP *ksp) 626 { 627 PC_Telescope red = (PC_Telescope)pc->data; 628 PetscFunctionBegin; 629 if (ksp) *ksp = red->ksp; 630 PetscFunctionReturn(0); 631 } 632 633 static PetscErrorCode PCTelescopeGetSubcommType_Telescope(PC pc,PetscSubcommType *subcommtype) 634 { 635 PC_Telescope red = (PC_Telescope)pc->data; 636 PetscFunctionBegin; 637 if (subcommtype) *subcommtype = red->subcommtype; 638 PetscFunctionReturn(0); 639 } 640 641 static PetscErrorCode PCTelescopeSetSubcommType_Telescope(PC pc,PetscSubcommType subcommtype) 642 { 643 PC_Telescope red = (PC_Telescope)pc->data; 644 645 PetscFunctionBegin; 646 if (pc->setupcalled) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"You cannot change the subcommunicator type for PCTelescope after it has been set up."); 647 red->subcommtype = subcommtype; 648 PetscFunctionReturn(0); 649 } 650 651 static PetscErrorCode PCTelescopeGetReductionFactor_Telescope(PC pc,PetscInt *fact) 652 { 653 PC_Telescope red = (PC_Telescope)pc->data; 654 PetscFunctionBegin; 655 if (fact) *fact = red->redfactor; 656 PetscFunctionReturn(0); 657 } 658 659 static PetscErrorCode PCTelescopeSetReductionFactor_Telescope(PC pc,PetscInt fact) 660 { 661 PC_Telescope red = (PC_Telescope)pc->data; 662 PetscMPIInt size; 663 PetscErrorCode ierr; 664 665 PetscFunctionBegin; 666 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);CHKERRQ(ierr); 667 if (fact <= 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Reduction factor of telescoping PC %D must be positive",fact); 668 if (fact > size) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Reduction factor of telescoping PC %D must be <= comm.size",fact); 669 red->redfactor = fact; 670 PetscFunctionReturn(0); 671 } 672 673 static PetscErrorCode PCTelescopeGetIgnoreDM_Telescope(PC pc,PetscBool *v) 674 { 675 PC_Telescope red = (PC_Telescope)pc->data; 676 PetscFunctionBegin; 677 if (v) *v = red->ignore_dm; 678 PetscFunctionReturn(0); 679 } 680 681 static PetscErrorCode PCTelescopeSetIgnoreDM_Telescope(PC pc,PetscBool v) 682 { 683 PC_Telescope red = (PC_Telescope)pc->data; 684 PetscFunctionBegin; 685 red->ignore_dm = v; 686 PetscFunctionReturn(0); 687 } 688 689 static PetscErrorCode PCTelescopeGetIgnoreKSPComputeOperators_Telescope(PC pc,PetscBool *v) 690 { 691 PC_Telescope red = (PC_Telescope)pc->data; 692 PetscFunctionBegin; 693 if (v) *v = red->ignore_kspcomputeoperators; 694 PetscFunctionReturn(0); 695 } 696 697 static PetscErrorCode PCTelescopeSetIgnoreKSPComputeOperators_Telescope(PC pc,PetscBool v) 698 { 699 PC_Telescope red = (PC_Telescope)pc->data; 700 PetscFunctionBegin; 701 red->ignore_kspcomputeoperators = v; 702 PetscFunctionReturn(0); 703 } 704 705 static PetscErrorCode PCTelescopeGetDM_Telescope(PC pc,DM *dm) 706 { 707 PC_Telescope red = (PC_Telescope)pc->data; 708 PetscFunctionBegin; 709 *dm = private_PCTelescopeGetSubDM(red); 710 PetscFunctionReturn(0); 711 } 712 713 /*@ 714 PCTelescopeGetKSP - Gets the KSP created by the telescoping PC. 715 716 Not Collective 717 718 Input Parameter: 719 . pc - the preconditioner context 720 721 Output Parameter: 722 . subksp - the KSP defined the smaller set of processes 723 724 Level: advanced 725 726 .keywords: PC, telescopting solve 727 @*/ 728 PetscErrorCode PCTelescopeGetKSP(PC pc,KSP *subksp) 729 { 730 PetscErrorCode ierr; 731 PetscFunctionBegin; 732 ierr = PetscUseMethod(pc,"PCTelescopeGetKSP_C",(PC,KSP*),(pc,subksp));CHKERRQ(ierr); 733 PetscFunctionReturn(0); 734 } 735 736 /*@ 737 PCTelescopeGetReductionFactor - Gets the factor by which the original number of processes has been reduced by. 738 739 Not Collective 740 741 Input Parameter: 742 . pc - the preconditioner context 743 744 Output Parameter: 745 . fact - the reduction factor 746 747 Level: advanced 748 749 .keywords: PC, telescoping solve 750 @*/ 751 PetscErrorCode PCTelescopeGetReductionFactor(PC pc,PetscInt *fact) 752 { 753 PetscErrorCode ierr; 754 PetscFunctionBegin; 755 ierr = PetscUseMethod(pc,"PCTelescopeGetReductionFactor_C",(PC,PetscInt*),(pc,fact));CHKERRQ(ierr); 756 PetscFunctionReturn(0); 757 } 758 759 /*@ 760 PCTelescopeSetReductionFactor - Sets the factor by which the original number of processes has been reduced by. 761 762 Not Collective 763 764 Input Parameter: 765 . pc - the preconditioner context 766 767 Output Parameter: 768 . fact - the reduction factor 769 770 Level: advanced 771 772 .keywords: PC, telescoping solve 773 @*/ 774 PetscErrorCode PCTelescopeSetReductionFactor(PC pc,PetscInt fact) 775 { 776 PetscErrorCode ierr; 777 PetscFunctionBegin; 778 ierr = PetscTryMethod(pc,"PCTelescopeSetReductionFactor_C",(PC,PetscInt),(pc,fact));CHKERRQ(ierr); 779 PetscFunctionReturn(0); 780 } 781 782 /*@ 783 PCTelescopeGetIgnoreDM - Get the flag indicating if any DM attached to the PC will be used. 784 785 Not Collective 786 787 Input Parameter: 788 . pc - the preconditioner context 789 790 Output Parameter: 791 . v - the flag 792 793 Level: advanced 794 795 .keywords: PC, telescoping solve 796 @*/ 797 PetscErrorCode PCTelescopeGetIgnoreDM(PC pc,PetscBool *v) 798 { 799 PetscErrorCode ierr; 800 PetscFunctionBegin; 801 ierr = PetscUseMethod(pc,"PCTelescopeGetIgnoreDM_C",(PC,PetscBool*),(pc,v));CHKERRQ(ierr); 802 PetscFunctionReturn(0); 803 } 804 805 /*@ 806 PCTelescopeSetIgnoreDM - Set a flag to ignore any DM attached to the PC. 807 808 Not Collective 809 810 Input Parameter: 811 . pc - the preconditioner context 812 813 Output Parameter: 814 . v - Use PETSC_TRUE to ignore any DM 815 816 Level: advanced 817 818 .keywords: PC, telescoping solve 819 @*/ 820 PetscErrorCode PCTelescopeSetIgnoreDM(PC pc,PetscBool v) 821 { 822 PetscErrorCode ierr; 823 PetscFunctionBegin; 824 ierr = PetscTryMethod(pc,"PCTelescopeSetIgnoreDM_C",(PC,PetscBool),(pc,v));CHKERRQ(ierr); 825 PetscFunctionReturn(0); 826 } 827 828 /*@ 829 PCTelescopeGetIgnoreKSPComputeOperators - Get the flag indicating if KSPComputeOperators will be used. 830 831 Not Collective 832 833 Input Parameter: 834 . pc - the preconditioner context 835 836 Output Parameter: 837 . v - the flag 838 839 Level: advanced 840 841 .keywords: PC, telescoping solve 842 @*/ 843 PetscErrorCode PCTelescopeGetIgnoreKSPComputeOperators(PC pc,PetscBool *v) 844 { 845 PetscErrorCode ierr; 846 PetscFunctionBegin; 847 ierr = PetscUseMethod(pc,"PCTelescopeGetIgnoreKSPComputeOperators_C",(PC,PetscBool*),(pc,v));CHKERRQ(ierr); 848 PetscFunctionReturn(0); 849 } 850 851 /*@ 852 PCTelescopeSetIgnoreKSPComputeOperators - Set a flag to ignore KSPComputeOperators. 853 854 Not Collective 855 856 Input Parameter: 857 . pc - the preconditioner context 858 859 Output Parameter: 860 . v - Use PETSC_TRUE to ignore the method (if defined) set via KSPSetComputeOperators on pc 861 862 Level: advanced 863 864 .keywords: PC, telescoping solve 865 @*/ 866 PetscErrorCode PCTelescopeSetIgnoreKSPComputeOperators(PC pc,PetscBool v) 867 { 868 PetscErrorCode ierr; 869 PetscFunctionBegin; 870 ierr = PetscTryMethod(pc,"PCTelescopeSetIgnoreKSPComputeOperators_C",(PC,PetscBool),(pc,v));CHKERRQ(ierr); 871 PetscFunctionReturn(0); 872 } 873 874 /*@ 875 PCTelescopeGetDM - Get the re-partitioned DM attached to the sub KSP. 876 877 Not Collective 878 879 Input Parameter: 880 . pc - the preconditioner context 881 882 Output Parameter: 883 . subdm - The re-partitioned DM 884 885 Level: advanced 886 887 .keywords: PC, telescoping solve 888 @*/ 889 PetscErrorCode PCTelescopeGetDM(PC pc,DM *subdm) 890 { 891 PetscErrorCode ierr; 892 PetscFunctionBegin; 893 ierr = PetscUseMethod(pc,"PCTelescopeGetDM_C",(PC,DM*),(pc,subdm));CHKERRQ(ierr); 894 PetscFunctionReturn(0); 895 } 896 897 /*@ 898 PCTelescopeSetSubcommType - set subcommunicator type (interlaced or contiguous) 899 900 Logically Collective 901 902 Input Parameter: 903 + pc - the preconditioner context 904 - subcommtype - the subcommunicator type (see PetscSubcommType) 905 906 Level: advanced 907 908 .keywords: PC, telescoping solve 909 910 .seealso: PetscSubcommType, PetscSubcomm, PCTELESCOPE 911 @*/ 912 PetscErrorCode PCTelescopeSetSubcommType(PC pc, PetscSubcommType subcommtype) 913 { 914 PetscErrorCode ierr; 915 PetscFunctionBegin; 916 ierr = PetscTryMethod(pc,"PCTelescopeSetSubcommType_C",(PC,PetscSubcommType),(pc,subcommtype));CHKERRQ(ierr); 917 PetscFunctionReturn(0); 918 } 919 920 /*@ 921 PCTelescopeGetSubcommType - Get the subcommunicator type (interlaced or contiguous) 922 923 Not Collective 924 925 Input Parameter: 926 . pc - the preconditioner context 927 928 Output Parameter: 929 . subcommtype - the subcommunicator type (see PetscSubcommType) 930 931 Level: advanced 932 933 .keywords: PC, telescoping solve 934 935 .seealso: PetscSubcomm, PetscSubcommType, PCTELESCOPE 936 @*/ 937 PetscErrorCode PCTelescopeGetSubcommType(PC pc, PetscSubcommType *subcommtype) 938 { 939 PetscErrorCode ierr; 940 PetscFunctionBegin; 941 ierr = PetscUseMethod(pc,"PCTelescopeGetSubcommType_C",(PC,PetscSubcommType*),(pc,subcommtype));CHKERRQ(ierr); 942 PetscFunctionReturn(0); 943 } 944 945 /* -------------------------------------------------------------------------------------*/ 946 /*MC 947 PCTELESCOPE - Runs a KSP solver on a sub-group of processors. MPI processes not in the sub-communicator are idle during the solve. 948 949 Options Database: 950 + -pc_telescope_reduction_factor <r> - factor to use communicator size by. e.g. with 64 MPI processes and r=4, the new sub-communicator will have 64/4 = 16 ranks. 951 - -pc_telescope_ignore_dm - flag to indicate whether an attached DM should be ignored 952 - -pc_telescope_subcomm_type <interlaced,contiguous> - how to define the reduced communicator. see PetscSubcomm for more. 953 954 Level: advanced 955 956 Notes: 957 The preconditioner is deemed telescopic as it only calls KSPSolve() on a single 958 sub-communicator, in contrast with PCREDUNDANT which calls KSPSolve() on N sub-communicators. 959 This means there will be MPI processes which will be idle during the application of this preconditioner. 960 961 The default KSP is PREONLY. If a DM is attached to the PC, it is re-partitioned on the sub-communicator. 962 Both the Bmat operator and the right hand side vector are permuted into the new DOF ordering defined by the re-partitioned DM. 963 Currently only support for re-partitioning a DMDA is provided. 964 Any nullspace attached to the original Bmat operator is extracted, re-partitioned and set on the repartitioned Bmat operator. 965 KSPSetComputeOperators() is not propagated to the sub KSP. 966 Currently there is no support for the flag -pc_use_amat 967 968 Assuming that the parent preconditioner (PC) is defined on a communicator c, this implementation 969 creates a child sub-communicator (c') containing fewer MPI processes than the original parent preconditioner (PC). 970 971 Developer Notes: 972 During PCSetup, the B operator is scattered onto c'. 973 Within PCApply, the RHS vector (x) is scattered into a redundant vector, xred (defined on c'). 974 Then, KSPSolve() is executed on the c' communicator. 975 976 The communicator used within the telescoping preconditioner is defined by a PetscSubcomm using the INTERLACED 977 creation routine by default (this can be changed with -pc_telescope_subcomm_type). We run the sub KSP on only the ranks within the communicator which have a color equal to zero. 978 979 The telescoping preconditioner is aware of nullspaces and near nullspaces which are attached to the B operator. 980 In the case where B has a (near) nullspace attached, the (near) nullspace vectors are extracted from B and mapped into 981 a new (near) nullspace, defined on the sub-communicator, which is attached to B' (the B operator which was scattered to c') 982 983 The telescoping preconditioner is aware of an attached DM. In the event that the DM is of type DMDA (2D or 3D - 984 1D support for 1D DMDAs is not provided), a new DMDA is created on c' (e.g. it is re-partitioned), and this new DM 985 is attached the sub KSPSolve(). The design of telescope is such that it should be possible to extend support 986 for re-partitioning other to DM's (e.g. DMPLEX). The user can supply a flag to ignore attached DMs. 987 988 By default, B' is defined by simply fusing rows from different MPI processes 989 990 When a DMDA is attached to the parent preconditioner, B' is defined by: (i) performing a symmetric permutation of B 991 into the ordering defined by the DMDA on c', (ii) extracting the local chunks via MatCreateSubMatrices(), (iii) fusing the 992 locally (sequential) matrices defined on the ranks common to c and c' into B' using MatCreateMPIMatConcatenateSeqMat() 993 994 Limitations/improvements include the following. 995 VecPlaceArray() could be used within PCApply() to improve efficiency and reduce memory usage. 996 997 The symmetric permutation used when a DMDA is encountered is performed via explicitly assmbleming a permutation matrix P, 998 and performing P^T.A.P. Possibly it might be more efficient to use MatPermute(). We opted to use P^T.A.P as it appears 999 VecPermute() does not supported for the use case required here. By computing P, one can permute both the operator and RHS in a 1000 consistent manner. 1001 1002 Mapping of vectors is performed in the following way. 1003 Suppose the parent comm size was 4, and we set a reduction factor of 2; this would give a comm size on c' of 2. 1004 Using the interlaced creation routine, the ranks in c with color = 0 will be rank 0 and 2. 1005 We perform the scatter to the sub-comm in the following way. 1006 [1] Given a vector x defined on comm c 1007 1008 rank(c) : _________ 0 ______ ________ 1 _______ ________ 2 _____________ ___________ 3 __________ 1009 x : [0, 1, 2, 3, 4, 5] [6, 7, 8, 9, 10, 11] [12, 13, 14, 15, 16, 17] [18, 19, 20, 21, 22, 23] 1010 1011 scatter to xtmp defined also on comm c so that we have the following values 1012 1013 rank(c) : ___________________ 0 ________________ _1_ ______________________ 2 _______________________ __3_ 1014 xtmp : [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] [ ] [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23] [ ] 1015 1016 The entries on rank 1 and 3 (ranks which do not have a color = 0 in c') have no values 1017 1018 1019 [2] Copy the value from rank 0, 2 (indices with respect to comm c) into the vector xred which is defined on communicator c'. 1020 Ranks 0 and 2 are the only ranks in the subcomm which have a color = 0. 1021 1022 rank(c') : ___________________ 0 _______________ ______________________ 1 _____________________ 1023 xred : [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23] 1024 1025 1026 Contributed by Dave May 1027 1028 Reference: 1029 Dave A. May, Patrick Sanan, Karl Rupp, Matthew G. Knepley, and Barry F. Smith, "Extreme-Scale Multigrid Components within PETSc". 2016. In Proceedings of the Platform for Advanced Scientific Computing Conference (PASC '16). DOI: 10.1145/2929908.2929913 1030 1031 .seealso: PCTelescopeGetKSP(), PCTelescopeGetDM(), PCTelescopeGetReductionFactor(), PCTelescopeSetReductionFactor(), PCTelescopeGetIgnoreDM(), PCTelescopeSetIgnoreDM(), PCREDUNDANT 1032 M*/ 1033 PETSC_EXTERN PetscErrorCode PCCreate_Telescope(PC pc) 1034 { 1035 PetscErrorCode ierr; 1036 struct _PC_Telescope *sred; 1037 1038 PetscFunctionBegin; 1039 ierr = PetscNewLog(pc,&sred);CHKERRQ(ierr); 1040 sred->psubcomm = NULL; 1041 sred->subcommtype = PETSC_SUBCOMM_INTERLACED; 1042 sred->subcomm = MPI_COMM_NULL; 1043 sred->redfactor = 1; 1044 sred->ignore_dm = PETSC_FALSE; 1045 sred->ignore_kspcomputeoperators = PETSC_FALSE; 1046 pc->data = (void*)sred; 1047 1048 pc->ops->apply = PCApply_Telescope; 1049 pc->ops->applytranspose = NULL; 1050 pc->ops->applyrichardson = PCApplyRichardson_Telescope; 1051 pc->ops->setup = PCSetUp_Telescope; 1052 pc->ops->destroy = PCDestroy_Telescope; 1053 pc->ops->reset = PCReset_Telescope; 1054 pc->ops->setfromoptions = PCSetFromOptions_Telescope; 1055 pc->ops->view = PCView_Telescope; 1056 1057 sred->pctelescope_setup_type = PCTelescopeSetUp_default; 1058 sred->pctelescope_matcreate_type = PCTelescopeMatCreate_default; 1059 sred->pctelescope_matnullspacecreate_type = PCTelescopeMatNullSpaceCreate_default; 1060 sred->pctelescope_reset_type = NULL; 1061 1062 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeGetKSP_C",PCTelescopeGetKSP_Telescope);CHKERRQ(ierr); 1063 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeGetSubcommType_C",PCTelescopeGetSubcommType_Telescope);CHKERRQ(ierr); 1064 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeSetSubcommType_C",PCTelescopeSetSubcommType_Telescope);CHKERRQ(ierr); 1065 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeGetReductionFactor_C",PCTelescopeGetReductionFactor_Telescope);CHKERRQ(ierr); 1066 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeSetReductionFactor_C",PCTelescopeSetReductionFactor_Telescope);CHKERRQ(ierr); 1067 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeGetIgnoreDM_C",PCTelescopeGetIgnoreDM_Telescope);CHKERRQ(ierr); 1068 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeSetIgnoreDM_C",PCTelescopeSetIgnoreDM_Telescope);CHKERRQ(ierr); 1069 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeGetIgnoreKSPComputeOperators_C",PCTelescopeGetIgnoreKSPComputeOperators_Telescope);CHKERRQ(ierr); 1070 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeSetIgnoreKSPComputeOperators_C",PCTelescopeSetIgnoreKSPComputeOperators_Telescope);CHKERRQ(ierr); 1071 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeGetDM_C",PCTelescopeGetDM_Telescope);CHKERRQ(ierr); 1072 PetscFunctionReturn(0); 1073 } 1074