Actual source code: aomapping.c


  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>

 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;

 18: PetscErrorCode AODestroy_Mapping(AO ao)
 19: {
 20:   AO_Mapping *aomap = (AO_Mapping *)ao->data;

 22:   PetscFunctionBegin;
 23:   PetscCall(PetscFree4(aomap->app, aomap->appPerm, aomap->petsc, aomap->petscPerm));
 24:   PetscCall(PetscFree(aomap));
 25:   PetscFunctionReturn(PETSC_SUCCESS);
 26: }

 28: PetscErrorCode AOView_Mapping(AO ao, PetscViewer viewer)
 29: {
 30:   AO_Mapping *aomap = (AO_Mapping *)ao->data;
 31:   PetscMPIInt rank;
 32:   PetscInt    i;
 33:   PetscBool   iascii;

 35:   PetscFunctionBegin;
 36:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ao), &rank));
 37:   if (rank) PetscFunctionReturn(PETSC_SUCCESS);
 38:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
 39:   if (iascii) {
 40:     PetscCall(PetscViewerASCIIPrintf(viewer, "Number of elements in ordering %" PetscInt_FMT "\n", aomap->N));
 41:     PetscCall(PetscViewerASCIIPrintf(viewer, "   App.   PETSc\n"));
 42:     for (i = 0; i < aomap->N; i++) PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "   %" PetscInt_FMT "    %" PetscInt_FMT "\n", i, aomap->app[i], aomap->petsc[aomap->appPerm[i]]));
 43:   }
 44:   PetscFunctionReturn(PETSC_SUCCESS);
 45: }

 47: PetscErrorCode AOPetscToApplication_Mapping(AO ao, PetscInt n, PetscInt *ia)
 48: {
 49:   AO_Mapping *aomap = (AO_Mapping *)ao->data;
 50:   PetscInt   *app   = aomap->app;
 51:   PetscInt   *petsc = aomap->petsc;
 52:   PetscInt   *perm  = aomap->petscPerm;
 53:   PetscInt    N     = aomap->N;
 54:   PetscInt    low, high, mid = 0;
 55:   PetscInt    idex;
 56:   PetscInt    i;

 58:   /* It would be possible to use a single bisection search, which
 59:      recursively divided the indices to be converted, and searched
 60:      partitions which contained an index. This would result in
 61:      better running times if indices are clustered.
 62:   */
 63:   PetscFunctionBegin;
 64:   for (i = 0; i < n; i++) {
 65:     idex = ia[i];
 66:     if (idex < 0) continue;
 67:     /* Use bisection since the array is sorted */
 68:     low  = 0;
 69:     high = N - 1;
 70:     while (low <= high) {
 71:       mid = (low + high) / 2;
 72:       if (idex == petsc[mid]) break;
 73:       else if (idex < petsc[mid]) high = mid - 1;
 74:       else low = mid + 1;
 75:     }
 76:     if (low > high) ia[i] = -1;
 77:     else ia[i] = app[perm[mid]];
 78:   }
 79:   PetscFunctionReturn(PETSC_SUCCESS);
 80: }

 82: PetscErrorCode AOApplicationToPetsc_Mapping(AO ao, PetscInt n, PetscInt *ia)
 83: {
 84:   AO_Mapping *aomap = (AO_Mapping *)ao->data;
 85:   PetscInt   *app   = aomap->app;
 86:   PetscInt   *petsc = aomap->petsc;
 87:   PetscInt   *perm  = aomap->appPerm;
 88:   PetscInt    N     = aomap->N;
 89:   PetscInt    low, high, mid = 0;
 90:   PetscInt    idex;
 91:   PetscInt    i;

 93:   /* It would be possible to use a single bisection search, which
 94:      recursively divided the indices to be converted, and searched
 95:      partitions which contained an index. This would result in
 96:      better running times if indices are clustered.
 97:   */
 98:   PetscFunctionBegin;
 99:   for (i = 0; i < n; i++) {
100:     idex = ia[i];
101:     if (idex < 0) continue;
102:     /* Use bisection since the array is sorted */
103:     low  = 0;
104:     high = N - 1;
105:     while (low <= high) {
106:       mid = (low + high) / 2;
107:       if (idex == app[mid]) break;
108:       else if (idex < app[mid]) high = mid - 1;
109:       else low = mid + 1;
110:     }
111:     if (low > high) ia[i] = -1;
112:     else ia[i] = petsc[perm[mid]];
113:   }
114:   PetscFunctionReturn(PETSC_SUCCESS);
115: }

117: static struct _AOOps AOps = {
118:   PetscDesignatedInitializer(view, AOView_Mapping),
119:   PetscDesignatedInitializer(destroy, AODestroy_Mapping),
120:   PetscDesignatedInitializer(petsctoapplication, AOPetscToApplication_Mapping),
121:   PetscDesignatedInitializer(applicationtopetsc, AOApplicationToPetsc_Mapping),
122: };

