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