Actual source code: aobasic.c

petsc-3.4.4 2014-03-13
  2: /*
  3:     The most basic AO application ordering routines. These store the
  4:   entire orderings on each processor.
  5: */

  7: #include <../src/vec/is/ao/aoimpl.h>          /*I  "petscao.h"   I*/

  9: typedef struct {
 10:   PetscInt *app;     /* app[i] is the partner for the ith PETSc slot */
 11:   PetscInt *petsc;   /* petsc[j] is the partner for the jth app slot */
 12: } AO_Basic;

 14: /*
 15:        All processors have the same data so processor 1 prints it
 16: */
 19: PetscErrorCode AOView_Basic(AO ao,PetscViewer viewer)
 20: {
 22:   PetscMPIInt    rank;
 23:   PetscInt       i;
 24:   AO_Basic       *aobasic = (AO_Basic*)ao->data;
 25:   PetscBool      iascii;

 28:   MPI_Comm_rank(PetscObjectComm((PetscObject)ao),&rank);
 29:   if (!rank) {
 30:     PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
 31:     if (iascii) {
 32:       PetscViewerASCIIPrintf(viewer,"Number of elements in ordering %D\n",ao->N);
 33:       PetscViewerASCIIPrintf(viewer,  "PETSc->App  App->PETSc\n");
 34:       for (i=0; i<ao->N; i++) {
 35:         PetscViewerASCIIPrintf(viewer,"%3D  %3D    %3D  %3D\n",i,aobasic->app[i],i,aobasic->petsc[i]);
 36:       }
 37:     }
 38:   }
 39:   PetscViewerFlush(viewer);
 40:   return(0);
 41: }

 45: PetscErrorCode AODestroy_Basic(AO ao)
 46: {
 47:   AO_Basic       *aobasic = (AO_Basic*)ao->data;

 51:   PetscFree2(aobasic->app,aobasic->petsc);
 52:   PetscFree(aobasic);
 53:   return(0);
 54: }

 58: PetscErrorCode AOBasicGetIndices_Private(AO ao,PetscInt **app,PetscInt **petsc)
 59: {
 60:   AO_Basic *basic = (AO_Basic*)ao->data;

 63:   if (app)   *app   = basic->app;
 64:   if (petsc) *petsc = basic->petsc;
 65:   return(0);
 66: }

 70: PetscErrorCode AOPetscToApplication_Basic(AO ao,PetscInt n,PetscInt *ia)
 71: {
 72:   PetscInt i,N=ao->N;
 73:   AO_Basic *aobasic = (AO_Basic*)ao->data;

 76:   for (i=0; i<n; i++) {
 77:     if (ia[i] >= 0 && ia[i] < N) {
 78:       ia[i] = aobasic->app[ia[i]];
 79:     } else {
 80:       ia[i] = -1;
 81:     }
 82:   }
 83:   return(0);
 84: }

 88: PetscErrorCode AOApplicationToPetsc_Basic(AO ao,PetscInt n,PetscInt *ia)
 89: {
 90:   PetscInt i,N=ao->N;
 91:   AO_Basic *aobasic = (AO_Basic*)ao->data;

 94:   for (i=0; i<n; i++) {
 95:     if (ia[i] >= 0 && ia[i] < N) {
 96:       ia[i] = aobasic->petsc[ia[i]];
 97:     } else {
 98:       ia[i] = -1;
 99:     }
100:   }
101:   return(0);
102: }

106: PetscErrorCode AOPetscToApplicationPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array)
107: {
108:   AO_Basic       *aobasic = (AO_Basic*) ao->data;
109:   PetscInt       *temp;
110:   PetscInt       i, j;

114:   PetscMalloc(ao->N*block * sizeof(PetscInt), &temp);
115:   for (i = 0; i < ao->N; i++) {
116:     for (j = 0; j < block; j++) temp[i*block+j] = array[aobasic->petsc[i]*block+j];
117:   }
118:   PetscMemcpy(array, temp, ao->N*block * sizeof(PetscInt));
119:   PetscFree(temp);
120:   return(0);
121: }

125: PetscErrorCode AOApplicationToPetscPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array)
126: {
127:   AO_Basic       *aobasic = (AO_Basic*) ao->data;
128:   PetscInt       *temp;
129:   PetscInt       i, j;

133:   PetscMalloc(ao->N*block * sizeof(PetscInt), &temp);
134:   for (i = 0; i < ao->N; i++) {
135:     for (j = 0; j < block; j++) temp[i*block+j] = array[aobasic->app[i]*block+j];
136:   }
137:   PetscMemcpy(array, temp, ao->N*block * sizeof(PetscInt));
138:   PetscFree(temp);
139:   return(0);
140: }

