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 /* 9 PCTelescopeSetUp_default() 10 PCTelescopeMatCreate_default() 11 12 default 13 14 // scatter in 15 x(comm) -> xtmp(comm) 16 17 xred(subcomm) <- xtmp 18 yred(subcomm) 19 20 yred(subcomm) --> xtmp 21 22 // scatter out 23 xtmp(comm) -> y(comm) 24 */ 25 26 PetscBool isActiveRank(PetscSubcomm scomm) 27 { 28 if (scomm->color == 0) { return PETSC_TRUE; } 29 else { return PETSC_FALSE; } 30 } 31 32 #undef __FUNCT__ 33 #define __FUNCT__ "private_PCTelescopeGetSubDM" 34 DM private_PCTelescopeGetSubDM(PC_Telescope sred) 35 { 36 DM subdm = NULL; 37 38 if (!isActiveRank(sred->psubcomm)) { subdm = NULL; } 39 else { 40 switch (sred->sr_type) { 41 case TELESCOPE_DEFAULT: subdm = NULL; 42 break; 43 case TELESCOPE_DMDA: subdm = ((PC_Telescope_DMDACtx*)sred->dm_ctx)->dmrepart; 44 break; 45 case TELESCOPE_DMPLEX: subdm = NULL; 46 break; 47 } 48 } 49 return(subdm); 50 } 51 52 #undef __FUNCT__ 53 #define __FUNCT__ "PCTelescopeSetUp_default" 54 PetscErrorCode PCTelescopeSetUp_default(PC pc,PC_Telescope sred) 55 { 56 PetscErrorCode ierr; 57 PetscInt m,M,bs,st,ed; 58 Vec x,xred,yred,xtmp; 59 Mat B; 60 MPI_Comm comm,subcomm; 61 VecScatter scatter; 62 IS isin; 63 64 PetscFunctionBegin; 65 ierr = PetscInfo(pc,"PCTelescope: setup (default)\n");CHKERRQ(ierr); 66 comm = PetscSubcommParent(sred->psubcomm); 67 subcomm = PetscSubcommChild(sred->psubcomm); 68 69 ierr = PCGetOperators(pc,NULL,&B);CHKERRQ(ierr); 70 ierr = MatGetSize(B,&M,NULL);CHKERRQ(ierr); 71 ierr = MatGetBlockSize(B,&bs);CHKERRQ(ierr); 72 ierr = MatCreateVecs(B,&x,NULL);CHKERRQ(ierr); 73 74 xred = NULL; 75 m = 0; 76 if (isActiveRank(sred->psubcomm)) { 77 ierr = VecCreate(subcomm,&xred);CHKERRQ(ierr); 78 ierr = VecSetSizes(xred,PETSC_DECIDE,M);CHKERRQ(ierr); 79 ierr = VecSetBlockSize(xred,bs);CHKERRQ(ierr); 80 ierr = VecSetFromOptions(xred);CHKERRQ(ierr); 81 ierr = VecGetLocalSize(xred,&m);CHKERRQ(ierr); 82 } 83 84 yred = NULL; 85 if (isActiveRank(sred->psubcomm)) { 86 ierr = VecDuplicate(xred,&yred);CHKERRQ(ierr); 87 } 88 89 ierr = VecCreate(comm,&xtmp);CHKERRQ(ierr); 90 ierr = VecSetSizes(xtmp,m,PETSC_DECIDE);CHKERRQ(ierr); 91 ierr = VecSetBlockSize(xtmp,bs);CHKERRQ(ierr); 92 ierr = VecSetType(xtmp,((PetscObject)x)->type_name);CHKERRQ(ierr); 93 94 if (isActiveRank(sred->psubcomm)) { 95 ierr = VecGetOwnershipRange(xred,&st,&ed);CHKERRQ(ierr); 96 ierr = ISCreateStride(comm,(ed-st),st,1,&isin);CHKERRQ(ierr); 97 } else { 98 ierr = VecGetOwnershipRange(x,&st,&ed);CHKERRQ(ierr); 99 ierr = ISCreateStride(comm,0,st,1,&isin);CHKERRQ(ierr); 100 } 101 ierr = ISSetBlockSize(isin,bs);CHKERRQ(ierr); 102 103 ierr = VecScatterCreate(x,isin,xtmp,NULL,&scatter);CHKERRQ(ierr); 104 105 sred->isin = isin; 106 sred->scatter = scatter; 107 sred->xred = xred; 108 sred->yred = yred; 109 sred->xtmp = xtmp; 110 ierr = VecDestroy(&x);CHKERRQ(ierr); 111 PetscFunctionReturn(0); 112 } 113 114 #undef __FUNCT__ 115 #define __FUNCT__ "PCTelescopeMatCreate_default" 116 PetscErrorCode PCTelescopeMatCreate_default(PC pc,PC_Telescope sred,MatReuse reuse,Mat *A) 117 { 118 PetscErrorCode ierr; 119 MPI_Comm comm,subcomm; 120 Mat Bred,B; 121 PetscInt nr,nc; 122 IS isrow,iscol; 123 Mat Blocal,*_Blocal; 124 125 PetscFunctionBegin; 126 ierr = PetscInfo(pc,"PCTelescope: updating the redundant preconditioned operator (default)\n");CHKERRQ(ierr); 127 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 128 subcomm = PetscSubcommChild(sred->psubcomm); 129 ierr = PCGetOperators(pc,NULL,&B);CHKERRQ(ierr); 130 ierr = MatGetSize(B,&nr,&nc);CHKERRQ(ierr); 131 isrow = sred->isin; 132 ierr = ISCreateStride(comm,nc,0,1,&iscol);CHKERRQ(ierr); 133 ierr = MatGetSubMatrices(B,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&_Blocal);CHKERRQ(ierr); 134 Blocal = *_Blocal; 135 ierr = PetscFree(_Blocal);CHKERRQ(ierr); 136 Bred = NULL; 137 if (isActiveRank(sred->psubcomm)) { 138 PetscInt mm; 139 140 if (reuse != MAT_INITIAL_MATRIX) { Bred = *A; } 141 142 ierr = MatGetSize(Blocal,&mm,NULL);CHKERRQ(ierr); 143 ierr = MatCreateMPIMatConcatenateSeqMat(subcomm,Blocal,mm,reuse,&Bred);CHKERRQ(ierr); 144 } 145 *A = Bred; 146 ierr = ISDestroy(&iscol);CHKERRQ(ierr); 147 ierr = MatDestroy(&Blocal);CHKERRQ(ierr); 148 PetscFunctionReturn(0); 149 } 150 151 #undef __FUNCT__ 152 #define __FUNCT__ "PCTelescopeMatNullSpaceCreate_default" 153 PetscErrorCode PCTelescopeMatNullSpaceCreate_default(PC pc,PC_Telescope sred,Mat sub_mat) 154 { 155 PetscErrorCode ierr; 156 MatNullSpace nullspace,sub_nullspace; 157 Mat A,B; 158 PetscBool has_const; 159 PetscInt i,k,n = 0; 160 const Vec *vecs; 161 Vec *sub_vecs = NULL; 162 MPI_Comm subcomm; 163 164 PetscFunctionBegin; 165 ierr = PCGetOperators(pc,&A,&B);CHKERRQ(ierr); 166 ierr = MatGetNullSpace(B,&nullspace);CHKERRQ(ierr); 167 if (!nullspace) PetscFunctionReturn(0); 168 169 ierr = PetscInfo(pc,"PCTelescope: generating nullspace (default)\n");CHKERRQ(ierr); 170 subcomm = PetscSubcommChild(sred->psubcomm); 171 ierr = MatNullSpaceGetVecs(nullspace,&has_const,&n,&vecs);CHKERRQ(ierr); 172 173 if (isActiveRank(sred->psubcomm)) { 174 if (n) { 175 ierr = VecDuplicateVecs(sred->xred,n,&sub_vecs);CHKERRQ(ierr); 176 } 177 } 178 179 /* copy entries */ 180 for (k=0; k<n; k++) { 181 const PetscScalar *x_array; 182 PetscScalar *LA_sub_vec; 183 PetscInt st,ed; 184 185 /* pull in vector x->xtmp */ 186 ierr = VecScatterBegin(sred->scatter,vecs[k],sred->xtmp,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 187 ierr = VecScatterEnd(sred->scatter,vecs[k],sred->xtmp,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 188 if (sub_vecs) { 189 /* copy vector entries into xred */ 190 ierr = VecGetArrayRead(sred->xtmp,&x_array);CHKERRQ(ierr); 191 if (sub_vecs[k]) { 192 ierr = VecGetOwnershipRange(sub_vecs[k],&st,&ed);CHKERRQ(ierr); 193 ierr = VecGetArray(sub_vecs[k],&LA_sub_vec);CHKERRQ(ierr); 194 for (i=0; i<ed-st; i++) { 195 LA_sub_vec[i] = x_array[i]; 196 } 197 ierr = VecRestoreArray(sub_vecs[k],&LA_sub_vec);CHKERRQ(ierr); 198 } 199 ierr = VecRestoreArrayRead(sred->xtmp,&x_array);CHKERRQ(ierr); 200 } 201 } 202 203 if (isActiveRank(sred->psubcomm)) { 204 /* create new nullspace for redundant object */ 205 ierr = MatNullSpaceCreate(subcomm,has_const,n,sub_vecs,&sub_nullspace);CHKERRQ(ierr); 206 sub_nullspace->remove = nullspace->remove; 207 sub_nullspace->rmctx = nullspace->rmctx; 208 209 /* attach redundant nullspace to Bred */ 210 ierr = MatSetNullSpace(sub_mat,sub_nullspace);CHKERRQ(ierr); 211 ierr = VecDestroyVecs(n,&sub_vecs);CHKERRQ(ierr); 212 } 213 PetscFunctionReturn(0); 214 } 215 216 #undef __FUNCT__ 217 #define __FUNCT__ "PCView_Telescope" 218 static PetscErrorCode PCView_Telescope(PC pc,PetscViewer viewer) 219 { 220 PC_Telescope sred = (PC_Telescope)pc->data; 221 PetscErrorCode ierr; 222 PetscBool iascii,isstring; 223 PetscViewer subviewer; 224 225 PetscFunctionBegin; 226 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 227 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr); 228 if (iascii) { 229 if (!sred->psubcomm) { 230 ierr = PetscViewerASCIIPrintf(viewer," Telescope: preconditioner not yet setup\n");CHKERRQ(ierr); 231 } else { 232 MPI_Comm comm,subcomm; 233 PetscMPIInt comm_size,subcomm_size; 234 DM dm,subdm; 235 236 ierr = PCGetDM(pc,&dm);CHKERRQ(ierr); 237 subdm = private_PCTelescopeGetSubDM(sred); 238 comm = PetscSubcommParent(sred->psubcomm); 239 subcomm = PetscSubcommChild(sred->psubcomm); 240 ierr = MPI_Comm_size(comm,&comm_size);CHKERRQ(ierr); 241 ierr = MPI_Comm_size(subcomm,&subcomm_size);CHKERRQ(ierr); 242 243 ierr = PetscViewerASCIIPrintf(viewer," Telescope: parent comm size reduction factor = %D\n",sred->redfactor);CHKERRQ(ierr); 244 ierr = PetscViewerASCIIPrintf(viewer," Telescope: comm_size = %d , subcomm_size = %d\n",(int)comm_size,(int)subcomm_size);CHKERRQ(ierr); 245 switch (sred->subcommtype) { 246 case PETSC_SUBCOMM_INTERLACED : 247 ierr = PetscViewerASCIIPrintf(viewer," Telescope: subcomm type: interlaced\n",sred->subcommtype);CHKERRQ(ierr); 248 break; 249 case PETSC_SUBCOMM_CONTIGUOUS : 250 ierr = PetscViewerASCIIPrintf(viewer," Telescope: subcomm type: contiguous\n",sred->subcommtype);CHKERRQ(ierr); 251 break; 252 default : 253 SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"General subcomm type not supported by PCTelescope"); 254 } 255 ierr = PetscViewerGetSubViewer(viewer,subcomm,&subviewer);CHKERRQ(ierr); 256 if (isActiveRank(sred->psubcomm)) { 257 ierr = PetscViewerASCIIPushTab(subviewer);CHKERRQ(ierr); 258 259 if (dm && sred->ignore_dm) { 260 ierr = PetscViewerASCIIPrintf(subviewer," Telescope: ignoring DM\n");CHKERRQ(ierr); 261 } 262 if (sred->ignore_kspcomputeoperators) { 263 ierr = PetscViewerASCIIPrintf(subviewer," Telescope: ignoring KSPComputeOperators\n");CHKERRQ(ierr); 264 } 265 switch (sred->sr_type) { 266 case TELESCOPE_DEFAULT: 267 ierr = PetscViewerASCIIPrintf(subviewer," Telescope: using default setup\n");CHKERRQ(ierr); 268 break; 269 case TELESCOPE_DMDA: 270 ierr = PetscViewerASCIIPrintf(subviewer," Telescope: DMDA detected\n");CHKERRQ(ierr); 271 ierr = DMView_DMDAShort(subdm,subviewer);CHKERRQ(ierr); 272 break; 273 case TELESCOPE_DMPLEX: 274 ierr = PetscViewerASCIIPrintf(subviewer," Telescope: DMPLEX detected\n");CHKERRQ(ierr); 275 break; 276 } 277 278 ierr = KSPView(sred->ksp,subviewer);CHKERRQ(ierr); 279 ierr = PetscViewerASCIIPopTab(subviewer);CHKERRQ(ierr); 280 } 281 ierr = PetscViewerRestoreSubViewer(viewer,subcomm,&subviewer);CHKERRQ(ierr); 282 } 283 } 284 PetscFunctionReturn(0); 285 } 286 287 #undef __FUNCT__ 288 #define __FUNCT__ "PCSetUp_Telescope" 289 static PetscErrorCode PCSetUp_Telescope(PC pc) 290 { 291 PC_Telescope sred = (PC_Telescope)pc->data; 292 PetscErrorCode ierr; 293 MPI_Comm comm,subcomm=0; 294 PCTelescopeType sr_type; 295 296 PetscFunctionBegin; 297 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 298 299 /* subcomm definition */ 300 if (!pc->setupcalled) { 301 if (!sred->psubcomm) { 302 ierr = PetscSubcommCreate(comm,&sred->psubcomm);CHKERRQ(ierr); 303 ierr = PetscSubcommSetNumber(sred->psubcomm,sred->redfactor);CHKERRQ(ierr); 304 ierr = PetscSubcommSetType(sred->psubcomm,sred->subcommtype);CHKERRQ(ierr); 305 ierr = PetscLogObjectMemory((PetscObject)pc,sizeof(PetscSubcomm));CHKERRQ(ierr); 306 } 307 } 308 subcomm = PetscSubcommChild(sred->psubcomm); 309 310 /* internal KSP */ 311 if (!pc->setupcalled) { 312 const char *prefix; 313 314 if (isActiveRank(sred->psubcomm)) { 315 ierr = KSPCreate(subcomm,&sred->ksp);CHKERRQ(ierr); 316 ierr = KSPSetErrorIfNotConverged(sred->ksp,pc->erroriffailure);CHKERRQ(ierr); 317 ierr = PetscObjectIncrementTabLevel((PetscObject)sred->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 318 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)sred->ksp);CHKERRQ(ierr); 319 ierr = KSPSetType(sred->ksp,KSPPREONLY);CHKERRQ(ierr); 320 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr); 321 ierr = KSPSetOptionsPrefix(sred->ksp,prefix);CHKERRQ(ierr); 322 ierr = KSPAppendOptionsPrefix(sred->ksp,"telescope_");CHKERRQ(ierr); 323 } 324 } 325 /* Determine type of setup/update */ 326 if (!pc->setupcalled) { 327 PetscBool has_dm,same; 328 DM dm; 329 330 sr_type = TELESCOPE_DEFAULT; 331 has_dm = PETSC_FALSE; 332 ierr = PCGetDM(pc,&dm);CHKERRQ(ierr); 333 if (dm) { has_dm = PETSC_TRUE; } 334 if (has_dm) { 335 /* check for dmda */ 336 ierr = PetscObjectTypeCompare((PetscObject)dm,DMDA,&same);CHKERRQ(ierr); 337 if (same) { 338 ierr = PetscInfo(pc,"PCTelescope: found DMDA\n");CHKERRQ(ierr); 339 sr_type = TELESCOPE_DMDA; 340 } 341 /* check for dmplex */ 342 ierr = PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&same);CHKERRQ(ierr); 343 if (same) { 344 PetscInfo(pc,"PCTelescope: found DMPLEX\n"); 345 sr_type = TELESCOPE_DMPLEX; 346 } 347 } 348 349 if (sred->ignore_dm) { 350 PetscInfo(pc,"PCTelescope: ignore DM\n"); 351 sr_type = TELESCOPE_DEFAULT; 352 } 353 sred->sr_type = sr_type; 354 } else { 355 sr_type = sred->sr_type; 356 } 357 358 /* set function pointers for repartition setup, matrix creation/update, matrix nullspace and reset functionality */ 359 switch (sr_type) { 360 case TELESCOPE_DEFAULT: 361 sred->pctelescope_setup_type = PCTelescopeSetUp_default; 362 sred->pctelescope_matcreate_type = PCTelescopeMatCreate_default; 363 sred->pctelescope_matnullspacecreate_type = PCTelescopeMatNullSpaceCreate_default; 364 sred->pctelescope_reset_type = NULL; 365 break; 366 case TELESCOPE_DMDA: 367 pc->ops->apply = PCApply_Telescope_dmda; 368 pc->ops->applyrichardson = PCApplyRichardson_Telescope_dmda; 369 sred->pctelescope_setup_type = PCTelescopeSetUp_dmda; 370 sred->pctelescope_matcreate_type = PCTelescopeMatCreate_dmda; 371 sred->pctelescope_matnullspacecreate_type = PCTelescopeMatNullSpaceCreate_dmda; 372 sred->pctelescope_reset_type = PCReset_Telescope_dmda; 373 break; 374 case TELESCOPE_DMPLEX: SETERRQ(comm,PETSC_ERR_SUP,"Supprt for DMPLEX is currently not available"); 375 break; 376 default: SETERRQ(comm,PETSC_ERR_SUP,"Only supprt for repartitioning DMDA is provided"); 377 default: SETERRQ(comm,PETSC_ERR_SUP,"Only support for repartitioning DMDA is provided"); 378 break; 379 } 380 381 /* setup */ 382 if (sred->pctelescope_setup_type) { 383 ierr = sred->pctelescope_setup_type(pc,sred);CHKERRQ(ierr); 384 } 385 /* update */ 386 if (!pc->setupcalled) { 387 if (sred->pctelescope_matcreate_type) { 388 ierr = sred->pctelescope_matcreate_type(pc,sred,MAT_INITIAL_MATRIX,&sred->Bred);CHKERRQ(ierr); 389 } 390 if (sred->pctelescope_matnullspacecreate_type) { 391 ierr = sred->pctelescope_matnullspacecreate_type(pc,sred,sred->Bred);CHKERRQ(ierr); 392 } 393 } else { 394 if (sred->pctelescope_matcreate_type) { 395 ierr = sred->pctelescope_matcreate_type(pc,sred,MAT_REUSE_MATRIX,&sred->Bred);CHKERRQ(ierr); 396 } 397 } 398 399 /* common - no construction */ 400 if (isActiveRank(sred->psubcomm)) { 401 ierr = KSPSetOperators(sred->ksp,sred->Bred,sred->Bred);CHKERRQ(ierr); 402 if (pc->setfromoptionscalled && !pc->setupcalled){ 403 ierr = KSPSetFromOptions(sred->ksp);CHKERRQ(ierr); 404 } 405 } 406 PetscFunctionReturn(0); 407 } 408 409 #undef __FUNCT__ 410 #define __FUNCT__ "PCApply_Telescope" 411 static PetscErrorCode PCApply_Telescope(PC pc,Vec x,Vec y) 412 { 413 PC_Telescope sred = (PC_Telescope)pc->data; 414 PetscErrorCode ierr; 415 Vec xtmp,xred,yred; 416 PetscInt i,st,ed; 417 VecScatter scatter; 418 PetscScalar *array; 419 const PetscScalar *x_array; 420 421 PetscFunctionBegin; 422 xtmp = sred->xtmp; 423 scatter = sred->scatter; 424 xred = sred->xred; 425 yred = sred->yred; 426 427 /* pull in vector x->xtmp */ 428 ierr = VecScatterBegin(scatter,x,xtmp,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 429 ierr = VecScatterEnd(scatter,x,xtmp,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 430 431 /* copy vector entires into xred */ 432 ierr = VecGetArrayRead(xtmp,&x_array);CHKERRQ(ierr); 433 if (xred) { 434 PetscScalar *LA_xred; 435 ierr = VecGetOwnershipRange(xred,&st,&ed);CHKERRQ(ierr); 436 ierr = VecGetArray(xred,&LA_xred);CHKERRQ(ierr); 437 for (i=0; i<ed-st; i++) { 438 LA_xred[i] = x_array[i]; 439 } 440 ierr = VecRestoreArray(xred,&LA_xred);CHKERRQ(ierr); 441 } 442 ierr = VecRestoreArrayRead(xtmp,&x_array);CHKERRQ(ierr); 443 /* solve */ 444 if (isActiveRank(sred->psubcomm)) { 445 ierr = KSPSolve(sred->ksp,xred,yred);CHKERRQ(ierr); 446 } 447 /* return vector */ 448 ierr = VecGetArray(xtmp,&array);CHKERRQ(ierr); 449 if (yred) { 450 const PetscScalar *LA_yred; 451 ierr = VecGetOwnershipRange(yred,&st,&ed);CHKERRQ(ierr); 452 ierr = VecGetArrayRead(yred,&LA_yred);CHKERRQ(ierr); 453 for (i=0; i<ed-st; i++) { 454 array[i] = LA_yred[i]; 455 } 456 ierr = VecRestoreArrayRead(yred,&LA_yred);CHKERRQ(ierr); 457 } 458 ierr = VecRestoreArray(xtmp,&array);CHKERRQ(ierr); 459 ierr = VecScatterBegin(scatter,xtmp,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 460 ierr = VecScatterEnd(scatter,xtmp,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 461 PetscFunctionReturn(0); 462 } 463 464 #undef __FUNCT__ 465 #define __FUNCT__ "PCApplyRichardson_Telescope" 466 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) 467 { 468 PC_Telescope sred = (PC_Telescope)pc->data; 469 PetscErrorCode ierr; 470 Vec xtmp,yred; 471 PetscInt i,st,ed; 472 VecScatter scatter; 473 const PetscScalar *x_array; 474 PetscBool default_init_guess_value; 475 476 PetscFunctionBegin; 477 xtmp = sred->xtmp; 478 scatter = sred->scatter; 479 yred = sred->yred; 480 481 if (its > 1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PCApplyRichardson_Telescope only supports max_it = 1"); 482 *reason = (PCRichardsonConvergedReason)0; 483 484 if (!zeroguess) { 485 ierr = PetscInfo(pc,"PCTelescope: Scattering y for non-zero initial guess\n");CHKERRQ(ierr); 486 /* pull in vector y->xtmp */ 487 ierr = VecScatterBegin(scatter,y,xtmp,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 488 ierr = VecScatterEnd(scatter,y,xtmp,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 489 490 /* copy vector entires into xred */ 491 ierr = VecGetArrayRead(xtmp,&x_array);CHKERRQ(ierr); 492 if (yred) { 493 PetscScalar *LA_yred; 494 ierr = VecGetOwnershipRange(yred,&st,&ed);CHKERRQ(ierr); 495 ierr = VecGetArray(yred,&LA_yred);CHKERRQ(ierr); 496 for (i=0; i<ed-st; i++) { 497 LA_yred[i] = x_array[i]; 498 } 499 ierr = VecRestoreArray(yred,&LA_yred);CHKERRQ(ierr); 500 } 501 ierr = VecRestoreArrayRead(xtmp,&x_array);CHKERRQ(ierr); 502 } 503 504 if (isActiveRank(sred->psubcomm)) { 505 ierr = KSPGetInitialGuessNonzero(sred->ksp,&default_init_guess_value);CHKERRQ(ierr); 506 if (!zeroguess) ierr = KSPSetInitialGuessNonzero(sred->ksp,PETSC_TRUE);CHKERRQ(ierr); 507 } 508 509 ierr = PCApply_Telescope(pc,x,y);CHKERRQ(ierr); 510 511 if (isActiveRank(sred->psubcomm)) { 512 ierr = KSPSetInitialGuessNonzero(sred->ksp,default_init_guess_value);CHKERRQ(ierr); 513 } 514 515 if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS; 516 *outits = 1; 517 PetscFunctionReturn(0); 518 } 519 520 #undef __FUNCT__ 521 #define __FUNCT__ "PCReset_Telescope" 522 static PetscErrorCode PCReset_Telescope(PC pc) 523 { 524 PC_Telescope sred = (PC_Telescope)pc->data; 525 PetscErrorCode ierr; 526 527 ierr = ISDestroy(&sred->isin);CHKERRQ(ierr); 528 ierr = VecScatterDestroy(&sred->scatter);CHKERRQ(ierr); 529 ierr = VecDestroy(&sred->xred);CHKERRQ(ierr); 530 ierr = VecDestroy(&sred->yred);CHKERRQ(ierr); 531 ierr = VecDestroy(&sred->xtmp);CHKERRQ(ierr); 532 ierr = MatDestroy(&sred->Bred);CHKERRQ(ierr); 533 ierr = KSPReset(sred->ksp);CHKERRQ(ierr); 534 if (sred->pctelescope_reset_type) { 535 ierr = sred->pctelescope_reset_type(pc);CHKERRQ(ierr); 536 } 537 PetscFunctionReturn(0); 538 } 539 540 #undef __FUNCT__ 541 #define __FUNCT__ "PCDestroy_Telescope" 542 static PetscErrorCode PCDestroy_Telescope(PC pc) 543 { 544 PC_Telescope sred = (PC_Telescope)pc->data; 545 PetscErrorCode ierr; 546 547 PetscFunctionBegin; 548 ierr = PCReset_Telescope(pc);CHKERRQ(ierr); 549 ierr = KSPDestroy(&sred->ksp);CHKERRQ(ierr); 550 ierr = PetscSubcommDestroy(&sred->psubcomm);CHKERRQ(ierr); 551 ierr = PetscFree(sred->dm_ctx);CHKERRQ(ierr); 552 ierr = PetscFree(pc->data);CHKERRQ(ierr); 553 PetscFunctionReturn(0); 554 } 555 556 #undef __FUNCT__ 557 #define __FUNCT__ "PCSetFromOptions_Telescope" 558 static PetscErrorCode PCSetFromOptions_Telescope(PetscOptionItems *PetscOptionsObject,PC pc) 559 { 560 PC_Telescope sred = (PC_Telescope)pc->data; 561 PetscErrorCode ierr; 562 MPI_Comm comm; 563 PetscMPIInt size; 564 PetscBool flg; 565 PetscSubcommType subcommtype; 566 567 PetscFunctionBegin; 568 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 569 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 570 ierr = PetscOptionsHead(PetscOptionsObject,"Telescope options");CHKERRQ(ierr); 571 ierr = PetscOptionsEnum("-pc_telescope_subcomm_type","Subcomm type (interlaced or contiguous)","PCTelescopeSetSubcommType",PetscSubcommTypes,(PetscEnum)sred->subcommtype,(PetscEnum*)&subcommtype,&flg);CHKERRQ(ierr); 572 if (flg) { 573 ierr = PCTelescopeSetSubcommType(pc,subcommtype);CHKERRQ(ierr); 574 } 575 ierr = PetscOptionsInt("-pc_telescope_reduction_factor","Factor to reduce comm size by","PCTelescopeSetReductionFactor",sred->redfactor,&sred->redfactor,0);CHKERRQ(ierr); 576 if (sred->redfactor > size) SETERRQ(comm,PETSC_ERR_ARG_WRONG,"-pc_telescope_reduction_factor <= comm size"); 577 ierr = PetscOptionsBool("-pc_telescope_ignore_dm","Ignore any DM attached to the PC","PCTelescopeSetIgnoreDM",sred->ignore_dm,&sred->ignore_dm,0);CHKERRQ(ierr); 578 ierr = PetscOptionsBool("-pc_telescope_ignore_kspcomputeoperators","Ignore method used to compute A","PCTelescopeSetIgnoreKSPComputeOperators",sred->ignore_kspcomputeoperators,&sred->ignore_kspcomputeoperators,0);CHKERRQ(ierr); 579 ierr = PetscOptionsTail();CHKERRQ(ierr); 580 PetscFunctionReturn(0); 581 } 582 583 /* PC simplementation specific API's */ 584 585 #undef __FUNCT__ 586 #define __FUNCT__ "PCTelescopeGetKSP_Telescope" 587 static PetscErrorCode PCTelescopeGetKSP_Telescope(PC pc,KSP *ksp) 588 { 589 PC_Telescope red = (PC_Telescope)pc->data; 590 PetscFunctionBegin; 591 if (ksp) *ksp = red->ksp; 592 PetscFunctionReturn(0); 593 } 594 595 #undef __FUNCT__ 596 #define __FUNCT__ "PCTelescopeGetSubcommType_Telescope" 597 static PetscErrorCode PCTelescopeGetSubcommType_Telescope(PC pc,PetscSubcommType *subcommtype) 598 { 599 PC_Telescope red = (PC_Telescope)pc->data; 600 PetscFunctionBegin; 601 if (subcommtype) *subcommtype = red->subcommtype; 602 PetscFunctionReturn(0); 603 } 604 605 #undef __FUNCT__ 606 #define __FUNCT__ "PCTelescopeSetSubcommType_Telescope" 607 static PetscErrorCode PCTelescopeSetSubcommType_Telescope(PC pc,PetscSubcommType subcommtype) 608 { 609 PC_Telescope red = (PC_Telescope)pc->data; 610 611 PetscFunctionBegin; 612 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."); 613 red->subcommtype = subcommtype; 614 PetscFunctionReturn(0); 615 } 616 617 #undef __FUNCT__ 618 #define __FUNCT__ "PCTelescopeGetReductionFactor_Telescope" 619 static PetscErrorCode PCTelescopeGetReductionFactor_Telescope(PC pc,PetscInt *fact) 620 { 621 PC_Telescope red = (PC_Telescope)pc->data; 622 PetscFunctionBegin; 623 if (fact) *fact = red->redfactor; 624 PetscFunctionReturn(0); 625 } 626 627 #undef __FUNCT__ 628 #define __FUNCT__ "PCTelescopeSetReductionFactor_Telescope" 629 static PetscErrorCode PCTelescopeSetReductionFactor_Telescope(PC pc,PetscInt fact) 630 { 631 PC_Telescope red = (PC_Telescope)pc->data; 632 PetscMPIInt size; 633 PetscErrorCode ierr; 634 635 PetscFunctionBegin; 636 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);CHKERRQ(ierr); 637 if (fact <= 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Reduction factor of telescoping PC %D must be positive",fact); 638 if (fact > size) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Reduction factor of telescoping PC %D must be <= comm.size",fact); 639 red->redfactor = fact; 640 PetscFunctionReturn(0); 641 } 642 643 #undef __FUNCT__ 644 #define __FUNCT__ "PCTelescopeGetIgnoreDM_Telescope" 645 static PetscErrorCode PCTelescopeGetIgnoreDM_Telescope(PC pc,PetscBool *v) 646 { 647 PC_Telescope red = (PC_Telescope)pc->data; 648 PetscFunctionBegin; 649 if (v) *v = red->ignore_dm; 650 PetscFunctionReturn(0); 651 } 652 653 #undef __FUNCT__ 654 #define __FUNCT__ "PCTelescopeSetIgnoreDM_Telescope" 655 static PetscErrorCode PCTelescopeSetIgnoreDM_Telescope(PC pc,PetscBool v) 656 { 657 PC_Telescope red = (PC_Telescope)pc->data; 658 PetscFunctionBegin; 659 red->ignore_dm = v; 660 PetscFunctionReturn(0); 661 } 662 663 #undef __FUNCT__ 664 #define __FUNCT__ "PCTelescopeGetIgnoreKSPComputeOperators_Telescope" 665 static PetscErrorCode PCTelescopeGetIgnoreKSPComputeOperators_Telescope(PC pc,PetscBool *v) 666 { 667 PC_Telescope red = (PC_Telescope)pc->data; 668 PetscFunctionBegin; 669 if (v) *v = red->ignore_kspcomputeoperators; 670 PetscFunctionReturn(0); 671 } 672 673 #undef __FUNCT__ 674 #define __FUNCT__ "PCTelescopeSetIgnoreKSPComputeOperators_Telescope" 675 static PetscErrorCode PCTelescopeSetIgnoreKSPComputeOperators_Telescope(PC pc,PetscBool v) 676 { 677 PC_Telescope red = (PC_Telescope)pc->data; 678 PetscFunctionBegin; 679 red->ignore_kspcomputeoperators = v; 680 PetscFunctionReturn(0); 681 } 682 683 #undef __FUNCT__ 684 #define __FUNCT__ "PCTelescopeGetDM_Telescope" 685 static PetscErrorCode PCTelescopeGetDM_Telescope(PC pc,DM *dm) 686 { 687 PC_Telescope red = (PC_Telescope)pc->data; 688 PetscFunctionBegin; 689 *dm = private_PCTelescopeGetSubDM(red); 690 PetscFunctionReturn(0); 691 } 692 693 #undef __FUNCT__ 694 #define __FUNCT__ "PCTelescopeGetKSP" 695 /*@ 696 PCTelescopeGetKSP - Gets the KSP created by the telescoping PC. 697 698 Not Collective 699 700 Input Parameter: 701 . pc - the preconditioner context 702 703 Output Parameter: 704 . subksp - the KSP defined the smaller set of processes 705 706 Level: advanced 707 708 .keywords: PC, telescopting solve 709 @*/ 710 PetscErrorCode PCTelescopeGetKSP(PC pc,KSP *subksp) 711 { 712 PetscErrorCode ierr; 713 PetscFunctionBegin; 714 ierr = PetscUseMethod(pc,"PCTelescopeGetKSP_C",(PC,KSP*),(pc,subksp));CHKERRQ(ierr); 715 PetscFunctionReturn(0); 716 } 717 718 #undef __FUNCT__ 719 #define __FUNCT__ "PCTelescopeGetReductionFactor" 720 /*@ 721 PCTelescopeGetReductionFactor - Gets the factor by which the original number of processes has been reduced by. 722 723 Not Collective 724 725 Input Parameter: 726 . pc - the preconditioner context 727 728 Output Parameter: 729 . fact - the reduction factor 730 731 Level: advanced 732 733 .keywords: PC, telescoping solve 734 @*/ 735 PetscErrorCode PCTelescopeGetReductionFactor(PC pc,PetscInt *fact) 736 { 737 PetscErrorCode ierr; 738 PetscFunctionBegin; 739 ierr = PetscUseMethod(pc,"PCTelescopeGetReductionFactor_C",(PC,PetscInt*),(pc,fact));CHKERRQ(ierr); 740 PetscFunctionReturn(0); 741 } 742 743 #undef __FUNCT__ 744 #define __FUNCT__ "PCTelescopeSetReductionFactor" 745 /*@ 746 PCTelescopeSetReductionFactor - Sets the factor by which the original number of processes has been reduced by. 747 748 Not Collective 749 750 Input Parameter: 751 . pc - the preconditioner context 752 753 Output Parameter: 754 . fact - the reduction factor 755 756 Level: advanced 757 758 .keywords: PC, telescoping solve 759 @*/ 760 PetscErrorCode PCTelescopeSetReductionFactor(PC pc,PetscInt fact) 761 { 762 PetscErrorCode ierr; 763 PetscFunctionBegin; 764 ierr = PetscTryMethod(pc,"PCTelescopeSetReductionFactor_C",(PC,PetscInt),(pc,fact));CHKERRQ(ierr); 765 PetscFunctionReturn(0); 766 } 767 768 #undef __FUNCT__ 769 #define __FUNCT__ "PCTelescopeGetIgnoreDM" 770 /*@ 771 PCTelescopeGetIgnoreDM - Get the flag indicating if any DM attached to the PC will be used. 772 773 Not Collective 774 775 Input Parameter: 776 . pc - the preconditioner context 777 778 Output Parameter: 779 . v - the flag 780 781 Level: advanced 782 783 .keywords: PC, telescoping solve 784 @*/ 785 PetscErrorCode PCTelescopeGetIgnoreDM(PC pc,PetscBool *v) 786 { 787 PetscErrorCode ierr; 788 PetscFunctionBegin; 789 ierr = PetscUseMethod(pc,"PCTelescopeGetIgnoreDM_C",(PC,PetscBool*),(pc,v));CHKERRQ(ierr); 790 PetscFunctionReturn(0); 791 } 792 793 #undef __FUNCT__ 794 #define __FUNCT__ "PCTelescopeSetIgnoreDM" 795 /*@ 796 PCTelescopeSetIgnoreDM - Set a flag to ignore any DM attached to the PC. 797 798 Not Collective 799 800 Input Parameter: 801 . pc - the preconditioner context 802 803 Output Parameter: 804 . v - Use PETSC_TRUE to ignore any DM 805 806 Level: advanced 807 808 .keywords: PC, telescoping solve 809 @*/ 810 PetscErrorCode PCTelescopeSetIgnoreDM(PC pc,PetscBool v) 811 { 812 PetscErrorCode ierr; 813 PetscFunctionBegin; 814 ierr = PetscTryMethod(pc,"PCTelescopeSetIgnoreDM_C",(PC,PetscBool),(pc,v));CHKERRQ(ierr); 815 PetscFunctionReturn(0); 816 } 817 818 #undef __FUNCT__ 819 #define __FUNCT__ "PCTelescopeGetIgnoreKSPComputeOperators" 820 /*@ 821 PCTelescopeGetIgnoreKSPComputeOperators - Get the flag indicating if KSPComputeOperators will be used. 822 823 Not Collective 824 825 Input Parameter: 826 . pc - the preconditioner context 827 828 Output Parameter: 829 . v - the flag 830 831 Level: advanced 832 833 .keywords: PC, telescoping solve 834 @*/ 835 PetscErrorCode PCTelescopeGetIgnoreKSPComputeOperators(PC pc,PetscBool *v) 836 { 837 PetscErrorCode ierr; 838 PetscFunctionBegin; 839 ierr = PetscUseMethod(pc,"PCTelescopeGetIgnoreKSPComputeOperators_C",(PC,PetscBool*),(pc,v));CHKERRQ(ierr); 840 PetscFunctionReturn(0); 841 } 842 843 #undef __FUNCT__ 844 #define __FUNCT__ "PCTelescopeSetIgnoreKSPComputeOperators" 845 /*@ 846 PCTelescopeSetIgnoreKSPComputeOperators - Set a flag to ignore KSPComputeOperators. 847 848 Not Collective 849 850 Input Parameter: 851 . pc - the preconditioner context 852 853 Output Parameter: 854 . v - Use PETSC_TRUE to ignore the method (if defined) set via KSPSetComputeOperators on pc 855 856 Level: advanced 857 858 .keywords: PC, telescoping solve 859 @*/ 860 PetscErrorCode PCTelescopeSetIgnoreKSPComputeOperators(PC pc,PetscBool v) 861 { 862 PetscErrorCode ierr; 863 PetscFunctionBegin; 864 ierr = PetscTryMethod(pc,"PCTelescopeSetIgnoreKSPComputeOperators_C",(PC,PetscBool),(pc,v));CHKERRQ(ierr); 865 PetscFunctionReturn(0); 866 } 867 868 #undef __FUNCT__ 869 #define __FUNCT__ "PCTelescopeGetDM" 870 /*@ 871 PCTelescopeGetDM - Get the re-partitioned DM attached to the sub KSP. 872 873 Not Collective 874 875 Input Parameter: 876 . pc - the preconditioner context 877 878 Output Parameter: 879 . subdm - The re-partitioned DM 880 881 Level: advanced 882 883 .keywords: PC, telescoping solve 884 @*/ 885 PetscErrorCode PCTelescopeGetDM(PC pc,DM *subdm) 886 { 887 PetscErrorCode ierr; 888 PetscFunctionBegin; 889 ierr = PetscUseMethod(pc,"PCTelescopeGetDM_C",(PC,DM*),(pc,subdm));CHKERRQ(ierr); 890 PetscFunctionReturn(0); 891 } 892 893 #undef __FUNCT__ 894 #define __FUNCT__ "PCTelescopeSetSubcommType" 895 /*@ 896 PCTelescopeSetSubcommType - set subcommunicator type (interlaced or contiguous) 897 898 Logically Collective 899 900 Input Parameter: 901 + pc - the preconditioner context 902 - subcommtype - the subcommunicator type (see PetscSubcommType) 903 904 Level: advanced 905 906 .keywords: PC, telescoping solve 907 908 .seealso: PetscSubcommType, PetscSubcomm, PCTELESCOPE 909 @*/ 910 PetscErrorCode PCTelescopeSetSubcommType(PC pc, PetscSubcommType subcommtype) 911 { 912 PetscErrorCode ierr; 913 PetscFunctionBegin; 914 ierr = PetscTryMethod(pc,"PCTelescopeSetSubcommType_C",(PC,PetscSubcommType),(pc,subcommtype));CHKERRQ(ierr); 915 PetscFunctionReturn(0); 916 } 917 918 #undef __FUNCT__ 919 #define __FUNCT__ "PCTelescopeGetSubcommType" 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 within c 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 B mat 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 are 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 less 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 which are attached to the only B operator. 980 In case where B has a n nullspace attached, these nullspaces vectors are extract from B and mapped into 981 a new 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 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 permuting of B 991 into the ordering defined by the DMDA on c', (ii) extracting the local chunks via MatGetSubMatrices(), (iii) fusing the 992 locally (sequential) matrices defined on the ranks common to c and c' into B' using MatCreateMPIMatConcatenateSeqMat() 993 994 Limitations/improvements 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(). I opted to use P^T.A.P as it appears 999 VecPermute() does not supported for the use case required here. By computing P, I can permute both the operator and RHS in a 1000 consistent manner. 1001 1002 Mapping of vectors is performed this way 1003 Suppose the parent comm size was 4, and we set a reduction factor of 2, thus 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 .seealso: PCTelescopeGetKSP(), PCTelescopeGetDM(), PCTelescopeGetReductionFactor(), PCTelescopeSetReductionFactor(), PCTelescopeGetIgnoreDM(), PCTelescopeSetIgnoreDM(), PCREDUNDANT 1029 M*/ 1030 #undef __FUNCT__ 1031 #define __FUNCT__ "PCCreate_Telescope" 1032 PETSC_EXTERN PetscErrorCode PCCreate_Telescope(PC pc) 1033 { 1034 PetscErrorCode ierr; 1035 struct _PC_Telescope *sred; 1036 1037 PetscFunctionBegin; 1038 ierr = PetscNewLog(pc,&sred);CHKERRQ(ierr); 1039 sred->subcommtype = PETSC_SUBCOMM_INTERLACED; 1040 sred->redfactor = 1; 1041 sred->ignore_dm = PETSC_FALSE; 1042 sred->ignore_kspcomputeoperators = PETSC_FALSE; 1043 pc->data = (void*)sred; 1044 1045 pc->ops->apply = PCApply_Telescope; 1046 pc->ops->applytranspose = NULL; 1047 pc->ops->applyrichardson = PCApplyRichardson_Telescope; 1048 pc->ops->setup = PCSetUp_Telescope; 1049 pc->ops->destroy = PCDestroy_Telescope; 1050 pc->ops->reset = PCReset_Telescope; 1051 pc->ops->setfromoptions = PCSetFromOptions_Telescope; 1052 pc->ops->view = PCView_Telescope; 1053 1054 sred->pctelescope_setup_type = PCTelescopeSetUp_default; 1055 sred->pctelescope_matcreate_type = PCTelescopeMatCreate_default; 1056 sred->pctelescope_matnullspacecreate_type = PCTelescopeMatNullSpaceCreate_default; 1057 sred->pctelescope_reset_type = NULL; 1058 1059 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeGetKSP_C",PCTelescopeGetKSP_Telescope);CHKERRQ(ierr); 1060 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeGetSubcommType_C",PCTelescopeGetSubcommType_Telescope);CHKERRQ(ierr); 1061 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeSetSubcommType_C",PCTelescopeSetSubcommType_Telescope);CHKERRQ(ierr); 1062 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeGetReductionFactor_C",PCTelescopeGetReductionFactor_Telescope);CHKERRQ(ierr); 1063 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeSetReductionFactor_C",PCTelescopeSetReductionFactor_Telescope);CHKERRQ(ierr); 1064 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeGetIgnoreDM_C",PCTelescopeGetIgnoreDM_Telescope);CHKERRQ(ierr); 1065 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeSetIgnoreDM_C",PCTelescopeSetIgnoreDM_Telescope);CHKERRQ(ierr); 1066 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeGetIgnoreKSPComputeOperators_C",PCTelescopeGetIgnoreKSPComputeOperators_Telescope);CHKERRQ(ierr); 1067 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeSetIgnoreKSPComputeOperators_C",PCTelescopeSetIgnoreKSPComputeOperators_Telescope);CHKERRQ(ierr); 1068 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeGetDM_C",PCTelescopeGetDM_Telescope);CHKERRQ(ierr); 1069 PetscFunctionReturn(0); 1070 } 1071