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:   PetscFree4(aomap->app, aomap->appPerm, aomap->petsc, aomap->petscPerm);
 23:   PetscFree(aomap);
 24:   return 0;
 25: }

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

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

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

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

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

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

113: static struct _AOOps AOps = {
114:   PetscDesignatedInitializer(view, AOView_Mapping),
115:   PetscDesignatedInitializer(destroy, AODestroy_Mapping),
116:   PetscDesignatedInitializer(petsctoapplication, AOPetscToApplication_Mapping),
117:   PetscDesignatedInitializer(applicationtopetsc, AOApplicationToPetsc_Mapping),
118: };

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

123:   Not Collective

125:   Input Parameters:
126: + ao       - The `AO`
127: - index    - The application index

129:   Output Parameter:
130: . hasIndex - Flag is `PETSC_TRUE` if the index exists

132:   Level: intermediate

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

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

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

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

166:   Not Collective

168:   Input Parameters:
169: + ao       - The `AO`
170: - index    - The petsc index

172:   Output Parameter:
173: . hasIndex - Flag is `PETSC_TRUE` if the index exists

175:   Level: intermediate

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

180: .seealso: [](sec_ao), `AOMappingHasApplicationIndex()`, `AOCreateMapping()`
181: @*/
182: PetscErrorCode AOMappingHasPetscIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
183: {
184:   AO_Mapping *aomap;
185:   PetscInt   *petsc;
186:   PetscInt    low, high, mid;

190:   aomap = (AO_Mapping *)ao->data;
191:   petsc = aomap->petsc;
192:   /* Use bisection since the array is sorted */
193:   low  = 0;
194:   high = aomap->N - 1;
195:   while (low <= high) {
196:     mid = (low + high) / 2;
197:     if (idex == petsc[mid]) break;
198:     else if (idex < petsc[mid]) high = mid - 1;
199:     else low = mid + 1;
200:   }
201:   if (low > high) *hasIndex = PETSC_FALSE;
202:   else *hasIndex = PETSC_TRUE;
203:   return 0;
204: }

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

209:   Input Parameters:
210: + comm    - MPI communicator that is to share the `AO`
211: . napp    - size of integer arrays
212: . myapp   - integer array that defines an ordering
213: - mypetsc - integer array that defines another ordering (may be NULL to indicate the identity ordering)

215:   Output Parameter:
216: . aoout   - the new application mapping

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

221:   Level: beginner

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

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

241:   *aoout = NULL;
242:   AOInitializePackage();

244:   PetscHeaderCreate(ao, AO_CLASSID, "AO", "Application Ordering", "AO", comm, AODestroy, AOView);
245:   PetscNew(&aomap);
246:   PetscMemcpy(ao->ops, &AOps, sizeof(AOps));
247:   ao->data = (void *)aomap;

249:   /* transmit all lengths to all processors */
250:   MPI_Comm_size(comm, &size);
251:   MPI_Comm_rank(comm, &rank);
252:   PetscMalloc2(size, &lens, size, &disp);
253:   nnapp = napp;
254:   MPI_Allgather(&nnapp, 1, MPI_INT, lens, 1, MPI_INT, comm);
255:   N = 0;
256:   for (i = 0; i < size; i++) {
257:     disp[i] = N;
258:     N += lens[i];
259:   }
260:   aomap->N = N;
261:   ao->N    = N;
262:   ao->n    = N;

264:   /* If mypetsc is 0 then use "natural" numbering */
265:   if (!mypetsc) {
266:     start = disp[rank];
267:     PetscMalloc1(napp + 1, &petsc);
268:     for (i = 0; i < napp; i++) petsc[i] = start + i;
269:   } else {
270:     petsc = (PetscInt *)mypetsc;
271:   }

273:   /* get all indices on all processors */
274:   PetscMalloc4(N, &allapp, N, &appPerm, N, &allpetsc, N, &petscPerm);
275:   MPI_Allgatherv((void *)myapp, napp, MPIU_INT, allapp, lens, disp, MPIU_INT, comm);
276:   MPI_Allgatherv((void *)petsc, napp, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm);
277:   PetscFree2(lens, disp);

279:   /* generate a list of application and PETSc node numbers */
280:   PetscMalloc4(N, &aomap->app, N, &aomap->appPerm, N, &aomap->petsc, N, &aomap->petscPerm);
281:   for (i = 0; i < N; i++) {
282:     appPerm[i]   = i;
283:     petscPerm[i] = i;
284:   }
285:   PetscSortIntWithPermutation(N, allpetsc, petscPerm);
286:   PetscSortIntWithPermutation(N, allapp, appPerm);
287:   /* Form sorted arrays of indices */
288:   for (i = 0; i < N; i++) {
289:     aomap->app[i]   = allapp[appPerm[i]];
290:     aomap->petsc[i] = allpetsc[petscPerm[i]];
291:   }
292:   /* Invert petscPerm[] into aomap->petscPerm[] */
293:   for (i = 0; i < N; i++) aomap->petscPerm[petscPerm[i]] = i;

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

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

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

304:   if (PetscDefined(USE_DEBUG)) {
305:     /* Check that the permutations are complementary */
307:   }
308:   /* Cleanup */
309:   if (!mypetsc) PetscFree(petsc);
310:   PetscFree4(allapp, appPerm, allpetsc, petscPerm);

312:   AOViewFromOptions(ao, NULL, "-ao_view");

314:   *aoout = ao;
315:   return 0;
316: }

318: /*@
319:   AOCreateMappingIS - Creates an application mapping using two index sets.

321:   Input Parameters:
322: + comm    - MPI communicator that is to share `AO`
323: . isapp   - index set that defines an ordering
324: - ispetsc - index set that defines another ordering, maybe NULL for identity `IS`

326:   Output Parameter:
327: . aoout   - the new application ordering

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

332:   Level: beginner

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

338: .seealso: [](sec_ao), [](sec_scatter), `AOCreateBasic()`, `AOCreateMapping()`, `AODestroy()`
339: @*/
340: PetscErrorCode AOCreateMappingIS(IS isapp, IS ispetsc, AO *aoout)
341: {
342:   MPI_Comm        comm;
343:   const PetscInt *mypetsc, *myapp;
344:   PetscInt        napp, npetsc;

346:   PetscObjectGetComm((PetscObject)isapp, &comm);
347:   ISGetLocalSize(isapp, &napp);
348:   if (ispetsc) {
349:     ISGetLocalSize(ispetsc, &npetsc);
351:     ISGetIndices(ispetsc, &mypetsc);
352:   } else {
353:     mypetsc = NULL;
354:   }
355:   ISGetIndices(isapp, &myapp);

357:   AOCreateMapping(comm, napp, myapp, mypetsc, aoout);

359:   ISRestoreIndices(isapp, &myapp);
360:   if (ispetsc) ISRestoreIndices(ispetsc, &mypetsc);
361:   return 0;
362: }