Actual source code: aomapping.c
petsc-3.7.5 2017-01-01
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> /*I "petscao.h" I*/
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;
20: PetscErrorCode AODestroy_Mapping(AO ao)
21: {
22: AO_Mapping *aomap = (AO_Mapping*) ao->data;
26: PetscFree4(aomap->app,aomap->appPerm,aomap->petsc,aomap->petscPerm);
27: PetscFree(aomap);
28: return(0);
29: }
33: PetscErrorCode AOView_Mapping(AO ao, PetscViewer viewer)
34: {
35: AO_Mapping *aomap = (AO_Mapping*) ao->data;
36: PetscMPIInt rank;
37: PetscInt i;
38: PetscBool iascii;
42: MPI_Comm_rank(PetscObjectComm((PetscObject)ao), &rank);
43: if (rank) return(0);
44: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
45: if (iascii) {
46: PetscViewerASCIIPrintf(viewer, "Number of elements in ordering %D\n", aomap->N);
47: PetscViewerASCIIPrintf(viewer, " App. PETSc\n");
48: for (i = 0; i < aomap->N; i++) {
49: PetscViewerASCIIPrintf(viewer, "%D %D %D\n", i, aomap->app[i], aomap->petsc[aomap->appPerm[i]]);
50: }
51: }
52: return(0);
53: }
57: PetscErrorCode AOPetscToApplication_Mapping(AO ao, PetscInt n, PetscInt *ia)
58: {
59: AO_Mapping *aomap = (AO_Mapping*) ao->data;
60: PetscInt *app = aomap->app;
61: PetscInt *petsc = aomap->petsc;
62: PetscInt *perm = aomap->petscPerm;
63: PetscInt N = aomap->N;
64: PetscInt low, high, mid=0;
65: PetscInt idex;
66: PetscInt i;
68: /* It would be possible to use a single bisection search, which
69: recursively divided the indices to be converted, and searched
70: partitions which contained an index. This would result in
71: better running times if indices are clustered.
72: */
74: for (i = 0; i < n; i++) {
75: idex = ia[i];
76: if (idex < 0) continue;
77: /* Use bisection since the array is sorted */
78: low = 0;
79: high = N - 1;
80: while (low <= high) {
81: mid = (low + high)/2;
82: if (idex == petsc[mid]) break;
83: else if (idex < petsc[mid]) high = mid - 1;
84: else low = mid + 1;
85: }
86: if (low > high) ia[i] = -1;
87: else ia[i] = app[perm[mid]];
88: }
89: return(0);
90: }
94: PetscErrorCode AOApplicationToPetsc_Mapping(AO ao, PetscInt n, PetscInt *ia)
95: {
96: AO_Mapping *aomap = (AO_Mapping*) ao->data;
97: PetscInt *app = aomap->app;
98: PetscInt *petsc = aomap->petsc;
99: PetscInt *perm = aomap->appPerm;
100: PetscInt N = aomap->N;
101: PetscInt low, high, mid=0;
102: PetscInt idex;
103: PetscInt i;
105: /* It would be possible to use a single bisection search, which
106: recursively divided the indices to be converted, and searched
107: partitions which contained an index. This would result in
108: better running times if indices are clustered.
109: */
111: for (i = 0; i < n; i++) {
112: idex = ia[i];
113: if (idex < 0) continue;
114: /* Use bisection since the array is sorted */
115: low = 0;
116: high = N - 1;
117: while (low <= high) {
118: mid = (low + high)/2;
119: if (idex == app[mid]) break;
120: else if (idex < app[mid]) high = mid - 1;
121: else low = mid + 1;
122: }
123: if (low > high) ia[i] = -1;
124: else ia[i] = petsc[perm[mid]];
125: }
126: return(0);
127: }
129: static struct _AOOps AOps = {AOView_Mapping,
130: AODestroy_Mapping,
131: AOPetscToApplication_Mapping,
132: AOApplicationToPetsc_Mapping,
133: NULL,
134: NULL,
135: NULL,
136: NULL};
140: /*@C
141: AOMappingHasApplicationIndex - Searches for the supplied application index.
143: Input Parameters:
144: + ao - The AOMapping
145: - index - The application index
147: Output Parameter:
148: . hasIndex - Flag is PETSC_TRUE if the index exists
150: Level: intermediate
152: .keywords: AO, index
153: .seealso: AOMappingHasPetscIndex(), AOCreateMapping()
154: @*/
155: PetscErrorCode AOMappingHasApplicationIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
156: {
157: AO_Mapping *aomap;
158: PetscInt *app;
159: PetscInt low, high, mid;
164: aomap = (AO_Mapping*) ao->data;
165: app = aomap->app;
166: /* Use bisection since the array is sorted */
167: low = 0;
168: high = aomap->N - 1;
169: while (low <= high) {
170: mid = (low + high)/2;
171: if (idex == app[mid]) break;
172: else if (idex < app[mid]) high = mid - 1;
173: else low = mid + 1;
174: }
175: if (low > high) *hasIndex = PETSC_FALSE;
176: else *hasIndex = PETSC_TRUE;
177: return(0);
178: }
182: /*@
183: AOMappingHasPetscIndex - Searches for the supplied petsc index.
185: Input Parameters:
186: + ao - The AOMapping
187: - index - The petsc index
189: Output Parameter:
190: . hasIndex - Flag is PETSC_TRUE if the index exists
192: Level: intermediate
194: .keywords: AO, index
195: .seealso: AOMappingHasApplicationIndex(), AOCreateMapping()
196: @*/
197: PetscErrorCode AOMappingHasPetscIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
198: {
199: AO_Mapping *aomap;
200: PetscInt *petsc;
201: PetscInt low, high, mid;
206: aomap = (AO_Mapping*) ao->data;
207: petsc = aomap->petsc;
208: /* Use bisection since the array is sorted */
209: low = 0;
210: high = aomap->N - 1;
211: while (low <= high) {
212: mid = (low + high)/2;
213: if (idex == petsc[mid]) break;
214: else if (idex < petsc[mid]) high = mid - 1;
215: else low = mid + 1;
216: }
217: if (low > high) *hasIndex = PETSC_FALSE;
218: else *hasIndex = PETSC_TRUE;
219: return(0);
220: }
224: /*@C
225: AOCreateMapping - Creates a basic application mapping using two integer arrays.
227: Input Parameters:
228: + comm - MPI communicator that is to share AO
229: . napp - size of integer arrays
230: . myapp - integer array that defines an ordering
231: - mypetsc - integer array that defines another ordering (may be NULL to indicate the identity ordering)
233: Output Parameter:
234: . aoout - the new application mapping
236: Options Database Key:
237: . -ao_view : call AOView() at the conclusion of AOCreateMapping()
239: Level: beginner
241: Notes: 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.
242: Use AOCreateBasic() or AOCreateBasicIS() if they do not have holes for better performance.
244: .keywords: AO, create
245: .seealso: AOCreateBasic(), AOCreateBasic(), AOCreateMappingIS(), AODestroy()
246: @*/
247: PetscErrorCode AOCreateMapping(MPI_Comm comm,PetscInt napp,const PetscInt myapp[],const PetscInt mypetsc[],AO *aoout)
248: {
249: AO ao;
250: AO_Mapping *aomap;
251: PetscInt *allpetsc, *allapp;
252: PetscInt *petscPerm, *appPerm;
253: PetscInt *petsc;
254: PetscMPIInt size, rank,*lens, *disp,nnapp;
255: PetscInt N, start;
256: PetscInt i;
261: *aoout = 0;
262: AOInitializePackage();
264: PetscHeaderCreate(ao, AO_CLASSID, "AO", "Application Ordering", "AO", comm, AODestroy, AOView);
265: PetscNewLog(ao,&aomap);
266: PetscMemcpy(ao->ops, &AOps, sizeof(AOps));
267: ao->data = (void*) aomap;
269: /* transmit all lengths to all processors */
270: MPI_Comm_size(comm, &size);
271: MPI_Comm_rank(comm, &rank);
272: PetscMalloc2(size, &lens,size,&disp);
273: nnapp = napp;
274: MPI_Allgather(&nnapp, 1, MPI_INT, lens, 1, MPI_INT, comm);
275: N = 0;
276: for (i = 0; i < size; i++) {
277: disp[i] = N;
278: N += lens[i];
279: }
280: aomap->N = N;
281: ao->N = N;
282: ao->n = N;
284: /* If mypetsc is 0 then use "natural" numbering */
285: if (!mypetsc) {
286: start = disp[rank];
287: PetscMalloc1(napp+1, &petsc);
288: for (i = 0; i < napp; i++) petsc[i] = start + i;
289: } else {
290: petsc = (PetscInt*)mypetsc;
291: }
293: /* get all indices on all processors */
294: PetscMalloc4(N, &allapp,N,&appPerm,N,&allpetsc,N,&petscPerm);
295: MPI_Allgatherv((void*)myapp, napp, MPIU_INT, allapp, lens, disp, MPIU_INT, comm);
296: MPI_Allgatherv((void*)petsc, napp, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm);
297: PetscFree2(lens,disp);
299: /* generate a list of application and PETSc node numbers */
300: PetscMalloc4(N, &aomap->app,N,&aomap->appPerm,N,&aomap->petsc,N,&aomap->petscPerm);
301: PetscLogObjectMemory((PetscObject)ao, 4*N * sizeof(PetscInt));
302: for (i = 0; i < N; i++) {
303: appPerm[i] = i;
304: petscPerm[i] = i;
305: }
306: PetscSortIntWithPermutation(N, allpetsc, petscPerm);
307: PetscSortIntWithPermutation(N, allapp, appPerm);
308: /* Form sorted arrays of indices */
309: for (i = 0; i < N; i++) {
310: aomap->app[i] = allapp[appPerm[i]];
311: aomap->petsc[i] = allpetsc[petscPerm[i]];
312: }
313: /* Invert petscPerm[] into aomap->petscPerm[] */
314: for (i = 0; i < N; i++) aomap->petscPerm[petscPerm[i]] = i;
316: /* Form map between aomap->app[] and aomap->petsc[] */
317: for (i = 0; i < N; i++) aomap->appPerm[i] = aomap->petscPerm[appPerm[i]];
319: /* Invert appPerm[] into allapp[] */
320: for (i = 0; i < N; i++) allapp[appPerm[i]] = i;
322: /* Form map between aomap->petsc[] and aomap->app[] */
323: for (i = 0; i < N; i++) aomap->petscPerm[i] = allapp[petscPerm[i]];
325: #if defined(PETSC_USE_DEBUG)
326: /* Check that the permutations are complementary */
327: for (i = 0; i < N; i++) {
328: if (i != aomap->appPerm[aomap->petscPerm[i]]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid ordering");
329: }
330: #endif
331: /* Cleanup */
332: if (!mypetsc) {
333: PetscFree(petsc);
334: }
335: PetscFree4(allapp,appPerm,allpetsc,petscPerm);
337: AOViewFromOptions(ao,NULL,"-ao_view");
339: *aoout = ao;
340: return(0);
341: }
345: /*@C
346: AOCreateMappingIS - Creates a basic application ordering using two index sets.
348: Input Parameters:
349: + comm - MPI communicator that is to share AO
350: . isapp - index set that defines an ordering
351: - ispetsc - index set that defines another ordering, maybe NULL for identity IS
353: Output Parameter:
354: . aoout - the new application ordering
356: Options Database Key:
357: . -ao_view : call AOView() at the conclusion of AOCreateMappingIS()
359: Level: beginner
361: Notes: 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.
362: Use AOCreateBasic() or AOCreateBasicIS() if they do not have holes for better performance.
364: .keywords: AO, create
365: .seealso: AOCreateBasic(), AOCreateMapping(), AODestroy()
366: @*/
367: PetscErrorCode AOCreateMappingIS(IS isapp, IS ispetsc, AO *aoout)
368: {
369: MPI_Comm comm;
370: const PetscInt *mypetsc, *myapp;
371: PetscInt napp, npetsc;
375: PetscObjectGetComm((PetscObject) isapp, &comm);
376: ISGetLocalSize(isapp, &napp);
377: if (ispetsc) {
378: ISGetLocalSize(ispetsc, &npetsc);
379: if (napp != npetsc) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Local IS lengths must match");
380: ISGetIndices(ispetsc, &mypetsc);
381: } else {
382: mypetsc = NULL;
383: }
384: ISGetIndices(isapp, &myapp);
386: AOCreateMapping(comm, napp, myapp, mypetsc, aoout);
388: ISRestoreIndices(isapp, &myapp);
389: if (ispetsc) {
390: ISRestoreIndices(ispetsc, &mypetsc);
391: }
392: return(0);
393: }