Actual source code: aomapping.c

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

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

 11:  #include src/dm/ao/aoimpl.h
 12:  #include petscsys.h

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

 22: #undef __FUNCT__  
 24: int AODestroy_Mapping(AO ao)
 25: {
 26:   AO_Mapping *aomap = (AO_Mapping *) ao->data;
 27:   int         ierr;

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

 35: #undef __FUNCT__  
 37: int AOView_Mapping(AO ao, PetscViewer viewer)
 38: {
 39:   AO_Mapping *aomap = (AO_Mapping *) ao->data;
 40:   int         rank, i;
 41:   PetscTruth  isascii;
 42:   int         ierr;

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

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

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

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

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

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

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

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

154: EXTERN_C_BEGIN
155: #undef __FUNCT__  
157: int AOSerialize_Mapping(MPI_Comm comm, AO *ao, PetscViewer viewer, PetscTruth store)
158: {
159:   AO          a;
160:   AO_Mapping *m;
161:   int         fd;
162:   int         ierr;

165:   PetscViewerBinaryGetDescriptor(viewer, &fd);
166:   if (store) {
167:     a = *ao;
168:     m = (AO_Mapping *) a->data;
169:     PetscBinaryWrite(fd, &a->N,         1,    PETSC_INT, 0);
170:     PetscBinaryWrite(fd, &a->n,         1,    PETSC_INT, 0);
171:     PetscBinaryWrite(fd, &m->N,         1,    PETSC_INT, 0);
172:     PetscBinaryWrite(fd,  m->app,       m->N, PETSC_INT, 0);
173:     PetscBinaryWrite(fd,  m->appPerm,   m->N, PETSC_INT, 0);
174:     PetscBinaryWrite(fd,  m->petsc,     m->N, PETSC_INT, 0);
175:     PetscBinaryWrite(fd,  m->petscPerm, m->N, PETSC_INT, 0);
176:   } else {
177:     PetscHeaderCreate(a, _p_AO, struct _AOOps, AO_COOKIE, AO_MAPPING, "AO", comm, AODestroy, AOView);
178:     PetscLogObjectCreate(a);
179:     PetscLogObjectMemory(a, sizeof(AO_Mapping) + sizeof(struct _p_AO));
180:     PetscNew(AO_Mapping, &m);
181:     a->data = (void *) m;
182:     PetscMemcpy(a->ops, &AOps, sizeof(AOps));

184:     PetscBinaryRead(fd, &a->N,     1,    PETSC_INT);
185:     PetscBinaryRead(fd, &a->n,     1,    PETSC_INT);
186:     PetscBinaryRead(fd, &m->N,     1,    PETSC_INT);
187:     PetscMalloc(m->N*4 * sizeof(int), &m->app);
188:     m->appPerm   = m->app     + m->N;
189:     m->petsc     = m->appPerm + m->N;
190:     m->petscPerm = m->petsc   + m->N;
191:     PetscLogObjectMemory(a, m->N*4 * sizeof(int));
192:     PetscBinaryRead(fd,  m->app,       m->N, PETSC_INT);
193:     PetscBinaryRead(fd,  m->appPerm,   m->N, PETSC_INT);
194:     PetscBinaryRead(fd,  m->petsc,     m->N, PETSC_INT);
195:     PetscBinaryRead(fd,  m->petscPerm, m->N, PETSC_INT);
196:     *ao = a;
197:   }

199:   return(0);
200: }
201: EXTERN_C_END

203: #undef __FUNCT__  
205: /*@C
206:   AOMappingHasApplicationIndex - Searches for the supplied application index.

208:   Input Parameters:
209: + ao       - The AOMapping
210: - index    - The application index

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

215:   Level: intermediate

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

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

252: #undef __FUNCT__  
254: /*@C
255:   AOMappingHasPetscIndex - Searches for the supplied petsc index.

257:   Input Parameters:
258: + ao       - The AOMapping
259: - index    - The petsc index

261:   Output Parameter:
262: . hasIndex - Flag is PETSC_TRUE if the index exists

264:   Level: intermediate

266: .keywords: AO, index
267: .seealso: AOMappingHasApplicationIndex(), AOCreateMapping()
268: @*/
269: int AOMappingHasPetscIndex(AO ao, int idex, PetscTruth *hasIndex)
270: {
271:   AO_Mapping *aomap;
272:   int        *petsc;
273:   int         low, high, mid;

278:   aomap = (AO_Mapping *) ao->data;
279:   petsc = aomap->petsc;
280:   /* Use bisection since the array is sorted */
281:   low  = 0;
282:   high = aomap->N - 1;
283:   while (low <= high) {
284:     mid = (low + high)/2;
285:     if (idex == petsc[mid]) {
286:       break;
287:     } else if (idex < petsc[mid]) {
288:       high = mid - 1;
289:     } else {
290:       low  = mid + 1;
291:     }
292:   }
293:   if (low > high) {
294:     *hasIndex = PETSC_FALSE;
295:   } else {
296:     *hasIndex = PETSC_TRUE;
297:   }
298:   return(0);
299: }

301: #undef __FUNCT__  
303: /*@C
304:   AOCreateMapping - Creates a basic application mapping using two integer arrays.

306:   Input Parameters:
307: + comm    - MPI communicator that is to share AO
308: . napp    - size of integer arrays
309: . myapp   - integer array that defines an ordering
310: - mypetsc - integer array that defines another ordering

312:   Output Parameter:
313: . aoout   - the new application mapping

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

318:   Level: beginner

320: .keywords: AO, create
321: .seealso: AOCreateDebug(), AOCreateBasic(), AOCreateMappingIS(), AODestroy()
322: @*/
323: int AOCreateMapping(MPI_Comm comm, int napp, int *myapp, int *mypetsc, AO *aoout)
324: {
325:   AO          ao;
326:   AO_Mapping *aomap;
327:   int        *lens, *disp;
328:   int        *allpetsc,  *allapp;
329:   int        *petscPerm, *appPerm;
330:   int        *petsc;
331:   int         size, rank;
332:   int         N, start;
333:   int         i;
334:   PetscTruth  opt;
335:   int         ierr;

339:   *aoout = 0;
340: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
341:   DMInitializePackage(PETSC_NULL);
342: #endif

344:   PetscHeaderCreate(ao, _p_AO, struct _AOOps, AO_COOKIE, AO_MAPPING, "AO", comm, AODestroy, AOView);
345:   PetscLogObjectCreate(ao);
346:   PetscNew(AO_Mapping, &aomap);
347:   PetscLogObjectMemory(ao, sizeof(struct _p_AO) + sizeof(AO_Mapping));
348:   PetscMemcpy(ao->ops, &AOps, sizeof(AOps));
349:   PetscStrallocpy(AO_SER_MAPPING_BINARY, &ao->serialize_name);
350:   ao->data = (void *) aomap;

352:   /* transmit all lengths to all processors */
353:   MPI_Comm_size(comm, &size);
354:   MPI_Comm_rank(comm, &rank);
355:   PetscMalloc(2*size * sizeof(int), &lens);
356:   disp = lens + size;
357:   MPI_Allgather(&napp, 1, MPI_INT, lens, 1, MPI_INT, comm);
358:   N    = 0;
359:   for(i = 0; i < size; i++) {
360:     disp[i] = N;
361:     N += lens[i];
362:   }
363:   aomap->N = N;
364:   ao->N    = N;
365:   ao->n    = N;

367:   /* If mypetsc is 0 then use "natural" numbering */
368:   if (!mypetsc) {
369:     start = disp[rank];
370:     ierr  = PetscMalloc((napp+1) * sizeof(int), &petsc);
371:     for(i = 0; i < napp; i++) {
372:       petsc[i] = start + i;
373:     }
374:   } else {
375:     petsc = mypetsc;
376:   }

378:   /* get all indices on all processors */
379:   PetscMalloc(N*4 * sizeof(int), &allapp);
380:   appPerm   = allapp   + N;
381:   allpetsc  = appPerm  + N;
382:   petscPerm = allpetsc + N;
383:   MPI_Allgatherv(myapp,   napp, MPI_INT, allapp,   lens, disp, MPI_INT, comm);
384:   MPI_Allgatherv(mypetsc, napp, MPI_INT, allpetsc, lens, disp, MPI_INT, comm);
385:   PetscFree(lens);

387:   /* generate a list of application and PETSc node numbers */
388:   PetscMalloc(N*4 * sizeof(int), &aomap->app);
389:   PetscLogObjectMemory(ao, 4*N * sizeof(int));
390:   aomap->appPerm   = aomap->app     + N;
391:   aomap->petsc     = aomap->appPerm + N;
392:   aomap->petscPerm = aomap->petsc   + N;
393:   for(i = 0; i < N; i++) {
394:     appPerm[i]   = i;
395:     petscPerm[i] = i;
396:   }
397:   PetscSortIntWithPermutation(N, allpetsc, petscPerm);
398:   PetscSortIntWithPermutation(N, allapp,   appPerm);
399:   /* Form sorted arrays of indices */
400:   for(i = 0; i < N; i++) {
401:     aomap->app[i]   = allapp[appPerm[i]];
402:     aomap->petsc[i] = allpetsc[petscPerm[i]];
403:   }
404:   /* Invert petscPerm[] into aomap->petscPerm[] */
405:   for(i = 0; i < N; i++) {
406:     aomap->petscPerm[petscPerm[i]] = i;
407:   }
408:   /* Form map between aomap->app[] and aomap->petsc[] */
409:   for(i = 0; i < N; i++) {
410:     aomap->appPerm[i] = aomap->petscPerm[appPerm[i]];
411:   }
412:   /* Invert appPerm[] into allapp[] */
413:   for(i = 0; i < N; i++) {
414:     allapp[appPerm[i]] = i;
415:   }
416:   /* Form map between aomap->petsc[] and aomap->app[] */
417:   for(i = 0; i < N; i++) {
418:     aomap->petscPerm[i] = allapp[petscPerm[i]];
419:   }
420: #ifdef PETSC_USE_BOPT_g
421:   /* Check that the permutations are complementary */
422:   for(i = 0; i < N; i++) {
423:     if (i != aomap->appPerm[aomap->petscPerm[i]])
424:       SETERRQ(PETSC_ERR_PLIB, "Invalid ordering");
425:   }
426: #endif
427:   /* Cleanup */
428:   if (!mypetsc) {
429:     PetscFree(petsc);
430:   }
431:   PetscFree(allapp);

433:   PetscOptionsHasName(PETSC_NULL, "-ao_view", &opt);
434:   if (opt == PETSC_TRUE) {
435:     AOView(ao, PETSC_VIEWER_STDOUT_SELF);
436:   }

438:   *aoout = ao;
439:   return(0);
440: }

442: #undef __FUNCT__  
444: /*@C
445:   AOCreateMappingIS - Creates a basic application ordering using two index sets.

447:   Input Parameters:
448: + comm    - MPI communicator that is to share AO
449: . isapp   - index set that defines an ordering
450: - ispetsc - index set that defines another ordering

452:   Output Parameter:
453: . aoout   - the new application ordering

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

458:   Level: beginner

460: .keywords: AO, create
461: .seealso: AOCreateBasic(), AOCreateMapping(), AODestroy()
462: @*/
463: int AOCreateMappingIS(IS isapp, IS ispetsc, AO *aoout)
464: {
465:   MPI_Comm comm;
466:   int     *mypetsc, *myapp;
467:   int      napp, npetsc;
468:   int      ierr;

471:   PetscObjectGetComm((PetscObject) isapp, &comm);
472:   ISGetSize(isapp, &napp);
473:   if (ispetsc) {
474:     ISGetSize(ispetsc, &npetsc);
475:     if (napp != npetsc) SETERRQ(PETSC_ERR_ARG_SIZ, "Local IS lengths must match");
476:     ISGetIndices(ispetsc, &mypetsc);
477:   }
478:   ISGetIndices(isapp, &myapp);

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

482:   ISRestoreIndices(isapp, &myapp);
483:   if (ispetsc) {
484:     ISRestoreIndices(ispetsc, &mypetsc);
485:   }
486:   return(0);
487: }