124: /*@C
125:   AOMappingHasApplicationIndex - Checks if an `AO` has a requested application index.

127:   Not Collective

129:   Input Parameters:
130: + ao       - The `AO`
131: - index    - The application index

133:   Output Parameter:
134: . hasIndex - Flag is `PETSC_TRUE` if the index exists

136:   Level: intermediate

138:   Developer Note:
139:   The name of the function is wrong, it should be `AOHasApplicationIndex`

141: .seealso: [](sec_ao), `AOMappingHasPetscIndex()`, `AOCreateMapping()`, `AO`
142: @*/
143: PetscErrorCode AOMappingHasApplicationIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
144: {
145:   AO_Mapping *aomap;
146:   PetscInt   *app;
147:   PetscInt    low, high, mid;

149:   PetscFunctionBegin;
152:   aomap = (AO_Mapping *)ao->data;
153:   app   = aomap->app;
154:   /* Use bisection since the array is sorted */
155:   low  = 0;
156:   high = aomap->N - 1;
157:   while (low <= high) {
158:     mid = (low + high) / 2;
159:     if (idex == app[mid]) break;
160:     else if (idex < app[mid]) high = mid - 1;
161:     else low = mid + 1;
162:   }
163:   if (low > high) *hasIndex = PETSC_FALSE;
164:   else *hasIndex = PETSC_TRUE;
165:   PetscFunctionReturn(PETSC_SUCCESS);
166: }

168: /*@
169:   AOMappingHasPetscIndex - checks if an `AO` has a requested petsc index.

171:   Not Collective

173:   Input Parameters:
174: + ao       - The `AO`
175: - index    - The petsc index

177:   Output Parameter:
178: . hasIndex - Flag is `PETSC_TRUE` if the index exists

180:   Level: intermediate

182:   Developer Note:
183:   The name of the function is wrong, it should be `AOHasPetscIndex`

185: .seealso: [](sec_ao), `AOMappingHasApplicationIndex()`, `AOCreateMapping()`
186: @*/
187: PetscErrorCode AOMappingHasPetscIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
188: {
189:   AO_Mapping *aomap;
190:   PetscInt   *petsc;
191:   PetscInt    low, high, mid;

193:   PetscFunctionBegin;
196:   aomap = (AO_Mapping *)ao->data;
197:   petsc = aomap->petsc;
198:   /* Use bisection since the array is sorted */
199:   low  = 0;
200:   high = aomap->N - 1;
201:   while (low <= high) {
202:     mid = (low + high) / 2;
203:     if (idex == petsc[mid]) break;
204:     else if (idex < petsc[mid]) high = mid - 1;
205:     else low = mid + 1;
206:   }
207:   if (low > high) *hasIndex = PETSC_FALSE;
208:   else *hasIndex = PETSC_TRUE;
209:   PetscFunctionReturn(PETSC_SUCCESS);
210: }

