Actual source code: isltog.c


  2: #include <petsc/private/isimpl.h>
  3: #include <petsc/private/hashmapi.h>
  4: #include <petscsf.h>
  5: #include <petscviewer.h>

  7: PetscClassId          IS_LTOGM_CLASSID;
  8: static PetscErrorCode ISLocalToGlobalMappingGetBlockInfo_Private(ISLocalToGlobalMapping, PetscInt *, PetscInt **, PetscInt **, PetscInt ***);

 10: typedef struct {
 11:   PetscInt *globals;
 12: } ISLocalToGlobalMapping_Basic;

 14: typedef struct {
 15:   PetscHMapI globalht;
 16: } ISLocalToGlobalMapping_Hash;

 18: /*@C
 19:   ISGetPointRange - Returns a description of the points in an `IS` suitable for traversal

 21:   Not Collective

 23:   Input Parameter:
 24: . pointIS - The `IS` object

 26:   Output Parameters:
 27: + pStart - The first index, see notes
 28: . pEnd   - One past the last index, see notes
 29: - points - The indices, see notes

 31:   Level: intermediate

 33:   Notes:
 34:   If the `IS` contains contiguous indices in an `ISSTRIDE`, then the indices are contained in [pStart, pEnd) and points = `NULL`.
 35:   Otherwise, `pStart = 0`, `pEnd = numIndices`, and points is an array of the indices. This supports the following pattern
 36: .vb
 37:   ISGetPointRange(is, &pStart, &pEnd, &points);
 38:   for (p = pStart; p < pEnd; ++p) {
 39:     const PetscInt point = points ? points[p] : p;
 40:     // use point
 41:   }
 42:   ISRestorePointRange(is, &pstart, &pEnd, &points);
 43: .ve
 44:   Hence the same code can be written for `pointIS` being a `ISSTRIDE` or `ISGENERAL`

 46: .seealso: [](sec_scatter), `IS`, `ISRestorePointRange()`, `ISGetPointSubrange()`, `ISGetIndices()`, `ISCreateStride()`
 47: @*/
 48: PetscErrorCode ISGetPointRange(IS pointIS, PetscInt *pStart, PetscInt *pEnd, const PetscInt **points)
 49: {
 50:   PetscInt  numCells, step = 1;
 51:   PetscBool isStride;

 53:   PetscFunctionBeginHot;
 54:   *pStart = 0;
 55:   *points = NULL;
 56:   PetscCall(ISGetLocalSize(pointIS, &numCells));
 57:   PetscCall(PetscObjectTypeCompare((PetscObject)pointIS, ISSTRIDE, &isStride));
 58:   if (isStride) PetscCall(ISStrideGetInfo(pointIS, pStart, &step));
 59:   *pEnd = *pStart + numCells;
 60:   if (!isStride || step != 1) PetscCall(ISGetIndices(pointIS, points));
 61:   PetscFunctionReturn(PETSC_SUCCESS);
 62: }

 64: /*@C
 65:   ISRestorePointRange - Destroys the traversal description created with `ISGetPointRange()`

 67:   Not Collective

 69:   Input Parameters:
 70: + pointIS - The `IS` object
 71: . pStart  - The first index, from `ISGetPointRange()`
 72: . pEnd    - One past the last index, from `ISGetPointRange()`
 73: - points  - The indices, from `ISGetPointRange()`

 75:   Level: intermediate

 77: .seealso: [](sec_scatter), `IS`, `ISGetPointRange()`, `ISGetPointSubrange()`, `ISGetIndices()`, `ISCreateStride()`
 78: @*/
 79: PetscErrorCode ISRestorePointRange(IS pointIS, PetscInt *pStart, PetscInt *pEnd, const PetscInt **points)
 80: {
 81:   PetscInt  step = 1;
 82:   PetscBool isStride;

 84:   PetscFunctionBeginHot;
 85:   PetscCall(PetscObjectTypeCompare((PetscObject)pointIS, ISSTRIDE, &isStride));
 86:   if (isStride) PetscCall(ISStrideGetInfo(pointIS, pStart, &step));
 87:   if (!isStride || step != 1) PetscCall(ISGetIndices(pointIS, points));
 88:   PetscFunctionReturn(PETSC_SUCCESS);
 89: }

 91: /*@C
 92:   ISGetPointSubrange - Configures the input `IS` to be a subrange for the traversal information given

 94:   Not Collective

 96:   Input Parameters:
 97: + subpointIS - The `IS` object to be configured
 98: . pStart     - The first index of the subrange
 99: . pEnd       - One past the last index for the subrange
100: - points     - The indices for the entire range, from `ISGetPointRange()`

102:   Output Parameters:
103: . subpointIS - The `IS` object now configured to be a subrange

105:   Level: intermediate

107:   Note:
108:   The input `IS` will now respond properly to calls to `ISGetPointRange()` and return the subrange.

110: .seealso: [](sec_scatter), `IS`, `ISGetPointRange()`, `ISRestorePointRange()`, `ISGetIndices()`, `ISCreateStride()`
111: @*/
112: PetscErrorCode ISGetPointSubrange(IS subpointIS, PetscInt pStart, PetscInt pEnd, const PetscInt *points)
113: {
114:   PetscFunctionBeginHot;
115:   if (points) {
116:     PetscCall(ISSetType(subpointIS, ISGENERAL));
117:     PetscCall(ISGeneralSetIndices(subpointIS, pEnd - pStart, &points[pStart], PETSC_USE_POINTER));
118:   } else {
119:     PetscCall(ISSetType(subpointIS, ISSTRIDE));
120:     PetscCall(ISStrideSetStride(subpointIS, pEnd - pStart, pStart, 1));
121:   }
122:   PetscFunctionReturn(PETSC_SUCCESS);
123: }

125: /* -----------------------------------------------------------------------------------------*/

127: /*
128:     Creates the global mapping information in the ISLocalToGlobalMapping structure

130:     If the user has not selected how to handle the global to local mapping then use HASH for "large" problems
131: */
132: static PetscErrorCode ISGlobalToLocalMappingSetUp(ISLocalToGlobalMapping mapping)
133: {
134:   PetscInt i, *idx = mapping->indices, n = mapping->n, end, start;

136:   PetscFunctionBegin;
137:   if (mapping->data) PetscFunctionReturn(PETSC_SUCCESS);
138:   end   = 0;
139:   start = PETSC_MAX_INT;

141:   for (i = 0; i < n; i++) {
142:     if (idx[i] < 0) continue;
143:     if (idx[i] < start) start = idx[i];
144:     if (idx[i] > end) end = idx[i];
145:   }
146:   if (start > end) {
147:     start = 0;
148:     end   = -1;
149:   }
150:   mapping->globalstart = start;
151:   mapping->globalend   = end;
152:   if (!((PetscObject)mapping)->type_name) {
153:     if ((end - start) > PetscMax(4 * n, 1000000)) {
154:       PetscCall(ISLocalToGlobalMappingSetType(mapping, ISLOCALTOGLOBALMAPPINGHASH));
155:     } else {
156:       PetscCall(ISLocalToGlobalMappingSetType(mapping, ISLOCALTOGLOBALMAPPINGBASIC));
157:     }
158:   }
159:   PetscTryTypeMethod(mapping, globaltolocalmappingsetup);
160:   PetscFunctionReturn(PETSC_SUCCESS);
161: }

163: static PetscErrorCode ISGlobalToLocalMappingSetUp_Basic(ISLocalToGlobalMapping mapping)
164: {
165:   PetscInt                      i, *idx = mapping->indices, n = mapping->n, end, start, *globals;
166:   ISLocalToGlobalMapping_Basic *map;

168:   PetscFunctionBegin;
169:   start = mapping->globalstart;
170:   end   = mapping->globalend;
171:   PetscCall(PetscNew(&map));
172:   PetscCall(PetscMalloc1(end - start + 2, &globals));
173:   map->globals = globals;
174:   for (i = 0; i < end - start + 1; i++) globals[i] = -1;
175:   for (i = 0; i < n; i++) {
176:     if (idx[i] < 0) continue;
177:     globals[idx[i] - start] = i;
178:   }
179:   mapping->data = (void *)map;
180:   PetscFunctionReturn(PETSC_SUCCESS);
181: }

183: static PetscErrorCode ISGlobalToLocalMappingSetUp_Hash(ISLocalToGlobalMapping mapping)
184: {
185:   PetscInt                     i, *idx = mapping->indices, n = mapping->n;
186:   ISLocalToGlobalMapping_Hash *map;

188:   PetscFunctionBegin;
189:   PetscCall(PetscNew(&map));
190:   PetscCall(PetscHMapICreate(&map->globalht));
191:   for (i = 0; i < n; i++) {
192:     if (idx[i] < 0) continue;
193:     PetscCall(PetscHMapISet(map->globalht, idx[i], i));
194:   }
195:   mapping->data = (void *)map;
196:   PetscFunctionReturn(PETSC_SUCCESS);
197: }

199: static PetscErrorCode ISLocalToGlobalMappingDestroy_Basic(ISLocalToGlobalMapping mapping)
200: {
201:   ISLocalToGlobalMapping_Basic *map = (ISLocalToGlobalMapping_Basic *)mapping->data;

203:   PetscFunctionBegin;
204:   if (!map) PetscFunctionReturn(PETSC_SUCCESS);
205:   PetscCall(PetscFree(map->globals));
206:   PetscCall(PetscFree(mapping->data));
207:   PetscFunctionReturn(PETSC_SUCCESS);
208: }

210: static PetscErrorCode ISLocalToGlobalMappingDestroy_Hash(ISLocalToGlobalMapping mapping)
211: {
212:   ISLocalToGlobalMapping_Hash *map = (ISLocalToGlobalMapping_Hash *)mapping->data;

214:   PetscFunctionBegin;
215:   if (!map) PetscFunctionReturn(PETSC_SUCCESS);
216:   PetscCall(PetscHMapIDestroy(&map->globalht));
217:   PetscCall(PetscFree(mapping->data));
218:   PetscFunctionReturn(PETSC_SUCCESS);
219: }

221: #define GTOLTYPE _Basic
222: #define GTOLNAME _Basic
223: #define GTOLBS   mapping->bs
224: #define GTOL(g, local) \
225:   do { \
226:     local = map->globals[g / bs - start]; \
227:     if (local >= 0) local = bs * local + (g % bs); \
228:   } while (0)

230: #include <../src/vec/is/utils/isltog.h>

232: #define GTOLTYPE _Basic
233: #define GTOLNAME Block_Basic
234: #define GTOLBS   1
235: #define GTOL(g, local) \
236:   do { \
237:     local = map->globals[g - start]; \
238:   } while (0)
239: #include <../src/vec/is/utils/isltog.h>

241: #define GTOLTYPE _Hash
242: #define GTOLNAME _Hash
243: #define GTOLBS   mapping->bs
244: #define GTOL(g, local) \
245:   do { \
246:     (void)PetscHMapIGet(map->globalht, g / bs, &local); \
247:     if (local >= 0) local = bs * local + (g % bs); \
248:   } while (0)
249: #include <../src/vec/is/utils/isltog.h>

251: #define GTOLTYPE _Hash
252: #define GTOLNAME Block_Hash
253: #define GTOLBS   1
254: #define GTOL(g, local) \
255:   do { \
256:     (void)PetscHMapIGet(map->globalht, g, &local); \
257:   } while (0)
258: #include <../src/vec/is/utils/isltog.h>

