Actual source code: isltog.c

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

  7: PetscClassId          IS_LTOGM_CLASSID;
  8: static PetscErrorCode ISLocalToGlobalMappingSetUpBlockInfo_Private(ISLocalToGlobalMapping);

 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: /*@
 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_INT_MAX;

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: static PetscErrorCode ISLocalToGlobalMappingResetBlockInfo_Private(ISLocalToGlobalMapping mapping)
222: {
223:   PetscFunctionBegin;
224:   PetscCall(PetscFree(mapping->info_procs));
225:   PetscCall(PetscFree(mapping->info_numprocs));
226:   if (mapping->info_indices) {
227:     for (PetscInt i = 0; i < mapping->info_nproc; i++) PetscCall(PetscFree(mapping->info_indices[i]));
228:     PetscCall(PetscFree(mapping->info_indices));
229:   }
230:   if (mapping->info_nodei) PetscCall(PetscFree(mapping->info_nodei[0]));
231:   PetscCall(PetscFree2(mapping->info_nodec, mapping->info_nodei));
232:   PetscCall(PetscSFDestroy(&mapping->multileaves_sf));
233:   PetscFunctionReturn(PETSC_SUCCESS);
234: }

236: #define GTOLTYPE _Basic
237: #define GTOLNAME _Basic
238: #define GTOLBS   mapping->bs
239: #define GTOL(g, local) \
240:   do { \
241:     local = map->globals[g / bs - start]; \
242:     if (local >= 0) local = bs * local + (g % bs); \
243:   } while (0)

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

247: #define GTOLTYPE _Basic
248: #define GTOLNAME Block_Basic
249: #define GTOLBS   1
250: #define GTOL(g, local) \
251:   do { \
252:     local = map->globals[g - start]; \
253:   } while (0)
254: #include <../src/vec/is/utils/isltog.h>

256: #define GTOLTYPE _Hash
257: #define GTOLNAME _Hash
258: #define GTOLBS   mapping->bs
259: #define GTOL(g, local) \
260:   do { \
261:     (void)PetscHMapIGet(map->globalht, g / bs, &local); \
262:     if (local >= 0) local = bs * local + (g % bs); \
263:   } while (0)
264: #include <../src/vec/is/utils/isltog.h>

266: #define GTOLTYPE _Hash
267: #define GTOLNAME Block_Hash
268: #define GTOLBS   1
269: #define GTOL(g, local) \
270:   do { \
271:     (void)PetscHMapIGet(map->globalht, g, &local); \
272:   } while (0)
273: #include <../src/vec/is/utils/isltog.h>

275: /*@
276:   ISLocalToGlobalMappingDuplicate - Duplicates the local to global mapping object

278:   Not Collective

280:   Input Parameter:
281: . ltog - local to global mapping

283:   Output Parameter:
284: . nltog - the duplicated local to global mapping

286:   Level: advanced

288: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreate()`
289: @*/
290: PetscErrorCode ISLocalToGlobalMappingDuplicate(ISLocalToGlobalMapping ltog, ISLocalToGlobalMapping *nltog)
291: {
292:   ISLocalToGlobalMappingType l2gtype;

294:   PetscFunctionBegin;
296:   PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)ltog), ltog->bs, ltog->n, ltog->indices, PETSC_COPY_VALUES, nltog));
297:   PetscCall(ISLocalToGlobalMappingGetType(ltog, &l2gtype));
298:   PetscCall(ISLocalToGlobalMappingSetType(*nltog, l2gtype));
299:   PetscFunctionReturn(PETSC_SUCCESS);
300: }

302: /*@
303:   ISLocalToGlobalMappingGetSize - Gets the local size of a local to global mapping

305:   Not Collective

307:   Input Parameter:
308: . mapping - local to global mapping

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

313:   Level: advanced

315: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreate()`
316: @*/
317: PetscErrorCode ISLocalToGlobalMappingGetSize(ISLocalToGlobalMapping mapping, PetscInt *n)
318: {
319:   PetscFunctionBegin;
321:   PetscAssertPointer(n, 2);
322:   *n = mapping->bs * mapping->n;
323:   PetscFunctionReturn(PETSC_SUCCESS);
324: }

326: /*@
327:   ISLocalToGlobalMappingViewFromOptions - View an `ISLocalToGlobalMapping` based on values in the options database

329:   Collective

331:   Input Parameters:
332: + A    - the local to global mapping object
333: . obj  - Optional object that provides the options prefix used for the options database query
334: - name - command line option

336:   Level: intermediate

338:   Note:
339:   See `PetscObjectViewFromOptions()` for the available `PetscViewer` and `PetscViewerFormat`

341: .seealso: [](sec_scatter), `PetscViewer`, `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingView`, `PetscObjectViewFromOptions()`, `ISLocalToGlobalMappingCreate()`
342: @*/
343: PetscErrorCode ISLocalToGlobalMappingViewFromOptions(ISLocalToGlobalMapping A, PetscObject obj, const char name[])
344: {
345:   PetscFunctionBegin;
347:   PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
348:   PetscFunctionReturn(PETSC_SUCCESS);
349: }

351: /*@
352:   ISLocalToGlobalMappingView - View a local to global mapping

354:   Collective on viewer

356:   Input Parameters:
357: + mapping - local to global mapping
358: - viewer  - viewer

360:   Level: intermediate

362: .seealso: [](sec_scatter), `PetscViewer`, `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreate()`
363: @*/
364: PetscErrorCode ISLocalToGlobalMappingView(ISLocalToGlobalMapping mapping, PetscViewer viewer)
365: {
366:   PetscBool         iascii, isbinary;
367:   PetscViewerFormat format;

369:   PetscFunctionBegin;
371:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mapping), &viewer));

374:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
375:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
376:   PetscCall(PetscViewerGetFormat(viewer, &format));
377:   if (iascii) {
378:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
379:       const PetscInt *idxs;
380:       IS              is;
381:       const char     *name = ((PetscObject)mapping)->name;
382:       char            iname[PETSC_MAX_PATH_LEN];

384:       PetscCall(PetscSNPrintf(iname, sizeof(iname), "%sl2g", name ? name : ""));
385:       PetscCall(ISLocalToGlobalMappingGetIndices(mapping, &idxs));
386:       PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)viewer), mapping->n * mapping->bs, idxs, PETSC_USE_POINTER, &is));
387:       PetscCall(PetscObjectSetName((PetscObject)is, iname));
388:       PetscCall(ISView(is, viewer));
389:       PetscCall(ISLocalToGlobalMappingRestoreIndices(mapping, &idxs));
390:       PetscCall(ISDestroy(&is));
391:     } else {
392:       PetscMPIInt rank;

394:       PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)mapping), &rank));
395:       PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)mapping, viewer));
396:       PetscCall(PetscViewerASCIIPushSynchronized(viewer));
397:       for (PetscInt i = 0; i < mapping->n; i++) {
398:         PetscInt bs = mapping->bs, g = mapping->indices[i];
399:         if (bs == 1) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] %" PetscInt_FMT " %" PetscInt_FMT "\n", rank, i, g));
400:         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));
401:       }
402:       PetscCall(PetscViewerFlush(viewer));
403:       PetscCall(PetscViewerASCIIPopSynchronized(viewer));
404:     }
405:   } else if (isbinary) {
406:     PetscBool skipHeader;

408:     PetscCall(PetscViewerSetUp(viewer));
409:     PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
410:     if (!skipHeader) {
411:       PetscMPIInt size;
412:       PetscInt    tr[3];

414:       PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)viewer), &size));
415:       tr[0] = IS_LTOGM_FILE_CLASSID;
416:       tr[1] = mapping->bs;
417:       tr[2] = size;
418:       PetscCall(PetscViewerBinaryWrite(viewer, tr, 3, PETSC_INT));
419:       PetscCall(PetscViewerBinaryWriteAll(viewer, &mapping->n, 1, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
420:     }
421:     /* write block indices */
422:     PetscCall(PetscViewerBinaryWriteAll(viewer, mapping->indices, mapping->n, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
423:   }
424:   PetscFunctionReturn(PETSC_SUCCESS);
425: }

427: /*@
428:   ISLocalToGlobalMappingLoad - Loads a local-to-global mapping that has been stored in binary format.

430:   Collective on viewer

432:   Input Parameters:
433: + mapping - the newly loaded map, this needs to have been created with `ISLocalToGlobalMappingCreate()` or some related function before a call to `ISLocalToGlobalMappingLoad()`
434: - viewer  - binary file viewer, obtained from `PetscViewerBinaryOpen()`

436:   Level: intermediate

438: .seealso: [](sec_scatter), `PetscViewer`, `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingView()`, `ISLocalToGlobalMappingCreate()`
439: @*/
440: PetscErrorCode ISLocalToGlobalMappingLoad(ISLocalToGlobalMapping mapping, PetscViewer viewer)
441: {
442:   PetscBool isbinary, skipHeader;

444:   PetscFunctionBegin;
447:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
448:   PetscCheck(isbinary, PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Invalid viewer of type %s", ((PetscObject)viewer)->type_name);

450:   /* reset previous data */
451:   PetscCall(ISLocalToGlobalMappingResetBlockInfo_Private(mapping));

453:   PetscCall(PetscViewerSetUp(viewer));
454:   PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));

456:   /* When skipping header, it assumes bs and n have been already set */
457:   if (!skipHeader) {
458:     MPI_Comm comm = PetscObjectComm((PetscObject)viewer);
459:     PetscInt tr[3], nold = mapping->n, *sizes, nmaps = PETSC_DECIDE, st = 0;

461:     PetscCall(PetscViewerBinaryRead(viewer, tr, 3, NULL, PETSC_INT));
462:     PetscCheck(tr[0] == IS_LTOGM_FILE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a local-to-global map next in file");

464:     mapping->bs = tr[1];
465:     PetscCall(PetscMalloc1(tr[2], &sizes));
466:     PetscCall(PetscViewerBinaryRead(viewer, sizes, tr[2], NULL, PETSC_INT));

468:     /* consume the input, read multiple maps per process if needed */
469:     PetscCall(PetscSplitOwnership(comm, &nmaps, &tr[2]));
470:     PetscCallMPI(MPI_Exscan(&nmaps, &st, 1, MPIU_INT, MPI_SUM, comm));
471:     mapping->n = 0;
472:     for (PetscInt i = st; i < st + nmaps; i++) mapping->n += sizes[i];
473:     PetscCall(PetscFree(sizes));

475:     if (nold != mapping->n) {
476:       if (mapping->dealloc_indices) PetscCall(PetscFree(mapping->indices));
477:       mapping->indices = NULL;
478:     }
479:   }

481:   /* read indices */
482:   if (mapping->n && !mapping->indices) {
483:     PetscCall(PetscMalloc1(mapping->n, &mapping->indices));
484:     mapping->dealloc_indices = PETSC_TRUE;
485:   }
486:   PetscCall(PetscViewerBinaryReadAll(viewer, mapping->indices, mapping->n, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
487:   PetscFunctionReturn(PETSC_SUCCESS);
488: }

490: /*@
491:   ISLocalToGlobalMappingCreateIS - Creates a mapping between a local (0 to n)
492:   ordering and a global parallel ordering.

494:   Not Collective

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

499:   Output Parameter:
500: . mapping - new mapping data structure

502:   Level: advanced

504:   Note:
505:   the block size of the `IS` determines the block size of the mapping

507: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingSetFromOptions()`
508: @*/
509: PetscErrorCode ISLocalToGlobalMappingCreateIS(IS is, ISLocalToGlobalMapping *mapping)
510: {
511:   PetscInt        n, bs;
512:   const PetscInt *indices;
513:   MPI_Comm        comm;
514:   PetscBool       isblock;

516:   PetscFunctionBegin;
518:   PetscAssertPointer(mapping, 2);

520:   PetscCall(PetscObjectGetComm((PetscObject)is, &comm));
521:   PetscCall(ISGetLocalSize(is, &n));
522:   PetscCall(PetscObjectTypeCompare((PetscObject)is, ISBLOCK, &isblock));
523:   if (!isblock) {
524:     PetscCall(ISGetIndices(is, &indices));
525:     PetscCall(ISLocalToGlobalMappingCreate(comm, 1, n, indices, PETSC_COPY_VALUES, mapping));
526:     PetscCall(ISRestoreIndices(is, &indices));
527:   } else {
528:     PetscCall(ISGetBlockSize(is, &bs));
529:     PetscCall(ISBlockGetIndices(is, &indices));
530:     PetscCall(ISLocalToGlobalMappingCreate(comm, bs, n / bs, indices, PETSC_COPY_VALUES, mapping));
531:     PetscCall(ISBlockRestoreIndices(is, &indices));
532:   }
533:   PetscFunctionReturn(PETSC_SUCCESS);
534: }

536: /*@
537:   ISLocalToGlobalMappingCreateSF - Creates a mapping between a local (0 to n) ordering and a global parallel ordering induced by a star forest.

539:   Collective

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

545:   Output Parameter:
546: . mapping - new mapping data structure

548:   Level: advanced

550:   Note:
551:   If a process calls this function with `start` = `PETSC_DECIDE` then all processes must, otherwise the program will hang.

553: .seealso: [](sec_scatter), `PetscSF`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingSetFromOptions()`
554: @*/
555: PetscErrorCode ISLocalToGlobalMappingCreateSF(PetscSF sf, PetscInt start, ISLocalToGlobalMapping *mapping)
556: {
557:   PetscInt i, maxlocal, nroots, nleaves, *globals, *ltog;
558:   MPI_Comm comm;

560:   PetscFunctionBegin;
562:   PetscAssertPointer(mapping, 3);
563:   PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
564:   PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL));
565:   if (start == PETSC_DECIDE) {
566:     start = 0;
567:     PetscCallMPI(MPI_Exscan(&nroots, &start, 1, MPIU_INT, MPI_SUM, comm));
568:   } else PetscCheck(start >= 0, comm, PETSC_ERR_ARG_OUTOFRANGE, "start must be nonnegative or PETSC_DECIDE");
569:   PetscCall(PetscSFGetLeafRange(sf, NULL, &maxlocal));
570:   ++maxlocal;
571:   PetscCall(PetscMalloc1(nroots, &globals));
572:   PetscCall(PetscMalloc1(maxlocal, &ltog));
573:   for (i = 0; i < nroots; i++) globals[i] = start + i;
574:   for (i = 0; i < maxlocal; i++) ltog[i] = -1;
575:   PetscCall(PetscSFBcastBegin(sf, MPIU_INT, globals, ltog, MPI_REPLACE));
576:   PetscCall(PetscSFBcastEnd(sf, MPIU_INT, globals, ltog, MPI_REPLACE));
577:   PetscCall(ISLocalToGlobalMappingCreate(comm, 1, maxlocal, ltog, PETSC_OWN_POINTER, mapping));
578:   PetscCall(PetscFree(globals));
579:   PetscFunctionReturn(PETSC_SUCCESS);
580: }

582: /*@
583:   ISLocalToGlobalMappingSetBlockSize - Sets the blocksize of the mapping

585:   Not Collective

587:   Input Parameters:
588: + mapping - mapping data structure
589: - bs      - the blocksize

591:   Level: advanced

593: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`
594: @*/
595: PetscErrorCode ISLocalToGlobalMappingSetBlockSize(ISLocalToGlobalMapping mapping, PetscInt bs)
596: {
597:   PetscInt       *nid;
598:   const PetscInt *oid;
599:   PetscInt        i, cn, on, obs, nn;

601:   PetscFunctionBegin;
603:   PetscCheck(bs >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid block size %" PetscInt_FMT, bs);
604:   if (bs == mapping->bs) PetscFunctionReturn(PETSC_SUCCESS);
605:   on  = mapping->n;
606:   obs = mapping->bs;
607:   oid = mapping->indices;
608:   nn  = (on * obs) / bs;
609:   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);

611:   PetscCall(PetscMalloc1(nn, &nid));
612:   PetscCall(ISLocalToGlobalMappingGetIndices(mapping, &oid));
613:   for (i = 0; i < nn; i++) {
614:     PetscInt j;
615:     for (j = 0, cn = 0; j < bs - 1; j++) {
616:       if (oid[i * bs + j] < 0) {
617:         cn++;
618:         continue;
619:       }
620:       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]);
621:     }
622:     if (oid[i * bs + j] < 0) cn++;
623:     if (cn) {
624:       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);
625:       nid[i] = -1;
626:     } else {
627:       nid[i] = oid[i * bs] / bs;
628:     }
629:   }
630:   PetscCall(ISLocalToGlobalMappingRestoreIndices(mapping, &oid));

632:   mapping->n  = nn;
633:   mapping->bs = bs;
634:   PetscCall(PetscFree(mapping->indices));
635:   mapping->indices     = nid;
636:   mapping->globalstart = 0;
637:   mapping->globalend   = 0;

639:   /* reset the cached information */
640:   PetscCall(ISLocalToGlobalMappingResetBlockInfo_Private(mapping));
641:   PetscTryTypeMethod(mapping, destroy);
642:   PetscFunctionReturn(PETSC_SUCCESS);
643: }

645: /*@
646:   ISLocalToGlobalMappingGetBlockSize - Gets the blocksize of the mapping
647:   ordering and a global parallel ordering.

649:   Not Collective

651:   Input Parameter:
652: . mapping - mapping data structure

654:   Output Parameter:
655: . bs - the blocksize

657:   Level: advanced

659: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`
660: @*/
661: PetscErrorCode ISLocalToGlobalMappingGetBlockSize(ISLocalToGlobalMapping mapping, PetscInt *bs)
662: {
663:   PetscFunctionBegin;
665:   *bs = mapping->bs;
666:   PetscFunctionReturn(PETSC_SUCCESS);
667: }

