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