Actual source code: aobasic.c

  1: #define PETSCDM_DLL

  3: /*
  4:     The most basic AO application ordering routines. These store the 
  5:   entire orderings on each processor.
  6: */

 8:  #include src/dm/ao/aoimpl.h
 9:  #include petscsys.h

 11: typedef struct {
 12:   PetscInt N;
 13:   PetscInt *app,*petsc;  /* app[i] is the partner for the ith PETSc slot */
 14:                          /* petsc[j] is the partner for the jth app slot */
 15: } AO_Basic;

 17: /*
 18:        All processors have the same data so processor 1 prints it
 19: */
 22: PetscErrorCode AOView_Basic(AO ao,PetscViewer viewer)
 23: {
 25:   PetscMPIInt    rank;
 26:   PetscInt       i;
 27:   AO_Basic       *aodebug = (AO_Basic*)ao->data;
 28:   PetscTruth     iascii;

 31:   MPI_Comm_rank(ao->comm,&rank);
 32:   if (!rank){
 33:     PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
 34:     if (iascii) {
 35:       PetscViewerASCIIPrintf(viewer,"Number of elements in ordering %D\n",aodebug->N);
 36:       PetscViewerASCIIPrintf(viewer,  "PETSc->App  App->PETSc\n");
 37:       for (i=0; i<aodebug->N; i++) {
 38:         PetscViewerASCIIPrintf(viewer,"%3D  %3D    %3D  %3D\n",i,aodebug->app[i],i,aodebug->petsc[i]);
 39:       }
 40:     } else {
 41:       SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for AOData basic",((PetscObject)viewer)->type_name);
 42:     }
 43:   }
 44:   PetscViewerFlush(viewer);
 45:   return(0);
 46: }

 50: PetscErrorCode AODestroy_Basic(AO ao)
 51: {
 52:   AO_Basic       *aodebug = (AO_Basic*)ao->data;

 56:   PetscFree(aodebug->app);
 57:   PetscFree(ao->data);
 58:   return(0);
 59: }

 63: PetscErrorCode AOBasicGetIndices_Private(AO ao,PetscInt **app,PetscInt **petsc)
 64: {
 65:   AO_Basic *basic = (AO_Basic*)ao->data;

 68:   if (app)   *app   = basic->app;
 69:   if (petsc) *petsc = basic->petsc;
 70:   return(0);
 71: }

 75: PetscErrorCode AOPetscToApplication_Basic(AO ao,PetscInt n,PetscInt *ia)
 76: {
 77:   PetscInt i;
 78:   AO_Basic *aodebug = (AO_Basic*)ao->data;

 81:   for (i=0; i<n; i++) {
 82:     if (ia[i] >= 0) {ia[i] = aodebug->app[ia[i]];}
 83:   }
 84:   return(0);
 85: }

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

 95:   for (i=0; i<n; i++) {
 96:     if (ia[i] >= 0) {ia[i] = aodebug->petsc[ia[i]];}
 97:   }
 98:   return(0);
 99: }

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

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

122: PetscErrorCode AOApplicationToPetscPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array)
123: {
124:   AO_Basic       *aodebug = (AO_Basic *) ao->data;
125:   PetscInt       *temp;
126:   PetscInt       i, j;

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

141: PetscErrorCode AOPetscToApplicationPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array)
142: {
143:   AO_Basic       *aodebug = (AO_Basic *) ao->data;
144:   PetscReal      *temp;
145:   PetscInt       i, j;

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

160: PetscErrorCode AOApplicationToPetscPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array)
161: {
162:   AO_Basic       *aodebug = (AO_Basic *) ao->data;
163:   PetscReal      *temp;
164:   PetscInt       i, j;

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

177: static struct _AOOps AOops = {AOView_Basic,
178:                               AODestroy_Basic,
179:                               AOPetscToApplication_Basic,
180:                               AOApplicationToPetsc_Basic,
181:                               AOPetscToApplicationPermuteInt_Basic,
182:                               AOApplicationToPetscPermuteInt_Basic,
183:                               AOPetscToApplicationPermuteReal_Basic,
184:                               AOApplicationToPetscPermuteReal_Basic};

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

191:    Collective on MPI_Comm

193:    Input Parameters:
194: +  comm - MPI communicator that is to share AO
195: .  napp - size of integer arrays
196: .  myapp - integer array that defines an ordering
197: -  mypetsc - integer array that defines another ordering (may be PETSC_NULL to 
198:              indicate the natural ordering)

200:    Output Parameter:
201: .  aoout - the new application ordering

203:    Options Database Key:
204: .   -ao_view - call AOView() at the conclusion of AOCreateBasic()

206:    Level: beginner

208: .keywords: AO, create

210: .seealso: AOCreateBasicIS(), AODestroy()
211: @*/
212: PetscErrorCode  AOCreateBasic(MPI_Comm comm,PetscInt napp,const PetscInt myapp[],const PetscInt mypetsc[],AO *aoout)
213: {
214:   AO_Basic       *aobasic;
215:   AO             ao;
216:   PetscMPIInt    *lens,size,rank,nnapp,*disp;
217:   PetscInt       *allpetsc,*allapp,ip,ia,N,i,*petsc,start;
218:   PetscTruth     opt;

223:   *aoout = 0;
224: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
225:   DMInitializePackage(PETSC_NULL);
226: #endif

228:   PetscHeaderCreate(ao, _p_AO, struct _AOOps, AO_COOKIE, AO_BASIC, "AO", comm, AODestroy, AOView);
229:   PetscNew(AO_Basic, &aobasic);
230:   PetscLogObjectMemory(ao, sizeof(struct _p_AO) + sizeof(AO_Basic));

232:   PetscMemcpy(ao->ops, &AOops, sizeof(AOops));
233:   ao->data = (void*) aobasic;

235:   /* transmit all lengths to all processors */
236:   MPI_Comm_size(comm, &size);
237:   MPI_Comm_rank(comm, &rank);
238:   PetscMalloc(2*size * sizeof(PetscMPIInt), &lens);
239:   disp  = lens + size;
240:   nnapp = napp;
241:   MPI_Allgather(&nnapp, 1, MPI_INT, lens, 1, MPI_INT, comm);
242:   N    =  0;
243:   for(i = 0; i < size; i++) {
244:     disp[i] = N;
245:     N += lens[i];
246:   }
247:   aobasic->N = N;

249:   /*
250:      If mypetsc is 0 then use "natural" numbering 
251:   */
252:   if (napp && !mypetsc) {
253:     start = disp[rank];
254:     PetscMalloc((napp+1) * sizeof(PetscInt), &petsc);
255:     for (i=0; i<napp; i++) {
256:       petsc[i] = start + i;
257:     }
258:   } else {
259:     petsc = (PetscInt*)mypetsc;
260:   }

262:   /* get all indices on all processors */
263:   PetscMalloc(2*N * sizeof(PetscInt), &allpetsc);
264:   allapp = allpetsc + N;
265:   MPI_Allgatherv(petsc, napp, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm);
266:   MPI_Allgatherv((void*)myapp, napp, MPIU_INT, allapp, lens, disp, MPIU_INT, comm);
267:   PetscFree(lens);

269: #if defined(PETSC_USE_DEBUG)
270:   {
271:     PetscInt *sorted;
272:     PetscMalloc((N+1)*sizeof(PetscInt),&sorted);

274:     PetscMemcpy(sorted,allapp,N*sizeof(PetscInt));
275:     PetscSortInt(N,sorted);
276:     for (i=0; i<N; i++) {
277:       if (sorted[i] != i) SETERRQ2(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]);
278:     }

280:     PetscMemcpy(sorted,allapp,N*sizeof(PetscInt));
281:     PetscSortInt(N,sorted);
282:     for (i=0; i<N; i++) {
283:       if (sorted[i] != i) SETERRQ2(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]);
284:     }

286:     PetscFree(sorted);
287:   }
288: #endif

