Actual source code: aobasic.c
petsc-3.10.5 2019-03-28
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: {
20: PetscMPIInt rank;
21: PetscInt i;
22: AO_Basic *aobasic = (AO_Basic*)ao->data;
23: PetscBool iascii;
26: MPI_Comm_rank(PetscObjectComm((PetscObject)ao),&rank);
27: if (!rank) {
28: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
29: if (iascii) {
30: PetscViewerASCIIPrintf(viewer,"Number of elements in ordering %D\n",ao->N);
31: PetscViewerASCIIPrintf(viewer, "PETSc->App App->PETSc\n");
32: for (i=0; i<ao->N; i++) {
33: PetscViewerASCIIPrintf(viewer,"%3D %3D %3D %3D\n",i,aobasic->app[i],i,aobasic->petsc[i]);
34: }
35: }
36: }
37: PetscViewerFlush(viewer);
38: return(0);
39: }
41: PetscErrorCode AODestroy_Basic(AO ao)
42: {
43: AO_Basic *aobasic = (AO_Basic*)ao->data;
47: PetscFree2(aobasic->app,aobasic->petsc);
48: PetscFree(aobasic);
49: return(0);
50: }
52: PetscErrorCode AOBasicGetIndices_Private(AO ao,PetscInt **app,PetscInt **petsc)
53: {
54: AO_Basic *basic = (AO_Basic*)ao->data;
57: if (app) *app = basic->app;
58: if (petsc) *petsc = basic->petsc;
59: return(0);
60: }
62: PetscErrorCode AOPetscToApplication_Basic(AO ao,PetscInt n,PetscInt *ia)
63: {
64: PetscInt i,N=ao->N;
65: AO_Basic *aobasic = (AO_Basic*)ao->data;
68: for (i=0; i<n; i++) {
69: if (ia[i] >= 0 && ia[i] < N) {
70: ia[i] = aobasic->app[ia[i]];
71: } else {
72: ia[i] = -1;
73: }
74: }
75: return(0);
76: }
78: PetscErrorCode AOApplicationToPetsc_Basic(AO ao,PetscInt n,PetscInt *ia)
79: {
80: PetscInt i,N=ao->N;
81: AO_Basic *aobasic = (AO_Basic*)ao->data;
84: for (i=0; i<n; i++) {
85: if (ia[i] >= 0 && ia[i] < N) {
86: ia[i] = aobasic->petsc[ia[i]];
87: } else {
88: ia[i] = -1;
89: }
90: }
91: return(0);
92: }
94: PetscErrorCode AOPetscToApplicationPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array)
95: {
96: AO_Basic *aobasic = (AO_Basic*) ao->data;
97: PetscInt *temp;
98: PetscInt i, j;
102: PetscMalloc1(ao->N*block, &temp);
103: for (i = 0; i < ao->N; i++) {
104: for (j = 0; j < block; j++) temp[i*block+j] = array[aobasic->petsc[i]*block+j];
105: }
106: PetscMemcpy(array, temp, ao->N*block * sizeof(PetscInt));
107: PetscFree(temp);
108: return(0);
109: }
111: PetscErrorCode AOApplicationToPetscPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array)
112: {
113: AO_Basic *aobasic = (AO_Basic*) ao->data;
114: PetscInt *temp;
115: PetscInt i, j;
119: PetscMalloc1(ao->N*block, &temp);
120: for (i = 0; i < ao->N; i++) {
121: for (j = 0; j < block; j++) temp[i*block+j] = array[aobasic->app[i]*block+j];
122: }
123: PetscMemcpy(array, temp, ao->N*block * sizeof(PetscInt));
124: PetscFree(temp);
125: return(0);
126: }
128: PetscErrorCode AOPetscToApplicationPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array)
129: {
130: AO_Basic *aobasic = (AO_Basic*) ao->data;
131: PetscReal *temp;
132: PetscInt i, j;
136: PetscMalloc1(ao->N*block, &temp);
137: for (i = 0; i < ao->N; i++) {
138: for (j = 0; j < block; j++) temp[i*block+j] = array[aobasic->petsc[i]*block+j];
139: }
140: PetscMemcpy(array, temp, ao->N*block * sizeof(PetscReal));
141: PetscFree(temp);
142: return(0);
143: }
145: PetscErrorCode AOApplicationToPetscPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array)
146: {
147: AO_Basic *aobasic = (AO_Basic*) ao->data;
148: PetscReal *temp;
149: PetscInt i, j;
153: PetscMalloc1(ao->N*block, &temp);
154: for (i = 0; i < ao->N; i++) {
155: for (j = 0; j < block; j++) temp[i*block+j] = array[aobasic->app[i]*block+j];
156: }
157: PetscMemcpy(array, temp, ao->N*block * sizeof(PetscReal));
158: PetscFree(temp);
159: return(0);
160: }
162: static struct _AOOps AOOps_Basic = {
163: AOView_Basic,
164: AODestroy_Basic,
165: AOPetscToApplication_Basic,
166: AOApplicationToPetsc_Basic,
167: AOPetscToApplicationPermuteInt_Basic,
168: AOApplicationToPetscPermuteInt_Basic,
169: AOPetscToApplicationPermuteReal_Basic,
170: AOApplicationToPetscPermuteReal_Basic
171: };
173: PETSC_EXTERN PetscErrorCode AOCreate_Basic(AO ao)
174: {
175: AO_Basic *aobasic;
176: PetscMPIInt size,rank,count,*lens,*disp;
177: PetscInt napp,*allpetsc,*allapp,ip,ia,N,i,*petsc=NULL,start;
179: IS isapp=ao->isapp,ispetsc=ao->ispetsc;
180: MPI_Comm comm;
181: const PetscInt *myapp,*mypetsc=NULL;
184: /* create special struct aobasic */
185: PetscNewLog(ao,&aobasic);
186: ao->data = (void*) aobasic;
187: PetscMemcpy(ao->ops,&AOOps_Basic,sizeof(struct _AOOps));
188: PetscObjectChangeTypeName((PetscObject)ao,AOBASIC);
190: ISGetLocalSize(isapp,&napp);
191: ISGetIndices(isapp,&myapp);
193: PetscMPIIntCast(napp,&count);
195: /* transmit all lengths to all processors */
196: PetscObjectGetComm((PetscObject)isapp,&comm);
197: MPI_Comm_size(comm, &size);
198: MPI_Comm_rank(comm, &rank);
199: PetscMalloc2(size, &lens,size,&disp);
200: MPI_Allgather(&count, 1, MPI_INT, lens, 1, MPI_INT, comm);
201: N = 0;
202: for (i = 0; i < size; i++) {
203: PetscMPIIntCast(N,disp+i); /* = sum(lens[j]), j< i */
204: N += lens[i];
205: }
206: ao->N = N;
207: ao->n = N;
209: /* If mypetsc is 0 then use "natural" numbering */
210: if (napp) {
211: if (!ispetsc) {
212: start = disp[rank];
213: PetscMalloc1(napp+1, &petsc);
214: for (i=0; i<napp; i++) petsc[i] = start + i;
215: } else {
216: ISGetIndices(ispetsc,&mypetsc);
217: petsc = (PetscInt*)mypetsc;
218: }
219: }
221: /* get all indices on all processors */
222: PetscMalloc2(N,&allpetsc,N,&allapp);
223: MPI_Allgatherv(petsc, count, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm);
224: MPI_Allgatherv((void*)myapp, count, MPIU_INT, allapp, lens, disp, MPIU_INT, comm);
225: PetscFree2(lens,disp);
227: #if defined(PETSC_USE_DEBUG)
228: {
229: PetscInt *sorted;
230: PetscMalloc1(N,&sorted);
232: PetscMemcpy(sorted,allpetsc,N*sizeof(PetscInt));
233: PetscSortInt(N,sorted);
234: for (i=0; i<N; i++) {
235: 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]);
236: }
238: PetscMemcpy(sorted,allapp,N*sizeof(PetscInt));
239: PetscSortInt(N,sorted);
240: for (i=0; i<N; i++) {
241: 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]);
242: }
244: PetscFree(sorted);
245: }
246: #endif
248: /* generate a list of application and PETSc node numbers */
249: PetscMalloc2(N, &aobasic->app,N,&aobasic->petsc);
250: PetscLogObjectMemory((PetscObject)ao,2*N*sizeof(PetscInt));
251: PetscMemzero(aobasic->app, N*sizeof(PetscInt));
252: PetscMemzero(aobasic->petsc, N*sizeof(PetscInt));
253: for (i = 0; i < N; i++) {
254: ip = allpetsc[i];
255: ia = allapp[i];
256: /* check there are no duplicates */
257: 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);
258: aobasic->app[ip] = ia + 1;
259: 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);
260: aobasic->petsc[ia] = ip + 1;
261: }
262: if (napp && !mypetsc) {
263: PetscFree(petsc);
264: }
265: PetscFree2(allpetsc,allapp);
266: /* shift indices down by one */
267: for (i = 0; i < N; i++) {
268: aobasic->app[i]--;
269: aobasic->petsc[i]--;
270: }
272: ISRestoreIndices(isapp,&myapp);
273: if (napp) {
274: if (ispetsc) {
275: ISRestoreIndices(ispetsc,&mypetsc);
276: } else {
277: PetscFree(petsc);
278: }
279: }
280: return(0);
281: }
283: /*@C
284: AOCreateBasic - Creates a basic application ordering using two integer arrays.
286: Collective on MPI_Comm
288: Input Parameters:
289: + comm - MPI communicator that is to share AO
290: . napp - size of integer arrays
291: . myapp - integer array that defines an ordering
292: - mypetsc - integer array that defines another ordering (may be NULL to
293: indicate the natural ordering, that is 0,1,2,3,...)
295: Output Parameter:
296: . aoout - the new application ordering
298: Level: beginner
300: Notes:
301: 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"
302: in the indices. Use AOCreateMapping() or AOCreateMappingIS() if you wish to have "holes" in the indices.
304: .keywords: AO, create
306: .seealso: AOCreateBasicIS(), AODestroy(), AOPetscToApplication(), AOApplicationToPetsc()
307: @*/
308: PetscErrorCode AOCreateBasic(MPI_Comm comm,PetscInt napp,const PetscInt myapp[],const PetscInt mypetsc[],AO *aoout)
309: {
311: IS isapp,ispetsc;
312: const PetscInt *app=myapp,*petsc=mypetsc;
315: ISCreateGeneral(comm,napp,app,PETSC_USE_POINTER,&isapp);
316: if (mypetsc) {
317: ISCreateGeneral(comm,napp,petsc,PETSC_USE_POINTER,&ispetsc);
318: } else {
319: ispetsc = NULL;
320: }
321: AOCreateBasicIS(isapp,ispetsc,aoout);
322: ISDestroy(&isapp);
323: if (mypetsc) {
324: ISDestroy(&ispetsc);
325: }
326: return(0);
327: }
329: /*@C
330: AOCreateBasicIS - Creates a basic application ordering using two index sets.
332: Collective on IS
334: Input Parameters:
335: + isapp - index set that defines an ordering
336: - ispetsc - index set that defines another ordering (may be NULL to use the
337: natural ordering)
339: Output Parameter:
340: . aoout - the new application ordering
342: Level: beginner
344: Notes:
345: 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;
346: that is there cannot be any "holes"
348: .keywords: AO, create
350: .seealso: AOCreateBasic(), AODestroy()
351: @*/
352: PetscErrorCode AOCreateBasicIS(IS isapp,IS ispetsc,AO *aoout)
353: {
355: MPI_Comm comm;
356: AO ao;
359: PetscObjectGetComm((PetscObject)isapp,&comm);
360: AOCreate(comm,&ao);
361: AOSetIS(ao,isapp,ispetsc);
362: AOSetType(ao,AOBASIC);
363: AOViewFromOptions(ao,NULL,"-ao_view");
364: *aoout = ao;
365: return(0);
366: }