1 2 #include <petsc/private/pcimpl.h> 3 #include <petscksp.h> /*I "petscksp.h" I*/ 4 #include <petscdm.h> 5 #include <petscdmda.h> 6 7 #include "telescope.h" 8 9 #undef __FUNCT__ 10 #define __FUNCT__ "_DMDADetermineRankFromGlobalIJK" 11 PetscErrorCode _DMDADetermineRankFromGlobalIJK(PetscInt dim,PetscInt i,PetscInt j,PetscInt k, 12 PetscInt Mp,PetscInt Np,PetscInt Pp, 13 PetscInt start_i[],PetscInt start_j[],PetscInt start_k[], 14 PetscInt span_i[],PetscInt span_j[],PetscInt span_k[], 15 PetscMPIInt *_pi,PetscMPIInt *_pj,PetscMPIInt *_pk,PetscMPIInt *rank_re) 16 { 17 PetscInt pi,pj,pk,n; 18 19 PetscFunctionBegin; 20 pi = pj = pk = -1; 21 if (_pi) { 22 for (n=0; n<Mp; n++) { 23 if ( (i >= start_i[n]) && (i < start_i[n]+span_i[n]) ) { 24 pi = n; 25 break; 26 } 27 } 28 if (pi == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"[dmda-ijk] pi cannot be determined : range %D, val %D",Mp,i); 29 *_pi = pi; 30 } 31 32 if (_pj) { 33 for (n=0; n<Np; n++) { 34 if ( (j >= start_j[n]) && (j < start_j[n]+span_j[n]) ) { 35 pj = n; 36 break; 37 } 38 } 39 if (pj == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"[dmda-ijk] pj cannot be determined : range %D, val %D",Np,j); 40 *_pj = pj; 41 } 42 43 if (_pk) { 44 for (n=0; n<Pp; n++) { 45 if ( (k >= start_k[n]) && (k < start_k[n]+span_k[n]) ) { 46 pk = n; 47 break; 48 } 49 } 50 if (pk == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"[dmda-ijk] pk cannot be determined : range %D, val %D",Pp,k); 51 *_pk = pk; 52 } 53 54 switch (dim) { 55 case 1: 56 *rank_re = pi; 57 break; 58 case 2: 59 *rank_re = pi + pj * Mp; 60 break; 61 case 3: 62 *rank_re = pi + pj * Mp + pk * (Mp*Np); 63 break; 64 } 65 PetscFunctionReturn(0); 66 } 67 68 #undef __FUNCT__ 69 #define __FUNCT__ "_DMDADetermineGlobalS0" 70 PetscErrorCode _DMDADetermineGlobalS0(PetscInt dim,PetscMPIInt rank_re,PetscInt Mp_re,PetscInt Np_re,PetscInt Pp_re, 71 PetscInt range_i_re[],PetscInt range_j_re[],PetscInt range_k_re[],PetscInt *s0) 72 { 73 PetscInt i,j,k,start_IJK = 0; 74 PetscInt rank_ijk; 75 76 PetscFunctionBegin; 77 switch (dim) { 78 case 1: 79 for (i=0; i<Mp_re; i++) { 80 rank_ijk = i; 81 if (rank_ijk < rank_re) { 82 start_IJK += range_i_re[i]; 83 } 84 } 85 break; 86 case 2: 87 for (j=0; j<Np_re; j++) { 88 for (i=0; i<Mp_re; i++) { 89 rank_ijk = i + j*Mp_re; 90 if (rank_ijk < rank_re) { 91 start_IJK += range_i_re[i]*range_j_re[j]; 92 } 93 } 94 } 95 break; 96 case 3: 97 for (k=0; k<Pp_re; k++) { 98 for (j=0; j<Np_re; j++) { 99 for (i=0; i<Mp_re; i++) { 100 rank_ijk = i + j*Mp_re + k*Mp_re*Np_re; 101 if (rank_ijk < rank_re) { 102 start_IJK += range_i_re[i]*range_j_re[j]*range_k_re[k]; 103 } 104 } 105 } 106 } 107 break; 108 } 109 *s0 = start_IJK; 110 PetscFunctionReturn(0); 111 } 112 113 #undef __FUNCT__ 114 #define __FUNCT__ "PCTelescopeSetUp_dmda_repart_coors2d" 115 PetscErrorCode PCTelescopeSetUp_dmda_repart_coors2d(PetscSubcomm psubcomm,DM dm,DM subdm) 116 { 117 PetscErrorCode ierr; 118 DM cdm; 119 Vec coor,coor_natural,perm_coors; 120 PetscInt i,j,si,sj,ni,nj,M,N,Ml,Nl,c,nidx; 121 PetscInt *fine_indices; 122 IS is_fine,is_local; 123 VecScatter sctx; 124 125 PetscFunctionBegin; 126 ierr = DMGetCoordinates(dm,&coor);CHKERRQ(ierr); 127 if (!coor) return(0); 128 if (isActiveRank(psubcomm)) { 129 ierr = DMDASetUniformCoordinates(subdm,0.0,1.0,0.0,1.0,0.0,1.0);CHKERRQ(ierr); 130 } 131 /* Get the coordinate vector from the distributed array */ 132 ierr = DMGetCoordinateDM(dm,&cdm);CHKERRQ(ierr); 133 ierr = DMDACreateNaturalVector(cdm,&coor_natural);CHKERRQ(ierr); 134 135 ierr = DMDAGlobalToNaturalBegin(cdm,coor,INSERT_VALUES,coor_natural);CHKERRQ(ierr); 136 ierr = DMDAGlobalToNaturalEnd(cdm,coor,INSERT_VALUES,coor_natural);CHKERRQ(ierr); 137 138 /* get indices of the guys I want to grab */ 139 ierr = DMDAGetInfo(dm,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 140 if (isActiveRank(psubcomm)) { 141 ierr = DMDAGetCorners(subdm,&si,&sj,NULL,&ni,&nj,NULL);CHKERRQ(ierr); 142 Ml = ni; 143 Nl = nj; 144 } else { 145 Ml = Nl = 0; 146 } 147 148 ierr = PetscMalloc1(Ml*Nl*2,&fine_indices);CHKERRQ(ierr); 149 c = 0; 150 if (isActiveRank(psubcomm)) { 151 for (j=sj; j<sj+nj; j++) { 152 for (i=si; i<si+ni; i++) { 153 nidx = (i) + (j)*M; 154 fine_indices[c ] = 2 * nidx ; 155 fine_indices[c+1] = 2 * nidx + 1 ; 156 c = c + 2; 157 } 158 } 159 if (c != Ml*Nl*2) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"c %D should equal 2 * Ml %D * Nl %D",c,Ml,Nl); 160 } 161 162 /* generate scatter */ 163 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm),Ml*Nl*2,fine_indices,PETSC_USE_POINTER,&is_fine);CHKERRQ(ierr); 164 ierr = ISCreateStride(PETSC_COMM_SELF,Ml*Nl*2,0,1,&is_local);CHKERRQ(ierr); 165 166 /* scatter */ 167 ierr = VecCreate(PETSC_COMM_SELF,&perm_coors);CHKERRQ(ierr); 168 ierr = VecSetSizes(perm_coors,PETSC_DECIDE,Ml*Nl*2);CHKERRQ(ierr); 169 ierr = VecSetType(perm_coors,VECSEQ);CHKERRQ(ierr); 170 171 ierr = VecScatterCreate(coor_natural,is_fine,perm_coors,is_local,&sctx);CHKERRQ(ierr); 172 ierr = VecScatterBegin(sctx,coor_natural,perm_coors,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 173 ierr = VecScatterEnd( sctx,coor_natural,perm_coors,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 174 /* access */ 175 if (isActiveRank(psubcomm)) { 176 Vec _coors; 177 const PetscScalar *LA_perm; 178 PetscScalar *LA_coors; 179 180 ierr = DMGetCoordinates(subdm,&_coors);CHKERRQ(ierr); 181 ierr = VecGetArrayRead(perm_coors,&LA_perm);CHKERRQ(ierr); 182 ierr = VecGetArray(_coors,&LA_coors);CHKERRQ(ierr); 183 for (i=0; i<Ml*Nl*2; i++) { 184 LA_coors[i] = LA_perm[i]; 185 } 186 ierr = VecRestoreArray(_coors,&LA_coors);CHKERRQ(ierr); 187 ierr = VecRestoreArrayRead(perm_coors,&LA_perm);CHKERRQ(ierr); 188 } 189 190 /* update local coords */ 191 if (isActiveRank(psubcomm)) { 192 DM _dmc; 193 Vec _coors,_coors_local; 194 ierr = DMGetCoordinateDM(subdm,&_dmc);CHKERRQ(ierr); 195 ierr = DMGetCoordinates(subdm,&_coors);CHKERRQ(ierr); 196 ierr = DMGetCoordinatesLocal(subdm,&_coors_local);CHKERRQ(ierr); 197 ierr = DMGlobalToLocalBegin(_dmc,_coors,INSERT_VALUES,_coors_local);CHKERRQ(ierr); 198 ierr = DMGlobalToLocalEnd(_dmc,_coors,INSERT_VALUES,_coors_local);CHKERRQ(ierr); 199 } 200 ierr = VecScatterDestroy(&sctx);CHKERRQ(ierr); 201 ierr = ISDestroy(&is_fine);CHKERRQ(ierr); 202 ierr = PetscFree(fine_indices);CHKERRQ(ierr); 203 ierr = ISDestroy(&is_local);CHKERRQ(ierr); 204 ierr = VecDestroy(&perm_coors);CHKERRQ(ierr); 205 ierr = VecDestroy(&coor_natural);CHKERRQ(ierr); 206 PetscFunctionReturn(0); 207 } 208 209 #undef __FUNCT__ 210 #define __FUNCT__ "PCTelescopeSetUp_dmda_repart_coors3d" 211 PetscErrorCode PCTelescopeSetUp_dmda_repart_coors3d(PetscSubcomm psubcomm,DM dm,DM subdm) 212 { 213 PetscErrorCode ierr; 214 DM cdm; 215 Vec coor,coor_natural,perm_coors; 216 PetscInt i,j,k,si,sj,sk,ni,nj,nk,M,N,P,Ml,Nl,Pl,c,nidx; 217 PetscInt *fine_indices; 218 IS is_fine,is_local; 219 VecScatter sctx; 220 221 PetscFunctionBegin; 222 ierr = DMGetCoordinates(dm,&coor);CHKERRQ(ierr); 223 if (!coor) PetscFunctionReturn(0); 224 225 if (isActiveRank(psubcomm)) { 226 ierr = DMDASetUniformCoordinates(subdm,0.0,1.0,0.0,1.0,0.0,1.0);CHKERRQ(ierr); 227 } 228 229 /* Get the coordinate vector from the distributed array */ 230 ierr = DMGetCoordinateDM(dm,&cdm);CHKERRQ(ierr); 231 ierr = DMDACreateNaturalVector(cdm,&coor_natural);CHKERRQ(ierr); 232 ierr = DMDAGlobalToNaturalBegin(cdm,coor,INSERT_VALUES,coor_natural);CHKERRQ(ierr); 233 ierr = DMDAGlobalToNaturalEnd(cdm,coor,INSERT_VALUES,coor_natural);CHKERRQ(ierr); 234 235 /* get indices of the guys I want to grab */ 236 ierr = DMDAGetInfo(dm,NULL,&M,&N,&P,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 237 238 if (isActiveRank(psubcomm)) { 239 ierr = DMDAGetCorners(subdm,&si,&sj,&sk,&ni,&nj,&nk);CHKERRQ(ierr); 240 241 Ml = ni - si; 242 Nl = nj - sj; 243 Pl = nk - sk; 244 } else { 245 Ml = Nl = Pl = 0; 246 } 247 248 ierr = PetscMalloc1(Ml*Nl*Pl*3,&fine_indices);CHKERRQ(ierr); 249 250 c = 0; 251 if (isActiveRank(psubcomm)) { 252 for (k=sk; k<sk+nk; k++) { 253 for (j=sj; j<sj+nj; j++) { 254 for (i=si; i<si+ni; i++) { 255 nidx = (i) + (j)*M + (k)*M*N; 256 fine_indices[c ] = 3 * nidx ; 257 fine_indices[c+1] = 3 * nidx + 1 ; 258 fine_indices[c+2] = 3 * nidx + 2 ; 259 c = c + 3; 260 } 261 } 262 } 263 } 264 265 /* generate scatter */ 266 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm),Ml*Nl*Pl*3,fine_indices,PETSC_USE_POINTER,&is_fine);CHKERRQ(ierr); 267 ierr = ISCreateStride(PETSC_COMM_SELF,Ml*Nl*Pl*3,0,1,&is_local);CHKERRQ(ierr); 268 269 /* scatter */ 270 ierr = VecCreate(PETSC_COMM_SELF,&perm_coors);CHKERRQ(ierr); 271 ierr = VecSetSizes(perm_coors,PETSC_DECIDE,Ml*Nl*Pl*3);CHKERRQ(ierr); 272 ierr = VecSetType(perm_coors,VECSEQ);CHKERRQ(ierr); 273 ierr = VecScatterCreate(coor_natural,is_fine,perm_coors,is_local,&sctx);CHKERRQ(ierr); 274 ierr = VecScatterBegin(sctx,coor_natural,perm_coors,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 275 ierr = VecScatterEnd( sctx,coor_natural,perm_coors,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 276 277 /* access */ 278 if (isActiveRank(psubcomm)) { 279 Vec _coors; 280 const PetscScalar *LA_perm; 281 PetscScalar *LA_coors; 282 283 ierr = DMGetCoordinates(subdm,&_coors);CHKERRQ(ierr); 284 ierr = VecGetArrayRead(perm_coors,&LA_perm);CHKERRQ(ierr); 285 ierr = VecGetArray(_coors,&LA_coors);CHKERRQ(ierr); 286 for (i=0; i<Ml*Nl*Pl*3; i++) { 287 LA_coors[i] = LA_perm[i]; 288 } 289 ierr = VecRestoreArray(_coors,&LA_coors);CHKERRQ(ierr); 290 ierr = VecRestoreArrayRead(perm_coors,&LA_perm);CHKERRQ(ierr); 291 } 292 293 /* update local coords */ 294 if (isActiveRank(psubcomm)) { 295 DM _dmc; 296 Vec _coors,_coors_local; 297 298 ierr = DMGetCoordinateDM(subdm,&_dmc);CHKERRQ(ierr); 299 ierr = DMGetCoordinates(subdm,&_coors);CHKERRQ(ierr); 300 ierr = DMGetCoordinatesLocal(subdm,&_coors_local);CHKERRQ(ierr); 301 ierr = DMGlobalToLocalBegin(_dmc,_coors,INSERT_VALUES,_coors_local);CHKERRQ(ierr); 302 ierr = DMGlobalToLocalEnd(_dmc,_coors,INSERT_VALUES,_coors_local);CHKERRQ(ierr); 303 } 304 305 ierr = VecScatterDestroy(&sctx);CHKERRQ(ierr); 306 ierr = ISDestroy(&is_fine);CHKERRQ(ierr); 307 ierr = PetscFree(fine_indices);CHKERRQ(ierr); 308 ierr = ISDestroy(&is_local);CHKERRQ(ierr); 309 ierr = VecDestroy(&perm_coors);CHKERRQ(ierr); 310 ierr = VecDestroy(&coor_natural);CHKERRQ(ierr); 311 PetscFunctionReturn(0); 312 } 313 314 #undef __FUNCT__ 315 #define __FUNCT__ "PCTelescopeSetUp_dmda_repart_coors" 316 PetscErrorCode PCTelescopeSetUp_dmda_repart_coors(PC pc,PC_Telescope sred,PC_Telescope_DMDACtx *ctx) 317 { 318 PetscInt dim; 319 DM dm,subdm; 320 PetscSubcomm psubcomm; 321 MPI_Comm comm; 322 Vec coor; 323 PetscErrorCode ierr; 324 325 PetscFunctionBegin; 326 ierr = PCGetDM(pc,&dm);CHKERRQ(ierr); 327 ierr = DMGetCoordinates(dm,&coor);CHKERRQ(ierr); 328 if (!coor) PetscFunctionReturn(0); 329 psubcomm = sred->psubcomm; 330 comm = PetscSubcommParent(psubcomm); 331 subdm = ctx->dmrepart; 332 333 334 ierr = PetscInfo(pc,"PCTelescope: setting up the coordinates (DMDA)\n");CHKERRQ(ierr); 335 ierr = DMDAGetInfo(dm,&dim,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 336 switch (dim) { 337 case 1: SETERRQ(comm,PETSC_ERR_SUP,"Telescope: DMDA (1D) repartitioning not provided"); 338 break; 339 case 2: PCTelescopeSetUp_dmda_repart_coors2d(psubcomm,dm,subdm); 340 break; 341 case 3: PCTelescopeSetUp_dmda_repart_coors3d(psubcomm,dm,subdm); 342 break; 343 } 344 PetscFunctionReturn(0); 345 } 346 347 /* setup repartitioned dm */ 348 #undef __FUNCT__ 349 #define __FUNCT__ "PCTelescopeSetUp_dmda_repart" 350 PetscErrorCode PCTelescopeSetUp_dmda_repart(PC pc,PC_Telescope sred,PC_Telescope_DMDACtx *ctx) 351 { 352 PetscErrorCode ierr; 353 DM dm; 354 PetscInt dim,nx,ny,nz,ndof,nsw,sum,k; 355 DMBoundaryType bx,by,bz; 356 DMDAStencilType stencil; 357 const PetscInt *_range_i_re; 358 const PetscInt *_range_j_re; 359 const PetscInt *_range_k_re; 360 DMDAInterpolationType itype; 361 PetscInt refine_x,refine_y,refine_z; 362 MPI_Comm comm,subcomm; 363 const char *prefix; 364 365 PetscFunctionBegin; 366 comm = PetscSubcommParent(sred->psubcomm); 367 subcomm = PetscSubcommChild(sred->psubcomm); 368 ierr = PCGetDM(pc,&dm);CHKERRQ(ierr); 369 370 ierr = DMDAGetInfo(dm,&dim,&nx,&ny,&nz,NULL,NULL,NULL,&ndof,&nsw,&bx,&by,&bz,&stencil);CHKERRQ(ierr); 371 ierr = DMDAGetInterpolationType(dm,&itype);CHKERRQ(ierr); 372 ierr = DMDAGetRefinementFactor(dm,&refine_x,&refine_y,&refine_z);CHKERRQ(ierr); 373 374 ctx->dmrepart = NULL; 375 _range_i_re = _range_j_re = _range_k_re = NULL; 376 /* Create DMDA on the child communicator */ 377 if (isActiveRank(sred->psubcomm)) { 378 switch (dim) { 379 case 1: 380 ierr = PetscInfo(pc,"PCTelescope: setting up the DMDA on comm subset (1D)\n");CHKERRQ(ierr); 381 /*ierr = DMDACreate1d(subcomm,bx,nx,ndof,nsw,NULL,&ctx->dmrepart);CHKERRQ(ierr);*/ 382 ny = nz = 1; 383 by = bz = DM_BOUNDARY_NONE; 384 break; 385 case 2: 386 ierr = PetscInfo(pc,"PCTelescope: setting up the DMDA on comm subset (2D)\n");CHKERRQ(ierr); 387 /*ierr = DMDACreate2d(subcomm,bx,by,stencil,nx,ny, PETSC_DECIDE,PETSC_DECIDE, ndof,nsw, NULL,NULL,&ctx->dmrepart);CHKERRQ(ierr);*/ 388 nz = 1; 389 bz = DM_BOUNDARY_NONE; 390 break; 391 case 3: 392 ierr = PetscInfo(pc,"PCTelescope: setting up the DMDA on comm subset (3D)\n");CHKERRQ(ierr); 393 /*ierr = DMDACreate3d(subcomm,bx,by,bz,stencil,nx,ny,nz, PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, ndof,nsw, NULL,NULL,NULL,&ctx->dmrepart);CHKERRQ(ierr);*/ 394 break; 395 } 396 /* 397 The API DMDACreate1d(), DMDACreate2d(), DMDACreate3d() does not allow us to set/append 398 a unique option prefix for the DM, thus I prefer to expose the contents of these API's here. 399 This allows users to control the partitioning of the subDM. 400 */ 401 ierr = DMDACreate(subcomm,&ctx->dmrepart);CHKERRQ(ierr); 402 /* Set unique option prefix name */ 403 ierr = DMGetOptionsPrefix(dm,&prefix);CHKERRQ(ierr); 404 ierr = DMSetOptionsPrefix(ctx->dmrepart,prefix);CHKERRQ(ierr); 405 ierr = DMAppendOptionsPrefix(ctx->dmrepart,"repart_");CHKERRQ(ierr); 406 /* standard setup from DMDACreate{1,2,3}d() */ 407 ierr = DMSetDimension(ctx->dmrepart,dim);CHKERRQ(ierr); 408 ierr = DMDASetSizes(ctx->dmrepart,nx,ny,nz);CHKERRQ(ierr); 409 ierr = DMDASetNumProcs(ctx->dmrepart,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 410 ierr = DMDASetBoundaryType(ctx->dmrepart,bx,by,bz);CHKERRQ(ierr); 411 ierr = DMDASetDof(ctx->dmrepart,ndof);CHKERRQ(ierr); 412 ierr = DMDASetStencilType(ctx->dmrepart,stencil);CHKERRQ(ierr); 413 ierr = DMDASetStencilWidth(ctx->dmrepart,nsw);CHKERRQ(ierr); 414 ierr = DMDASetOwnershipRanges(ctx->dmrepart,NULL,NULL,NULL);CHKERRQ(ierr); 415 ierr = DMSetFromOptions(ctx->dmrepart);CHKERRQ(ierr); 416 ierr = DMSetUp(ctx->dmrepart);CHKERRQ(ierr); 417 /* Set refinement factors and interpolation type from the partent */ 418 ierr = DMDASetRefinementFactor(ctx->dmrepart,refine_x,refine_y,refine_z);CHKERRQ(ierr); 419 ierr = DMDASetInterpolationType(ctx->dmrepart,itype);CHKERRQ(ierr); 420 421 ierr = DMDAGetInfo(ctx->dmrepart,NULL,NULL,NULL,NULL,&ctx->Mp_re,&ctx->Np_re,&ctx->Pp_re,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 422 ierr = DMDAGetOwnershipRanges(ctx->dmrepart,&_range_i_re,&_range_j_re,&_range_k_re);CHKERRQ(ierr); 423 } 424 425 /* generate ranges for repartitioned dm */ 426 /* note - assume rank 0 always participates */ 427 ierr = MPI_Bcast(&ctx->Mp_re,1,MPIU_INT,0,comm);CHKERRQ(ierr); 428 ierr = MPI_Bcast(&ctx->Np_re,1,MPIU_INT,0,comm);CHKERRQ(ierr); 429 ierr = MPI_Bcast(&ctx->Pp_re,1,MPIU_INT,0,comm);CHKERRQ(ierr); 430 431 ierr = PetscCalloc1(ctx->Mp_re,&ctx->range_i_re);CHKERRQ(ierr); 432 ierr = PetscCalloc1(ctx->Np_re,&ctx->range_j_re);CHKERRQ(ierr); 433 ierr = PetscCalloc1(ctx->Pp_re,&ctx->range_k_re);CHKERRQ(ierr); 434 435 if (_range_i_re) {ierr = PetscMemcpy(ctx->range_i_re,_range_i_re,sizeof(PetscInt)*ctx->Mp_re);CHKERRQ(ierr);} 436 if (_range_j_re) {ierr = PetscMemcpy(ctx->range_j_re,_range_j_re,sizeof(PetscInt)*ctx->Np_re);CHKERRQ(ierr);} 437 if (_range_k_re) {ierr = PetscMemcpy(ctx->range_k_re,_range_k_re,sizeof(PetscInt)*ctx->Pp_re);CHKERRQ(ierr);} 438 439 ierr = MPI_Bcast(ctx->range_i_re,ctx->Mp_re,MPIU_INT,0,comm);CHKERRQ(ierr); 440 ierr = MPI_Bcast(ctx->range_j_re,ctx->Np_re,MPIU_INT,0,comm);CHKERRQ(ierr); 441 ierr = MPI_Bcast(ctx->range_k_re,ctx->Pp_re,MPIU_INT,0,comm);CHKERRQ(ierr); 442 443 ierr = PetscMalloc1(ctx->Mp_re,&ctx->start_i_re);CHKERRQ(ierr); 444 ierr = PetscMalloc1(ctx->Np_re,&ctx->start_j_re);CHKERRQ(ierr); 445 ierr = PetscMalloc1(ctx->Pp_re,&ctx->start_k_re);CHKERRQ(ierr); 446 447 sum = 0; 448 for (k=0; k<ctx->Mp_re; k++) { 449 ctx->start_i_re[k] = sum; 450 sum += ctx->range_i_re[k]; 451 } 452 453 sum = 0; 454 for (k=0; k<ctx->Np_re; k++) { 455 ctx->start_j_re[k] = sum; 456 sum += ctx->range_j_re[k]; 457 } 458 459 sum = 0; 460 for (k=0; k<ctx->Pp_re; k++) { 461 ctx->start_k_re[k] = sum; 462 sum += ctx->range_k_re[k]; 463 } 464 465 /* attach dm to ksp on sub communicator */ 466 if (isActiveRank(sred->psubcomm)) { 467 ierr = KSPSetDM(sred->ksp,ctx->dmrepart);CHKERRQ(ierr); 468 ierr = KSPSetDMActive(sred->ksp,PETSC_FALSE);CHKERRQ(ierr); 469 } 470 PetscFunctionReturn(0); 471 } 472 473 #undef __FUNCT__ 474 #define __FUNCT__ "PCTelescopeSetUp_dmda_permutation_3d" 475 PetscErrorCode PCTelescopeSetUp_dmda_permutation_3d(PC pc,PC_Telescope sred,PC_Telescope_DMDACtx *ctx) 476 { 477 PetscErrorCode ierr; 478 DM dm; 479 MPI_Comm comm; 480 Mat Pscalar,P; 481 PetscInt ndof; 482 PetscInt i,j,k,location,startI[3],endI[3],lenI[3],nx,ny,nz; 483 PetscInt sr,er,Mr; 484 Vec V; 485 486 PetscFunctionBegin; 487 ierr = PetscInfo(pc,"PCTelescope: setting up the permutation matrix (DMDA-3D)\n");CHKERRQ(ierr); 488 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 489 490 ierr = PCGetDM(pc,&dm);CHKERRQ(ierr); 491 ierr = DMDAGetInfo(dm,NULL,&nx,&ny,&nz,NULL,NULL,NULL,&ndof,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 492 493 ierr = DMGetGlobalVector(dm,&V);CHKERRQ(ierr); 494 ierr = VecGetSize(V,&Mr);CHKERRQ(ierr); 495 ierr = VecGetOwnershipRange(V,&sr,&er);CHKERRQ(ierr); 496 ierr = DMRestoreGlobalVector(dm,&V);CHKERRQ(ierr); 497 sr = sr / ndof; 498 er = er / ndof; 499 Mr = Mr / ndof; 500 501 ierr = MatCreate(comm,&Pscalar);CHKERRQ(ierr); 502 ierr = MatSetSizes(Pscalar,(er-sr),(er-sr),Mr,Mr);CHKERRQ(ierr); 503 ierr = MatSetType(Pscalar,MATAIJ);CHKERRQ(ierr); 504 ierr = MatSeqAIJSetPreallocation(Pscalar,2,NULL);CHKERRQ(ierr); 505 ierr = MatMPIAIJSetPreallocation(Pscalar,2,NULL,2,NULL);CHKERRQ(ierr); 506 507 ierr = DMDAGetCorners(dm,NULL,NULL,NULL,&lenI[0],&lenI[1],&lenI[2]);CHKERRQ(ierr); 508 ierr = DMDAGetCorners(dm,&startI[0],&startI[1],&startI[2],&endI[0],&endI[1],&endI[2]);CHKERRQ(ierr); 509 endI[0] += startI[0]; 510 endI[1] += startI[1]; 511 endI[2] += startI[2]; 512 513 for (k=startI[2]; k<endI[2]; k++) { 514 for (j=startI[1]; j<endI[1]; j++) { 515 for (i=startI[0]; i<endI[0]; i++) { 516 PetscMPIInt rank_ijk_re,rank_reI[3]; 517 PetscInt s0_re; 518 PetscInt ii,jj,kk,local_ijk_re,mapped_ijk; 519 PetscInt lenI_re[3]; 520 521 location = (i - startI[0]) + (j - startI[1])*lenI[0] + (k - startI[2])*lenI[0]*lenI[1]; 522 ierr = _DMDADetermineRankFromGlobalIJK(3,i,j,k, ctx->Mp_re,ctx->Np_re,ctx->Pp_re, 523 ctx->start_i_re,ctx->start_j_re,ctx->start_k_re, 524 ctx->range_i_re,ctx->range_j_re,ctx->range_k_re, 525 &rank_reI[0],&rank_reI[1],&rank_reI[2],&rank_ijk_re);CHKERRQ(ierr); 526 ierr = _DMDADetermineGlobalS0(3,rank_ijk_re, ctx->Mp_re,ctx->Np_re,ctx->Pp_re, ctx->range_i_re,ctx->range_j_re,ctx->range_k_re, &s0_re);CHKERRQ(ierr); 527 ii = i - ctx->start_i_re[ rank_reI[0] ]; 528 if (ii < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"[dmdarepart-perm3d] index error ii"); 529 jj = j - ctx->start_j_re[ rank_reI[1] ]; 530 if (jj < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"[dmdarepart-perm3d] index error jj"); 531 kk = k - ctx->start_k_re[ rank_reI[2] ]; 532 if (kk < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"[dmdarepart-perm3d] index error kk"); 533 lenI_re[0] = ctx->range_i_re[ rank_reI[0] ]; 534 lenI_re[1] = ctx->range_j_re[ rank_reI[1] ]; 535 lenI_re[2] = ctx->range_k_re[ rank_reI[2] ]; 536 local_ijk_re = ii + jj * lenI_re[0] + kk * lenI_re[0] * lenI_re[1]; 537 mapped_ijk = s0_re + local_ijk_re; 538 ierr = MatSetValue(Pscalar,sr+location,mapped_ijk,1.0,INSERT_VALUES);CHKERRQ(ierr); 539 } 540 } 541 } 542 ierr = MatAssemblyBegin(Pscalar,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 543 ierr = MatAssemblyEnd(Pscalar,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 544 ierr = MatCreateMAIJ(Pscalar,ndof,&P);CHKERRQ(ierr); 545 ierr = MatDestroy(&Pscalar);CHKERRQ(ierr); 546 ctx->permutation = P; 547 PetscFunctionReturn(0); 548 } 549 550 #undef __FUNCT__ 551 #define __FUNCT__ "PCTelescopeSetUp_dmda_permutation_2d" 552 PetscErrorCode PCTelescopeSetUp_dmda_permutation_2d(PC pc,PC_Telescope sred,PC_Telescope_DMDACtx *ctx) 553 { 554 PetscErrorCode ierr; 555 DM dm; 556 MPI_Comm comm; 557 Mat Pscalar,P; 558 PetscInt ndof; 559 PetscInt i,j,location,startI[2],endI[2],lenI[2],nx,ny,nz; 560 PetscInt sr,er,Mr; 561 Vec V; 562 563 PetscFunctionBegin; 564 ierr = PetscInfo(pc,"PCTelescope: setting up the permutation matrix (DMDA-2D)\n");CHKERRQ(ierr); 565 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 566 ierr = PCGetDM(pc,&dm);CHKERRQ(ierr); 567 ierr = DMDAGetInfo(dm,NULL,&nx,&ny,&nz,NULL,NULL,NULL,&ndof,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 568 ierr = DMGetGlobalVector(dm,&V);CHKERRQ(ierr); 569 ierr = VecGetSize(V,&Mr);CHKERRQ(ierr); 570 ierr = VecGetOwnershipRange(V,&sr,&er);CHKERRQ(ierr); 571 ierr = DMRestoreGlobalVector(dm,&V);CHKERRQ(ierr); 572 sr = sr / ndof; 573 er = er / ndof; 574 Mr = Mr / ndof; 575 576 ierr = MatCreate(comm,&Pscalar);CHKERRQ(ierr); 577 ierr = MatSetSizes(Pscalar,(er-sr),(er-sr),Mr,Mr);CHKERRQ(ierr); 578 ierr = MatSetType(Pscalar,MATAIJ);CHKERRQ(ierr); 579 ierr = MatSeqAIJSetPreallocation(Pscalar,2,NULL);CHKERRQ(ierr); 580 ierr = MatMPIAIJSetPreallocation(Pscalar,2,NULL,2,NULL);CHKERRQ(ierr); 581 582 ierr = DMDAGetCorners(dm,NULL,NULL,NULL,&lenI[0],&lenI[1],NULL);CHKERRQ(ierr); 583 ierr = DMDAGetCorners(dm,&startI[0],&startI[1],NULL,&endI[0],&endI[1],NULL);CHKERRQ(ierr); 584 endI[0] += startI[0]; 585 endI[1] += startI[1]; 586 587 for (j=startI[1]; j<endI[1]; j++) { 588 for (i=startI[0]; i<endI[0]; i++) { 589 PetscMPIInt rank_ijk_re,rank_reI[3]; 590 PetscInt s0_re; 591 PetscInt ii,jj,local_ijk_re,mapped_ijk; 592 PetscInt lenI_re[3]; 593 594 location = (i - startI[0]) + (j - startI[1])*lenI[0]; 595 ierr = _DMDADetermineRankFromGlobalIJK(2,i,j,0, ctx->Mp_re,ctx->Np_re,ctx->Pp_re, 596 ctx->start_i_re,ctx->start_j_re,ctx->start_k_re, 597 ctx->range_i_re,ctx->range_j_re,ctx->range_k_re, 598 &rank_reI[0],&rank_reI[1],NULL,&rank_ijk_re);CHKERRQ(ierr); 599 600 ierr = _DMDADetermineGlobalS0(2,rank_ijk_re, ctx->Mp_re,ctx->Np_re,ctx->Pp_re, ctx->range_i_re,ctx->range_j_re,ctx->range_k_re, &s0_re);CHKERRQ(ierr); 601 602 ii = i - ctx->start_i_re[ rank_reI[0] ]; 603 if (ii < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"[dmdarepart-perm2d] index error ii"); 604 jj = j - ctx->start_j_re[ rank_reI[1] ]; 605 if (jj < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"[dmdarepart-perm2d] index error jj"); 606 607 lenI_re[0] = ctx->range_i_re[ rank_reI[0] ]; 608 lenI_re[1] = ctx->range_j_re[ rank_reI[1] ]; 609 local_ijk_re = ii + jj * lenI_re[0]; 610 mapped_ijk = s0_re + local_ijk_re; 611 ierr = MatSetValue(Pscalar,sr+location,mapped_ijk,1.0,INSERT_VALUES);CHKERRQ(ierr); 612 } 613 } 614 ierr = MatAssemblyBegin(Pscalar,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 615 ierr = MatAssemblyEnd(Pscalar,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 616 ierr = MatCreateMAIJ(Pscalar,ndof,&P);CHKERRQ(ierr); 617 ierr = MatDestroy(&Pscalar);CHKERRQ(ierr); 618 ctx->permutation = P; 619 PetscFunctionReturn(0); 620 } 621 622 #undef __FUNCT__ 623 #define __FUNCT__ "PCTelescopeSetUp_dmda_scatters" 624 PetscErrorCode PCTelescopeSetUp_dmda_scatters(PC pc,PC_Telescope sred,PC_Telescope_DMDACtx *ctx) 625 { 626 PetscErrorCode ierr; 627 Vec xred,yred,xtmp,x,xp; 628 VecScatter scatter; 629 IS isin; 630 Mat B; 631 PetscInt m,bs,st,ed; 632 MPI_Comm comm; 633 634 PetscFunctionBegin; 635 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 636 ierr = PCGetOperators(pc,NULL,&B);CHKERRQ(ierr); 637 ierr = MatCreateVecs(B,&x,NULL);CHKERRQ(ierr); 638 ierr = MatGetBlockSize(B,&bs);CHKERRQ(ierr); 639 ierr = VecDuplicate(x,&xp);CHKERRQ(ierr); 640 m = 0; 641 xred = NULL; 642 yred = NULL; 643 if (isActiveRank(sred->psubcomm)) { 644 ierr = DMCreateGlobalVector(ctx->dmrepart,&xred);CHKERRQ(ierr); 645 ierr = VecDuplicate(xred,&yred);CHKERRQ(ierr); 646 ierr = VecGetOwnershipRange(xred,&st,&ed);CHKERRQ(ierr); 647 ierr = ISCreateStride(comm,ed-st,st,1,&isin);CHKERRQ(ierr); 648 ierr = VecGetLocalSize(xred,&m); 649 } else { 650 ierr = VecGetOwnershipRange(x,&st,&ed);CHKERRQ(ierr); 651 ierr = ISCreateStride(comm,0,st,1,&isin);CHKERRQ(ierr); 652 } 653 ierr = ISSetBlockSize(isin,bs);CHKERRQ(ierr); 654 ierr = VecCreate(comm,&xtmp);CHKERRQ(ierr); 655 ierr = VecSetSizes(xtmp,m,PETSC_DECIDE);CHKERRQ(ierr); 656 ierr = VecSetBlockSize(xtmp,bs);CHKERRQ(ierr); 657 ierr = VecSetType(xtmp,((PetscObject)x)->type_name);CHKERRQ(ierr); 658 ierr = VecScatterCreate(x,isin,xtmp,NULL,&scatter);CHKERRQ(ierr); 659 sred->xred = xred; 660 sred->yred = yred; 661 sred->isin = isin; 662 sred->scatter = scatter; 663 sred->xtmp = xtmp; 664 665 ctx->xp = xp; 666 ierr = VecDestroy(&x);CHKERRQ(ierr); 667 PetscFunctionReturn(0); 668 } 669 670 #undef __FUNCT__ 671 #define __FUNCT__ "PCTelescopeSetUp_dmda" 672 PetscErrorCode PCTelescopeSetUp_dmda(PC pc,PC_Telescope sred) 673 { 674 PC_Telescope_DMDACtx *ctx; 675 PetscInt dim; 676 DM dm; 677 MPI_Comm comm; 678 PetscErrorCode ierr; 679 680 PetscFunctionBegin; 681 ierr = PetscInfo(pc,"PCTelescope: setup (DMDA)\n");CHKERRQ(ierr); 682 PetscMalloc1(1,&ctx); 683 PetscMemzero(ctx,sizeof(PC_Telescope_DMDACtx)); 684 sred->dm_ctx = (void*)ctx; 685 686 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 687 ierr = PCGetDM(pc,&dm);CHKERRQ(ierr); 688 ierr = DMDAGetInfo(dm,&dim,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 689 690 PCTelescopeSetUp_dmda_repart(pc,sred,ctx); 691 PCTelescopeSetUp_dmda_repart_coors(pc,sred,ctx); 692 switch (dim) { 693 case 1: SETERRQ(comm,PETSC_ERR_SUP,"Telescope: DMDA (1D) repartitioning not provided"); 694 break; 695 case 2: ierr = PCTelescopeSetUp_dmda_permutation_2d(pc,sred,ctx);CHKERRQ(ierr); 696 break; 697 case 3: ierr = PCTelescopeSetUp_dmda_permutation_3d(pc,sred,ctx);CHKERRQ(ierr); 698 break; 699 } 700 ierr = PCTelescopeSetUp_dmda_scatters(pc,sred,ctx);CHKERRQ(ierr); 701 PetscFunctionReturn(0); 702 } 703 704 #undef __FUNCT__ 705 #define __FUNCT__ "PCTelescopeMatCreate_dmda" 706 PetscErrorCode PCTelescopeMatCreate_dmda(PC pc,PC_Telescope sred,MatReuse reuse,Mat *A) 707 { 708 PetscErrorCode ierr; 709 PC_Telescope_DMDACtx *ctx; 710 MPI_Comm comm,subcomm; 711 Mat Bperm,Bred,B,P; 712 PetscInt nr,nc; 713 IS isrow,iscol; 714 Mat Blocal,*_Blocal; 715 716 PetscFunctionBegin; 717 ierr = PetscInfo(pc,"PCTelescope: updating the redundant preconditioned operator (DMDA)\n");CHKERRQ(ierr); 718 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 719 subcomm = PetscSubcommChild(sred->psubcomm); 720 ctx = (PC_Telescope_DMDACtx*)sred->dm_ctx; 721 722 ierr = PCGetOperators(pc,NULL,&B);CHKERRQ(ierr); 723 ierr = MatGetSize(B,&nr,&nc);CHKERRQ(ierr); 724 725 P = ctx->permutation; 726 ierr = MatPtAP(B,P,MAT_INITIAL_MATRIX,1.1,&Bperm);CHKERRQ(ierr); 727 728 /* Get submatrices */ 729 isrow = sred->isin; 730 ierr = ISCreateStride(comm,nc,0,1,&iscol);CHKERRQ(ierr); 731 732 ierr = MatGetSubMatrices(Bperm,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&_Blocal);CHKERRQ(ierr); 733 Blocal = *_Blocal; 734 Bred = NULL; 735 if (isActiveRank(sred->psubcomm)) { 736 PetscInt mm; 737 738 if (reuse != MAT_INITIAL_MATRIX) {Bred = *A;} 739 ierr = MatGetSize(Blocal,&mm,NULL);CHKERRQ(ierr); 740 /* ierr = MatCreateMPIMatConcatenateSeqMat(subcomm,Blocal,PETSC_DECIDE,reuse,&Bred);CHKERRQ(ierr); */ 741 ierr = MatCreateMPIMatConcatenateSeqMat(subcomm,Blocal,mm,reuse,&Bred);CHKERRQ(ierr); 742 } 743 *A = Bred; 744 745 ierr = ISDestroy(&iscol);CHKERRQ(ierr); 746 ierr = MatDestroy(&Bperm);CHKERRQ(ierr); 747 ierr = MatDestroyMatrices(1,&_Blocal);CHKERRQ(ierr); 748 PetscFunctionReturn(0); 749 } 750 751 #undef __FUNCT__ 752 #define __FUNCT__ "PCTelescopeMatNullSpaceCreate_dmda" 753 PetscErrorCode PCTelescopeMatNullSpaceCreate_dmda(PC pc,PC_Telescope sred,Mat sub_mat) 754 { 755 PetscErrorCode ierr; 756 MatNullSpace nullspace,sub_nullspace; 757 Mat A,B; 758 PetscBool has_const; 759 PetscInt i,k,n; 760 const Vec *vecs; 761 Vec *sub_vecs; 762 MPI_Comm subcomm; 763 PC_Telescope_DMDACtx *ctx; 764 765 PetscFunctionBegin; 766 ierr = PCGetOperators(pc,&A,&B);CHKERRQ(ierr); 767 ierr = MatGetNullSpace(B,&nullspace);CHKERRQ(ierr); 768 if (!nullspace) return(0); 769 770 ierr = PetscInfo(pc,"PCTelescope: generating nullspace (DMDA)\n");CHKERRQ(ierr); 771 ctx = (PC_Telescope_DMDACtx*)sred->dm_ctx; 772 subcomm = PetscSubcommChild(sred->psubcomm); 773 ierr = MatNullSpaceGetVecs(nullspace,&has_const,&n,&vecs);CHKERRQ(ierr); 774 775 if (isActiveRank(sred->psubcomm)) { 776 /* create new vectors */ 777 if (n) { 778 ierr = VecDuplicateVecs(sred->xred,n,&sub_vecs);CHKERRQ(ierr); 779 } 780 } 781 782 /* copy entries */ 783 for (k=0; k<n; k++) { 784 const PetscScalar *x_array; 785 PetscScalar *LA_sub_vec; 786 PetscInt st,ed,bs; 787 788 /* permute vector into ordering associated with re-partitioned dmda */ 789 ierr = MatMultTranspose(ctx->permutation,vecs[k],ctx->xp);CHKERRQ(ierr); 790 791 /* pull in vector x->xtmp */ 792 ierr = VecScatterBegin(sred->scatter,ctx->xp,sred->xtmp,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 793 ierr = VecScatterEnd(sred->scatter,ctx->xp,sred->xtmp,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 794 795 /* copy vector entires into xred */ 796 ierr = VecGetBlockSize(sred->xtmp,&bs);CHKERRQ(ierr); 797 ierr = VecGetArrayRead(sred->xtmp,&x_array);CHKERRQ(ierr); 798 if (sub_vecs[k]) { 799 ierr = VecGetOwnershipRange(sub_vecs[k],&st,&ed);CHKERRQ(ierr); 800 ierr = VecGetArray(sub_vecs[k],&LA_sub_vec);CHKERRQ(ierr); 801 for (i=0; i<ed-st; i++) { 802 LA_sub_vec[i] = x_array[i]; 803 } 804 ierr = VecRestoreArray(sub_vecs[k],&LA_sub_vec);CHKERRQ(ierr); 805 } 806 ierr = VecRestoreArrayRead(sred->xtmp,&x_array);CHKERRQ(ierr); 807 } 808 809 if (isActiveRank(sred->psubcomm)) { 810 /* create new nullspace for redundant object */ 811 ierr = MatNullSpaceCreate(subcomm,has_const,n,sub_vecs,&sub_nullspace);CHKERRQ(ierr); 812 813 /* attach redundant nullspace to Bred */ 814 ierr = MatSetNullSpace(sub_mat,sub_nullspace);CHKERRQ(ierr); 815 ierr = VecDestroyVecs(n,&sub_vecs);CHKERRQ(ierr); 816 } 817 PetscFunctionReturn(0); 818 } 819 820 #undef __FUNCT__ 821 #define __FUNCT__ "PCApply_Telescope_dmda" 822 PetscErrorCode PCApply_Telescope_dmda(PC pc,Vec x,Vec y) 823 { 824 PC_Telescope sred = (PC_Telescope)pc->data; 825 PetscErrorCode ierr; 826 Mat perm; 827 Vec xtmp,xp,xred,yred; 828 PetscInt i,st,ed,bs; 829 VecScatter scatter; 830 PetscScalar *array; 831 const PetscScalar *x_array; 832 PC_Telescope_DMDACtx *ctx; 833 834 ctx = (PC_Telescope_DMDACtx*)sred->dm_ctx; 835 xtmp = sred->xtmp; 836 scatter = sred->scatter; 837 xred = sred->xred; 838 yred = sred->yred; 839 perm = ctx->permutation; 840 xp = ctx->xp; 841 842 PetscFunctionBegin; 843 /* permute vector into ordering associated with re-partitioned dmda */ 844 ierr = MatMultTranspose(perm,x,xp);CHKERRQ(ierr); 845 846 /* pull in vector x->xtmp */ 847 ierr = VecScatterBegin(scatter,xp,xtmp,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 848 ierr = VecScatterEnd(scatter,xp,xtmp,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 849 850 /* copy vector entires into xred */ 851 ierr = VecGetBlockSize(xtmp,&bs);CHKERRQ(ierr); 852 ierr = VecGetArrayRead(xtmp,&x_array);CHKERRQ(ierr); 853 if (xred) { 854 PetscScalar *LA_xred; 855 ierr = VecGetOwnershipRange(xred,&st,&ed);CHKERRQ(ierr); 856 857 ierr = VecGetArray(xred,&LA_xred);CHKERRQ(ierr); 858 for (i=0; i<ed-st; i++) { 859 LA_xred[i] = x_array[i]; 860 } 861 ierr = VecRestoreArray(xred,&LA_xred);CHKERRQ(ierr); 862 } 863 ierr = VecRestoreArrayRead(xtmp,&x_array);CHKERRQ(ierr); 864 865 /* solve */ 866 if (isActiveRank(sred->psubcomm)) { 867 ierr = KSPSolve(sred->ksp,xred,yred);CHKERRQ(ierr); 868 } 869 870 /* return vector */ 871 ierr = VecGetBlockSize(xtmp,&bs);CHKERRQ(ierr); 872 ierr = VecGetArray(xtmp,&array);CHKERRQ(ierr); 873 if (yred) { 874 const PetscScalar *LA_yred; 875 ierr = VecGetOwnershipRange(yred,&st,&ed);CHKERRQ(ierr); 876 ierr = VecGetArrayRead(yred,&LA_yred);CHKERRQ(ierr); 877 for (i=0; i<ed-st; i++) { 878 array[i] = LA_yred[i]; 879 } 880 ierr = VecRestoreArrayRead(yred,&LA_yred);CHKERRQ(ierr); 881 } 882 ierr = VecRestoreArray(xtmp,&array);CHKERRQ(ierr); 883 ierr = VecScatterBegin(scatter,xtmp,xp,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 884 ierr = VecScatterEnd(scatter,xtmp,xp,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 885 ierr = MatMult(perm,xp,y);CHKERRQ(ierr); 886 PetscFunctionReturn(0); 887 } 888 889 #undef __FUNCT__ 890 #define __FUNCT__ "PCReset_Telescope_dmda" 891 PetscErrorCode PCReset_Telescope_dmda(PC pc) 892 { 893 PetscErrorCode ierr; 894 PC_Telescope sred = (PC_Telescope)pc->data; 895 PC_Telescope_DMDACtx *ctx; 896 897 PetscFunctionBegin; 898 ctx = (PC_Telescope_DMDACtx*)sred->dm_ctx; 899 ierr = VecDestroy(&ctx->xp);CHKERRQ(ierr); 900 ierr = MatDestroy(&ctx->permutation);CHKERRQ(ierr); 901 ierr = DMDestroy(&ctx->dmrepart);CHKERRQ(ierr); 902 ierr = PetscFree(ctx->range_i_re);CHKERRQ(ierr); 903 ierr = PetscFree(ctx->range_j_re);CHKERRQ(ierr); 904 ierr = PetscFree(ctx->range_k_re);CHKERRQ(ierr); 905 ierr = PetscFree(ctx->start_i_re);CHKERRQ(ierr); 906 ierr = PetscFree(ctx->start_j_re);CHKERRQ(ierr); 907 ierr = PetscFree(ctx->start_k_re);CHKERRQ(ierr); 908 PetscFunctionReturn(0); 909 } 910 911 #undef __FUNCT__ 912 #define __FUNCT__ "DMView_DMDAShort_3d" 913 PetscErrorCode DMView_DMDAShort_3d(DM dm,PetscViewer v) 914 { 915 PetscInt M,N,P,m,n,p,ndof,nsw; 916 MPI_Comm comm; 917 PetscMPIInt size; 918 const char* prefix; 919 PetscErrorCode ierr; 920 921 PetscFunctionBegin; 922 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 923 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 924 ierr = DMGetOptionsPrefix(dm,&prefix);CHKERRQ(ierr); 925 ierr = DMDAGetInfo(dm,NULL,&M,&N,&P,&m,&n,&p,&ndof,&nsw,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 926 if (prefix) {ierr = PetscViewerASCIIPrintf(v,"DMDA Object: (%s) %d MPI processes\n",prefix,size);CHKERRQ(ierr);} 927 else {ierr = PetscViewerASCIIPrintf(v,"DMDA Object: %d MPI processes\n",size);CHKERRQ(ierr);} 928 ierr = PetscViewerASCIIPrintf(v," M %D N %D P %D m %D n %D p %D dof %D overlap %D\n",M,N,P,m,n,p,ndof,nsw);CHKERRQ(ierr); 929 PetscFunctionReturn(0); 930 } 931 932 #undef __FUNCT__ 933 #define __FUNCT__ "DMView_DMDAShort_2d" 934 PetscErrorCode DMView_DMDAShort_2d(DM dm,PetscViewer v) 935 { 936 PetscInt M,N,m,n,ndof,nsw; 937 MPI_Comm comm; 938 PetscMPIInt size; 939 const char* prefix; 940 PetscErrorCode ierr; 941 942 PetscFunctionBegin; 943 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 944 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 945 ierr = DMGetOptionsPrefix(dm,&prefix);CHKERRQ(ierr); 946 ierr = DMDAGetInfo(dm,NULL,&M,&N,NULL,&m,&n,NULL,&ndof,&nsw,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 947 if (prefix) {PetscViewerASCIIPrintf(v,"DMDA Object: (%s) %d MPI processes\n",prefix,size);CHKERRQ(ierr);} 948 else {ierr = PetscViewerASCIIPrintf(v,"DMDA Object: %d MPI processes\n",size);CHKERRQ(ierr);} 949 ierr = PetscViewerASCIIPrintf(v," M %D N %D m %D n %D dof %D overlap %D\n",M,N,m,n,ndof,nsw);CHKERRQ(ierr); 950 PetscFunctionReturn(0); 951 } 952 953 #undef __FUNCT__ 954 #define __FUNCT__ "DMView_DMDAShort" 955 PetscErrorCode DMView_DMDAShort(DM dm,PetscViewer v) 956 { 957 PetscErrorCode ierr; 958 PetscInt dim; 959 960 PetscFunctionBegin; 961 ierr = DMDAGetInfo(dm,&dim,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 962 switch (dim) { 963 case 2: ierr = DMView_DMDAShort_2d(dm,v);CHKERRQ(ierr); 964 break; 965 case 3: ierr = DMView_DMDAShort_3d(dm,v);CHKERRQ(ierr); 966 break; 967 } 968 PetscFunctionReturn(0); 969 } 970 971