1 /*
2 The most basic AO application ordering routines. These store the
3 entire orderings on each processor to be efficient but can require excessive memory
4 */
5
6 #include <../src/vec/is/ao/aoimpl.h> /*I "petscao.h" I*/
7
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;
12
13 /*
14 All processors have the same data so processor 1 prints it
15 */
AOView_Basic(AO ao,PetscViewer viewer)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 isascii;
22
23 PetscFunctionBegin;
24 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ao), &rank));
25 if (rank == 0) {
26 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
27 if (isascii) {
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 }
36
AODestroy_Basic(AO ao)37 static PetscErrorCode AODestroy_Basic(AO ao)
38 {
39 AO_Basic *aobasic = (AO_Basic *)ao->data;
40
41 PetscFunctionBegin;
42 PetscCall(PetscFree2(aobasic->app, aobasic->petsc));
43 PetscCall(PetscFree(aobasic));
44 PetscFunctionReturn(PETSC_SUCCESS);
45 }
46
AOPetscToApplication_Basic(AO ao,PetscInt n,PetscInt * ia)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;
51
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 }
62
AOApplicationToPetsc_Basic(AO ao,PetscInt n,PetscInt * ia)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;
67
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 }
78
AOPetscToApplicationPermuteInt_Basic(AO ao,PetscInt block,PetscInt * array)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;
84
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 }
94
AOApplicationToPetscPermuteInt_Basic(AO ao,PetscInt block,PetscInt * array)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;
100
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 }
110
AOPetscToApplicationPermuteReal_Basic(AO ao,PetscInt block,PetscReal * array)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;
116
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 }
126
AOApplicationToPetscPermuteReal_Basic(AO ao,PetscInt block,PetscReal * array)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;
132
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 }
142
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 };
153
AOCreate_Basic(AO ao)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;
162
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));
169
170 PetscCall(ISGetLocalSize(isapp, &napp));
171 PetscCall(ISGetIndices(isapp, &myapp));
172
173 PetscCall(PetscMPIIntCast(napp, &count));
174
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;
188
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 }
200
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));
206
207 if (PetscDefined(USE_DEBUG)) {
208 PetscInt *sorted;
209 PetscCall(PetscMalloc1(N, &sorted));
210
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]);
214
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]);
218
219 PetscCall(PetscFree(sorted));
220 }
221
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 }
240
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 }
251
252 /*@
253 AOCreateBasic - Creates a basic application ordering using two integer arrays.
254
255 Collective
256
257 Input Parameters:
258 + comm - MPI communicator that is to share `AO`
259 . napp - size of `myapp` and `mypetsc`
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,...)
263
264 Output Parameter:
265 . aoout - the new application ordering
266
267 Level: beginner
268
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.
272
273 .seealso: [](sec_ao), [](sec_scatter), `AO`, `AOCreateBasicIS()`, `AODestroy()`, `AOPetscToApplication()`, `AOApplicationToPetsc()`
274 @*/
AOCreateBasic(MPI_Comm comm,PetscInt napp,const PetscInt myapp[],const PetscInt mypetsc[],AO * aoout)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;
279
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 }
292
293 /*@
294 AOCreateBasicIS - Creates a basic application ordering using two `IS` index sets.
295
296 Collective
297
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 natural ordering)
301
302 Output Parameter:
303 . aoout - the new application ordering
304
305 Level: beginner
306
307 Note:
308 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;
309 that is there cannot be any "holes"
310
311 .seealso: [](sec_ao), [](sec_scatter), `IS`, `AO`, `AOCreateBasic()`, `AODestroy()`
312 @*/
AOCreateBasicIS(IS isapp,IS ispetsc,AO * aoout)313 PetscErrorCode AOCreateBasicIS(IS isapp, IS ispetsc, AO *aoout)
314 {
315 MPI_Comm comm;
316 AO ao;
317
318 PetscFunctionBegin;
319 PetscCall(PetscObjectGetComm((PetscObject)isapp, &comm));
320 PetscCall(AOCreate(comm, &ao));
321 PetscCall(AOSetIS(ao, isapp, ispetsc));
322 PetscCall(AOSetType(ao, AOBASIC));
323 PetscCall(AOViewFromOptions(ao, NULL, "-ao_view"));
324 *aoout = ao;
325 PetscFunctionReturn(PETSC_SUCCESS);
326 }
327