260: /*@
261:     ISLocalToGlobalMappingDuplicate - Duplicates the local to global mapping object

263:     Not Collective

265:     Input Parameter:
266: .   ltog - local to global mapping

268:     Output Parameter:
269: .   nltog - the duplicated local to global mapping

271:     Level: advanced

273: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreate()`
274: @*/
275: PetscErrorCode ISLocalToGlobalMappingDuplicate(ISLocalToGlobalMapping ltog, ISLocalToGlobalMapping *nltog)
276: {
277:   ISLocalToGlobalMappingType l2gtype;

279:   PetscFunctionBegin;
281:   PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)ltog), ltog->bs, ltog->n, ltog->indices, PETSC_COPY_VALUES, nltog));
282:   PetscCall(ISLocalToGlobalMappingGetType(ltog, &l2gtype));
283:   PetscCall(ISLocalToGlobalMappingSetType(*nltog, l2gtype));
284:   PetscFunctionReturn(PETSC_SUCCESS);
285: }

287: /*@
288:     ISLocalToGlobalMappingGetSize - Gets the local size of a local to global mapping

290:     Not Collective

292:     Input Parameter:
293: .   ltog - local to global mapping

295:     Output Parameter:
296: .   n - the number of entries in the local mapping, `ISLocalToGlobalMappingGetIndices()` returns an array of this length

298:     Level: advanced

300: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreate()`
301: @*/
302: PetscErrorCode ISLocalToGlobalMappingGetSize(ISLocalToGlobalMapping mapping, PetscInt *n)
303: {
304:   PetscFunctionBegin;
307:   *n = mapping->bs * mapping->n;
308:   PetscFunctionReturn(PETSC_SUCCESS);
309: }

311: /*@C
312:    ISLocalToGlobalMappingViewFromOptions - View an `ISLocalToGlobalMapping` based on values in the options database

314:    Collective

316:    Input Parameters:
317: +  A - the local to global mapping object
318: .  obj - Optional object that provides the options prefix used for the options database query
319: -  name - command line option

321:    Level: intermediate

323:    Note:
324:    See `PetscObjectViewFromOptions()` for the available `PetscViewer` and `PetscViewerFormat`

326: .seealso: [](sec_scatter), `PetscViewer`, ``ISLocalToGlobalMapping`, `ISLocalToGlobalMappingView`, `PetscObjectViewFromOptions()`, `ISLocalToGlobalMappingCreate()`
327: @*/
328: PetscErrorCode ISLocalToGlobalMappingViewFromOptions(ISLocalToGlobalMapping A, PetscObject obj, const char name[])
329: {
330:   PetscFunctionBegin;
332:   PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
333:   PetscFunctionReturn(PETSC_SUCCESS);
334: }

336: /*@C
337:     ISLocalToGlobalMappingView - View a local to global mapping

339:     Not Collective

341:     Input Parameters:
342: +   ltog - local to global mapping
343: -   viewer - viewer

345:     Level: advanced

347: .seealso: [](sec_scatter), `PetscViewer`, `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreate()`
348: @*/
349: PetscErrorCode ISLocalToGlobalMappingView(ISLocalToGlobalMapping mapping, PetscViewer viewer)
350: {
351:   PetscInt    i;
352:   PetscMPIInt rank;
353:   PetscBool   iascii;

355:   PetscFunctionBegin;
357:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mapping), &viewer));

360:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)mapping), &rank));
361:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
362:   if (iascii) {
363:     PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)mapping, viewer));
364:     PetscCall(PetscViewerASCIIPushSynchronized(viewer));
365:     for (i = 0; i < mapping->n; i++) {
366:       PetscInt bs = mapping->bs, g = mapping->indices[i];
367:       if (bs == 1) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] %" PetscInt_FMT " %" PetscInt_FMT "\n", rank, i, g));
368:       else PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] %" PetscInt_FMT ":%" PetscInt_FMT " %" PetscInt_FMT ":%" PetscInt_FMT "\n", rank, i * bs, (i + 1) * bs, g * bs, (g + 1) * bs));
369:     }
370:     PetscCall(PetscViewerFlush(viewer));
371:     PetscCall(PetscViewerASCIIPopSynchronized(viewer));
372:   }
373:   PetscFunctionReturn(PETSC_SUCCESS);
374: }

376: /*@
377:     ISLocalToGlobalMappingCreateIS - Creates a mapping between a local (0 to n)
378:     ordering and a global parallel ordering.

380:     Not Collective

382:     Input Parameter:
383: .   is - index set containing the global numbers for each local number

385:     Output Parameter:
386: .   mapping - new mapping data structure

388:     Level: advanced

390:     Note:
391:     the block size of the `IS` determines the block size of the mapping

393: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingSetFromOptions()`
394: @*/
395: PetscErrorCode ISLocalToGlobalMappingCreateIS(IS is, ISLocalToGlobalMapping *mapping)
396: {
397:   PetscInt        n, bs;
398:   const PetscInt *indices;
399:   MPI_Comm        comm;
400:   PetscBool       isblock;

402:   PetscFunctionBegin;

406:   PetscCall(PetscObjectGetComm((PetscObject)is, &comm));
407:   PetscCall(ISGetLocalSize(is, &n));
408:   PetscCall(PetscObjectTypeCompare((PetscObject)is, ISBLOCK, &isblock));
409:   if (!isblock) {
410:     PetscCall(ISGetIndices(is, &indices));
411:     PetscCall(ISLocalToGlobalMappingCreate(comm, 1, n, indices, PETSC_COPY_VALUES, mapping));
412:     PetscCall(ISRestoreIndices(is, &indices));
413:   } else {
414:     PetscCall(ISGetBlockSize(is, &bs));
415:     PetscCall(ISBlockGetIndices(is, &indices));
416:     PetscCall(ISLocalToGlobalMappingCreate(comm, bs, n / bs, indices, PETSC_COPY_VALUES, mapping));
417:     PetscCall(ISBlockRestoreIndices(is, &indices));
418:   }
419:   PetscFunctionReturn(PETSC_SUCCESS);
420: }

422: /*@C
423:     ISLocalToGlobalMappingCreateSF - Creates a mapping between a local (0 to n)
424:     ordering and a global parallel ordering.

426:     Collective

428:     Input Parameters:
429: +   sf - star forest mapping contiguous local indices to (rank, offset)
430: -   start - first global index on this process, or `PETSC_DECIDE` to compute contiguous global numbering automatically

432:     Output Parameter:
433: .   mapping - new mapping data structure

435:     Level: advanced

437:     Note:
438:     If any processor calls this with `start` = `PETSC_DECIDE` then all processors must, otherwise the program will hang.

440: .seealso: [](sec_scatter), `PetscSF`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingSetFromOptions()`
441: @*/
442: PetscErrorCode ISLocalToGlobalMappingCreateSF(PetscSF sf, PetscInt start, ISLocalToGlobalMapping *mapping)
443: {
444:   PetscInt i, maxlocal, nroots, nleaves, *globals, *ltog;
445:   MPI_Comm comm;

447:   PetscFunctionBegin;
450:   PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
451:   PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL));
452:   if (start == PETSC_DECIDE) {
453:     start = 0;
454:     PetscCallMPI(MPI_Exscan(&nroots, &start, 1, MPIU_INT, MPI_SUM, comm));
455:   } else PetscCheck(start >= 0, comm, PETSC_ERR_ARG_OUTOFRANGE, "start must be nonnegative or PETSC_DECIDE");
456:   PetscCall(PetscSFGetLeafRange(sf, NULL, &maxlocal));
457:   ++maxlocal;
458:   PetscCall(PetscMalloc1(nroots, &globals));
459:   PetscCall(PetscMalloc1(maxlocal, &ltog));
460:   for (i = 0; i < nroots; i++) globals[i] = start + i;
461:   for (i = 0; i < maxlocal; i++) ltog[i] = -1;
462:   PetscCall(PetscSFBcastBegin(sf, MPIU_INT, globals, ltog, MPI_REPLACE));
463:   PetscCall(PetscSFBcastEnd(sf, MPIU_INT, globals, ltog, MPI_REPLACE));
464:   PetscCall(ISLocalToGlobalMappingCreate(comm, 1, maxlocal, ltog, PETSC_OWN_POINTER, mapping));
465:   PetscCall(PetscFree(globals));
466:   PetscFunctionReturn(PETSC_SUCCESS);
467: }