212: /*@C
213:   AOCreateMapping - Creates an application mapping using two integer arrays.

215:   Input Parameters:
216: + comm    - MPI communicator that is to share the `AO`
217: . napp    - size of integer arrays
218: . myapp   - integer array that defines an ordering
219: - mypetsc - integer array that defines another ordering (may be NULL to indicate the identity ordering)

221:   Output Parameter:
222: . aoout   - the new application mapping

224:   Options Database Key:
225: . -ao_view - call `AOView()` at the conclusion of `AOCreateMapping()`

227:   Level: beginner

229:     Note:
230:     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.
231:     Use `AOCreateBasic()` or `AOCreateBasicIS()` if they do not have holes for better performance.

233: .seealso: [](sec_ao), `AOCreateBasic()`, `AOCreateBasic()`, `AOCreateMappingIS()`, `AODestroy()`
234: @*/
235: PetscErrorCode AOCreateMapping(MPI_Comm comm, PetscInt napp, const PetscInt myapp[], const PetscInt mypetsc[], AO *aoout)
236: {
237:   AO          ao;
238:   AO_Mapping *aomap;
239:   PetscInt   *allpetsc, *allapp;
240:   PetscInt   *petscPerm, *appPerm;
241:   PetscInt   *petsc;
242:   PetscMPIInt size, rank, *lens, *disp, nnapp;
243:   PetscInt    N, start;
244:   PetscInt    i;

246:   PetscFunctionBegin;
248:   *aoout = NULL;
249:   PetscCall(AOInitializePackage());

251:   PetscCall(PetscHeaderCreate(ao, AO_CLASSID, "AO", "Application Ordering", "AO", comm, AODestroy, AOView));
252:   PetscCall(PetscNew(&aomap));
253:   PetscCall(PetscMemcpy(ao->ops, &AOps, sizeof(AOps)));
254:   ao->data = (void *)aomap;

256:   /* transmit all lengths to all processors */
257:   PetscCallMPI(MPI_Comm_size(comm, &size));
258:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
259:   PetscCall(PetscMalloc2(size, &lens, size, &disp));
260:   nnapp = napp;
261:   PetscCallMPI(MPI_Allgather(&nnapp, 1, MPI_INT, lens, 1, MPI_INT, comm));
262:   N = 0;
263:   for (i = 0; i < size; i++) {
264:     disp[i] = N;
265:     N += lens[i];
266:   }
267:   aomap->N = N;
268:   ao->N    = N;
269:   ao->n    = N;

271:   /* If mypetsc is 0 then use "natural" numbering */
272:   if (!mypetsc) {
273:     start = disp[rank];
274:     PetscCall(PetscMalloc1(napp + 1, &petsc));
275:     for (i = 0; i < napp; i++) petsc[i] = start + i;
276:   } else {
277:     petsc = (PetscInt *)mypetsc;
278:   }

280:   /* get all indices on all processors */
281:   PetscCall(PetscMalloc4(N, &allapp, N, &appPerm, N, &allpetsc, N, &petscPerm));
282:   PetscCallMPI(MPI_Allgatherv((void *)myapp, napp, MPIU_INT, allapp, lens, disp, MPIU_INT, comm));
283:   PetscCallMPI(MPI_Allgatherv((void *)petsc, napp, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm));
284:   PetscCall(PetscFree2(lens, disp));

286:   /* generate a list of application and PETSc node numbers */
287:   PetscCall(PetscMalloc4(N, &aomap->app, N, &aomap->appPerm, N, &aomap->petsc, N, &aomap->petscPerm));
288:   for (i = 0; i < N; i++) {
289:     appPerm[i]   = i;
290:     petscPerm[i] = i;
291:   }
292:   PetscCall(PetscSortIntWithPermutation(N, allpetsc, petscPerm));
293:   PetscCall(PetscSortIntWithPermutation(N, allapp, appPerm));
294:   /* Form sorted arrays of indices */
295:   for (i = 0; i < N; i++) {
296:     aomap->app[i]   = allapp[appPerm[i]];
297:     aomap->petsc[i] = allpetsc[petscPerm[i]];
298:   }
299:   /* Invert petscPerm[] into aomap->petscPerm[] */
300:   for (i = 0; i < N; i++) aomap->petscPerm[petscPerm[i]] = i;

302:   /* Form map between aomap->app[] and aomap->petsc[] */
303:   for (i = 0; i < N; i++) aomap->appPerm[i] = aomap->petscPerm[appPerm[i]];

305:   /* Invert appPerm[] into allapp[] */
306:   for (i = 0; i < N; i++) allapp[appPerm[i]] = i;

308:   /* Form map between aomap->petsc[] and aomap->app[] */
309:   for (i = 0; i < N; i++) aomap->petscPerm[i] = allapp[petscPerm[i]];

311:   if (PetscDefined(USE_DEBUG)) {
312:     /* Check that the permutations are complementary */
313:     for (i = 0; i < N; i++) PetscCheck(i == aomap->appPerm[aomap->petscPerm[i]], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid ordering");
314:   }
315:   /* Cleanup */
316:   if (!mypetsc) PetscCall(PetscFree(petsc));
317:   PetscCall(PetscFree4(allapp, appPerm, allpetsc, petscPerm));

319:   PetscCall(AOViewFromOptions(ao, NULL, "-ao_view"));

321:   *aoout = ao;
322:   PetscFunctionReturn(PETSC_SUCCESS);
323: }

325: /*@
326:   AOCreateMappingIS - Creates an application mapping using two index sets.

328:   Input Parameters:
329: + comm    - MPI communicator that is to share `AO`
330: . isapp   - index set that defines an ordering
331: - ispetsc - index set that defines another ordering, maybe NULL for identity `IS`

333:   Output Parameter:
334: . aoout   - the new application ordering

336:   Options Database Key:
337: . -ao_view - call `AOView()` at the conclusion of `AOCreateMappingIS()`

339:   Level: beginner

341:   Note:
342:   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.
343:   Use `AOCreateBasic()` or `AOCreateBasicIS()` if they do not have holes for better performance.

345: .seealso: [](sec_ao), [](sec_scatter), `AOCreateBasic()`, `AOCreateMapping()`, `AODestroy()`
346: @*/
347: PetscErrorCode AOCreateMappingIS(IS isapp, IS ispetsc, AO *aoout)
348: {
349:   MPI_Comm        comm;
350:   const PetscInt *mypetsc, *myapp;
351:   PetscInt        napp, npetsc;

353:   PetscFunctionBegin;
354:   PetscCall(PetscObjectGetComm((PetscObject)isapp, &comm));
355:   PetscCall(ISGetLocalSize(isapp, &napp));
356:   if (ispetsc) {
357:     PetscCall(ISGetLocalSize(ispetsc, &npetsc));
358:     PetscCheck(napp == npetsc, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local IS lengths must match");
359:     PetscCall(ISGetIndices(ispetsc, &mypetsc));
360:   } else {
361:     mypetsc = NULL;
362:   }
363:   PetscCall(ISGetIndices(isapp, &myapp));

365:   PetscCall(AOCreateMapping(comm, napp, myapp, mypetsc, aoout));

367:   PetscCall(ISRestoreIndices(isapp, &myapp));
368:   if (ispetsc) PetscCall(ISRestoreIndices(ispetsc, &mypetsc));
369:   PetscFunctionReturn(PETSC_SUCCESS);
370: }