Actual source code: aobasic.c


  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>

  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: */
 17: PetscErrorCode AOView_Basic(AO ao,PetscViewer viewer)
 18: {
 19:   PetscMPIInt    rank;
 20:   PetscInt       i;
 21:   AO_Basic       *aobasic = (AO_Basic*)ao->data;
 22:   PetscBool      iascii;

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

 39: PetscErrorCode AODestroy_Basic(AO ao)
 40: {
 41:   AO_Basic       *aobasic = (AO_Basic*)ao->data;

 43:   PetscFree2(aobasic->app,aobasic->petsc);
 44:   PetscFree(aobasic);
 45:   return 0;
 46: }

 48: PetscErrorCode AOBasicGetIndices_Private(AO ao,PetscInt **app,PetscInt **petsc)
 49: {
 50:   AO_Basic *basic = (AO_Basic*)ao->data;

 52:   if (app)   *app   = basic->app;
 53:   if (petsc) *petsc = basic->petsc;
 54:   return 0;
 55: }

 57: PetscErrorCode AOPetscToApplication_Basic(AO ao,PetscInt n,PetscInt *ia)
 58: {
 59:   PetscInt i,N=ao->N;
 60:   AO_Basic *aobasic = (AO_Basic*)ao->data;

 62:   for (i=0; i<n; i++) {
 63:     if (ia[i] >= 0 && ia[i] < N) {
 64:       ia[i] = aobasic->app[ia[i]];
 65:     } else {
 66:       ia[i] = -1;
 67:     }
 68:   }
 69:   return 0;
 70: }

 72: PetscErrorCode AOApplicationToPetsc_Basic(AO ao,PetscInt n,PetscInt *ia)
 73: {
 74:   PetscInt i,N=ao->N;
 75:   AO_Basic *aobasic = (AO_Basic*)ao->data;

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

 87: PetscErrorCode AOPetscToApplicationPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array)
 88: {
 89:   AO_Basic       *aobasic = (AO_Basic*) ao->data;
 90:   PetscInt       *temp;
 91:   PetscInt       i, j;

 93:   PetscMalloc1(ao->N*block, &temp);
 94:   for (i = 0; i < ao->N; i++) {
 95:     for (j = 0; j < block; j++) temp[i*block+j] = array[aobasic->petsc[i]*block+j];
 96:   }
 97:   PetscArraycpy(array, temp, ao->N*block);
 98:   PetscFree(temp);
 99:   return 0;
100: }

102: PetscErrorCode AOApplicationToPetscPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array)
103: {
104:   AO_Basic       *aobasic = (AO_Basic*) ao->data;
105:   PetscInt       *temp;
106:   PetscInt       i, j;

108:   PetscMalloc1(ao->N*block, &temp);
109:   for (i = 0; i < ao->N; i++) {
110:     for (j = 0; j < block; j++) temp[i*block+j] = array[aobasic->app[i]*block+j];
111:   }
112:   PetscArraycpy(array, temp, ao->N*block);
113:   PetscFree(temp);
114:   return 0;
115: }

117: PetscErrorCode AOPetscToApplicationPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array)
118: {
119:   AO_Basic       *aobasic = (AO_Basic*) ao->data;
120:   PetscReal      *temp;
121:   PetscInt       i, j;

123:   PetscMalloc1(ao->N*block, &temp);
124:   for (i = 0; i < ao->N; i++) {
125:     for (j = 0; j < block; j++) temp[i*block+j] = array[aobasic->petsc[i]*block+j];
126:   }
127:   PetscArraycpy(array, temp, ao->N*block);
128:   PetscFree(temp);
129:   return 0;
130: }

132: PetscErrorCode AOApplicationToPetscPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array)
133: {
134:   AO_Basic       *aobasic = (AO_Basic*) ao->data;
135:   PetscReal      *temp;
136:   PetscInt       i, j;

138:   PetscMalloc1(ao->N*block, &temp);
139:   for (i = 0; i < ao->N; i++) {
140:     for (j = 0; j < block; j++) temp[i*block+j] = array[aobasic->app[i]*block+j];
141:   }
142:   PetscArraycpy(array, temp, ao->N*block);
143:   PetscFree(temp);
144:   return 0;
145: }