469: /*@
470:     ISLocalToGlobalMappingSetBlockSize - Sets the blocksize of the mapping

472:     Not Collective

474:     Input Parameters:
475: +   mapping - mapping data structure
476: -   bs - the blocksize

478:     Level: advanced

480: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`
481: @*/
482: PetscErrorCode ISLocalToGlobalMappingSetBlockSize(ISLocalToGlobalMapping mapping, PetscInt bs)
483: {
484:   PetscInt       *nid;
485:   const PetscInt *oid;
486:   PetscInt        i, cn, on, obs, nn;

488:   PetscFunctionBegin;
490:   PetscCheck(bs >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid block size %" PetscInt_FMT, bs);
491:   if (bs == mapping->bs) PetscFunctionReturn(PETSC_SUCCESS);
492:   on  = mapping->n;
493:   obs = mapping->bs;
494:   oid = mapping->indices;
495:   nn  = (on * obs) / bs;
496:   PetscCheck((on * obs) % bs == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Block size %" PetscInt_FMT " is inconsistent with block size %" PetscInt_FMT " and number of block indices %" PetscInt_FMT, bs, obs, on);

498:   PetscCall(PetscMalloc1(nn, &nid));
499:   PetscCall(ISLocalToGlobalMappingGetIndices(mapping, &oid));
500:   for (i = 0; i < nn; i++) {
501:     PetscInt j;
502:     for (j = 0, cn = 0; j < bs - 1; j++) {
503:       if (oid[i * bs + j] < 0) {
504:         cn++;
505:         continue;
506:       }
507:       PetscCheck(oid[i * bs + j] == oid[i * bs + j + 1] - 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Block sizes %" PetscInt_FMT " and %" PetscInt_FMT " are incompatible with the block indices: non consecutive indices %" PetscInt_FMT " %" PetscInt_FMT, bs, obs, oid[i * bs + j], oid[i * bs + j + 1]);
508:     }
509:     if (oid[i * bs + j] < 0) cn++;
510:     if (cn) {
511:       PetscCheck(cn == bs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Block sizes %" PetscInt_FMT " and %" PetscInt_FMT " are incompatible with the block indices: invalid number of negative entries in block %" PetscInt_FMT, bs, obs, cn);
512:       nid[i] = -1;
513:     } else {
514:       nid[i] = oid[i * bs] / bs;
515:     }
516:   }
517:   PetscCall(ISLocalToGlobalMappingRestoreIndices(mapping, &oid));

519:   mapping->n  = nn;
520:   mapping->bs = bs;
521:   PetscCall(PetscFree(mapping->indices));
522:   mapping->indices     = nid;
523:   mapping->globalstart = 0;
524:   mapping->globalend   = 0;

526:   /* reset the cached information */
527:   PetscCall(PetscFree(mapping->info_procs));
528:   PetscCall(PetscFree(mapping->info_numprocs));
529:   if (mapping->info_indices) {
530:     PetscInt i;

532:     PetscCall(PetscFree((mapping->info_indices)[0]));
533:     for (i = 1; i < mapping->info_nproc; i++) PetscCall(PetscFree(mapping->info_indices[i]));
534:     PetscCall(PetscFree(mapping->info_indices));
535:   }
536:   mapping->info_cached = PETSC_FALSE;

538:   PetscTryTypeMethod(mapping, destroy);
539:   PetscFunctionReturn(PETSC_SUCCESS);
540: }

542: /*@
543:     ISLocalToGlobalMappingGetBlockSize - Gets the blocksize of the mapping
544:     ordering and a global parallel ordering.

546:     Not Collective

548:     Input Parameter:
549: .   mapping - mapping data structure

551:     Output Parameter:
552: .   bs - the blocksize

554:     Level: advanced

556: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`
557: @*/
558: PetscErrorCode ISLocalToGlobalMappingGetBlockSize(ISLocalToGlobalMapping mapping, PetscInt *bs)
559: {
560:   PetscFunctionBegin;
562:   *bs = mapping->bs;
563:   PetscFunctionReturn(PETSC_SUCCESS);
564: }

566: /*@
567:     ISLocalToGlobalMappingCreate - Creates a mapping between a local (0 to n)
568:     ordering and a global parallel ordering.

570:     Not Collective, but communicator may have more than one process

572:     Input Parameters:
573: +   comm - MPI communicator
574: .   bs - the block size
575: .   n - the number of local elements divided by the block size, or equivalently the number of block indices
576: .   indices - the global index for each local element, these do not need to be in increasing order (sorted), these values should not be scaled (i.e. multiplied) by the blocksize bs
577: -   mode - see PetscCopyMode

579:     Output Parameter:
580: .   mapping - new mapping data structure

582:     Level: advanced

584:     Notes:
585:     There is one integer value in indices per block and it represents the actual indices bs*idx + j, where j=0,..,bs-1

587:     For "small" problems when using `ISGlobalToLocalMappingApply()` and `ISGlobalToLocalMappingApplyBlock()`, the `ISLocalToGlobalMappingType`
588:     of `ISLOCALTOGLOBALMAPPINGBASIC` will be used; this uses more memory but is faster; this approach is not scalable for extremely large mappings.

590:     For large problems `ISLOCALTOGLOBALMAPPINGHASH` is used, this is scalable.
591:     Use `ISLocalToGlobalMappingSetType()` or call `ISLocalToGlobalMappingSetFromOptions()` with the option
592:     `-islocaltoglobalmapping_type` <`basic`,`hash`> to control which is used.

594: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingSetFromOptions()`,
595:           `ISLOCALTOGLOBALMAPPINGBASIC`, `ISLOCALTOGLOBALMAPPINGHASH`
596:           `ISLocalToGlobalMappingSetType()`, `ISLocalToGlobalMappingType`
597: @*/
598: PetscErrorCode ISLocalToGlobalMappingCreate(MPI_Comm comm, PetscInt bs, PetscInt n, const PetscInt indices[], PetscCopyMode mode, ISLocalToGlobalMapping *mapping)
599: {
600:   PetscInt *in;

602:   PetscFunctionBegin;

606:   *mapping = NULL;
607:   PetscCall(ISInitializePackage());

609:   PetscCall(PetscHeaderCreate(*mapping, IS_LTOGM_CLASSID, "ISLocalToGlobalMapping", "Local to global mapping", "IS", comm, ISLocalToGlobalMappingDestroy, ISLocalToGlobalMappingView));
610:   (*mapping)->n  = n;
611:   (*mapping)->bs = bs;
612:   if (mode == PETSC_COPY_VALUES) {
613:     PetscCall(PetscMalloc1(n, &in));
614:     PetscCall(PetscArraycpy(in, indices, n));
615:     (*mapping)->indices         = in;
616:     (*mapping)->dealloc_indices = PETSC_TRUE;
617:   } else if (mode == PETSC_OWN_POINTER) {
618:     (*mapping)->indices         = (PetscInt *)indices;
619:     (*mapping)->dealloc_indices = PETSC_TRUE;
620:   } else if (mode == PETSC_USE_POINTER) {
621:     (*mapping)->indices = (PetscInt *)indices;
622:   } else SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid mode %d", mode);
623:   PetscFunctionReturn(PETSC_SUCCESS);
624: }

626: PetscFunctionList ISLocalToGlobalMappingList = NULL;

628: /*@
629:    ISLocalToGlobalMappingSetFromOptions - Set mapping options from the options database.

631:    Not Collective

633:    Input Parameter:
634: .  mapping - mapping data structure

636:    Options Database Key:
637: .  -islocaltoglobalmapping_type - <basic,hash> nonscalable and scalable versions

639:    Level: advanced

641: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingSetFromOptions()`,
642:           `ISLOCALTOGLOBALMAPPINGBASIC`, `ISLOCALTOGLOBALMAPPINGHASH`
643:           `ISLocalToGlobalMappingSetType()`, `ISLocalToGlobalMappingType`
644: @*/
645: PetscErrorCode ISLocalToGlobalMappingSetFromOptions(ISLocalToGlobalMapping mapping)
646: {
647:   char                       type[256];
648:   ISLocalToGlobalMappingType defaulttype = "Not set";
649:   PetscBool                  flg;

651:   PetscFunctionBegin;
653:   PetscCall(ISLocalToGlobalMappingRegisterAll());
654:   PetscObjectOptionsBegin((PetscObject)mapping);
655:   PetscCall(PetscOptionsFList("-islocaltoglobalmapping_type", "ISLocalToGlobalMapping method", "ISLocalToGlobalMappingSetType", ISLocalToGlobalMappingList, (char *)(((PetscObject)mapping)->type_name) ? ((PetscObject)mapping)->type_name : defaulttype, type, 256, &flg));
656:   if (flg) PetscCall(ISLocalToGlobalMappingSetType(mapping, type));
657:   PetscOptionsEnd();
658:   PetscFunctionReturn(PETSC_SUCCESS);
659: }

661: /*@
662:    ISLocalToGlobalMappingDestroy - Destroys a mapping between a local (0 to n)
663:    ordering and a global parallel ordering.

665:    Not Collective

667:    Input Parameter:
668: .  mapping - mapping data structure

670:    Level: advanced

672: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingCreate()`
673: @*/
674: PetscErrorCode ISLocalToGlobalMappingDestroy(ISLocalToGlobalMapping *mapping)
675: {
676:   PetscFunctionBegin;
677:   if (!*mapping) PetscFunctionReturn(PETSC_SUCCESS);
679:   if (--((PetscObject)(*mapping))->refct > 0) {
680:     *mapping = NULL;
681:     PetscFunctionReturn(PETSC_SUCCESS);
682:   }
683:   if ((*mapping)->dealloc_indices) PetscCall(PetscFree((*mapping)->indices));
684:   PetscCall(PetscFree((*mapping)->info_procs));
685:   PetscCall(PetscFree((*mapping)->info_numprocs));
686:   if ((*mapping)->info_indices) {
687:     PetscInt i;

689:     PetscCall(PetscFree(((*mapping)->info_indices)[0]));
690:     for (i = 1; i < (*mapping)->info_nproc; i++) PetscCall(PetscFree(((*mapping)->info_indices)[i]));
691:     PetscCall(PetscFree((*mapping)->info_indices));
692:   }
693:   if ((*mapping)->info_nodei) PetscCall(PetscFree(((*mapping)->info_nodei)[0]));
694:   PetscCall(PetscFree2((*mapping)->info_nodec, (*mapping)->info_nodei));
695:   if ((*mapping)->ops->destroy) PetscCall((*(*mapping)->ops->destroy)(*mapping));
696:   PetscCall(PetscHeaderDestroy(mapping));
697:   *mapping = NULL;
698:   PetscFunctionReturn(PETSC_SUCCESS);
699: }

701: /*@
702:     ISLocalToGlobalMappingApplyIS - Creates from an `IS` in the local numbering
703:     a new index set using the global numbering defined in an `ISLocalToGlobalMapping`
704:     context.

706:     Collective

708:     Input Parameters:
709: +   mapping - mapping between local and global numbering
710: -   is - index set in local numbering

712:     Output Parameter:
713: .   newis - index set in global numbering

715:     Level: advanced

717:     Note:
718:     The output `IS` will have the same communicator as the input `IS`.

720: .seealso: [](sec_scatter), `ISLocalToGlobalMappingApply()`, `ISLocalToGlobalMappingCreate()`,
721:           `ISLocalToGlobalMappingDestroy()`, `ISGlobalToLocalMappingApply()`
722: @*/
723: PetscErrorCode ISLocalToGlobalMappingApplyIS(ISLocalToGlobalMapping mapping, IS is, IS *newis)
724: {
725:   PetscInt        n, *idxout;
726:   const PetscInt *idxin;

728:   PetscFunctionBegin;

733:   PetscCall(ISGetLocalSize(is, &n));
734:   PetscCall(ISGetIndices(is, &idxin));
735:   PetscCall(PetscMalloc1(n, &idxout));
736:   PetscCall(ISLocalToGlobalMappingApply(mapping, n, idxin, idxout));
737:   PetscCall(ISRestoreIndices(is, &idxin));
738:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is), n, idxout, PETSC_OWN_POINTER, newis));
739:   PetscFunctionReturn(PETSC_SUCCESS);
740: }

742: /*@
743:    ISLocalToGlobalMappingApply - Takes a list of integers in a local numbering
744:    and converts them to the global numbering.

746:    Not Collective

748:    Input Parameters:
749: +  mapping - the local to global mapping context
750: .  N - number of integers
751: -  in - input indices in local numbering

753:    Output Parameter:
754: .  out - indices in global numbering

756:    Level: advanced

758:    Note:
759:    The `in` and `out` array parameters may be identical.

761: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingApplyBlock()`, `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingDestroy()`,
762:           `ISLocalToGlobalMappingApplyIS()`, `AOCreateBasic()`, `AOApplicationToPetsc()`,
763:           `AOPetscToApplication()`, `ISGlobalToLocalMappingApply()`
764: @*/
765: PetscErrorCode ISLocalToGlobalMappingApply(ISLocalToGlobalMapping mapping, PetscInt N, const PetscInt in[], PetscInt out[])
766: {
767:   PetscInt i, bs, Nmax;

769:   PetscFunctionBegin;
771:   bs   = mapping->bs;
772:   Nmax = bs * mapping->n;
773:   if (bs == 1) {
774:     const PetscInt *idx = mapping->indices;
775:     for (i = 0; i < N; i++) {
776:       if (in[i] < 0) {
777:         out[i] = in[i];
778:         continue;
779:       }
780:       PetscCheck(in[i] < Nmax, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local index %" PetscInt_FMT " too large %" PetscInt_FMT " (max) at %" PetscInt_FMT, in[i], Nmax - 1, i);
781:       out[i] = idx[in[i]];
782:     }
783:   } else {
784:     const PetscInt *idx = mapping->indices;
785:     for (i = 0; i < N; i++) {
786:       if (in[i] < 0) {
787:         out[i] = in[i];
788:         continue;
789:       }
790:       PetscCheck(in[i] < Nmax, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local index %" PetscInt_FMT " too large %" PetscInt_FMT " (max) at %" PetscInt_FMT, in[i], Nmax - 1, i);
791:       out[i] = idx[in[i] / bs] * bs + (in[i] % bs);
792:     }
793:   }
794:   PetscFunctionReturn(PETSC_SUCCESS);
795: }

797: /*@
798:    ISLocalToGlobalMappingApplyBlock - Takes a list of integers in a local block numbering and converts them to the global block numbering

800:    Not Collective

802:    Input Parameters:
803: +  mapping - the local to global mapping context
804: .  N - number of integers
805: -  in - input indices in local block numbering

807:    Output Parameter:
808: .  out - indices in global block numbering

810:    Example:
811:      If the index values are {0,1,6,7} set with a call to `ISLocalToGlobalMappingCreate`(`PETSC_COMM_SELF`,2,2,{0,3}) then the mapping applied to 0
812:      (the first block) would produce 0 and the mapping applied to 1 (the second block) would produce 3.

814:    Level: advanced

816:    Note:
817:    The `in` and `out` array parameters may be identical.

819: .seealso: [](sec_scatter), `ISLocalToGlobalMappingApply()`, `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingDestroy()`,
820:           `ISLocalToGlobalMappingApplyIS()`, `AOCreateBasic()`, `AOApplicationToPetsc()`,
821:           `AOPetscToApplication()`, `ISGlobalToLocalMappingApply()`
822: @*/
823: PetscErrorCode ISLocalToGlobalMappingApplyBlock(ISLocalToGlobalMapping mapping, PetscInt N, const PetscInt in[], PetscInt out[])
824: {
825:   PetscInt        i, Nmax;
826:   const PetscInt *idx;

828:   PetscFunctionBegin;
830:   Nmax = mapping->n;
831:   idx  = mapping->indices;
832:   for (i = 0; i < N; i++) {
833:     if (in[i] < 0) {
834:       out[i] = in[i];
835:       continue;
836:     }
837:     PetscCheck(in[i] < Nmax, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local block index %" PetscInt_FMT " too large %" PetscInt_FMT " (max) at %" PetscInt_FMT, in[i], Nmax - 1, i);
838:     out[i] = idx[in[i]];
839:   }
840:   PetscFunctionReturn(PETSC_SUCCESS);
841: }

843: /*@
844:     ISGlobalToLocalMappingApply - Provides the local numbering for a list of integers
845:     specified with a global numbering.

847:     Not Collective

849:     Input Parameters:
850: +   mapping - mapping between local and global numbering
851: .   type - `IS_GTOLM_MASK` - maps global indices with no local value to -1 in the output list (i.e., mask them)
852:            `IS_GTOLM_DROP` - drops the indices with no local value from the output list
853: .   n - number of global indices to map
854: -   idx - global indices to map

856:     Output Parameters:
857: +   nout - number of indices in output array (if type == `IS_GTOLM_MASK` then nout = n)
858: -   idxout - local index of each global index, one must pass in an array long enough
859:              to hold all the indices. You can call `ISGlobalToLocalMappingApply()` with
860:              idxout == NULL to determine the required length (returned in nout)
861:              and then allocate the required space and call `ISGlobalToLocalMappingApply()`
862:              a second time to set the values.

864:     Level: advanced

866:     Notes:
867:     Either `nout` or `idxout` may be `NULL`. `idx` and `idxout` may be identical.

869:     For "small" problems when using `ISGlobalToLocalMappingApply()` and `ISGlobalToLocalMappingApplyBlock()`, the `ISLocalToGlobalMappingType` of
870:     `ISLOCALTOGLOBALMAPPINGBASIC` will be used;
871:     this uses more memory but is faster; this approach is not scalable for extremely large mappings. For large problems `ISLOCALTOGLOBALMAPPINGHASH` is used, this is scalable.
872:     Use `ISLocalToGlobalMappingSetType()` or call `ISLocalToGlobalMappingSetFromOptions()` with the option -islocaltoglobalmapping_type <basic,hash> to control which is used.

874:     Developer Note:
875:     The manual page states that `idx` and `idxout` may be identical but the calling
876:     sequence declares `idx` as const so it cannot be the same as `idxout`.

878: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingApply()`, `ISGlobalToLocalMappingApplyBlock()`, `ISLocalToGlobalMappingCreate()`,
879:           `ISLocalToGlobalMappingDestroy()`
880: @*/
881: PetscErrorCode ISGlobalToLocalMappingApply(ISLocalToGlobalMapping mapping, ISGlobalToLocalMappingMode type, PetscInt n, const PetscInt idx[], PetscInt *nout, PetscInt idxout[])
882: {
883:   PetscFunctionBegin;
885:   if (!mapping->data) PetscCall(ISGlobalToLocalMappingSetUp(mapping));
886:   PetscUseTypeMethod(mapping, globaltolocalmappingapply, type, n, idx, nout, idxout);
887:   PetscFunctionReturn(PETSC_SUCCESS);
888: }

890: /*@
891:     ISGlobalToLocalMappingApplyIS - Creates from an `IS` in the global numbering
892:     a new index set using the local numbering defined in an `ISLocalToGlobalMapping`
893:     context.

895:     Not Collective

897:     Input Parameters:
898: +   mapping - mapping between local and global numbering
899: .   type - `IS_GTOLM_MASK` - maps global indices with no local value to -1 in the output list (i.e., mask them)
900:            `IS_GTOLM_DROP` - drops the indices with no local value from the output list
901: -   is - index set in global numbering

903:     Output Parameter:
904: .   newis - index set in local numbering

906:     Level: advanced

908:     Note:
909:     The output `IS` will be sequential, as it encodes a purely local operation

911: .seealso: [](sec_scatter), `ISGlobalToLocalMapping`, `ISGlobalToLocalMappingApply()`, `ISLocalToGlobalMappingCreate()`,
912:           `ISLocalToGlobalMappingDestroy()`
913: @*/
914: PetscErrorCode ISGlobalToLocalMappingApplyIS(ISLocalToGlobalMapping mapping, ISGlobalToLocalMappingMode type, IS is, IS *newis)
915: {
916:   PetscInt        n, nout, *idxout;
917:   const PetscInt *idxin;

919:   PetscFunctionBegin;

924:   PetscCall(ISGetLocalSize(is, &n));
925:   PetscCall(ISGetIndices(is, &idxin));
926:   if (type == IS_GTOLM_MASK) {
927:     PetscCall(PetscMalloc1(n, &idxout));
928:   } else {
929:     PetscCall(ISGlobalToLocalMappingApply(mapping, type, n, idxin, &nout, NULL));
930:     PetscCall(PetscMalloc1(nout, &idxout));
931:   }
932:   PetscCall(ISGlobalToLocalMappingApply(mapping, type, n, idxin, &nout, idxout));
933:   PetscCall(ISRestoreIndices(is, &idxin));
934:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nout, idxout, PETSC_OWN_POINTER, newis));
935:   PetscFunctionReturn(PETSC_SUCCESS);
936: }

938: /*@
939:     ISGlobalToLocalMappingApplyBlock - Provides the local block numbering for a list of integers
940:     specified with a block global numbering.

942:     Not Collective

944:     Input Parameters:
945: +   mapping - mapping between local and global numbering
946: .   type - `IS_GTOLM_MASK` - maps global indices with no local value to -1 in the output list (i.e., mask them)
947:            `IS_GTOLM_DROP` - drops the indices with no local value from the output list
948: .   n - number of global indices to map
949: -   idx - global indices to map

951:     Output Parameters:
952: +   nout - number of indices in output array (if type == `IS_GTOLM_MASK` then nout = n)
953: -   idxout - local index of each global index, one must pass in an array long enough
954:              to hold all the indices. You can call `ISGlobalToLocalMappingApplyBlock()` with
955:              idxout == NULL to determine the required length (returned in nout)
956:              and then allocate the required space and call `ISGlobalToLocalMappingApplyBlock()`
957:              a second time to set the values.

959:     Level: advanced

961:     Notes:
962:     Either `nout` or `idxout` may be `NULL`. `idx` and `idxout` may be identical.

964:     For "small" problems when using `ISGlobalToLocalMappingApply()` and `ISGlobalToLocalMappingApplyBlock()`, the `ISLocalToGlobalMappingType` of
965:     `ISLOCALTOGLOBALMAPPINGBASIC` will be used;
966:     this uses more memory but is faster; this approach is not scalable for extremely large mappings. For large problems `ISLOCALTOGLOBALMAPPINGHASH` is used, this is scalable.
967:     Use `ISLocalToGlobalMappingSetType()` or call `ISLocalToGlobalMappingSetFromOptions()` with the option -islocaltoglobalmapping_type <basic,hash> to control which is used.

969:     Developer Note:
970:     The manual page states that `idx` and `idxout` may be identical but the calling
971:     sequence declares `idx` as const so it cannot be the same as `idxout`.

973: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingApply()`, `ISGlobalToLocalMappingApply()`, `ISLocalToGlobalMappingCreate()`,
974:           `ISLocalToGlobalMappingDestroy()`
975: @*/
976: PetscErrorCode ISGlobalToLocalMappingApplyBlock(ISLocalToGlobalMapping mapping, ISGlobalToLocalMappingMode type, PetscInt n, const PetscInt idx[], PetscInt *nout, PetscInt idxout[])
977: {
978:   PetscFunctionBegin;
980:   if (!mapping->data) PetscCall(ISGlobalToLocalMappingSetUp(mapping));
981:   PetscUseTypeMethod(mapping, globaltolocalmappingapplyblock, type, n, idx, nout, idxout);
982:   PetscFunctionReturn(PETSC_SUCCESS);
983: }

985: /*@C
986:     ISLocalToGlobalMappingGetBlockInfo - Gets the neighbor information for each processor and
987:      each index shared by more than one processor

989:     Collective

991:     Input Parameter:
992: .   mapping - the mapping from local to global indexing

994:     Output Parameters:
995: +   nproc - number of processors that are connected to this one
996: .   proc - neighboring processors
997: .   numproc - number of indices for each subdomain (processor)
998: -   indices - indices of nodes (in local numbering) shared with neighbors (sorted by global numbering)

1000:     Level: advanced

1002:     Fortran Usage:
1003: .vb
1004:         PetscInt indices[nproc][numprocmax],ierr)
1005:         ISLocalToGlobalMpngGetInfoSize(ISLocalToGlobalMapping,PetscInt nproc,PetscInt numprocmax,ierr) followed by
1006:         ISLocalToGlobalMappingGetInfo(ISLocalToGlobalMapping,PetscInt nproc, PetscInt procs[nproc],PetscInt numprocs[nproc],
1007: .ve

1009:    Fortran Note:
1010:         There is no `ISLocalToGlobalMappingRestoreInfo()` in Fortran. You must make sure that `procs`[], `numprocs`[] and
1011:         `indices`[][] are large enough arrays, either by allocating them dynamically or defining static ones large enough.

1013: .seealso: [](sec_scatter), `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
1014:           `ISLocalToGlobalMappingRestoreInfo()`
1015: @*/
1016: PetscErrorCode ISLocalToGlobalMappingGetBlockInfo(ISLocalToGlobalMapping mapping, PetscInt *nproc, PetscInt *procs[], PetscInt *numprocs[], PetscInt **indices[])
1017: {
1018:   PetscFunctionBegin;
1020:   if (mapping->info_cached) {
1021:     *nproc    = mapping->info_nproc;
1022:     *procs    = mapping->info_procs;
1023:     *numprocs = mapping->info_numprocs;
1024:     *indices  = mapping->info_indices;
1025:   } else {
1026:     PetscCall(ISLocalToGlobalMappingGetBlockInfo_Private(mapping, nproc, procs, numprocs, indices));
1027:   }
1028:   PetscFunctionReturn(PETSC_SUCCESS);
1029: }

1031: static PetscErrorCode ISLocalToGlobalMappingGetBlockInfo_Private(ISLocalToGlobalMapping mapping, PetscInt *nproc, PetscInt *procs[], PetscInt *numprocs[], PetscInt **indices[])
1032: {
1033:   PetscMPIInt  size, rank, tag1, tag2, tag3, *len, *source, imdex;
1034:   PetscInt     i, n = mapping->n, Ng, ng, max = 0, *lindices = mapping->indices;
1035:   PetscInt    *nprocs, *owner, nsends, *sends, j, *starts, nmax, nrecvs, *recvs, proc;
1036:   PetscInt     cnt, scale, *ownedsenders, *nownedsenders, rstart;
1037:   PetscInt     node, nownedm, nt, *sends2, nsends2, *starts2, *lens2, *dest, nrecvs2, *starts3, *recvs2, k, *bprocs, *tmp;
1038:   PetscInt     first_procs, first_numprocs, *first_indices;
1039:   MPI_Request *recv_waits, *send_waits;
1040:   MPI_Status   recv_status, *send_status, *recv_statuses;
1041:   MPI_Comm     comm;
1042:   PetscBool    debug = PETSC_FALSE;

1044:   PetscFunctionBegin;
1045:   PetscCall(PetscObjectGetComm((PetscObject)mapping, &comm));
1046:   PetscCallMPI(MPI_Comm_size(comm, &size));
1047:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
1048:   if (size == 1) {
1049:     *nproc = 0;
1050:     *procs = NULL;
1051:     PetscCall(PetscNew(numprocs));
1052:     (*numprocs)[0] = 0;
1053:     PetscCall(PetscNew(indices));
1054:     (*indices)[0] = NULL;
1055:     /* save info for reuse */
1056:     mapping->info_nproc    = *nproc;
1057:     mapping->info_procs    = *procs;
1058:     mapping->info_numprocs = *numprocs;
1059:     mapping->info_indices  = *indices;
1060:     mapping->info_cached   = PETSC_TRUE;
1061:     PetscFunctionReturn(PETSC_SUCCESS);
1062:   }

1064:   PetscCall(PetscOptionsGetBool(((PetscObject)mapping)->options, NULL, "-islocaltoglobalmappinggetinfo_debug", &debug, NULL));

1066:   /*
1067:     Notes on ISLocalToGlobalMappingGetBlockInfo

1069:     globally owned node - the nodes that have been assigned to this processor in global
1070:            numbering, just for this routine.

1072:     nontrivial globally owned node - node assigned to this processor that is on a subdomain
1073:            boundary (i.e. is has more than one local owner)

1075:     locally owned node - node that exists on this processors subdomain

1077:     nontrivial locally owned node - node that is not in the interior (i.e. has more than one
1078:            local subdomain
1079:   */
1080:   PetscCall(PetscObjectGetNewTag((PetscObject)mapping, &tag1));
1081:   PetscCall(PetscObjectGetNewTag((PetscObject)mapping, &tag2));
1082:   PetscCall(PetscObjectGetNewTag((PetscObject)mapping, &tag3));

1084:   for (i = 0; i < n; i++) {
1085:     if (lindices[i] > max) max = lindices[i];
1086:   }
1087:   PetscCall(MPIU_Allreduce(&max, &Ng, 1, MPIU_INT, MPI_MAX, comm));
1088:   Ng++;
1089:   PetscCallMPI(MPI_Comm_size(comm, &size));
1090:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
1091:   scale = Ng / size + 1;
1092:   ng    = scale;
1093:   if (rank == size - 1) ng = Ng - scale * (size - 1);
1094:   ng     = PetscMax(1, ng);
1095:   rstart = scale * rank;

1097:   /* determine ownership ranges of global indices */
1098:   PetscCall(PetscMalloc1(2 * size, &nprocs));
1099:   PetscCall(PetscArrayzero(nprocs, 2 * size));

1101:   /* determine owners of each local node  */
1102:   PetscCall(PetscMalloc1(n, &owner));
1103:   for (i = 0; i < n; i++) {
1104:     proc                 = lindices[i] / scale; /* processor that globally owns this index */
1105:     nprocs[2 * proc + 1] = 1;                   /* processor globally owns at least one of ours */
1106:     owner[i]             = proc;
1107:     nprocs[2 * proc]++; /* count of how many that processor globally owns of ours */
1108:   }
1109:   nsends = 0;
1110:   for (i = 0; i < size; i++) nsends += nprocs[2 * i + 1];
1111:   PetscCall(PetscInfo(mapping, "Number of global owners for my local data %" PetscInt_FMT "\n", nsends));

1113:   /* inform other processors of number of messages and max length*/
1114:   PetscCall(PetscMaxSum(comm, nprocs, &nmax, &nrecvs));
1115:   PetscCall(PetscInfo(mapping, "Number of local owners for my global data %" PetscInt_FMT "\n", nrecvs));

1117:   /* post receives for owned rows */
1118:   PetscCall(PetscMalloc1((2 * nrecvs + 1) * (nmax + 1), &recvs));
1119:   PetscCall(PetscMalloc1(nrecvs + 1, &recv_waits));
1120:   for (i = 0; i < nrecvs; i++) PetscCallMPI(MPI_Irecv(recvs + 2 * nmax * i, 2 * nmax, MPIU_INT, MPI_ANY_SOURCE, tag1, comm, recv_waits + i));

1122:   /* pack messages containing lists of local nodes to owners */
1123:   PetscCall(PetscMalloc1(2 * n + 1, &sends));
1124:   PetscCall(PetscMalloc1(size + 1, &starts));
1125:   starts[0] = 0;
1126:   for (i = 1; i < size; i++) starts[i] = starts[i - 1] + 2 * nprocs[2 * i - 2];
1127:   for (i = 0; i < n; i++) {
1128:     sends[starts[owner[i]]++] = lindices[i];
1129:     sends[starts[owner[i]]++] = i;
1130:   }
1131:   PetscCall(PetscFree(owner));
1132:   starts[0] = 0;
1133:   for (i = 1; i < size; i++) starts[i] = starts[i - 1] + 2 * nprocs[2 * i - 2];

1135:   /* send the messages */
1136:   PetscCall(PetscMalloc1(nsends + 1, &send_waits));
1137:   PetscCall(PetscMalloc1(nsends + 1, &dest));
1138:   cnt = 0;
1139:   for (i = 0; i < size; i++) {
1140:     if (nprocs[2 * i]) {
1141:       PetscCallMPI(MPI_Isend(sends + starts[i], 2 * nprocs[2 * i], MPIU_INT, i, tag1, comm, send_waits + cnt));
1142:       dest[cnt] = i;
1143:       cnt++;
1144:     }
1145:   }
1146:   PetscCall(PetscFree(starts));

1148:   /* wait on receives */
1149:   PetscCall(PetscMalloc1(nrecvs + 1, &source));
1150:   PetscCall(PetscMalloc1(nrecvs + 1, &len));
1151:   cnt = nrecvs;
1152:   PetscCall(PetscCalloc1(ng + 1, &nownedsenders));
1153:   while (cnt) {
1154:     PetscCallMPI(MPI_Waitany(nrecvs, recv_waits, &imdex, &recv_status));
1155:     /* unpack receives into our local space */
1156:     PetscCallMPI(MPI_Get_count(&recv_status, MPIU_INT, &len[imdex]));
1157:     source[imdex] = recv_status.MPI_SOURCE;
1158:     len[imdex]    = len[imdex] / 2;
1159:     /* count how many local owners for each of my global owned indices */
1160:     for (i = 0; i < len[imdex]; i++) nownedsenders[recvs[2 * imdex * nmax + 2 * i] - rstart]++;
1161:     cnt--;
1162:   }
1163:   PetscCall(PetscFree(recv_waits));

1165:   /* count how many globally owned indices are on an edge multiplied by how many processors own them. */
1166:   nownedm = 0;
1167:   for (i = 0; i < ng; i++) {
1168:     if (nownedsenders[i] > 1) nownedm += nownedsenders[i];
1169:   }

1171:   /* create single array to contain rank of all local owners of each globally owned index */
1172:   PetscCall(PetscMalloc1(nownedm + 1, &ownedsenders));
1173:   PetscCall(PetscMalloc1(ng + 1, &starts));
1174:   starts[0] = 0;
1175:   for (i = 1; i < ng; i++) {
1176:     if (nownedsenders[i - 1] > 1) starts[i] = starts[i - 1] + nownedsenders[i - 1];
1177:     else starts[i] = starts[i - 1];
1178:   }

1180:   /* for each nontrivial globally owned node list all arriving processors */
1181:   for (i = 0; i < nrecvs; i++) {
1182:     for (j = 0; j < len[i]; j++) {
1183:       node = recvs[2 * i * nmax + 2 * j] - rstart;
1184:       if (nownedsenders[node] > 1) ownedsenders[starts[node]++] = source[i];
1185:     }
1186:   }

1188:   if (debug) { /* -----------------------------------  */
1189:     starts[0] = 0;
1190:     for (i = 1; i < ng; i++) {
1191:       if (nownedsenders[i - 1] > 1) starts[i] = starts[i - 1] + nownedsenders[i - 1];
1192:       else starts[i] = starts[i - 1];
1193:     }
1194:     for (i = 0; i < ng; i++) {
1195:       if (nownedsenders[i] > 1) {
1196:         PetscCall(PetscSynchronizedPrintf(comm, "[%d] global node %" PetscInt_FMT " local owner processors: ", rank, i + rstart));
1197:         for (j = 0; j < nownedsenders[i]; j++) PetscCall(PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", ownedsenders[starts[i] + j]));
1198:         PetscCall(PetscSynchronizedPrintf(comm, "\n"));
1199:       }
1200:     }
1201:     PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));
1202:   } /* -----------------------------------  */

1204:   /* wait on original sends */
1205:   if (nsends) {
1206:     PetscCall(PetscMalloc1(nsends, &send_status));
1207:     PetscCallMPI(MPI_Waitall(nsends, send_waits, send_status));
1208:     PetscCall(PetscFree(send_status));
1209:   }
1210:   PetscCall(PetscFree(send_waits));
1211:   PetscCall(PetscFree(sends));
1212:   PetscCall(PetscFree(nprocs));

1214:   /* pack messages to send back to local owners */
1215:   starts[0] = 0;
1216:   for (i = 1; i < ng; i++) {
1217:     if (nownedsenders[i - 1] > 1) starts[i] = starts[i - 1] + nownedsenders[i - 1];
1218:     else starts[i] = starts[i - 1];
1219:   }
1220:   nsends2 = nrecvs;
1221:   PetscCall(PetscMalloc1(nsends2 + 1, &nprocs)); /* length of each message */
1222:   for (i = 0; i < nrecvs; i++) {
1223:     nprocs[i] = 1;
1224:     for (j = 0; j < len[i]; j++) {
1225:       node = recvs[2 * i * nmax + 2 * j] - rstart;
1226:       if (nownedsenders[node] > 1) nprocs[i] += 2 + nownedsenders[node];
1227:     }
1228:   }
1229:   nt = 0;
1230:   for (i = 0; i < nsends2; i++) nt += nprocs[i];

1232:   PetscCall(PetscMalloc1(nt + 1, &sends2));
1233:   PetscCall(PetscMalloc1(nsends2 + 1, &starts2));

1235:   starts2[0] = 0;
1236:   for (i = 1; i < nsends2; i++) starts2[i] = starts2[i - 1] + nprocs[i - 1];
1237:   /*
1238:      Each message is 1 + nprocs[i] long, and consists of
1239:        (0) the number of nodes being sent back
1240:        (1) the local node number,
1241:        (2) the number of processors sharing it,
1242:        (3) the processors sharing it
1243:   */
1244:   for (i = 0; i < nsends2; i++) {
1245:     cnt                = 1;
1246:     sends2[starts2[i]] = 0;
1247:     for (j = 0; j < len[i]; j++) {
1248:       node = recvs[2 * i * nmax + 2 * j] - rstart;
1249:       if (nownedsenders[node] > 1) {
1250:         sends2[starts2[i]]++;
1251:         sends2[starts2[i] + cnt++] = recvs[2 * i * nmax + 2 * j + 1];
1252:         sends2[starts2[i] + cnt++] = nownedsenders[node];
1253:         PetscCall(PetscArraycpy(&sends2[starts2[i] + cnt], &ownedsenders[starts[node]], nownedsenders[node]));
1254:         cnt += nownedsenders[node];
1255:       }
1256:     }
1257:   }

1259:   /* receive the message lengths */
1260:   nrecvs2 = nsends;
1261:   PetscCall(PetscMalloc1(nrecvs2 + 1, &lens2));
1262:   PetscCall(PetscMalloc1(nrecvs2 + 1, &starts3));
1263:   PetscCall(PetscMalloc1(nrecvs2 + 1, &recv_waits));
1264:   for (i = 0; i < nrecvs2; i++) PetscCallMPI(MPI_Irecv(&lens2[i], 1, MPIU_INT, dest[i], tag2, comm, recv_waits + i));

1266:   /* send the message lengths */
1267:   for (i = 0; i < nsends2; i++) PetscCallMPI(MPI_Send(&nprocs[i], 1, MPIU_INT, source[i], tag2, comm));

1269:   /* wait on receives of lens */
1270:   if (nrecvs2) {
1271:     PetscCall(PetscMalloc1(nrecvs2, &recv_statuses));
1272:     PetscCallMPI(MPI_Waitall(nrecvs2, recv_waits, recv_statuses));
1273:     PetscCall(PetscFree(recv_statuses));
1274:   }
1275:   PetscCall(PetscFree(recv_waits));

1277:   starts3[0] = 0;
1278:   nt         = 0;
1279:   for (i = 0; i < nrecvs2 - 1; i++) {
1280:     starts3[i + 1] = starts3[i] + lens2[i];
1281:     nt += lens2[i];
1282:   }
1283:   if (nrecvs2) nt += lens2[nrecvs2 - 1];

1285:   PetscCall(PetscMalloc1(nt + 1, &recvs2));
1286:   PetscCall(PetscMalloc1(nrecvs2 + 1, &recv_waits));
1287:   for (i = 0; i < nrecvs2; i++) PetscCallMPI(MPI_Irecv(recvs2 + starts3[i], lens2[i], MPIU_INT, dest[i], tag3, comm, recv_waits + i));

1289:   /* send the messages */
1290:   PetscCall(PetscMalloc1(nsends2 + 1, &send_waits));
1291:   for (i = 0; i < nsends2; i++) PetscCallMPI(MPI_Isend(sends2 + starts2[i], nprocs[i], MPIU_INT, source[i], tag3, comm, send_waits + i));

1293:   /* wait on receives */
1294:   if (nrecvs2) {
1295:     PetscCall(PetscMalloc1(nrecvs2, &recv_statuses));
1296:     PetscCallMPI(MPI_Waitall(nrecvs2, recv_waits, recv_statuses));
1297:     PetscCall(PetscFree(recv_statuses));
1298:   }
1299:   PetscCall(PetscFree(recv_waits));
1300:   PetscCall(PetscFree(nprocs));

1302:   if (debug) { /* -----------------------------------  */
1303:     cnt = 0;
1304:     for (i = 0; i < nrecvs2; i++) {
1305:       nt = recvs2[cnt++];
1306:       for (j = 0; j < nt; j++) {
1307:         PetscCall(PetscSynchronizedPrintf(comm, "[%d] local node %" PetscInt_FMT " number of subdomains %" PetscInt_FMT ": ", rank, recvs2[cnt], recvs2[cnt + 1]));
1308:         for (k = 0; k < recvs2[cnt + 1]; k++) PetscCall(PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", recvs2[cnt + 2 + k]));
1309:         cnt += 2 + recvs2[cnt + 1];
1310:         PetscCall(PetscSynchronizedPrintf(comm, "\n"));
1311:       }
1312:     }
1313:     PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));
1314:   } /* -----------------------------------  */

1316:   /* count number subdomains for each local node */
1317:   PetscCall(PetscCalloc1(size, &nprocs));
1318:   cnt = 0;
1319:   for (i = 0; i < nrecvs2; i++) {
1320:     nt = recvs2[cnt++];
1321:     for (j = 0; j < nt; j++) {
1322:       for (k = 0; k < recvs2[cnt + 1]; k++) nprocs[recvs2[cnt + 2 + k]]++;
1323:       cnt += 2 + recvs2[cnt + 1];
1324:     }
1325:   }
1326:   nt = 0;
1327:   for (i = 0; i < size; i++) nt += (nprocs[i] > 0);
1328:   *nproc = nt;
1329:   PetscCall(PetscMalloc1(nt + 1, procs));
1330:   PetscCall(PetscMalloc1(nt + 1, numprocs));
1331:   PetscCall(PetscMalloc1(nt + 1, indices));
1332:   for (i = 0; i < nt + 1; i++) (*indices)[i] = NULL;
1333:   PetscCall(PetscMalloc1(size, &bprocs));
1334:   cnt = 0;
1335:   for (i = 0; i < size; i++) {
1336:     if (nprocs[i] > 0) {
1337:       bprocs[i]        = cnt;
1338:       (*procs)[cnt]    = i;
1339:       (*numprocs)[cnt] = nprocs[i];
1340:       PetscCall(PetscMalloc1(nprocs[i], &(*indices)[cnt]));
1341:       cnt++;
1342:     }
1343:   }

1345:   /* make the list of subdomains for each nontrivial local node */
1346:   PetscCall(PetscArrayzero(*numprocs, nt));
1347:   cnt = 0;
1348:   for (i = 0; i < nrecvs2; i++) {
1349:     nt = recvs2[cnt++];
1350:     for (j = 0; j < nt; j++) {
1351:       for (k = 0; k < recvs2[cnt + 1]; k++) (*indices)[bprocs[recvs2[cnt + 2 + k]]][(*numprocs)[bprocs[recvs2[cnt + 2 + k]]]++] = recvs2[cnt];
1352:       cnt += 2 + recvs2[cnt + 1];
1353:     }
1354:   }
1355:   PetscCall(PetscFree(bprocs));
1356:   PetscCall(PetscFree(recvs2));

1358:   /* sort the node indexing by their global numbers */
1359:   nt = *nproc;
1360:   for (i = 0; i < nt; i++) {
1361:     PetscCall(PetscMalloc1((*numprocs)[i], &tmp));
1362:     for (j = 0; j < (*numprocs)[i]; j++) tmp[j] = lindices[(*indices)[i][j]];
1363:     PetscCall(PetscSortIntWithArray((*numprocs)[i], tmp, (*indices)[i]));
1364:     PetscCall(PetscFree(tmp));
1365:   }

1367:   if (debug) { /* -----------------------------------  */
1368:     nt = *nproc;
1369:     for (i = 0; i < nt; i++) {
1370:       PetscCall(PetscSynchronizedPrintf(comm, "[%d] subdomain %" PetscInt_FMT " number of indices %" PetscInt_FMT ": ", rank, (*procs)[i], (*numprocs)[i]));
1371:       for (j = 0; j < (*numprocs)[i]; j++) PetscCall(PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", (*indices)[i][j]));
1372:       PetscCall(PetscSynchronizedPrintf(comm, "\n"));
1373:     }
1374:     PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));
1375:   } /* -----------------------------------  */

