Actual source code: aomapping.c
petsc-3.13.6 2020-09-29
2: /*
3: These AO Section 1.5 Writing Application Codes with PETSc 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;
24: PetscFree4(aomap->app,aomap->appPerm,aomap->petsc,aomap->petscPerm);
25: PetscFree(aomap);
26: return(0);
27: }
29: PetscErrorCode AOView_Mapping(AO ao, PetscViewer viewer)
30: {
31: AO_Mapping *aomap = (AO_Mapping*) ao->data;
32: PetscMPIInt rank;
33: PetscInt i;
34: PetscBool iascii;
38: MPI_Comm_rank(PetscObjectComm((PetscObject)ao), &rank);
39: if (rank) return(0);
40: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
41: if (iascii) {
42: PetscViewerASCIIPrintf(viewer, "Number of elements in ordering %D\n", aomap->N);
43: PetscViewerASCIIPrintf(viewer, " App. PETSc\n");
44: for (i = 0; i < aomap->N; i++) {
45: PetscViewerASCIIPrintf(viewer, "%D %D %D\n", i, aomap->app[i], aomap->petsc[aomap->appPerm[i]]);
46: }
47: }
48: return(0);
49: }
51: PetscErrorCode AOPetscToApplication_Mapping(AO ao, PetscInt n, PetscInt *ia)
52: {
53: AO_Mapping *aomap = (AO_Mapping*) ao->data;
54: PetscInt *app = aomap->app;
55: PetscInt *petsc = aomap->petsc;
56: PetscInt *perm = aomap->petscPerm;
57: PetscInt N = aomap->N;
58: PetscInt low, high, mid=0;
59: PetscInt idex;
60: PetscInt i;
62: /* It would be possible to use a single bisection search, which
63: recursively divided the indices to be converted, and searched
64: partitions which contained an index. This would result in
65: better running times if indices are clustered.
66: */
68: for (i = 0; i < n; i++) {
69: idex = ia[i];
70: if (idex < 0) continue;
71: /* Use bisection since the array is sorted */
72: low = 0;
73: high = N - 1;
74: while (low <= high) {
75: mid = (low + high)/2;
76: if (idex == petsc[mid]) break;
77: else if (idex < petsc[mid]) high = mid - 1;
78: else low = mid + 1;
79: }
80: if (low > high) ia[i] = -1;
81: else ia[i] = app[perm[mid]];
82: }
83: return(0);
84: }
86: PetscErrorCode AOApplicationToPetsc_Mapping(AO ao, PetscInt n, PetscInt *ia)
87: {
88: AO_Mapping *aomap = (AO_Mapping*) ao->data;
89: PetscInt *app = aomap->app;
90: PetscInt *petsc = aomap->petsc;
91: PetscInt *perm = aomap->appPerm;
92: PetscInt N = aomap->N;
93: PetscInt low, high, mid=0;
94: PetscInt idex;
95: PetscInt i;
97: /* It would be possible to use a single bisection search, which
98: recursively divided the indices to be converted, and searched
99: partitions which contained an index. This would result in
100: better running times if indices are clustered.
101: */
103: for (i = 0; i < n; i++) {
104: idex = ia[i];
105: if (idex < 0) continue;
106: /* Use bisection since the array is sorted */
107: low = 0;
108: high = N - 1;
109: while (low <= high) {
110: mid = (low + high)/2;
111: if (idex == app[mid]) break;
112: else if (idex < app[mid]) high = mid - 1;
113: else low = mid + 1;
114: }
115: if (low > high) ia[i] = -1;
116: else ia[i] = petsc[perm[mid]];
117: }
118: return(0);
119: }
121: static struct _AOOps AOps = {AOView_Mapping,
122: AODestroy_Mapping,
123: AOPetscToApplication_Mapping,
124: AOApplicationToPetsc_Mapping,
125: NULL,
126: NULL,
127: NULL,
128: NULL};
130: /*@C
131: AOMappingHasApplicationIndex - Searches for the supplied Section 1.5 Writing Application Codes with PETSc index.
133: Input Parameters:
134: + ao - The AOMapping
135: - index - The Section 1.5 Writing Application Codes with PETSc index
137: Output Parameter:
138: . hasIndex - Flag is PETSC_TRUE if the index exists
140: Level: intermediate
142: .seealso: AOMappingHasPetscIndex(), AOCreateMapping()
143: @*/
144: PetscErrorCode AOMappingHasApplicationIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
145: {
146: AO_Mapping *aomap;
147: PetscInt *app;
148: PetscInt low, high, mid;
153: aomap = (AO_Mapping*) ao->data;
154: app = aomap->app;
155: /* Use bisection since the array is sorted */
156: low = 0;
157: high = aomap->N - 1;
158: while (low <= high) {
159: mid = (low + high)/2;
160: if (idex == app[mid]) break;
161: else if (idex < app[mid]) high = mid - 1;
162: else low = mid + 1;
163: }
164: if (low > high) *hasIndex = PETSC_FALSE;
165: else *hasIndex = PETSC_TRUE;
166: return(0);
167: }
169: /*@
170: AOMappingHasPetscIndex - Searches for the supplied petsc index.
172: Input Parameters:
173: + ao - The AOMapping
174: - index - The petsc index
176: Output Parameter:
177: . hasIndex - Flag is PETSC_TRUE if the index exists
179: Level: intermediate
181: .seealso: AOMappingHasApplicationIndex(), AOCreateMapping()
182: @*/
183: PetscErrorCode AOMappingHasPetscIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
184: {
185: AO_Mapping *aomap;
186: PetscInt *petsc;
187: PetscInt low, high, mid;
192: aomap = (AO_Mapping*) ao->data;
193: petsc = aomap->petsc;
194: /* Use bisection since the array is sorted */
195: low = 0;
196: high = aomap->N - 1;
197: while (low <= high) {
198: mid = (low + high)/2;
199: if (idex == petsc[mid]) break;
200: else if (idex < petsc[mid]) high = mid - 1;
201: else low = mid + 1;
202: }
203: if (low > high) *hasIndex = PETSC_FALSE;
204: else *hasIndex = PETSC_TRUE;
205: return(0);
206: }
208: /*@C
209: AOCreateMapping - Creates a basic Section 1.5 Writing Application Codes with PETSc mapping using two integer arrays.
211: Input Parameters:
212: + comm - MPI communicator that is to share AO
213: . napp - size of integer arrays
214: . myapp - integer array that defines an ordering
215: - mypetsc - integer array that defines another ordering (may be NULL to indicate the identity ordering)
217: Output Parameter:
218: . aoout - the new Section 1.5 Writing Application Codes with PETSc mapping
220: Options Database Key:
221: . -ao_view : call AOView() at the conclusion of AOCreateMapping()
223: Level: beginner
225: Notes:
226: 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.
227: Use AOCreateBasic() or AOCreateBasicIS() if they do not have holes for better performance.
229: .seealso: AOCreateBasic(), AOCreateBasic(), AOCreateMappingIS(), AODestroy()
230: @*/
231: PetscErrorCode AOCreateMapping(MPI_Comm comm,PetscInt napp,const PetscInt myapp[],const PetscInt mypetsc[],AO *aoout)
232: {
233: AO ao;
234: AO_Mapping *aomap;
235: PetscInt *allpetsc, *allapp;
236: PetscInt *petscPerm, *appPerm;
237: PetscInt *petsc;
238: PetscMPIInt size, rank,*lens, *disp,nnapp;
239: PetscInt N, start;
240: PetscInt i;
245: *aoout = 0;
246: AOInitializePackage();
248: PetscHeaderCreate(ao, AO_CLASSID, "AO", "Application Ordering", "AO", comm, AODestroy, AOView);
249: PetscNewLog(ao,&aomap);
250: PetscMemcpy(ao->ops, &AOps, sizeof(AOps));
251: ao->data = (void*) aomap;
253: /* transmit all lengths to all processors */
254: MPI_Comm_size(comm, &size);
255: MPI_Comm_rank(comm, &rank);
256: PetscMalloc2(size, &lens,size,&disp);
257: nnapp = napp;
258: MPI_Allgather(&nnapp, 1, MPI_INT, lens, 1, MPI_INT, comm);
259: N = 0;
260: for (i = 0; i < size; i++) {
261: disp[i] = N;
262: N += lens[i];
263: }
264: aomap->N = N;
265: ao->N = N;
266: ao->n = N;
268: /* If mypetsc is 0 then use "natural" numbering */
269: if (!mypetsc) {
270: start = disp[rank];
271: PetscMalloc1(napp+1, &petsc);
272: for (i = 0; i < napp; i++) petsc[i] = start + i;
273: } else {
274: petsc = (PetscInt*)mypetsc;
275: }
277: /* get all indices on all processors */
278: PetscMalloc4(N, &allapp,N,&appPerm,N,&allpetsc,N,&petscPerm);
279: MPI_Allgatherv((void*)myapp, napp, MPIU_INT, allapp, lens, disp, MPIU_INT, comm);
280: MPI_Allgatherv((void*)petsc, napp, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm);
281: PetscFree2(lens,disp);
283: /* generate a list of Section 1.5 Writing Application Codes with PETSc and PETSc node numbers */
284: PetscMalloc4(N, &aomap->app,N,&aomap->appPerm,N,&aomap->petsc,N,&aomap->petscPerm);
285: PetscLogObjectMemory((PetscObject)ao, 4*N * sizeof(PetscInt));
286: for (i = 0; i < N; i++) {
287: appPerm[i] = i;
288: petscPerm[i] = i;
289: }
290: PetscSortIntWithPermutation(N, allpetsc, petscPerm);
291: PetscSortIntWithPermutation(N, allapp, appPerm);
292: /* Form sorted arrays of indices */
293: for (i = 0; i < N; i++) {
294: aomap->app[i] = allapp[appPerm[i]];
295: aomap->petsc[i] = allpetsc[petscPerm[i]];
296: }
297: /* Invert petscPerm[] into aomap->petscPerm[] */
298: for (i = 0; i < N; i++) aomap->petscPerm[petscPerm[i]] = i;
300: /* Form map between aomap->app[] and aomap->petsc[] */
301: for (i = 0; i < N; i++) aomap->appPerm[i] = aomap->petscPerm[appPerm[i]];
303: /* Invert appPerm[] into allapp[] */
304: for (i = 0; i < N; i++) allapp[appPerm[i]] = i;
306: /* Form map between aomap->petsc[] and aomap->app[] */
307: for (i = 0; i < N; i++) aomap->petscPerm[i] = allapp[petscPerm[i]];
309: #if defined(PETSC_USE_DEBUG)
310: /* Check that the permutations are complementary */
311: for (i = 0; i < N; i++) {
312: if (i != aomap->appPerm[aomap->petscPerm[i]]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid ordering");
313: }
314: #endif
315: /* Cleanup */
316: if (!mypetsc) {
317: PetscFree(petsc);
318: }
319: PetscFree4(allapp,appPerm,allpetsc,petscPerm);
321: AOViewFromOptions(ao,NULL,"-ao_view");
323: *aoout = ao;
324: return(0);
325: }
327: /*@
328: AOCreateMappingIS - Creates a basic Section 1.5 Writing Application Codes with PETSc ordering using two index sets.
330: Input Parameters:
331: + comm - MPI communicator that is to share AO
332: . isapp - index set that defines an ordering
333: - ispetsc - index set that defines another ordering, maybe NULL for identity IS
335: Output Parameter:
336: . aoout - the new Section 1.5 Writing Application Codes with PETSc ordering
338: Options Database Key:
339: . -ao_view : call AOView() at the conclusion of AOCreateMappingIS()
341: Level: beginner
343: Notes:
344: 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.
345: Use AOCreateBasic() or AOCreateBasicIS() if they do not have holes for better performance.
347: .seealso: AOCreateBasic(), AOCreateMapping(), AODestroy()
348: @*/
349: PetscErrorCode AOCreateMappingIS(IS isapp, IS ispetsc, AO *aoout)
350: {
351: MPI_Comm comm;
352: const PetscInt *mypetsc, *myapp;
353: PetscInt napp, npetsc;
357: PetscObjectGetComm((PetscObject) isapp, &comm);
358: ISGetLocalSize(isapp, &napp);
359: if (ispetsc) {
360: ISGetLocalSize(ispetsc, &npetsc);
361: if (napp != npetsc) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Local IS lengths must match");
362: ISGetIndices(ispetsc, &mypetsc);
363: } else {
364: mypetsc = NULL;
365: }
366: ISGetIndices(isapp, &myapp);
368: AOCreateMapping(comm, napp, myapp, mypetsc, aoout);
370: ISRestoreIndices(isapp, &myapp);
371: if (ispetsc) {
372: ISRestoreIndices(ispetsc, &mypetsc);
373: }
374: return(0);
375: }