144: PetscErrorCode AOPetscToApplicationPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array)
145: {
146:   AO_Basic       *aobasic = (AO_Basic*) ao->data;
147:   PetscReal      *temp;
148:   PetscInt       i, j;

152:   PetscMalloc(ao->N*block * sizeof(PetscReal), &temp);
153:   for (i = 0; i < ao->N; i++) {
154:     for (j = 0; j < block; j++) temp[i*block+j] = array[aobasic->petsc[i]*block+j];
155:   }
156:   PetscMemcpy(array, temp, ao->N*block * sizeof(PetscReal));
157:   PetscFree(temp);
158:   return(0);
159: }

163: PetscErrorCode AOApplicationToPetscPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array)
164: {
165:   AO_Basic       *aobasic = (AO_Basic*) ao->data;
166:   PetscReal      *temp;
167:   PetscInt       i, j;

171:   PetscMalloc(ao->N*block * sizeof(PetscReal), &temp);
172:   for (i = 0; i < ao->N; i++) {
173:     for (j = 0; j < block; j++) temp[i*block+j] = array[aobasic->app[i]*block+j];
174:   }
175:   PetscMemcpy(array, temp, ao->N*block * sizeof(PetscReal));
176:   PetscFree(temp);
177:   return(0);
178: }

180: static struct _AOOps AOOps_Basic = {
181:   AOView_Basic,
182:   AODestroy_Basic,
183:   AOPetscToApplication_Basic,
184:   AOApplicationToPetsc_Basic,
185:   AOPetscToApplicationPermuteInt_Basic,
186:   AOApplicationToPetscPermuteInt_Basic,
187:   AOPetscToApplicationPermuteReal_Basic,
188:   AOApplicationToPetscPermuteReal_Basic
189: };

