Actual source code: aomapping.c

  1: #define PETSCDM_DLL

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

 9:  #include src/dm/ao/aoimpl.h
 10:  #include petscsys.h

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

 22: PetscErrorCode AODestroy_Mapping(AO ao)
 23: {
 24:   AO_Mapping     *aomap = (AO_Mapping *) ao->data;

 28:   PetscFree(aomap->app);
 29:   PetscFree(ao->data);
 30:   return(0);
 31: }

 35: PetscErrorCode AOView_Mapping(AO ao, PetscViewer viewer)
 36: {
 37:   AO_Mapping     *aomap = (AO_Mapping *) ao->data;
 38:   PetscMPIInt    rank;
 39:   PetscInt       i;
 40:   PetscTruth     iascii;

 44:   MPI_Comm_rank(ao->comm, &rank);
 45:   if (rank) return(0);

 47:   if (!viewer) {
 48:     viewer = PETSC_VIEWER_STDOUT_SELF;
 49:   }

 51:   PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_ASCII, &iascii);
 52:   if (iascii) {
 53:     PetscViewerASCIIPrintf(viewer, "Number of elements in ordering %D\n", aomap->N);
 54:     PetscViewerASCIIPrintf(viewer, "   App.   PETSc\n");
 55:     for(i = 0; i < aomap->N; i++) {
 56:       PetscViewerASCIIPrintf(viewer, "%D   %D    %D\n", i, aomap->app[i], aomap->petsc[aomap->appPerm[i]]);
 57:     }
 58:   }
 59:   return(0);
 60: }

 64: PetscErrorCode AOPetscToApplication_Mapping(AO ao, PetscInt n, PetscInt *ia)
 65: {
 66:   AO_Mapping *aomap = (AO_Mapping *) ao->data;
 67:   PetscInt   *app   = aomap->app;
 68:   PetscInt   *petsc = aomap->petsc;
 69:   PetscInt   *perm  = aomap->petscPerm;
 70:   PetscInt   N     = aomap->N;
 71:   PetscInt   low, high, mid=0;
 72:   PetscInt   idex;
 73:   PetscInt   i;

 75:   /* It would be possible to use a single bisection search, which
 76:      recursively divided the indices to be converted, and searched
 77:      partitions which contained an index. This would result in
 78:      better running times if indices are clustered.
 79:   */
 81:   for(i = 0; i < n; i++) {
 82:     idex = ia[i];
 83:     if (idex < 0) continue;
 84:     /* Use bisection since the array is sorted */
 85:     low  = 0;
 86:     high = N - 1;
 87:     while (low <= high) {
 88:       mid = (low + high)/2;
 89:       if (idex == petsc[mid]) {
 90:         break;
 91:       } else if (idex < petsc[mid]) {
 92:         high = mid - 1;
 93:       } else {
 94:         low  = mid + 1;
 95:       }
 96:     }
 97:     if (low > high) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Invalid input index %D", idex);
 98:     ia[i] = app[perm[mid]];
 99:   }
100:   return(0);
101: }

105: PetscErrorCode AOApplicationToPetsc_Mapping(AO ao, PetscInt n, PetscInt *ia)
106: {
107:   AO_Mapping *aomap = (AO_Mapping *) ao->data;
108:   PetscInt   *app   = aomap->app;
109:   PetscInt   *petsc = aomap->petsc;
110:   PetscInt   *perm  = aomap->appPerm;
111:   PetscInt   N     = aomap->N;
112:   PetscInt   low, high, mid=0;
113:   PetscInt   idex;
114:   PetscInt   i;

116:   /* It would be possible to use a single bisection search, which
117:      recursively divided the indices to be converted, and searched
118:      partitions which contained an index. This would result in
119:      better running times if indices are clustered.
120:   */
122:   for(i = 0; i < n; i++) {
123:     idex = ia[i];
124:     if (idex < 0) continue;
125:     /* Use bisection since the array is sorted */
126:     low  = 0;
127:     high = N - 1;
128:     while (low <= high) {
129:       mid = (low + high)/2;
130:       if (idex == app[mid]) {
131:         break;
132:       } else if (idex < app[mid]) {
133:         high = mid - 1;
134:       } else {
135:         low  = mid + 1;
136:       }
137:     }
138:     if (low > high) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Invalid input index %D", idex);
139:     ia[i] = petsc[perm[mid]];
140:   }
141:   return(0);
142: }

144: static struct _AOOps AOps = {AOView_Mapping,
145:                              AODestroy_Mapping,
146:                              AOPetscToApplication_Mapping,
147:                              AOApplicationToPetsc_Mapping,
148:                              PETSC_NULL,
149:                              PETSC_NULL,
150:                              PETSC_NULL,
151:                              PETSC_NULL};