1377:   /* wait on sends */
1378:   if (nsends2) {
1379:     PetscCall(PetscMalloc1(nsends2, &send_status));
1380:     PetscCallMPI(MPI_Waitall(nsends2, send_waits, send_status));
1381:     PetscCall(PetscFree(send_status));
1382:   }

1384:   PetscCall(PetscFree(starts3));
1385:   PetscCall(PetscFree(dest));
1386:   PetscCall(PetscFree(send_waits));

1388:   PetscCall(PetscFree(nownedsenders));
1389:   PetscCall(PetscFree(ownedsenders));
1390:   PetscCall(PetscFree(starts));
1391:   PetscCall(PetscFree(starts2));
1392:   PetscCall(PetscFree(lens2));

1394:   PetscCall(PetscFree(source));
1395:   PetscCall(PetscFree(len));
1396:   PetscCall(PetscFree(recvs));
1397:   PetscCall(PetscFree(nprocs));
1398:   PetscCall(PetscFree(sends2));

1400:   /* put the information about myself as the first entry in the list */
1401:   first_procs    = (*procs)[0];
1402:   first_numprocs = (*numprocs)[0];
1403:   first_indices  = (*indices)[0];
1404:   for (i = 0; i < *nproc; i++) {
1405:     if ((*procs)[i] == rank) {
1406:       (*procs)[0]    = (*procs)[i];
1407:       (*numprocs)[0] = (*numprocs)[i];
1408:       (*indices)[0]  = (*indices)[i];
1409:       (*procs)[i]    = first_procs;
1410:       (*numprocs)[i] = first_numprocs;
1411:       (*indices)[i]  = first_indices;
1412:       break;
1413:     }
1414:   }

