Actual source code: aomapping.c
petsc-3.4.2 2013-07-02
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);
45: if (!viewer) viewer = PETSC_VIEWER_STDOUT_SELF;
47: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
48: if (iascii) {
49: PetscViewerASCIIPrintf(viewer, "Number of elements in ordering %D\n", aomap->N);
50: PetscViewerASCIIPrintf(viewer, " App. PETSc\n");
51: for (i = 0; i < aomap->N; i++) {
52: PetscViewerASCIIPrintf(viewer, "%D %D %D\n", i, aomap->app[i], aomap->petsc[aomap->appPerm[i]]);
53: }
54: }
55: return(0);
56: }
60: PetscErrorCode AOPetscToApplication_Mapping(AO ao, PetscInt n, PetscInt *ia)
61: {
62: AO_Mapping *aomap = (AO_Mapping*) ao->data;
63: PetscInt *app = aomap->app;
64: PetscInt *petsc = aomap->petsc;
65: PetscInt *perm = aomap->petscPerm;
66: PetscInt N = aomap->N;
67: PetscInt low, high, mid=0;
68: PetscInt idex;
69: PetscInt i;
71: /* It would be possible to use a single bisection search, which
72: recursively divided the indices to be converted, and searched
73: partitions which contained an index. This would result in
74: better running times if indices are clustered.
75: */
77: for (i = 0; i < n; i++) {
78: idex = ia[i];
79: if (idex < 0) continue;
80: /* Use bisection since the array is sorted */
81: low = 0;
82: high = N - 1;
83: while (low <= high) {
84: mid = (low + high)/2;
85: if (idex == petsc[mid]) break;
86: else if (idex < petsc[mid]) high = mid - 1;
87: else low = mid + 1;
88: }
89: if (low > high) ia[i] = -1;
90: else ia[i] = app[perm[mid]];
91: }
92: return(0);
93: }
97: PetscErrorCode AOApplicationToPetsc_Mapping(AO ao, PetscInt n, PetscInt *ia)
98: {
99: AO_Mapping *aomap = (AO_Mapping*) ao->data;
100: PetscInt *app = aomap->app;
101: PetscInt *petsc = aomap->petsc;
102: PetscInt *perm = aomap->appPerm;
103: PetscInt N = aomap->N;
104: PetscInt low, high, mid=0;
105: PetscInt idex;
106: PetscInt i;
108: /* It would be possible to use a single bisection search, which
109: recursively divided the indices to be converted, and searched
110: partitions which contained an index. This would result in
111: better running times if indices are clustered.
112: */
114: for (i = 0; i < n; i++) {
115: idex = ia[i];
116: if (idex < 0) continue;
117: /* Use bisection since the array is sorted */
118: low = 0;
119: high = N - 1;
120: while (low <= high) {
121: mid = (low + high)/2;
122: if (idex == app[mid]) break;
123: else if (idex < app[mid]) high = mid - 1;
124: else low = mid + 1;
125: }
126: if (low > high) ia[i] = -1;
127: else ia[i] = petsc[perm[mid]];
128: }
129: return(0);
130: }
132: static struct _AOOps AOps = {AOView_Mapping,
133: AODestroy_Mapping,
134: AOPetscToApplication_Mapping,
135: AOApplicationToPetsc_Mapping,
136: NULL,
137: NULL,
138: NULL,
139: NULL};
143: /*@C
144: AOMappingHasApplicationIndex - Searches for the supplied application index.
146: Input Parameters:
147: + ao - The AOMapping
148: - index - The application index
150: Output Parameter:
151: . hasIndex - Flag is PETSC_TRUE if the index exists
153: Level: intermediate
155: .keywords: AO, index
156: .seealso: AOMappingHasPetscIndex(), AOCreateMapping()
157: @*/
158: PetscErrorCode AOMappingHasApplicationIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
159: {
160: AO_Mapping *aomap;
161: PetscInt *app;
162: PetscInt low, high, mid;
167: aomap = (AO_Mapping*) ao->data;
168: app = aomap->app;
169: /* Use bisection since the array is sorted */
170: low = 0;
171: high = aomap->N - 1;
172: while (low <= high) {
173: mid = (low + high)/2;
174: if (idex == app[mid]) break;
175: else if (idex < app[mid]) high = mid - 1;
176: else low = mid + 1;
177: }
178: if (low > high) *hasIndex = PETSC_FALSE;
179: else *hasIndex = PETSC_TRUE;
180: return(0);
181: }
185: /*@
186: AOMappingHasPetscIndex - Searches for the supplied petsc index.
188: Input Parameters:
189: + ao - The AOMapping
190: - index - The petsc index
192: Output Parameter:
193: . hasIndex - Flag is PETSC_TRUE if the index exists
195: Level: intermediate
197: .keywords: AO, index
198: .seealso: AOMappingHasApplicationIndex(), AOCreateMapping()
199: @*/
200: PetscErrorCode AOMappingHasPetscIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
201: {
202: AO_Mapping *aomap;
203: PetscInt *petsc;
204: PetscInt low, high, mid;
209: aomap = (AO_Mapping*) ao->data;
210: petsc = aomap->petsc;
211: /* Use bisection since the array is sorted */
212: low = 0;
213: high = aomap->N - 1;
214: while (low <= high) {
215: mid = (low + high)/2;
216: if (idex == petsc[mid]) break;
217: else if (idex < petsc[mid]) high = mid - 1;
218: else low = mid + 1;
219: }
220: if (low > high) *hasIndex = PETSC_FALSE;
221: else *hasIndex = PETSC_TRUE;
222: return(0);
223: }
227: /*@C
228: AOCreateMapping - Creates a basic application mapping using two integer arrays.
230: Input Parameters:
231: + comm - MPI communicator that is to share AO
232: . napp - size of integer arrays
233: . myapp - integer array that defines an ordering
234: - mypetsc - integer array that defines another ordering (may be NULL to indicate the identity ordering)
236: Output Parameter:
237: . aoout - the new application mapping
239: Options Database Key:
240: $ -ao_view : call AOView() at the conclusion of AOCreateMapping()
242: Level: beginner
244: 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.
245: Use AOCreateBasic() or AOCreateBasicIS() if they do not have holes for better performance.
247: .keywords: AO, create
248: .seealso: AOCreateBasic(), AOCreateBasic(), AOCreateMappingIS(), AODestroy()
249: @*/
250: PetscErrorCode AOCreateMapping(MPI_Comm comm,PetscInt napp,const PetscInt myapp[],const PetscInt mypetsc[],AO *aoout)
251: {
252: AO ao;
253: AO_Mapping *aomap;
254: PetscInt *allpetsc, *allapp;
255: PetscInt *petscPerm, *appPerm;
256: PetscInt *petsc;
257: PetscMPIInt size, rank,*lens, *disp,nnapp;
258: PetscInt N, start;
259: PetscInt i;
260: PetscBool opt;
265: *aoout = 0;
266: #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
267: AOInitializePackage();
268: #endif
270: PetscHeaderCreate(ao, _p_AO, struct _AOOps, AO_CLASSID, "AO", "Application Ordering", "AO", comm, AODestroy, AOView);
271: PetscNewLog(ao, AO_Mapping, &aomap);
272: PetscMemcpy(ao->ops, &AOps, sizeof(AOps));
273: ao->data = (void*) aomap;
275: /* transmit all lengths to all processors */
276: MPI_Comm_size(comm, &size);
277: MPI_Comm_rank(comm, &rank);
278: PetscMalloc2(size,PetscMPIInt, &lens,size,PetscMPIInt,&disp);
279: nnapp = napp;
280: MPI_Allgather(&nnapp, 1, MPI_INT, lens, 1, MPI_INT, comm);
281: N = 0;
282: for (i = 0; i < size; i++) {
283: disp[i] = N;
284: N += lens[i];
285: }
286: aomap->N = N;
287: ao->N = N;
288: ao->n = N;
290: /* If mypetsc is 0 then use "natural" numbering */
291: if (!mypetsc) {
292: start = disp[rank];
293: PetscMalloc((napp+1) * sizeof(PetscInt), &petsc);
294: for (i = 0; i < napp; i++) petsc[i] = start + i;
295: } else {
296: petsc = (PetscInt*)mypetsc;
297: }
299: /* get all indices on all processors */
300: PetscMalloc4(N,PetscInt, &allapp,N,PetscInt,&appPerm,N,PetscInt,&allpetsc,N,PetscInt,&petscPerm);
301: MPI_Allgatherv((void*)myapp, napp, MPIU_INT, allapp, lens, disp, MPIU_INT, comm);
302: MPI_Allgatherv((void*)petsc, napp, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm);
303: PetscFree2(lens,disp);
305: /* generate a list of application and PETSc node numbers */
306: PetscMalloc4(N,PetscInt, &aomap->app,N,PetscInt,&aomap->appPerm,N,PetscInt,&aomap->petsc,N,PetscInt,&aomap->petscPerm);
307: PetscLogObjectMemory(ao, 4*N * sizeof(PetscInt));
308: for (i = 0; i < N; i++) {
309: appPerm[i] = i;
310: petscPerm[i] = i;
311: }
312: PetscSortIntWithPermutation(N, allpetsc, petscPerm);
313: PetscSortIntWithPermutation(N, allapp, appPerm);
314: /* Form sorted arrays of indices */
315: for (i = 0; i < N; i++) {
316: aomap->app[i] = allapp[appPerm[i]];
317: aomap->petsc[i] = allpetsc[petscPerm[i]];
318: }
319: /* Invert petscPerm[] into aomap->petscPerm[] */
320: for (i = 0; i < N; i++) aomap->petscPerm[petscPerm[i]] = i;
322: /* Form map between aomap->app[] and aomap->petsc[] */
323: for (i = 0; i < N; i++) aomap->appPerm[i] = aomap->petscPerm[appPerm[i]];
325: /* Invert appPerm[] into allapp[] */
326: for (i = 0; i < N; i++) allapp[appPerm[i]] = i;
328: /* Form map between aomap->petsc[] and aomap->app[] */
329: for (i = 0; i < N; i++) aomap->petscPerm[i] = allapp[petscPerm[i]];
331: #if defined(PETSC_USE_DEBUG)
332: /* Check that the permutations are complementary */
333: for (i = 0; i < N; i++) {
334: if (i != aomap->appPerm[aomap->petscPerm[i]]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid ordering");
335: }
336: #endif
337: /* Cleanup */
338: if (!mypetsc) {
339: PetscFree(petsc);
340: }
341: PetscFree4(allapp,appPerm,allpetsc,petscPerm);
343: opt = PETSC_FALSE;
344: PetscOptionsGetBool(NULL, "-ao_view", &opt,NULL);
345: if (opt) {
346: AOView(ao, PETSC_VIEWER_STDOUT_SELF);
347: }
349: *aoout = ao;
350: return(0);
351: }
355: /*@C
356: AOCreateMappingIS - Creates a basic application ordering using two index sets.
358: Input Parameters:
359: + comm - MPI communicator that is to share AO
360: . isapp - index set that defines an ordering
361: - ispetsc - index set that defines another ordering, maybe NULL for identity IS
363: Output Parameter:
364: . aoout - the new application ordering
366: Options Database Key:
367: $ -ao_view : call AOView() at the conclusion of AOCreateMappingIS()
369: Level: beginner
371: 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.
372: Use AOCreateBasic() or AOCreateBasicIS() if they do not have holes for better performance.
374: .keywords: AO, create
375: .seealso: AOCreateBasic(), AOCreateMapping(), AODestroy()
376: @*/
377: PetscErrorCode AOCreateMappingIS(IS isapp, IS ispetsc, AO *aoout)
378: {
379: MPI_Comm comm;
380: const PetscInt *mypetsc, *myapp;
381: PetscInt napp, npetsc;
385: PetscObjectGetComm((PetscObject) isapp, &comm);
386: ISGetLocalSize(isapp, &napp);
387: if (ispetsc) {
388: ISGetLocalSize(ispetsc, &npetsc);
389: if (napp != npetsc) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Local IS lengths must match");
390: ISGetIndices(ispetsc, &mypetsc);
391: } else {
392: mypetsc = NULL;
393: }
394: ISGetIndices(isapp, &myapp);
396: AOCreateMapping(comm, napp, myapp, mypetsc, aoout);
398: ISRestoreIndices(isapp, &myapp);
399: if (ispetsc) {
400: ISRestoreIndices(ispetsc, &mypetsc);
401: }
402: return(0);
403: }