147: static struct _AOOps AOOps_Basic = {
148:   PetscDesignatedInitializer(view,AOView_Basic),
149:   PetscDesignatedInitializer(destroy,AODestroy_Basic),
150:   PetscDesignatedInitializer(petsctoapplication,AOPetscToApplication_Basic),
151:   PetscDesignatedInitializer(applicationtopetsc,AOApplicationToPetsc_Basic),
152:   PetscDesignatedInitializer(petsctoapplicationpermuteint,AOPetscToApplicationPermuteInt_Basic),
153:   PetscDesignatedInitializer(applicationtopetscpermuteint,AOApplicationToPetscPermuteInt_Basic),
154:   PetscDesignatedInitializer(petsctoapplicationpermutereal,AOPetscToApplicationPermuteReal_Basic),
155:   PetscDesignatedInitializer(applicationtopetscpermutereal,AOApplicationToPetscPermuteReal_Basic),
156: };

158: PETSC_EXTERN PetscErrorCode AOCreate_Basic(AO ao)
159: {
160:   AO_Basic       *aobasic;
161:   PetscMPIInt    size,rank,count,*lens,*disp;
162:   PetscInt       napp,*allpetsc,*allapp,ip,ia,N,i,*petsc=NULL,start;
163:   IS             isapp=ao->isapp,ispetsc=ao->ispetsc;
164:   MPI_Comm       comm;
165:   const PetscInt *myapp,*mypetsc=NULL;

167:   /* create special struct aobasic */
168:   PetscNewLog(ao,&aobasic);
169:   ao->data = (void*) aobasic;
170:   PetscMemcpy(ao->ops,&AOOps_Basic,sizeof(struct _AOOps));
171:   PetscObjectChangeTypeName((PetscObject)ao,AOBASIC);

173:   ISGetLocalSize(isapp,&napp);
174:   ISGetIndices(isapp,&myapp);

176:   PetscMPIIntCast(napp,&count);

178:   /* transmit all lengths to all processors */
179:   PetscObjectGetComm((PetscObject)isapp,&comm);
180:   MPI_Comm_size(comm, &size);
181:   MPI_Comm_rank(comm, &rank);
182:   PetscMalloc2(size, &lens,size,&disp);
183:   MPI_Allgather(&count, 1, MPI_INT, lens, 1, MPI_INT, comm);
184:   N    =  0;
185:   for (i = 0; i < size; i++) {
186:     PetscMPIIntCast(N,disp+i); /* = sum(lens[j]), j< i */
187:     N   += lens[i];
188:   }
189:   ao->N = N;
190:   ao->n = N;

192:   /* If mypetsc is 0 then use "natural" numbering */
193:   if (napp) {
194:     if (!ispetsc) {
195:       start = disp[rank];
196:       PetscMalloc1(napp+1, &petsc);
197:       for (i=0; i<napp; i++) petsc[i] = start + i;
198:     } else {
199:       ISGetIndices(ispetsc,&mypetsc);
200:       petsc = (PetscInt*)mypetsc;
201:     }
202:   }

204:   /* get all indices on all processors */
205:   PetscMalloc2(N,&allpetsc,N,&allapp);
206:   MPI_Allgatherv(petsc, count, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm);
207:   MPI_Allgatherv((void*)myapp, count, MPIU_INT, allapp, lens, disp, MPIU_INT, comm);
208:   PetscFree2(lens,disp);

210:   if (PetscDefined(USE_DEBUG)) {
211:     PetscInt *sorted;
212:     PetscMalloc1(N,&sorted);

214:     PetscArraycpy(sorted,allpetsc,N);
215:     PetscSortInt(N,sorted);
216:     for (i=0; i<N; i++) {
218:     }

220:     PetscArraycpy(sorted,allapp,N);
221:     PetscSortInt(N,sorted);
222:     for (i=0; i<N; i++) {
224:     }

226:     PetscFree(sorted);
227:   }

229:   /* generate a list of application and PETSc node numbers */
230:   PetscCalloc2(N, &aobasic->app,N,&aobasic->petsc);
231:   PetscLogObjectMemory((PetscObject)ao,2*N*sizeof(PetscInt));
232:   for (i = 0; i < N; i++) {
233:     ip = allpetsc[i];
234:     ia = allapp[i];
235:     /* check there are no duplicates */
237:     aobasic->app[ip] = ia + 1;
239:     aobasic->petsc[ia] = ip + 1;
240:   }
241:   if (napp && !mypetsc) {
242:     PetscFree(petsc);
243:   }
244:   PetscFree2(allpetsc,allapp);
245:   /* shift indices down by one */
246:   for (i = 0; i < N; i++) {
247:     aobasic->app[i]--;
248:     aobasic->petsc[i]--;
249:   }

251:   ISRestoreIndices(isapp,&myapp);
252:   if (napp) {
253:     if (ispetsc) {
254:       ISRestoreIndices(ispetsc,&mypetsc);
255:     } else {
256:       PetscFree(petsc);
257:     }
258:   }
259:   return 0;
260: }

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

