Actual source code: aomapping.c

petsc-3.7.2 2016-06-05
Report Typos and Errors
  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>          /*I  "petscao.h" I*/

 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;

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

 26:   PetscFree4(aomap->app,aomap->appPerm,aomap->petsc,aomap->petscPerm);
 27:   PetscFree(aomap);
 28:   return(0);
 29: }

 33: PetscErrorCode AOView_Mapping(AO ao, PetscViewer viewer)
 34: {
 35:   AO_Mapping     *aomap = (AO_Mapping*) ao->data;
 36:   PetscMPIInt    rank;
 37:   PetscInt       i;
 38:   PetscBool      iascii;

 42:   MPI_Comm_rank(PetscObjectComm((PetscObject)ao), &rank);
 43:   if (rank) return(0);
 44:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
 45:   if (iascii) {
 46:     PetscViewerASCIIPrintf(viewer, "Number of elements in ordering %D\n", aomap->N);
 47:     PetscViewerASCIIPrintf(viewer, "   App.   PETSc\n");
 48:     for (i = 0; i < aomap->N; i++) {
 49:       PetscViewerASCIIPrintf(viewer, "%D   %D    %D\n", i, aomap->app[i], aomap->petsc[aomap->appPerm[i]]);
 50:     }
 51:   }
 52:   return(0);
 53: }

 57: PetscErrorCode AOPetscToApplication_Mapping(AO ao, PetscInt n, PetscInt *ia)
 58: {
 59:   AO_Mapping *aomap = (AO_Mapping*) ao->data;
 60:   PetscInt   *app   = aomap->app;
 61:   PetscInt   *petsc = aomap->petsc;
 62:   PetscInt   *perm  = aomap->petscPerm;
 63:   PetscInt   N      = aomap->N;
 64:   PetscInt   low, high, mid=0;
 65:   PetscInt   idex;
 66:   PetscInt   i;

 68:   /* It would be possible to use a single bisection search, which
 69:      recursively divided the indices to be converted, and searched
 70:      partitions which contained an index. This would result in
 71:      better running times if indices are clustered.
 72:   */
 74:   for (i = 0; i < n; i++) {
 75:     idex = ia[i];
 76:     if (idex < 0) continue;
 77:     /* Use bisection since the array is sorted */
 78:     low  = 0;
 79:     high = N - 1;
 80:     while (low <= high) {
 81:       mid = (low + high)/2;
 82:       if (idex == petsc[mid]) break;
 83:       else if (idex < petsc[mid]) high = mid - 1;
 84:       else low = mid + 1;
 85:     }
 86:     if (low > high) ia[i] = -1;
 87:     else            ia[i] = app[perm[mid]];
 88:   }
 89:   return(0);
 90: }

 94: PetscErrorCode AOApplicationToPetsc_Mapping(AO ao, PetscInt n, PetscInt *ia)
 95: {
 96:   AO_Mapping *aomap = (AO_Mapping*) ao->data;
 97:   PetscInt   *app   = aomap->app;
 98:   PetscInt   *petsc = aomap->petsc;
 99:   PetscInt   *perm  = aomap->appPerm;
100:   PetscInt   N      = aomap->N;
101:   PetscInt   low, high, mid=0;
102:   PetscInt   idex;
103:   PetscInt   i;

105:   /* It would be possible to use a single bisection search, which
106:      recursively divided the indices to be converted, and searched
107:      partitions which contained an index. This would result in
108:      better running times if indices are clustered.
109:   */
111:   for (i = 0; i < n; i++) {
112:     idex = ia[i];
113:     if (idex < 0) continue;
114:     /* Use bisection since the array is sorted */
115:     low  = 0;
116:     high = N - 1;
117:     while (low <= high) {
118:       mid = (low + high)/2;
119:       if (idex == app[mid]) break;
120:       else if (idex < app[mid]) high = mid - 1;
121:       else low = mid + 1;
122:     }
123:     if (low > high) ia[i] = -1;
124:     else            ia[i] = petsc[perm[mid]];
125:   }
126:   return(0);
127: }

129: static struct _AOOps AOps = {AOView_Mapping,
130:                              AODestroy_Mapping,
131:                              AOPetscToApplication_Mapping,
132:                              AOApplicationToPetsc_Mapping,
133:                              NULL,
134:                              NULL,
135:                              NULL,
136:                              NULL};

