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: }