Actual source code: aomemscalable.c


  2: /*
  3:     The memory scalable AO application ordering routines. These store the
  4:   orderings on each processor for that processor's range of values
  5: */

  7: #include <../src/vec/is/ao/aoimpl.h>

  9: typedef struct {
 10:   PetscInt   *app_loc;   /* app_loc[i] is the partner for the ith local PETSc slot */
 11:   PetscInt   *petsc_loc; /* petsc_loc[j] is the partner for the jth local app slot */
 12:   PetscLayout map;       /* determines the local sizes of ao */
 13: } AO_MemoryScalable;

 15: /*
 16:        All processors ship the data to process 0 to be printed; note that this is not scalable because
 17:        process 0 allocates space for all the orderings entry across all the processes
 18: */
 19: PetscErrorCode AOView_MemoryScalable(AO ao, PetscViewer viewer)
 20: {
 21:   PetscMPIInt        rank, size;
 22:   AO_MemoryScalable *aomems = (AO_MemoryScalable *)ao->data;
 23:   PetscBool          iascii;
 24:   PetscMPIInt        tag_app, tag_petsc;
 25:   PetscLayout        map = aomems->map;
 26:   PetscInt          *app, *app_loc, *petsc, *petsc_loc, len, i, j;
 27:   MPI_Status         status;

 29:   PetscFunctionBegin;
 30:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
 31:   PetscCheck(iascii, PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Viewer type %s not supported for AO MemoryScalable", ((PetscObject)viewer)->type_name);

 33:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ao), &rank));
 34:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ao), &size));

 36:   PetscCall(PetscObjectGetNewTag((PetscObject)ao, &tag_app));
 37:   PetscCall(PetscObjectGetNewTag((PetscObject)ao, &tag_petsc));

 39:   if (rank == 0) {
 40:     PetscCall(PetscViewerASCIIPrintf(viewer, "Number of elements in ordering %" PetscInt_FMT "\n", ao->N));
 41:     PetscCall(PetscViewerASCIIPrintf(viewer, "PETSc->App  App->PETSc\n"));

 43:     PetscCall(PetscMalloc2(map->N, &app, map->N, &petsc));
 44:     len = map->n;
 45:     /* print local AO */
 46:     PetscCall(PetscViewerASCIIPrintf(viewer, "Process [%d]\n", rank));
 47:     for (i = 0; i < len; i++) PetscCall(PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT "  %3" PetscInt_FMT "    %3" PetscInt_FMT "  %3" PetscInt_FMT "\n", i, aomems->app_loc[i], i, aomems->petsc_loc[i]));

 49:     /* recv and print off-processor's AO */
 50:     for (i = 1; i < size; i++) {
 51:       len       = map->range[i + 1] - map->range[i];
 52:       app_loc   = app + map->range[i];
 53:       petsc_loc = petsc + map->range[i];
 54:       PetscCallMPI(MPI_Recv(app_loc, (PetscMPIInt)len, MPIU_INT, i, tag_app, PetscObjectComm((PetscObject)ao), &status));
 55:       PetscCallMPI(MPI_Recv(petsc_loc, (PetscMPIInt)len, MPIU_INT, i, tag_petsc, PetscObjectComm((PetscObject)ao), &status));
 56:       PetscCall(PetscViewerASCIIPrintf(viewer, "Process [%" PetscInt_FMT "]\n", i));
 57:       for (j = 0; j < len; j++) PetscCall(PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT "  %3" PetscInt_FMT "    %3" PetscInt_FMT "  %3" PetscInt_FMT "\n", map->range[i] + j, app_loc[j], map->range[i] + j, petsc_loc[j]));
 58:     }
 59:     PetscCall(PetscFree2(app, petsc));

 61:   } else {
 62:     /* send values */
 63:     PetscCallMPI(MPI_Send((void *)aomems->app_loc, map->n, MPIU_INT, 0, tag_app, PetscObjectComm((PetscObject)ao)));
 64:     PetscCallMPI(MPI_Send((void *)aomems->petsc_loc, map->n, MPIU_INT, 0, tag_petsc, PetscObjectComm((PetscObject)ao)));
 65:   }
 66:   PetscCall(PetscViewerFlush(viewer));
 67:   PetscFunctionReturn(PETSC_SUCCESS);
 68: }

 70: PetscErrorCode AODestroy_MemoryScalable(AO ao)
 71: {
 72:   AO_MemoryScalable *aomems = (AO_MemoryScalable *)ao->data;

 74:   PetscFunctionBegin;
 75:   PetscCall(PetscFree2(aomems->app_loc, aomems->petsc_loc));
 76:   PetscCall(PetscLayoutDestroy(&aomems->map));
 77:   PetscCall(PetscFree(aomems));
 78:   PetscFunctionReturn(PETSC_SUCCESS);
 79: }

 81: /*
 82:    Input Parameters:
 83: +   ao - the application ordering context
 84: .   n  - the number of integers in ia[]
 85: .   ia - the integers; these are replaced with their mapped value
 86: -   maploc - app_loc or petsc_loc in struct "AO_MemoryScalable"

 88:    Output Parameter:
 89: .   ia - the mapped interges
 90:  */
 91: PetscErrorCode AOMap_MemoryScalable_private(AO ao, PetscInt n, PetscInt *ia, const PetscInt *maploc)
 92: {
 93:   AO_MemoryScalable *aomems = (AO_MemoryScalable *)ao->data;
 94:   MPI_Comm           comm;
 95:   PetscMPIInt        rank, size, tag1, tag2;
 96:   PetscInt          *owner, *start, *sizes, nsends, nreceives;
 97:   PetscInt           nmax, count, *sindices, *rindices, i, j, idx, lastidx, *sindices2, *rindices2;
 98:   const PetscInt    *owners = aomems->map->range;
 99:   MPI_Request       *send_waits, *recv_waits, *send_waits2, *recv_waits2;
100:   MPI_Status         recv_status;
101:   PetscMPIInt        nindices, source, widx;
102:   PetscInt          *rbuf, *sbuf;
103:   MPI_Status        *send_status, *send_status2;

105:   PetscFunctionBegin;
106:   PetscCall(PetscObjectGetComm((PetscObject)ao, &comm));
107:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
108:   PetscCallMPI(MPI_Comm_size(comm, &size));

110:   /*  first count number of contributors to each processor */
111:   PetscCall(PetscMalloc1(size, &start));
112:   PetscCall(PetscCalloc2(2 * size, &sizes, n, &owner));

114:   j       = 0;
115:   lastidx = -1;
116:   for (i = 0; i < n; i++) {
117:     if (ia[i] < 0) owner[i] = -1;      /* mark negative entries (which are not to be mapped) with a special negative value */
118:     if (ia[i] >= ao->N) owner[i] = -2; /* mark out of range entries with special negative value */
119:     else {
120:       /* if indices are NOT locally sorted, need to start search at the beginning */
121:       if (lastidx > (idx = ia[i])) j = 0;
122:       lastidx = idx;
123:       for (; j < size; j++) {
124:         if (idx >= owners[j] && idx < owners[j + 1]) {
125:           sizes[2 * j]++;       /* num of indices to be sent */
126:           sizes[2 * j + 1] = 1; /* send to proc[j] */
127:           owner[i]         = j;
128:           break;
129:         }
130:       }
131:     }
132:   }
133:   sizes[2 * rank] = sizes[2 * rank + 1] = 0; /* do not receive from self! */
134:   nsends                                = 0;
135:   for (i = 0; i < size; i++) nsends += sizes[2 * i + 1];

137:   /* inform other processors of number of messages and max length*/
138:   PetscCall(PetscMaxSum(comm, sizes, &nmax, &nreceives));

140:   /* allocate arrays */
141:   PetscCall(PetscObjectGetNewTag((PetscObject)ao, &tag1));
142:   PetscCall(PetscObjectGetNewTag((PetscObject)ao, &tag2));

144:   PetscCall(PetscMalloc2(nreceives * nmax, &rindices, nreceives, &recv_waits));
145:   PetscCall(PetscMalloc2(nsends * nmax, &rindices2, nsends, &recv_waits2));

147:   PetscCall(PetscMalloc3(n, &sindices, nsends, &send_waits, nsends, &send_status));
148:   PetscCall(PetscMalloc3(n, &sindices2, nreceives, &send_waits2, nreceives, &send_status2));

150:   /* post 1st receives: receive others requests
151:      since we don't know how long each individual message is we
152:      allocate the largest needed buffer for each receive. Potentially
153:      this is a lot of wasted space.
154:   */
155:   for (i = 0, count = 0; i < nreceives; i++) PetscCallMPI(MPI_Irecv(rindices + nmax * i, nmax, MPIU_INT, MPI_ANY_SOURCE, tag1, comm, recv_waits + count++));

157:   /* do 1st sends:
158:       1) starts[i] gives the starting index in svalues for stuff going to
159:          the ith processor
160:   */
161:   start[0] = 0;
162:   for (i = 1; i < size; i++) start[i] = start[i - 1] + sizes[2 * i - 2];
163:   for (i = 0; i < n; i++) {
164:     j = owner[i];
165:     if (j == -1) continue; /* do not remap negative entries in ia[] */
166:     else if (j == -2) { /* out of range entries get mapped to -1 */ ia[i] = -1;
167:       continue;
168:     } else if (j != rank) {
169:       sindices[start[j]++] = ia[i];
170:     } else { /* compute my own map */
171:       ia[i] = maploc[ia[i] - owners[rank]];
172:     }
173:   }

175:   start[0] = 0;
176:   for (i = 1; i < size; i++) start[i] = start[i - 1] + sizes[2 * i - 2];
177:   for (i = 0, count = 0; i < size; i++) {
178:     if (sizes[2 * i + 1]) {
179:       /* send my request to others */
180:       PetscCallMPI(MPI_Isend(sindices + start[i], sizes[2 * i], MPIU_INT, i, tag1, comm, send_waits + count));
181:       /* post receive for the answer of my request */
182:       PetscCallMPI(MPI_Irecv(sindices2 + start[i], sizes[2 * i], MPIU_INT, i, tag2, comm, recv_waits2 + count));
183:       count++;
184:     }
185:   }
186:   PetscCheck(nsends == count, comm, PETSC_ERR_SUP, "nsends %" PetscInt_FMT " != count %" PetscInt_FMT, nsends, count);

188:   /* wait on 1st sends */
189:   if (nsends) PetscCallMPI(MPI_Waitall(nsends, send_waits, send_status));

191:   /* 1st recvs: other's requests */
192:   for (j = 0; j < nreceives; j++) {
193:     PetscCallMPI(MPI_Waitany(nreceives, recv_waits, &widx, &recv_status)); /* idx: index of handle for operation that completed */
194:     PetscCallMPI(MPI_Get_count(&recv_status, MPIU_INT, &nindices));
195:     rbuf   = rindices + nmax * widx; /* global index */
196:     source = recv_status.MPI_SOURCE;

198:     /* compute mapping */
199:     sbuf = rbuf;
200:     for (i = 0; i < nindices; i++) sbuf[i] = maploc[rbuf[i] - owners[rank]];

202:     /* send mapping back to the sender */
203:     PetscCallMPI(MPI_Isend(sbuf, nindices, MPIU_INT, source, tag2, comm, send_waits2 + widx));
204:   }

206:   /* wait on 2nd sends */
207:   if (nreceives) PetscCallMPI(MPI_Waitall(nreceives, send_waits2, send_status2));

209:   /* 2nd recvs: for the answer of my request */
210:   for (j = 0; j < nsends; j++) {
211:     PetscCallMPI(MPI_Waitany(nsends, recv_waits2, &widx, &recv_status));
212:     PetscCallMPI(MPI_Get_count(&recv_status, MPIU_INT, &nindices));
213:     source = recv_status.MPI_SOURCE;
214:     /* pack output ia[] */
215:     rbuf  = sindices2 + start[source];
216:     count = 0;
217:     for (i = 0; i < n; i++) {
218:       if (source == owner[i]) ia[i] = rbuf[count++];
219:     }
220:   }

222:   /* free arrays */
223:   PetscCall(PetscFree(start));
224:   PetscCall(PetscFree2(sizes, owner));
225:   PetscCall(PetscFree2(rindices, recv_waits));
226:   PetscCall(PetscFree2(rindices2, recv_waits2));
227:   PetscCall(PetscFree3(sindices, send_waits, send_status));
228:   PetscCall(PetscFree3(sindices2, send_waits2, send_status2));
229:   PetscFunctionReturn(PETSC_SUCCESS);
230: }