193: PETSC_EXTERN PetscErrorCode AOCreate_Basic(AO ao)
194: {
195:   AO_Basic       *aobasic;
196:   PetscMPIInt    size,rank,count,*lens,*disp;
197:   PetscInt       napp,*allpetsc,*allapp,ip,ia,N,i,*petsc=NULL,start;
199:   IS             isapp=ao->isapp,ispetsc=ao->ispetsc;
200:   MPI_Comm       comm;
201:   const PetscInt *myapp,*mypetsc=NULL;

204:   /* create special struct aobasic */
205:   PetscNewLog(ao, AO_Basic, &aobasic);
206:   ao->data = (void*) aobasic;
207:   PetscMemcpy(ao->ops,&AOOps_Basic,sizeof(struct _AOOps));
208:   PetscObjectChangeTypeName((PetscObject)ao,AOBASIC);

210:   ISGetLocalSize(isapp,&napp);
211:   ISGetIndices(isapp,&myapp);

213:   PetscMPIIntCast(napp,&count);

215:   /* transmit all lengths to all processors */
216:   PetscObjectGetComm((PetscObject)isapp,&comm);
217:   MPI_Comm_size(comm, &size);
218:   MPI_Comm_rank(comm, &rank);
219:   PetscMalloc2(size,PetscMPIInt, &lens,size,PetscMPIInt,&disp);
220:   MPI_Allgather(&count, 1, MPI_INT, lens, 1, MPI_INT, comm);
221:   N    =  0;
222:   for (i = 0; i < size; i++) {
223:     PetscMPIIntCast(N,disp+i); /* = sum(lens[j]), j< i */
224:     N   += lens[i];
225:   }
226:   ao->N = N;
227:   ao->n = N;

229:   /* If mypetsc is 0 then use "natural" numbering */
230:   if (napp) {
231:     if (!ispetsc) {
232:       start = disp[rank];
233:       PetscMalloc((napp+1) * sizeof(PetscInt), &petsc);
234:       for (i=0; i<napp; i++) petsc[i] = start + i;
235:     } else {
236:       ISGetIndices(ispetsc,&mypetsc);
237:       petsc = (PetscInt*)mypetsc;
238:     }
239:   }

241:   /* get all indices on all processors */
242:   PetscMalloc2(N,PetscInt,&allpetsc,N,PetscInt,&allapp);
243:   MPI_Allgatherv(petsc, count, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm);
244:   MPI_Allgatherv((void*)myapp, count, MPIU_INT, allapp, lens, disp, MPIU_INT, comm);
245:   PetscFree2(lens,disp);

247: #if defined(PETSC_USE_DEBUG)
248:   {
249:     PetscInt *sorted;
250:     PetscMalloc(N*sizeof(PetscInt),&sorted);

252:     PetscMemcpy(sorted,allpetsc,N*sizeof(PetscInt));
253:     PetscSortInt(N,sorted);
254:     for (i=0; i<N; i++) {
255:       if (sorted[i] != i) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PETSc ordering requires a permutation of numbers 0 to N-1\n it is missing %D has %D",i,sorted[i]);
256:     }

258:     PetscMemcpy(sorted,allapp,N*sizeof(PetscInt));
259:     PetscSortInt(N,sorted);
260:     for (i=0; i<N; i++) {
261:       if (sorted[i] != i) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Application ordering requires a permutation of numbers 0 to N-1\n it is missing %D has %D",i,sorted[i]);
262:     }

264:     PetscFree(sorted);
265:   }
266: #endif

268:   /* generate a list of application and PETSc node numbers */
269:   PetscMalloc2(N,PetscInt, &aobasic->app,N,PetscInt,&aobasic->petsc);
270:   PetscLogObjectMemory(ao,2*N*sizeof(PetscInt));
271:   PetscMemzero(aobasic->app, N*sizeof(PetscInt));
272:   PetscMemzero(aobasic->petsc, N*sizeof(PetscInt));
273:   for (i = 0; i < N; i++) {
274:     ip = allpetsc[i];
275:     ia = allapp[i];
276:     /* check there are no duplicates */
277:     if (aobasic->app[ip]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Duplicate in PETSc ordering at position %d. Already mapped to %d, not %d.", i, aobasic->app[ip]-1, ia);
278:     aobasic->app[ip] = ia + 1;
279:     if (aobasic->petsc[ia]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Duplicate in Application ordering at position %d. Already mapped to %d, not %d.", i, aobasic->petsc[ia]-1, ip);
280:     aobasic->petsc[ia] = ip + 1;
281:   }
282:   if (napp && !mypetsc) {
283:     PetscFree(petsc);
284:   }
285:   PetscFree2(allpetsc,allapp);
286:   /* shift indices down by one */
287:   for (i = 0; i < N; i++) {
288:     aobasic->app[i]--;
289:     aobasic->petsc[i]--;
290:   }

292:   ISRestoreIndices(isapp,&myapp);
293:   if (napp) {
294:     if (ispetsc) {
295:       ISRestoreIndices(ispetsc,&mypetsc);
296:     } else {
297:       PetscFree(petsc);
298:     }
299:   }
300:   return(0);
301: }

305: /*@C
306:    AOCreateBasic - Creates a basic application ordering using two integer arrays.

308:    Collective on MPI_Comm

310:    Input Parameters:
311: +  comm - MPI communicator that is to share AO
312: .  napp - size of integer arrays
313: .  myapp - integer array that defines an ordering
314: -  mypetsc - integer array that defines another ordering (may be NULL to
315:              indicate the natural ordering, that is 0,1,2,3,...)

317:    Output Parameter:
318: .  aoout - the new application ordering

320:    Level: beginner

322:     Notes: the arrays myapp and mypetsc must contain the all the integers 0 to napp-1 with no duplicates; that is there cannot be any "holes"
323:            in the indices. Use AOCreateMapping() or AOCreateMappingIS() if you wish to have "holes" in the indices.

325: .keywords: AO, create

327: .seealso: AOCreateBasicIS(), AODestroy(), AOPetscToApplication(), AOApplicationToPetsc()
328: @*/
329: PetscErrorCode  AOCreateBasic(MPI_Comm comm,PetscInt napp,const PetscInt myapp[],const PetscInt mypetsc[],AO *aoout)
330: {
332:   IS             isapp,ispetsc;
333:   const PetscInt *app=myapp,*petsc=mypetsc;

336:   ISCreateGeneral(comm,napp,app,PETSC_USE_POINTER,&isapp);
337:   if (mypetsc) {
338:     ISCreateGeneral(comm,napp,petsc,PETSC_USE_POINTER,&ispetsc);
339:   } else {
340:     ispetsc = NULL;
341:   }
342:   AOCreateBasicIS(isapp,ispetsc,aoout);
343:   ISDestroy(&isapp);
344:   if (mypetsc) {
345:     ISDestroy(&ispetsc);
346:   }
347:   return(0);
348: }

352: /*@C
353:    AOCreateBasicIS - Creates a basic application ordering using two index sets.

355:    Collective on IS

357:    Input Parameters:
358: +  isapp - index set that defines an ordering
359: -  ispetsc - index set that defines another ordering (may be NULL to use the
360:              natural ordering)

362:    Output Parameter:
363: .  aoout - the new application ordering

365:    Level: beginner

367:     Notes: the index sets isapp and ispetsc must contain the all the integers 0 to napp-1 (where napp is the length of the index sets) with no duplicates;
368:            that is there cannot be any "holes"

370: .keywords: AO, create

372: .seealso: AOCreateBasic(),  AODestroy()
373: @*/
374: PetscErrorCode AOCreateBasicIS(IS isapp,IS ispetsc,AO *aoout)
375: {
377:   MPI_Comm       comm;
378:   AO             ao;

381:   PetscObjectGetComm((PetscObject)isapp,&comm);
382:   AOCreate(comm,&ao);
383:   AOSetIS(ao,isapp,ispetsc);
384:   AOSetType(ao,AOBASIC);
385:   AOViewFromOptions(ao,NULL,"-ao_view");
386:   *aoout = ao;
387:   return(0);
388: }