669: /*@
670:   ISLocalToGlobalMappingCreate - Creates a mapping between a local (0 to n)
671:   ordering and a global parallel ordering.

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

675:   Input Parameters:
676: + comm    - MPI communicator
677: . bs      - the block size
678: . n       - the number of local elements divided by the block size, or equivalently the number of block indices
679: . 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
680: - mode    - see PetscCopyMode

682:   Output Parameter:
683: . mapping - new mapping data structure

685:   Level: advanced

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

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

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

697: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingSetFromOptions()`,
698:           `ISLOCALTOGLOBALMAPPINGBASIC`, `ISLOCALTOGLOBALMAPPINGHASH`
699:           `ISLocalToGlobalMappingSetType()`, `ISLocalToGlobalMappingType`
700: @*/
701: PetscErrorCode ISLocalToGlobalMappingCreate(MPI_Comm comm, PetscInt bs, PetscInt n, const PetscInt indices[], PetscCopyMode mode, ISLocalToGlobalMapping *mapping)
702: {
703:   PetscInt *in;

705:   PetscFunctionBegin;
706:   if (n) PetscAssertPointer(indices, 4);
707:   PetscAssertPointer(mapping, 6);

709:   *mapping = NULL;
710:   PetscCall(ISInitializePackage());

712:   PetscCall(PetscHeaderCreate(*mapping, IS_LTOGM_CLASSID, "ISLocalToGlobalMapping", "Local to global mapping", "IS", comm, ISLocalToGlobalMappingDestroy, ISLocalToGlobalMappingView));
713:   (*mapping)->n  = n;
714:   (*mapping)->bs = bs;
715:   if (mode == PETSC_COPY_VALUES) {
716:     PetscCall(PetscMalloc1(n, &in));
717:     PetscCall(PetscArraycpy(in, indices, n));
718:     (*mapping)->indices         = in;
719:     (*mapping)->dealloc_indices = PETSC_TRUE;
720:   } else if (mode == PETSC_OWN_POINTER) {
721:     (*mapping)->indices         = (PetscInt *)indices;
722:     (*mapping)->dealloc_indices = PETSC_TRUE;
723:   } else if (mode == PETSC_USE_POINTER) {
724:     (*mapping)->indices = (PetscInt *)indices;
725:   } else SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid mode %d", mode);
726:   PetscFunctionReturn(PETSC_SUCCESS);
727: }

729: PetscFunctionList ISLocalToGlobalMappingList = NULL;

731: /*@
732:   ISLocalToGlobalMappingSetFromOptions - Set mapping options from the options database.

734:   Not Collective

736:   Input Parameter:
737: . mapping - mapping data structure

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

742:   Level: advanced

744: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingDestroy()`,
745: `ISLocalToGlobalMappingCreateIS()`, `ISLOCALTOGLOBALMAPPINGBASIC`,
746: `ISLOCALTOGLOBALMAPPINGHASH`, `ISLocalToGlobalMappingSetType()`, `ISLocalToGlobalMappingType`
747: @*/
748: PetscErrorCode ISLocalToGlobalMappingSetFromOptions(ISLocalToGlobalMapping mapping)
749: {
750:   char                       type[256];
751:   ISLocalToGlobalMappingType defaulttype = "Not set";
752:   PetscBool                  flg;

754:   PetscFunctionBegin;
756:   PetscCall(ISLocalToGlobalMappingRegisterAll());
757:   PetscObjectOptionsBegin((PetscObject)mapping);
758:   PetscCall(PetscOptionsFList("-islocaltoglobalmapping_type", "ISLocalToGlobalMapping method", "ISLocalToGlobalMappingSetType", ISLocalToGlobalMappingList, ((PetscObject)mapping)->type_name ? ((PetscObject)mapping)->type_name : defaulttype, type, 256, &flg));
759:   if (flg) PetscCall(ISLocalToGlobalMappingSetType(mapping, type));
760:   PetscOptionsEnd();
761:   PetscFunctionReturn(PETSC_SUCCESS);
762: }

764: /*@
765:   ISLocalToGlobalMappingDestroy - Destroys a mapping between a local (0 to n)
766:   ordering and a global parallel ordering.

768:   Not Collective

770:   Input Parameter:
771: . mapping - mapping data structure

773:   Level: advanced

775: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingCreate()`
776: @*/
777: PetscErrorCode ISLocalToGlobalMappingDestroy(ISLocalToGlobalMapping *mapping)
778: {
779:   PetscFunctionBegin;
780:   if (!*mapping) PetscFunctionReturn(PETSC_SUCCESS);
782:   if (--((PetscObject)*mapping)->refct > 0) {
783:     *mapping = NULL;
784:     PetscFunctionReturn(PETSC_SUCCESS);
785:   }
786:   if ((*mapping)->dealloc_indices) PetscCall(PetscFree((*mapping)->indices));
787:   PetscCall(ISLocalToGlobalMappingResetBlockInfo_Private(*mapping));
788:   PetscTryTypeMethod(*mapping, destroy);
789:   PetscCall(PetscHeaderDestroy(mapping));
790:   PetscFunctionReturn(PETSC_SUCCESS);
791: }

793: /*@
794:   ISLocalToGlobalMappingApplyIS - Creates from an `IS` in the local numbering
795:   a new index set using the global numbering defined in an `ISLocalToGlobalMapping`
796:   context.

798:   Collective

800:   Input Parameters:
801: + mapping - mapping between local and global numbering
802: - is      - index set in local numbering

804:   Output Parameter:
805: . newis - index set in global numbering

807:   Level: advanced

809:   Note:
810:   The output `IS` will have the same communicator as the input `IS`.

812: .seealso: [](sec_scatter), `ISLocalToGlobalMappingApply()`, `ISLocalToGlobalMappingCreate()`,
813:           `ISLocalToGlobalMappingDestroy()`, `ISGlobalToLocalMappingApply()`
814: @*/
815: PetscErrorCode ISLocalToGlobalMappingApplyIS(ISLocalToGlobalMapping mapping, IS is, IS *newis)
816: {
817:   PetscInt        n, *idxout;
818:   const PetscInt *idxin;

820:   PetscFunctionBegin;
823:   PetscAssertPointer(newis, 3);

825:   PetscCall(ISGetLocalSize(is, &n));
826:   PetscCall(ISGetIndices(is, &idxin));
827:   PetscCall(PetscMalloc1(n, &idxout));
828:   PetscCall(ISLocalToGlobalMappingApply(mapping, n, idxin, idxout));
829:   PetscCall(ISRestoreIndices(is, &idxin));
830:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is), n, idxout, PETSC_OWN_POINTER, newis));
831:   PetscFunctionReturn(PETSC_SUCCESS);
832: }

834: /*@
835:   ISLocalToGlobalMappingApply - Takes a list of integers in a local numbering
836:   and converts them to the global numbering.

838:   Not Collective

840:   Input Parameters:
841: + mapping - the local to global mapping context
842: . N       - number of integers
843: - in      - input indices in local numbering

845:   Output Parameter:
846: . out - indices in global numbering

848:   Level: advanced

850:   Note:
851:   The `in` and `out` array parameters may be identical.

853: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingApplyBlock()`, `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingDestroy()`,
854:           `ISLocalToGlobalMappingApplyIS()`, `AOCreateBasic()`, `AOApplicationToPetsc()`,
855:           `AOPetscToApplication()`, `ISGlobalToLocalMappingApply()`
856: @*/
857: PetscErrorCode ISLocalToGlobalMappingApply(ISLocalToGlobalMapping mapping, PetscInt N, const PetscInt in[], PetscInt out[])
858: {
859:   PetscInt i, bs, Nmax;

861:   PetscFunctionBegin;
863:   bs   = mapping->bs;
864:   Nmax = bs * mapping->n;
865:   if (bs == 1) {
866:     const PetscInt *idx = mapping->indices;
867:     for (i = 0; i < N; i++) {
868:       if (in[i] < 0) {
869:         out[i] = in[i];
870:         continue;
871:       }
872:       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);
873:       out[i] = idx[in[i]];
874:     }
875:   } else {
876:     const PetscInt *idx = mapping->indices;
877:     for (i = 0; i < N; i++) {
878:       if (in[i] < 0) {
879:         out[i] = in[i];
880:         continue;
881:       }
882:       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);
883:       out[i] = idx[in[i] / bs] * bs + (in[i] % bs);
884:     }
885:   }
886:   PetscFunctionReturn(PETSC_SUCCESS);
887: }

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

892:   Not Collective

894:   Input Parameters:
895: + mapping - the local to global mapping context
896: . N       - number of integers
897: - in      - input indices in local block numbering

899:   Output Parameter:
900: . out - indices in global block numbering

902:   Example:
903:   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
904:   (the first block) would produce 0 and the mapping applied to 1 (the second block) would produce 3.

906:   Level: advanced

908:   Note:
909:   The `in` and `out` array parameters may be identical.

911: .seealso: [](sec_scatter), `ISLocalToGlobalMappingApply()`, `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingDestroy()`,
912:           `ISLocalToGlobalMappingApplyIS()`, `AOCreateBasic()`, `AOApplicationToPetsc()`,
913:           `AOPetscToApplication()`, `ISGlobalToLocalMappingApply()`
914: @*/
915: PetscErrorCode ISLocalToGlobalMappingApplyBlock(ISLocalToGlobalMapping mapping, PetscInt N, const PetscInt in[], PetscInt out[])
916: {
917:   PetscInt        i, Nmax;
918:   const PetscInt *idx;

920:   PetscFunctionBegin;
922:   Nmax = mapping->n;
923:   idx  = mapping->indices;
924:   for (i = 0; i < N; i++) {
925:     if (in[i] < 0) {
926:       out[i] = in[i];
927:       continue;
928:     }
929:     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);
930:     out[i] = idx[in[i]];
931:   }
932:   PetscFunctionReturn(PETSC_SUCCESS);
933: }

