Actual source code: aomapping.c
petsc-3.3-p7 2013-05-11
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/dm/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(((PetscObject)ao)->comm, &rank);
43: if (rank) return(0);
45: if (!viewer) {
46: viewer = PETSC_VIEWER_STDOUT_SELF;
47: }
49: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
50: if (iascii) {
51: PetscViewerASCIIPrintf(viewer, "Number of elements in ordering %D\n", aomap->N);
52: PetscViewerASCIIPrintf(viewer, " App. PETSc\n");
53: for(i = 0; i < aomap->N; i++) {
54: PetscViewerASCIIPrintf(viewer, "%D %D %D\n", i, aomap->app[i], aomap->petsc[aomap->appPerm[i]]);
55: }
56: }
57: return(0);
58: }
62: PetscErrorCode AOPetscToApplication_Mapping(AO ao, PetscInt n, PetscInt *ia)
63: {
64: AO_Mapping *aomap = (AO_Mapping *) ao->data;
65: PetscInt *app = aomap->app;
66: PetscInt *petsc = aomap->petsc;
67: PetscInt *perm = aomap->petscPerm;
68: PetscInt N = aomap->N;
69: PetscInt low, high, mid=0;
70: PetscInt idex;
71: PetscInt i;
73: /* It would be possible to use a single bisection search, which
74: recursively divided the indices to be converted, and searched
75: partitions which contained an index. This would result in
76: better running times if indices are clustered.
77: */
79: for(i = 0; i < n; i++) {
80: idex = ia[i];
81: if (idex < 0) continue;
82: /* Use bisection since the array is sorted */
83: low = 0;
84: high = N - 1;
85: while (low <= high) {
86: mid = (low + high)/2;
87: if (idex == petsc[mid]) {
88: break;
89: } else if (idex < petsc[mid]) {
90: high = mid - 1;
91: } else {
92: low = mid + 1;
93: }
94: }
95: if (low > high) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Invalid input index %D", idex);
96: ia[i] = app[perm[mid]];
97: }
98: return(0);
99: }
103: PetscErrorCode AOApplicationToPetsc_Mapping(AO ao, PetscInt n, PetscInt *ia)
104: {
105: AO_Mapping *aomap = (AO_Mapping *) ao->data;
106: PetscInt *app = aomap->app;
107: PetscInt *petsc = aomap->petsc;
108: PetscInt *perm = aomap->appPerm;
109: PetscInt N = aomap->N;
110: PetscInt low, high, mid=0;
111: PetscInt idex;
112: PetscInt i;
114: /* It would be possible to use a single bisection search, which
115: recursively divided the indices to be converted, and searched
116: partitions which contained an index. This would result in
117: better running times if indices are clustered.
118: */
120: for(i = 0; i < n; i++) {
121: idex = ia[i];
122: if (idex < 0) continue;
123: /* Use bisection since the array is sorted */
124: low = 0;
125: high = N - 1;
126: while (low <= high) {
127: mid = (low + high)/2;
128: if (idex == app[mid]) {
129: break;
130: } else if (idex < app[mid]) {
131: high = mid - 1;
132: } else {
133: low = mid + 1;
134: }
135: }
136: if (low > high) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Invalid input index %D", idex);
137: ia[i] = petsc[perm[mid]];
138: }
139: return(0);
140: }
142: static struct _AOOps AOps = {AOView_Mapping,
143: AODestroy_Mapping,
144: AOPetscToApplication_Mapping,
145: AOApplicationToPetsc_Mapping,
146: PETSC_NULL,
147: PETSC_NULL,
148: PETSC_NULL,
149: PETSC_NULL};
153: /*@C
154: AOMappingHasApplicationIndex - Searches for the supplied application index.
156: Input Parameters:
157: + ao - The AOMapping
158: - index - The application index
160: Output Parameter:
161: . hasIndex - Flag is PETSC_TRUE if the index exists
163: Level: intermediate
165: .keywords: AO, index
166: .seealso: AOMappingHasPetscIndex(), AOCreateMapping()
167: @*/
168: PetscErrorCode AOMappingHasApplicationIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
169: {
170: AO_Mapping *aomap;
171: PetscInt *app;
172: PetscInt low, high, mid;
177: aomap = (AO_Mapping *) ao->data;
178: app = aomap->app;
179: /* Use bisection since the array is sorted */
180: low = 0;
181: high = aomap->N - 1;
182: while (low <= high) {
183: mid = (low + high)/2;
184: if (idex == app[mid]) {
185: break;
186: } else if (idex < app[mid]) {
187: high = mid - 1;
188: } else {
189: low = mid + 1;
190: }
191: }
192: if (low > high) {
193: *hasIndex = PETSC_FALSE;
194: } else {
195: *hasIndex = PETSC_TRUE;
196: }
197: return(0);
198: }
202: /*@
203: AOMappingHasPetscIndex - Searches for the supplied petsc index.
205: Input Parameters:
206: + ao - The AOMapping
207: - index - The petsc index
209: Output Parameter:
210: . hasIndex - Flag is PETSC_TRUE if the index exists
212: Level: intermediate
214: .keywords: AO, index
215: .seealso: AOMappingHasApplicationIndex(), AOCreateMapping()
216: @*/
217: PetscErrorCode AOMappingHasPetscIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
218: {
219: AO_Mapping *aomap;
220: PetscInt *petsc;
221: PetscInt low, high, mid;
226: aomap = (AO_Mapping *) ao->data;
227: petsc = aomap->petsc;
228: /* Use bisection since the array is sorted */
229: low = 0;
230: high = aomap->N - 1;
231: while (low <= high) {
232: mid = (low + high)/2;
233: if (idex == petsc[mid]) {
234: break;
235: } else if (idex < petsc[mid]) {
236: high = mid - 1;
237: } else {
238: low = mid + 1;
239: }
240: }
241: if (low > high) {
242: *hasIndex = PETSC_FALSE;
243: } else {
244: *hasIndex = PETSC_TRUE;
245: }
246: return(0);
247: }
251: /*@C
252: AOCreateMapping - Creates a basic application mapping using two integer arrays.
254: Input Parameters:
255: + comm - MPI communicator that is to share AO
256: . napp - size of integer arrays
257: . myapp - integer array that defines an ordering
258: - mypetsc - integer array that defines another ordering (may be PETSC_NULL to indicate the identity ordering)
260: Output Parameter:
261: . aoout - the new application mapping
263: Options Database Key:
264: $ -ao_view : call AOView() at the conclusion of AOCreateMapping()
266: Level: beginner
268: 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.
269: Use AOCreateBasic() or AOCreateBasicIS() if they do not have holes for better performance.
271: .keywords: AO, create
272: .seealso: AOCreateBasic(), AOCreateBasic(), AOCreateMappingIS(), AODestroy()
273: @*/
274: PetscErrorCode AOCreateMapping(MPI_Comm comm,PetscInt napp,const PetscInt myapp[],const PetscInt mypetsc[],AO *aoout)
275: {
276: AO ao;
277: AO_Mapping *aomap;
278: PetscInt *allpetsc, *allapp;
279: PetscInt *petscPerm, *appPerm;
280: PetscInt *petsc;
281: PetscMPIInt size, rank,*lens, *disp,nnapp;
282: PetscInt N, start;
283: PetscInt i;
284: PetscBool opt;
289: *aoout = 0;
290: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
291: AOInitializePackage(PETSC_NULL);
292: #endif
294: PetscHeaderCreate(ao, _p_AO, struct _AOOps, AO_CLASSID, -1, "AO", "Application Ordering", "AO", comm, AODestroy, AOView);
295: PetscNewLog(ao, AO_Mapping, &aomap);
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: PetscMalloc2(size,PetscMPIInt, &lens,size,PetscMPIInt,&disp);
303: nnapp = napp;
304: MPI_Allgather(&nnapp, 1, MPI_INT, lens, 1, MPI_INT, comm);
305: N = 0;
306: for(i = 0; i < size; i++) {
307: disp[i] = N;
308: N += lens[i];
309: }
310: aomap->N = N;
311: ao->N = N;
312: ao->n = N;
314: /* If mypetsc is 0 then use "natural" numbering */
315: if (!mypetsc) {
316: start = disp[rank];
317: PetscMalloc((napp+1) * sizeof(PetscInt), &petsc);
318: for(i = 0; i < napp; i++) {
319: petsc[i] = start + i;
320: }
321: } else {
322: petsc = (PetscInt*)mypetsc;
323: }
325: /* get all indices on all processors */
326: PetscMalloc4(N,PetscInt, &allapp,N,PetscInt,&appPerm,N,PetscInt,&allpetsc,N,PetscInt,&petscPerm);
327: MPI_Allgatherv((void*)myapp, napp, MPIU_INT, allapp, lens, disp, MPIU_INT, comm);
328: MPI_Allgatherv((void*)petsc, napp, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm);
329: PetscFree2(lens,disp);
331: /* generate a list of application and PETSc node numbers */
332: PetscMalloc4(N,PetscInt, &aomap->app,N,PetscInt,&aomap->appPerm,N,PetscInt,&aomap->petsc,N,PetscInt,&aomap->petscPerm);
333: PetscLogObjectMemory(ao, 4*N * sizeof(PetscInt));
334: for(i = 0; i < N; i++) {
335: appPerm[i] = i;
336: petscPerm[i] = i;
337: }
338: PetscSortIntWithPermutation(N, allpetsc, petscPerm);
339: PetscSortIntWithPermutation(N, allapp, appPerm);
340: /* Form sorted arrays of indices */
341: for(i = 0; i < N; i++) {
342: aomap->app[i] = allapp[appPerm[i]];
343: aomap->petsc[i] = allpetsc[petscPerm[i]];
344: }
345: /* Invert petscPerm[] into aomap->petscPerm[] */
346: for(i = 0; i < N; i++) {
347: aomap->petscPerm[petscPerm[i]] = i;
348: }
349: /* Form map between aomap->app[] and aomap->petsc[] */
350: for(i = 0; i < N; i++) {
351: aomap->appPerm[i] = aomap->petscPerm[appPerm[i]];
352: }
353: /* Invert appPerm[] into allapp[] */
354: for(i = 0; i < N; i++) {
355: allapp[appPerm[i]] = i;
356: }
357: /* Form map between aomap->petsc[] and aomap->app[] */
358: for(i = 0; i < N; i++) {
359: aomap->petscPerm[i] = allapp[petscPerm[i]];
360: }
361: #ifdef PETSC_USE_DEBUG
362: /* Check that the permutations are complementary */
363: for(i = 0; i < N; i++) {
364: if (i != aomap->appPerm[aomap->petscPerm[i]]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid ordering");
365: }
366: #endif
367: /* Cleanup */
368: if (!mypetsc) {
369: PetscFree(petsc);
370: }
371: PetscFree4(allapp,appPerm,allpetsc,petscPerm);
373: opt = PETSC_FALSE;
374: PetscOptionsGetBool(PETSC_NULL, "-ao_view", &opt,PETSC_NULL);
375: if (opt) {
376: AOView(ao, PETSC_VIEWER_STDOUT_SELF);
377: }
379: *aoout = ao;
380: return(0);
381: }
385: /*@C
386: AOCreateMappingIS - Creates a basic application ordering using two index sets.
388: Input Parameters:
389: + comm - MPI communicator that is to share AO
390: . isapp - index set that defines an ordering
391: - ispetsc - index set that defines another ordering, maybe PETSC_NULL for identity IS
393: Output Parameter:
394: . aoout - the new application ordering
396: Options Database Key:
397: $ -ao_view : call AOView() at the conclusion of AOCreateMappingIS()
399: Level: beginner
401: 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.
402: Use AOCreateBasic() or AOCreateBasicIS() if they do not have holes for better performance.
404: .keywords: AO, create
405: .seealso: AOCreateBasic(), AOCreateMapping(), AODestroy()
406: @*/
407: PetscErrorCode AOCreateMappingIS(IS isapp, IS ispetsc, AO *aoout)
408: {
409: MPI_Comm comm;
410: const PetscInt *mypetsc, *myapp;
411: PetscInt napp, npetsc;
412: PetscErrorCode ierr;
415: PetscObjectGetComm((PetscObject) isapp, &comm);
416: ISGetLocalSize(isapp, &napp);
417: if (ispetsc) {
418: ISGetLocalSize(ispetsc, &npetsc);
419: if (napp != npetsc) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Local IS lengths must match");
420: ISGetIndices(ispetsc, &mypetsc);
421: } else {
422: mypetsc = NULL;
423: }
424: ISGetIndices(isapp, &myapp);
426: AOCreateMapping(comm, napp, myapp, mypetsc, aoout);
428: ISRestoreIndices(isapp, &myapp);
429: if (ispetsc) {
430: ISRestoreIndices(ispetsc, &mypetsc);
431: }
432: return(0);
433: }