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:   PetscFunctionBegin;
 25:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ao), &rank));
 26:   if (rank == 0) {
 27:     PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
 28:     if (iascii) {
 29:       PetscCall(PetscViewerASCIIPrintf(viewer, "Number of elements in ordering %" PetscInt_FMT "\n", ao->N));
 30:       PetscCall(PetscViewerASCIIPrintf(viewer, "PETSc->App  App->PETSc\n"));
 31:       for (i = 0; i < ao->N; i++) PetscCall(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:   PetscCall(PetscViewerFlush(viewer));
 35:   PetscFunctionReturn(PETSC_SUCCESS);
 36: }

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

 42:   PetscFunctionBegin;
 43:   PetscCall(PetscFree2(aobasic->app, aobasic->petsc));
 44:   PetscCall(PetscFree(aobasic));
 45:   PetscFunctionReturn(PETSC_SUCCESS);
 46: }

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

 52:   PetscFunctionBegin;
 53:   if (app) *app = basic->app;
 54:   if (petsc) *petsc = basic->petsc;
 55:   PetscFunctionReturn(PETSC_SUCCESS);
 56: }

 58: PetscErrorCode AOPetscToApplication_Basic(AO ao, PetscInt n, PetscInt *ia)
 59: {
 60:   PetscInt  i, N = ao->N;
 61:   AO_Basic *aobasic = (AO_Basic *)ao->data;

 63:   PetscFunctionBegin;
 64:   for (i = 0; i < n; i++) {
 65:     if (ia[i] >= 0 && ia[i] < N) {
 66:       ia[i] = aobasic->app[ia[i]];
 67:     } else {
 68:       ia[i] = -1;
 69:     }
 70:   }
 71:   PetscFunctionReturn(PETSC_SUCCESS);
 72: }

 74: PetscErrorCode AOApplicationToPetsc_Basic(AO ao, PetscInt n, PetscInt *ia)
 75: {
 76:   PetscInt  i, N = ao->N;
 77:   AO_Basic *aobasic = (AO_Basic *)ao->data;

 79:   PetscFunctionBegin;
 80:   for (i = 0; i < n; i++) {
 81:     if (ia[i] >= 0 && ia[i] < N) {
 82:       ia[i] = aobasic->petsc[ia[i]];
 83:     } else {
 84:       ia[i] = -1;
 85:     }
 86:   }
 87:   PetscFunctionReturn(PETSC_SUCCESS);
 88: }

 90: PetscErrorCode AOPetscToApplicationPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array)
 91: {
 92:   AO_Basic *aobasic = (AO_Basic *)ao->data;
 93:   PetscInt *temp;
 94:   PetscInt  i, j;

 96:   PetscFunctionBegin;
 97:   PetscCall(PetscMalloc1(ao->N * block, &temp));
 98:   for (i = 0; i < ao->N; i++) {
 99:     for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->petsc[i] * block + j];
100:   }
101:   PetscCall(PetscArraycpy(array, temp, ao->N * block));
102:   PetscCall(PetscFree(temp));
103:   PetscFunctionReturn(PETSC_SUCCESS);
104: }

106: PetscErrorCode AOApplicationToPetscPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array)
107: {
108:   AO_Basic *aobasic = (AO_Basic *)ao->data;
109:   PetscInt *temp;
110:   PetscInt  i, j;

112:   PetscFunctionBegin;
113:   PetscCall(PetscMalloc1(ao->N * block, &temp));
114:   for (i = 0; i < ao->N; i++) {
115:     for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->app[i] * block + j];
116:   }
117:   PetscCall(PetscArraycpy(array, temp, ao->N * block));
118:   PetscCall(PetscFree(temp));
119:   PetscFunctionReturn(PETSC_SUCCESS);
120: }

122: PetscErrorCode AOPetscToApplicationPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array)
123: {
124:   AO_Basic  *aobasic = (AO_Basic *)ao->data;
125:   PetscReal *temp;
126:   PetscInt   i, j;

128:   PetscFunctionBegin;
129:   PetscCall(PetscMalloc1(ao->N * block, &temp));
130:   for (i = 0; i < ao->N; i++) {
131:     for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->petsc[i] * block + j];
132:   }
133:   PetscCall(PetscArraycpy(array, temp, ao->N * block));
134:   PetscCall(PetscFree(temp));
135:   PetscFunctionReturn(PETSC_SUCCESS);
136: }

