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++) PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT "  %3" PetscInt_FMT "    %3" PetscInt_FMT "  %3" PetscInt_FMT "\n", i, aobasic->app[i], i, aobasic->petsc[i]);
 31:     }
 32:   }
 33:   PetscViewerFlush(viewer);
 34:   return 0;
 35: }

 37: PetscErrorCode AODestroy_Basic(AO ao)
 38: {
 39:   AO_Basic *aobasic = (AO_Basic *)ao->data;

 41:   PetscFree2(aobasic->app, aobasic->petsc);
 42:   PetscFree(aobasic);
 43:   return 0;
 44: }

 46: PetscErrorCode AOBasicGetIndices_Private(AO ao, PetscInt **app, PetscInt **petsc)
 47: {
 48:   AO_Basic *basic = (AO_Basic *)ao->data;

 50:   if (app) *app = basic->app;
 51:   if (petsc) *petsc = basic->petsc;
 52:   return 0;
 53: }

 55: PetscErrorCode AOPetscToApplication_Basic(AO ao, PetscInt n, PetscInt *ia)
 56: {
 57:   PetscInt  i, N = ao->N;
 58:   AO_Basic *aobasic = (AO_Basic *)ao->data;

 60:   for (i = 0; i < n; i++) {
 61:     if (ia[i] >= 0 && ia[i] < N) {
 62:       ia[i] = aobasic->app[ia[i]];
 63:     } else {
 64:       ia[i] = -1;
 65:     }
 66:   }
 67:   return 0;
 68: }

 70: PetscErrorCode AOApplicationToPetsc_Basic(AO ao, PetscInt n, PetscInt *ia)
 71: {
 72:   PetscInt  i, N = ao->N;
 73:   AO_Basic *aobasic = (AO_Basic *)ao->data;

 75:   for (i = 0; i < n; i++) {
 76:     if (ia[i] >= 0 && ia[i] < N) {
 77:       ia[i] = aobasic->petsc[ia[i]];
 78:     } else {
 79:       ia[i] = -1;
 80:     }
 81:   }
 82:   return 0;
 83: }

 85: PetscErrorCode AOPetscToApplicationPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array)
 86: {
 87:   AO_Basic *aobasic = (AO_Basic *)ao->data;
 88:   PetscInt *temp;
 89:   PetscInt  i, j;

 91:   PetscMalloc1(ao->N * block, &temp);
 92:   for (i = 0; i < ao->N; i++) {
 93:     for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->petsc[i] * block + j];
 94:   }
 95:   PetscArraycpy(array, temp, ao->N * block);
 96:   PetscFree(temp);
 97:   return 0;
 98: }

100: PetscErrorCode AOApplicationToPetscPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array)
101: {
102:   AO_Basic *aobasic = (AO_Basic *)ao->data;
103:   PetscInt *temp;
104:   PetscInt  i, j;

106:   PetscMalloc1(ao->N * block, &temp);
107:   for (i = 0; i < ao->N; i++) {
108:     for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->app[i] * block + j];
109:   }
110:   PetscArraycpy(array, temp, ao->N * block);
111:   PetscFree(temp);
112:   return 0;
113: }

115: PetscErrorCode AOPetscToApplicationPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array)
116: {
117:   AO_Basic  *aobasic = (AO_Basic *)ao->data;
118:   PetscReal *temp;
119:   PetscInt   i, j;

121:   PetscMalloc1(ao->N * block, &temp);
122:   for (i = 0; i < ao->N; i++) {
123:     for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->petsc[i] * block + j];
124:   }
125:   PetscArraycpy(array, temp, ao->N * block);
126:   PetscFree(temp);
127:   return 0;
128: }

130: PetscErrorCode AOApplicationToPetscPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array)
131: {
132:   AO_Basic  *aobasic = (AO_Basic *)ao->data;
133:   PetscReal *temp;
134:   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->app[i] * block + j];
139:   }
140:   PetscArraycpy(array, temp, ao->N * block);
141:   PetscFree(temp);
142:   return 0;
143: }

