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++) PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", i, aomap->app[i], aomap->petsc[aomap->appPerm[i]]);
41: }
42: return 0;
43: }
45: PetscErrorCode AOPetscToApplication_Mapping(AO ao, PetscInt n, PetscInt *ia)
46: {
47: AO_Mapping *aomap = (AO_Mapping *)ao->data;
48: PetscInt *app = aomap->app;
49: PetscInt *petsc = aomap->petsc;
50: PetscInt *perm = aomap->petscPerm;
51: PetscInt N = aomap->N;
52: PetscInt low, high, mid = 0;
53: PetscInt idex;
54: PetscInt i;
56: /* It would be possible to use a single bisection search, which
57: recursively divided the indices to be converted, and searched
58: partitions which contained an index. This would result in
59: better running times if indices are clustered.
60: */
61: for (i = 0; i < n; i++) {
62: idex = ia[i];
63: if (idex < 0) continue;
64: /* Use bisection since the array is sorted */
65: low = 0;
66: high = N - 1;
67: while (low <= high) {
68: mid = (low + high) / 2;
69: if (idex == petsc[mid]) break;
70: else if (idex < petsc[mid]) high = mid - 1;
71: else low = mid + 1;
72: }
73: if (low > high) ia[i] = -1;
74: else ia[i] = app[perm[mid]];
75: }
76: return 0;
77: }
79: PetscErrorCode AOApplicationToPetsc_Mapping(AO ao, PetscInt n, PetscInt *ia)
80: {
81: AO_Mapping *aomap = (AO_Mapping *)ao->data;
82: PetscInt *app = aomap->app;
83: PetscInt *petsc = aomap->petsc;
84: PetscInt *perm = aomap->appPerm;
85: PetscInt N = aomap->N;
86: PetscInt low, high, mid = 0;
87: PetscInt idex;
88: PetscInt i;
90: /* It would be possible to use a single bisection search, which
91: recursively divided the indices to be converted, and searched
92: partitions which contained an index. This would result in
93: better running times if indices are clustered.
94: */
95: for (i = 0; i < n; i++) {
96: idex = ia[i];
97: if (idex < 0) continue;
98: /* Use bisection since the array is sorted */
99: low = 0;
100: high = N - 1;
101: while (low <= high) {
102: mid = (low + high) / 2;
103: if (idex == app[mid]) break;
104: else if (idex < app[mid]) high = mid - 1;
105: else low = mid + 1;
106: }
107: if (low > high) ia[i] = -1;
108: else ia[i] = petsc[perm[mid]];
109: }
110: return 0;
111: }
113: static struct _AOOps AOps = {
114: PetscDesignatedInitializer(view, AOView_Mapping),
115: PetscDesignatedInitializer(destroy, AODestroy_Mapping),
116: PetscDesignatedInitializer(petsctoapplication, AOPetscToApplication_Mapping),
117: PetscDesignatedInitializer(applicationtopetsc, AOApplicationToPetsc_Mapping),
118: };
120: /*@C
121: AOMappingHasApplicationIndex - Checks if an `AO` has a requested application index.
123: Not Collective
125: Input Parameters:
126: + ao - The `AO`
127: - index - The application index
129: Output Parameter:
130: . hasIndex - Flag is `PETSC_TRUE` if the index exists
132: Level: intermediate
134: Developer Note:
135: The name of the function is wrong, it should be `AOHasApplicationIndex`
137: .seealso: [](sec_ao), `AOMappingHasPetscIndex()`, `AOCreateMapping()`, `AO`
138: @*/
139: PetscErrorCode AOMappingHasApplicationIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
140: {
141: AO_Mapping *aomap;
142: PetscInt *app;
143: PetscInt low, high, mid;
147: aomap = (AO_Mapping *)ao->data;
148: app = aomap->app;
149: /* Use bisection since the array is sorted */
150: low = 0;
151: high = aomap->N - 1;
152: while (low <= high) {
153: mid = (low + high) / 2;
154: if (idex == app[mid]) break;
155: else if (idex < app[mid]) high = mid - 1;
156: else low = mid + 1;
157: }
158: if (low > high) *hasIndex = PETSC_FALSE;
159: else *hasIndex = PETSC_TRUE;
160: return 0;
161: }
163: /*@
164: AOMappingHasPetscIndex - checks if an `AO` has a requested petsc index.
166: Not Collective
168: Input Parameters:
169: + ao - The `AO`
170: - index - The petsc index
172: Output Parameter:
173: . hasIndex - Flag is `PETSC_TRUE` if the index exists
175: Level: intermediate
177: Developer Note:
178: The name of the function is wrong, it should be `AOHasPetscIndex`
180: .seealso: [](sec_ao), `AOMappingHasApplicationIndex()`, `AOCreateMapping()`
181: @*/
182: PetscErrorCode AOMappingHasPetscIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
183: {
184: AO_Mapping *aomap;
185: PetscInt *petsc;
186: PetscInt low, high, mid;
190: aomap = (AO_Mapping *)ao->data;
191: petsc = aomap->petsc;
192: /* Use bisection since the array is sorted */
193: low = 0;
194: high = aomap->N - 1;
195: while (low <= high) {
196: mid = (low + high) / 2;
197: if (idex == petsc[mid]) break;
198: else if (idex < petsc[mid]) high = mid - 1;
199: else low = mid + 1;
200: }
201: if (low > high) *hasIndex = PETSC_FALSE;
202: else *hasIndex = PETSC_TRUE;
203: return 0;
204: }
206: /*@C
207: AOCreateMapping - Creates an application mapping using two integer arrays.
209: Input Parameters:
210: + comm - MPI communicator that is to share the `AO`
211: . napp - size of integer arrays
212: . myapp - integer array that defines an ordering
213: - mypetsc - integer array that defines another ordering (may be NULL to indicate the identity ordering)
215: Output Parameter:
216: . aoout - the new application mapping
218: Options Database Key:
219: . -ao_view - call `AOView()` at the conclusion of `AOCreateMapping()`
221: Level: beginner
223: Note:
224: 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.
225: Use `AOCreateBasic()` or `AOCreateBasicIS()` if they do not have holes for better performance.
227: .seealso: [](sec_ao), `AOCreateBasic()`, `AOCreateBasic()`, `AOCreateMappingIS()`, `AODestroy()`
228: @*/
229: PetscErrorCode AOCreateMapping(MPI_Comm comm, PetscInt napp, const PetscInt myapp[], const PetscInt mypetsc[], AO *aoout)
230: {
231: AO ao;
232: AO_Mapping *aomap;
233: PetscInt *allpetsc, *allapp;
234: PetscInt *petscPerm, *appPerm;
235: PetscInt *petsc;
236: PetscMPIInt size, rank, *lens, *disp, nnapp;
237: PetscInt N, start;
238: PetscInt i;
241: *aoout = NULL;
242: AOInitializePackage();
244: PetscHeaderCreate(ao, AO_CLASSID, "AO", "Application Ordering", "AO", comm, AODestroy, AOView);
245: PetscNew(&aomap);
246: PetscMemcpy(ao->ops, &AOps, sizeof(AOps));
247: ao->data = (void *)aomap;
249: /* transmit all lengths to all processors */
250: MPI_Comm_size(comm, &size);
251: MPI_Comm_rank(comm, &rank);
252: PetscMalloc2(size, &lens, size, &disp);
253: nnapp = napp;
254: MPI_Allgather(&nnapp, 1, MPI_INT, lens, 1, MPI_INT, comm);
255: N = 0;
256: for (i = 0; i < size; i++) {
257: disp[i] = N;
258: N += lens[i];
259: }
260: aomap->N = N;
261: ao->N = N;
262: ao->n = N;
264: /* If mypetsc is 0 then use "natural" numbering */
265: if (!mypetsc) {
266: start = disp[rank];
267: PetscMalloc1(napp + 1, &petsc);
268: for (i = 0; i < napp; i++) petsc[i] = start + i;
269: } else {
270: petsc = (PetscInt *)mypetsc;
271: }
273: /* get all indices on all processors */
274: PetscMalloc4(N, &allapp, N, &appPerm, N, &allpetsc, N, &petscPerm);
275: MPI_Allgatherv((void *)myapp, napp, MPIU_INT, allapp, lens, disp, MPIU_INT, comm);
276: MPI_Allgatherv((void *)petsc, napp, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm);
277: PetscFree2(lens, disp);
279: /* generate a list of application and PETSc node numbers */
280: PetscMalloc4(N, &aomap->app, N, &aomap->appPerm, N, &aomap->petsc, N, &aomap->petscPerm);
281: for (i = 0; i < N; i++) {
282: appPerm[i] = i;
283: petscPerm[i] = i;
284: }
285: PetscSortIntWithPermutation(N, allpetsc, petscPerm);
286: PetscSortIntWithPermutation(N, allapp, appPerm);
287: /* Form sorted arrays of indices */
288: for (i = 0; i < N; i++) {
289: aomap->app[i] = allapp[appPerm[i]];
290: aomap->petsc[i] = allpetsc[petscPerm[i]];
291: }
292: /* Invert petscPerm[] into aomap->petscPerm[] */
293: for (i = 0; i < N; i++) aomap->petscPerm[petscPerm[i]] = i;
295: /* Form map between aomap->app[] and aomap->petsc[] */
296: for (i = 0; i < N; i++) aomap->appPerm[i] = aomap->petscPerm[appPerm[i]];
298: /* Invert appPerm[] into allapp[] */
299: for (i = 0; i < N; i++) allapp[appPerm[i]] = i;
301: /* Form map between aomap->petsc[] and aomap->app[] */
302: for (i = 0; i < N; i++) aomap->petscPerm[i] = allapp[petscPerm[i]];
304: if (PetscDefined(USE_DEBUG)) {
305: /* Check that the permutations are complementary */
307: }
308: /* Cleanup */
309: if (!mypetsc) PetscFree(petsc);
310: PetscFree4(allapp, appPerm, allpetsc, petscPerm);
312: AOViewFromOptions(ao, NULL, "-ao_view");
314: *aoout = ao;
315: return 0;
316: }
318: /*@
319: AOCreateMappingIS - Creates an application mapping using two index sets.
321: Input Parameters:
322: + comm - MPI communicator that is to share `AO`
323: . isapp - index set that defines an ordering
324: - ispetsc - index set that defines another ordering, maybe NULL for identity `IS`
326: Output Parameter:
327: . aoout - the new application ordering
329: Options Database Key:
330: . -ao_view - call `AOView()` at the conclusion of `AOCreateMappingIS()`
332: Level: beginner
334: Note:
335: 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.
336: Use `AOCreateBasic()` or `AOCreateBasicIS()` if they do not have holes for better performance.
338: .seealso: [](sec_ao), [](sec_scatter), `AOCreateBasic()`, `AOCreateMapping()`, `AODestroy()`
339: @*/
340: PetscErrorCode AOCreateMappingIS(IS isapp, IS ispetsc, AO *aoout)
341: {
342: MPI_Comm comm;
343: const PetscInt *mypetsc, *myapp;
344: PetscInt napp, npetsc;
346: PetscObjectGetComm((PetscObject)isapp, &comm);
347: ISGetLocalSize(isapp, &napp);
348: if (ispetsc) {
349: ISGetLocalSize(ispetsc, &npetsc);
351: ISGetIndices(ispetsc, &mypetsc);
352: } else {
353: mypetsc = NULL;
354: }
355: ISGetIndices(isapp, &myapp);
357: AOCreateMapping(comm, napp, myapp, mypetsc, aoout);
359: ISRestoreIndices(isapp, &myapp);
360: if (ispetsc) ISRestoreIndices(ispetsc, &mypetsc);
361: return 0;
362: }