935: /*@
936:   ISGlobalToLocalMappingApply - Provides the local numbering for a list of integers
937:   specified with a global numbering.

939:   Not Collective

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

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

956:   Level: advanced

958:   Notes:
959:   Either `nout` or `idxout` may be `NULL`. `idx` and `idxout` may be identical.

961:   For "small" problems when using `ISGlobalToLocalMappingApply()` and `ISGlobalToLocalMappingApplyBlock()`, the `ISLocalToGlobalMappingType` of
962:   `ISLOCALTOGLOBALMAPPINGBASIC` will be used;
963:   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.
964:   Use `ISLocalToGlobalMappingSetType()` or call `ISLocalToGlobalMappingSetFromOptions()` with the option -islocaltoglobalmapping_type <basic,hash> to control which is used.

966:   Developer Notes:
967:   The manual page states that `idx` and `idxout` may be identical but the calling
968:   sequence declares `idx` as const so it cannot be the same as `idxout`.

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

982: /*@
983:   ISGlobalToLocalMappingApplyIS - Creates from an `IS` in the global numbering
984:   a new index set using the local numbering defined in an `ISLocalToGlobalMapping`
985:   context.

987:   Not Collective

989:   Input Parameters:
990: + mapping - mapping between local and global numbering
991: . type    - `IS_GTOLM_MASK` - maps global indices with no local value to -1 in the output list (i.e., mask them)
992:            `IS_GTOLM_DROP` - drops the indices with no local value from the output list
993: - is      - index set in global numbering

995:   Output Parameter:
996: . newis - index set in local numbering

998:   Level: advanced

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

1003: .seealso: [](sec_scatter), `ISGlobalToLocalMapping`, `ISGlobalToLocalMappingApply()`, `ISLocalToGlobalMappingCreate()`,
1004:           `ISLocalToGlobalMappingDestroy()`
1005: @*/
1006: PetscErrorCode ISGlobalToLocalMappingApplyIS(ISLocalToGlobalMapping mapping, ISGlobalToLocalMappingMode type, IS is, IS *newis)
1007: {
1008:   PetscInt        n, nout, *idxout;
1009:   const PetscInt *idxin;

1011:   PetscFunctionBegin;
1014:   PetscAssertPointer(newis, 4);

1016:   PetscCall(ISGetLocalSize(is, &n));
1017:   PetscCall(ISGetIndices(is, &idxin));
1018:   if (type == IS_GTOLM_MASK) {
1019:     PetscCall(PetscMalloc1(n, &idxout));
1020:   } else {
1021:     PetscCall(ISGlobalToLocalMappingApply(mapping, type, n, idxin, &nout, NULL));
1022:     PetscCall(PetscMalloc1(nout, &idxout));
1023:   }
1024:   PetscCall(ISGlobalToLocalMappingApply(mapping, type, n, idxin, &nout, idxout));
1025:   PetscCall(ISRestoreIndices(is, &idxin));
1026:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nout, idxout, PETSC_OWN_POINTER, newis));
1027:   PetscFunctionReturn(PETSC_SUCCESS);
1028: }

1030: /*@
1031:   ISGlobalToLocalMappingApplyBlock - Provides the local block numbering for a list of integers
1032:   specified with a block global numbering.

1034:   Not Collective

1036:   Input Parameters:
1037: + mapping - mapping between local and global numbering
1038: . type    - `IS_GTOLM_MASK` - maps global indices with no local value to -1 in the output list (i.e., mask them)
1039:            `IS_GTOLM_DROP` - drops the indices with no local value from the output list
1040: . n       - number of global indices to map
1041: - idx     - global indices to map

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

1051:   Level: advanced

1053:   Notes:
1054:   Either `nout` or `idxout` may be `NULL`. `idx` and `idxout` may be identical.

1056:   For "small" problems when using `ISGlobalToLocalMappingApply()` and `ISGlobalToLocalMappingApplyBlock()`, the `ISLocalToGlobalMappingType` of
1057:   `ISLOCALTOGLOBALMAPPINGBASIC` will be used;
1058:   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.
1059:   Use `ISLocalToGlobalMappingSetType()` or call `ISLocalToGlobalMappingSetFromOptions()` with the option -islocaltoglobalmapping_type <basic,hash> to control which is used.

1061:   Developer Notes:
1062:   The manual page states that `idx` and `idxout` may be identical but the calling
1063:   sequence declares `idx` as const so it cannot be the same as `idxout`.

1065: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingApply()`, `ISGlobalToLocalMappingApply()`, `ISLocalToGlobalMappingCreate()`,
1066:           `ISLocalToGlobalMappingDestroy()`
1067: @*/
1068: PetscErrorCode ISGlobalToLocalMappingApplyBlock(ISLocalToGlobalMapping mapping, ISGlobalToLocalMappingMode type, PetscInt n, const PetscInt idx[], PetscInt *nout, PetscInt idxout[])
1069: {
1070:   PetscFunctionBegin;
1072:   if (!mapping->data) PetscCall(ISGlobalToLocalMappingSetUp(mapping));
1073:   PetscUseTypeMethod(mapping, globaltolocalmappingapplyblock, type, n, idx, nout, idxout);
1074:   PetscFunctionReturn(PETSC_SUCCESS);
1075: }

1077: /*@C
1078:   ISLocalToGlobalMappingGetBlockInfo - Gets the neighbor information

1080:   Collective the first time it is called

1082:   Input Parameter:
1083: . mapping - the mapping from local to global indexing

1085:   Output Parameters:
1086: + nproc    - number of processes that are connected to the calling process
1087: . procs    - neighboring processes
1088: . numprocs - number of block indices for each process
1089: - indices  - block indices (in local numbering) shared with neighbors (sorted by global numbering)

1091:   Level: advanced

1093: .seealso: [](sec_scatter), `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
1094:           `ISLocalToGlobalMappingRestoreBlockInfo()`, `ISLocalToGlobalMappingGetBlockMultiLeavesSF()`
1095: @*/
1096: PetscErrorCode ISLocalToGlobalMappingGetBlockInfo(ISLocalToGlobalMapping mapping, PetscInt *nproc, PetscInt *procs[], PetscInt *numprocs[], PetscInt **indices[])
1097: {
1098:   PetscFunctionBegin;
1100:   PetscCall(ISLocalToGlobalMappingSetUpBlockInfo_Private(mapping));
1101:   if (nproc) *nproc = mapping->info_nproc;
1102:   if (procs) *procs = mapping->info_procs;
1103:   if (numprocs) *numprocs = mapping->info_numprocs;
1104:   if (indices) *indices = mapping->info_indices;
1105:   PetscFunctionReturn(PETSC_SUCCESS);
1106: }

1108: /*@C
1109:   ISLocalToGlobalMappingGetBlockNodeInfo - Gets the neighbor information for each local block index

1111:   Collective the first time it is called

1113:   Input Parameter:
1114: . mapping - the mapping from local to global indexing

1116:   Output Parameters:
1117: + n       - number of local block nodes
1118: . n_procs - an array storing the number of processes for each local block node (including self)
1119: - procs   - the processes' rank for each local block node (sorted, self is first)

1121:   Level: advanced

1123:   Notes:
1124:   The user needs to call `ISLocalToGlobalMappingRestoreBlockNodeInfo()` when the data is no longer needed.
1125:   The information returned by this function complements that of `ISLocalToGlobalMappingGetBlockInfo()`.
1126:   The latter only provides local information, and the neighboring information
1127:   cannot be inferred in the general case, unless the mapping is locally one-to-one on each process.

1129: .seealso: `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
1130:           `ISLocalToGlobalMappingGetBlockInfo()`, `ISLocalToGlobalMappingRestoreBlockNodeInfo()`, `ISLocalToGlobalMappingGetNodeInfo()`
1131: @*/
1132: PetscErrorCode ISLocalToGlobalMappingGetBlockNodeInfo(ISLocalToGlobalMapping mapping, PetscInt *n, PetscInt *n_procs[], PetscInt **procs[])
1133: {
1134:   PetscFunctionBegin;
1136:   PetscCall(ISLocalToGlobalMappingSetUpBlockInfo_Private(mapping));
1137:   if (n) *n = mapping->n;
1138:   if (n_procs) *n_procs = mapping->info_nodec;
1139:   if (procs) *procs = mapping->info_nodei;
1140:   PetscFunctionReturn(PETSC_SUCCESS);
1141: }