232: PetscErrorCode AOPetscToApplication_MemoryScalable(AO ao, PetscInt n, PetscInt *ia)
233: {
234:   AO_MemoryScalable *aomems  = (AO_MemoryScalable *)ao->data;
235:   PetscInt          *app_loc = aomems->app_loc;

237:   PetscFunctionBegin;
238:   PetscCall(AOMap_MemoryScalable_private(ao, n, ia, app_loc));
239:   PetscFunctionReturn(PETSC_SUCCESS);
240: }

242: PetscErrorCode AOApplicationToPetsc_MemoryScalable(AO ao, PetscInt n, PetscInt *ia)
243: {
244:   AO_MemoryScalable *aomems    = (AO_MemoryScalable *)ao->data;
245:   PetscInt          *petsc_loc = aomems->petsc_loc;

247:   PetscFunctionBegin;
248:   PetscCall(AOMap_MemoryScalable_private(ao, n, ia, petsc_loc));
249:   PetscFunctionReturn(PETSC_SUCCESS);
250: }

252: static struct _AOOps AOOps_MemoryScalable = {
253:   PetscDesignatedInitializer(view, AOView_MemoryScalable),
254:   PetscDesignatedInitializer(destroy, AODestroy_MemoryScalable),
255:   PetscDesignatedInitializer(petsctoapplication, AOPetscToApplication_MemoryScalable),
256:   PetscDesignatedInitializer(applicationtopetsc, AOApplicationToPetsc_MemoryScalable),
257: };

