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