Actual source code: aobasic.c
petsc-3.3-p7 2013-05-11
2: /*
3: The most basic AO application ordering routines. These store the
4: entire orderings on each processor.
5: */
7: #include <../src/dm/ao/aoimpl.h> /*I "petscao.h" I*/
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: */
19: PetscErrorCode AOView_Basic(AO ao,PetscViewer viewer)
20: {
22: PetscMPIInt rank;
23: PetscInt i;
24: AO_Basic *aobasic = (AO_Basic*)ao->data;
25: PetscBool iascii;
28: MPI_Comm_rank(((PetscObject)ao)->comm,&rank);
29: if (!rank){
30: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
31: if (iascii) {
32: PetscViewerASCIIPrintf(viewer,"Number of elements in ordering %D\n",ao->N);
33: PetscViewerASCIIPrintf(viewer, "PETSc->App App->PETSc\n");
34: for (i=0; i<ao->N; i++) {
35: PetscViewerASCIIPrintf(viewer,"%3D %3D %3D %3D\n",i,aobasic->app[i],i,aobasic->petsc[i]);
36: }
37: } else SETERRQ1(((PetscObject)viewer)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for AO basic",((PetscObject)viewer)->type_name);
38: }
39: PetscViewerFlush(viewer);
40: return(0);
41: }
45: PetscErrorCode AODestroy_Basic(AO ao)
46: {
47: AO_Basic *aobasic = (AO_Basic*)ao->data;
51: PetscFree2(aobasic->app,aobasic->petsc);
52: PetscFree(aobasic);
53: return(0);
54: }
58: PetscErrorCode AOBasicGetIndices_Private(AO ao,PetscInt **app,PetscInt **petsc)
59: {
60: AO_Basic *basic = (AO_Basic*)ao->data;
63: if (app) *app = basic->app;
64: if (petsc) *petsc = basic->petsc;
65: return(0);
66: }
70: PetscErrorCode AOPetscToApplication_Basic(AO ao,PetscInt n,PetscInt *ia)
71: {
72: PetscInt i,N=ao->N;
73: AO_Basic *aobasic = (AO_Basic*)ao->data;
76: for (i=0; i<n; i++) {
77: if (ia[i] >= 0 && ia[i] < N ) {
78: ia[i] = aobasic->app[ia[i]];
79: } else {
80: ia[i] = -1;
81: }
82: }
83: return(0);
84: }
88: PetscErrorCode AOApplicationToPetsc_Basic(AO ao,PetscInt n,PetscInt *ia)
89: {
90: PetscInt i,N=ao->N;
91: AO_Basic *aobasic = (AO_Basic*)ao->data;
94: for (i=0; i<n; i++) {
95: if (ia[i] >= 0 && ia[i] < N) {
96: ia[i] = aobasic->petsc[ia[i]];
97: } else {
98: ia[i] = -1;
99: }
100: }
101: return(0);
102: }
106: PetscErrorCode AOPetscToApplicationPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array)
107: {
108: AO_Basic *aobasic = (AO_Basic *) ao->data;
109: PetscInt *temp;
110: PetscInt i, j;
114: PetscMalloc(ao->N*block * sizeof(PetscInt), &temp);
115: for(i = 0; i < ao->N; i++) {
116: for(j = 0; j < block; j++) temp[i*block+j] = array[aobasic->petsc[i]*block+j];
117: }
118: PetscMemcpy(array, temp, ao->N*block * sizeof(PetscInt));
119: PetscFree(temp);
120: return(0);
121: }
125: PetscErrorCode AOApplicationToPetscPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array)
126: {
127: AO_Basic *aobasic = (AO_Basic *) ao->data;
128: PetscInt *temp;
129: PetscInt i, j;
133: PetscMalloc(ao->N*block * sizeof(PetscInt), &temp);
134: for(i = 0; i < ao->N; i++) {
135: for(j = 0; j < block; j++) temp[i*block+j] = array[aobasic->app[i]*block+j];
136: }
137: PetscMemcpy(array, temp, ao->N*block * sizeof(PetscInt));
138: PetscFree(temp);
139: return(0);
140: }
144: PetscErrorCode AOPetscToApplicationPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array)
145: {
146: AO_Basic *aobasic = (AO_Basic *) ao->data;
147: PetscReal *temp;
148: PetscInt i, j;
152: PetscMalloc(ao->N*block * sizeof(PetscReal), &temp);
153: for(i = 0; i < ao->N; i++) {
154: for(j = 0; j < block; j++) temp[i*block+j] = array[aobasic->petsc[i]*block+j];
155: }
156: PetscMemcpy(array, temp, ao->N*block * sizeof(PetscReal));
157: PetscFree(temp);
158: return(0);
159: }
163: PetscErrorCode AOApplicationToPetscPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array)
164: {
165: AO_Basic *aobasic = (AO_Basic *) ao->data;
166: PetscReal *temp;
167: PetscInt i, j;
171: PetscMalloc(ao->N*block * sizeof(PetscReal), &temp);
172: for(i = 0; i < ao->N; i++) {
173: for(j = 0; j < block; j++) temp[i*block+j] = array[aobasic->app[i]*block+j];
174: }
175: PetscMemcpy(array, temp, ao->N*block * sizeof(PetscReal));
176: PetscFree(temp);
177: return(0);
178: }
180: static struct _AOOps AOOps_Basic = {
181: AOView_Basic,
182: AODestroy_Basic,
183: AOPetscToApplication_Basic,
184: AOApplicationToPetsc_Basic,
185: AOPetscToApplicationPermuteInt_Basic,
186: AOApplicationToPetscPermuteInt_Basic,
187: AOPetscToApplicationPermuteReal_Basic,
188: AOApplicationToPetscPermuteReal_Basic};
190: EXTERN_C_BEGIN
193: PetscErrorCode AOCreate_Basic(AO ao)
194: {
195: AO_Basic *aobasic;
196: PetscMPIInt size,rank,count,*lens,*disp;
197: PetscInt napp,*allpetsc,*allapp,ip,ia,N,i,*petsc=PETSC_NULL,start;
199: IS isapp=ao->isapp,ispetsc=ao->ispetsc;
200: MPI_Comm comm;
201: const PetscInt *myapp,*mypetsc=PETSC_NULL;
204: /* create special struct aobasic */
205: PetscNewLog(ao, AO_Basic, &aobasic);
206: ao->data = (void*) aobasic;
207: PetscMemcpy(ao->ops,&AOOps_Basic,sizeof(struct _AOOps));
208: PetscObjectChangeTypeName((PetscObject)ao,AOBASIC);
210: ISGetLocalSize(isapp,&napp);
211: ISGetIndices(isapp,&myapp);
213: count = PetscMPIIntCast(napp);
215: /* transmit all lengths to all processors */
216: PetscObjectGetComm((PetscObject)isapp,&comm);
217: MPI_Comm_size(comm, &size);
218: MPI_Comm_rank(comm, &rank);
219: PetscMalloc2(size,PetscMPIInt, &lens,size,PetscMPIInt,&disp);
220: MPI_Allgather(&count, 1, MPI_INT, lens, 1, MPI_INT, comm);
221: N = 0;
222: for(i = 0; i < size; i++) {
223: disp[i] = PetscMPIIntCast(N); /* = sum(lens[j]), j< i */
224: N += lens[i];
225: }
226: ao->N = N;
227: ao->n = N;
229: /* If mypetsc is 0 then use "natural" numbering */
230: if (napp){
231: if (!ispetsc) {
232: start = disp[rank];
233: PetscMalloc((napp+1) * sizeof(PetscInt), &petsc);
234: for (i=0; i<napp; i++) petsc[i] = start + i;
235: } else {
236: ISGetIndices(ispetsc,&mypetsc);
237: petsc = (PetscInt*)mypetsc;
238: }
239: }
241: /* get all indices on all processors */
242: PetscMalloc2(N,PetscInt,&allpetsc,N,PetscInt,&allapp);
243: MPI_Allgatherv(petsc, count, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm);
244: MPI_Allgatherv((void*)myapp, count, MPIU_INT, allapp, lens, disp, MPIU_INT, comm);
245: PetscFree2(lens,disp);
247: #if defined(PETSC_USE_DEBUG)
248: {
249: PetscInt *sorted;
250: PetscMalloc(N*sizeof(PetscInt),&sorted);
252: PetscMemcpy(sorted,allpetsc,N*sizeof(PetscInt));
253: PetscSortInt(N,sorted);
254: for (i=0; i<N; i++) {
255: 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]);
256: }
258: PetscMemcpy(sorted,allapp,N*sizeof(PetscInt));
259: PetscSortInt(N,sorted);
260: for (i=0; i<N; i++) {
261: 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]);
262: }
264: PetscFree(sorted);
265: }
266: #endif
268: /* generate a list of application and PETSc node numbers */
269: PetscMalloc2(N,PetscInt, &aobasic->app,N,PetscInt,&aobasic->petsc);
270: PetscLogObjectMemory(ao,2*N*sizeof(PetscInt));
271: PetscMemzero(aobasic->app, N*sizeof(PetscInt));
272: PetscMemzero(aobasic->petsc, N*sizeof(PetscInt));
273: for(i = 0; i < N; i++) {
274: ip = allpetsc[i];
275: ia = allapp[i];
276: /* check there are no duplicates */
277: 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);
278: aobasic->app[ip] = ia + 1;
279: 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);
280: aobasic->petsc[ia] = ip + 1;
281: }
282: if (napp && !mypetsc) {
283: PetscFree(petsc);
284: }
285: PetscFree2(allpetsc,allapp);
286: /* shift indices down by one */
287: for(i = 0; i < N; i++) {
288: aobasic->app[i]--;
289: aobasic->petsc[i]--;
290: }
292: ISRestoreIndices(isapp,&myapp);
293: if (napp){
294: if (ispetsc){
295: ISRestoreIndices(ispetsc,&mypetsc);
296: } else {
297: PetscFree(petsc);
298: }
299: }
300: return(0);
301: }
302: EXTERN_C_END
306: /*@C
307: AOCreateBasic - Creates a basic application ordering using two integer arrays.
309: Collective on MPI_Comm
311: Input Parameters:
312: + comm - MPI communicator that is to share AO
313: . napp - size of integer arrays
314: . myapp - integer array that defines an ordering
315: - mypetsc - integer array that defines another ordering (may be PETSC_NULL to
316: indicate the natural ordering, that is 0,1,2,3,...)
318: Output Parameter:
319: . aoout - the new application ordering
321: Level: beginner
323: 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"
324: in the indices. Use AOCreateMapping() or AOCreateMappingIS() if you wish to have "holes" in the indices.
326: .keywords: AO, create
328: .seealso: AOCreateBasicIS(), AODestroy(), AOPetscToApplication(), AOApplicationToPetsc()
329: @*/
330: PetscErrorCode AOCreateBasic(MPI_Comm comm,PetscInt napp,const PetscInt myapp[],const PetscInt mypetsc[],AO *aoout)
331: {
333: IS isapp,ispetsc;
334: const PetscInt *app=myapp,*petsc=mypetsc;
337: ISCreateGeneral(comm,napp,app,PETSC_USE_POINTER,&isapp);
338: if (mypetsc){
339: ISCreateGeneral(comm,napp,petsc,PETSC_USE_POINTER,&ispetsc);
340: } else {
341: ispetsc = PETSC_NULL;
342: }
343: AOCreateBasicIS(isapp,ispetsc,aoout);
344: ISDestroy(&isapp);
345: if (mypetsc){
346: ISDestroy(&ispetsc);
347: }
348: return(0);
349: }
353: /*@C
354: AOCreateBasicIS - Creates a basic application ordering using two index sets.
356: Collective on IS
358: Input Parameters:
359: + isapp - index set that defines an ordering
360: - ispetsc - index set that defines another ordering (may be PETSC_NULL to use the
361: natural ordering)
363: Output Parameter:
364: . aoout - the new application ordering
366: Level: beginner
368: 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;
369: that is there cannot be any "holes"
371: .keywords: AO, create
373: .seealso: AOCreateBasic(), AODestroy()
374: @*/
375: PetscErrorCode AOCreateBasicIS(IS isapp,IS ispetsc,AO *aoout)
376: {
378: MPI_Comm comm;
379: AO ao;
382: PetscObjectGetComm((PetscObject)isapp,&comm);
383: AOCreate(comm,&ao);
384: AOSetIS(ao,isapp,ispetsc);
385: AOSetType(ao,AOBASIC);
386: *aoout = ao;
387: return(0);
388: }