259: PetscErrorCode AOCreateMemoryScalable_private(MPI_Comm comm, PetscInt napp, const PetscInt from_array[], const PetscInt to_array[], AO ao, PetscInt *aomap_loc)
260: {
261:   AO_MemoryScalable *aomems  = (AO_MemoryScalable *)ao->data;
262:   PetscLayout        map     = aomems->map;
263:   PetscInt           n_local = map->n, i, j;
264:   PetscMPIInt        rank, size, tag;
265:   PetscInt          *owner, *start, *sizes, nsends, nreceives;
266:   PetscInt           nmax, count, *sindices, *rindices, idx, lastidx;
267:   PetscInt          *owners = aomems->map->range;
268:   MPI_Request       *send_waits, *recv_waits;
269:   MPI_Status         recv_status;
270:   PetscMPIInt        nindices, widx;
271:   PetscInt          *rbuf;
272:   PetscInt           n = napp, ip, ia;
273:   MPI_Status        *send_status;

275:   PetscFunctionBegin;
276:   PetscCall(PetscArrayzero(aomap_loc, n_local));

278:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
279:   PetscCallMPI(MPI_Comm_size(comm, &size));

281:   /*  first count number of contributors (of from_array[]) to each processor */
282:   PetscCall(PetscCalloc1(2 * size, &sizes));
283:   PetscCall(PetscMalloc1(n, &owner));

285:   j       = 0;
286:   lastidx = -1;
287:   for (i = 0; i < n; i++) {
288:     /* if indices are NOT locally sorted, need to start search at the beginning */
289:     if (lastidx > (idx = from_array[i])) j = 0;
290:     lastidx = idx;
291:     for (; j < size; j++) {
292:       if (idx >= owners[j] && idx < owners[j + 1]) {
293:         sizes[2 * j] += 2;    /* num of indices to be sent - in pairs (ip,ia) */
294:         sizes[2 * j + 1] = 1; /* send to proc[j] */
295:         owner[i]         = j;
296:         break;
297:       }
298:     }
299:   }
300:   sizes[2 * rank] = sizes[2 * rank + 1] = 0; /* do not receive from self! */
301:   nsends                                = 0;
302:   for (i = 0; i < size; i++) nsends += sizes[2 * i + 1];

304:   /* inform other processors of number of messages and max length*/
305:   PetscCall(PetscMaxSum(comm, sizes, &nmax, &nreceives));

307:   /* allocate arrays */
308:   PetscCall(PetscObjectGetNewTag((PetscObject)ao, &tag));
309:   PetscCall(PetscMalloc2(nreceives * nmax, &rindices, nreceives, &recv_waits));
310:   PetscCall(PetscMalloc3(2 * n, &sindices, nsends, &send_waits, nsends, &send_status));
311:   PetscCall(PetscMalloc1(size, &start));

313:   /* post receives: */
314:   for (i = 0; i < nreceives; i++) PetscCallMPI(MPI_Irecv(rindices + nmax * i, nmax, MPIU_INT, MPI_ANY_SOURCE, tag, comm, recv_waits + i));

316:   /* do sends:
317:       1) starts[i] gives the starting index in svalues for stuff going to
318:          the ith processor
319:   */
320:   start[0] = 0;
321:   for (i = 1; i < size; i++) start[i] = start[i - 1] + sizes[2 * i - 2];
322:   for (i = 0; i < n; i++) {
323:     j = owner[i];
324:     if (j != rank) {
325:       ip                   = from_array[i];
326:       ia                   = to_array[i];
327:       sindices[start[j]++] = ip;
328:       sindices[start[j]++] = ia;
329:     } else { /* compute my own map */
330:       ip            = from_array[i] - owners[rank];
331:       ia            = to_array[i];
332:       aomap_loc[ip] = ia;
333:     }
334:   }

336:   start[0] = 0;
337:   for (i = 1; i < size; i++) start[i] = start[i - 1] + sizes[2 * i - 2];
338:   for (i = 0, count = 0; i < size; i++) {
339:     if (sizes[2 * i + 1]) {
340:       PetscCallMPI(MPI_Isend(sindices + start[i], sizes[2 * i], MPIU_INT, i, tag, comm, send_waits + count));
341:       count++;
342:     }
343:   }
344:   PetscCheck(nsends == count, comm, PETSC_ERR_SUP, "nsends %" PetscInt_FMT " != count %" PetscInt_FMT, nsends, count);

346:   /* wait on sends */
347:   if (nsends) PetscCallMPI(MPI_Waitall(nsends, send_waits, send_status));

349:   /* recvs */
350:   count = 0;
351:   for (j = nreceives; j > 0; j--) {
352:     PetscCallMPI(MPI_Waitany(nreceives, recv_waits, &widx, &recv_status));
353:     PetscCallMPI(MPI_Get_count(&recv_status, MPIU_INT, &nindices));
354:     rbuf = rindices + nmax * widx; /* global index */

356:     /* compute local mapping */
357:     for (i = 0; i < nindices; i += 2) {       /* pack aomap_loc */
358:       ip            = rbuf[i] - owners[rank]; /* local index */
359:       ia            = rbuf[i + 1];
360:       aomap_loc[ip] = ia;
361:     }
362:     count++;
363:   }

365:   PetscCall(PetscFree(start));
366:   PetscCall(PetscFree3(sindices, send_waits, send_status));
367:   PetscCall(PetscFree2(rindices, recv_waits));
368:   PetscCall(PetscFree(owner));
369:   PetscCall(PetscFree(sizes));
370:   PetscFunctionReturn(PETSC_SUCCESS);
371: }