138: PetscErrorCode AOApplicationToPetscPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array)
139: {
140:   AO_Basic  *aobasic = (AO_Basic *)ao->data;
141:   PetscReal *temp;
142:   PetscInt   i, j;

144:   PetscFunctionBegin;
145:   PetscCall(PetscMalloc1(ao->N * block, &temp));
146:   for (i = 0; i < ao->N; i++) {
147:     for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->app[i] * block + j];
148:   }
149:   PetscCall(PetscArraycpy(array, temp, ao->N * block));
150:   PetscCall(PetscFree(temp));
151:   PetscFunctionReturn(PETSC_SUCCESS);
152: }

154: static struct _AOOps AOOps_Basic = {
155:   PetscDesignatedInitializer(view, AOView_Basic),
156:   PetscDesignatedInitializer(destroy, AODestroy_Basic),
157:   PetscDesignatedInitializer(petsctoapplication, AOPetscToApplication_Basic),
158:   PetscDesignatedInitializer(applicationtopetsc, AOApplicationToPetsc_Basic),
159:   PetscDesignatedInitializer(petsctoapplicationpermuteint, AOPetscToApplicationPermuteInt_Basic),
160:   PetscDesignatedInitializer(applicationtopetscpermuteint, AOApplicationToPetscPermuteInt_Basic),
161:   PetscDesignatedInitializer(petsctoapplicationpermutereal, AOPetscToApplicationPermuteReal_Basic),
162:   PetscDesignatedInitializer(applicationtopetscpermutereal, AOApplicationToPetscPermuteReal_Basic),
163: };