140: /*@C
141:   AOMappingHasApplicationIndex - Searches for the supplied application index.

143:   Input Parameters:
144: + ao       - The AOMapping
145: - index    - The application index

147:   Output Parameter:
148: . hasIndex - Flag is PETSC_TRUE if the index exists

150:   Level: intermediate

152: .keywords: AO, index
153: .seealso: AOMappingHasPetscIndex(), AOCreateMapping()
154: @*/
155: PetscErrorCode  AOMappingHasApplicationIndex(AO ao, PetscInt idex, PetscBool  *hasIndex)
156: {
157:   AO_Mapping *aomap;
158:   PetscInt   *app;
159:   PetscInt   low, high, mid;

164:   aomap = (AO_Mapping*) ao->data;
165:   app   = aomap->app;
166:   /* Use bisection since the array is sorted */
167:   low  = 0;
168:   high = aomap->N - 1;
169:   while (low <= high) {
170:     mid = (low + high)/2;
171:     if (idex == app[mid]) break;
172:     else if (idex < app[mid]) high = mid - 1;
173:     else low = mid + 1;
174:   }
175:   if (low > high) *hasIndex = PETSC_FALSE;
176:   else *hasIndex = PETSC_TRUE;
177:   return(0);
178: }

182: /*@
183:   AOMappingHasPetscIndex - Searches for the supplied petsc index.

185:   Input Parameters:
186: + ao       - The AOMapping
187: - index    - The petsc index

189:   Output Parameter:
190: . hasIndex - Flag is PETSC_TRUE if the index exists

192:   Level: intermediate

194: .keywords: AO, index
195: .seealso: AOMappingHasApplicationIndex(), AOCreateMapping()
196: @*/
197: PetscErrorCode  AOMappingHasPetscIndex(AO ao, PetscInt idex, PetscBool  *hasIndex)
198: {
199:   AO_Mapping *aomap;
200:   PetscInt   *petsc;
201:   PetscInt   low, high, mid;

206:   aomap = (AO_Mapping*) ao->data;
207:   petsc = aomap->petsc;
208:   /* Use bisection since the array is sorted */
209:   low  = 0;
210:   high = aomap->N - 1;
211:   while (low <= high) {
212:     mid = (low + high)/2;
213:     if (idex == petsc[mid]) break;
214:     else if (idex < petsc[mid]) high = mid - 1;
215:     else low = mid + 1;
216:   }
217:   if (low > high) *hasIndex = PETSC_FALSE;
218:   else *hasIndex = PETSC_TRUE;
219:   return(0);
220: }