373: PETSC_EXTERN PetscErrorCode AOCreate_MemoryScalable(AO ao)
374: {
375:   IS                 isapp = ao->isapp, ispetsc = ao->ispetsc;
376:   const PetscInt    *mypetsc, *myapp;
377:   PetscInt           napp, n_local, N, i, start, *petsc, *lens, *disp;
378:   MPI_Comm           comm;
379:   AO_MemoryScalable *aomems;
380:   PetscLayout        map;
381:   PetscMPIInt        size, rank;

383:   PetscFunctionBegin;
384:   PetscCheck(isapp, PetscObjectComm((PetscObject)ao), PETSC_ERR_ARG_WRONGSTATE, "AOSetIS() must be called before AOSetType()");
385:   /* create special struct aomems */
386:   PetscCall(PetscNew(&aomems));
387:   ao->data = (void *)aomems;
388:   PetscCall(PetscMemcpy(ao->ops, &AOOps_MemoryScalable, sizeof(struct _AOOps)));
389:   PetscCall(PetscObjectChangeTypeName((PetscObject)ao, AOMEMORYSCALABLE));

391:   /* transmit all local lengths of isapp to all processors */
392:   PetscCall(PetscObjectGetComm((PetscObject)isapp, &comm));
393:   PetscCallMPI(MPI_Comm_size(comm, &size));
394:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
395:   PetscCall(PetscMalloc2(size, &lens, size, &disp));
396:   PetscCall(ISGetLocalSize(isapp, &napp));
397:   PetscCallMPI(MPI_Allgather(&napp, 1, MPIU_INT, lens, 1, MPIU_INT, comm));

399:   N = 0;
400:   for (i = 0; i < size; i++) {
401:     disp[i] = N;
402:     N += lens[i];
403:   }

405:   /* If ispetsc is 0 then use "natural" numbering */
406:   if (napp) {
407:     if (!ispetsc) {
408:       start = disp[rank];
409:       PetscCall(PetscMalloc1(napp + 1, &petsc));
410:       for (i = 0; i < napp; i++) petsc[i] = start + i;
411:     } else {
412:       PetscCall(ISGetIndices(ispetsc, &mypetsc));
413:       petsc = (PetscInt *)mypetsc;
414:     }
415:   } else {
416:     petsc = NULL;
417:   }

419:   /* create a map with global size N - used to determine the local sizes of ao - shall we use local napp instead of N? */
420:   PetscCall(PetscLayoutCreate(comm, &map));
421:   map->bs = 1;
422:   map->N  = N;
423:   PetscCall(PetscLayoutSetUp(map));

425:   ao->N       = N;
426:   ao->n       = map->n;
427:   aomems->map = map;

429:   /* create distributed indices app_loc: petsc->app and petsc_loc: app->petsc */
430:   n_local = map->n;
431:   PetscCall(PetscCalloc2(n_local, &aomems->app_loc, n_local, &aomems->petsc_loc));
432:   PetscCall(ISGetIndices(isapp, &myapp));

434:   PetscCall(AOCreateMemoryScalable_private(comm, napp, petsc, myapp, ao, aomems->app_loc));
435:   PetscCall(AOCreateMemoryScalable_private(comm, napp, myapp, petsc, ao, aomems->petsc_loc));

437:   PetscCall(ISRestoreIndices(isapp, &myapp));
438:   if (napp) {
439:     if (ispetsc) {
440:       PetscCall(ISRestoreIndices(ispetsc, &mypetsc));
441:     } else {
442:       PetscCall(PetscFree(petsc));
443:     }
444:   }
445:   PetscCall(PetscFree2(lens, disp));
446:   PetscFunctionReturn(PETSC_SUCCESS);
447: }