1143: /*@C
1144:   ISLocalToGlobalMappingRestoreBlockNodeInfo - Frees the memory allocated by `ISLocalToGlobalMappingGetBlockNodeInfo()`

1146:   Not Collective

1148:   Input Parameters:
1149: + mapping - the mapping from local to global indexing
1150: . n       - number of local block nodes
1151: . n_procs - an array storing the number of processes for each local block nodes (including self)
1152: - procs   - the processes' rank for each local block node (sorted, self is first)

1154:   Level: advanced

1156: .seealso: `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
1157:           `ISLocalToGlobalMappingGetBlockNodeInfo()`
1158: @*/
1159: PetscErrorCode ISLocalToGlobalMappingRestoreBlockNodeInfo(ISLocalToGlobalMapping mapping, PetscInt *n, PetscInt *n_procs[], PetscInt **procs[])
1160: {
1161:   PetscFunctionBegin;
1163:   if (n) *n = 0;
1164:   if (n_procs) *n_procs = NULL;
1165:   if (procs) *procs = NULL;
1166:   PetscFunctionReturn(PETSC_SUCCESS);
1167: }

1169: /*@C
1170:   ISLocalToGlobalMappingGetBlockMultiLeavesSF - Get the star-forest to communicate multi-leaf block data

1172:   Collective the first time it is called

1174:   Input Parameter:
1175: . mapping - the mapping from local to global indexing

1177:   Output Parameter:
1178: . mlsf - the `PetscSF`

1180:   Level: advanced

1182:   Notes:
1183:   The returned star forest is suitable to exchange local information with other processes sharing the same global block index.
1184:   For example, suppose a mapping with two processes has been created with
1185: .vb
1186:     rank 0 global block indices: [0, 1, 2]
1187:     rank 1 global block indices: [2, 3, 4]
1188: .ve
1189:   and we want to share the local information
1190: .vb
1191:     rank 0 data: [-1, -2, -3]
1192:     rank 1 data: [1, 2, 3]
1193: .ve
1194:   then, the broadcasting action of `mlsf` will allow to collect
1195: .vb
1196:     rank 0 mlleafdata: [-1, -2, -3, 3]
1197:     rank 1 mlleafdata: [-3, 3, 1, 2]
1198: .ve
1199:   Use ``ISLocalToGlobalMappingGetBlockNodeInfo()`` to index into the multi-leaf data.

1201: .seealso: [](sec_scatter), `ISLocalToGlobalMappingGetBlockNodeInfo()`, `PetscSF`
1202: @*/
1203: PetscErrorCode ISLocalToGlobalMappingGetBlockMultiLeavesSF(ISLocalToGlobalMapping mapping, PetscSF *mlsf)
1204: {
1205:   PetscFunctionBegin;
1207:   PetscAssertPointer(mlsf, 2);
1208:   PetscCall(ISLocalToGlobalMappingSetUpBlockInfo_Private(mapping));
1209:   *mlsf = mapping->multileaves_sf;
1210:   PetscFunctionReturn(PETSC_SUCCESS);
1211: }

1213: static PetscErrorCode ISLocalToGlobalMappingSetUpBlockInfo_Private(ISLocalToGlobalMapping mapping)
1214: {
1215:   PetscSF            sf, sf2, imsf, msf;
1216:   MPI_Comm           comm;
1217:   const PetscSFNode *sfnode;
1218:   PetscSFNode       *newsfnode;
1219:   PetscLayout        layout;
1220:   PetscHMapI         neighs;
1221:   PetscHashIter      iter;
1222:   PetscBool          missing;
1223:   const PetscInt    *gidxs, *rootdegree;
1224:   PetscInt          *mask, *mrootdata, *leafdata, *newleafdata, *leafrd, *tmpg;
1225:   PetscInt           nroots, nleaves, newnleaves, bs, i, j, m, mnroots, p;
1226:   PetscMPIInt        rank, size;

1228:   PetscFunctionBegin;
1229:   if (mapping->multileaves_sf) PetscFunctionReturn(PETSC_SUCCESS);
1230:   PetscCall(PetscObjectGetComm((PetscObject)mapping, &comm));
1231:   PetscCallMPI(MPI_Comm_size(comm, &size));
1232:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

1234:   /* Get mapping indices */
1235:   PetscCall(ISLocalToGlobalMappingGetBlockSize(mapping, &bs));
1236:   PetscCall(ISLocalToGlobalMappingGetBlockIndices(mapping, &gidxs));
1237:   PetscCall(ISLocalToGlobalMappingGetSize(mapping, &nleaves));
1238:   nleaves /= bs;

1240:   /* Create layout for global indices */
1241:   for (i = 0, m = 0; i < nleaves; i++) m = PetscMax(m, gidxs[i]);
1242:   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &m, 1, MPIU_INT, MPI_MAX, comm));
1243:   PetscCall(PetscLayoutCreate(comm, &layout));
1244:   PetscCall(PetscLayoutSetSize(layout, m + 1));
1245:   PetscCall(PetscLayoutSetUp(layout));

1247:   /* Create SF to share global indices */
1248:   PetscCall(PetscSFCreate(comm, &sf));
1249:   PetscCall(PetscSFSetGraphLayout(sf, layout, nleaves, NULL, PETSC_OWN_POINTER, gidxs));
1250:   PetscCall(PetscSFSetUp(sf));
1251:   PetscCall(PetscLayoutDestroy(&layout));

1253:   /* communicate root degree to leaves */
1254:   PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, &sfnode));
1255:   PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
1256:   PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
1257:   for (i = 0, mnroots = 0; i < nroots; i++) mnroots += rootdegree[i];
1258:   PetscCall(PetscMalloc3(2 * PetscMax(mnroots, nroots), &mrootdata, 2 * nleaves, &leafdata, nleaves, &leafrd));
1259:   for (i = 0, m = 0; i < nroots; i++) {
1260:     mrootdata[2 * i + 0] = rootdegree[i];
1261:     mrootdata[2 * i + 1] = m;
1262:     m += rootdegree[i];
1263:   }
1264:   PetscCall(PetscSFBcastBegin(sf, MPIU_2INT, mrootdata, leafdata, MPI_REPLACE));
1265:   PetscCall(PetscSFBcastEnd(sf, MPIU_2INT, mrootdata, leafdata, MPI_REPLACE));

1267:   /* allocate enough space to store ranks */
1268:   for (i = 0, newnleaves = 0; i < nleaves; i++) {
1269:     newnleaves += leafdata[2 * i];
1270:     leafrd[i] = leafdata[2 * i];
1271:   }

1273:   /* create new SF nodes to collect multi-root data at leaves */
1274:   PetscCall(PetscMalloc1(newnleaves, &newsfnode));
1275:   for (i = 0, m = 0; i < nleaves; i++) {
1276:     for (j = 0; j < leafrd[i]; j++) {
1277:       newsfnode[m].rank  = sfnode[i].rank;
1278:       newsfnode[m].index = leafdata[2 * i + 1] + j;
1279:       m++;
1280:     }
1281:   }

1283:   /* gather ranks at multi roots */
1284:   for (i = 0; i < mnroots; i++) mrootdata[i] = -1;
1285:   for (i = 0; i < nleaves; i++) leafdata[i] = (PetscInt)rank;

1287:   PetscCall(PetscSFGatherBegin(sf, MPIU_INT, leafdata, mrootdata));
1288:   PetscCall(PetscSFGatherEnd(sf, MPIU_INT, leafdata, mrootdata));

1290:   /* from multi-roots to multi-leaves */
1291:   PetscCall(PetscSFCreate(comm, &sf2));
1292:   PetscCall(PetscSFSetGraph(sf2, mnroots, newnleaves, NULL, PETSC_OWN_POINTER, newsfnode, PETSC_OWN_POINTER));
1293:   PetscCall(PetscSFSetUp(sf2));

1295:   /* broadcast multi-root data to multi-leaves */
1296:   PetscCall(PetscMalloc1(newnleaves, &newleafdata));
1297:   PetscCall(PetscSFBcastBegin(sf2, MPIU_INT, mrootdata, newleafdata, MPI_REPLACE));
1298:   PetscCall(PetscSFBcastEnd(sf2, MPIU_INT, mrootdata, newleafdata, MPI_REPLACE));

1300:   /* sort sharing ranks */
1301:   for (i = 0, m = 0; i < nleaves; i++) {
1302:     PetscCall(PetscSortInt(leafrd[i], newleafdata + m));
1303:     m += leafrd[i];
1304:   }

1306:   /* Number of neighbors and their ranks */
1307:   PetscCall(PetscHMapICreate(&neighs));
1308:   for (i = 0; i < newnleaves; i++) PetscCall(PetscHMapIPut(neighs, newleafdata[i], &iter, &missing));
1309:   PetscCall(PetscHMapIGetSize(neighs, &mapping->info_nproc));
1310:   PetscCall(PetscMalloc1(mapping->info_nproc + 1, &mapping->info_procs));
1311:   PetscCall(PetscHMapIGetKeys(neighs, (i = 0, &i), mapping->info_procs));
1312:   for (i = 0; i < mapping->info_nproc; i++) { /* put info for self first */
1313:     if (mapping->info_procs[i] == rank) {
1314:       PetscInt newr = mapping->info_procs[0];

1316:       mapping->info_procs[0] = rank;
1317:       mapping->info_procs[i] = newr;
1318:       break;
1319:     }
1320:   }
1321:   if (mapping->info_nproc) PetscCall(PetscSortInt(mapping->info_nproc - 1, mapping->info_procs + 1));
1322:   PetscCall(PetscHMapIDestroy(&neighs));

