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: PetscFunctionBegin;
25: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ao), &rank));
26: if (rank == 0) {
27: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
28: if (iascii) {
29: PetscCall(PetscViewerASCIIPrintf(viewer, "Number of elements in ordering %" PetscInt_FMT "\n", ao->N));
30: PetscCall(PetscViewerASCIIPrintf(viewer, "PETSc->App App->PETSc\n"));
31: for (i = 0; i < ao->N; i++) PetscCall(PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT " %3" PetscInt_FMT " %3" PetscInt_FMT " %3" PetscInt_FMT "\n", i, aobasic->app[i], i, aobasic->petsc[i]));
32: }
33: }
34: PetscCall(PetscViewerFlush(viewer));
35: PetscFunctionReturn(PETSC_SUCCESS);
36: }
38: PetscErrorCode AODestroy_Basic(AO ao)
39: {
40: AO_Basic *aobasic = (AO_Basic *)ao->data;
42: PetscFunctionBegin;
43: PetscCall(PetscFree2(aobasic->app, aobasic->petsc));
44: PetscCall(PetscFree(aobasic));
45: PetscFunctionReturn(PETSC_SUCCESS);
46: }
48: PetscErrorCode AOBasicGetIndices_Private(AO ao, PetscInt **app, PetscInt **petsc)
49: {
50: AO_Basic *basic = (AO_Basic *)ao->data;
52: PetscFunctionBegin;
53: if (app) *app = basic->app;
54: if (petsc) *petsc = basic->petsc;
55: PetscFunctionReturn(PETSC_SUCCESS);
56: }
58: PetscErrorCode AOPetscToApplication_Basic(AO ao, PetscInt n, PetscInt *ia)
59: {
60: PetscInt i, N = ao->N;
61: AO_Basic *aobasic = (AO_Basic *)ao->data;
63: PetscFunctionBegin;
64: for (i = 0; i < n; i++) {
65: if (ia[i] >= 0 && ia[i] < N) {
66: ia[i] = aobasic->app[ia[i]];
67: } else {
68: ia[i] = -1;
69: }
70: }
71: PetscFunctionReturn(PETSC_SUCCESS);
72: }
74: PetscErrorCode AOApplicationToPetsc_Basic(AO ao, PetscInt n, PetscInt *ia)
75: {
76: PetscInt i, N = ao->N;
77: AO_Basic *aobasic = (AO_Basic *)ao->data;
79: PetscFunctionBegin;
80: for (i = 0; i < n; i++) {
81: if (ia[i] >= 0 && ia[i] < N) {
82: ia[i] = aobasic->petsc[ia[i]];
83: } else {
84: ia[i] = -1;
85: }
86: }
87: PetscFunctionReturn(PETSC_SUCCESS);
88: }
90: PetscErrorCode AOPetscToApplicationPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array)
91: {
92: AO_Basic *aobasic = (AO_Basic *)ao->data;
93: PetscInt *temp;
94: PetscInt i, j;
96: PetscFunctionBegin;
97: PetscCall(PetscMalloc1(ao->N * block, &temp));
98: for (i = 0; i < ao->N; i++) {
99: for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->petsc[i] * block + j];
100: }
101: PetscCall(PetscArraycpy(array, temp, ao->N * block));
102: PetscCall(PetscFree(temp));
103: PetscFunctionReturn(PETSC_SUCCESS);
104: }
106: PetscErrorCode AOApplicationToPetscPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array)
107: {
108: AO_Basic *aobasic = (AO_Basic *)ao->data;
109: PetscInt *temp;
110: PetscInt i, j;
112: PetscFunctionBegin;
113: PetscCall(PetscMalloc1(ao->N * block, &temp));
114: for (i = 0; i < ao->N; i++) {
115: for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->app[i] * block + j];
116: }
117: PetscCall(PetscArraycpy(array, temp, ao->N * block));
118: PetscCall(PetscFree(temp));
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: PetscErrorCode AOPetscToApplicationPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array)
123: {
124: AO_Basic *aobasic = (AO_Basic *)ao->data;
125: PetscReal *temp;
126: PetscInt i, j;
128: PetscFunctionBegin;
129: PetscCall(PetscMalloc1(ao->N * block, &temp));
130: for (i = 0; i < ao->N; i++) {
131: for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->petsc[i] * block + j];
132: }
133: PetscCall(PetscArraycpy(array, temp, ao->N * block));
134: PetscCall(PetscFree(temp));
135: PetscFunctionReturn(PETSC_SUCCESS);
136: }
138: PetscErrorCode AOApplicationToPetscPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array)
139: {
140: AO_Basic *aobasic = (AO_Basic *)ao->data;
141: PetscReal *temp;
142: PetscInt i, j;
144: PetscFunctionBegin;
145: PetscCall(PetscMalloc1(ao->N * block, &temp));
146: for (i = 0; i < ao->N; i++) {
147: for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->app[i] * block + j];
148: }
149: PetscCall(PetscArraycpy(array, temp, ao->N * block));
150: PetscCall(PetscFree(temp));
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: static struct _AOOps AOOps_Basic = {
155: PetscDesignatedInitializer(view, AOView_Basic),
156: PetscDesignatedInitializer(destroy, AODestroy_Basic),
157: PetscDesignatedInitializer(petsctoapplication, AOPetscToApplication_Basic),
158: PetscDesignatedInitializer(applicationtopetsc, AOApplicationToPetsc_Basic),
159: PetscDesignatedInitializer(petsctoapplicationpermuteint, AOPetscToApplicationPermuteInt_Basic),
160: PetscDesignatedInitializer(applicationtopetscpermuteint, AOApplicationToPetscPermuteInt_Basic),
161: PetscDesignatedInitializer(petsctoapplicationpermutereal, AOPetscToApplicationPermuteReal_Basic),
162: PetscDesignatedInitializer(applicationtopetscpermutereal, AOApplicationToPetscPermuteReal_Basic),
163: };
165: PETSC_EXTERN PetscErrorCode AOCreate_Basic(AO ao)
166: {
167: AO_Basic *aobasic;
168: PetscMPIInt size, rank, count, *lens, *disp;
169: PetscInt napp, *allpetsc, *allapp, ip, ia, N, i, *petsc = NULL, start;
170: IS isapp = ao->isapp, ispetsc = ao->ispetsc;
171: MPI_Comm comm;
172: const PetscInt *myapp, *mypetsc = NULL;
174: PetscFunctionBegin;
175: /* create special struct aobasic */
176: PetscCall(PetscNew(&aobasic));
177: ao->data = (void *)aobasic;
178: PetscCall(PetscMemcpy(ao->ops, &AOOps_Basic, sizeof(struct _AOOps)));
179: PetscCall(PetscObjectChangeTypeName((PetscObject)ao, AOBASIC));
181: PetscCall(ISGetLocalSize(isapp, &napp));
182: PetscCall(ISGetIndices(isapp, &myapp));
184: PetscCall(PetscMPIIntCast(napp, &count));
186: /* transmit all lengths to all processors */
187: PetscCall(PetscObjectGetComm((PetscObject)isapp, &comm));
188: PetscCallMPI(MPI_Comm_size(comm, &size));
189: PetscCallMPI(MPI_Comm_rank(comm, &rank));
190: PetscCall(PetscMalloc2(size, &lens, size, &disp));
191: PetscCallMPI(MPI_Allgather(&count, 1, MPI_INT, lens, 1, MPI_INT, comm));
192: N = 0;
193: for (i = 0; i < size; i++) {
194: PetscCall(PetscMPIIntCast(N, disp + i)); /* = sum(lens[j]), j< i */
195: N += lens[i];
196: }
197: ao->N = N;
198: ao->n = N;
200: /* If mypetsc is 0 then use "natural" numbering */
201: if (napp) {
202: if (!ispetsc) {
203: start = disp[rank];
204: PetscCall(PetscMalloc1(napp + 1, &petsc));
205: for (i = 0; i < napp; i++) petsc[i] = start + i;
206: } else {
207: PetscCall(ISGetIndices(ispetsc, &mypetsc));
208: petsc = (PetscInt *)mypetsc;
209: }
210: }
212: /* get all indices on all processors */
213: PetscCall(PetscMalloc2(N, &allpetsc, N, &allapp));
214: PetscCallMPI(MPI_Allgatherv(petsc, count, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm));
215: PetscCallMPI(MPI_Allgatherv((void *)myapp, count, MPIU_INT, allapp, lens, disp, MPIU_INT, comm));
216: PetscCall(PetscFree2(lens, disp));
218: if (PetscDefined(USE_DEBUG)) {
219: PetscInt *sorted;
220: PetscCall(PetscMalloc1(N, &sorted));
222: PetscCall(PetscArraycpy(sorted, allpetsc, N));
223: PetscCall(PetscSortInt(N, sorted));
224: for (i = 0; i < N; i++) PetscCheck(sorted[i] == i, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PETSc ordering requires a permutation of numbers 0 to N-1\n it is missing %" PetscInt_FMT " has %" PetscInt_FMT, i, sorted[i]);
226: PetscCall(PetscArraycpy(sorted, allapp, N));
227: PetscCall(PetscSortInt(N, sorted));
228: for (i = 0; i < N; i++) PetscCheck(sorted[i] == i, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Application ordering requires a permutation of numbers 0 to N-1\n it is missing %" PetscInt_FMT " has %" PetscInt_FMT, i, sorted[i]);
230: PetscCall(PetscFree(sorted));
231: }
233: /* generate a list of application and PETSc node numbers */
234: PetscCall(PetscCalloc2(N, &aobasic->app, N, &aobasic->petsc));
235: for (i = 0; i < N; i++) {
236: ip = allpetsc[i];
237: ia = allapp[i];
238: /* check there are no duplicates */
239: PetscCheck(!aobasic->app[ip], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Duplicate in PETSc ordering at position %" PetscInt_FMT ". Already mapped to %" PetscInt_FMT ", not %" PetscInt_FMT ".", i, aobasic->app[ip] - 1, ia);
240: aobasic->app[ip] = ia + 1;
241: PetscCheck(!aobasic->petsc[ia], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Duplicate in Application ordering at position %" PetscInt_FMT ". Already mapped to %" PetscInt_FMT ", not %" PetscInt_FMT ".", i, aobasic->petsc[ia] - 1, ip);
242: aobasic->petsc[ia] = ip + 1;
243: }
244: if (napp && !mypetsc) PetscCall(PetscFree(petsc));
245: PetscCall(PetscFree2(allpetsc, allapp));
246: /* shift indices down by one */
247: for (i = 0; i < N; i++) {
248: aobasic->app[i]--;
249: aobasic->petsc[i]--;
250: }
252: PetscCall(ISRestoreIndices(isapp, &myapp));
253: if (napp) {
254: if (ispetsc) {
255: PetscCall(ISRestoreIndices(ispetsc, &mypetsc));
256: } else {
257: PetscCall(PetscFree(petsc));
258: }
259: }
260: PetscFunctionReturn(PETSC_SUCCESS);
261: }
263: /*@C
264: AOCreateBasic - Creates a basic application ordering using two integer arrays.
266: Collective
268: Input Parameters:
269: + comm - MPI communicator that is to share `AO`
270: . napp - size of integer arrays
271: . myapp - integer array that defines an ordering
272: - mypetsc - integer array that defines another ordering (may be NULL to
273: indicate the natural ordering, that is 0,1,2,3,...)
275: Output Parameter:
276: . aoout - the new application ordering
278: Level: beginner
280: Note:
281: 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"
282: in the indices. Use `AOCreateMapping()` or `AOCreateMappingIS()` if you wish to have "holes" in the indices.
284: .seealso: [](sec_ao), [](sec_scatter), `AO`, `AOCreateBasicIS()`, `AODestroy()`, `AOPetscToApplication()`, `AOApplicationToPetsc()`
285: @*/
286: PetscErrorCode AOCreateBasic(MPI_Comm comm, PetscInt napp, const PetscInt myapp[], const PetscInt mypetsc[], AO *aoout)
287: {
288: IS isapp, ispetsc;
289: const PetscInt *app = myapp, *petsc = mypetsc;
291: PetscFunctionBegin;
292: PetscCall(ISCreateGeneral(comm, napp, app, PETSC_USE_POINTER, &isapp));
293: if (mypetsc) {
294: PetscCall(ISCreateGeneral(comm, napp, petsc, PETSC_USE_POINTER, &ispetsc));
295: } else {
296: ispetsc = NULL;
297: }
298: PetscCall(AOCreateBasicIS(isapp, ispetsc, aoout));
299: PetscCall(ISDestroy(&isapp));
300: if (mypetsc) PetscCall(ISDestroy(&ispetsc));
301: PetscFunctionReturn(PETSC_SUCCESS);
302: }
304: /*@C
305: AOCreateBasicIS - Creates a basic application ordering using two `IS` index sets.
307: Collective
309: Input Parameters:
310: + isapp - index set that defines an ordering
311: - ispetsc - index set that defines another ordering (may be NULL to use the
312: natural ordering)
314: Output Parameter:
315: . aoout - the new application ordering
317: Level: beginner
319: Note:
320: 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;
321: that is there cannot be any "holes"
323: .seealso: [](sec_ao), [](sec_scatter), `IS`, `AO`, `AOCreateBasic()`, `AODestroy()`
324: @*/
325: PetscErrorCode AOCreateBasicIS(IS isapp, IS ispetsc, AO *aoout)
326: {
327: MPI_Comm comm;
328: AO ao;
330: PetscFunctionBegin;
331: PetscCall(PetscObjectGetComm((PetscObject)isapp, &comm));
332: PetscCall(AOCreate(comm, &ao));
333: PetscCall(AOSetIS(ao, isapp, ispetsc));
334: PetscCall(AOSetType(ao, AOBASIC));
335: PetscCall(AOViewFromOptions(ao, NULL, "-ao_view"));
336: *aoout = ao;
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }