Actual source code: aomapping.c

petsc-3.11.3 2019-06-26
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>

 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;

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

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

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

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

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

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

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

121: static struct _AOOps AOps = {AOView_Mapping,
122:                              AODestroy_Mapping,
123:                              AOPetscToApplication_Mapping,
124:                              AOApplicationToPetsc_Mapping,
125:                              NULL,
126:                              NULL,
127:                              NULL,
128:                              NULL};

130: /*@C
131:   AOMappingHasApplicationIndex - Searches for the supplied application index.

133:   Input Parameters:
134: + ao       - The AOMapping
135: - index    - The application index

137:   Output Parameter:
138: . hasIndex - Flag is PETSC_TRUE if the index exists

140:   Level: intermediate

142: .keywords: AO, index
143: .seealso: AOMappingHasPetscIndex(), AOCreateMapping()
144: @*/
145: PetscErrorCode  AOMappingHasApplicationIndex(AO ao, PetscInt idex, PetscBool  *hasIndex)
146: {
147:   AO_Mapping *aomap;
148:   PetscInt   *app;
149:   PetscInt   low, high, mid;

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

170: /*@
171:   AOMappingHasPetscIndex - Searches for the supplied petsc index.

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

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

180:   Level: intermediate

182: .keywords: AO, index
183: .seealso: AOMappingHasApplicationIndex(), AOCreateMapping()
184: @*/
185: PetscErrorCode  AOMappingHasPetscIndex(AO ao, PetscInt idex, PetscBool  *hasIndex)
186: {
187:   AO_Mapping *aomap;
188:   PetscInt   *petsc;
189:   PetscInt   low, high, mid;

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

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

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

219:   Output Parameter:
220: . aoout   - the new application mapping

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

225:   Level: beginner

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

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

248:   *aoout = 0;
249:   AOInitializePackage();

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

256:   /* transmit all lengths to all processors */
257:   MPI_Comm_size(comm, &size);
258:   MPI_Comm_rank(comm, &rank);
259:   PetscMalloc2(size, &lens,size,&disp);
260:   nnapp = napp;
261:   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:     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:   PetscMalloc4(N, &allapp,N,&appPerm,N,&allpetsc,N,&petscPerm);
282:   MPI_Allgatherv((void*)myapp, napp, MPIU_INT, allapp,   lens, disp, MPIU_INT, comm);
283:   MPI_Allgatherv((void*)petsc, napp, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm);
284:   PetscFree2(lens,disp);

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

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

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

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

312: #if defined(PETSC_USE_DEBUG)
313:   /* Check that the permutations are complementary */
314:   for (i = 0; i < N; i++) {
315:     if (i != aomap->appPerm[aomap->petscPerm[i]]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid ordering");
316:   }
317: #endif
318:   /* Cleanup */
319:   if (!mypetsc) {
320:     PetscFree(petsc);
321:   }
322:   PetscFree4(allapp,appPerm,allpetsc,petscPerm);

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

326:   *aoout = ao;
327:   return(0);
328: }

330: /*@
331:   AOCreateMappingIS - Creates a basic application ordering using two index sets.

333:   Input Parameters:
334: + comm    - MPI communicator that is to share AO
335: . isapp   - index set that defines an ordering
336: - ispetsc - index set that defines another ordering, maybe NULL for identity IS

338:   Output Parameter:
339: . aoout   - the new application ordering

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

344:   Level: beginner

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

350: .keywords: AO, create
351: .seealso: AOCreateBasic(), AOCreateMapping(), AODestroy()
352: @*/
353: PetscErrorCode  AOCreateMappingIS(IS isapp, IS ispetsc, AO *aoout)
354: {
355:   MPI_Comm       comm;
356:   const PetscInt *mypetsc, *myapp;
357:   PetscInt       napp, npetsc;

361:   PetscObjectGetComm((PetscObject) isapp, &comm);
362:   ISGetLocalSize(isapp, &napp);
363:   if (ispetsc) {
364:     ISGetLocalSize(ispetsc, &npetsc);
365:     if (napp != npetsc) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Local IS lengths must match");
366:     ISGetIndices(ispetsc, &mypetsc);
367:   } else {
368:     mypetsc = NULL;
369:   }
370:   ISGetIndices(isapp, &myapp);

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

374:   ISRestoreIndices(isapp, &myapp);
375:   if (ispetsc) {
376:     ISRestoreIndices(ispetsc, &mypetsc);
377:   }
378:   return(0);
379: }