1324:   /* collect info data */
1325:   PetscCall(PetscMalloc1(mapping->info_nproc, &mapping->info_numprocs));
1326:   PetscCall(PetscMalloc1(mapping->info_nproc, &mapping->info_indices));
1327:   for (i = 0; i < mapping->info_nproc; i++) mapping->info_indices[i] = NULL;

1329:   PetscCall(PetscMalloc1(nleaves, &mask));
1330:   PetscCall(PetscMalloc1(nleaves, &tmpg));
1331:   for (p = 0; p < mapping->info_nproc; p++) {
1332:     PetscInt *tmp, trank = mapping->info_procs[p];

1334:     PetscCall(PetscMemzero(mask, nleaves * sizeof(*mask)));
1335:     for (i = 0, m = 0; i < nleaves; i++) {
1336:       for (j = 0; j < leafrd[i]; j++) {
1337:         if (newleafdata[m] == trank) mask[i]++;
1338:         if (!p && newleafdata[m] != rank) mask[i]++;
1339:         m++;
1340:       }
1341:     }
1342:     for (i = 0, m = 0; i < nleaves; i++)
1343:       if (mask[i] > (!p ? 1 : 0)) m++;

1345:     PetscCall(PetscMalloc1(m, &tmp));
1346:     for (i = 0, m = 0; i < nleaves; i++)
1347:       if (mask[i] > (!p ? 1 : 0)) {
1348:         tmp[m]  = i;
1349:         tmpg[m] = gidxs[i];
1350:         m++;
1351:       }
1352:     PetscCall(PetscSortIntWithArray(m, tmpg, tmp));
1353:     mapping->info_indices[p]  = tmp;
1354:     mapping->info_numprocs[p] = m;
1355:   }

1357:   /* Node info */
1358:   PetscCall(PetscMalloc2(nleaves, &mapping->info_nodec, nleaves + 1, &mapping->info_nodei));
1359:   PetscCall(PetscArraycpy(mapping->info_nodec, leafrd, nleaves));
1360:   PetscCall(PetscMalloc1(newnleaves, &mapping->info_nodei[0]));
1361:   for (i = 0; i < nleaves - 1; i++) mapping->info_nodei[i + 1] = mapping->info_nodei[i] + mapping->info_nodec[i];
1362:   PetscCall(PetscArraycpy(mapping->info_nodei[0], newleafdata, newnleaves));

1364:   /* Create SF from leaves to multi-leaves */
1365:   PetscCall(PetscSFGetMultiSF(sf, &msf));
1366:   PetscCall(PetscSFCreateInverseSF(msf, &imsf));
1367:   PetscCall(PetscSFCompose(imsf, sf2, &mapping->multileaves_sf));
1368:   PetscCall(PetscSFDestroy(&imsf));
1369:   PetscCall(PetscSFDestroy(&sf));
1370:   PetscCall(PetscSFDestroy(&sf2));

1372:   PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(mapping, &gidxs));
1373:   PetscCall(PetscFree(tmpg));
1374:   PetscCall(PetscFree(mask));
1375:   PetscCall(PetscFree3(mrootdata, leafdata, leafrd));
1376:   PetscCall(PetscFree(newleafdata));
1377:   PetscFunctionReturn(PETSC_SUCCESS);
1378: }

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

1383:   Not Collective

1385:   Input Parameters:
1386: + mapping  - the mapping from local to global indexing
1387: . nproc    - number of processes that are connected to the calling process
1388: . procs    - neighboring processes
1389: . numprocs - number of block indices for each process
1390: - indices  - block indices (in local numbering) shared with neighbors (sorted by global numbering)

1392:   Level: advanced

1394: .seealso: [](sec_scatter), `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
1395:           `ISLocalToGlobalMappingGetInfo()`
1396: @*/
1397: PetscErrorCode ISLocalToGlobalMappingRestoreBlockInfo(ISLocalToGlobalMapping mapping, PetscInt *nproc, PetscInt *procs[], PetscInt *numprocs[], PetscInt **indices[])
1398: {
1399:   PetscFunctionBegin;
1401:   if (nproc) *nproc = 0;
1402:   if (procs) *procs = NULL;
1403:   if (numprocs) *numprocs = NULL;
1404:   if (indices) *indices = NULL;
1405:   PetscFunctionReturn(PETSC_SUCCESS);
1406: }

1408: /*@C
1409:   ISLocalToGlobalMappingGetInfo - Gets the neighbor information for each process

1411:   Collective the first time it is called

1413:   Input Parameter:
1414: . mapping - the mapping from local to global indexing

1416:   Output Parameters:
1417: + nproc    - number of processes that are connected to the calling process
1418: . procs    - neighboring processes
1419: . numprocs - number of indices for each process
1420: - indices  - indices (in local numbering) shared with neighbors (sorted by global numbering)

1422:   Level: advanced

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

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

1432: .seealso: [](sec_scatter), `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
1433:           `ISLocalToGlobalMappingRestoreInfo()`, `ISLocalToGlobalMappingGetNodeInfo()`
1434: @*/
1435: PetscErrorCode ISLocalToGlobalMappingGetInfo(ISLocalToGlobalMapping mapping, PetscInt *nproc, PetscInt *procs[], PetscInt *numprocs[], PetscInt **indices[])
1436: {
1437:   PetscInt **bindices = NULL, *bnumprocs = NULL, bs, i, j, k, n, *bprocs;

1439:   PetscFunctionBegin;
1441:   bs = mapping->bs;
1442:   PetscCall(ISLocalToGlobalMappingGetBlockInfo(mapping, &n, &bprocs, &bnumprocs, &bindices));
1443:   if (bs > 1) { /* we need to expand the cached info */
1444:     if (indices) PetscCall(PetscCalloc1(n, indices));
1445:     if (numprocs) PetscCall(PetscCalloc1(n, numprocs));
1446:     if (indices || numprocs) {
1447:       for (i = 0; i < n; i++) {
1448:         if (indices) {
1449:           PetscCall(PetscMalloc1(bs * bnumprocs[i], &(*indices)[i]));
1450:           for (j = 0; j < bnumprocs[i]; j++) {
1451:             for (k = 0; k < bs; k++) (*indices)[i][j * bs + k] = bs * bindices[i][j] + k;
1452:           }
1453:         }
1454:         if (numprocs) (*numprocs)[i] = bnumprocs[i] * bs;
1455:       }
1456:     }
1457:   } else {
1458:     if (numprocs) *numprocs = bnumprocs;
1459:     if (indices) *indices = bindices;
1460:   }
1461:   if (nproc) *nproc = n;
1462:   if (procs) *procs = bprocs;
1463:   PetscFunctionReturn(PETSC_SUCCESS);
1464: }

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

1469:   Not Collective

1471:   Input Parameters:
1472: + mapping  - the mapping from local to global indexing
1473: . nproc    - number of processes that are connected to the calling process
1474: . procs    - neighboring processes
1475: . numprocs - number of indices for each process
1476: - indices  - indices (in local numbering) shared with neighbors (sorted by global numbering)

1478:   Level: advanced

1480: .seealso: [](sec_scatter), `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
1481:           `ISLocalToGlobalMappingGetInfo()`
1482: @*/
1483: PetscErrorCode ISLocalToGlobalMappingRestoreInfo(ISLocalToGlobalMapping mapping, PetscInt *nproc, PetscInt *procs[], PetscInt *numprocs[], PetscInt **indices[])
1484: {
1485:   PetscFunctionBegin;
1487:   if (mapping->bs > 1) {
1488:     if (numprocs) PetscCall(PetscFree(*numprocs));
1489:     if (indices) {
1490:       if (*indices)
1491:         for (PetscInt i = 0; i < *nproc; i++) PetscCall(PetscFree((*indices)[i]));
1492:       PetscCall(PetscFree(*indices));
1493:     }
1494:   }
1495:   PetscFunctionReturn(PETSC_SUCCESS);
1496: }

1498: /*@C
1499:   ISLocalToGlobalMappingGetNodeInfo - Gets the neighbor information of local nodes

1501:   Collective the first time it is called

1503:   Input Parameter:
1504: . mapping - the mapping from local to global indexing

1506:   Output Parameters:
1507: + n       - number of local nodes
1508: . n_procs - an array storing the number of processes for each local node (including self)
1509: - procs   - the processes' rank for each local node (sorted, self is first)

1511:   Level: advanced

1513:   Note:
1514:   The user needs to call `ISLocalToGlobalMappingRestoreNodeInfo()` when the data is no longer needed.

1516: .seealso: [](sec_scatter), `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
1517:           `ISLocalToGlobalMappingGetInfo()`, `ISLocalToGlobalMappingRestoreNodeInfo()`, `ISLocalToGlobalMappingGetBlockNodeInfo()`
1518: @*/
1519: PetscErrorCode ISLocalToGlobalMappingGetNodeInfo(ISLocalToGlobalMapping mapping, PetscInt *n, PetscInt *n_procs[], PetscInt **procs[])
1520: {
1521:   PetscInt **bprocs = NULL, *bn_procs = NULL, bs, i, j, k, bn;

1523:   PetscFunctionBegin;
1525:   bs = mapping->bs;
1526:   PetscCall(ISLocalToGlobalMappingGetBlockNodeInfo(mapping, &bn, &bn_procs, &bprocs));
1527:   if (bs > 1) { /* we need to expand the cached info */
1528:     PetscInt *tn_procs;
1529:     PetscInt  c;

1531:     PetscCall(PetscMalloc1(bn * bs, &tn_procs));
1532:     for (i = 0, c = 0; i < bn; i++) {
1533:       for (k = 0; k < bs; k++) tn_procs[i * bs + k] = bn_procs[i];
1534:       c += bs * bn_procs[i];
1535:     }
1536:     if (n) *n = bn * bs;
1537:     if (procs) {
1538:       PetscInt **tprocs;
1539:       PetscInt   tn = bn * bs;

1541:       PetscCall(PetscMalloc1(tn, &tprocs));
1542:       if (tn) PetscCall(PetscMalloc1(c, &tprocs[0]));
1543:       for (i = 0; i < tn - 1; i++) tprocs[i + 1] = tprocs[i] + tn_procs[i];
1544:       for (i = 0; i < bn; i++) {
1545:         for (k = 0; k < bs; k++) {
1546:           for (j = 0; j < bn_procs[i]; j++) tprocs[i * bs + k][j] = bprocs[i][j];
1547:         }
1548:       }
1549:       *procs = tprocs;
1550:     }
1551:     if (n_procs) *n_procs = tn_procs;
1552:     else PetscCall(PetscFree(tn_procs));
1553:   } else {
1554:     if (n) *n = bn;
1555:     if (n_procs) *n_procs = bn_procs;
1556:     if (procs) *procs = bprocs;
1557:   }
1558:   PetscFunctionReturn(PETSC_SUCCESS);
1559: }

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