1416:   /* save info for reuse */
1417:   mapping->info_nproc    = *nproc;
1418:   mapping->info_procs    = *procs;
1419:   mapping->info_numprocs = *numprocs;
1420:   mapping->info_indices  = *indices;
1421:   mapping->info_cached   = PETSC_TRUE;
1422:   PetscFunctionReturn(PETSC_SUCCESS);
1423: }

1425: /*@C
1426:     ISLocalToGlobalMappingRestoreBlockInfo - Frees the memory allocated by `ISLocalToGlobalMappingGetBlockInfo()`

1428:     Collective

1430:     Input Parameter:
1431: .   mapping - the mapping from local to global indexing

1433:     Output Parameters:
1434: +   nproc - number of processors that are connected to this one
1435: .   proc - neighboring processors
1436: .   numproc - number of indices for each processor
1437: -   indices - indices of local nodes shared with neighbor (sorted by global numbering)

1439:     Level: advanced

1441: .seealso: [](sec_scatter), `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
1442:           `ISLocalToGlobalMappingGetInfo()`
1443: @*/
1444: PetscErrorCode ISLocalToGlobalMappingRestoreBlockInfo(ISLocalToGlobalMapping mapping, PetscInt *nproc, PetscInt *procs[], PetscInt *numprocs[], PetscInt **indices[])
1445: {
1446:   PetscFunctionBegin;
1448:   if (mapping->info_free) {
1449:     PetscCall(PetscFree(*numprocs));
1450:     if (*indices) {
1451:       PetscInt i;

1453:       PetscCall(PetscFree((*indices)[0]));
1454:       for (i = 1; i < *nproc; i++) PetscCall(PetscFree((*indices)[i]));
1455:       PetscCall(PetscFree(*indices));
1456:     }
1457:   }
1458:   *nproc    = 0;
1459:   *procs    = NULL;
1460:   *numprocs = NULL;
1461:   *indices  = NULL;
1462:   PetscFunctionReturn(PETSC_SUCCESS);
1463: }

