Actual source code: aobasic.c
2: /*
3: The most basic AO application ordering routines. These store the
4: entire orderings on each processor.
5: */
7: #include <../src/vec/is/ao/aoimpl.h>
9: typedef struct {
10: PetscInt *app; /* app[i] is the partner for the ith PETSc slot */
11: PetscInt *petsc; /* petsc[j] is the partner for the jth app slot */
12: } AO_Basic;
14: /*
15: All processors have the same data so processor 1 prints it
16: */
17: PetscErrorCode AOView_Basic(AO ao, PetscViewer viewer)
18: {
19: PetscMPIInt rank;
20: PetscInt i;
21: AO_Basic *aobasic = (AO_Basic *)ao->data;
22: PetscBool iascii;
24: MPI_Comm_rank(PetscObjectComm((PetscObject)ao), &rank);
25: if (rank == 0) {
26: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
27: if (iascii) {
28: PetscViewerASCIIPrintf(viewer, "Number of elements in ordering %" PetscInt_FMT "\n", ao->N);
29: PetscViewerASCIIPrintf(viewer, "PETSc->App App->PETSc\n");
30: for (i = 0; i < ao->N; i++) PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT " %3" PetscInt_FMT " %3" PetscInt_FMT " %3" PetscInt_FMT "\n", i, aobasic->app[i], i, aobasic->petsc[i]);
31: }
32: }
33: PetscViewerFlush(viewer);
34: return 0;
35: }
37: PetscErrorCode AODestroy_Basic(AO ao)
38: {
39: AO_Basic *aobasic = (AO_Basic *)ao->data;
41: PetscFree2(aobasic->app, aobasic->petsc);
42: PetscFree(aobasic);
43: return 0;
44: }
46: PetscErrorCode AOBasicGetIndices_Private(AO ao, PetscInt **app, PetscInt **petsc)
47: {
48: AO_Basic *basic = (AO_Basic *)ao->data;
50: if (app) *app = basic->app;
51: if (petsc) *petsc = basic->petsc;
52: return 0;
53: }
55: PetscErrorCode AOPetscToApplication_Basic(AO ao, PetscInt n, PetscInt *ia)
56: {
57: PetscInt i, N = ao->N;
58: AO_Basic *aobasic = (AO_Basic *)ao->data;
60: for (i = 0; i < n; i++) {
61: if (ia[i] >= 0 && ia[i] < N) {
62: ia[i] = aobasic->app[ia[i]];
63: } else {
64: ia[i] = -1;
65: }
66: }
67: return 0;
68: }
70: PetscErrorCode AOApplicationToPetsc_Basic(AO ao, PetscInt n, PetscInt *ia)
71: {
72: PetscInt i, N = ao->N;
73: AO_Basic *aobasic = (AO_Basic *)ao->data;
75: for (i = 0; i < n; i++) {
76: if (ia[i] >= 0 && ia[i] < N) {
77: ia[i] = aobasic->petsc[ia[i]];
78: } else {
79: ia[i] = -1;
80: }
81: }
82: return 0;
83: }
85: PetscErrorCode AOPetscToApplicationPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array)
86: {
87: AO_Basic *aobasic = (AO_Basic *)ao->data;
88: PetscInt *temp;
89: PetscInt i, j;
91: PetscMalloc1(ao->N * block, &temp);
92: for (i = 0; i < ao->N; i++) {
93: for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->petsc[i] * block + j];
94: }
95: PetscArraycpy(array, temp, ao->N * block);
96: PetscFree(temp);
97: return 0;
98: }
100: PetscErrorCode AOApplicationToPetscPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array)
101: {
102: AO_Basic *aobasic = (AO_Basic *)ao->data;
103: PetscInt *temp;
104: PetscInt i, j;
106: PetscMalloc1(ao->N * block, &temp);
107: for (i = 0; i < ao->N; i++) {
108: for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->app[i] * block + j];
109: }
110: PetscArraycpy(array, temp, ao->N * block);
111: PetscFree(temp);
112: return 0;
113: }
115: PetscErrorCode AOPetscToApplicationPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array)
116: {
117: AO_Basic *aobasic = (AO_Basic *)ao->data;
118: PetscReal *temp;
119: PetscInt i, j;
121: PetscMalloc1(ao->N * block, &temp);
122: for (i = 0; i < ao->N; i++) {
123: for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->petsc[i] * block + j];
124: }
125: PetscArraycpy(array, temp, ao->N * block);
126: PetscFree(temp);
127: return 0;
128: }
130: PetscErrorCode AOApplicationToPetscPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array)
131: {
132: AO_Basic *aobasic = (AO_Basic *)ao->data;
133: PetscReal *temp;
134: PetscInt i, j;
136: PetscMalloc1(ao->N * block, &temp);
137: for (i = 0; i < ao->N; i++) {
138: for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->app[i] * block + j];
139: }
140: PetscArraycpy(array, temp, ao->N * block);
141: PetscFree(temp);
142: return 0;
143: }
145: static struct _AOOps AOOps_Basic = {
146: PetscDesignatedInitializer(view, AOView_Basic),
147: PetscDesignatedInitializer(destroy, AODestroy_Basic),
148: PetscDesignatedInitializer(petsctoapplication, AOPetscToApplication_Basic),
149: PetscDesignatedInitializer(applicationtopetsc, AOApplicationToPetsc_Basic),
150: PetscDesignatedInitializer(petsctoapplicationpermuteint, AOPetscToApplicationPermuteInt_Basic),
151: PetscDesignatedInitializer(applicationtopetscpermuteint, AOApplicationToPetscPermuteInt_Basic),
152: PetscDesignatedInitializer(petsctoapplicationpermutereal, AOPetscToApplicationPermuteReal_Basic),
153: PetscDesignatedInitializer(applicationtopetscpermutereal, AOApplicationToPetscPermuteReal_Basic),
154: };
156: PETSC_EXTERN PetscErrorCode AOCreate_Basic(AO ao)
157: {
158: AO_Basic *aobasic;
159: PetscMPIInt size, rank, count, *lens, *disp;
160: PetscInt napp, *allpetsc, *allapp, ip, ia, N, i, *petsc = NULL, start;
161: IS isapp = ao->isapp, ispetsc = ao->ispetsc;
162: MPI_Comm comm;
163: const PetscInt *myapp, *mypetsc = NULL;
165: /* create special struct aobasic */
166: PetscNew(&aobasic);
167: ao->data = (void *)aobasic;
168: PetscMemcpy(ao->ops, &AOOps_Basic, sizeof(struct _AOOps));
169: PetscObjectChangeTypeName((PetscObject)ao, AOBASIC);
171: ISGetLocalSize(isapp, &napp);
172: ISGetIndices(isapp, &myapp);
174: PetscMPIIntCast(napp, &count);
176: /* transmit all lengths to all processors */
177: PetscObjectGetComm((PetscObject)isapp, &comm);
178: MPI_Comm_size(comm, &size);
179: MPI_Comm_rank(comm, &rank);
180: PetscMalloc2(size, &lens, size, &disp);
181: MPI_Allgather(&count, 1, MPI_INT, lens, 1, MPI_INT, comm);
182: N = 0;
183: for (i = 0; i < size; i++) {
184: PetscMPIIntCast(N, disp + i); /* = sum(lens[j]), j< i */
185: N += lens[i];
186: }
187: ao->N = N;
188: ao->n = N;
190: /* If mypetsc is 0 then use "natural" numbering */
191: if (napp) {
192: if (!ispetsc) {
193: start = disp[rank];
194: PetscMalloc1(napp + 1, &petsc);
195: for (i = 0; i < napp; i++) petsc[i] = start + i;
196: } else {
197: ISGetIndices(ispetsc, &mypetsc);
198: petsc = (PetscInt *)mypetsc;
199: }
200: }
202: /* get all indices on all processors */
203: PetscMalloc2(N, &allpetsc, N, &allapp);
204: MPI_Allgatherv(petsc, count, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm);
205: MPI_Allgatherv((void *)myapp, count, MPIU_INT, allapp, lens, disp, MPIU_INT, comm);
206: PetscFree2(lens, disp);
208: if (PetscDefined(USE_DEBUG)) {
209: PetscInt *sorted;
210: PetscMalloc1(N, &sorted);
212: PetscArraycpy(sorted, allpetsc, N);
213: PetscSortInt(N, sorted);
216: PetscArraycpy(sorted, allapp, N);
217: PetscSortInt(N, sorted);
220: PetscFree(sorted);
221: }
223: /* generate a list of application and PETSc node numbers */
224: PetscCalloc2(N, &aobasic->app, N, &aobasic->petsc);
225: for (i = 0; i < N; i++) {
226: ip = allpetsc[i];
227: ia = allapp[i];
228: /* check there are no duplicates */
230: aobasic->app[ip] = ia + 1;
232: aobasic->petsc[ia] = ip + 1;
233: }
234: if (napp && !mypetsc) PetscFree(petsc);
235: PetscFree2(allpetsc, allapp);
236: /* shift indices down by one */
237: for (i = 0; i < N; i++) {
238: aobasic->app[i]--;
239: aobasic->petsc[i]--;
240: }
242: ISRestoreIndices(isapp, &myapp);
243: if (napp) {
244: if (ispetsc) {
245: ISRestoreIndices(ispetsc, &mypetsc);
246: } else {
247: PetscFree(petsc);
248: }
249: }
250: return 0;
251: }
253: /*@C
254: AOCreateBasic - Creates a basic application ordering using two integer arrays.
256: Collective
258: Input Parameters:
259: + comm - MPI communicator that is to share `AO`
260: . napp - size of integer arrays
261: . myapp - integer array that defines an ordering
262: - mypetsc - integer array that defines another ordering (may be NULL to
263: indicate the natural ordering, that is 0,1,2,3,...)
265: Output Parameter:
266: . aoout - the new application ordering
268: Level: beginner
270: Note:
271: The arrays myapp and mypetsc must contain the all the integers 0 to napp-1 with no duplicates; that is there cannot be any "holes"
272: in the indices. Use `AOCreateMapping()` or `AOCreateMappingIS()` if you wish to have "holes" in the indices.
274: .seealso: [](sec_ao), [](sec_scatter), `AO`, `AOCreateBasicIS()`, `AODestroy()`, `AOPetscToApplication()`, `AOApplicationToPetsc()`
275: @*/
276: PetscErrorCode AOCreateBasic(MPI_Comm comm, PetscInt napp, const PetscInt myapp[], const PetscInt mypetsc[], AO *aoout)
277: {
278: IS isapp, ispetsc;
279: const PetscInt *app = myapp, *petsc = mypetsc;
281: ISCreateGeneral(comm, napp, app, PETSC_USE_POINTER, &isapp);
282: if (mypetsc) {
283: ISCreateGeneral(comm, napp, petsc, PETSC_USE_POINTER, &ispetsc);
284: } else {
285: ispetsc = NULL;
286: }
287: AOCreateBasicIS(isapp, ispetsc, aoout);
288: ISDestroy(&isapp);
289: if (mypetsc) ISDestroy(&ispetsc);
290: return 0;
291: }
293: /*@C
294: AOCreateBasicIS - Creates a basic application ordering using two `IS` index sets.
296: Collective
298: Input Parameters:
299: + isapp - index set that defines an ordering
300: - ispetsc - index set that defines another ordering (may be NULL to use the
301: natural ordering)
303: Output Parameter:
304: . aoout - the new application ordering
306: Level: beginner
308: Note:
309: The index sets isapp and ispetsc must contain the all the integers 0 to napp-1 (where napp is the length of the index sets) with no duplicates;
310: that is there cannot be any "holes"
312: .seealso: [](sec_ao), [](sec_scatter), `IS`, `AO`, `AOCreateBasic()`, `AODestroy()`
313: @*/
314: PetscErrorCode AOCreateBasicIS(IS isapp, IS ispetsc, AO *aoout)
315: {
316: MPI_Comm comm;
317: AO ao;
319: PetscObjectGetComm((PetscObject)isapp, &comm);
320: AOCreate(comm, &ao);
321: AOSetIS(ao, isapp, ispetsc);
322: AOSetType(ao, AOBASIC);
323: AOViewFromOptions(ao, NULL, "-ao_view");
324: *aoout = ao;
325: return 0;
326: }