1564:   Not Collective

1566:   Input Parameters:
1567: + mapping - the mapping from local to global indexing
1568: . n       - number of local nodes
1569: . n_procs - an array storing the number of processes for each local node (including self)
1570: - procs   - the processes' rank for each local node (sorted, self is first)

1572:   Level: advanced

1574: .seealso: [](sec_scatter), `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
1575:           `ISLocalToGlobalMappingGetInfo()`
1576: @*/
1577: PetscErrorCode ISLocalToGlobalMappingRestoreNodeInfo(ISLocalToGlobalMapping mapping, PetscInt *n, PetscInt *n_procs[], PetscInt **procs[])
1578: {
1579:   PetscFunctionBegin;
1581:   if (mapping->bs > 1) {
1582:     if (n_procs) PetscCall(PetscFree(*n_procs));
1583:     if (procs) {
1584:       if (*procs) PetscCall(PetscFree((*procs)[0]));
1585:       PetscCall(PetscFree(*procs));
1586:     }
1587:   }
1588:   PetscCall(ISLocalToGlobalMappingRestoreBlockNodeInfo(mapping, n, n_procs, procs));
1589:   PetscFunctionReturn(PETSC_SUCCESS);
1590: }

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

1595:     Synopsis:
1596:     ISLocalToGlobalMappingGetIndicesF90(ISLocalToGlobalMapping ltog, PetscInt, pointer :: array(:)}, integer ierr)

1598:     Not Collective

1600:     Input Parameter:
1601: .   A - the matrix

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

1606:     Level: advanced

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

1611: .seealso: [](sec_fortranarrays), [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingGetIndices()`, `ISLocalToGlobalMappingRestoreIndices()`,
1612:           `ISLocalToGlobalMappingRestoreIndicesF90()`
1613: M*/

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

1618:     Synopsis:
1619:     ISLocalToGlobalMappingRestoreIndicesF90(ISLocalToGlobalMapping ltog, PetscInt, pointer :: array(:)}, integer ierr)

1621:     Not Collective

1623:     Input Parameters:
1624: +   A - the matrix
1625: -   array - array of indices, the length of this array may be obtained with `ISLocalToGlobalMappingGetSize()`

1627:     Level: advanced

1629: .seealso: [](sec_fortranarrays), [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingGetIndices()`, `ISLocalToGlobalMappingRestoreIndices()`,
1630:           `ISLocalToGlobalMappingGetIndicesF90()`
1631: M*/

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

1636:   Not Collective

1638:   Input Parameter:
1639: . ltog - local to global mapping

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

1644:   Level: advanced

1646:   Note:
1647:   `ISLocalToGlobalMappingGetSize()` returns the length the this array

1649: .seealso: [](sec_scatter), `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingApply()`, `ISLocalToGlobalMappingRestoreIndices()`,
1650:           `ISLocalToGlobalMappingGetBlockIndices()`, `ISLocalToGlobalMappingRestoreBlockIndices()`
1651: @*/
1652: PetscErrorCode ISLocalToGlobalMappingGetIndices(ISLocalToGlobalMapping ltog, const PetscInt *array[])
1653: {
1654:   PetscFunctionBegin;
1656:   PetscAssertPointer(array, 2);
1657:   if (ltog->bs == 1) {
1658:     *array = ltog->indices;
1659:   } else {
1660:     PetscInt       *jj, k, i, j, n = ltog->n, bs = ltog->bs;
1661:     const PetscInt *ii;

1663:     PetscCall(PetscMalloc1(bs * n, &jj));
1664:     *array = jj;
1665:     k      = 0;
1666:     ii     = ltog->indices;
1667:     for (i = 0; i < n; i++)
1668:       for (j = 0; j < bs; j++) jj[k++] = bs * ii[i] + j;
1669:   }
1670:   PetscFunctionReturn(PETSC_SUCCESS);
1671: }

1673: /*@C
1674:   ISLocalToGlobalMappingRestoreIndices - Restore indices obtained with `ISLocalToGlobalMappingGetIndices()`

1676:   Not Collective

1678:   Input Parameters:
1679: + ltog  - local to global mapping
1680: - array - array of indices

1682:   Level: advanced

1684: .seealso: [](sec_scatter), `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingApply()`, `ISLocalToGlobalMappingGetIndices()`
1685: @*/
1686: PetscErrorCode ISLocalToGlobalMappingRestoreIndices(ISLocalToGlobalMapping ltog, const PetscInt *array[])
1687: {
1688:   PetscFunctionBegin;
1690:   PetscAssertPointer(array, 2);
1691:   PetscCheck(ltog->bs != 1 || *array == ltog->indices, PETSC_COMM_SELF, PETSC_ERR_ARG_BADPTR, "Trying to return mismatched pointer");
1692:   if (ltog->bs > 1) PetscCall(PetscFree(*(void **)array));
1693:   PetscFunctionReturn(PETSC_SUCCESS);
1694: }

1696: /*@C
1697:   ISLocalToGlobalMappingGetBlockIndices - Get global indices for every local block

1699:   Not Collective

1701:   Input Parameter:
1702: . ltog - local to global mapping

1704:   Output Parameter:
1705: . array - array of indices

1707:   Level: advanced

1709: .seealso: [](sec_scatter), `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingApply()`, `ISLocalToGlobalMappingRestoreBlockIndices()`
1710: @*/
1711: PetscErrorCode ISLocalToGlobalMappingGetBlockIndices(ISLocalToGlobalMapping ltog, const PetscInt *array[])
1712: {
1713:   PetscFunctionBegin;
1715:   PetscAssertPointer(array, 2);
1716:   *array = ltog->indices;
1717:   PetscFunctionReturn(PETSC_SUCCESS);
1718: }

1720: /*@C
1721:   ISLocalToGlobalMappingRestoreBlockIndices - Restore indices obtained with `ISLocalToGlobalMappingGetBlockIndices()`

1723:   Not Collective

1725:   Input Parameters:
1726: + ltog  - local to global mapping
1727: - array - array of indices

1729:   Level: advanced

1731: .seealso: [](sec_scatter), `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingApply()`, `ISLocalToGlobalMappingGetIndices()`
1732: @*/
1733: PetscErrorCode ISLocalToGlobalMappingRestoreBlockIndices(ISLocalToGlobalMapping ltog, const PetscInt *array[])
1734: {
1735:   PetscFunctionBegin;
1737:   PetscAssertPointer(array, 2);
1738:   PetscCheck(*array == ltog->indices, PETSC_COMM_SELF, PETSC_ERR_ARG_BADPTR, "Trying to return mismatched pointer");
1739:   *array = NULL;
1740:   PetscFunctionReturn(PETSC_SUCCESS);
1741: }

1743: /*@
1744:   ISLocalToGlobalMappingConcatenate - Create a new mapping that concatenates a list of mappings

1746:   Not Collective

1748:   Input Parameters:
1749: + comm  - communicator for the new mapping, must contain the communicator of every mapping to concatenate
1750: . n     - number of mappings to concatenate
1751: - ltogs - local to global mappings

1753:   Output Parameter:
1754: . ltogcat - new mapping

1756:   Level: advanced

1758:   Note:
1759:   This currently always returns a mapping with block size of 1

1761:   Developer Notes:
1762:   If all the input mapping have the same block size we could easily handle that as a special case

1764: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingCreate()`
1765: @*/
1766: PetscErrorCode ISLocalToGlobalMappingConcatenate(MPI_Comm comm, PetscInt n, const ISLocalToGlobalMapping ltogs[], ISLocalToGlobalMapping *ltogcat)
1767: {
1768:   PetscInt i, cnt, m, *idx;

1770:   PetscFunctionBegin;
1771:   PetscCheck(n >= 0, comm, PETSC_ERR_ARG_OUTOFRANGE, "Must have a non-negative number of mappings, given %" PetscInt_FMT, n);
1772:   if (n > 0) PetscAssertPointer(ltogs, 3);
1774:   PetscAssertPointer(ltogcat, 4);
1775:   for (cnt = 0, i = 0; i < n; i++) {
1776:     PetscCall(ISLocalToGlobalMappingGetSize(ltogs[i], &m));
1777:     cnt += m;
1778:   }
1779:   PetscCall(PetscMalloc1(cnt, &idx));
1780:   for (cnt = 0, i = 0; i < n; i++) {
1781:     const PetscInt *subidx;
1782:     PetscCall(ISLocalToGlobalMappingGetSize(ltogs[i], &m));
1783:     PetscCall(ISLocalToGlobalMappingGetIndices(ltogs[i], &subidx));
1784:     PetscCall(PetscArraycpy(&idx[cnt], subidx, m));
1785:     PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogs[i], &subidx));
1786:     cnt += m;
1787:   }
1788:   PetscCall(ISLocalToGlobalMappingCreate(comm, 1, cnt, idx, PETSC_OWN_POINTER, ltogcat));
1789:   PetscFunctionReturn(PETSC_SUCCESS);
1790: }

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