224: /*@C
225:   AOCreateMapping - Creates a basic application mapping using two integer arrays.

227:   Input Parameters:
228: + comm    - MPI communicator that is to share AO
229: . napp    - size of integer arrays
230: . myapp   - integer array that defines an ordering
231: - mypetsc - integer array that defines another ordering (may be NULL to indicate the identity ordering)

233:   Output Parameter:
234: . aoout   - the new application mapping

236:   Options Database Key:
237: . -ao_view : call AOView() at the conclusion of AOCreateMapping()

239:   Level: beginner

241:     Notes: 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.
242:        Use AOCreateBasic() or AOCreateBasicIS() if they do not have holes for better performance.

244: .keywords: AO, create
245: .seealso: AOCreateBasic(), AOCreateBasic(), AOCreateMappingIS(), AODestroy()
246: @*/
247: PetscErrorCode  AOCreateMapping(MPI_Comm comm,PetscInt napp,const PetscInt myapp[],const PetscInt mypetsc[],AO *aoout)
248: {
249:   AO             ao;
250:   AO_Mapping     *aomap;
251:   PetscInt       *allpetsc,  *allapp;
252:   PetscInt       *petscPerm, *appPerm;
253:   PetscInt       *petsc;
254:   PetscMPIInt    size, rank,*lens, *disp,nnapp;
255:   PetscInt       N, start;
256:   PetscInt       i;

261:   *aoout = 0;
262:   AOInitializePackage();

264:   PetscHeaderCreate(ao, AO_CLASSID, "AO", "Application Ordering", "AO", comm, AODestroy, AOView);
265:   PetscNewLog(ao,&aomap);
266:   PetscMemcpy(ao->ops, &AOps, sizeof(AOps));
267:   ao->data = (void*) aomap;

269:   /* transmit all lengths to all processors */
270:   MPI_Comm_size(comm, &size);
271:   MPI_Comm_rank(comm, &rank);
272:   PetscMalloc2(size, &lens,size,&disp);
273:   nnapp = napp;
274:   MPI_Allgather(&nnapp, 1, MPI_INT, lens, 1, MPI_INT, comm);
275:   N     = 0;
276:   for (i = 0; i < size; i++) {
277:     disp[i] = N;
278:     N      += lens[i];
279:   }
280:   aomap->N = N;
281:   ao->N    = N;
282:   ao->n    = N;

284:   /* If mypetsc is 0 then use "natural" numbering */
285:   if (!mypetsc) {
286:     start = disp[rank];
287:     PetscMalloc1(napp+1, &petsc);
288:     for (i = 0; i < napp; i++) petsc[i] = start + i;
289:   } else {
290:     petsc = (PetscInt*)mypetsc;
291:   }

293:   /* get all indices on all processors */
294:   PetscMalloc4(N, &allapp,N,&appPerm,N,&allpetsc,N,&petscPerm);
295:   MPI_Allgatherv((void*)myapp, napp, MPIU_INT, allapp,   lens, disp, MPIU_INT, comm);
296:   MPI_Allgatherv((void*)petsc, napp, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm);
297:   PetscFree2(lens,disp);

299:   /* generate a list of application and PETSc node numbers */
300:   PetscMalloc4(N, &aomap->app,N,&aomap->appPerm,N,&aomap->petsc,N,&aomap->petscPerm);
301:   PetscLogObjectMemory((PetscObject)ao, 4*N * sizeof(PetscInt));
302:   for (i = 0; i < N; i++) {
303:     appPerm[i]   = i;
304:     petscPerm[i] = i;
305:   }
306:   PetscSortIntWithPermutation(N, allpetsc, petscPerm);
307:   PetscSortIntWithPermutation(N, allapp,   appPerm);
308:   /* Form sorted arrays of indices */
309:   for (i = 0; i < N; i++) {
310:     aomap->app[i]   = allapp[appPerm[i]];
311:     aomap->petsc[i] = allpetsc[petscPerm[i]];
312:   }
313:   /* Invert petscPerm[] into aomap->petscPerm[] */
314:   for (i = 0; i < N; i++) aomap->petscPerm[petscPerm[i]] = i;

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

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

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

325: #if defined(PETSC_USE_DEBUG)
326:   /* Check that the permutations are complementary */
327:   for (i = 0; i < N; i++) {
328:     if (i != aomap->appPerm[aomap->petscPerm[i]]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid ordering");
329:   }
330: #endif
331:   /* Cleanup */
332:   if (!mypetsc) {
333:     PetscFree(petsc);
334:   }
335:   PetscFree4(allapp,appPerm,allpetsc,petscPerm);

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

339:   *aoout = ao;
340:   return(0);
341: }

345: /*@C
346:   AOCreateMappingIS - Creates a basic application ordering using two index sets.

348:   Input Parameters:
349: + comm    - MPI communicator that is to share AO
350: . isapp   - index set that defines an ordering
351: - ispetsc - index set that defines another ordering, maybe NULL for identity IS

353:   Output Parameter:
354: . aoout   - the new application ordering

356:   Options Database Key:
357: . -ao_view : call AOView() at the conclusion of AOCreateMappingIS()

359:   Level: beginner

361:     Notes: 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.
362:        Use AOCreateBasic() or AOCreateBasicIS() if they do not have holes for better performance.

364: .keywords: AO, create
365: .seealso: AOCreateBasic(), AOCreateMapping(), AODestroy()
366: @*/
367: PetscErrorCode  AOCreateMappingIS(IS isapp, IS ispetsc, AO *aoout)
368: {
369:   MPI_Comm       comm;
370:   const PetscInt *mypetsc, *myapp;
371:   PetscInt       napp, npetsc;

375:   PetscObjectGetComm((PetscObject) isapp, &comm);
376:   ISGetLocalSize(isapp, &napp);
377:   if (ispetsc) {
378:     ISGetLocalSize(ispetsc, &npetsc);
379:     if (napp != npetsc) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Local IS lengths must match");
380:     ISGetIndices(ispetsc, &mypetsc);
381:   } else {
382:     mypetsc = NULL;
383:   }
384:   ISGetIndices(isapp, &myapp);

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

388:   ISRestoreIndices(isapp, &myapp);
389:   if (ispetsc) {
390:     ISRestoreIndices(ispetsc, &mypetsc);
391:   }
392:   return(0);
393: }