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