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: PetscFree4(aomap->app,aomap->appPerm,aomap->petsc,aomap->petscPerm);
23: PetscFree(aomap);
24: return 0;
25: }
27: 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: MPI_Comm_rank(PetscObjectComm((PetscObject)ao), &rank);
35: if (rank) return 0;
36: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
37: if (iascii) {
38: PetscViewerASCIIPrintf(viewer, "Number of elements in ordering %" PetscInt_FMT "\n", aomap->N);
39: PetscViewerASCIIPrintf(viewer, " App. PETSc\n");
40: for (i = 0; i < aomap->N; i++) {
41: PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", i, aomap->app[i], aomap->petsc[aomap->appPerm[i]]);
42: }
43: }
44: return 0;
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: 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: return 0;
79: }
81: 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: for (i = 0; i < n; i++) {
98: idex = ia[i];
99: if (idex < 0) continue;
100: /* Use bisection since the array is sorted */
101: low = 0;
102: high = N - 1;
103: while (low <= high) {
104: mid = (low + high)/2;
105: if (idex == app[mid]) break;
106: else if (idex < app[mid]) high = mid - 1;
107: else low = mid + 1;
108: }
109: if (low > high) ia[i] = -1;
110: else ia[i] = petsc[perm[mid]];
111: }
112: return 0;
113: }
115: static struct _AOOps AOps = {
116: PetscDesignatedInitializer(view,AOView_Mapping),
117: PetscDesignatedInitializer(destroy,AODestroy_Mapping),
118: PetscDesignatedInitializer(petsctoapplication,AOPetscToApplication_Mapping),
119: PetscDesignatedInitializer(applicationtopetsc,AOApplicationToPetsc_Mapping),
120: };
122: /*@C
123: AOMappingHasApplicationIndex - Searches for the supplied application index.
125: Input Parameters:
126: + ao - The AOMapping
127: - index - The application index
129: Output Parameter:
130: . hasIndex - Flag is PETSC_TRUE if the index exists
132: Level: intermediate
134: .seealso: AOMappingHasPetscIndex(), AOCreateMapping()
135: @*/
136: PetscErrorCode AOMappingHasApplicationIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
137: {
138: AO_Mapping *aomap;
139: PetscInt *app;
140: PetscInt low, high, mid;
144: aomap = (AO_Mapping*) ao->data;
145: app = aomap->app;
146: /* Use bisection since the array is sorted */
147: low = 0;
148: high = aomap->N - 1;
149: while (low <= high) {
150: mid = (low + high)/2;
151: if (idex == app[mid]) break;
152: else if (idex < app[mid]) high = mid - 1;
153: else low = mid + 1;
154: }
155: if (low > high) *hasIndex = PETSC_FALSE;
156: else *hasIndex = PETSC_TRUE;
157: return 0;
158: }
160: /*@
161: AOMappingHasPetscIndex - Searches for the supplied petsc index.
163: Input Parameters:
164: + ao - The AOMapping
165: - index - The petsc index
167: Output Parameter:
168: . hasIndex - Flag is PETSC_TRUE if the index exists
170: Level: intermediate
172: .seealso: AOMappingHasApplicationIndex(), AOCreateMapping()
173: @*/
174: PetscErrorCode AOMappingHasPetscIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
175: {
176: AO_Mapping *aomap;
177: PetscInt *petsc;
178: PetscInt low, high, mid;
182: aomap = (AO_Mapping*) ao->data;
183: petsc = aomap->petsc;
184: /* Use bisection since the array is sorted */
185: low = 0;
186: high = aomap->N - 1;
187: while (low <= high) {
188: mid = (low + high)/2;
189: if (idex == petsc[mid]) break;
190: else if (idex < petsc[mid]) high = mid - 1;
191: else low = mid + 1;
192: }
193: if (low > high) *hasIndex = PETSC_FALSE;
194: else *hasIndex = PETSC_TRUE;
195: return 0;
196: }
198: /*@C
199: AOCreateMapping - Creates a basic application mapping using two integer arrays.
201: Input Parameters:
202: + comm - MPI communicator that is to share AO
203: . napp - size of integer arrays
204: . myapp - integer array that defines an ordering
205: - mypetsc - integer array that defines another ordering (may be NULL to indicate the identity ordering)
207: Output Parameter:
208: . aoout - the new application mapping
210: Options Database Key:
211: . -ao_view - call AOView() at the conclusion of AOCreateMapping()
213: Level: beginner
215: Notes:
216: 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.
217: Use AOCreateBasic() or AOCreateBasicIS() if they do not have holes for better performance.
219: .seealso: AOCreateBasic(), AOCreateBasic(), AOCreateMappingIS(), AODestroy()
220: @*/
221: PetscErrorCode AOCreateMapping(MPI_Comm comm,PetscInt napp,const PetscInt myapp[],const PetscInt mypetsc[],AO *aoout)
222: {
223: AO ao;
224: AO_Mapping *aomap;
225: PetscInt *allpetsc, *allapp;
226: PetscInt *petscPerm, *appPerm;
227: PetscInt *petsc;
228: PetscMPIInt size, rank,*lens, *disp,nnapp;
229: PetscInt N, start;
230: PetscInt i;
233: *aoout = NULL;
234: AOInitializePackage();
236: PetscHeaderCreate(ao, AO_CLASSID, "AO", "Application Ordering", "AO", comm, AODestroy, AOView);
237: PetscNewLog(ao,&aomap);
238: PetscMemcpy(ao->ops, &AOps, sizeof(AOps));
239: ao->data = (void*) aomap;
241: /* transmit all lengths to all processors */
242: MPI_Comm_size(comm, &size);
243: MPI_Comm_rank(comm, &rank);
244: PetscMalloc2(size, &lens,size,&disp);
245: nnapp = napp;
246: MPI_Allgather(&nnapp, 1, MPI_INT, lens, 1, MPI_INT, comm);
247: N = 0;
248: for (i = 0; i < size; i++) {
249: disp[i] = N;
250: N += lens[i];
251: }
252: aomap->N = N;
253: ao->N = N;
254: ao->n = N;
256: /* If mypetsc is 0 then use "natural" numbering */
257: if (!mypetsc) {
258: start = disp[rank];
259: PetscMalloc1(napp+1, &petsc);
260: for (i = 0; i < napp; i++) petsc[i] = start + i;
261: } else {
262: petsc = (PetscInt*)mypetsc;
263: }
265: /* get all indices on all processors */
266: PetscMalloc4(N, &allapp,N,&appPerm,N,&allpetsc,N,&petscPerm);
267: MPI_Allgatherv((void*)myapp, napp, MPIU_INT, allapp, lens, disp, MPIU_INT, comm);
268: MPI_Allgatherv((void*)petsc, napp, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm);
269: PetscFree2(lens,disp);
271: /* generate a list of application and PETSc node numbers */
272: PetscMalloc4(N, &aomap->app,N,&aomap->appPerm,N,&aomap->petsc,N,&aomap->petscPerm);
273: PetscLogObjectMemory((PetscObject)ao, 4*N * sizeof(PetscInt));
274: for (i = 0; i < N; i++) {
275: appPerm[i] = i;
276: petscPerm[i] = i;
277: }
278: PetscSortIntWithPermutation(N, allpetsc, petscPerm);
279: PetscSortIntWithPermutation(N, allapp, appPerm);
280: /* Form sorted arrays of indices */
281: for (i = 0; i < N; i++) {
282: aomap->app[i] = allapp[appPerm[i]];
283: aomap->petsc[i] = allpetsc[petscPerm[i]];
284: }
285: /* Invert petscPerm[] into aomap->petscPerm[] */
286: for (i = 0; i < N; i++) aomap->petscPerm[petscPerm[i]] = i;
288: /* Form map between aomap->app[] and aomap->petsc[] */
289: for (i = 0; i < N; i++) aomap->appPerm[i] = aomap->petscPerm[appPerm[i]];
291: /* Invert appPerm[] into allapp[] */
292: for (i = 0; i < N; i++) allapp[appPerm[i]] = i;
294: /* Form map between aomap->petsc[] and aomap->app[] */
295: for (i = 0; i < N; i++) aomap->petscPerm[i] = allapp[petscPerm[i]];
297: if (PetscDefined(USE_DEBUG)) {
298: /* Check that the permutations are complementary */
299: for (i = 0; i < N; i++) {
301: }
302: }
303: /* Cleanup */
304: if (!mypetsc) {
305: PetscFree(petsc);
306: }
307: PetscFree4(allapp,appPerm,allpetsc,petscPerm);
309: AOViewFromOptions(ao,NULL,"-ao_view");
311: *aoout = ao;
312: return 0;
313: }
315: /*@
316: AOCreateMappingIS - Creates a basic application ordering using two index sets.
318: Input Parameters:
319: + comm - MPI communicator that is to share AO
320: . isapp - index set that defines an ordering
321: - ispetsc - index set that defines another ordering, maybe NULL for identity IS
323: Output Parameter:
324: . aoout - the new application ordering
326: Options Database Key:
327: . -ao_view - call AOView() at the conclusion of AOCreateMappingIS()
329: Level: beginner
331: Notes:
332: 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.
333: Use AOCreateBasic() or AOCreateBasicIS() if they do not have holes for better performance.
335: .seealso: AOCreateBasic(), AOCreateMapping(), AODestroy()
336: @*/
337: PetscErrorCode AOCreateMappingIS(IS isapp, IS ispetsc, AO *aoout)
338: {
339: MPI_Comm comm;
340: const PetscInt *mypetsc, *myapp;
341: PetscInt napp, npetsc;
343: PetscObjectGetComm((PetscObject) isapp, &comm);
344: ISGetLocalSize(isapp, &napp);
345: if (ispetsc) {
346: ISGetLocalSize(ispetsc, &npetsc);
348: ISGetIndices(ispetsc, &mypetsc);
349: } else {
350: mypetsc = NULL;
351: }
352: ISGetIndices(isapp, &myapp);
354: AOCreateMapping(comm, napp, myapp, mypetsc, aoout);
356: ISRestoreIndices(isapp, &myapp);
357: if (ispetsc) {
358: ISRestoreIndices(ispetsc, &mypetsc);
359: }
360: return 0;
361: }