145: static struct _AOOps AOOps_Basic = {
146:   PetscDesignatedInitializer(view, AOView_Basic),
147:   PetscDesignatedInitializer(destroy, AODestroy_Basic),
148:   PetscDesignatedInitializer(petsctoapplication, AOPetscToApplication_Basic),
149:   PetscDesignatedInitializer(applicationtopetsc, AOApplicationToPetsc_Basic),
150:   PetscDesignatedInitializer(petsctoapplicationpermuteint, AOPetscToApplicationPermuteInt_Basic),
151:   PetscDesignatedInitializer(applicationtopetscpermuteint, AOApplicationToPetscPermuteInt_Basic),
152:   PetscDesignatedInitializer(petsctoapplicationpermutereal, AOPetscToApplicationPermuteReal_Basic),
153:   PetscDesignatedInitializer(applicationtopetscpermutereal, AOApplicationToPetscPermuteReal_Basic),
154: };

156: PETSC_EXTERN PetscErrorCode AOCreate_Basic(AO ao)
157: {
158:   AO_Basic       *aobasic;
159:   PetscMPIInt     size, rank, count, *lens, *disp;
160:   PetscInt        napp, *allpetsc, *allapp, ip, ia, N, i, *petsc = NULL, start;
161:   IS              isapp = ao->isapp, ispetsc = ao->ispetsc;
162:   MPI_Comm        comm;
163:   const PetscInt *myapp, *mypetsc = NULL;

165:   /* create special struct aobasic */
166:   PetscNew(&aobasic);
167:   ao->data = (void *)aobasic;
168:   PetscMemcpy(ao->ops, &AOOps_Basic, sizeof(struct _AOOps));
169:   PetscObjectChangeTypeName((PetscObject)ao, AOBASIC);

171:   ISGetLocalSize(isapp, &napp);
172:   ISGetIndices(isapp, &myapp);

174:   PetscMPIIntCast(napp, &count);

176:   /* transmit all lengths to all processors */
177:   PetscObjectGetComm((PetscObject)isapp, &comm);
178:   MPI_Comm_size(comm, &size);
179:   MPI_Comm_rank(comm, &rank);
180:   PetscMalloc2(size, &lens, size, &disp);
181:   MPI_Allgather(&count, 1, MPI_INT, lens, 1, MPI_INT, comm);
182:   N = 0;
183:   for (i = 0; i < size; i++) {
184:     PetscMPIIntCast(N, disp + i); /* = sum(lens[j]), j< i */
185:     N += lens[i];
186:   }
187:   ao->N = N;
188:   ao->n = N;

190:   /* If mypetsc is 0 then use "natural" numbering */
191:   if (napp) {
192:     if (!ispetsc) {
193:       start = disp[rank];
194:       PetscMalloc1(napp + 1, &petsc);
195:       for (i = 0; i < napp; i++) petsc[i] = start + i;
196:     } else {
197:       ISGetIndices(ispetsc, &mypetsc);
198:       petsc = (PetscInt *)mypetsc;
199:     }
200:   }

202:   /* get all indices on all processors */
203:   PetscMalloc2(N, &allpetsc, N, &allapp);
204:   MPI_Allgatherv(petsc, count, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm);
205:   MPI_Allgatherv((void *)myapp, count, MPIU_INT, allapp, lens, disp, MPIU_INT, comm);
206:   PetscFree2(lens, disp);

208:   if (PetscDefined(USE_DEBUG)) {
209:     PetscInt *sorted;
210:     PetscMalloc1(N, &sorted);

212:     PetscArraycpy(sorted, allpetsc, N);
213:     PetscSortInt(N, sorted);

216:     PetscArraycpy(sorted, allapp, N);
217:     PetscSortInt(N, sorted);

220:     PetscFree(sorted);
221:   }

223:   /* generate a list of application and PETSc node numbers */
224:   PetscCalloc2(N, &aobasic->app, N, &aobasic->petsc);
225:   for (i = 0; i < N; i++) {
226:     ip = allpetsc[i];
227:     ia = allapp[i];
228:     /* check there are no duplicates */
230:     aobasic->app[ip] = ia + 1;
232:     aobasic->petsc[ia] = ip + 1;
233:   }
234:   if (napp && !mypetsc) PetscFree(petsc);
235:   PetscFree2(allpetsc, allapp);
236:   /* shift indices down by one */
237:   for (i = 0; i < N; i++) {
238:     aobasic->app[i]--;
239:     aobasic->petsc[i]--;
240:   }

242:   ISRestoreIndices(isapp, &myapp);
243:   if (napp) {
244:     if (ispetsc) {
245:       ISRestoreIndices(ispetsc, &mypetsc);
246:     } else {
247:       PetscFree(petsc);
248:     }
249:   }
250:   return 0;
251: }

