Actual source code: aobasic.c
1: /*
2: The most basic AO application ordering routines. These store the
3: entire orderings on each processor.
4: */
6: #include <../src/vec/is/ao/aoimpl.h>
8: typedef struct {
9: PetscInt *app; /* app[i] is the partner for the ith PETSc slot */
10: PetscInt *petsc; /* petsc[j] is the partner for the jth app slot */
11: } AO_Basic;
13: /*
14: All processors have the same data so processor 1 prints it
15: */
16: static PetscErrorCode AOView_Basic(AO ao, PetscViewer viewer)
17: {
18: PetscMPIInt rank;
19: PetscInt i;
20: AO_Basic *aobasic = (AO_Basic *)ao->data;
21: PetscBool iascii;
23: PetscFunctionBegin;
24: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ao), &rank));
25: if (rank == 0) {
26: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
27: if (iascii) {
28: PetscCall(PetscViewerASCIIPrintf(viewer, "Number of elements in ordering %" PetscInt_FMT "\n", ao->N));
29: PetscCall(PetscViewerASCIIPrintf(viewer, "PETSc->App App->PETSc\n"));
30: 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]));
31: }
32: }
33: PetscCall(PetscViewerFlush(viewer));
34: PetscFunctionReturn(PETSC_SUCCESS);
35: }
37: static PetscErrorCode AODestroy_Basic(AO ao)
38: {
39: AO_Basic *aobasic = (AO_Basic *)ao->data;
41: PetscFunctionBegin;
42: PetscCall(PetscFree2(aobasic->app, aobasic->petsc));
43: PetscCall(PetscFree(aobasic));
44: PetscFunctionReturn(PETSC_SUCCESS);
45: }
47: static PetscErrorCode AOPetscToApplication_Basic(AO ao, PetscInt n, PetscInt *ia)
48: {
49: PetscInt i, N = ao->N;
50: AO_Basic *aobasic = (AO_Basic *)ao->data;
52: PetscFunctionBegin;
53: for (i = 0; i < n; i++) {
54: if (ia[i] >= 0 && ia[i] < N) {
55: ia[i] = aobasic->app[ia[i]];
56: } else {
57: ia[i] = -1;
58: }
59: }
60: PetscFunctionReturn(PETSC_SUCCESS);
61: }
63: static PetscErrorCode AOApplicationToPetsc_Basic(AO ao, PetscInt n, PetscInt *ia)
64: {
65: PetscInt i, N = ao->N;
66: AO_Basic *aobasic = (AO_Basic *)ao->data;
68: PetscFunctionBegin;
69: for (i = 0; i < n; i++) {
70: if (ia[i] >= 0 && ia[i] < N) {
71: ia[i] = aobasic->petsc[ia[i]];
72: } else {
73: ia[i] = -1;
74: }
75: }
76: PetscFunctionReturn(PETSC_SUCCESS);
77: }
79: static PetscErrorCode AOPetscToApplicationPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array)
80: {
81: AO_Basic *aobasic = (AO_Basic *)ao->data;
82: PetscInt *temp;
83: PetscInt i, j;
85: PetscFunctionBegin;
86: PetscCall(PetscMalloc1(ao->N * block, &temp));
87: for (i = 0; i < ao->N; i++) {
88: for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->petsc[i] * block + j];
89: }
90: PetscCall(PetscArraycpy(array, temp, ao->N * block));
91: PetscCall(PetscFree(temp));
92: PetscFunctionReturn(PETSC_SUCCESS);
93: }
95: static PetscErrorCode AOApplicationToPetscPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array)
96: {
97: AO_Basic *aobasic = (AO_Basic *)ao->data;
98: PetscInt *temp;
99: PetscInt i, j;
101: PetscFunctionBegin;
102: PetscCall(PetscMalloc1(ao->N * block, &temp));
103: for (i = 0; i < ao->N; i++) {
104: for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->app[i] * block + j];
105: }
106: PetscCall(PetscArraycpy(array, temp, ao->N * block));
107: PetscCall(PetscFree(temp));
108: PetscFunctionReturn(PETSC_SUCCESS);
109: }
111: static PetscErrorCode AOPetscToApplicationPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array)
112: {
113: AO_Basic *aobasic = (AO_Basic *)ao->data;
114: PetscReal *temp;
115: PetscInt i, j;
117: PetscFunctionBegin;
118: PetscCall(PetscMalloc1(ao->N * block, &temp));
119: for (i = 0; i < ao->N; i++) {
120: for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->petsc[i] * block + j];
121: }
122: PetscCall(PetscArraycpy(array, temp, ao->N * block));
123: PetscCall(PetscFree(temp));
124: PetscFunctionReturn(PETSC_SUCCESS);
125: }
127: static PetscErrorCode AOApplicationToPetscPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array)
128: {
129: AO_Basic *aobasic = (AO_Basic *)ao->data;
130: PetscReal *temp;
131: PetscInt i, j;
133: PetscFunctionBegin;
134: PetscCall(PetscMalloc1(ao->N * block, &temp));
135: for (i = 0; i < ao->N; i++) {
136: for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->app[i] * block + j];
137: }
138: PetscCall(PetscArraycpy(array, temp, ao->N * block));
139: PetscCall(PetscFree(temp));
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
143: static const struct _AOOps AOOps_Basic = {
144: PetscDesignatedInitializer(view, AOView_Basic),
145: PetscDesignatedInitializer(destroy, AODestroy_Basic),
146: PetscDesignatedInitializer(petsctoapplication, AOPetscToApplication_Basic),
147: PetscDesignatedInitializer(applicationtopetsc, AOApplicationToPetsc_Basic),
148: PetscDesignatedInitializer(petsctoapplicationpermuteint, AOPetscToApplicationPermuteInt_Basic),
149: PetscDesignatedInitializer(applicationtopetscpermuteint, AOApplicationToPetscPermuteInt_Basic),
150: PetscDesignatedInitializer(petsctoapplicationpermutereal, AOPetscToApplicationPermuteReal_Basic),
151: PetscDesignatedInitializer(applicationtopetscpermutereal, AOApplicationToPetscPermuteReal_Basic),
152: };
154: PETSC_INTERN PetscErrorCode AOCreate_Basic(AO ao)
155: {
156: AO_Basic *aobasic;
157: PetscMPIInt size, rank, count, *lens, *disp;
158: PetscInt napp, *allpetsc, *allapp, ip, ia, N, i, *petsc = NULL, start;
159: IS isapp = ao->isapp, ispetsc = ao->ispetsc;
160: MPI_Comm comm;
161: const PetscInt *myapp, *mypetsc = NULL;
163: PetscFunctionBegin;
164: /* create special struct aobasic */
165: PetscCall(PetscNew(&aobasic));
166: ao->data = (void *)aobasic;
167: ao->ops[0] = AOOps_Basic;
168: PetscCall(PetscObjectChangeTypeName((PetscObject)ao, AOBASIC));
170: PetscCall(ISGetLocalSize(isapp, &napp));
171: PetscCall(ISGetIndices(isapp, &myapp));
173: PetscCall(PetscMPIIntCast(napp, &count));
175: /* transmit all lengths to all processors */
176: PetscCall(PetscObjectGetComm((PetscObject)isapp, &comm));
177: PetscCallMPI(MPI_Comm_size(comm, &size));
178: PetscCallMPI(MPI_Comm_rank(comm, &rank));
179: PetscCall(PetscMalloc2(size, &lens, size, &disp));
180: PetscCallMPI(MPI_Allgather(&count, 1, MPI_INT, lens, 1, MPI_INT, comm));
181: N = 0;
182: for (i = 0; i < size; i++) {
183: PetscCall(PetscMPIIntCast(N, disp + i)); /* = sum(lens[j]), j< i */
184: N += lens[i];
185: }
186: ao->N = N;
187: ao->n = N;
189: /* If mypetsc is 0 then use "natural" numbering */
190: if (napp) {
191: if (!ispetsc) {
192: start = disp[rank];
193: PetscCall(PetscMalloc1(napp + 1, &petsc));
194: for (i = 0; i < napp; i++) petsc[i] = start + i;
195: } else {
196: PetscCall(ISGetIndices(ispetsc, &mypetsc));
197: petsc = (PetscInt *)mypetsc;
198: }
199: }
201: /* get all indices on all processors */
202: PetscCall(PetscMalloc2(N, &allpetsc, N, &allapp));
203: PetscCallMPI(MPI_Allgatherv(petsc, count, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm));
204: PetscCallMPI(MPI_Allgatherv((void *)myapp, count, MPIU_INT, allapp, lens, disp, MPIU_INT, comm));
205: PetscCall(PetscFree2(lens, disp));
207: if (PetscDefined(USE_DEBUG)) {
208: PetscInt *sorted;
209: PetscCall(PetscMalloc1(N, &sorted));
211: PetscCall(PetscArraycpy(sorted, allpetsc, N));
212: PetscCall(PetscSortInt(N, sorted));
213: 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, it is missing %" PetscInt_FMT " has %" PetscInt_FMT, i, sorted[i]);
215: PetscCall(PetscArraycpy(sorted, allapp, N));
216: PetscCall(PetscSortInt(N, sorted));
217: 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, it is missing %" PetscInt_FMT " has %" PetscInt_FMT, i, sorted[i]);
219: PetscCall(PetscFree(sorted));
220: }
222: /* generate a list of application and PETSc node numbers */
223: PetscCall(PetscCalloc2(N, &aobasic->app, N, &aobasic->petsc));
224: for (i = 0; i < N; i++) {
225: ip = allpetsc[i];
226: ia = allapp[i];
227: /* check there are no duplicates */
228: 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);
229: aobasic->app[ip] = ia + 1;
230: 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);
231: aobasic->petsc[ia] = ip + 1;
232: }
233: if (napp && !mypetsc) PetscCall(PetscFree(petsc));
234: PetscCall(PetscFree2(allpetsc, allapp));
235: /* shift indices down by one */
236: for (i = 0; i < N; i++) {
237: aobasic->app[i]--;
238: aobasic->petsc[i]--;
239: }
241: PetscCall(ISRestoreIndices(isapp, &myapp));
242: if (napp) {
243: if (ispetsc) {
244: PetscCall(ISRestoreIndices(ispetsc, &mypetsc));
245: } else {
246: PetscCall(PetscFree(petsc));
247: }
248: }
249: PetscFunctionReturn(PETSC_SUCCESS);
250: }
252: /*@C
253: AOCreateBasic - Creates a basic application ordering using two integer arrays.
255: Collective
257: Input Parameters:
258: + comm - MPI communicator that is to share `AO`
259: . napp - size of integer arrays
260: . myapp - integer array that defines an ordering
261: - mypetsc - integer array that defines another ordering (may be `NULL` to
262: indicate the natural ordering, that is 0,1,2,3,...)
264: Output Parameter:
265: . aoout - the new application ordering
267: Level: beginner
269: Note:
270: 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"
271: in the indices. Use `AOCreateMapping()` or `AOCreateMappingIS()` if you wish to have "holes" in the indices.
273: .seealso: [](sec_ao), [](sec_scatter), `AO`, `AOCreateBasicIS()`, `AODestroy()`, `AOPetscToApplication()`, `AOApplicationToPetsc()`
274: @*/
275: PetscErrorCode AOCreateBasic(MPI_Comm comm, PetscInt napp, const PetscInt myapp[], const PetscInt mypetsc[], AO *aoout)
276: {
277: IS isapp, ispetsc;
278: const PetscInt *app = myapp, *petsc = mypetsc;
280: PetscFunctionBegin;
281: PetscCall(ISCreateGeneral(comm, napp, app, PETSC_USE_POINTER, &isapp));
282: if (mypetsc) {
283: PetscCall(ISCreateGeneral(comm, napp, petsc, PETSC_USE_POINTER, &ispetsc));
284: } else {
285: ispetsc = NULL;
286: }
287: PetscCall(AOCreateBasicIS(isapp, ispetsc, aoout));
288: PetscCall(ISDestroy(&isapp));
289: if (mypetsc) PetscCall(ISDestroy(&ispetsc));
290: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
320: PetscCall(PetscObjectGetComm((PetscObject)isapp, &comm));
321: PetscCall(AOCreate(comm, &ao));
322: PetscCall(AOSetIS(ao, isapp, ispetsc));
323: PetscCall(AOSetType(ao, AOBASIC));
324: PetscCall(AOViewFromOptions(ao, NULL, "-ao_view"));
325: *aoout = ao;
326: PetscFunctionReturn(PETSC_SUCCESS);
327: }