449: /*@C
450:    AOCreateMemoryScalable - Creates a memory scalable application ordering using two integer arrays.

452:    Collective

454:    Input Parameters:
455: +  comm - MPI communicator that is to share the `AO`
456: .  napp - size of integer arrays
457: .  myapp - integer array that defines an ordering
458: -  mypetsc - integer array that defines another ordering (may be NULL to
459:              indicate the natural ordering, that is 0,1,2,3,...)

461:    Output Parameter:
462: .  aoout - the new application ordering

464:    Level: beginner

466:     Note:
467:     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"
468:     in the indices. Use `AOCreateMapping()` or `AOCreateMappingIS()` if you wish to have "holes" in the indices.
469:     Comparing with `AOCreateBasic()`, this routine trades memory with message communication.

471: .seealso: [](sec_ao), [](sec_scatter), `AO`, `AOCreateMemoryScalableIS()`, `AODestroy()`, `AOPetscToApplication()`, `AOApplicationToPetsc()`
472: @*/
473: PetscErrorCode AOCreateMemoryScalable(MPI_Comm comm, PetscInt napp, const PetscInt myapp[], const PetscInt mypetsc[], AO *aoout)
474: {
475:   IS              isapp, ispetsc;
476:   const PetscInt *app = myapp, *petsc = mypetsc;

478:   PetscFunctionBegin;
479:   PetscCall(ISCreateGeneral(comm, napp, app, PETSC_USE_POINTER, &isapp));
480:   if (mypetsc) {
481:     PetscCall(ISCreateGeneral(comm, napp, petsc, PETSC_USE_POINTER, &ispetsc));
482:   } else {
483:     ispetsc = NULL;
484:   }
485:   PetscCall(AOCreateMemoryScalableIS(isapp, ispetsc, aoout));
486:   PetscCall(ISDestroy(&isapp));
487:   if (mypetsc) PetscCall(ISDestroy(&ispetsc));
488:   PetscFunctionReturn(PETSC_SUCCESS);
489: }

