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: PetscDesignatedInitializer(petsctoapplicationpermuteint, NULL),
122: PetscDesignatedInitializer(applicationtopetscpermuteint, NULL),
123: PetscDesignatedInitializer(petsctoapplicationpermutereal, NULL),
124: PetscDesignatedInitializer(applicationtopetscpermutereal, NULL),
125: };
127: /*@
128: AOMappingHasApplicationIndex - Checks if an `AO` has a requested application index.
130: Not Collective
132: Input Parameters:
133: + ao - The `AO`
134: - idex - The application index
136: Output Parameter:
137: . hasIndex - Flag is `PETSC_TRUE` if the index exists
139: Level: intermediate
141: Developer Notes:
142: The name of the function is wrong, it should be `AOHasApplicationIndex`
144: .seealso: [](sec_ao), `AOMappingHasPetscIndex()`, `AOCreateMapping()`, `AO`
145: @*/
146: PetscErrorCode AOMappingHasApplicationIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
147: {
148: AO_Mapping *aomap;
149: PetscInt *app;
150: PetscInt low, high, mid;
152: PetscFunctionBegin;
154: PetscAssertPointer(hasIndex, 3);
155: aomap = (AO_Mapping *)ao->data;
156: app = aomap->app;
157: /* Use bisection since the array is sorted */
158: low = 0;
159: high = aomap->N - 1;
160: while (low <= high) {
161: mid = (low + high) / 2;
162: if (idex == app[mid]) break;
163: else if (idex < app[mid]) high = mid - 1;
164: else low = mid + 1;
165: }
166: if (low > high) *hasIndex = PETSC_FALSE;
167: else *hasIndex = PETSC_TRUE;
168: PetscFunctionReturn(PETSC_SUCCESS);
169: }
171: /*@
172: AOMappingHasPetscIndex - checks if an `AO` has a requested petsc index.
174: Not Collective
176: Input Parameters:
177: + ao - The `AO`
178: - idex - The petsc index
180: Output Parameter:
181: . hasIndex - Flag is `PETSC_TRUE` if the index exists
183: Level: intermediate
185: Developer Notes:
186: The name of the function is wrong, it should be `AOHasPetscIndex`
188: .seealso: [](sec_ao), `AOMappingHasApplicationIndex()`, `AOCreateMapping()`
189: @*/
190: PetscErrorCode AOMappingHasPetscIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
191: {
192: AO_Mapping *aomap;
193: PetscInt *petsc;
194: PetscInt low, high, mid;
196: PetscFunctionBegin;
198: PetscAssertPointer(hasIndex, 3);
199: aomap = (AO_Mapping *)ao->data;
200: petsc = aomap->petsc;
201: /* Use bisection since the array is sorted */
202: low = 0;
203: high = aomap->N - 1;
204: while (low <= high) {
205: mid = (low + high) / 2;
206: if (idex == petsc[mid]) break;
207: else if (idex < petsc[mid]) high = mid - 1;
208: else low = mid + 1;
209: }
210: if (low > high) *hasIndex = PETSC_FALSE;
211: else *hasIndex = PETSC_TRUE;
212: PetscFunctionReturn(PETSC_SUCCESS);
213: }
215: /*@
216: AOCreateMapping - Creates an application mapping using two integer arrays.
218: Input Parameters:
219: + comm - MPI communicator that is to share the `AO`
220: . napp - size of integer arrays
221: . myapp - integer array that defines an ordering
222: - mypetsc - integer array that defines another ordering (may be `NULL` to indicate the identity ordering)
224: Output Parameter:
225: . aoout - the new application mapping
227: Options Database Key:
228: . -ao_view - call `AOView()` at the conclusion of `AOCreateMapping()`
230: Level: beginner
232: Note:
233: 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.
234: Use `AOCreateBasic()` or `AOCreateBasicIS()` if they do not have holes for better performance.
236: .seealso: [](sec_ao), `AOCreateBasic()`, `AOCreateMappingIS()`, `AODestroy()`
237: @*/
238: PetscErrorCode AOCreateMapping(MPI_Comm comm, PetscInt napp, const PetscInt myapp[], const PetscInt mypetsc[], AO *aoout)
239: {
240: AO ao;
241: AO_Mapping *aomap;
242: PetscInt *allpetsc, *allapp;
243: PetscInt *petscPerm, *appPerm;
244: PetscInt *petsc;
245: PetscMPIInt size, rank, *lens, *disp, nnapp;
246: PetscCount N;
247: PetscInt start;
249: PetscFunctionBegin;
250: PetscAssertPointer(aoout, 5);
251: *aoout = NULL;
252: PetscCall(AOInitializePackage());
254: PetscCall(PetscHeaderCreate(ao, AO_CLASSID, "AO", "Application Ordering", "AO", comm, AODestroy, AOView));
255: PetscCall(PetscNew(&aomap));
256: ao->ops[0] = AOps;
257: ao->data = (void *)aomap;
259: /* transmit all lengths to all processors */
260: PetscCallMPI(MPI_Comm_size(comm, &size));
261: PetscCallMPI(MPI_Comm_rank(comm, &rank));
262: PetscCall(PetscMalloc2(size, &lens, size, &disp));
263: PetscCall(PetscMPIIntCast(napp, &nnapp));
264: PetscCallMPI(MPI_Allgather(&nnapp, 1, MPI_INT, lens, 1, MPI_INT, comm));
265: N = 0;
266: for (PetscMPIInt i = 0; i < size; i++) {
267: PetscCall(PetscMPIIntCast(N, &disp[i]));
268: N += lens[i];
269: }
270: PetscCall(PetscIntCast(N, &aomap->N));
271: ao->N = (PetscInt)N;
272: ao->n = (PetscInt)N;
274: /* If mypetsc is 0 then use "natural" numbering */
275: if (!mypetsc) {
276: start = disp[rank];
277: PetscCall(PetscMalloc1(napp + 1, &petsc));
278: for (PetscInt i = 0; i < napp; i++) petsc[i] = start + i;
279: } else {
280: petsc = (PetscInt *)mypetsc;
281: }
283: /* get all indices on all processors */
284: PetscCall(PetscMalloc4(N, &allapp, N, &appPerm, N, &allpetsc, N, &petscPerm));
285: PetscCallMPI(MPI_Allgatherv((void *)myapp, nnapp, MPIU_INT, allapp, lens, disp, MPIU_INT, comm));
286: PetscCallMPI(MPI_Allgatherv((void *)petsc, nnapp, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm));
287: PetscCall(PetscFree2(lens, disp));
289: /* generate a list of application and PETSc node numbers */
290: PetscCall(PetscMalloc4(N, &aomap->app, N, &aomap->appPerm, N, &aomap->petsc, N, &aomap->petscPerm));
291: for (PetscInt i = 0; i < N; i++) {
292: appPerm[i] = i;
293: petscPerm[i] = i;
294: }
295: PetscCall(PetscSortIntWithPermutation((PetscInt)N, allpetsc, petscPerm));
296: PetscCall(PetscSortIntWithPermutation((PetscInt)N, allapp, appPerm));
297: /* Form sorted arrays of indices */
298: for (PetscInt i = 0; i < N; i++) {
299: aomap->app[i] = allapp[appPerm[i]];
300: aomap->petsc[i] = allpetsc[petscPerm[i]];
301: }
302: /* Invert petscPerm[] into aomap->petscPerm[] */
303: for (PetscInt i = 0; i < N; i++) aomap->petscPerm[petscPerm[i]] = i;
305: /* Form map between aomap->app[] and aomap->petsc[] */
306: for (PetscInt i = 0; i < N; i++) aomap->appPerm[i] = aomap->petscPerm[appPerm[i]];
308: /* Invert appPerm[] into allapp[] */
309: for (PetscInt i = 0; i < N; i++) allapp[appPerm[i]] = i;
311: /* Form map between aomap->petsc[] and aomap->app[] */
312: for (PetscInt i = 0; i < N; i++) aomap->petscPerm[i] = allapp[petscPerm[i]];
314: if (PetscDefined(USE_DEBUG)) {
315: /* Check that the permutations are complementary */
316: for (PetscInt i = 0; i < N; i++) PetscCheck(i == aomap->appPerm[aomap->petscPerm[i]], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid ordering");
317: }
318: /* Cleanup */
319: if (!mypetsc) PetscCall(PetscFree(petsc));
320: PetscCall(PetscFree4(allapp, appPerm, allpetsc, petscPerm));
321: PetscCall(AOViewFromOptions(ao, NULL, "-ao_view"));
322: *aoout = ao;
323: PetscFunctionReturn(PETSC_SUCCESS);
324: }
326: /*@
327: AOCreateMappingIS - Creates an application mapping using two index sets.
329: Input Parameters:
330: + isapp - index set that defines an ordering
331: - ispetsc - index set that defines another ordering, maybe `NULL` for identity `IS`
333: Output Parameter:
334: . aoout - the new application ordering
336: Options Database Key:
337: . -ao_view - call `AOView()` at the conclusion of `AOCreateMappingIS()`
339: Level: beginner
341: Note:
342: 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.
343: Use `AOCreateBasic()` or `AOCreateBasicIS()` if they do not have holes for better performance.
345: .seealso: [](sec_ao), [](sec_scatter), `AOCreateBasic()`, `AOCreateMapping()`, `AODestroy()`
346: @*/
347: PetscErrorCode AOCreateMappingIS(IS isapp, IS ispetsc, AO *aoout)
348: {
349: MPI_Comm comm;
350: const PetscInt *mypetsc, *myapp;
351: PetscInt napp, npetsc;
353: PetscFunctionBegin;
354: PetscCall(PetscObjectGetComm((PetscObject)isapp, &comm));
355: PetscCall(ISGetLocalSize(isapp, &napp));
356: if (ispetsc) {
357: PetscCall(ISGetLocalSize(ispetsc, &npetsc));
358: PetscCheck(napp == npetsc, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local IS lengths must match");
359: PetscCall(ISGetIndices(ispetsc, &mypetsc));
360: } else {
361: mypetsc = NULL;
362: }
363: PetscCall(ISGetIndices(isapp, &myapp));
365: PetscCall(AOCreateMapping(comm, napp, myapp, mypetsc, aoout));
367: PetscCall(ISRestoreIndices(isapp, &myapp));
368: if (ispetsc) PetscCall(ISRestoreIndices(ispetsc, &mypetsc));
369: PetscFunctionReturn(PETSC_SUCCESS);
370: }