Actual source code: aomapping.c
1: /*
2: These AO application ordering routines do not require that the input
3: be a permutation, but merely a 1-1 mapping. This implementation still
4: keeps the entire ordering on each processor.
5: */
7: #include <../src/vec/is/ao/aoimpl.h>
9: typedef struct {
10: PetscInt N;
11: PetscInt *app; /* app[i] is the partner for petsc[appPerm[i]] */
12: PetscInt *appPerm;
13: PetscInt *petsc; /* petsc[j] is the partner for app[petscPerm[j]] */
14: PetscInt *petscPerm;
15: } AO_Mapping;
17: static PetscErrorCode AODestroy_Mapping(AO ao)
18: {
19: AO_Mapping *aomap = (AO_Mapping *)ao->data;
21: PetscFunctionBegin;
22: PetscCall(PetscFree4(aomap->app, aomap->appPerm, aomap->petsc, aomap->petscPerm));
23: PetscCall(PetscFree(aomap));
24: PetscFunctionReturn(PETSC_SUCCESS);
25: }
27: static 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: PetscFunctionBegin;
35: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ao), &rank));
36: if (rank) PetscFunctionReturn(PETSC_SUCCESS);
37: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
38: if (iascii) {
39: PetscCall(PetscViewerASCIIPrintf(viewer, "Number of elements in ordering %" PetscInt_FMT "\n", aomap->N));
40: PetscCall(PetscViewerASCIIPrintf(viewer, " App. PETSc\n"));
41: for (i = 0; i < aomap->N; i++) PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", i, aomap->app[i], aomap->petsc[aomap->appPerm[i]]));
42: }
43: PetscFunctionReturn(PETSC_SUCCESS);
44: }
46: static PetscErrorCode AOPetscToApplication_Mapping(AO ao, PetscInt n, PetscInt *ia)
47: {
48: AO_Mapping *aomap = (AO_Mapping *)ao->data;
49: PetscInt *app = aomap->app;
50: PetscInt *petsc = aomap->petsc;
51: PetscInt *perm = aomap->petscPerm;
52: PetscInt N = aomap->N;
53: PetscInt low, high, mid = 0;
54: PetscInt idex;
55: PetscInt i;
57: /* It would be possible to use a single bisection search, which
58: recursively divided the indices to be converted, and searched
59: partitions which contained an index. This would result in
60: better running times if indices are clustered.
61: */
62: PetscFunctionBegin;
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: PetscFunctionReturn(PETSC_SUCCESS);
79: }
81: static 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: PetscFunctionBegin;
98: for (i = 0; i < n; i++) {
99: idex = ia[i];
100: if (idex < 0) continue;
101: /* Use bisection since the array is sorted */
102: low = 0;
103: high = N - 1;
104: while (low <= high) {
105: mid = (low + high) / 2;
106: if (idex == app[mid]) break;
107: else if (idex < app[mid]) high = mid - 1;
108: else low = mid + 1;
109: }
110: if (low > high) ia[i] = -1;
111: else ia[i] = petsc[perm[mid]];
112: }
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: static const struct _AOOps AOps = {
117: PetscDesignatedInitializer(view, AOView_Mapping),
118: PetscDesignatedInitializer(destroy, AODestroy_Mapping),
119: PetscDesignatedInitializer(petsctoapplication, AOPetscToApplication_Mapping),
120: PetscDesignatedInitializer(applicationtopetsc, AOApplicationToPetsc_Mapping),
121: };
123: /*@C
124: AOMappingHasApplicationIndex - Checks if an `AO` has a requested application index.
126: Not Collective
128: Input Parameters:
129: + ao - The `AO`
130: - idex - The application index
132: Output Parameter:
133: . hasIndex - Flag is `PETSC_TRUE` if the index exists
135: Level: intermediate
137: Developer Notes:
138: The name of the function is wrong, it should be `AOHasApplicationIndex`
140: .seealso: [](sec_ao), `AOMappingHasPetscIndex()`, `AOCreateMapping()`, `AO`
141: @*/
142: PetscErrorCode AOMappingHasApplicationIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
143: {
144: AO_Mapping *aomap;
145: PetscInt *app;
146: PetscInt low, high, mid;
148: PetscFunctionBegin;
150: PetscAssertPointer(hasIndex, 3);
151: aomap = (AO_Mapping *)ao->data;
152: app = aomap->app;
153: /* Use bisection since the array is sorted */
154: low = 0;
155: high = aomap->N - 1;
156: while (low <= high) {
157: mid = (low + high) / 2;
158: if (idex == app[mid]) break;
159: else if (idex < app[mid]) high = mid - 1;
160: else low = mid + 1;
161: }
162: if (low > high) *hasIndex = PETSC_FALSE;
163: else *hasIndex = PETSC_TRUE;
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }
167: /*@
168: AOMappingHasPetscIndex - checks if an `AO` has a requested petsc index.
170: Not Collective
172: Input Parameters:
173: + ao - The `AO`
174: - idex - The petsc index
176: Output Parameter:
177: . hasIndex - Flag is `PETSC_TRUE` if the index exists
179: Level: intermediate
181: Developer Notes:
182: The name of the function is wrong, it should be `AOHasPetscIndex`
184: .seealso: [](sec_ao), `AOMappingHasApplicationIndex()`, `AOCreateMapping()`
185: @*/
186: PetscErrorCode AOMappingHasPetscIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
187: {
188: AO_Mapping *aomap;
189: PetscInt *petsc;
190: PetscInt low, high, mid;
192: PetscFunctionBegin;
194: PetscAssertPointer(hasIndex, 3);
195: aomap = (AO_Mapping *)ao->data;
196: petsc = aomap->petsc;
197: /* Use bisection since the array is sorted */
198: low = 0;
199: high = aomap->N - 1;
200: while (low <= high) {
201: mid = (low + high) / 2;
202: if (idex == petsc[mid]) break;
203: else if (idex < petsc[mid]) high = mid - 1;
204: else low = mid + 1;
205: }
206: if (low > high) *hasIndex = PETSC_FALSE;
207: else *hasIndex = PETSC_TRUE;
208: PetscFunctionReturn(PETSC_SUCCESS);
209: }
211: /*@C
212: AOCreateMapping - Creates an application mapping using two integer arrays.
214: Input Parameters:
215: + comm - MPI communicator that is to share the `AO`
216: . napp - size of integer arrays
217: . myapp - integer array that defines an ordering
218: - mypetsc - integer array that defines another ordering (may be `NULL` to indicate the identity ordering)
220: Output Parameter:
221: . aoout - the new application mapping
223: Options Database Key:
224: . -ao_view - call `AOView()` at the conclusion of `AOCreateMapping()`
226: Level: beginner
228: Note:
229: 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.
230: Use `AOCreateBasic()` or `AOCreateBasicIS()` if they do not have holes for better performance.
232: .seealso: [](sec_ao), `AOCreateBasic()`, `AOCreateMappingIS()`, `AODestroy()`
233: @*/
234: PetscErrorCode AOCreateMapping(MPI_Comm comm, PetscInt napp, const PetscInt myapp[], const PetscInt mypetsc[], AO *aoout)
235: {
236: AO ao;
237: AO_Mapping *aomap;
238: PetscInt *allpetsc, *allapp;
239: PetscInt *petscPerm, *appPerm;
240: PetscInt *petsc;
241: PetscMPIInt size, rank, *lens, *disp, nnapp;
242: PetscInt N, start;
243: PetscInt i;
245: PetscFunctionBegin;
246: PetscAssertPointer(aoout, 5);
247: *aoout = NULL;
248: PetscCall(AOInitializePackage());
250: PetscCall(PetscHeaderCreate(ao, AO_CLASSID, "AO", "Application Ordering", "AO", comm, AODestroy, AOView));
251: PetscCall(PetscNew(&aomap));
252: ao->ops[0] = AOps;
253: ao->data = (void *)aomap;
255: /* transmit all lengths to all processors */
256: PetscCallMPI(MPI_Comm_size(comm, &size));
257: PetscCallMPI(MPI_Comm_rank(comm, &rank));
258: PetscCall(PetscMalloc2(size, &lens, size, &disp));
259: nnapp = napp;
260: PetscCallMPI(MPI_Allgather(&nnapp, 1, MPI_INT, lens, 1, MPI_INT, comm));
261: N = 0;
262: for (i = 0; i < size; i++) {
263: disp[i] = N;
264: N += lens[i];
265: }
266: aomap->N = N;
267: ao->N = N;
268: ao->n = N;
270: /* If mypetsc is 0 then use "natural" numbering */
271: if (!mypetsc) {
272: start = disp[rank];
273: PetscCall(PetscMalloc1(napp + 1, &petsc));
274: for (i = 0; i < napp; i++) petsc[i] = start + i;
275: } else {
276: petsc = (PetscInt *)mypetsc;
277: }
279: /* get all indices on all processors */
280: PetscCall(PetscMalloc4(N, &allapp, N, &appPerm, N, &allpetsc, N, &petscPerm));
281: PetscCallMPI(MPI_Allgatherv((void *)myapp, napp, MPIU_INT, allapp, lens, disp, MPIU_INT, comm));
282: PetscCallMPI(MPI_Allgatherv((void *)petsc, napp, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm));
283: PetscCall(PetscFree2(lens, disp));
285: /* generate a list of application and PETSc node numbers */
286: PetscCall(PetscMalloc4(N, &aomap->app, N, &aomap->appPerm, N, &aomap->petsc, N, &aomap->petscPerm));
287: for (i = 0; i < N; i++) {
288: appPerm[i] = i;
289: petscPerm[i] = i;
290: }
291: PetscCall(PetscSortIntWithPermutation(N, allpetsc, petscPerm));
292: PetscCall(PetscSortIntWithPermutation(N, allapp, appPerm));
293: /* Form sorted arrays of indices */
294: for (i = 0; i < N; i++) {
295: aomap->app[i] = allapp[appPerm[i]];
296: aomap->petsc[i] = allpetsc[petscPerm[i]];
297: }
298: /* Invert petscPerm[] into aomap->petscPerm[] */
299: for (i = 0; i < N; i++) aomap->petscPerm[petscPerm[i]] = i;
301: /* Form map between aomap->app[] and aomap->petsc[] */
302: for (i = 0; i < N; i++) aomap->appPerm[i] = aomap->petscPerm[appPerm[i]];
304: /* Invert appPerm[] into allapp[] */
305: for (i = 0; i < N; i++) allapp[appPerm[i]] = i;
307: /* Form map between aomap->petsc[] and aomap->app[] */
308: for (i = 0; i < N; i++) aomap->petscPerm[i] = allapp[petscPerm[i]];
310: if (PetscDefined(USE_DEBUG)) {
311: /* Check that the permutations are complementary */
312: for (i = 0; i < N; i++) PetscCheck(i == aomap->appPerm[aomap->petscPerm[i]], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid ordering");
313: }
314: /* Cleanup */
315: if (!mypetsc) PetscCall(PetscFree(petsc));
316: PetscCall(PetscFree4(allapp, appPerm, allpetsc, petscPerm));
318: PetscCall(AOViewFromOptions(ao, NULL, "-ao_view"));
320: *aoout = ao;
321: PetscFunctionReturn(PETSC_SUCCESS);
322: }
324: /*@
325: AOCreateMappingIS - Creates an application mapping using two index sets.
327: Input Parameters:
328: + isapp - index set that defines an ordering
329: - ispetsc - index set that defines another ordering, maybe `NULL` for identity `IS`
331: Output Parameter:
332: . aoout - the new application ordering
334: Options Database Key:
335: . -ao_view - call `AOView()` at the conclusion of `AOCreateMappingIS()`
337: Level: beginner
339: Note:
340: 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.
341: Use `AOCreateBasic()` or `AOCreateBasicIS()` if they do not have holes for better performance.
343: .seealso: [](sec_ao), [](sec_scatter), `AOCreateBasic()`, `AOCreateMapping()`, `AODestroy()`
344: @*/
345: PetscErrorCode AOCreateMappingIS(IS isapp, IS ispetsc, AO *aoout)
346: {
347: MPI_Comm comm;
348: const PetscInt *mypetsc, *myapp;
349: PetscInt napp, npetsc;
351: PetscFunctionBegin;
352: PetscCall(PetscObjectGetComm((PetscObject)isapp, &comm));
353: PetscCall(ISGetLocalSize(isapp, &napp));
354: if (ispetsc) {
355: PetscCall(ISGetLocalSize(ispetsc, &npetsc));
356: PetscCheck(napp == npetsc, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local IS lengths must match");
357: PetscCall(ISGetIndices(ispetsc, &mypetsc));
358: } else {
359: mypetsc = NULL;
360: }
361: PetscCall(ISGetIndices(isapp, &myapp));
363: PetscCall(AOCreateMapping(comm, napp, myapp, mypetsc, aoout));
365: PetscCall(ISRestoreIndices(isapp, &myapp));
366: if (ispetsc) PetscCall(ISRestoreIndices(ispetsc, &mypetsc));
367: PetscFunctionReturn(PETSC_SUCCESS);
368: }