Actual source code: aobasic.c

  1: /*$Id: aobasic.c,v 1.60 2001/03/23 23:24:52 balay Exp $*/

  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:   int N;
 13:   int *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: */
 20: #undef __FUNCT__  
 22: int AOView_Basic(AO ao,PetscViewer viewer)
 23: {
 24:   int        rank,ierr,i;
 25:   AO_Basic   *aodebug = (AO_Basic*)ao->data;
 26:   PetscTruth isascii;

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

 46: #undef __FUNCT__  
 48: int AODestroy_Basic(AO ao)
 49: {
 50:   AO_Basic *aodebug = (AO_Basic*)ao->data;
 51:   int      ierr;

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

 59: #undef __FUNCT__  
 61: int AOBasicGetIndices_Private(AO ao,int **app,int **petsc)
 62: {
 63:   AO_Basic *basic = (AO_Basic*)ao->data;

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

 71: #undef __FUNCT__  
 73: int AOPetscToApplication_Basic(AO ao,int n,int *ia)
 74: {
 75:   int      i;
 76:   AO_Basic *aodebug = (AO_Basic*)ao->data;

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

 85: #undef __FUNCT__  
 87: int AOApplicationToPetsc_Basic(AO ao,int n,int *ia)
 88: {
 89:   int      i;
 90:   AO_Basic *aodebug = (AO_Basic*)ao->data;

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

 99: #undef __FUNCT__  
101: int AOPetscToApplicationPermuteInt_Basic(AO ao, int block, int *array)
102: {
103:   AO_Basic *aodebug = (AO_Basic *) ao->data;
104:   int      *temp;
105:   int       i, j;
106:   int       ierr;

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

118: #undef __FUNCT__  
120: int AOApplicationToPetscPermuteInt_Basic(AO ao, int block, int *array)
121: {
122:   AO_Basic *aodebug = (AO_Basic *) ao->data;
123:   int      *temp;
124:   int       i, j;
125:   int       ierr;

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

137: #undef __FUNCT__  
139: int AOPetscToApplicationPermuteReal_Basic(AO ao, int block, double *array)
140: {
141:   AO_Basic *aodebug = (AO_Basic *) ao->data;
142:   double   *temp;
143:   int       i, j;
144:   int       ierr;

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

156: #undef __FUNCT__  
158: int AOApplicationToPetscPermuteReal_Basic(AO ao, int block, double *array)
159: {
160:   AO_Basic *aodebug = (AO_Basic *) ao->data;
161:   double   *temp;
162:   int       i, j;
163:   int       ierr;

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

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

184: EXTERN_C_BEGIN
185: #undef __FUNCT__  
187: int AOSerialize_Basic(MPI_Comm comm, AO *ao, PetscViewer viewer, PetscTruth store)
188: {
189:   AO        a;
190:   AO_Basic *d;
191:   int       fd;
192:   int       ierr;

195:   PetscViewerBinaryGetDescriptor(viewer, &fd);
196:   if (store) {
197:     a = *ao;
198:     d = (AO_Basic *) a->data;
199:     PetscBinaryWrite(fd, &a->N,     1,    PETSC_INT, 0);
200:     PetscBinaryWrite(fd, &a->n,     1,    PETSC_INT, 0);
201:     PetscBinaryWrite(fd, &d->N,     1,    PETSC_INT, 0);
202:     PetscBinaryWrite(fd,  d->app,   d->N, PETSC_INT, 0);
203:     PetscBinaryWrite(fd,  d->petsc, d->N, PETSC_INT, 0);
204:   } else {
205:     PetscHeaderCreate(a, _p_AO, struct _AOOps, AO_COOKIE, AO_BASIC, "AO", comm, AODestroy, AOView);
206:     PetscLogObjectCreate(a);
207:     PetscLogObjectMemory(a, sizeof(AO_Basic) + sizeof(struct _p_AO));
208:     PetscNew(AO_Basic, &d);
209:     a->data = (void *) d;
210:     PetscMemcpy(a->ops, &AOops, sizeof(AOops));

212:     PetscBinaryRead(fd, &a->N,     1,    PETSC_INT);
213:     PetscBinaryRead(fd, &a->n,     1,    PETSC_INT);
214:     PetscBinaryRead(fd, &d->N,     1,    PETSC_INT);
215:     PetscMalloc(d->N*2 * sizeof(int), &d->app);
216:     d->petsc = d->app + d->N;
217:     PetscLogObjectMemory(a, d->N*2 * sizeof(int));
218:     PetscBinaryRead(fd,  d->app,   d->N, PETSC_INT);
219:     PetscBinaryRead(fd,  d->petsc, d->N, PETSC_INT);
220:     *ao = a;
221:   }

223:   return(0);
224: }
225: EXTERN_C_END

227: #undef __FUNCT__  
229: /*@C
230:    AOCreateBasic - Creates a basic application ordering using two integer arrays.

232:    Collective on MPI_Comm

234:    Input Parameters:
235: +  comm - MPI communicator that is to share AO
236: .  napp - size of integer arrays
237: .  myapp - integer array that defines an ordering
238: -  mypetsc - integer array that defines another ordering (may be PETSC_NULL to 
239:              indicate the natural ordering)

241:    Output Parameter:
242: .  aoout - the new application ordering

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

247:    Level: beginner

249: .keywords: AO, create

251: .seealso: AOCreateBasicIS(), AODestroy()
252: @*/
253: int AOCreateBasic(MPI_Comm comm,int napp,int *myapp,int *mypetsc,AO *aoout)
254: {
255:   AO_Basic   *aobasic;
256:   AO         ao;
257:   int        *lens,size,rank,N,i,*petsc,start;
258:   int        *allpetsc,*allapp,*disp,ip,ia;
259:   PetscTruth opt;
260:   int        ierr;

264:   *aoout = 0;
265: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
266:   DMInitializePackage(PETSC_NULL);
267: #endif

269:   PetscHeaderCreate(ao, _p_AO, struct _AOOps, AO_COOKIE, AO_BASIC, "AO", comm, AODestroy, AOView);
270:   PetscLogObjectCreate(ao);
271:   PetscNew(AO_Basic, &aobasic);
272:   PetscLogObjectMemory(ao, sizeof(struct _p_AO) + sizeof(AO_Basic));

274:   PetscMemcpy(ao->ops, &AOops, sizeof(AOops));
275:   PetscStrallocpy(AO_SER_BASIC_BINARY, &ao->serialize_name);
276:   ao->data = (void *) aobasic;

278:   /* transmit all lengths to all processors */
279:   MPI_Comm_size(comm, &size);
280:   MPI_Comm_rank(comm, &rank);
281:   PetscMalloc(2*size * sizeof(int), &lens);
282:   disp = lens + size;
283:   MPI_Allgather(&napp, 1, MPI_INT, lens, 1, MPI_INT, comm);
284:   N    =  0;
285:   for(i = 0; i < size; i++) {
286:     disp[i] = N;
287:     N += lens[i];
288:   }
289:   aobasic->N = N;

291:   /*
292:      If mypetsc is 0 then use "natural" numbering 
293:   */
294:   if (!mypetsc) {
295:     start = disp[rank];
296:     ierr  = PetscMalloc((napp+1) * sizeof(int), &petsc);
297:     for (i=0; i<napp; i++) {
298:       petsc[i] = start + i;
299:     }
300:   } else {
301:     petsc = mypetsc;
302:   }

304:   /* get all indices on all processors */
305:   ierr   = PetscMalloc(2*N * sizeof(int), &allpetsc);
306:   allapp = allpetsc + N;
307:   ierr   = MPI_Allgatherv(petsc, napp, MPI_INT, allpetsc, lens, disp, MPI_INT, comm);
308:   ierr   = MPI_Allgatherv(myapp, napp, MPI_INT, allapp, lens, disp, MPI_INT, comm);
309:   ierr   = PetscFree(lens);

311:   /* generate a list of application and PETSc node numbers */
312:   PetscMalloc(2*N * sizeof(int), &aobasic->app);
313:   PetscLogObjectMemory(ao,2*N*sizeof(int));
314:   aobasic->petsc = aobasic->app + N;
315:   PetscMemzero(aobasic->app, 2*N*sizeof(int));
316:   for(i = 0; i < N; i++) {
317:     ip = allpetsc[i];
318:     ia = allapp[i];
319:     /* check there are no duplicates */
320:     if (aobasic->app[ip]) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Duplicate in PETSc ordering");
321:     aobasic->app[ip] = ia + 1;
322:     if (aobasic->petsc[ia]) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Duplicate in Application ordering");
323:     aobasic->petsc[ia] = ip + 1;
324:   }
325:   if (!mypetsc) {
326:     PetscFree(petsc);
327:   }
328:   PetscFree(allpetsc);
329:   /* shift indices down by one */
330:   for(i = 0; i < N; i++) {
331:     aobasic->app[i]--;
332:     aobasic->petsc[i]--;
333:   }

335:   PetscOptionsHasName(PETSC_NULL, "-ao_view", &opt);
336:   if (opt) {
337:     AOView(ao, PETSC_VIEWER_STDOUT_SELF);
338:   }

340:   *aoout = ao;
341:   return(0);
342: }

344: #undef __FUNCT__  
346: /*@C
347:    AOCreateBasicIS - Creates a basic application ordering using two index sets.

349:    Collective on IS

351:    Input Parameters:
352: +  isapp - index set that defines an ordering
353: -  ispetsc - index set that defines another ordering (may be PETSC_NULL to use the
354:              natural ordering)

356:    Output Parameter:
357: .  aoout - the new application ordering

359:    Options Database Key:
360: -   -ao_view - call AOView() at the conclusion of AOCreateBasicIS()

362:    Level: beginner

364: .keywords: AO, create

366: .seealso: AOCreateBasic(),  AODestroy()
367: @*/
368: int AOCreateBasicIS(IS isapp,IS ispetsc,AO *aoout)
369: {
370:   int       *mypetsc = 0,*myapp,ierr,napp,npetsc;
371:   MPI_Comm  comm;

374:   PetscObjectGetComm((PetscObject)isapp,&comm);
375:   ISGetLocalSize(isapp,&napp);
376:   if (ispetsc) {
377:     ISGetLocalSize(ispetsc,&npetsc);
378:     if (napp != npetsc) SETERRQ(PETSC_ERR_ARG_SIZ,"Local IS lengths must match");
379:     ISGetIndices(ispetsc,&mypetsc);
380:   }
381:   ISGetIndices(isapp,&myapp);

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

385:   ISRestoreIndices(isapp,&myapp);
386:   if (ispetsc) {
387:     ISRestoreIndices(ispetsc,&mypetsc);
388:   }
389:   return(0);
390: }