1 2 /* 3 Code for manipulating distributed regular arrays in parallel. 4 */ 5 6 #include <petsc-private/daimpl.h> /*I "petscdmda.h" I*/ 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "DMDAGetRay" 10 /*@C 11 DMDAGetRay - Returns a vector on process zero that contains a row or column of the values in a DMDA vector 12 13 Collective on DMDA 14 15 Input Parameters: 16 + da - the distributed array 17 . vec - the vector 18 . dir - Cartesian direction, either DMDA_X, DMDA_Y, or DMDA_Z 19 - gp - global grid point number in this direction 20 21 Output Parameters: 22 + newvec - the new vector that can hold the values (size zero on all processes except process 0) 23 - scatter - the VecScatter that will map from the original vector to the slice 24 25 Level: advanced 26 27 Notes: 28 All processors that share the DMDA must call this with the same gp value 29 30 .keywords: distributed array, get, processor subset 31 @*/ 32 PetscErrorCode DMDAGetRay(DM da,DMDADirection dir,PetscInt gp,Vec *newvec,VecScatter *scatter) 33 { 34 PetscMPIInt rank; 35 DM_DA *dd = (DM_DA*)da->data; 36 PetscErrorCode ierr; 37 IS is; 38 AO ao; 39 Vec vec; 40 PetscInt *indices,i,j; 41 42 PetscFunctionBegin; 43 if (dd->dim == 1) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"Cannot get slice from 1d DMDA"); 44 if (dd->dim == 3) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"Cannot get slice from 3d DMDA"); 45 ierr = DMDAGetAO(da,&ao);CHKERRQ(ierr); 46 ierr = MPI_Comm_rank(((PetscObject)da)->comm,&rank);CHKERRQ(ierr); 47 if (!rank) { 48 if (dir == DMDA_Y) { 49 ierr = PetscMalloc(dd->w*dd->M*sizeof(PetscInt),&indices);CHKERRQ(ierr); 50 indices[0] = gp*dd->M*dd->w; 51 for (i=1; i<dd->M*dd->w; i++) {indices[i] = indices[i-1] + 1;} 52 ierr = AOApplicationToPetsc(ao,dd->M*dd->w,indices);CHKERRQ(ierr); 53 ierr = VecCreate(PETSC_COMM_SELF,newvec);CHKERRQ(ierr); 54 ierr = VecSetBlockSize(*newvec,dd->w);CHKERRQ(ierr); 55 ierr = VecSetSizes(*newvec,dd->M*dd->w,PETSC_DETERMINE);CHKERRQ(ierr); 56 ierr = VecSetType(*newvec,VECSEQ);CHKERRQ(ierr); 57 ierr = ISCreateGeneral(PETSC_COMM_SELF,dd->w*dd->M,indices,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 58 } else if (dir == DMDA_X) { 59 ierr = PetscMalloc(dd->w*dd->N*sizeof(PetscInt),&indices);CHKERRQ(ierr); 60 indices[0] = dd->w*gp; 61 for (j=1; j<dd->w; j++) indices[j] = indices[j-1] + 1; 62 for (i=1; i<dd->N; i++) { 63 indices[i*dd->w] = indices[i*dd->w-1] + dd->w*dd->M - dd->w + 1; 64 for (j=1; j<dd->w; j++) indices[i*dd->w + j] = indices[i*dd->w + j - 1] + 1; 65 } 66 ierr = AOApplicationToPetsc(ao,dd->w*dd->N,indices);CHKERRQ(ierr); 67 ierr = VecCreate(PETSC_COMM_SELF,newvec);CHKERRQ(ierr); 68 ierr = VecSetBlockSize(*newvec,dd->w);CHKERRQ(ierr); 69 ierr = VecSetSizes(*newvec,dd->N*dd->w,PETSC_DETERMINE);CHKERRQ(ierr); 70 ierr = VecSetType(*newvec,VECSEQ);CHKERRQ(ierr); 71 ierr = ISCreateGeneral(PETSC_COMM_SELF,dd->w*dd->N,indices,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 72 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unknown DMDADirection"); 73 } else { 74 ierr = VecCreateSeq(PETSC_COMM_SELF,0,newvec);CHKERRQ(ierr); 75 ierr = ISCreateGeneral(PETSC_COMM_SELF,0,0,PETSC_COPY_VALUES,&is);CHKERRQ(ierr); 76 } 77 ierr = DMGetGlobalVector(da,&vec);CHKERRQ(ierr); 78 ierr = VecScatterCreate(vec,is,*newvec,PETSC_NULL,scatter);CHKERRQ(ierr); 79 ierr = DMRestoreGlobalVector(da,&vec);CHKERRQ(ierr); 80 ierr = ISDestroy(&is);CHKERRQ(ierr); 81 PetscFunctionReturn(0); 82 } 83 84 #undef __FUNCT__ 85 #define __FUNCT__ "DMDAGetProcessorSubset" 86 /*@C 87 DMDAGetProcessorSubset - Returns a communicator consisting only of the 88 processors in a DMDA that own a particular global x, y, or z grid point 89 (corresponding to a logical plane in a 3D grid or a line in a 2D grid). 90 91 Collective on DMDA 92 93 Input Parameters: 94 + da - the distributed array 95 . dir - Cartesian direction, either DMDA_X, DMDA_Y, or DMDA_Z 96 - gp - global grid point number in this direction 97 98 Output Parameters: 99 . comm - new communicator 100 101 Level: advanced 102 103 Notes: 104 All processors that share the DMDA must call this with the same gp value 105 106 This routine is particularly useful to compute boundary conditions 107 or other application-specific calculations that require manipulating 108 sets of data throughout a logical plane of grid points. 109 110 .keywords: distributed array, get, processor subset 111 @*/ 112 PetscErrorCode DMDAGetProcessorSubset(DM da,DMDADirection dir,PetscInt gp,MPI_Comm *comm) 113 { 114 MPI_Group group,subgroup; 115 PetscErrorCode ierr; 116 PetscInt i,ict,flag,*owners,xs,xm,ys,ym,zs,zm; 117 PetscMPIInt size,*ranks = PETSC_NULL; 118 DM_DA *dd = (DM_DA*)da->data; 119 120 PetscFunctionBegin; 121 PetscValidHeaderSpecific(da,DM_CLASSID,1); 122 flag = 0; 123 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr); 124 ierr = MPI_Comm_size(((PetscObject)da)->comm,&size);CHKERRQ(ierr); 125 if (dir == DMDA_Z) { 126 if (dd->dim < 3) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_OUTOFRANGE,"DMDA_Z invalid for DMDA dim < 3"); 127 if (gp < 0 || gp > dd->P) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"invalid grid point"); 128 if (gp >= zs && gp < zs+zm) flag = 1; 129 } else if (dir == DMDA_Y) { 130 if (dd->dim == 1) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_OUTOFRANGE,"DMDA_Y invalid for DMDA dim = 1"); 131 if (gp < 0 || gp > dd->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"invalid grid point"); 132 if (gp >= ys && gp < ys+ym) flag = 1; 133 } else if (dir == DMDA_X) { 134 if (gp < 0 || gp > dd->M) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"invalid grid point"); 135 if (gp >= xs && gp < xs+xm) flag = 1; 136 } else SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Invalid direction"); 137 138 ierr = PetscMalloc2(size,PetscInt,&owners,size,PetscMPIInt,&ranks);CHKERRQ(ierr); 139 ierr = MPI_Allgather(&flag,1,MPIU_INT,owners,1,MPIU_INT,((PetscObject)da)->comm);CHKERRQ(ierr); 140 ict = 0; 141 ierr = PetscInfo2(da,"DMDAGetProcessorSubset: dim=%D, direction=%d, procs: ",dd->dim,(int)dir);CHKERRQ(ierr); 142 for (i=0; i<size; i++) { 143 if (owners[i]) { 144 ranks[ict] = i; ict++; 145 ierr = PetscInfo1(da,"%D ",i);CHKERRQ(ierr); 146 } 147 } 148 ierr = PetscInfo(da,"\n");CHKERRQ(ierr); 149 ierr = MPI_Comm_group(((PetscObject)da)->comm,&group);CHKERRQ(ierr); 150 ierr = MPI_Group_incl(group,ict,ranks,&subgroup);CHKERRQ(ierr); 151 ierr = MPI_Comm_create(((PetscObject)da)->comm,subgroup,comm);CHKERRQ(ierr); 152 ierr = MPI_Group_free(&subgroup);CHKERRQ(ierr); 153 ierr = MPI_Group_free(&group);CHKERRQ(ierr); 154 ierr = PetscFree2(owners,ranks);CHKERRQ(ierr); 155 PetscFunctionReturn(0); 156 } 157 158 #undef __FUNCT__ 159 #define __FUNCT__ "DMDAGetProcessorSubsets" 160 /*@C 161 DMDAGetProcessorSubsets - Returns communicators consisting only of the 162 processors in a DMDA adjacent in a particular dimension, 163 corresponding to a logical plane in a 3D grid or a line in a 2D grid. 164 165 Collective on DMDA 166 167 Input Parameters: 168 + da - the distributed array 169 - dir - Cartesian direction, either DMDA_X, DMDA_Y, or DMDA_Z 170 171 Output Parameters: 172 . subcomm - new communicator 173 174 Level: advanced 175 176 Notes: 177 This routine is useful for distributing one-dimensional data in a tensor product grid. 178 179 .keywords: distributed array, get, processor subset 180 @*/ 181 PetscErrorCode DMDAGetProcessorSubsets(DM da, DMDADirection dir, MPI_Comm *subcomm) 182 { 183 MPI_Comm comm; 184 MPI_Group group, subgroup; 185 PetscInt subgroupSize = 0; 186 PetscInt *firstPoints; 187 PetscMPIInt size, *subgroupRanks = PETSC_NULL; 188 PetscInt xs, xm, ys, ym, zs, zm, firstPoint, p; 189 PetscErrorCode ierr; 190 DM_DA *dd = (DM_DA*)da->data; 191 192 PetscFunctionBegin; 193 PetscValidHeaderSpecific(da, DM_CLASSID, 1); 194 comm = ((PetscObject) da)->comm; 195 ierr = DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm);CHKERRQ(ierr); 196 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 197 if (dir == DMDA_Z) { 198 if (dd->dim < 3) SETERRQ(comm,PETSC_ERR_ARG_OUTOFRANGE,"DMDA_Z invalid for DMDA dim < 3"); 199 firstPoint = zs; 200 } else if (dir == DMDA_Y) { 201 if (dd->dim == 1) SETERRQ(comm,PETSC_ERR_ARG_OUTOFRANGE,"DMDA_Y invalid for DMDA dim = 1"); 202 firstPoint = ys; 203 } else if (dir == DMDA_X) { 204 firstPoint = xs; 205 } else SETERRQ(comm,PETSC_ERR_ARG_OUTOFRANGE,"Invalid direction"); 206 207 ierr = PetscMalloc2(size, PetscInt, &firstPoints, size, PetscMPIInt, &subgroupRanks);CHKERRQ(ierr); 208 ierr = MPI_Allgather(&firstPoint, 1, MPIU_INT, firstPoints, 1, MPIU_INT, comm);CHKERRQ(ierr); 209 ierr = PetscInfo2(da,"DMDAGetProcessorSubset: dim=%D, direction=%d, procs: ",dd->dim,(int)dir);CHKERRQ(ierr); 210 for (p = 0; p < size; ++p) { 211 if (firstPoints[p] == firstPoint) { 212 subgroupRanks[subgroupSize++] = p; 213 ierr = PetscInfo1(da, "%D ", p);CHKERRQ(ierr); 214 } 215 } 216 ierr = PetscInfo(da, "\n");CHKERRQ(ierr); 217 ierr = MPI_Comm_group(comm, &group);CHKERRQ(ierr); 218 ierr = MPI_Group_incl(group, subgroupSize, subgroupRanks, &subgroup);CHKERRQ(ierr); 219 ierr = MPI_Comm_create(comm, subgroup, subcomm);CHKERRQ(ierr); 220 ierr = MPI_Group_free(&subgroup);CHKERRQ(ierr); 221 ierr = MPI_Group_free(&group);CHKERRQ(ierr); 222 ierr = PetscFree2(firstPoints, subgroupRanks);CHKERRQ(ierr); 223 PetscFunctionReturn(0); 224 } 225