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