1465: /*@C
1466:     ISLocalToGlobalMappingGetInfo - Gets the neighbor information for each processor and
1467:      each index shared by more than one processor

1469:     Collective

1471:     Input Parameter:
1472: .   mapping - the mapping from local to global indexing

1474:     Output Parameters:
1475: +   nproc - number of processors that are connected to this one
1476: .   proc - neighboring processors
1477: .   numproc - number of indices for each subdomain (processor)
1478: -   indices - indices of nodes (in local numbering) shared with neighbors (sorted by global numbering)

1480:     Level: advanced

1482:     Note:
1483:     The user needs to call `ISLocalToGlobalMappingRestoreInfo()` when the data is no longer needed.

1485:     Fortran Usage:
1486: .vb
1487:         PetscInt indices[nproc][numprocmax],ierr)
1488:         ISLocalToGlobalMpngGetInfoSize(ISLocalToGlobalMapping,PetscInt nproc,PetscInt numprocmax,ierr) followed by
1489:         ISLocalToGlobalMappingGetInfo(ISLocalToGlobalMapping,PetscInt nproc, PetscInt procs[nproc],PetscInt numprocs[nproc],
1490: .ve

1492:     Fortran Note:
1493:         There is no `ISLocalToGlobalMappingRestoreInfo()` in Fortran. You must make sure that `procs`[], `numprocs`[] and
1494:         `indices`[][] are large enough arrays, either by allocating them dynamically or defining static ones large enough.

1496: .seealso: [](sec_scatter), `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
1497:           `ISLocalToGlobalMappingRestoreInfo()`
1498: @*/
1499: PetscErrorCode ISLocalToGlobalMappingGetInfo(ISLocalToGlobalMapping mapping, PetscInt *nproc, PetscInt *procs[], PetscInt *numprocs[], PetscInt **indices[])
1500: {
1501:   PetscInt **bindices = NULL, *bnumprocs = NULL, bs, i, j, k;

1503:   PetscFunctionBegin;
1505:   bs = mapping->bs;
1506:   PetscCall(ISLocalToGlobalMappingGetBlockInfo(mapping, nproc, procs, &bnumprocs, &bindices));
1507:   if (bs > 1) { /* we need to expand the cached info */
1508:     PetscCall(PetscCalloc1(*nproc, &*indices));
1509:     PetscCall(PetscCalloc1(*nproc, &*numprocs));
1510:     for (i = 0; i < *nproc; i++) {
1511:       PetscCall(PetscMalloc1(bs * bnumprocs[i], &(*indices)[i]));
1512:       for (j = 0; j < bnumprocs[i]; j++) {
1513:         for (k = 0; k < bs; k++) (*indices)[i][j * bs + k] = bs * bindices[i][j] + k;
1514:       }
1515:       (*numprocs)[i] = bnumprocs[i] * bs;
1516:     }
1517:     mapping->info_free = PETSC_TRUE;
1518:   } else {
1519:     *numprocs = bnumprocs;
1520:     *indices  = bindices;
1521:   }
1522:   PetscFunctionReturn(PETSC_SUCCESS);
1523: }

