Actual source code: ghost.c

  1: #ifdef PETSC_RCS_HEADER
  2: static char vcid[] = "$Id: ghost.c,v 1.3 2000/01/10 03:27:01 knepley Exp $";
  3: #endif

 5:  #include petscsys.h

  7: #undef  __FUNCT__
  9: /*@C
 10:   PetscGhostExchange - This functions transfers data between local and ghost storage without a predefined mapping.

 12:   Collective on MPI_Comm

 14:   Input Parameters:
 15: + comm         - The communicator
 16: . numGhosts    - The number of ghosts in this domain
 17: . ghostProcs   - The processor from which to obtain each ghost
 18: . ghostIndices - The global index for each ghost
 19: . dataType     - The type of the variables
 20: . firstVar     - The first variable on each processor
 21: . addv         - The insert mode, INSERT_VALUES or ADD_VALUES
 22: - mode         - The direction of the transfer, SCATTER_FORWARD or SCATTER_REVERSE

 24:   Output Parameters:
 25: + locVars      - The local variable array
 26: - ghostVars    - The ghost variables

 28:   Note:
 29:   The data in ghostVars is assumed contiguous and implicitly indexed by the order of
 30:   ghostProcs and ghostIndices. The SCATTER_FORWARD mode will take the requested data
 31:   from locVars and copy it to ghostVars in the order specified by ghostIndices. The
 32:   SCATTER_REVERSE mode will take data from ghostVars and copy it to locVars.

 34:   Level: advanced

 36: .keywords: ghost, exchange
 37: .seealso: GridGlobalToLocal(), GridLocalToGlobal()
 38: @*/
 39: int PetscGhostExchange(MPI_Comm comm, int numGhosts, int *ghostProcs, int *ghostIndices, PetscDataType dataType,
 40:                       int *firstVar, InsertMode addv, ScatterMode mode, void *locVars, void *ghostVars)
 41: {
 42:   int         *numSendGhosts; /* The number of ghosts from each domain */
 43:   int         *numRecvGhosts; /* The number of local variables which are ghosts in each domain */
 44:   int         *sumSendGhosts; /* The prefix sums of numSendGhosts */
 45:   int         *sumRecvGhosts; /* The prefix sums of numRecvGhosts */
 46:   int         *offsets;       /* The offset into the send array for each domain */
 47:   int          totSendGhosts; /* The number of ghosts to request variables for */
 48:   int          totRecvGhosts; /* The number of nodes to provide class info about */
 49:   int         *sendIndices;   /* The canonical indices of ghosts in this domain */
 50:   int         *recvIndices;   /* The canonical indices of ghosts to return variables for */
 51:   char        *tempVars;      /* The variables of the requested or submitted ghosts */
 52:   char        *locBytes   = (char *) locVars;
 53:   MPI_Datatype MPIType;
 54:   int          typeSize;
 55: #ifdef PETSC_USE_BOPT_g
 56:   int          numLocVars;
 57: #endif
 58:   int          numProcs, rank;
 59:   int          proc, ghost, locIndex, byte;
 60:   int          ierr;

 63:   /* Initialize communication */
 64:   MPI_Comm_size(comm, &numProcs);
 65:   MPI_Comm_rank(comm, &rank);
 66:   PetscMalloc(numProcs * sizeof(int), &numSendGhosts);
 67:   PetscMalloc(numProcs * sizeof(int), &numRecvGhosts);
 68:   PetscMalloc(numProcs * sizeof(int), &sumSendGhosts);
 69:   PetscMalloc(numProcs * sizeof(int), &sumRecvGhosts);
 70:   PetscMalloc(numProcs * sizeof(int), &offsets);
 71:   PetscMemzero(numSendGhosts,  numProcs * sizeof(int));
 72:   PetscMemzero(numRecvGhosts,  numProcs * sizeof(int));
 73:   PetscMemzero(sumSendGhosts,  numProcs * sizeof(int));
 74:   PetscMemzero(sumRecvGhosts,  numProcs * sizeof(int));
 75:   PetscMemzero(offsets,        numProcs * sizeof(int));
 76: #ifdef PETSC_USE_BOPT_g
 77:   numLocVars = firstVar[rank+1] - firstVar[rank];
 78: #endif

 80:   /* Get number of ghosts needed from each processor */
 81:   for(ghost = 0; ghost < numGhosts; ghost++) {
 82:     numSendGhosts[ghostProcs[ghost]]++;
 83:   }

 85:   /* Get number of ghosts to provide variables for */
 86:   MPI_Alltoall(numSendGhosts, 1, MPI_INT, numRecvGhosts, 1, MPI_INT, comm);
 87:   for(proc = 1; proc < numProcs; proc++) {
 88:     sumSendGhosts[proc] = sumSendGhosts[proc-1] + numSendGhosts[proc-1];
 89:     sumRecvGhosts[proc] = sumRecvGhosts[proc-1] + numRecvGhosts[proc-1];
 90:     offsets[proc]       = sumSendGhosts[proc];
 91:   }
 92:   totSendGhosts = sumSendGhosts[numProcs-1] + numSendGhosts[numProcs-1];
 93:   totRecvGhosts = sumRecvGhosts[numProcs-1] + numRecvGhosts[numProcs-1];
 94:   if (numGhosts != totSendGhosts) {
 95:     SETERRQ2(PETSC_ERR_PLIB, "Invalid number of ghosts %d in send, should be %d", totSendGhosts, numGhosts);
 96:   }

 98:   PetscDataTypeGetSize(dataType, &typeSize);
 99:   if (totSendGhosts) {
100:     PetscMalloc(totSendGhosts * sizeof(int), &sendIndices);
101:   }
102:   if (totRecvGhosts) {
103:     PetscMalloc(totRecvGhosts * sizeof(int), &recvIndices);
104:     PetscMalloc(totRecvGhosts * typeSize,    &tempVars);
105:   }

107:   /* Must order ghosts by processor */
108:   for(ghost = 0; ghost < numGhosts; ghost++) {
109:     sendIndices[offsets[ghostProcs[ghost]]++] = ghostIndices[ghost];
110:   }

112:   /* Get canonical indices of ghosts to provide variables for */
113:   MPI_Alltoallv(sendIndices, numSendGhosts, sumSendGhosts, MPI_INT,
114:                        recvIndices, numRecvGhosts, sumRecvGhosts, MPI_INT, comm);
115: 

117:   switch(mode)
118:   {
119:   case SCATTER_FORWARD:
120:     /* Get ghost variables */
121:     if (addv == INSERT_VALUES) {
122:       for(ghost = 0; ghost < totRecvGhosts; ghost++) {
123:         locIndex = recvIndices[ghost] - firstVar[rank];
124: #ifdef PETSC_USE_BOPT_g
125:         if ((locIndex < 0) || (locIndex >= numLocVars)) {
126:           SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid ghost index %d, not in [0,%d)", locIndex, numLocVars);
127:         }
128: #endif
129:         for(byte = 0; byte < typeSize; byte++) {
130:           tempVars[ghost*typeSize+byte] = locBytes[locIndex*typeSize+byte];
131:         }
132:       }
133:     } else {
134:       for(ghost = 0; ghost < totRecvGhosts; ghost++) {
135:         locIndex = recvIndices[ghost] - firstVar[rank];
136: #ifdef PETSC_USE_BOPT_g
137:         if ((locIndex < 0) || (locIndex >= numLocVars)) {
138:           SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid ghost index %d, not in [0,%d)", locIndex, numLocVars);
139:         }
140: #endif
141:         for(byte = 0; byte < typeSize; byte++) {
142:           tempVars[ghost*typeSize+byte] += locBytes[locIndex*typeSize+byte];
143:         }
144:       }
145:     }

147:     /* Communicate local variables to ghost storage */
148:     PetscDataTypeToMPIDataType(dataType, &MPIType);
149:     MPI_Alltoallv(tempVars,  numRecvGhosts, sumRecvGhosts, MPIType,
150:                          ghostVars, numSendGhosts, sumSendGhosts, MPIType, comm);
151: 
152:     break;
153:   case SCATTER_REVERSE:
154:     /* Communicate ghost variables to local storage */
155:     PetscDataTypeToMPIDataType(dataType, &MPIType);
156:     MPI_Alltoallv(ghostVars, numSendGhosts, sumSendGhosts, MPIType,
157:                          tempVars,  numRecvGhosts, sumRecvGhosts, MPIType, comm);
158: 

160:     /* Get ghost variables */
161:     if (addv == INSERT_VALUES) {
162:       for(ghost = 0; ghost < totRecvGhosts; ghost++) {
163:         locIndex = recvIndices[ghost] - firstVar[rank];
164: #ifdef PETSC_USE_BOPT_g
165:         if ((locIndex < 0) || (locIndex >= numLocVars)) {
166:           SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid ghost index %d, not in [0,%d)", locIndex, numLocVars);
167:         }
168: #endif
169:         for(byte = 0; byte < typeSize; byte++) {
170:           locBytes[locIndex*typeSize+byte] = tempVars[ghost*typeSize+byte];
171:         }
172:       }
173:     } else {
174:       /* There must be a better way to do this -- Ask Bill */
175:       if (typeSize == sizeof(int)) {
176:         int *tempInt = (int *) tempVars;
177:         int *locInt  = (int *) locVars;

179:         for(ghost = 0; ghost < totRecvGhosts; ghost++) {
180:           locIndex = recvIndices[ghost] - firstVar[rank];
181: #ifdef PETSC_USE_BOPT_g
182:           if ((locIndex < 0) || (locIndex >= numLocVars)) {
183:             SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid ghost index %d, not in [0,%d)", locIndex, numLocVars);
184:           }
185: #endif
186:           locInt[locIndex] += tempInt[ghost];
187:         }
188:       } else if (typeSize == sizeof(long int)) {
189:         long int *tempLongInt = (long int *) tempVars;
190:         long int *locLongInt  = (long int *) locVars;

192:         for(ghost = 0; ghost < totRecvGhosts; ghost++) {
193:           locIndex = recvIndices[ghost] - firstVar[rank];
194: #ifdef PETSC_USE_BOPT_g
195:           if ((locIndex < 0) || (locIndex >= numLocVars)) {
196:             SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid ghost index %d, not in [0,%d)", locIndex, numLocVars);
197:           }
198: #endif
199:           locLongInt[locIndex] += tempLongInt[ghost];
200:         }
201:       } else {
202:         for(ghost = 0; ghost < totRecvGhosts; ghost++) {
203:           locIndex = recvIndices[ghost] - firstVar[rank];
204: #ifdef PETSC_USE_BOPT_g
205:           if ((locIndex < 0) || (locIndex >= numLocVars)) {
206:             SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid ghost index %d, not in [0,%d)", locIndex, numLocVars);
207:           }
208: #endif
209:           for(byte = 0; byte < typeSize; byte++) {
210:             locBytes[locIndex*typeSize+byte] += tempVars[ghost*typeSize+byte];
211:           }
212:         }
213:       }
214:     }
215:     break;
216:   default:
217:     SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid scatter mode %d", mode);
218:   }

220:   /* Cleanup */
221:   PetscFree(numSendGhosts);
222:   PetscFree(numRecvGhosts);
223:   PetscFree(sumSendGhosts);
224:   PetscFree(sumRecvGhosts);
225:   PetscFree(offsets);
226:   if (totSendGhosts) {
227:     PetscFree(sendIndices);
228:   }
229:   if (totRecvGhosts) {
230:     PetscFree(recvIndices);
231:     PetscFree(tempVars);
232:   }
233:   return(0);
234: }