Actual source code: aobasic.c

petsc-3.8.3 2017-12-09
Report Typos and Errors

  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: {
 20:   PetscMPIInt    rank;
 21:   PetscInt       i;
 22:   AO_Basic       *aobasic = (AO_Basic*)ao->data;
 23:   PetscBool      iascii;

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

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

 47:   PetscFree2(aobasic->app,aobasic->petsc);
 48:   PetscFree(aobasic);
 49:   return(0);
 50: }

 52: PetscErrorCode AOBasicGetIndices_Private(AO ao,PetscInt **app,PetscInt **petsc)
 53: {
 54:   AO_Basic *basic = (AO_Basic*)ao->data;

 57:   if (app)   *app   = basic->app;
 58:   if (petsc) *petsc = basic->petsc;
 59:   return(0);
 60: }

 62: PetscErrorCode AOPetscToApplication_Basic(AO ao,PetscInt n,PetscInt *ia)
 63: {
 64:   PetscInt i,N=ao->N;
 65:   AO_Basic *aobasic = (AO_Basic*)ao->data;

 68:   for (i=0; i<n; i++) {
 69:     if (ia[i] >= 0 && ia[i] < N) {
 70:       ia[i] = aobasic->app[ia[i]];
 71:     } else {
 72:       ia[i] = -1;
 73:     }
 74:   }
 75:   return(0);
 76: }

 78: PetscErrorCode AOApplicationToPetsc_Basic(AO ao,PetscInt n,PetscInt *ia)
 79: {
 80:   PetscInt i,N=ao->N;
 81:   AO_Basic *aobasic = (AO_Basic*)ao->data;

 84:   for (i=0; i<n; i++) {
 85:     if (ia[i] >= 0 && ia[i] < N) {
 86:       ia[i] = aobasic->petsc[ia[i]];
 87:     } else {
 88:       ia[i] = -1;
 89:     }
 90:   }
 91:   return(0);
 92: }

 94: PetscErrorCode AOPetscToApplicationPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array)
 95: {
 96:   AO_Basic       *aobasic = (AO_Basic*) ao->data;
 97:   PetscInt       *temp;
 98:   PetscInt       i, j;

102:   PetscMalloc1(ao->N*block, &temp);
103:   for (i = 0; i < ao->N; i++) {
104:     for (j = 0; j < block; j++) temp[i*block+j] = array[aobasic->petsc[i]*block+j];
105:   }
106:   PetscMemcpy(array, temp, ao->N*block * sizeof(PetscInt));
107:   PetscFree(temp);
108:   return(0);
109: }

111: PetscErrorCode AOApplicationToPetscPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array)
112: {
113:   AO_Basic       *aobasic = (AO_Basic*) ao->data;
114:   PetscInt       *temp;
115:   PetscInt       i, j;

119:   PetscMalloc1(ao->N*block, &temp);
120:   for (i = 0; i < ao->N; i++) {
121:     for (j = 0; j < block; j++) temp[i*block+j] = array[aobasic->app[i]*block+j];
122:   }
123:   PetscMemcpy(array, temp, ao->N*block * sizeof(PetscInt));
124:   PetscFree(temp);
125:   return(0);
126: }

128: PetscErrorCode AOPetscToApplicationPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array)
129: {
130:   AO_Basic       *aobasic = (AO_Basic*) ao->data;
131:   PetscReal      *temp;
132:   PetscInt       i, j;

136:   PetscMalloc1(ao->N*block, &temp);
137:   for (i = 0; i < ao->N; i++) {
138:     for (j = 0; j < block; j++) temp[i*block+j] = array[aobasic->petsc[i]*block+j];
139:   }
140:   PetscMemcpy(array, temp, ao->N*block * sizeof(PetscReal));
141:   PetscFree(temp);
142:   return(0);
143: }

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

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

162: static struct _AOOps AOOps_Basic = {
163:   AOView_Basic,
164:   AODestroy_Basic,
165:   AOPetscToApplication_Basic,
166:   AOApplicationToPetsc_Basic,
167:   AOPetscToApplicationPermuteInt_Basic,
168:   AOApplicationToPetscPermuteInt_Basic,
169:   AOPetscToApplicationPermuteReal_Basic,
170:   AOApplicationToPetscPermuteReal_Basic
171: };

