Actual source code: aomapping.c


  2: /*
  3:   These AO application ordering routines do not require that the input
  4:   be a permutation, but merely a 1-1 mapping. This implementation still
  5:   keeps the entire ordering on each processor.
  6: */

  8: #include <../src/vec/is/ao/aoimpl.h>

 10: typedef struct {
 11:   PetscInt N;
 12:   PetscInt *app;       /* app[i] is the partner for petsc[appPerm[i]] */
 13:   PetscInt *appPerm;
 14:   PetscInt *petsc;     /* petsc[j] is the partner for app[petscPerm[j]] */
 15:   PetscInt *petscPerm;
 16: } AO_Mapping;

 18: PetscErrorCode AODestroy_Mapping(AO ao)
 19: {
 20:   AO_Mapping     *aomap = (AO_Mapping*) ao->data;

 22:   PetscFree4(aomap->app,aomap->appPerm,aomap->petsc,aomap->petscPerm);
 23:   PetscFree(aomap);
 24:   return 0;
 25: }

 27: PetscErrorCode AOView_Mapping(AO ao, PetscViewer viewer)
 28: {
 29:   AO_Mapping     *aomap = (AO_Mapping*) ao->data;
 30:   PetscMPIInt    rank;
 31:   PetscInt       i;
 32:   PetscBool      iascii;

 34:   MPI_Comm_rank(PetscObjectComm((PetscObject)ao), &rank);
 35:   if (rank) return 0;
 36:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
 37:   if (iascii) {
 38:     PetscViewerASCIIPrintf(viewer, "Number of elements in ordering %" PetscInt_FMT "\n", aomap->N);
 39:     PetscViewerASCIIPrintf(viewer, "   App.   PETSc\n");
 40:     for (i = 0; i < aomap->N; i++) {
 41:       PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "   %" PetscInt_FMT "    %" PetscInt_FMT "\n", i, aomap->app[i], aomap->petsc[aomap->appPerm[i]]);
 42:     }
 43:   }
 44:   return 0;
 45: }

 47: PetscErrorCode AOPetscToApplication_Mapping(AO ao, PetscInt n, PetscInt *ia)
 48: {
 49:   AO_Mapping *aomap = (AO_Mapping*) ao->data;
 50:   PetscInt   *app   = aomap->app;
 51:   PetscInt   *petsc = aomap->petsc;
 52:   PetscInt   *perm  = aomap->petscPerm;
 53:   PetscInt   N      = aomap->N;
 54:   PetscInt   low, high, mid=0;
 55:   PetscInt   idex;
 56:   PetscInt   i;

 58:   /* It would be possible to use a single bisection search, which
 59:      recursively divided the indices to be converted, and searched
 60:      partitions which contained an index. This would result in
 61:      better running times if indices are clustered.
 62:   */
 63:   for (i = 0; i < n; i++) {
 64:     idex = ia[i];
 65:     if (idex < 0) continue;
 66:     /* Use bisection since the array is sorted */
 67:     low  = 0;
 68:     high = N - 1;
 69:     while (low <= high) {
 70:       mid = (low + high)/2;
 71:       if (idex == petsc[mid]) break;
 72:       else if (idex < petsc[mid]) high = mid - 1;
 73:       else low = mid + 1;
 74:     }
 75:     if (low > high) ia[i] = -1;
 76:     else            ia[i] = app[perm[mid]];
 77:   }
 78:   return 0;
 79: }

 81: PetscErrorCode AOApplicationToPetsc_Mapping(AO ao, PetscInt n, PetscInt *ia)
 82: {
 83:   AO_Mapping *aomap = (AO_Mapping*) ao->data;
 84:   PetscInt   *app   = aomap->app;
 85:   PetscInt   *petsc = aomap->petsc;
 86:   PetscInt   *perm  = aomap->appPerm;
 87:   PetscInt   N      = aomap->N;
 88:   PetscInt   low, high, mid=0;
 89:   PetscInt   idex;
 90:   PetscInt   i;

 92:   /* It would be possible to use a single bisection search, which
 93:      recursively divided the indices to be converted, and searched
 94:      partitions which contained an index. This would result in
 95:      better running times if indices are clustered.
 96:   */
 97:   for (i = 0; i < n; i++) {
 98:     idex = ia[i];
 99:     if (idex < 0) continue;
100:     /* Use bisection since the array is sorted */
101:     low  = 0;
102:     high = N - 1;
103:     while (low <= high) {
104:       mid = (low + high)/2;
105:       if (idex == app[mid]) break;
106:       else if (idex < app[mid]) high = mid - 1;
107:       else low = mid + 1;
108:     }
109:     if (low > high) ia[i] = -1;
110:     else            ia[i] = petsc[perm[mid]];
111:   }
112:   return 0;
113: }