1525: /*@C
1526:     ISLocalToGlobalMappingRestoreInfo - Frees the memory allocated by `ISLocalToGlobalMappingGetInfo()`

1528:     Collective

1530:     Input Parameter:
1531: .   mapping - the mapping from local to global indexing

1533:     Output Parameters:
1534: +   nproc - number of processors that are connected to this one
1535: .   proc - neighboring processors
1536: .   numproc - number of indices for each processor
1537: -   indices - indices of local nodes shared with neighbor (sorted by global numbering)

1539:     Level: advanced

1541: .seealso: [](sec_scatter), `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
1542:           `ISLocalToGlobalMappingGetInfo()`
1543: @*/
1544: PetscErrorCode ISLocalToGlobalMappingRestoreInfo(ISLocalToGlobalMapping mapping, PetscInt *nproc, PetscInt *procs[], PetscInt *numprocs[], PetscInt **indices[])
1545: {
1546:   PetscFunctionBegin;
1547:   PetscCall(ISLocalToGlobalMappingRestoreBlockInfo(mapping, nproc, procs, numprocs, indices));
1548:   PetscFunctionReturn(PETSC_SUCCESS);
1549: }

1551: /*@C
1552:     ISLocalToGlobalMappingGetNodeInfo - Gets the neighbor information for each MPI rank

1554:     Collective

1556:     Input Parameter:
1557: .   mapping - the mapping from local to global indexing

1559:     Output Parameters:
1560: +   nnodes - number of local nodes (same `ISLocalToGlobalMappingGetSize()`)
1561: .   count - number of neighboring processors per node
1562: -   indices - indices of processes sharing the node (sorted)

1564:     Level: advanced

1566:     Note:
1567:     The user needs to call `ISLocalToGlobalMappingRestoreInfo()` when the data is no longer needed.

1569: .seealso: [](sec_scatter), `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
1570:           `ISLocalToGlobalMappingGetInfo()`, `ISLocalToGlobalMappingRestoreNodeInfo()`
1571: @*/
1572: PetscErrorCode ISLocalToGlobalMappingGetNodeInfo(ISLocalToGlobalMapping mapping, PetscInt *nnodes, PetscInt *count[], PetscInt **indices[])
1573: {
1574:   PetscInt n;

1576:   PetscFunctionBegin;
1578:   PetscCall(ISLocalToGlobalMappingGetSize(mapping, &n));
1579:   if (!mapping->info_nodec) {
1580:     PetscInt i, m, n_neigh, *neigh, *n_shared, **shared;

1582:     PetscCall(PetscMalloc2(n + 1, &mapping->info_nodec, n, &mapping->info_nodei));
1583:     PetscCall(ISLocalToGlobalMappingGetInfo(mapping, &n_neigh, &neigh, &n_shared, &shared));
1584:     for (i = 0; i < n; i++) mapping->info_nodec[i] = 1;
1585:     m                      = n;
1586:     mapping->info_nodec[n] = 0;
1587:     for (i = 1; i < n_neigh; i++) {
1588:       PetscInt j;

1590:       m += n_shared[i];
1591:       for (j = 0; j < n_shared[i]; j++) mapping->info_nodec[shared[i][j]] += 1;
1592:     }
1593:     if (n) PetscCall(PetscMalloc1(m, &mapping->info_nodei[0]));
1594:     for (i = 1; i < n; i++) mapping->info_nodei[i] = mapping->info_nodei[i - 1] + mapping->info_nodec[i - 1];
1595:     PetscCall(PetscArrayzero(mapping->info_nodec, n));
1596:     for (i = 0; i < n; i++) {
1597:       mapping->info_nodec[i]    = 1;
1598:       mapping->info_nodei[i][0] = neigh[0];
1599:     }
1600:     for (i = 1; i < n_neigh; i++) {
1601:       PetscInt j;

1603:       for (j = 0; j < n_shared[i]; j++) {
1604:         PetscInt k = shared[i][j];

1606:         mapping->info_nodei[k][mapping->info_nodec[k]] = neigh[i];
1607:         mapping->info_nodec[k] += 1;
1608:       }
1609:     }
1610:     for (i = 0; i < n; i++) PetscCall(PetscSortRemoveDupsInt(&mapping->info_nodec[i], mapping->info_nodei[i]));
1611:     PetscCall(ISLocalToGlobalMappingRestoreInfo(mapping, &n_neigh, &neigh, &n_shared, &shared));
1612:   }
1613:   if (nnodes) *nnodes = n;
1614:   if (count) *count = mapping->info_nodec;
1615:   if (indices) *indices = mapping->info_nodei;
1616:   PetscFunctionReturn(PETSC_SUCCESS);
1617: }

1619: /*@C
1620:     ISLocalToGlobalMappingRestoreNodeInfo - Frees the memory allocated by `ISLocalToGlobalMappingGetNodeInfo()`

1622:     Collective

1624:     Input Parameter:
1625: .   mapping - the mapping from local to global indexing

1627:     Output Parameters:
1628: +   nnodes - number of local nodes
1629: .   count - number of neighboring processors per node
1630: -   indices - indices of processes sharing the node (sorted)

1632:     Level: advanced

1634: .seealso: [](sec_scatter), `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
1635:           `ISLocalToGlobalMappingGetInfo()`
1636: @*/
1637: PetscErrorCode ISLocalToGlobalMappingRestoreNodeInfo(ISLocalToGlobalMapping mapping, PetscInt *nnodes, PetscInt *count[], PetscInt **indices[])
1638: {
1639:   PetscFunctionBegin;
1641:   if (nnodes) *nnodes = 0;
1642:   if (count) *count = NULL;
1643:   if (indices) *indices = NULL;
1644:   PetscFunctionReturn(PETSC_SUCCESS);
1645: }

1647: /*MC
1648:     ISLocalToGlobalMappingGetIndicesF90 - Get global indices for every local point that is mapped

1650:     Synopsis:
1651:     ISLocalToGlobalMappingGetIndicesF90(ISLocalToGlobalMapping ltog, PetscInt, pointer :: array(:)}, integer ierr)

1653:     Not Collective

1655:     Input Parameter:
1656: .   A - the matrix

1658:     Output Parameter:
1659: .   array - array of indices, the length of this array may be obtained with `ISLocalToGlobalMappingGetSize()`

1661:     Level: advanced

1663:     Note:
1664:     Use  `ISLocalToGlobalMappingGetIndicesF90()` when you no longer need access to the data

1666: .seealso: [](sec_fortranarrays), [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingGetIndices()`, `ISLocalToGlobalMappingRestoreIndices()`,
1667:           `ISLocalToGlobalMappingRestoreIndicesF90()`
1668: M*/

1670: /*MC
1671:     ISLocalToGlobalMappingRestoreIndicesF90 - restores the global indices for every local point that is mapped obtained with `ISLocalToGlobalMappingGetIndicesF90()`

1673:     Synopsis:
1674:     ISLocalToGlobalMappingRestoreIndicesF90(ISLocalToGlobalMapping ltog, PetscInt, pointer :: array(:)}, integer ierr)

1676:     Not Collective

1678:     Input Parameters:
1679: +   A - the matrix
1680: -   array - array of indices, the length of this array may be obtained with `ISLocalToGlobalMappingGetSize()`

1682:     Level: advanced

1684: .seealso: [](sec_fortranarrays), [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingGetIndices()`, `ISLocalToGlobalMappingRestoreIndices()`,
1685:           `ISLocalToGlobalMappingGetIndicesF90()`
1686: M*/

1688: /*@C
1689:    ISLocalToGlobalMappingGetIndices - Get global indices for every local point that is mapped

1691:    Not Collective

1693:    Input Parameter:
1694: . ltog - local to global mapping

1696:    Output Parameter:
1697: . array - array of indices, the length of this array may be obtained with `ISLocalToGlobalMappingGetSize()`

1699:    Level: advanced

1701:    Note:
1702:     `ISLocalToGlobalMappingGetSize()` returns the length the this array

1704: .seealso: [](sec_scatter), `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingApply()`, `ISLocalToGlobalMappingRestoreIndices()`,
1705:           `ISLocalToGlobalMappingGetBlockIndices()`, `ISLocalToGlobalMappingRestoreBlockIndices()`
1706: @*/
1707: PetscErrorCode ISLocalToGlobalMappingGetIndices(ISLocalToGlobalMapping ltog, const PetscInt **array)
1708: {
1709:   PetscFunctionBegin;
1712:   if (ltog->bs == 1) {
1713:     *array = ltog->indices;
1714:   } else {
1715:     PetscInt       *jj, k, i, j, n = ltog->n, bs = ltog->bs;
1716:     const PetscInt *ii;

1718:     PetscCall(PetscMalloc1(bs * n, &jj));
1719:     *array = jj;
1720:     k      = 0;
1721:     ii     = ltog->indices;
1722:     for (i = 0; i < n; i++)
1723:       for (j = 0; j < bs; j++) jj[k++] = bs * ii[i] + j;
1724:   }
1725:   PetscFunctionReturn(PETSC_SUCCESS);
1726: }

1728: /*@C
1729:    ISLocalToGlobalMappingRestoreIndices - Restore indices obtained with `ISLocalToGlobalMappingGetIndices()`

1731:    Not Collective

1733:    Input Parameters:
1734: + ltog - local to global mapping
1735: - array - array of indices

1737:    Level: advanced

1739: .seealso: [](sec_scatter), `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingApply()`, `ISLocalToGlobalMappingGetIndices()`
1740: @*/
1741: PetscErrorCode ISLocalToGlobalMappingRestoreIndices(ISLocalToGlobalMapping ltog, const PetscInt **array)
1742: {
1743:   PetscFunctionBegin;
1746:   PetscCheck(ltog->bs != 1 || *array == ltog->indices, PETSC_COMM_SELF, PETSC_ERR_ARG_BADPTR, "Trying to return mismatched pointer");

1748:   if (ltog->bs > 1) PetscCall(PetscFree(*(void **)array));
1749:   PetscFunctionReturn(PETSC_SUCCESS);
1750: }

1752: /*@C
1753:    ISLocalToGlobalMappingGetBlockIndices - Get global indices for every local block

1755:    Not Collective

1757:    Input Parameter:
1758: . ltog - local to global mapping

1760:    Output Parameter:
1761: . array - array of indices

1763:    Level: advanced

1765: .seealso: [](sec_scatter), `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingApply()`, `ISLocalToGlobalMappingRestoreBlockIndices()`
1766: @*/
1767: PetscErrorCode ISLocalToGlobalMappingGetBlockIndices(ISLocalToGlobalMapping ltog, const PetscInt **array)
1768: {
1769:   PetscFunctionBegin;
1772:   *array = ltog->indices;
1773:   PetscFunctionReturn(PETSC_SUCCESS);
1774: }