173: PETSC_EXTERN PetscErrorCode AOCreate_Basic(AO ao)
174: {
175:   AO_Basic       *aobasic;
176:   PetscMPIInt    size,rank,count,*lens,*disp;
177:   PetscInt       napp,*allpetsc,*allapp,ip,ia,N,i,*petsc=NULL,start;
179:   IS             isapp=ao->isapp,ispetsc=ao->ispetsc;
180:   MPI_Comm       comm;
181:   const PetscInt *myapp,*mypetsc=NULL;

184:   /* create special struct aobasic */
185:   PetscNewLog(ao,&aobasic);
186:   ao->data = (void*) aobasic;
187:   PetscMemcpy(ao->ops,&AOOps_Basic,sizeof(struct _AOOps));
188:   PetscObjectChangeTypeName((PetscObject)ao,AOBASIC);

190:   ISGetLocalSize(isapp,&napp);
191:   ISGetIndices(isapp,&myapp);

193:   PetscMPIIntCast(napp,&count);

195:   /* transmit all lengths to all processors */
196:   PetscObjectGetComm((PetscObject)isapp,&comm);
197:   MPI_Comm_size(comm, &size);
198:   MPI_Comm_rank(comm, &rank);
199:   PetscMalloc2(size, &lens,size,&disp);
200:   MPI_Allgather(&count, 1, MPI_INT, lens, 1, MPI_INT, comm);
201:   N    =  0;
202:   for (i = 0; i < size; i++) {
203:     PetscMPIIntCast(N,disp+i); /* = sum(lens[j]), j< i */
204:     N   += lens[i];
205:   }
206:   ao->N = N;
207:   ao->n = N;

209:   /* If mypetsc is 0 then use "natural" numbering */
210:   if (napp) {
211:     if (!ispetsc) {
212:       start = disp[rank];
213:       PetscMalloc1(napp+1, &petsc);
214:       for (i=0; i<napp; i++) petsc[i] = start + i;
215:     } else {
216:       ISGetIndices(ispetsc,&mypetsc);
217:       petsc = (PetscInt*)mypetsc;
218:     }
219:   }

221:   /* get all indices on all processors */
222:   PetscMalloc2(N,&allpetsc,N,&allapp);
223:   MPI_Allgatherv(petsc, count, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm);
224:   MPI_Allgatherv((void*)myapp, count, MPIU_INT, allapp, lens, disp, MPIU_INT, comm);
225:   PetscFree2(lens,disp);

227: #if defined(PETSC_USE_DEBUG)
228:   {
229:     PetscInt *sorted;
230:     PetscMalloc1(N,&sorted);

232:     PetscMemcpy(sorted,allpetsc,N*sizeof(PetscInt));
233:     PetscSortInt(N,sorted);
234:     for (i=0; i<N; i++) {
235:       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]);
236:     }

238:     PetscMemcpy(sorted,allapp,N*sizeof(PetscInt));
239:     PetscSortInt(N,sorted);
240:     for (i=0; i<N; i++) {
241:       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]);
242:     }

244:     PetscFree(sorted);
245:   }
246: #endif

248:   /* generate a list of application and PETSc node numbers */
249:   PetscMalloc2(N, &aobasic->app,N,&aobasic->petsc);
250:   PetscLogObjectMemory((PetscObject)ao,2*N*sizeof(PetscInt));
251:   PetscMemzero(aobasic->app, N*sizeof(PetscInt));
252:   PetscMemzero(aobasic->petsc, N*sizeof(PetscInt));
253:   for (i = 0; i < N; i++) {
254:     ip = allpetsc[i];
255:     ia = allapp[i];
256:     /* check there are no duplicates */
257:     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);
258:     aobasic->app[ip] = ia + 1;
259:     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);
260:     aobasic->petsc[ia] = ip + 1;
261:   }
262:   if (napp && !mypetsc) {
263:     PetscFree(petsc);
264:   }
265:   PetscFree2(allpetsc,allapp);
266:   /* shift indices down by one */
267:   for (i = 0; i < N; i++) {
268:     aobasic->app[i]--;
269:     aobasic->petsc[i]--;
270:   }

272:   ISRestoreIndices(isapp,&myapp);
273:   if (napp) {
274:     if (ispetsc) {
275:       ISRestoreIndices(ispetsc,&mypetsc);
276:     } else {
277:       PetscFree(petsc);
278:     }
279:   }
280:   return(0);
281: }

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

286:    Collective on MPI_Comm

288:    Input Parameters:
289: +  comm - MPI communicator that is to share AO
290: .  napp - size of integer arrays
291: .  myapp - integer array that defines an ordering
292: -  mypetsc - integer array that defines another ordering (may be NULL to
293:              indicate the natural ordering, that is 0,1,2,3,...)

295:    Output Parameter:
296: .  aoout - the new application ordering

298:    Level: beginner

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

303: .keywords: AO, create

305: .seealso: AOCreateBasicIS(), AODestroy(), AOPetscToApplication(), AOApplicationToPetsc()
306: @*/
307: PetscErrorCode  AOCreateBasic(MPI_Comm comm,PetscInt napp,const PetscInt myapp[],const PetscInt mypetsc[],AO *aoout)
308: {
310:   IS             isapp,ispetsc;
311:   const PetscInt *app=myapp,*petsc=mypetsc;

314:   ISCreateGeneral(comm,napp,app,PETSC_USE_POINTER,&isapp);
315:   if (mypetsc) {
316:     ISCreateGeneral(comm,napp,petsc,PETSC_USE_POINTER,&ispetsc);
317:   } else {
318:     ispetsc = NULL;
319:   }
320:   AOCreateBasicIS(isapp,ispetsc,aoout);
321:   ISDestroy(&isapp);
322:   if (mypetsc) {
323:     ISDestroy(&ispetsc);
324:   }
325:   return(0);
326: }

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

331:    Collective on IS

333:    Input Parameters:
334: +  isapp - index set that defines an ordering
335: -  ispetsc - index set that defines another ordering (may be NULL to use the
336:              natural ordering)

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

341:    Level: beginner

343:     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;
344:            that is there cannot be any "holes"

346: .keywords: AO, create

348: .seealso: AOCreateBasic(),  AODestroy()
349: @*/
350: PetscErrorCode AOCreateBasicIS(IS isapp,IS ispetsc,AO *aoout)
351: {
353:   MPI_Comm       comm;
354:   AO             ao;

357:   PetscObjectGetComm((PetscObject)isapp,&comm);
358:   AOCreate(comm,&ao);
359:   AOSetIS(ao,isapp,ispetsc);
360:   AOSetType(ao,AOBASIC);
361:   AOViewFromOptions(ao,NULL,"-ao_view");
362:   *aoout = ao;
363:   return(0);
364: }