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