290:   /* generate a list of application and PETSc node numbers */
291:   PetscMalloc(2*N * sizeof(PetscInt), &aobasic->app);
292:   PetscLogObjectMemory(ao,2*N*sizeof(PetscInt));
293:   aobasic->petsc = aobasic->app + N;
294:   PetscMemzero(aobasic->app, 2*N*sizeof(PetscInt));
295:   for(i = 0; i < N; i++) {
296:     ip = allpetsc[i];
297:     ia = allapp[i];
298:     /* check there are no duplicates */
299:     if (aobasic->app[ip]) SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"Duplicate in PETSc ordering at position %d. Already mapped to %d, not %d.", i, aobasic->app[ip]-1, ia);
300:     aobasic->app[ip] = ia + 1;
301:     if (aobasic->petsc[ia]) SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"Duplicate in Application ordering at position %d. Already mapped to %d, not %d.", i, aobasic->petsc[ia]-1, ip);
302:     aobasic->petsc[ia] = ip + 1;
303:   }
304:   if (!mypetsc) {
305:     PetscFree(petsc);
306:   }
307:   PetscFree(allpetsc);
308:   /* shift indices down by one */
309:   for(i = 0; i < N; i++) {
310:     aobasic->app[i]--;
311:     aobasic->petsc[i]--;
312:   }

314:   PetscOptionsHasName(PETSC_NULL, "-ao_view", &opt);
315:   if (opt) {
316:     AOView(ao, PETSC_VIEWER_STDOUT_SELF);
317:   }

319:   *aoout = ao;
320:   return(0);
321: }

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

328:    Collective on IS

330:    Input Parameters:
331: +  isapp - index set that defines an ordering
332: -  ispetsc - index set that defines another ordering (may be PETSC_NULL to use the
333:              natural ordering)

335:    Output Parameter:
336: .  aoout - the new application ordering

338:    Options Database Key:
339: -   -ao_view - call AOView() at the conclusion of AOCreateBasicIS()

341:    Level: beginner

343: .keywords: AO, create

345: .seealso: AOCreateBasic(),  AODestroy()
346: @*/
347: PetscErrorCode  AOCreateBasicIS(IS isapp,IS ispetsc,AO *aoout)
348: {
350:   PetscInt       *mypetsc = 0,*myapp,napp,npetsc;
351:   MPI_Comm       comm;

354:   PetscObjectGetComm((PetscObject)isapp,&comm);
355:   ISGetLocalSize(isapp,&napp);
356:   if (ispetsc) {
357:     ISGetLocalSize(ispetsc,&npetsc);
358:     if (napp != npetsc) SETERRQ(PETSC_ERR_ARG_SIZ,"Local IS lengths must match");
359:     ISGetIndices(ispetsc,&mypetsc);
360:   }
361:   ISGetIndices(isapp,&myapp);

363:   AOCreateBasic(comm,napp,myapp,mypetsc,aoout);

365:   ISRestoreIndices(isapp,&myapp);
366:   if (ispetsc) {
367:     ISRestoreIndices(ispetsc,&mypetsc);
368:   }
369:   return(0);
370: }