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: }