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