115: static struct _AOOps AOps = {
116:   PetscDesignatedInitializer(view,AOView_Mapping),
117:   PetscDesignatedInitializer(destroy,AODestroy_Mapping),
118:   PetscDesignatedInitializer(petsctoapplication,AOPetscToApplication_Mapping),
119:   PetscDesignatedInitializer(applicationtopetsc,AOApplicationToPetsc_Mapping),
120: };

122: /*@C
123:   AOMappingHasApplicationIndex - Searches for the supplied application index.

125:   Input Parameters:
126: + ao       - The AOMapping
127: - index    - The application index

129:   Output Parameter:
130: . hasIndex - Flag is PETSC_TRUE if the index exists

132:   Level: intermediate

134: .seealso: AOMappingHasPetscIndex(), AOCreateMapping()
135: @*/
136: PetscErrorCode  AOMappingHasApplicationIndex(AO ao, PetscInt idex, PetscBool  *hasIndex)
137: {
138:   AO_Mapping *aomap;
139:   PetscInt   *app;
140:   PetscInt   low, high, mid;

144:   aomap = (AO_Mapping*) ao->data;
145:   app   = aomap->app;
146:   /* Use bisection since the array is sorted */
147:   low  = 0;
148:   high = aomap->N - 1;
149:   while (low <= high) {
150:     mid = (low + high)/2;
151:     if (idex == app[mid]) break;
152:     else if (idex < app[mid]) high = mid - 1;
153:     else low = mid + 1;
154:   }
155:   if (low > high) *hasIndex = PETSC_FALSE;
156:   else *hasIndex = PETSC_TRUE;
157:   return 0;
158: }

160: /*@
161:   AOMappingHasPetscIndex - Searches for the supplied petsc index.

163:   Input Parameters:
164: + ao       - The AOMapping
165: - index    - The petsc index

167:   Output Parameter:
168: . hasIndex - Flag is PETSC_TRUE if the index exists

170:   Level: intermediate

172: .seealso: AOMappingHasApplicationIndex(), AOCreateMapping()
173: @*/
174: PetscErrorCode  AOMappingHasPetscIndex(AO ao, PetscInt idex, PetscBool  *hasIndex)
175: {
176:   AO_Mapping *aomap;
177:   PetscInt   *petsc;
178:   PetscInt   low, high, mid;

182:   aomap = (AO_Mapping*) ao->data;
183:   petsc = aomap->petsc;
184:   /* Use bisection since the array is sorted */
185:   low  = 0;
186:   high = aomap->N - 1;
187:   while (low <= high) {
188:     mid = (low + high)/2;
189:     if (idex == petsc[mid]) break;
190:     else if (idex < petsc[mid]) high = mid - 1;
191:     else low = mid + 1;
192:   }
193:   if (low > high) *hasIndex = PETSC_FALSE;
194:   else *hasIndex = PETSC_TRUE;
195:   return 0;
196: }