165: PETSC_EXTERN PetscErrorCode AOCreate_Basic(AO ao)
166: {
167:   AO_Basic       *aobasic;
168:   PetscMPIInt     size, rank, count, *lens, *disp;
169:   PetscInt        napp, *allpetsc, *allapp, ip, ia, N, i, *petsc = NULL, start;
170:   IS              isapp = ao->isapp, ispetsc = ao->ispetsc;
171:   MPI_Comm        comm;
172:   const PetscInt *myapp, *mypetsc = NULL;

174:   PetscFunctionBegin;
175:   /* create special struct aobasic */
176:   PetscCall(PetscNew(&aobasic));
177:   ao->data = (void *)aobasic;
178:   PetscCall(PetscMemcpy(ao->ops, &AOOps_Basic, sizeof(struct _AOOps)));
179:   PetscCall(PetscObjectChangeTypeName((PetscObject)ao, AOBASIC));

181:   PetscCall(ISGetLocalSize(isapp, &napp));
182:   PetscCall(ISGetIndices(isapp, &myapp));

184:   PetscCall(PetscMPIIntCast(napp, &count));

186:   /* transmit all lengths to all processors */
187:   PetscCall(PetscObjectGetComm((PetscObject)isapp, &comm));
188:   PetscCallMPI(MPI_Comm_size(comm, &size));
189:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
190:   PetscCall(PetscMalloc2(size, &lens, size, &disp));
191:   PetscCallMPI(MPI_Allgather(&count, 1, MPI_INT, lens, 1, MPI_INT, comm));
192:   N = 0;
193:   for (i = 0; i < size; i++) {
194:     PetscCall(PetscMPIIntCast(N, disp + i)); /* = sum(lens[j]), j< i */
195:     N += lens[i];
196:   }
197:   ao->N = N;
198:   ao->n = N;

200:   /* If mypetsc is 0 then use "natural" numbering */
201:   if (napp) {
202:     if (!ispetsc) {
203:       start = disp[rank];
204:       PetscCall(PetscMalloc1(napp + 1, &petsc));
205:       for (i = 0; i < napp; i++) petsc[i] = start + i;
206:     } else {
207:       PetscCall(ISGetIndices(ispetsc, &mypetsc));
208:       petsc = (PetscInt *)mypetsc;
209:     }
210:   }

212:   /* get all indices on all processors */
213:   PetscCall(PetscMalloc2(N, &allpetsc, N, &allapp));
214:   PetscCallMPI(MPI_Allgatherv(petsc, count, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm));
215:   PetscCallMPI(MPI_Allgatherv((void *)myapp, count, MPIU_INT, allapp, lens, disp, MPIU_INT, comm));
216:   PetscCall(PetscFree2(lens, disp));

218:   if (PetscDefined(USE_DEBUG)) {
219:     PetscInt *sorted;
220:     PetscCall(PetscMalloc1(N, &sorted));

222:     PetscCall(PetscArraycpy(sorted, allpetsc, N));
223:     PetscCall(PetscSortInt(N, sorted));
224:     for (i = 0; i < N; i++) PetscCheck(sorted[i] == i, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PETSc ordering requires a permutation of numbers 0 to N-1\n it is missing %" PetscInt_FMT " has %" PetscInt_FMT, i, sorted[i]);

226:     PetscCall(PetscArraycpy(sorted, allapp, N));
227:     PetscCall(PetscSortInt(N, sorted));
228:     for (i = 0; i < N; i++) PetscCheck(sorted[i] == i, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Application ordering requires a permutation of numbers 0 to N-1\n it is missing %" PetscInt_FMT " has %" PetscInt_FMT, i, sorted[i]);

230:     PetscCall(PetscFree(sorted));
231:   }

233:   /* generate a list of application and PETSc node numbers */
234:   PetscCall(PetscCalloc2(N, &aobasic->app, N, &aobasic->petsc));
235:   for (i = 0; i < N; i++) {
236:     ip = allpetsc[i];
237:     ia = allapp[i];
238:     /* check there are no duplicates */
239:     PetscCheck(!aobasic->app[ip], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Duplicate in PETSc ordering at position %" PetscInt_FMT ". Already mapped to %" PetscInt_FMT ", not %" PetscInt_FMT ".", i, aobasic->app[ip] - 1, ia);
240:     aobasic->app[ip] = ia + 1;
241:     PetscCheck(!aobasic->petsc[ia], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Duplicate in Application ordering at position %" PetscInt_FMT ". Already mapped to %" PetscInt_FMT ", not %" PetscInt_FMT ".", i, aobasic->petsc[ia] - 1, ip);
242:     aobasic->petsc[ia] = ip + 1;
243:   }
244:   if (napp && !mypetsc) PetscCall(PetscFree(petsc));
245:   PetscCall(PetscFree2(allpetsc, allapp));
246:   /* shift indices down by one */
247:   for (i = 0; i < N; i++) {
248:     aobasic->app[i]--;
249:     aobasic->petsc[i]--;
250:   }

252:   PetscCall(ISRestoreIndices(isapp, &myapp));
253:   if (napp) {
254:     if (ispetsc) {
255:       PetscCall(ISRestoreIndices(ispetsc, &mypetsc));
256:     } else {
257:       PetscCall(PetscFree(petsc));
258:     }
259:   }
260:   PetscFunctionReturn(PETSC_SUCCESS);
261: }

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

266:    Collective

268:    Input Parameters:
269: +  comm - MPI communicator that is to share `AO`
270: .  napp - size of integer arrays
271: .  myapp - integer array that defines an ordering
272: -  mypetsc - integer array that defines another ordering (may be NULL to
273:              indicate the natural ordering, that is 0,1,2,3,...)

275:    Output Parameter:
276: .  aoout - the new application ordering

278:    Level: beginner

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

284: .seealso: [](sec_ao), [](sec_scatter), `AO`, `AOCreateBasicIS()`, `AODestroy()`, `AOPetscToApplication()`, `AOApplicationToPetsc()`
285: @*/
286: PetscErrorCode AOCreateBasic(MPI_Comm comm, PetscInt napp, const PetscInt myapp[], const PetscInt mypetsc[], AO *aoout)
287: {
288:   IS              isapp, ispetsc;
289:   const PetscInt *app = myapp, *petsc = mypetsc;

291:   PetscFunctionBegin;
292:   PetscCall(ISCreateGeneral(comm, napp, app, PETSC_USE_POINTER, &isapp));
293:   if (mypetsc) {
294:     PetscCall(ISCreateGeneral(comm, napp, petsc, PETSC_USE_POINTER, &ispetsc));
295:   } else {
296:     ispetsc = NULL;
297:   }
298:   PetscCall(AOCreateBasicIS(isapp, ispetsc, aoout));
299:   PetscCall(ISDestroy(&isapp));
300:   if (mypetsc) PetscCall(ISDestroy(&ispetsc));
301:   PetscFunctionReturn(PETSC_SUCCESS);
302: }

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

307:    Collective

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:     Note:
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: [](sec_ao), [](sec_scatter), `IS`, `AO`, `AOCreateBasic()`, `AODestroy()`
324: @*/
325: PetscErrorCode AOCreateBasicIS(IS isapp, IS ispetsc, AO *aoout)
326: {
327:   MPI_Comm comm;
328:   AO       ao;

330:   PetscFunctionBegin;
331:   PetscCall(PetscObjectGetComm((PetscObject)isapp, &comm));
332:   PetscCall(AOCreate(comm, &ao));
333:   PetscCall(AOSetIS(ao, isapp, ispetsc));
334:   PetscCall(AOSetType(ao, AOBASIC));
335:   PetscCall(AOViewFromOptions(ao, NULL, "-ao_view"));
336:   *aoout = ao;
337:   PetscFunctionReturn(PETSC_SUCCESS);
338: }