155: /*@C
156:   AOMappingHasApplicationIndex - Searches for the supplied application index.

158:   Input Parameters:
159: + ao       - The AOMapping
160: - index    - The application index

162:   Output Parameter:
163: . hasIndex - Flag is PETSC_TRUE if the index exists

165:   Level: intermediate

167: .keywords: AO, index
168: .seealso: AOMappingHasPetscIndex(), AOCreateMapping()
169: @*/
170: PetscErrorCode  AOMappingHasApplicationIndex(AO ao, PetscInt idex, PetscTruth *hasIndex)
171: {
172:   AO_Mapping *aomap;
173:   PetscInt   *app;
174:   PetscInt   low, high, mid;

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

204: /*@C
205:   AOMappingHasPetscIndex - Searches for the supplied petsc index.

207:   Input Parameters:
208: + ao       - The AOMapping
209: - index    - The petsc index

211:   Output Parameter:
212: . hasIndex - Flag is PETSC_TRUE if the index exists

214:   Level: intermediate

216: .keywords: AO, index
217: .seealso: AOMappingHasApplicationIndex(), AOCreateMapping()
218: @*/
219: PetscErrorCode  AOMappingHasPetscIndex(AO ao, PetscInt idex, PetscTruth *hasIndex)
220: {
221:   AO_Mapping *aomap;
222:   PetscInt   *petsc;
223:   PetscInt   low, high, mid;

228:   aomap = (AO_Mapping *) ao->data;
229:   petsc = aomap->petsc;
230:   /* Use bisection since the array is sorted */
231:   low  = 0;
232:   high = aomap->N - 1;
233:   while (low <= high) {
234:     mid = (low + high)/2;
235:     if (idex == petsc[mid]) {
236:       break;
237:     } else if (idex < petsc[mid]) {
238:       high = mid - 1;
239:     } else {
240:       low  = mid + 1;
241:     }
242:   }
243:   if (low > high) {
244:     *hasIndex = PETSC_FALSE;
245:   } else {
246:     *hasIndex = PETSC_TRUE;
247:   }
248:   return(0);
249: }

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

256:   Input Parameters:
257: + comm    - MPI communicator that is to share AO
258: . napp    - size of integer arrays
259: . myapp   - integer array that defines an ordering
260: - mypetsc - integer array that defines another ordering

262:   Output Parameter:
263: . aoout   - the new application mapping

265:   Options Database Key:
266: $ -ao_view : call AOView() at the conclusion of AOCreateMapping()

268:   Level: beginner

270: .keywords: AO, create
271: .seealso: AOCreateDebug(), AOCreateBasic(), AOCreateMappingIS(), AODestroy()
272: @*/
273: PetscErrorCode  AOCreateMapping(MPI_Comm comm,PetscInt napp,const PetscInt myapp[],const PetscInt mypetsc[],AO *aoout)
274: {
275:   AO             ao;
276:   AO_Mapping     *aomap;
277:   PetscInt       *allpetsc,  *allapp;
278:   PetscInt       *petscPerm, *appPerm;
279:   PetscInt       *petsc;
280:   PetscMPIInt    size, rank,*lens, *disp,nnapp;
281:   PetscInt       N, start;
282:   PetscInt       i;
283:   PetscTruth     opt;

288:   *aoout = 0;
289: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
290:   DMInitializePackage(PETSC_NULL);
291: #endif

293:   PetscHeaderCreate(ao, _p_AO, struct _AOOps, AO_COOKIE, AO_MAPPING, "AO", comm, AODestroy, AOView);
294:   PetscNew(AO_Mapping, &aomap);
295:   PetscLogObjectMemory(ao, sizeof(struct _p_AO) + sizeof(AO_Mapping));
296:   PetscMemcpy(ao->ops, &AOps, sizeof(AOps));
297:   ao->data = (void*) aomap;

299:   /* transmit all lengths to all processors */
300:   MPI_Comm_size(comm, &size);
301:   MPI_Comm_rank(comm, &rank);
302:   PetscMalloc(2*size * sizeof(PetscMPIInt), &lens);
303:   disp  = lens + size;
304:   nnapp = napp;
305:   MPI_Allgather(&nnapp, 1, MPI_INT, lens, 1, MPI_INT, comm);
306:   N    = 0;
307:   for(i = 0; i < size; i++) {
308:     disp[i] = N;
309:     N += lens[i];
310:   }
311:   aomap->N = N;
312:   ao->N    = N;
313:   ao->n    = N;

315:   /* If mypetsc is 0 then use "natural" numbering */
316:   if (!mypetsc) {
317:     start = disp[rank];
318:     PetscMalloc((napp+1) * sizeof(PetscInt), &petsc);
319:     for(i = 0; i < napp; i++) {
320:       petsc[i] = start + i;
321:     }
322:   } else {
323:     petsc = (PetscInt*)mypetsc;
324:   }

326:   /* get all indices on all processors */
327:   PetscMalloc(N*4 * sizeof(PetscInt), &allapp);
328:   appPerm   = allapp   + N;
329:   allpetsc  = appPerm  + N;
330:   petscPerm = allpetsc + N;
331:   MPI_Allgatherv((void*)myapp, napp, MPIU_INT, allapp,   lens, disp, MPIU_INT, comm);
332:   MPI_Allgatherv((void*)petsc, napp, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm);
333:   PetscFree(lens);

335:   /* generate a list of application and PETSc node numbers */
336:   PetscMalloc(N*4 * sizeof(PetscInt), &aomap->app);
337:   PetscLogObjectMemory(ao, 4*N * sizeof(PetscInt));
338:   aomap->appPerm   = aomap->app     + N;
339:   aomap->petsc     = aomap->appPerm + N;
340:   aomap->petscPerm = aomap->petsc   + N;
341:   for(i = 0; i < N; i++) {
342:     appPerm[i]   = i;
343:     petscPerm[i] = i;
344:   }
345:   PetscSortIntWithPermutation(N, allpetsc, petscPerm);
346:   PetscSortIntWithPermutation(N, allapp,   appPerm);
347:   /* Form sorted arrays of indices */
348:   for(i = 0; i < N; i++) {
349:     aomap->app[i]   = allapp[appPerm[i]];
350:     aomap->petsc[i] = allpetsc[petscPerm[i]];
351:   }
352:   /* Invert petscPerm[] into aomap->petscPerm[] */
353:   for(i = 0; i < N; i++) {
354:     aomap->petscPerm[petscPerm[i]] = i;
355:   }
356:   /* Form map between aomap->app[] and aomap->petsc[] */
357:   for(i = 0; i < N; i++) {
358:     aomap->appPerm[i] = aomap->petscPerm[appPerm[i]];
359:   }
360:   /* Invert appPerm[] into allapp[] */
361:   for(i = 0; i < N; i++) {
362:     allapp[appPerm[i]] = i;
363:   }
364:   /* Form map between aomap->petsc[] and aomap->app[] */
365:   for(i = 0; i < N; i++) {
366:     aomap->petscPerm[i] = allapp[petscPerm[i]];
367:   }
368: #ifdef PETSC_USE_DEBUG
369:   /* Check that the permutations are complementary */
370:   for(i = 0; i < N; i++) {
371:     if (i != aomap->appPerm[aomap->petscPerm[i]])
372:       SETERRQ(PETSC_ERR_PLIB, "Invalid ordering");
373:   }
374: #endif
375:   /* Cleanup */
376:   if (!mypetsc) {
377:     PetscFree(petsc);
378:   }
379:   PetscFree(allapp);

381:   PetscOptionsHasName(PETSC_NULL, "-ao_view", &opt);
382:   if (opt) {
383:     AOView(ao, PETSC_VIEWER_STDOUT_SELF);
384:   }

386:   *aoout = ao;
387:   return(0);
388: }

392: /*@C
393:   AOCreateMappingIS - Creates a basic application ordering using two index sets.

395:   Input Parameters:
396: + comm    - MPI communicator that is to share AO
397: . isapp   - index set that defines an ordering
398: - ispetsc - index set that defines another ordering

400:   Output Parameter:
401: . aoout   - the new application ordering

403:   Options Database Key:
404: $ -ao_view : call AOView() at the conclusion of AOCreateMappingIS()

406:   Level: beginner

408: .keywords: AO, create
409: .seealso: AOCreateBasic(), AOCreateMapping(), AODestroy()
410: @*/
411: PetscErrorCode  AOCreateMappingIS(IS isapp, IS ispetsc, AO *aoout)
412: {
413:   MPI_Comm       comm;
414:   PetscInt       *mypetsc, *myapp;
415:   PetscInt       napp, npetsc;

419:   PetscObjectGetComm((PetscObject) isapp, &comm);
420:   ISGetLocalSize(isapp, &napp);
421:   if (ispetsc) {
422:     ISGetLocalSize(ispetsc, &npetsc);
423:     if (napp != npetsc) SETERRQ(PETSC_ERR_ARG_SIZ, "Local IS lengths must match");
424:     ISGetIndices(ispetsc, &mypetsc);
425:   } else {
426:     mypetsc = NULL;
427:   }
428:   ISGetIndices(isapp, &myapp);

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

432:   ISRestoreIndices(isapp, &myapp);
433:   if (ispetsc) {
434:     ISRestoreIndices(ispetsc, &mypetsc);
435:   }
436:   return(0);
437: }