1796:    Options Database Key:
1797: .   -islocaltoglobalmapping_type basic - select this method

1799:    Level: beginner

1801:    Developer Note:
1802:    This stores all the mapping information on each MPI rank.

1804: .seealso: [](sec_scatter), `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingSetType()`, `ISLOCALTOGLOBALMAPPINGHASH`
1805: M*/
1806: PETSC_EXTERN PetscErrorCode ISLocalToGlobalMappingCreate_Basic(ISLocalToGlobalMapping ltog)
1807: {
1808:   PetscFunctionBegin;
1809:   ltog->ops->globaltolocalmappingapply      = ISGlobalToLocalMappingApply_Basic;
1810:   ltog->ops->globaltolocalmappingsetup      = ISGlobalToLocalMappingSetUp_Basic;
1811:   ltog->ops->globaltolocalmappingapplyblock = ISGlobalToLocalMappingApplyBlock_Basic;
1812:   ltog->ops->destroy                        = ISLocalToGlobalMappingDestroy_Basic;
1813:   PetscFunctionReturn(PETSC_SUCCESS);
1814: }

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

1820:    Options Database Key:
1821: .   -islocaltoglobalmapping_type hash - select this method

1823:    Level: beginner

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

1828: .seealso: [](sec_scatter), `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingSetType()`, `ISLOCALTOGLOBALMAPPINGBASIC`
1829: M*/
1830: PETSC_EXTERN PetscErrorCode ISLocalToGlobalMappingCreate_Hash(ISLocalToGlobalMapping ltog)
1831: {
1832:   PetscFunctionBegin;
1833:   ltog->ops->globaltolocalmappingapply      = ISGlobalToLocalMappingApply_Hash;
1834:   ltog->ops->globaltolocalmappingsetup      = ISGlobalToLocalMappingSetUp_Hash;
1835:   ltog->ops->globaltolocalmappingapplyblock = ISGlobalToLocalMappingApplyBlock_Hash;
1836:   ltog->ops->destroy                        = ISLocalToGlobalMappingDestroy_Hash;
1837:   PetscFunctionReturn(PETSC_SUCCESS);
1838: }

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

1843:   Not Collective, No Fortran Support

1845:   Input Parameters:
1846: + sname    - name of a new method
1847: - function - routine to create method context

1849:   Example Usage:
1850: .vb
1851:    ISLocalToGlobalMappingRegister("my_mapper", MyCreate);
1852: .ve

1854:   Then, your mapping can be chosen with the procedural interface via
1855: $     ISLocalToGlobalMappingSetType(ltog, "my_mapper")
1856:   or at runtime via the option
1857: $     -islocaltoglobalmapping_type my_mapper

1859:   Level: advanced

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

1864: .seealso: [](sec_scatter), `ISLocalToGlobalMappingRegisterAll()`, `ISLocalToGlobalMappingRegisterDestroy()`, `ISLOCALTOGLOBALMAPPINGBASIC`,
1865:           `ISLOCALTOGLOBALMAPPINGHASH`, `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingApply()`
1866: @*/
1867: PetscErrorCode ISLocalToGlobalMappingRegister(const char sname[], PetscErrorCode (*function)(ISLocalToGlobalMapping))
1868: {
1869:   PetscFunctionBegin;
1870:   PetscCall(ISInitializePackage());
1871:   PetscCall(PetscFunctionListAdd(&ISLocalToGlobalMappingList, sname, function));
1872:   PetscFunctionReturn(PETSC_SUCCESS);
1873: }

1875: /*@
1876:   ISLocalToGlobalMappingSetType - Sets the implementation type `ISLocalToGlobalMapping` will use

1878:   Logically Collective

1880:   Input Parameters:
1881: + ltog - the `ISLocalToGlobalMapping` object
1882: - type - a known method

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

1887:   Level: intermediate

1889:   Notes:
1890:   See `ISLocalToGlobalMappingType` for available methods

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

1896:   Developer Notes:
1897:   `ISLocalToGlobalMappingRegister()` is used to add new types to `ISLocalToGlobalMappingList` from which they
1898:   are accessed by `ISLocalToGlobalMappingSetType()`.

1900: .seealso: [](sec_scatter), `ISLocalToGlobalMappingType`, `ISLocalToGlobalMappingRegister()`, `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingGetType()`
1901: @*/
1902: PetscErrorCode ISLocalToGlobalMappingSetType(ISLocalToGlobalMapping ltog, ISLocalToGlobalMappingType type)
1903: {
1904:   PetscBool match;
1905:   PetscErrorCode (*r)(ISLocalToGlobalMapping) = NULL;

1907:   PetscFunctionBegin;
1909:   if (type) PetscAssertPointer(type, 2);

1911:   PetscCall(PetscObjectTypeCompare((PetscObject)ltog, type, &match));
1912:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

1914:   /* L2G maps defer type setup at globaltolocal calls, allow passing NULL here */
1915:   if (type) {
1916:     PetscCall(PetscFunctionListFind(ISLocalToGlobalMappingList, type, &r));
1917:     PetscCheck(r, PetscObjectComm((PetscObject)ltog), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unable to find requested ISLocalToGlobalMapping type %s", type);
1918:   }
1919:   /* Destroy the previous private LTOG context */
1920:   PetscTryTypeMethod(ltog, destroy);
1921:   ltog->ops->destroy = NULL;

1923:   PetscCall(PetscObjectChangeTypeName((PetscObject)ltog, type));
1924:   if (r) PetscCall((*r)(ltog));
1925:   PetscFunctionReturn(PETSC_SUCCESS);
1926: }

1928: /*@
1929:   ISLocalToGlobalMappingGetType - Get the type of the `ISLocalToGlobalMapping`

1931:   Not Collective

1933:   Input Parameter:
1934: . ltog - the `ISLocalToGlobalMapping` object

1936:   Output Parameter:
1937: . type - the type

1939:   Level: intermediate

1941: .seealso: [](sec_scatter), `ISLocalToGlobalMappingType`, `ISLocalToGlobalMappingRegister()`, `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingSetType()`
1942: @*/
1943: PetscErrorCode ISLocalToGlobalMappingGetType(ISLocalToGlobalMapping ltog, ISLocalToGlobalMappingType *type)
1944: {
1945:   PetscFunctionBegin;
1947:   PetscAssertPointer(type, 2);
1948:   *type = ((PetscObject)ltog)->type_name;
1949:   PetscFunctionReturn(PETSC_SUCCESS);
1950: }

1952: PetscBool ISLocalToGlobalMappingRegisterAllCalled = PETSC_FALSE;

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

1957:   Not Collective

1959:   Level: advanced

1961: .seealso: [](sec_scatter), `ISRegister()`, `ISLocalToGlobalRegister()`
1962: @*/
1963: PetscErrorCode ISLocalToGlobalMappingRegisterAll(void)
1964: {
1965:   PetscFunctionBegin;
1966:   if (ISLocalToGlobalMappingRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
1967:   ISLocalToGlobalMappingRegisterAllCalled = PETSC_TRUE;
1968:   PetscCall(ISLocalToGlobalMappingRegister(ISLOCALTOGLOBALMAPPINGBASIC, ISLocalToGlobalMappingCreate_Basic));
1969:   PetscCall(ISLocalToGlobalMappingRegister(ISLOCALTOGLOBALMAPPINGHASH, ISLocalToGlobalMappingCreate_Hash));
1970:   PetscFunctionReturn(PETSC_SUCCESS);
1971: }