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: PetscBool isascii;
33: PetscFunctionBegin;
34: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ao), &rank));
35: if (rank) PetscFunctionReturn(PETSC_SUCCESS);
36: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
37: if (isascii) {
38: PetscCall(PetscViewerASCIIPrintf(viewer, "Number of elements in ordering %" PetscInt_FMT "\n", aomap->N));
39: PetscCall(PetscViewerASCIIPrintf(viewer, " App. PETSc\n"));
40: for (PetscInt 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]]));
41: }
42: PetscFunctionReturn(PETSC_SUCCESS);
43: }
45: static 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;
55: /* It would be possible to use a single bisection search, which
56: recursively divided the indices to be converted, and searched
57: partitions which contained an index. This would result in
58: better running times if indices are clustered.
59: */
60: PetscFunctionBegin;
61: for (PetscInt 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: PetscFunctionReturn(PETSC_SUCCESS);
77: }
79: static 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;
89: /* It would be possible to use a single bisection search, which
90: recursively divided the indices to be converted, and searched
91: partitions which contained an index. This would result in
92: better running times if indices are clustered.
93: */
94: PetscFunctionBegin;
95: for (PetscInt 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: PetscFunctionReturn(PETSC_SUCCESS);
111: }
113: static const struct _AOOps AOps = {
114: PetscDesignatedInitializer(view, AOView_Mapping),
115: PetscDesignatedInitializer(destroy, AODestroy_Mapping),
116: PetscDesignatedInitializer(petsctoapplication, AOPetscToApplication_Mapping),
117: PetscDesignatedInitializer(applicationtopetsc, AOApplicationToPetsc_Mapping),
118: PetscDesignatedInitializer(petsctoapplicationpermuteint, NULL),
119: PetscDesignatedInitializer(applicationtopetscpermuteint, NULL),
120: PetscDesignatedInitializer(petsctoapplicationpermutereal, NULL),
121: PetscDesignatedInitializer(applicationtopetscpermutereal, NULL),
122: };
124: /*@
125: AOMappingHasApplicationIndex - Checks if an `AO` has a requested application index.
127: Not Collective
129: Input Parameters:
130: + ao - The `AO`
131: - idex - The application index
133: Output Parameter:
134: . hasIndex - Flag is `PETSC_TRUE` if the index exists
136: Level: intermediate
138: Developer Notes:
139: The name of the function is wrong, it should be `AOHasApplicationIndex`
141: .seealso: [](sec_ao), `AOMappingHasPetscIndex()`, `AOCreateMapping()`, `AO`
142: @*/
143: PetscErrorCode AOMappingHasApplicationIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
144: {
145: AO_Mapping *aomap;
146: PetscInt *app;
147: PetscInt low, high, mid;
149: PetscFunctionBegin;
151: PetscAssertPointer(hasIndex, 3);
152: aomap = (AO_Mapping *)ao->data;
153: app = aomap->app;
154: /* Use bisection since the array is sorted */
155: low = 0;
156: high = aomap->N - 1;
157: while (low <= high) {
158: mid = (low + high) / 2;
159: if (idex == app[mid]) break;
160: else if (idex < app[mid]) high = mid - 1;
161: else low = mid + 1;
162: }
163: if (low > high) *hasIndex = PETSC_FALSE;
164: else *hasIndex = PETSC_TRUE;
165: PetscFunctionReturn(PETSC_SUCCESS);
166: }
168: /*@
169: AOMappingHasPetscIndex - checks if an `AO` has a requested PETSc index.
171: Not Collective
173: Input Parameters:
174: + ao - The `AO`
175: - idex - The PETSc index
177: Output Parameter:
178: . hasIndex - Flag is `PETSC_TRUE` if the index exists
180: Level: intermediate
182: Developer Notes:
183: The name of the function is wrong, it should be `AOHasPetscIndex`
185: .seealso: [](sec_ao), `AOMappingHasApplicationIndex()`, `AOCreateMapping()`
186: @*/
187: PetscErrorCode AOMappingHasPetscIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
188: {
189: AO_Mapping *aomap;
190: PetscInt *petsc;
191: PetscInt low, high, mid;
193: PetscFunctionBegin;
195: PetscAssertPointer(hasIndex, 3);
196: aomap = (AO_Mapping *)ao->data;
197: petsc = aomap->petsc;
198: /* Use bisection since the array is sorted */
199: low = 0;
200: high = aomap->N - 1;
201: while (low <= high) {
202: mid = (low + high) / 2;
203: if (idex == petsc[mid]) break;
204: else if (idex < petsc[mid]) high = mid - 1;
205: else low = mid + 1;
206: }
207: if (low > high) *hasIndex = PETSC_FALSE;
208: else *hasIndex = PETSC_TRUE;
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
212: /*@
213: AOCreateMapping - Creates an application mapping using two integer arrays.
215: Input Parameters:
216: + comm - MPI communicator that is to share the `AO`
217: . napp - size of integer arrays
218: . myapp - integer array that defines an ordering
219: - mypetsc - integer array that defines another ordering (may be `NULL` to indicate the identity ordering)
221: Output Parameter:
222: . aoout - the new application mapping
224: Options Database Key:
225: . -ao_view - call `AOView()` at the conclusion of `AOCreateMapping()`
227: Level: beginner
229: Note:
230: 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.
231: Use `AOCreateBasic()` or `AOCreateBasicIS()` if they do not have holes for better performance.
233: .seealso: [](sec_ao), `AOCreateBasic()`, `AOCreateMappingIS()`, `AODestroy()`
234: @*/
235: PetscErrorCode AOCreateMapping(MPI_Comm comm, PetscInt napp, const PetscInt myapp[], const PetscInt mypetsc[], AO *aoout)
236: {
237: AO ao;
238: AO_Mapping *aomap;
239: PetscInt *allpetsc, *allapp;
240: PetscInt *petscPerm, *appPerm;
241: PetscInt *petsc;
242: PetscMPIInt size, rank, *lens, *disp, nnapp;
243: PetscCount N;
244: PetscInt start;
246: PetscFunctionBegin;
247: PetscAssertPointer(aoout, 5);
248: *aoout = NULL;
249: PetscCall(AOInitializePackage());
251: PetscCall(PetscHeaderCreate(ao, AO_CLASSID, "AO", "Application Ordering", "AO", comm, AODestroy, AOView));
252: PetscCall(PetscNew(&aomap));
253: ao->ops[0] = AOps;
254: ao->data = (void *)aomap;
256: /* transmit all lengths to all processors */
257: PetscCallMPI(MPI_Comm_size(comm, &size));
258: PetscCallMPI(MPI_Comm_rank(comm, &rank));
259: PetscCall(PetscMalloc2(size, &lens, size, &disp));
260: PetscCall(PetscMPIIntCast(napp, &nnapp));
261: PetscCallMPI(MPI_Allgather(&nnapp, 1, MPI_INT, lens, 1, MPI_INT, comm));
262: N = 0;
263: for (PetscMPIInt i = 0; i < size; i++) {
264: PetscCall(PetscMPIIntCast(N, &disp[i]));
265: N += lens[i];
266: }
267: PetscCall(PetscIntCast(N, &aomap->N));
268: ao->N = aomap->N;
269: ao->n = aomap->N;
271: /* If mypetsc is 0 then use "natural" numbering */
272: if (!mypetsc) {
273: start = disp[rank];
274: PetscCall(PetscMalloc1(napp + 1, &petsc));
275: for (PetscInt i = 0; i < napp; i++) petsc[i] = start + i;
276: } else {
277: petsc = (PetscInt *)mypetsc;
278: }
280: /* get all indices on all processors */
281: PetscCall(PetscMalloc4(N, &allapp, N, &appPerm, N, &allpetsc, N, &petscPerm));
282: PetscCallMPI(MPI_Allgatherv((void *)myapp, nnapp, MPIU_INT, allapp, lens, disp, MPIU_INT, comm));
283: PetscCallMPI(MPI_Allgatherv((void *)petsc, nnapp, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm));
284: PetscCall(PetscFree2(lens, disp));
286: /* generate a list of application and PETSc node numbers */
287: PetscCall(PetscMalloc4(N, &aomap->app, N, &aomap->appPerm, N, &aomap->petsc, N, &aomap->petscPerm));
288: for (PetscInt i = 0; i < N; i++) {
289: appPerm[i] = i;
290: petscPerm[i] = i;
291: }
292: PetscCall(PetscSortIntWithPermutation(aomap->N, allpetsc, petscPerm));
293: PetscCall(PetscSortIntWithPermutation(aomap->N, allapp, appPerm));
294: /* Form sorted arrays of indices */
295: for (PetscInt i = 0; i < N; i++) {
296: aomap->app[i] = allapp[appPerm[i]];
297: aomap->petsc[i] = allpetsc[petscPerm[i]];
298: }
299: /* Invert petscPerm[] into aomap->petscPerm[] */
300: for (PetscInt i = 0; i < N; i++) aomap->petscPerm[petscPerm[i]] = i;
302: /* Form map between aomap->app[] and aomap->petsc[] */
303: for (PetscInt i = 0; i < N; i++) aomap->appPerm[i] = aomap->petscPerm[appPerm[i]];
305: /* Invert appPerm[] into allapp[] */
306: for (PetscInt i = 0; i < N; i++) allapp[appPerm[i]] = i;
308: /* Form map between aomap->petsc[] and aomap->app[] */
309: for (PetscInt i = 0; i < N; i++) aomap->petscPerm[i] = allapp[petscPerm[i]];
311: if (PetscDefined(USE_DEBUG)) {
312: /* Check that the permutations are complementary */
313: for (PetscInt i = 0; i < N; i++) PetscCheck(i == aomap->appPerm[aomap->petscPerm[i]], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid ordering");
314: }
315: /* Cleanup */
316: if (!mypetsc) PetscCall(PetscFree(petsc));
317: PetscCall(PetscFree4(allapp, appPerm, allpetsc, petscPerm));
318: PetscCall(AOViewFromOptions(ao, NULL, "-ao_view"));
319: *aoout = ao;
320: PetscFunctionReturn(PETSC_SUCCESS);
321: }
323: /*@
324: AOCreateMappingIS - Creates an application mapping using two index sets.
326: Input Parameters:
327: + isapp - index set that defines an ordering
328: - ispetsc - index set that defines another ordering, maybe `NULL` for identity `IS`
330: Output Parameter:
331: . aoout - the new application ordering
333: Options Database Key:
334: . -ao_view - call `AOView()` at the conclusion of `AOCreateMappingIS()`
336: Level: beginner
338: Note:
339: 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.
340: Use `AOCreateBasic()` or `AOCreateBasicIS()` if they do not have holes for better performance.
342: .seealso: [](sec_ao), [](sec_scatter), `AOCreateBasic()`, `AOCreateMapping()`, `AODestroy()`
343: @*/
344: PetscErrorCode AOCreateMappingIS(IS isapp, IS ispetsc, AO *aoout)
345: {
346: MPI_Comm comm;
347: const PetscInt *mypetsc, *myapp;
348: PetscInt napp, npetsc;
350: PetscFunctionBegin;
351: PetscCall(PetscObjectGetComm((PetscObject)isapp, &comm));
352: PetscCall(ISGetLocalSize(isapp, &napp));
353: if (ispetsc) {
354: PetscCall(ISGetLocalSize(ispetsc, &npetsc));
355: PetscCheck(napp == npetsc, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local IS lengths must match");
356: PetscCall(ISGetIndices(ispetsc, &mypetsc));
357: } else {
358: mypetsc = NULL;
359: }
360: PetscCall(ISGetIndices(isapp, &myapp));
362: PetscCall(AOCreateMapping(comm, napp, myapp, mypetsc, aoout));
364: PetscCall(ISRestoreIndices(isapp, &myapp));
365: if (ispetsc) PetscCall(ISRestoreIndices(ispetsc, &mypetsc));
366: PetscFunctionReturn(PETSC_SUCCESS);
367: }