Actual source code: aomapping.c
1: #define PETSCDM_DLL
3: /*
4: These AO application ordering routines do not require that the input
5: be a permutation, but merely a 1-1 mapping. This implementation still
6: keeps the entire ordering on each processor.
7: */
9: #include src/dm/ao/aoimpl.h
10: #include petscsys.h
12: typedef struct {
13: PetscInt N;
14: PetscInt *app; /* app[i] is the partner for petsc[appPerm[i]] */
15: PetscInt *appPerm;
16: PetscInt *petsc; /* petsc[j] is the partner for app[petscPerm[j]] */
17: PetscInt *petscPerm;
18: } AO_Mapping;
22: PetscErrorCode AODestroy_Mapping(AO ao)
23: {
24: AO_Mapping *aomap = (AO_Mapping *) ao->data;
28: PetscFree(aomap->app);
29: PetscFree(ao->data);
30: return(0);
31: }
35: PetscErrorCode AOView_Mapping(AO ao, PetscViewer viewer)
36: {
37: AO_Mapping *aomap = (AO_Mapping *) ao->data;
38: PetscMPIInt rank;
39: PetscInt i;
40: PetscTruth iascii;
44: MPI_Comm_rank(ao->comm, &rank);
45: if (rank) return(0);
47: if (!viewer) {
48: viewer = PETSC_VIEWER_STDOUT_SELF;
49: }
51: PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_ASCII, &iascii);
52: if (iascii) {
53: PetscViewerASCIIPrintf(viewer, "Number of elements in ordering %D\n", aomap->N);
54: PetscViewerASCIIPrintf(viewer, " App. PETSc\n");
55: for(i = 0; i < aomap->N; i++) {
56: PetscViewerASCIIPrintf(viewer, "%D %D %D\n", i, aomap->app[i], aomap->petsc[aomap->appPerm[i]]);
57: }
58: }
59: return(0);
60: }
64: PetscErrorCode AOPetscToApplication_Mapping(AO ao, PetscInt n, PetscInt *ia)
65: {
66: AO_Mapping *aomap = (AO_Mapping *) ao->data;
67: PetscInt *app = aomap->app;
68: PetscInt *petsc = aomap->petsc;
69: PetscInt *perm = aomap->petscPerm;
70: PetscInt N = aomap->N;
71: PetscInt low, high, mid=0;
72: PetscInt idex;
73: PetscInt i;
75: /* It would be possible to use a single bisection search, which
76: recursively divided the indices to be converted, and searched
77: partitions which contained an index. This would result in
78: better running times if indices are clustered.
79: */
81: for(i = 0; i < n; i++) {
82: idex = ia[i];
83: if (idex < 0) continue;
84: /* Use bisection since the array is sorted */
85: low = 0;
86: high = N - 1;
87: while (low <= high) {
88: mid = (low + high)/2;
89: if (idex == petsc[mid]) {
90: break;
91: } else if (idex < petsc[mid]) {
92: high = mid - 1;
93: } else {
94: low = mid + 1;
95: }
96: }
97: if (low > high) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Invalid input index %D", idex);
98: ia[i] = app[perm[mid]];
99: }
100: return(0);
101: }
105: PetscErrorCode AOApplicationToPetsc_Mapping(AO ao, PetscInt n, PetscInt *ia)
106: {
107: AO_Mapping *aomap = (AO_Mapping *) ao->data;
108: PetscInt *app = aomap->app;
109: PetscInt *petsc = aomap->petsc;
110: PetscInt *perm = aomap->appPerm;
111: PetscInt N = aomap->N;
112: PetscInt low, high, mid=0;
113: PetscInt idex;
114: PetscInt i;
116: /* It would be possible to use a single bisection search, which
117: recursively divided the indices to be converted, and searched
118: partitions which contained an index. This would result in
119: better running times if indices are clustered.
120: */
122: for(i = 0; i < n; i++) {
123: idex = ia[i];
124: if (idex < 0) continue;
125: /* Use bisection since the array is sorted */
126: low = 0;
127: high = N - 1;
128: while (low <= high) {
129: mid = (low + high)/2;
130: if (idex == app[mid]) {
131: break;
132: } else if (idex < app[mid]) {
133: high = mid - 1;
134: } else {
135: low = mid + 1;
136: }
137: }
138: if (low > high) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Invalid input index %D", idex);
139: ia[i] = petsc[perm[mid]];
140: }
141: return(0);
142: }
144: static struct _AOOps AOps = {AOView_Mapping,
145: AODestroy_Mapping,
146: AOPetscToApplication_Mapping,
147: AOApplicationToPetsc_Mapping,
148: PETSC_NULL,
149: PETSC_NULL,
150: PETSC_NULL,
151: PETSC_NULL};
155: /*@C
156: AOMappingHasApplicationIndex - Searches for the supplied application index.
158: Input Parameters:
159: + ao - The AOMapping
160: - index - The application index
162: Output Parameter:
163: . hasIndex - Flag is PETSC_TRUE if the index exists
165: Level: intermediate
167: .keywords: AO, index
168: .seealso: AOMappingHasPetscIndex(), AOCreateMapping()
169: @*/
170: PetscErrorCode AOMappingHasApplicationIndex(AO ao, PetscInt idex, PetscTruth *hasIndex)
171: {
172: AO_Mapping *aomap;
173: PetscInt *app;
174: PetscInt low, high, mid;
179: aomap = (AO_Mapping *) ao->data;
180: app = aomap->app;
181: /* Use bisection since the array is sorted */
182: low = 0;
183: high = aomap->N - 1;
184: while (low <= high) {
185: mid = (low + high)/2;
186: if (idex == app[mid]) {
187: break;
188: } else if (idex < app[mid]) {
189: high = mid - 1;
190: } else {
191: low = mid + 1;
192: }
193: }
194: if (low > high) {
195: *hasIndex = PETSC_FALSE;
196: } else {
197: *hasIndex = PETSC_TRUE;
198: }
199: return(0);
200: }
204: /*@C
205: AOMappingHasPetscIndex - Searches for the supplied petsc index.
207: Input Parameters:
208: + ao - The AOMapping
209: - index - The petsc index
211: Output Parameter:
212: . hasIndex - Flag is PETSC_TRUE if the index exists
214: Level: intermediate
216: .keywords: AO, index
217: .seealso: AOMappingHasApplicationIndex(), AOCreateMapping()
218: @*/
219: PetscErrorCode AOMappingHasPetscIndex(AO ao, PetscInt idex, PetscTruth *hasIndex)
220: {
221: AO_Mapping *aomap;
222: PetscInt *petsc;
223: PetscInt low, high, mid;
228: aomap = (AO_Mapping *) ao->data;
229: petsc = aomap->petsc;
230: /* Use bisection since the array is sorted */
231: low = 0;
232: high = aomap->N - 1;
233: while (low <= high) {
234: mid = (low + high)/2;
235: if (idex == petsc[mid]) {
236: break;
237: } else if (idex < petsc[mid]) {
238: high = mid - 1;
239: } else {
240: low = mid + 1;
241: }
242: }
243: if (low > high) {
244: *hasIndex = PETSC_FALSE;
245: } else {
246: *hasIndex = PETSC_TRUE;
247: }
248: return(0);
249: }
253: /*@C
254: AOCreateMapping - Creates a basic application mapping using two integer arrays.
256: Input Parameters:
257: + comm - MPI communicator that is to share AO
258: . napp - size of integer arrays
259: . myapp - integer array that defines an ordering
260: - mypetsc - integer array that defines another ordering
262: Output Parameter:
263: . aoout - the new application mapping
265: Options Database Key:
266: $ -ao_view : call AOView() at the conclusion of AOCreateMapping()
268: Level: beginner
270: .keywords: AO, create
271: .seealso: AOCreateDebug(), AOCreateBasic(), AOCreateMappingIS(), AODestroy()
272: @*/
273: PetscErrorCode AOCreateMapping(MPI_Comm comm,PetscInt napp,const PetscInt myapp[],const PetscInt mypetsc[],AO *aoout)
274: {
275: AO ao;
276: AO_Mapping *aomap;
277: PetscInt *allpetsc, *allapp;
278: PetscInt *petscPerm, *appPerm;
279: PetscInt *petsc;
280: PetscMPIInt size, rank,*lens, *disp,nnapp;
281: PetscInt N, start;
282: PetscInt i;
283: PetscTruth opt;
288: *aoout = 0;
289: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
290: DMInitializePackage(PETSC_NULL);
291: #endif
293: PetscHeaderCreate(ao, _p_AO, struct _AOOps, AO_COOKIE, AO_MAPPING, "AO", comm, AODestroy, AOView);
294: PetscNew(AO_Mapping, &aomap);
295: PetscLogObjectMemory(ao, sizeof(struct _p_AO) + sizeof(AO_Mapping));
296: PetscMemcpy(ao->ops, &AOps, sizeof(AOps));
297: ao->data = (void*) aomap;
299: /* transmit all lengths to all processors */
300: MPI_Comm_size(comm, &size);
301: MPI_Comm_rank(comm, &rank);
302: PetscMalloc(2*size * sizeof(PetscMPIInt), &lens);
303: disp = lens + size;
304: nnapp = napp;
305: MPI_Allgather(&nnapp, 1, MPI_INT, lens, 1, MPI_INT, comm);
306: N = 0;
307: for(i = 0; i < size; i++) {
308: disp[i] = N;
309: N += lens[i];
310: }
311: aomap->N = N;
312: ao->N = N;
313: ao->n = N;
315: /* If mypetsc is 0 then use "natural" numbering */
316: if (!mypetsc) {
317: start = disp[rank];
318: PetscMalloc((napp+1) * sizeof(PetscInt), &petsc);
319: for(i = 0; i < napp; i++) {
320: petsc[i] = start + i;
321: }
322: } else {
323: petsc = (PetscInt*)mypetsc;
324: }
326: /* get all indices on all processors */
327: PetscMalloc(N*4 * sizeof(PetscInt), &allapp);
328: appPerm = allapp + N;
329: allpetsc = appPerm + N;
330: petscPerm = allpetsc + N;
331: MPI_Allgatherv((void*)myapp, napp, MPIU_INT, allapp, lens, disp, MPIU_INT, comm);
332: MPI_Allgatherv((void*)petsc, napp, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm);
333: PetscFree(lens);
335: /* generate a list of application and PETSc node numbers */
336: PetscMalloc(N*4 * sizeof(PetscInt), &aomap->app);
337: PetscLogObjectMemory(ao, 4*N * sizeof(PetscInt));
338: aomap->appPerm = aomap->app + N;
339: aomap->petsc = aomap->appPerm + N;
340: aomap->petscPerm = aomap->petsc + N;
341: for(i = 0; i < N; i++) {
342: appPerm[i] = i;
343: petscPerm[i] = i;
344: }
345: PetscSortIntWithPermutation(N, allpetsc, petscPerm);
346: PetscSortIntWithPermutation(N, allapp, appPerm);
347: /* Form sorted arrays of indices */
348: for(i = 0; i < N; i++) {
349: aomap->app[i] = allapp[appPerm[i]];
350: aomap->petsc[i] = allpetsc[petscPerm[i]];
351: }
352: /* Invert petscPerm[] into aomap->petscPerm[] */
353: for(i = 0; i < N; i++) {
354: aomap->petscPerm[petscPerm[i]] = i;
355: }
356: /* Form map between aomap->app[] and aomap->petsc[] */
357: for(i = 0; i < N; i++) {
358: aomap->appPerm[i] = aomap->petscPerm[appPerm[i]];
359: }
360: /* Invert appPerm[] into allapp[] */
361: for(i = 0; i < N; i++) {
362: allapp[appPerm[i]] = i;
363: }
364: /* Form map between aomap->petsc[] and aomap->app[] */
365: for(i = 0; i < N; i++) {
366: aomap->petscPerm[i] = allapp[petscPerm[i]];
367: }
368: #ifdef PETSC_USE_DEBUG
369: /* Check that the permutations are complementary */
370: for(i = 0; i < N; i++) {
371: if (i != aomap->appPerm[aomap->petscPerm[i]])
372: SETERRQ(PETSC_ERR_PLIB, "Invalid ordering");
373: }
374: #endif
375: /* Cleanup */
376: if (!mypetsc) {
377: PetscFree(petsc);
378: }
379: PetscFree(allapp);
381: PetscOptionsHasName(PETSC_NULL, "-ao_view", &opt);
382: if (opt) {
383: AOView(ao, PETSC_VIEWER_STDOUT_SELF);
384: }
386: *aoout = ao;
387: return(0);
388: }
392: /*@C
393: AOCreateMappingIS - Creates a basic application ordering using two index sets.
395: Input Parameters:
396: + comm - MPI communicator that is to share AO
397: . isapp - index set that defines an ordering
398: - ispetsc - index set that defines another ordering
400: Output Parameter:
401: . aoout - the new application ordering
403: Options Database Key:
404: $ -ao_view : call AOView() at the conclusion of AOCreateMappingIS()
406: Level: beginner
408: .keywords: AO, create
409: .seealso: AOCreateBasic(), AOCreateMapping(), AODestroy()
410: @*/
411: PetscErrorCode AOCreateMappingIS(IS isapp, IS ispetsc, AO *aoout)
412: {
413: MPI_Comm comm;
414: PetscInt *mypetsc, *myapp;
415: PetscInt napp, npetsc;
419: PetscObjectGetComm((PetscObject) isapp, &comm);
420: ISGetLocalSize(isapp, &napp);
421: if (ispetsc) {
422: ISGetLocalSize(ispetsc, &npetsc);
423: if (napp != npetsc) SETERRQ(PETSC_ERR_ARG_SIZ, "Local IS lengths must match");
424: ISGetIndices(ispetsc, &mypetsc);
425: } else {
426: mypetsc = NULL;
427: }
428: ISGetIndices(isapp, &myapp);
430: AOCreateMapping(comm, napp, myapp, mypetsc, aoout);
432: ISRestoreIndices(isapp, &myapp);
433: if (ispetsc) {
434: ISRestoreIndices(ispetsc, &mypetsc);
435: }
436: return(0);
437: }