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: PetscFunctionBegin;
23: PetscCall(PetscFree4(aomap->app, aomap->appPerm, aomap->petsc, aomap->petscPerm));
24: PetscCall(PetscFree(aomap));
25: PetscFunctionReturn(PETSC_SUCCESS);
26: }
28: PetscErrorCode AOView_Mapping(AO ao, PetscViewer viewer)
29: {
30: AO_Mapping *aomap = (AO_Mapping *)ao->data;
31: PetscMPIInt rank;
32: PetscInt i;
33: PetscBool iascii;
35: PetscFunctionBegin;
36: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ao), &rank));
37: if (rank) PetscFunctionReturn(PETSC_SUCCESS);
38: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
39: if (iascii) {
40: PetscCall(PetscViewerASCIIPrintf(viewer, "Number of elements in ordering %" PetscInt_FMT "\n", aomap->N));
41: PetscCall(PetscViewerASCIIPrintf(viewer, " App. PETSc\n"));
42: 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]]));
43: }
44: PetscFunctionReturn(PETSC_SUCCESS);
45: }
47: PetscErrorCode AOPetscToApplication_Mapping(AO ao, PetscInt n, PetscInt *ia)
48: {
49: AO_Mapping *aomap = (AO_Mapping *)ao->data;
50: PetscInt *app = aomap->app;
51: PetscInt *petsc = aomap->petsc;
52: PetscInt *perm = aomap->petscPerm;
53: PetscInt N = aomap->N;
54: PetscInt low, high, mid = 0;
55: PetscInt idex;
56: PetscInt i;
58: /* It would be possible to use a single bisection search, which
59: recursively divided the indices to be converted, and searched
60: partitions which contained an index. This would result in
61: better running times if indices are clustered.
62: */
63: PetscFunctionBegin;
64: for (i = 0; i < n; i++) {
65: idex = ia[i];
66: if (idex < 0) continue;
67: /* Use bisection since the array is sorted */
68: low = 0;
69: high = N - 1;
70: while (low <= high) {
71: mid = (low + high) / 2;
72: if (idex == petsc[mid]) break;
73: else if (idex < petsc[mid]) high = mid - 1;
74: else low = mid + 1;
75: }
76: if (low > high) ia[i] = -1;
77: else ia[i] = app[perm[mid]];
78: }
79: PetscFunctionReturn(PETSC_SUCCESS);
80: }
82: PetscErrorCode AOApplicationToPetsc_Mapping(AO ao, PetscInt n, PetscInt *ia)
83: {
84: AO_Mapping *aomap = (AO_Mapping *)ao->data;
85: PetscInt *app = aomap->app;
86: PetscInt *petsc = aomap->petsc;
87: PetscInt *perm = aomap->appPerm;
88: PetscInt N = aomap->N;
89: PetscInt low, high, mid = 0;
90: PetscInt idex;
91: PetscInt i;
93: /* It would be possible to use a single bisection search, which
94: recursively divided the indices to be converted, and searched
95: partitions which contained an index. This would result in
96: better running times if indices are clustered.
97: */
98: PetscFunctionBegin;
99: for (i = 0; i < n; i++) {
100: idex = ia[i];
101: if (idex < 0) continue;
102: /* Use bisection since the array is sorted */
103: low = 0;
104: high = N - 1;
105: while (low <= high) {
106: mid = (low + high) / 2;
107: if (idex == app[mid]) break;
108: else if (idex < app[mid]) high = mid - 1;
109: else low = mid + 1;
110: }
111: if (low > high) ia[i] = -1;
112: else ia[i] = petsc[perm[mid]];
113: }
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }
117: static struct _AOOps AOps = {
118: PetscDesignatedInitializer(view, AOView_Mapping),
119: PetscDesignatedInitializer(destroy, AODestroy_Mapping),
120: PetscDesignatedInitializer(petsctoapplication, AOPetscToApplication_Mapping),
121: PetscDesignatedInitializer(applicationtopetsc, AOApplicationToPetsc_Mapping),
122: };
124: /*@C
125: AOMappingHasApplicationIndex - Checks if an `AO` has a requested application index.
127: Not Collective
129: Input Parameters:
130: + ao - The `AO`
131: - index - The application index
133: Output Parameter:
134: . hasIndex - Flag is `PETSC_TRUE` if the index exists
136: Level: intermediate
138: Developer Note:
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;
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: - index - The petsc index
177: Output Parameter:
178: . hasIndex - Flag is `PETSC_TRUE` if the index exists
180: Level: intermediate
182: Developer Note:
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;
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: /*@C
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()`, `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: PetscInt N, start;
244: PetscInt i;
246: PetscFunctionBegin;
248: *aoout = NULL;
249: PetscCall(AOInitializePackage());
251: PetscCall(PetscHeaderCreate(ao, AO_CLASSID, "AO", "Application Ordering", "AO", comm, AODestroy, AOView));
252: PetscCall(PetscNew(&aomap));
253: PetscCall(PetscMemcpy(ao->ops, &AOps, sizeof(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: nnapp = napp;
261: PetscCallMPI(MPI_Allgather(&nnapp, 1, MPI_INT, lens, 1, MPI_INT, comm));
262: N = 0;
263: for (i = 0; i < size; i++) {
264: disp[i] = N;
265: N += lens[i];
266: }
267: aomap->N = N;
268: ao->N = N;
269: ao->n = 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 (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, napp, MPIU_INT, allapp, lens, disp, MPIU_INT, comm));
283: PetscCallMPI(MPI_Allgatherv((void *)petsc, napp, 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 (i = 0; i < N; i++) {
289: appPerm[i] = i;
290: petscPerm[i] = i;
291: }
292: PetscCall(PetscSortIntWithPermutation(N, allpetsc, petscPerm));
293: PetscCall(PetscSortIntWithPermutation(N, allapp, appPerm));
294: /* Form sorted arrays of indices */
295: for (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 (i = 0; i < N; i++) aomap->petscPerm[petscPerm[i]] = i;
302: /* Form map between aomap->app[] and aomap->petsc[] */
303: for (i = 0; i < N; i++) aomap->appPerm[i] = aomap->petscPerm[appPerm[i]];
305: /* Invert appPerm[] into allapp[] */
306: for (i = 0; i < N; i++) allapp[appPerm[i]] = i;
308: /* Form map between aomap->petsc[] and aomap->app[] */
309: for (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 (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));
319: PetscCall(AOViewFromOptions(ao, NULL, "-ao_view"));
321: *aoout = ao;
322: PetscFunctionReturn(PETSC_SUCCESS);
323: }
325: /*@
326: AOCreateMappingIS - Creates an application mapping using two index sets.
328: Input Parameters:
329: + comm - MPI communicator that is to share `AO`
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: }