265:    Collective

267:    Input Parameters:
268: +  comm - MPI communicator that is to share AO
269: .  napp - size of integer arrays
270: .  myapp - integer array that defines an ordering
271: -  mypetsc - integer array that defines another ordering (may be NULL to
272:              indicate the natural ordering, that is 0,1,2,3,...)

274:    Output Parameter:
275: .  aoout - the new application ordering

277:    Level: beginner

279:     Notes:
280:     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"
281:            in the indices. Use AOCreateMapping() or AOCreateMappingIS() if you wish to have "holes" in the indices.

283: .seealso: AOCreateBasicIS(), AODestroy(), AOPetscToApplication(), AOApplicationToPetsc()
284: @*/
285: PetscErrorCode  AOCreateBasic(MPI_Comm comm,PetscInt napp,const PetscInt myapp[],const PetscInt mypetsc[],AO *aoout)
286: {
287:   IS             isapp,ispetsc;
288:   const PetscInt *app=myapp,*petsc=mypetsc;

290:   ISCreateGeneral(comm,napp,app,PETSC_USE_POINTER,&isapp);
291:   if (mypetsc) {
292:     ISCreateGeneral(comm,napp,petsc,PETSC_USE_POINTER,&ispetsc);
293:   } else {
294:     ispetsc = NULL;
295:   }
296:   AOCreateBasicIS(isapp,ispetsc,aoout);
297:   ISDestroy(&isapp);
298:   if (mypetsc) {
299:     ISDestroy(&ispetsc);
300:   }
301:   return 0;
302: }

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

307:    Collective on IS

309:    Input Parameters:
310: +  isapp - index set that defines an ordering
311: -  ispetsc - index set that defines another ordering (may be NULL to use the
312:              natural ordering)

314:    Output Parameter:
315: .  aoout - the new application ordering

317:    Level: beginner

319:     Notes:
320:     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;
321:            that is there cannot be any "holes"

323: .seealso: AOCreateBasic(),  AODestroy()
324: @*/
325: PetscErrorCode AOCreateBasicIS(IS isapp,IS ispetsc,AO *aoout)
326: {
327:   MPI_Comm       comm;
328:   AO             ao;

330:   PetscObjectGetComm((PetscObject)isapp,&comm);
331:   AOCreate(comm,&ao);
332:   AOSetIS(ao,isapp,ispetsc);
333:   AOSetType(ao,AOBASIC);
334:   AOViewFromOptions(ao,NULL,"-ao_view");
335:   *aoout = ao;
336:   return 0;
337: }