491: /*@C
492:    AOCreateMemoryScalableIS - Creates a memory scalable application ordering using two index sets.

494:    Collective

496:    Input Parameters:
497: +  isapp - index set that defines an ordering
498: -  ispetsc - index set that defines another ordering (may be NULL to use the
499:              natural ordering)

501:    Output Parameter:
502: .  aoout - the new application ordering

504:    Level: beginner

506:     Notes:
507:     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;
508:     that is there cannot be any "holes".

510:     Comparing with `AOCreateBasicIS()`, this routine trades memory with message communication.

512: .seealso: [](sec_ao), [](sec_scatter), `AO`, `AOCreateBasicIS()`, `AOCreateMemoryScalable()`, `AODestroy()`
513: @*/
514: PetscErrorCode AOCreateMemoryScalableIS(IS isapp, IS ispetsc, AO *aoout)
515: {
516:   MPI_Comm comm;
517:   AO       ao;

519:   PetscFunctionBegin;
520:   PetscCall(PetscObjectGetComm((PetscObject)isapp, &comm));
521:   PetscCall(AOCreate(comm, &ao));
522:   PetscCall(AOSetIS(ao, isapp, ispetsc));
523:   PetscCall(AOSetType(ao, AOMEMORYSCALABLE));
524:   PetscCall(AOViewFromOptions(ao, NULL, "-ao_view"));
525:   *aoout = ao;
526:   PetscFunctionReturn(PETSC_SUCCESS);
527: }