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 */
6
7 #include <../src/vec/is/ao/aoimpl.h> /*I "petscao.h" I*/
8
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;
16
AODestroy_Mapping(AO ao)17 static PetscErrorCode AODestroy_Mapping(AO ao)
18 {
19 AO_Mapping *aomap = (AO_Mapping *)ao->data;
20
21 PetscFunctionBegin;
22 PetscCall(PetscFree4(aomap->app, aomap->appPerm, aomap->petsc, aomap->petscPerm));
23 PetscCall(PetscFree(aomap));
24 PetscFunctionReturn(PETSC_SUCCESS);
25 }
26
AOView_Mapping(AO ao,PetscViewer viewer)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 isascii;
33
34 PetscFunctionBegin;
35 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ao), &rank));
36 if (rank) PetscFunctionReturn(PETSC_SUCCESS);
37 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
38 if (isascii) {
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 }
45
AOPetscToApplication_Mapping(AO ao,PetscInt n,PetscInt * ia)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;
56
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 }
80
AOApplicationToPetsc_Mapping(AO ao,PetscInt n,PetscInt * ia)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;
91
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 }
115
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 };
126
127 /*@
128 AOMappingHasApplicationIndex - Checks if an `AO` has a requested application index.
129
130 Not Collective
131
132 Input Parameters:
133 + ao - The `AO`
134 - idex - The application index
135
136 Output Parameter:
137 . hasIndex - Flag is `PETSC_TRUE` if the index exists
138
139 Level: intermediate
140
141 Developer Notes:
142 The name of the function is wrong, it should be `AOHasApplicationIndex`
143
144 .seealso: [](sec_ao), `AOMappingHasPetscIndex()`, `AOCreateMapping()`, `AO`
145 @*/
AOMappingHasApplicationIndex(AO ao,PetscInt idex,PetscBool * hasIndex)146 PetscErrorCode AOMappingHasApplicationIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
147 {
148 AO_Mapping *aomap;
149 PetscInt *app;
150 PetscInt low, high, mid;
151
152 PetscFunctionBegin;
153 PetscValidHeaderSpecific(ao, AO_CLASSID, 1);
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 }
170
171 /*@
172 AOMappingHasPetscIndex - checks if an `AO` has a requested PETSc index.
173
174 Not Collective
175
176 Input Parameters:
177 + ao - The `AO`
178 - idex - The PETSc index
179
180 Output Parameter:
181 . hasIndex - Flag is `PETSC_TRUE` if the index exists
182
183 Level: intermediate
184
185 Developer Notes:
186 The name of the function is wrong, it should be `AOHasPetscIndex`
187
188 .seealso: [](sec_ao), `AOMappingHasApplicationIndex()`, `AOCreateMapping()`
189 @*/
AOMappingHasPetscIndex(AO ao,PetscInt idex,PetscBool * hasIndex)190 PetscErrorCode AOMappingHasPetscIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
191 {
192 AO_Mapping *aomap;
193 PetscInt *petsc;
194 PetscInt low, high, mid;
195
196 PetscFunctionBegin;
197 PetscValidHeaderSpecific(ao, AO_CLASSID, 1);
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 }
214
215 /*@
216 AOCreateMapping - Creates an application mapping using two integer arrays.
217
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)
223
224 Output Parameter:
225 . aoout - the new application mapping
226
227 Options Database Key:
228 . -ao_view - call `AOView()` at the conclusion of `AOCreateMapping()`
229
230 Level: beginner
231
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.
235
236 .seealso: [](sec_ao), `AOCreateBasic()`, `AOCreateMappingIS()`, `AODestroy()`
237 @*/
AOCreateMapping(MPI_Comm comm,PetscInt napp,const PetscInt myapp[],const PetscInt mypetsc[],AO * aoout)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;
248
249 PetscFunctionBegin;
250 PetscAssertPointer(aoout, 5);
251 *aoout = NULL;
252 PetscCall(AOInitializePackage());
253
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;
258
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 = aomap->N;
272 ao->n = aomap->N;
273
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 }
282
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));
288
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(aomap->N, allpetsc, petscPerm));
296 PetscCall(PetscSortIntWithPermutation(aomap->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;
304
305 /* Form map between aomap->app[] and aomap->petsc[] */
306 for (PetscInt i = 0; i < N; i++) aomap->appPerm[i] = aomap->petscPerm[appPerm[i]];
307
308 /* Invert appPerm[] into allapp[] */
309 for (PetscInt i = 0; i < N; i++) allapp[appPerm[i]] = i;
310
311 /* Form map between aomap->petsc[] and aomap->app[] */
312 for (PetscInt i = 0; i < N; i++) aomap->petscPerm[i] = allapp[petscPerm[i]];
313
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 }
325
326 /*@
327 AOCreateMappingIS - Creates an application mapping using two index sets.
328
329 Input Parameters:
330 + isapp - index set that defines an ordering
331 - ispetsc - index set that defines another ordering, maybe `NULL` for identity `IS`
332
333 Output Parameter:
334 . aoout - the new application ordering
335
336 Options Database Key:
337 . -ao_view - call `AOView()` at the conclusion of `AOCreateMappingIS()`
338
339 Level: beginner
340
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.
344
345 .seealso: [](sec_ao), [](sec_scatter), `AOCreateBasic()`, `AOCreateMapping()`, `AODestroy()`
346 @*/
AOCreateMappingIS(IS isapp,IS ispetsc,AO * aoout)347 PetscErrorCode AOCreateMappingIS(IS isapp, IS ispetsc, AO *aoout)
348 {
349 MPI_Comm comm;
350 const PetscInt *mypetsc, *myapp;
351 PetscInt napp, npetsc;
352
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));
364
365 PetscCall(AOCreateMapping(comm, napp, myapp, mypetsc, aoout));
366
367 PetscCall(ISRestoreIndices(isapp, &myapp));
368 if (ispetsc) PetscCall(ISRestoreIndices(ispetsc, &mypetsc));
369 PetscFunctionReturn(PETSC_SUCCESS);
370 }
371