Actual source code: aomapping.c

petsc-3.4.4 2014-03-13
  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);

 45:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_SELF;

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

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

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

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

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

132: static struct _AOOps AOps = {AOView_Mapping,
133:                              AODestroy_Mapping,
134:                              AOPetscToApplication_Mapping,
135:                              AOApplicationToPetsc_Mapping,
136:                              NULL,
137:                              NULL,
138:                              NULL,
139:                              NULL};

143: /*@C
144:   AOMappingHasApplicationIndex - Searches for the supplied application index.

146:   Input Parameters:
147: + ao       - The AOMapping
148: - index    - The application index

150:   Output Parameter:
151: . hasIndex - Flag is PETSC_TRUE if the index exists

153:   Level: intermediate

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

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

185: /*@
186:   AOMappingHasPetscIndex - Searches for the supplied petsc index.

188:   Input Parameters:
189: + ao       - The AOMapping
190: - index    - The petsc index

192:   Output Parameter:
193: . hasIndex - Flag is PETSC_TRUE if the index exists

195:   Level: intermediate

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

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

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

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

236:   Output Parameter:
237: . aoout   - the new application mapping

239:   Options Database Key:
240: $ -ao_view : call AOView() at the conclusion of AOCreateMapping()

242:   Level: beginner

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

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

264:   *aoout = 0;
265: #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
266:   AOInitializePackage();
267: #endif

269:   PetscHeaderCreate(ao, _p_AO, struct _AOOps, AO_CLASSID, "AO", "Application Ordering", "AO", comm, AODestroy, AOView);
270:   PetscNewLog(ao, AO_Mapping, &aomap);
271:   PetscMemcpy(ao->ops, &AOps, sizeof(AOps));
272:   ao->data = (void*) aomap;

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

289:   /* If mypetsc is 0 then use "natural" numbering */
290:   if (!mypetsc) {
291:     start = disp[rank];
292:     PetscMalloc((napp+1) * sizeof(PetscInt), &petsc);
293:     for (i = 0; i < napp; i++) petsc[i] = start + i;
294:   } else {
295:     petsc = (PetscInt*)mypetsc;
296:   }

298:   /* get all indices on all processors */
299:   PetscMalloc4(N,PetscInt, &allapp,N,PetscInt,&appPerm,N,PetscInt,&allpetsc,N,PetscInt,&petscPerm);
300:   MPI_Allgatherv((void*)myapp, napp, MPIU_INT, allapp,   lens, disp, MPIU_INT, comm);
301:   MPI_Allgatherv((void*)petsc, napp, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm);
302:   PetscFree2(lens,disp);

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

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

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

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

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

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

344:   *aoout = ao;
345:   return(0);
346: }

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

353:   Input Parameters:
354: + comm    - MPI communicator that is to share AO
355: . isapp   - index set that defines an ordering
356: - ispetsc - index set that defines another ordering, maybe NULL for identity IS

358:   Output Parameter:
359: . aoout   - the new application ordering

361:   Options Database Key:
362: $ -ao_view : call AOView() at the conclusion of AOCreateMappingIS()

364:   Level: beginner

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

369: .keywords: AO, create
370: .seealso: AOCreateBasic(), AOCreateMapping(), AODestroy()
371: @*/
372: PetscErrorCode  AOCreateMappingIS(IS isapp, IS ispetsc, AO *aoout)
373: {
374:   MPI_Comm       comm;
375:   const PetscInt *mypetsc, *myapp;
376:   PetscInt       napp, npetsc;

380:   PetscObjectGetComm((PetscObject) isapp, &comm);
381:   ISGetLocalSize(isapp, &napp);
382:   if (ispetsc) {
383:     ISGetLocalSize(ispetsc, &npetsc);
384:     if (napp != npetsc) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Local IS lengths must match");
385:     ISGetIndices(ispetsc, &mypetsc);
386:   } else {
387:     mypetsc = NULL;
388:   }
389:   ISGetIndices(isapp, &myapp);

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

393:   ISRestoreIndices(isapp, &myapp);
394:   if (ispetsc) {
395:     ISRestoreIndices(ispetsc, &mypetsc);
396:   }
397:   return(0);
398: }