198: /*@C
199:   AOCreateMapping - Creates a basic application mapping using two integer arrays.

201:   Input Parameters:
202: + comm    - MPI communicator that is to share AO
203: . napp    - size of integer arrays
204: . myapp   - integer array that defines an ordering
205: - mypetsc - integer array that defines another ordering (may be NULL to indicate the identity ordering)

207:   Output Parameter:
208: . aoout   - the new application mapping

210:   Options Database Key:
211: . -ao_view - call AOView() at the conclusion of AOCreateMapping()

213:   Level: beginner

215:     Notes:
216:     the arrays myapp and mypetsc need NOT contain the all the integers 0 to napp-1, that is there CAN be "holes"  in the indices.
217:        Use AOCreateBasic() or AOCreateBasicIS() if they do not have holes for better performance.

219: .seealso: AOCreateBasic(), AOCreateBasic(), AOCreateMappingIS(), AODestroy()
220: @*/
221: PetscErrorCode  AOCreateMapping(MPI_Comm comm,PetscInt napp,const PetscInt myapp[],const PetscInt mypetsc[],AO *aoout)
222: {
223:   AO             ao;
224:   AO_Mapping     *aomap;
225:   PetscInt       *allpetsc,  *allapp;
226:   PetscInt       *petscPerm, *appPerm;
227:   PetscInt       *petsc;
228:   PetscMPIInt    size, rank,*lens, *disp,nnapp;
229:   PetscInt       N, start;
230:   PetscInt       i;

233:   *aoout = NULL;
234:   AOInitializePackage();

236:   PetscHeaderCreate(ao, AO_CLASSID, "AO", "Application Ordering", "AO", comm, AODestroy, AOView);
237:   PetscNewLog(ao,&aomap);
238:   PetscMemcpy(ao->ops, &AOps, sizeof(AOps));
239:   ao->data = (void*) aomap;

241:   /* transmit all lengths to all processors */
242:   MPI_Comm_size(comm, &size);
243:   MPI_Comm_rank(comm, &rank);
244:   PetscMalloc2(size, &lens,size,&disp);
245:   nnapp = napp;
246:   MPI_Allgather(&nnapp, 1, MPI_INT, lens, 1, MPI_INT, comm);
247:   N     = 0;
248:   for (i = 0; i < size; i++) {
249:     disp[i] = N;
250:     N      += lens[i];
251:   }
252:   aomap->N = N;
253:   ao->N    = N;
254:   ao->n    = N;

256:   /* If mypetsc is 0 then use "natural" numbering */
257:   if (!mypetsc) {
258:     start = disp[rank];
259:     PetscMalloc1(napp+1, &petsc);
260:     for (i = 0; i < napp; i++) petsc[i] = start + i;
261:   } else {
262:     petsc = (PetscInt*)mypetsc;
263:   }

265:   /* get all indices on all processors */
266:   PetscMalloc4(N, &allapp,N,&appPerm,N,&allpetsc,N,&petscPerm);
267:   MPI_Allgatherv((void*)myapp, napp, MPIU_INT, allapp,   lens, disp, MPIU_INT, comm);
268:   MPI_Allgatherv((void*)petsc, napp, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm);
269:   PetscFree2(lens,disp);

271:   /* generate a list of application and PETSc node numbers */
272:   PetscMalloc4(N, &aomap->app,N,&aomap->appPerm,N,&aomap->petsc,N,&aomap->petscPerm);
273:   PetscLogObjectMemory((PetscObject)ao, 4*N * sizeof(PetscInt));
274:   for (i = 0; i < N; i++) {
275:     appPerm[i]   = i;
276:     petscPerm[i] = i;
277:   }
278:   PetscSortIntWithPermutation(N, allpetsc, petscPerm);
279:   PetscSortIntWithPermutation(N, allapp,   appPerm);
280:   /* Form sorted arrays of indices */
281:   for (i = 0; i < N; i++) {
282:     aomap->app[i]   = allapp[appPerm[i]];
283:     aomap->petsc[i] = allpetsc[petscPerm[i]];
284:   }
285:   /* Invert petscPerm[] into aomap->petscPerm[] */
286:   for (i = 0; i < N; i++) aomap->petscPerm[petscPerm[i]] = i;

288:   /* Form map between aomap->app[] and aomap->petsc[] */
289:   for (i = 0; i < N; i++) aomap->appPerm[i] = aomap->petscPerm[appPerm[i]];

291:   /* Invert appPerm[] into allapp[] */
292:   for (i = 0; i < N; i++) allapp[appPerm[i]] = i;

294:   /* Form map between aomap->petsc[] and aomap->app[] */
295:   for (i = 0; i < N; i++) aomap->petscPerm[i] = allapp[petscPerm[i]];

297:   if (PetscDefined(USE_DEBUG)) {
298:     /* Check that the permutations are complementary */
299:     for (i = 0; i < N; i++) {
301:     }
302:   }
303:   /* Cleanup */
304:   if (!mypetsc) {
305:     PetscFree(petsc);
306:   }
307:   PetscFree4(allapp,appPerm,allpetsc,petscPerm);

309:   AOViewFromOptions(ao,NULL,"-ao_view");

311:   *aoout = ao;
312:   return 0;
313: }

315: /*@
316:   AOCreateMappingIS - Creates a basic application ordering using two index sets.

318:   Input Parameters:
319: + comm    - MPI communicator that is to share AO
320: . isapp   - index set that defines an ordering
321: - ispetsc - index set that defines another ordering, maybe NULL for identity IS

323:   Output Parameter:
324: . aoout   - the new application ordering

326:   Options Database Key:
327: . -ao_view - call AOView() at the conclusion of AOCreateMappingIS()

329:   Level: beginner

331:     Notes:
332:     the index sets isapp and ispetsc need NOT contain the all the integers 0 to N-1, that is there CAN be "holes"  in the indices.
333:        Use AOCreateBasic() or AOCreateBasicIS() if they do not have holes for better performance.

335: .seealso: AOCreateBasic(), AOCreateMapping(), AODestroy()
336: @*/
337: PetscErrorCode  AOCreateMappingIS(IS isapp, IS ispetsc, AO *aoout)
338: {
339:   MPI_Comm       comm;
340:   const PetscInt *mypetsc, *myapp;
341:   PetscInt       napp, npetsc;

343:   PetscObjectGetComm((PetscObject) isapp, &comm);
344:   ISGetLocalSize(isapp, &napp);
345:   if (ispetsc) {
346:     ISGetLocalSize(ispetsc, &npetsc);
348:     ISGetIndices(ispetsc, &mypetsc);
349:   } else {
350:     mypetsc = NULL;
351:   }
352:   ISGetIndices(isapp, &myapp);

354:   AOCreateMapping(comm, napp, myapp, mypetsc, aoout);

356:   ISRestoreIndices(isapp, &myapp);
357:   if (ispetsc) {
358:     ISRestoreIndices(ispetsc, &mypetsc);
359:   }
360:   return 0;
361: }