253: /*@C
254:    AOCreateBasic - Creates a basic application ordering using two integer arrays.

256:    Collective

258:    Input Parameters:
259: +  comm - MPI communicator that is to share `AO`
260: .  napp - size of integer arrays
261: .  myapp - integer array that defines an ordering
262: -  mypetsc - integer array that defines another ordering (may be NULL to
263:              indicate the natural ordering, that is 0,1,2,3,...)

265:    Output Parameter:
266: .  aoout - the new application ordering

268:    Level: beginner

270:    Note:
271:    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"
272:    in the indices. Use `AOCreateMapping()` or `AOCreateMappingIS()` if you wish to have "holes" in the indices.

274: .seealso: [](sec_ao), [](sec_scatter), `AO`, `AOCreateBasicIS()`, `AODestroy()`, `AOPetscToApplication()`, `AOApplicationToPetsc()`
275: @*/
276: PetscErrorCode AOCreateBasic(MPI_Comm comm, PetscInt napp, const PetscInt myapp[], const PetscInt mypetsc[], AO *aoout)
277: {
278:   IS              isapp, ispetsc;
279:   const PetscInt *app = myapp, *petsc = mypetsc;

281:   ISCreateGeneral(comm, napp, app, PETSC_USE_POINTER, &isapp);
282:   if (mypetsc) {
283:     ISCreateGeneral(comm, napp, petsc, PETSC_USE_POINTER, &ispetsc);
284:   } else {
285:     ispetsc = NULL;
286:   }
287:   AOCreateBasicIS(isapp, ispetsc, aoout);
288:   ISDestroy(&isapp);
289:   if (mypetsc) ISDestroy(&ispetsc);
290:   return 0;
291: }

293: /*@C
294:    AOCreateBasicIS - Creates a basic application ordering using two `IS` index sets.

296:    Collective

298:    Input Parameters:
299: +  isapp - index set that defines an ordering
300: -  ispetsc - index set that defines another ordering (may be NULL to use the
301:              natural ordering)

303:    Output Parameter:
304: .  aoout - the new application ordering

306:    Level: beginner

308:     Note:
309:     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;
310:     that is there cannot be any "holes"

312: .seealso: [](sec_ao), [](sec_scatter), `IS`, `AO`, `AOCreateBasic()`, `AODestroy()`
313: @*/
314: PetscErrorCode AOCreateBasicIS(IS isapp, IS ispetsc, AO *aoout)
315: {
316:   MPI_Comm comm;
317:   AO       ao;

319:   PetscObjectGetComm((PetscObject)isapp, &comm);
320:   AOCreate(comm, &ao);
321:   AOSetIS(ao, isapp, ispetsc);
322:   AOSetType(ao, AOBASIC);
323:   AOViewFromOptions(ao, NULL, "-ao_view");
324:   *aoout = ao;
325:   return 0;
326: }