1776: /*@C
1777:    ISLocalToGlobalMappingRestoreBlockIndices - Restore indices obtained with `ISLocalToGlobalMappingGetBlockIndices()`

1779:    Not Collective

1781:    Input Parameters:
1782: + ltog - local to global mapping
1783: - array - array of indices

1785:    Level: advanced

1787: .seealso: [](sec_scatter), `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingApply()`, `ISLocalToGlobalMappingGetIndices()`
1788: @*/
1789: PetscErrorCode ISLocalToGlobalMappingRestoreBlockIndices(ISLocalToGlobalMapping ltog, const PetscInt **array)
1790: {
1791:   PetscFunctionBegin;
1794:   PetscCheck(*array == ltog->indices, PETSC_COMM_SELF, PETSC_ERR_ARG_BADPTR, "Trying to return mismatched pointer");
1795:   *array = NULL;
1796:   PetscFunctionReturn(PETSC_SUCCESS);
1797: }

1799: /*@C
1800:    ISLocalToGlobalMappingConcatenate - Create a new mapping that concatenates a list of mappings

1802:    Not Collective

1804:    Input Parameters:
1805: + comm - communicator for the new mapping, must contain the communicator of every mapping to concatenate
1806: . n - number of mappings to concatenate
1807: - ltogs - local to global mappings

1809:    Output Parameter:
1810: . ltogcat - new mapping

1812:    Level: advanced

1814:    Note:
1815:    This currently always returns a mapping with block size of 1

1817:    Developer Note:
1818:    If all the input mapping have the same block size we could easily handle that as a special case

1820: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingCreate()`
1821: @*/
1822: PetscErrorCode ISLocalToGlobalMappingConcatenate(MPI_Comm comm, PetscInt n, const ISLocalToGlobalMapping ltogs[], ISLocalToGlobalMapping *ltogcat)
1823: {
1824:   PetscInt i, cnt, m, *idx;

1826:   PetscFunctionBegin;
1827:   PetscCheck(n >= 0, comm, PETSC_ERR_ARG_OUTOFRANGE, "Must have a non-negative number of mappings, given %" PetscInt_FMT, n);
1831:   for (cnt = 0, i = 0; i < n; i++) {
1832:     PetscCall(ISLocalToGlobalMappingGetSize(ltogs[i], &m));
1833:     cnt += m;
1834:   }
1835:   PetscCall(PetscMalloc1(cnt, &idx));
1836:   for (cnt = 0, i = 0; i < n; i++) {
1837:     const PetscInt *subidx;
1838:     PetscCall(ISLocalToGlobalMappingGetSize(ltogs[i], &m));
1839:     PetscCall(ISLocalToGlobalMappingGetIndices(ltogs[i], &subidx));
1840:     PetscCall(PetscArraycpy(&idx[cnt], subidx, m));
1841:     PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogs[i], &subidx));
1842:     cnt += m;
1843:   }
1844:   PetscCall(ISLocalToGlobalMappingCreate(comm, 1, cnt, idx, PETSC_OWN_POINTER, ltogcat));
1845:   PetscFunctionReturn(PETSC_SUCCESS);
1846: }

1848: /*MC
1849:       ISLOCALTOGLOBALMAPPINGBASIC - basic implementation of the `ISLocalToGlobalMapping` object. When `ISGlobalToLocalMappingApply()` is
1850:                                     used this is good for only small and moderate size problems.

1852:    Options Database Key:
1853: .   -islocaltoglobalmapping_type basic - select this method

1855:    Level: beginner

1857:    Developer Note:
1858:    This stores all the mapping information on each MPI rank.

1860: .seealso: [](sec_scatter), `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingSetType()`, `ISLOCALTOGLOBALMAPPINGHASH`
1861: M*/
1862: PETSC_EXTERN PetscErrorCode ISLocalToGlobalMappingCreate_Basic(ISLocalToGlobalMapping ltog)
1863: {
1864:   PetscFunctionBegin;
1865:   ltog->ops->globaltolocalmappingapply      = ISGlobalToLocalMappingApply_Basic;
1866:   ltog->ops->globaltolocalmappingsetup      = ISGlobalToLocalMappingSetUp_Basic;
1867:   ltog->ops->globaltolocalmappingapplyblock = ISGlobalToLocalMappingApplyBlock_Basic;
1868:   ltog->ops->destroy                        = ISLocalToGlobalMappingDestroy_Basic;
1869:   PetscFunctionReturn(PETSC_SUCCESS);
1870: }

1872: /*MC
1873:       ISLOCALTOGLOBALMAPPINGHASH - hash implementation of the `ISLocalToGlobalMapping` object. When `ISGlobalToLocalMappingApply()` is
1874:                                     used this is good for large memory problems.

1876:    Options Database Key:
1877: .   -islocaltoglobalmapping_type hash - select this method

1879:    Level: beginner

1881:    Note:
1882:     This is selected automatically for large problems if the user does not set the type.

1884: .seealso: [](sec_scatter), `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingSetType()`, `ISLOCALTOGLOBALMAPPINGBASIC`
1885: M*/
1886: PETSC_EXTERN PetscErrorCode ISLocalToGlobalMappingCreate_Hash(ISLocalToGlobalMapping ltog)
1887: {
1888:   PetscFunctionBegin;
1889:   ltog->ops->globaltolocalmappingapply      = ISGlobalToLocalMappingApply_Hash;
1890:   ltog->ops->globaltolocalmappingsetup      = ISGlobalToLocalMappingSetUp_Hash;
1891:   ltog->ops->globaltolocalmappingapplyblock = ISGlobalToLocalMappingApplyBlock_Hash;
1892:   ltog->ops->destroy                        = ISLocalToGlobalMappingDestroy_Hash;
1893:   PetscFunctionReturn(PETSC_SUCCESS);
1894: }

1896: /*@C
1897:     ISLocalToGlobalMappingRegister -  Registers a method for applying a global to local mapping with an `ISLocalToGlobalMapping`

1899:    Not Collective

1901:    Input Parameters:
1902: +  sname - name of a new method
1903: -  function - routine to create method context

1905:    Sample usage:
1906: .vb
1907:    ISLocalToGlobalMappingRegister("my_mapper", MyCreate);
1908: .ve

1910:    Then, your mapping can be chosen with the procedural interface via
1911: $     ISLocalToGlobalMappingSetType(ltog, "my_mapper")
1912:    or at runtime via the option
1913: $     -islocaltoglobalmapping_type my_mapper

1915:    Level: advanced

1917:    Note:
1918:    `ISLocalToGlobalMappingRegister()` may be called multiple times to add several user-defined mappings.

1920: .seealso: [](sec_scatter), `ISLocalToGlobalMappingRegisterAll()`, `ISLocalToGlobalMappingRegisterDestroy()`, `ISLOCALTOGLOBALMAPPINGBASIC`,
1921:           `ISLOCALTOGLOBALMAPPINGHASH`, `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingApply()`
1922: @*/
1923: PetscErrorCode ISLocalToGlobalMappingRegister(const char sname[], PetscErrorCode (*function)(ISLocalToGlobalMapping))
1924: {
1925:   PetscFunctionBegin;
1926:   PetscCall(ISInitializePackage());
1927:   PetscCall(PetscFunctionListAdd(&ISLocalToGlobalMappingList, sname, function));
1928:   PetscFunctionReturn(PETSC_SUCCESS);
1929: }

1931: /*@C
1932:    ISLocalToGlobalMappingSetType - Sets the implementation type `ISLocalToGlobalMapping` will use

1934:    Logically Collective

1936:    Input Parameters:
1937: +  ltog - the `ISLocalToGlobalMapping` object
1938: -  type - a known method

1940:    Options Database Key:
1941: .  -islocaltoglobalmapping_type  <method> - Sets the method; use -help for a list of available methods (for instance, basic or hash)

1943:   Level: intermediate

1945:    Notes:
1946:    See `ISLocalToGlobalMappingType` for available methods

1948:   Normally, it is best to use the `ISLocalToGlobalMappingSetFromOptions()` command and
1949:   then set the `ISLocalToGlobalMappingType` from the options database rather than by using
1950:   this routine.

1952:   Developer Note:
1953:   `ISLocalToGlobalMappingRegister()` is used to add new types to `ISLocalToGlobalMappingList` from which they
1954:   are accessed by `ISLocalToGlobalMappingSetType()`.

1956: .seealso: [](sec_scatter), `ISLocalToGlobalMappingType`, `ISLocalToGlobalMappingType`, `ISLocalToGlobalMappingRegister()`, `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingGetType()`
1957: @*/
1958: PetscErrorCode ISLocalToGlobalMappingSetType(ISLocalToGlobalMapping ltog, ISLocalToGlobalMappingType type)
1959: {
1960:   PetscBool match;
1961:   PetscErrorCode (*r)(ISLocalToGlobalMapping) = NULL;

1963:   PetscFunctionBegin;

1967:   PetscCall(PetscObjectTypeCompare((PetscObject)ltog, type, &match));
1968:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

1970:   /* L2G maps defer type setup at globaltolocal calls, allow passing NULL here */
1971:   if (type) {
1972:     PetscCall(PetscFunctionListFind(ISLocalToGlobalMappingList, type, &r));
1973:     PetscCheck(r, PetscObjectComm((PetscObject)ltog), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unable to find requested ISLocalToGlobalMapping type %s", type);
1974:   }
1975:   /* Destroy the previous private LTOG context */
1976:   PetscTryTypeMethod(ltog, destroy);
1977:   ltog->ops->destroy = NULL;

1979:   PetscCall(PetscObjectChangeTypeName((PetscObject)ltog, type));
1980:   if (r) PetscCall((*r)(ltog));
1981:   PetscFunctionReturn(PETSC_SUCCESS);
1982: }

1984: /*@C
1985:    ISLocalToGlobalMappingGetType - Get the type of the `ISLocalToGlobalMapping`

1987:    Not Collective

1989:    Input Parameter:
1990: .  ltog - the `ISLocalToGlobalMapping` object

1992:    Output Parameter:
1993: .  type - the type

1995:    Level: intermediate

1997: .seealso: [](sec_scatter), `ISLocalToGlobalMappingType`, `ISLocalToGlobalMappingType`, `ISLocalToGlobalMappingRegister()`, `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingSetType()`
1998: @*/
1999: PetscErrorCode ISLocalToGlobalMappingGetType(ISLocalToGlobalMapping ltog, ISLocalToGlobalMappingType *type)
2000: {
2001:   PetscFunctionBegin;
2004:   *type = ((PetscObject)ltog)->type_name;
2005:   PetscFunctionReturn(PETSC_SUCCESS);
2006: }

2008: PetscBool ISLocalToGlobalMappingRegisterAllCalled = PETSC_FALSE;

2010: /*@C
2011:   ISLocalToGlobalMappingRegisterAll - Registers all of the local to global mapping components in the `IS` package.

2013:   Not Collective

2015:   Level: advanced

2017: .seealso: [](sec_scatter), `ISRegister()`, `ISLocalToGlobalRegister()`
2018: @*/
2019: PetscErrorCode ISLocalToGlobalMappingRegisterAll(void)
2020: {
2021:   PetscFunctionBegin;
2022:   if (ISLocalToGlobalMappingRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
2023:   ISLocalToGlobalMappingRegisterAllCalled = PETSC_TRUE;
2024:   PetscCall(ISLocalToGlobalMappingRegister(ISLOCALTOGLOBALMAPPINGBASIC, ISLocalToGlobalMappingCreate_Basic));
2025:   PetscCall(ISLocalToGlobalMappingRegister(ISLOCALTOGLOBALMAPPINGHASH, ISLocalToGlobalMappingCreate_Hash));
2026:   PetscFunctionReturn(PETSC_SUCCESS);
2027: }