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: /*
126:     Creates the global mapping information in the ISLocalToGlobalMapping structure

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

134:   PetscFunctionBegin;
135:   if (mapping->data) PetscFunctionReturn(PETSC_SUCCESS);
136:   end   = 0;
137:   start = PETSC_INT_MAX;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

273: /*@
274:   ISLocalToGlobalMappingDuplicate - Duplicates the local to global mapping object

276:   Not Collective

278:   Input Parameter:
279: . ltog - local to global mapping

281:   Output Parameter:
282: . nltog - the duplicated local to global mapping

284:   Level: advanced

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

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

300: /*@
301:   ISLocalToGlobalMappingGetSize - Gets the local size of a local to global mapping

303:   Not Collective

305:   Input Parameter:
306: . mapping - local to global mapping

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

311:   Level: advanced

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

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

327:   Collective

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

334:   Level: intermediate

336:   Note:
337:   See `PetscObjectViewFromOptions()` for the available `PetscViewer` and `PetscViewerFormat`

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

349: /*@
350:   ISLocalToGlobalMappingView - View a local to global mapping

352:   Collective on viewer

354:   Input Parameters:
355: + mapping - local to global mapping
356: - viewer  - viewer

358:   Level: intermediate

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

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

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

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

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

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

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

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

428:   Collective on viewer

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

434:   Level: intermediate

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

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

448:   /* reset previous data */
449:   PetscCall(ISLocalToGlobalMappingResetBlockInfo_Private(mapping));

451:   PetscCall(PetscViewerSetUp(viewer));
452:   PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));

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

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

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

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

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

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

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

492:   Not Collective

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

497:   Output Parameter:
498: . mapping - new mapping data structure

500:   Level: advanced

502:   Note:
503:   the block size of the `IS` determines the block size of the mapping

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

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

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

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

537:   Collective

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

543:   Output Parameter:
544: . mapping - new mapping data structure

546:   Level: advanced

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

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

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

580: /*@
581:   ISLocalToGlobalMappingSetBlockSize - Sets the blocksize of the mapping

583:   Not Collective

585:   Input Parameters:
586: + mapping - mapping data structure
587: - bs      - the blocksize

589:   Level: advanced

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

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

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

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

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

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

647:   Not Collective

649:   Input Parameter:
650: . mapping - mapping data structure

652:   Output Parameter:
653: . bs - the blocksize

655:   Level: advanced

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

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

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

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

680:   Output Parameter:
681: . mapping - new mapping data structure

683:   Level: advanced

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

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

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

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

703:   PetscFunctionBegin;
704:   if (n) PetscAssertPointer(indices, 4);
705:   PetscAssertPointer(mapping, 6);

707:   *mapping = NULL;
708:   PetscCall(ISInitializePackage());

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

727: PetscFunctionList ISLocalToGlobalMappingList = NULL;

729: /*@
730:   ISLocalToGlobalMappingSetFromOptions - Set mapping options from the options database.

732:   Not Collective

734:   Input Parameter:
735: . mapping - mapping data structure

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

740:   Level: advanced

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

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

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

766:   Not Collective

768:   Input Parameter:
769: . mapping - mapping data structure

771:   Level: advanced

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

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

796:   Collective

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

802:   Output Parameter:
803: . newis - index set in global numbering

805:   Level: advanced

807:   Note:
808:   The output `IS` will have the same communicator as the input `IS` as well as the same block size.

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

818:   PetscFunctionBegin;
821:   PetscAssertPointer(newis, 3);

823:   PetscCall(ISGetLocalSize(is, &n));
824:   PetscCall(ISGetBlockSize(is, &bs));
825:   PetscCall(ISGetIndices(is, &idxin));
826:   PetscCall(PetscMalloc1(n, &idxout));
827:   PetscCall(ISLocalToGlobalMappingApply(mapping, n, idxin, idxout));
828:   PetscCall(ISRestoreIndices(is, &idxin));
829:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is), n, idxout, PETSC_OWN_POINTER, newis));
830:   PetscCall(ISSetBlockSize(*newis, bs));
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:   Notes:
1001:   The output `IS` will be sequential, as it encodes a purely local operation

1003:   If `type` is `IS_GTOLM_MASK`, `newis` will have the same block size as `is`

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

1013:   PetscFunctionBegin;
1016:   PetscAssertPointer(newis, 4);

1018:   PetscCall(ISGetLocalSize(is, &n));
1019:   PetscCall(ISGetIndices(is, &idxin));
1020:   if (type == IS_GTOLM_MASK) {
1021:     PetscCall(PetscMalloc1(n, &idxout));
1022:   } else {
1023:     PetscCall(ISGlobalToLocalMappingApply(mapping, type, n, idxin, &nout, NULL));
1024:     PetscCall(PetscMalloc1(nout, &idxout));
1025:   }
1026:   PetscCall(ISGlobalToLocalMappingApply(mapping, type, n, idxin, &nout, idxout));
1027:   PetscCall(ISRestoreIndices(is, &idxin));
1028:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nout, idxout, PETSC_OWN_POINTER, newis));
1029:   if (type == IS_GTOLM_MASK) {
1030:     PetscCall(ISGetBlockSize(is, &bs));
1031:     PetscCall(ISSetBlockSize(*newis, bs));
1032:   }
1033:   PetscFunctionReturn(PETSC_SUCCESS);
1034: }

1036: /*@
1037:   ISGlobalToLocalMappingApplyBlock - Provides the local block numbering for a list of integers
1038:   specified with a block global numbering.

1040:   Not Collective

1042:   Input Parameters:
1043: + mapping - mapping between local and global numbering
1044: . type    - `IS_GTOLM_MASK` - maps global indices with no local value to -1 in the output list (i.e., mask them)
1045:            `IS_GTOLM_DROP` - drops the indices with no local value from the output list
1046: . n       - number of global indices to map
1047: - idx     - global indices to map

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

1057:   Level: advanced

1059:   Notes:
1060:   Either `nout` or `idxout` may be `NULL`. `idx` and `idxout` may be identical.

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

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

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

1083: /*@C
1084:   ISLocalToGlobalMappingGetBlockInfo - Gets the neighbor information

1086:   Collective the first time it is called

1088:   Input Parameter:
1089: . mapping - the mapping from local to global indexing

1091:   Output Parameters:
1092: + nproc    - number of processes that are connected to the calling process
1093: . procs    - neighboring processes
1094: . numprocs - number of block indices for each process
1095: - indices  - block indices (in local numbering) shared with neighbors (sorted by global numbering)

1097:   Level: advanced

1099: .seealso: [](sec_scatter), `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
1100:           `ISLocalToGlobalMappingRestoreBlockInfo()`, `ISLocalToGlobalMappingGetBlockMultiLeavesSF()`
1101: @*/
1102: PetscErrorCode ISLocalToGlobalMappingGetBlockInfo(ISLocalToGlobalMapping mapping, PetscInt *nproc, PetscInt *procs[], PetscInt *numprocs[], PetscInt **indices[])
1103: {
1104:   PetscFunctionBegin;
1106:   PetscCall(ISLocalToGlobalMappingSetUpBlockInfo_Private(mapping));
1107:   if (nproc) *nproc = mapping->info_nproc;
1108:   if (procs) *procs = mapping->info_procs;
1109:   if (numprocs) *numprocs = mapping->info_numprocs;
1110:   if (indices) *indices = mapping->info_indices;
1111:   PetscFunctionReturn(PETSC_SUCCESS);
1112: }

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

1117:   Collective the first time it is called

1119:   Input Parameter:
1120: . mapping - the mapping from local to global indexing

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

1127:   Level: advanced

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

1135: .seealso: `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
1136:           `ISLocalToGlobalMappingGetBlockInfo()`, `ISLocalToGlobalMappingRestoreBlockNodeInfo()`, `ISLocalToGlobalMappingGetNodeInfo()`
1137: @*/
1138: PetscErrorCode ISLocalToGlobalMappingGetBlockNodeInfo(ISLocalToGlobalMapping mapping, PetscInt *n, PetscInt *n_procs[], PetscInt **procs[])
1139: {
1140:   PetscFunctionBegin;
1142:   PetscCall(ISLocalToGlobalMappingSetUpBlockInfo_Private(mapping));
1143:   if (n) *n = mapping->n;
1144:   if (n_procs) *n_procs = mapping->info_nodec;
1145:   if (procs) *procs = mapping->info_nodei;
1146:   PetscFunctionReturn(PETSC_SUCCESS);
1147: }

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

1152:   Not Collective

1154:   Input Parameters:
1155: + mapping - the mapping from local to global indexing
1156: . n       - number of local block nodes
1157: . n_procs - an array storing the number of processes for each local block nodes (including self)
1158: - procs   - the processes' rank for each local block node (sorted, self is first)

1160:   Level: advanced

1162: .seealso: `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
1163:           `ISLocalToGlobalMappingGetBlockNodeInfo()`
1164: @*/
1165: PetscErrorCode ISLocalToGlobalMappingRestoreBlockNodeInfo(ISLocalToGlobalMapping mapping, PetscInt *n, PetscInt *n_procs[], PetscInt **procs[])
1166: {
1167:   PetscFunctionBegin;
1169:   if (n) *n = 0;
1170:   if (n_procs) *n_procs = NULL;
1171:   if (procs) *procs = NULL;
1172:   PetscFunctionReturn(PETSC_SUCCESS);
1173: }

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

1178:   Collective the first time it is called

1180:   Input Parameter:
1181: . mapping - the mapping from local to global indexing

1183:   Output Parameter:
1184: . mlsf - the `PetscSF`

1186:   Level: advanced

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

1207: .seealso: [](sec_scatter), `ISLocalToGlobalMappingGetBlockNodeInfo()`, `PetscSF`
1208: @*/
1209: PetscErrorCode ISLocalToGlobalMappingGetBlockMultiLeavesSF(ISLocalToGlobalMapping mapping, PetscSF *mlsf)
1210: {
1211:   PetscFunctionBegin;
1213:   PetscAssertPointer(mlsf, 2);
1214:   PetscCall(ISLocalToGlobalMappingSetUpBlockInfo_Private(mapping));
1215:   *mlsf = mapping->multileaves_sf;
1216:   PetscFunctionReturn(PETSC_SUCCESS);
1217: }

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

1234:   PetscFunctionBegin;
1235:   if (mapping->multileaves_sf) PetscFunctionReturn(PETSC_SUCCESS);
1236:   PetscCall(PetscObjectGetComm((PetscObject)mapping, &comm));
1237:   PetscCallMPI(MPI_Comm_size(comm, &size));
1238:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

1240:   /* Get mapping indices */
1241:   PetscCall(ISLocalToGlobalMappingGetBlockSize(mapping, &bs));
1242:   PetscCall(ISLocalToGlobalMappingGetBlockIndices(mapping, &gidxs));
1243:   PetscCall(ISLocalToGlobalMappingGetSize(mapping, &nleaves));
1244:   nleaves /= bs;

1246:   /* Create layout for global indices */
1247:   for (i = 0, m = 0; i < nleaves; i++) m = PetscMax(m, gidxs[i]);
1248:   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &m, 1, MPIU_INT, MPI_MAX, comm));
1249:   PetscCall(PetscLayoutCreate(comm, &layout));
1250:   PetscCall(PetscLayoutSetSize(layout, m + 1));
1251:   PetscCall(PetscLayoutSetUp(layout));

1253:   /* Create SF to share global indices */
1254:   PetscCall(PetscSFCreate(comm, &sf));
1255:   PetscCall(PetscSFSetGraphLayout(sf, layout, nleaves, NULL, PETSC_OWN_POINTER, gidxs));
1256:   PetscCall(PetscSFSetUp(sf));
1257:   PetscCall(PetscLayoutDestroy(&layout));

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

1273:   /* allocate enough space to store ranks */
1274:   for (i = 0, newnleaves = 0; i < nleaves; i++) {
1275:     newnleaves += leafdata[2 * i];
1276:     leafrd[i] = leafdata[2 * i];
1277:   }

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

1289:   /* gather ranks at multi roots */
1290:   for (i = 0; i < mnroots; i++) mrootdata[i] = -1;
1291:   for (i = 0; i < nleaves; i++) leafdata[i] = rank;

1293:   PetscCall(PetscSFGatherBegin(sf, MPIU_INT, leafdata, mrootdata));
1294:   PetscCall(PetscSFGatherEnd(sf, MPIU_INT, leafdata, mrootdata));

1296:   /* from multi-roots to multi-leaves */
1297:   PetscCall(PetscSFCreate(comm, &sf2));
1298:   PetscCall(PetscSFSetGraph(sf2, mnroots, newnleaves, NULL, PETSC_OWN_POINTER, newsfnode, PETSC_OWN_POINTER));
1299:   PetscCall(PetscSFSetUp(sf2));

1301:   /* broadcast multi-root data to multi-leaves */
1302:   PetscCall(PetscMalloc1(newnleaves, &newleafdata));
1303:   PetscCall(PetscSFBcastBegin(sf2, MPIU_INT, mrootdata, newleafdata, MPI_REPLACE));
1304:   PetscCall(PetscSFBcastEnd(sf2, MPIU_INT, mrootdata, newleafdata, MPI_REPLACE));

1306:   /* sort sharing ranks */
1307:   for (i = 0, m = 0; i < nleaves; i++) {
1308:     PetscCall(PetscSortInt(leafrd[i], newleafdata + m));
1309:     m += leafrd[i];
1310:   }

1312:   /* Number of neighbors and their ranks */
1313:   PetscCall(PetscHMapICreate(&neighs));
1314:   for (i = 0; i < newnleaves; i++) PetscCall(PetscHMapIPut(neighs, newleafdata[i], &iter, &missing));
1315:   PetscCall(PetscHMapIGetSize(neighs, &mapping->info_nproc));
1316:   PetscCall(PetscMalloc1(mapping->info_nproc, &mapping->info_procs));
1317:   PetscCall(PetscHMapIGetKeys(neighs, (i = 0, &i), mapping->info_procs));
1318:   for (i = 0; i < mapping->info_nproc; i++) { /* put info for self first */
1319:     if (mapping->info_procs[i] == rank) {
1320:       PetscInt newr = mapping->info_procs[0];

1322:       mapping->info_procs[0] = rank;
1323:       mapping->info_procs[i] = newr;
1324:       break;
1325:     }
1326:   }
1327:   if (mapping->info_nproc) PetscCall(PetscSortInt(mapping->info_nproc - 1, mapping->info_procs + 1));
1328:   PetscCall(PetscHMapIDestroy(&neighs));

1330:   /* collect info data */
1331:   PetscCall(PetscMalloc1(mapping->info_nproc, &mapping->info_numprocs));
1332:   PetscCall(PetscMalloc1(mapping->info_nproc, &mapping->info_indices));
1333:   for (i = 0; i < mapping->info_nproc; i++) mapping->info_indices[i] = NULL;

1335:   PetscCall(PetscMalloc1(nleaves, &mask));
1336:   PetscCall(PetscMalloc1(nleaves, &tmpg));
1337:   for (p = 0; p < mapping->info_nproc; p++) {
1338:     PetscInt *tmp, trank = mapping->info_procs[p];

1340:     PetscCall(PetscMemzero(mask, nleaves * sizeof(*mask)));
1341:     for (i = 0, m = 0; i < nleaves; i++) {
1342:       for (j = 0; j < leafrd[i]; j++) {
1343:         if (newleafdata[m] == trank) mask[i]++;
1344:         if (!p && newleafdata[m] != rank) mask[i]++;
1345:         m++;
1346:       }
1347:     }
1348:     for (i = 0, m = 0; i < nleaves; i++)
1349:       if (mask[i] > (!p ? 1 : 0)) m++;

1351:     PetscCall(PetscMalloc1(m, &tmp));
1352:     for (i = 0, m = 0; i < nleaves; i++)
1353:       if (mask[i] > (!p ? 1 : 0)) {
1354:         tmp[m]  = i;
1355:         tmpg[m] = gidxs[i];
1356:         m++;
1357:       }
1358:     PetscCall(PetscSortIntWithArray(m, tmpg, tmp));
1359:     mapping->info_indices[p]  = tmp;
1360:     mapping->info_numprocs[p] = m;
1361:   }

1363:   /* Node info */
1364:   PetscCall(PetscMalloc2(nleaves, &mapping->info_nodec, nleaves + 1, &mapping->info_nodei));
1365:   PetscCall(PetscArraycpy(mapping->info_nodec, leafrd, nleaves));
1366:   PetscCall(PetscMalloc1(newnleaves, &mapping->info_nodei[0]));
1367:   for (i = 0; i < nleaves - 1; i++) mapping->info_nodei[i + 1] = mapping->info_nodei[i] + mapping->info_nodec[i];
1368:   PetscCall(PetscArraycpy(mapping->info_nodei[0], newleafdata, newnleaves));

1370:   /* Create SF from leaves to multi-leaves */
1371:   PetscCall(PetscSFGetMultiSF(sf, &msf));
1372:   PetscCall(PetscSFCreateInverseSF(msf, &imsf));
1373:   PetscCall(PetscSFCompose(imsf, sf2, &mapping->multileaves_sf));
1374:   PetscCall(PetscSFDestroy(&imsf));
1375:   PetscCall(PetscSFDestroy(&sf));
1376:   PetscCall(PetscSFDestroy(&sf2));

1378:   PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(mapping, &gidxs));
1379:   PetscCall(PetscFree(tmpg));
1380:   PetscCall(PetscFree(mask));
1381:   PetscCall(PetscFree3(mrootdata, leafdata, leafrd));
1382:   PetscCall(PetscFree(newleafdata));
1383:   PetscFunctionReturn(PETSC_SUCCESS);
1384: }

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

1389:   Not Collective

1391:   Input Parameters:
1392: + mapping  - the mapping from local to global indexing
1393: . nproc    - number of processes that are connected to the calling process
1394: . procs    - neighboring processes
1395: . numprocs - number of block indices for each process
1396: - indices  - block indices (in local numbering) shared with neighbors (sorted by global numbering)

1398:   Level: advanced

1400: .seealso: [](sec_scatter), `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
1401:           `ISLocalToGlobalMappingGetInfo()`
1402: @*/
1403: PetscErrorCode ISLocalToGlobalMappingRestoreBlockInfo(ISLocalToGlobalMapping mapping, PetscInt *nproc, PetscInt *procs[], PetscInt *numprocs[], PetscInt **indices[])
1404: {
1405:   PetscFunctionBegin;
1407:   if (nproc) *nproc = 0;
1408:   if (procs) *procs = NULL;
1409:   if (numprocs) *numprocs = NULL;
1410:   if (indices) *indices = NULL;
1411:   PetscFunctionReturn(PETSC_SUCCESS);
1412: }

1414: /*@C
1415:   ISLocalToGlobalMappingGetInfo - Gets the neighbor information for each process

1417:   Collective the first time it is called

1419:   Input Parameter:
1420: . mapping - the mapping from local to global indexing

1422:   Output Parameters:
1423: + nproc    - number of processes that are connected to the calling process
1424: . procs    - neighboring processes
1425: . numprocs - number of indices for each process
1426: - indices  - indices (in local numbering) shared with neighbors (sorted by global numbering)

1428:   Level: advanced

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

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

1438: .seealso: [](sec_scatter), `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
1439:           `ISLocalToGlobalMappingRestoreInfo()`, `ISLocalToGlobalMappingGetNodeInfo()`
1440: @*/
1441: PetscErrorCode ISLocalToGlobalMappingGetInfo(ISLocalToGlobalMapping mapping, PetscInt *nproc, PetscInt *procs[], PetscInt *numprocs[], PetscInt **indices[])
1442: {
1443:   PetscInt **bindices = NULL, *bnumprocs = NULL, bs, i, j, k, n, *bprocs;

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

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

1475:   Not Collective

1477:   Input Parameters:
1478: + mapping  - the mapping from local to global indexing
1479: . nproc    - number of processes that are connected to the calling process
1480: . procs    - neighboring processes
1481: . numprocs - number of indices for each process
1482: - indices  - indices (in local numbering) shared with neighbors (sorted by global numbering)

1484:   Level: advanced

1486: .seealso: [](sec_scatter), `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
1487:           `ISLocalToGlobalMappingGetInfo()`
1488: @*/
1489: PetscErrorCode ISLocalToGlobalMappingRestoreInfo(ISLocalToGlobalMapping mapping, PetscInt *nproc, PetscInt *procs[], PetscInt *numprocs[], PetscInt **indices[])
1490: {
1491:   PetscFunctionBegin;
1493:   if (mapping->bs > 1) {
1494:     if (numprocs) PetscCall(PetscFree(*numprocs));
1495:     if (indices) {
1496:       if (*indices)
1497:         for (PetscInt i = 0; i < *nproc; i++) PetscCall(PetscFree((*indices)[i]));
1498:       PetscCall(PetscFree(*indices));
1499:     }
1500:   }
1501:   PetscFunctionReturn(PETSC_SUCCESS);
1502: }

1504: /*@C
1505:   ISLocalToGlobalMappingGetNodeInfo - Gets the neighbor information of local nodes

1507:   Collective the first time it is called

1509:   Input Parameter:
1510: . mapping - the mapping from local to global indexing

1512:   Output Parameters:
1513: + n       - number of local nodes
1514: . n_procs - an array storing the number of processes for each local node (including self)
1515: - procs   - the processes' rank for each local node (sorted, self is first)

1517:   Level: advanced

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

1522: .seealso: [](sec_scatter), `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
1523:           `ISLocalToGlobalMappingGetInfo()`, `ISLocalToGlobalMappingRestoreNodeInfo()`, `ISLocalToGlobalMappingGetBlockNodeInfo()`
1524: @*/
1525: PetscErrorCode ISLocalToGlobalMappingGetNodeInfo(ISLocalToGlobalMapping mapping, PetscInt *n, PetscInt *n_procs[], PetscInt **procs[])
1526: {
1527:   PetscInt **bprocs = NULL, *bn_procs = NULL, bs, i, j, k, bn;

1529:   PetscFunctionBegin;
1531:   bs = mapping->bs;
1532:   PetscCall(ISLocalToGlobalMappingGetBlockNodeInfo(mapping, &bn, &bn_procs, &bprocs));
1533:   if (bs > 1) { /* we need to expand the cached info */
1534:     PetscInt *tn_procs;
1535:     PetscInt  c;

1537:     PetscCall(PetscMalloc1(bn * bs, &tn_procs));
1538:     for (i = 0, c = 0; i < bn; i++) {
1539:       for (k = 0; k < bs; k++) tn_procs[i * bs + k] = bn_procs[i];
1540:       c += bs * bn_procs[i];
1541:     }
1542:     if (n) *n = bn * bs;
1543:     if (procs) {
1544:       PetscInt **tprocs;
1545:       PetscInt   tn = bn * bs;

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

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

1570:   Not Collective

1572:   Input Parameters:
1573: + mapping - the mapping from local to global indexing
1574: . n       - number of local nodes
1575: . n_procs - an array storing the number of processes for each local node (including self)
1576: - procs   - the processes' rank for each local node (sorted, self is first)

1578:   Level: advanced

1580: .seealso: [](sec_scatter), `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
1581:           `ISLocalToGlobalMappingGetInfo()`
1582: @*/
1583: PetscErrorCode ISLocalToGlobalMappingRestoreNodeInfo(ISLocalToGlobalMapping mapping, PetscInt *n, PetscInt *n_procs[], PetscInt **procs[])
1584: {
1585:   PetscFunctionBegin;
1587:   if (mapping->bs > 1) {
1588:     if (n_procs) PetscCall(PetscFree(*n_procs));
1589:     if (procs) {
1590:       if (*procs) PetscCall(PetscFree((*procs)[0]));
1591:       PetscCall(PetscFree(*procs));
1592:     }
1593:   }
1594:   PetscCall(ISLocalToGlobalMappingRestoreBlockNodeInfo(mapping, n, n_procs, procs));
1595:   PetscFunctionReturn(PETSC_SUCCESS);
1596: }

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

1601:   Not Collective

1603:   Input Parameter:
1604: . ltog - local to global mapping

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

1609:   Level: advanced

1611:   Note:
1612:   `ISLocalToGlobalMappingGetSize()` returns the length the this array

1614: .seealso: [](sec_scatter), `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingApply()`, `ISLocalToGlobalMappingRestoreIndices()`,
1615:           `ISLocalToGlobalMappingGetBlockIndices()`, `ISLocalToGlobalMappingRestoreBlockIndices()`
1616: @*/
1617: PetscErrorCode ISLocalToGlobalMappingGetIndices(ISLocalToGlobalMapping ltog, const PetscInt *array[])
1618: {
1619:   PetscFunctionBegin;
1621:   PetscAssertPointer(array, 2);
1622:   if (ltog->bs == 1) {
1623:     *array = ltog->indices;
1624:   } else {
1625:     PetscInt       *jj, k, i, j, n = ltog->n, bs = ltog->bs;
1626:     const PetscInt *ii;

1628:     PetscCall(PetscMalloc1(bs * n, &jj));
1629:     *array = jj;
1630:     k      = 0;
1631:     ii     = ltog->indices;
1632:     for (i = 0; i < n; i++)
1633:       for (j = 0; j < bs; j++) jj[k++] = bs * ii[i] + j;
1634:   }
1635:   PetscFunctionReturn(PETSC_SUCCESS);
1636: }

1638: /*@C
1639:   ISLocalToGlobalMappingRestoreIndices - Restore indices obtained with `ISLocalToGlobalMappingGetIndices()`

1641:   Not Collective

1643:   Input Parameters:
1644: + ltog  - local to global mapping
1645: - array - array of indices

1647:   Level: advanced

1649: .seealso: [](sec_scatter), `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingApply()`, `ISLocalToGlobalMappingGetIndices()`
1650: @*/
1651: PetscErrorCode ISLocalToGlobalMappingRestoreIndices(ISLocalToGlobalMapping ltog, const PetscInt *array[])
1652: {
1653:   PetscFunctionBegin;
1655:   PetscAssertPointer(array, 2);
1656:   PetscCheck(ltog->bs != 1 || *array == ltog->indices, PETSC_COMM_SELF, PETSC_ERR_ARG_BADPTR, "Trying to return mismatched pointer");
1657:   if (ltog->bs > 1) PetscCall(PetscFree(*(void **)array));
1658:   PetscFunctionReturn(PETSC_SUCCESS);
1659: }

1661: /*@C
1662:   ISLocalToGlobalMappingGetBlockIndices - Get global indices for every local block in a `ISLocalToGlobalMapping`

1664:   Not Collective

1666:   Input Parameter:
1667: . ltog - local to global mapping

1669:   Output Parameter:
1670: . array - array of indices

1672:   Level: advanced

1674: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingApply()`,
1675:           `ISLocalToGlobalMappingRestoreBlockIndices()`
1676: @*/
1677: PetscErrorCode ISLocalToGlobalMappingGetBlockIndices(ISLocalToGlobalMapping ltog, const PetscInt *array[])
1678: {
1679:   PetscFunctionBegin;
1681:   PetscAssertPointer(array, 2);
1682:   *array = ltog->indices;
1683:   PetscFunctionReturn(PETSC_SUCCESS);
1684: }

1686: /*@C
1687:   ISLocalToGlobalMappingRestoreBlockIndices - Restore indices obtained with `ISLocalToGlobalMappingGetBlockIndices()`

1689:   Not Collective

1691:   Input Parameters:
1692: + ltog  - local to global mapping
1693: - array - array of indices

1695:   Level: advanced

1697: .seealso: [](sec_scatter), `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingApply()`, `ISLocalToGlobalMappingGetIndices()`
1698: @*/
1699: PetscErrorCode ISLocalToGlobalMappingRestoreBlockIndices(ISLocalToGlobalMapping ltog, const PetscInt *array[])
1700: {
1701:   PetscFunctionBegin;
1703:   PetscAssertPointer(array, 2);
1704:   PetscCheck(*array == ltog->indices, PETSC_COMM_SELF, PETSC_ERR_ARG_BADPTR, "Trying to return mismatched pointer");
1705:   *array = NULL;
1706:   PetscFunctionReturn(PETSC_SUCCESS);
1707: }

1709: /*@
1710:   ISLocalToGlobalMappingConcatenate - Create a new mapping that concatenates a list of mappings

1712:   Not Collective

1714:   Input Parameters:
1715: + comm  - communicator for the new mapping, must contain the communicator of every mapping to concatenate
1716: . n     - number of mappings to concatenate
1717: - ltogs - local to global mappings

1719:   Output Parameter:
1720: . ltogcat - new mapping

1722:   Level: advanced

1724:   Note:
1725:   This currently always returns a mapping with block size of 1

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

1730: .seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingCreate()`
1731: @*/
1732: PetscErrorCode ISLocalToGlobalMappingConcatenate(MPI_Comm comm, PetscInt n, const ISLocalToGlobalMapping ltogs[], ISLocalToGlobalMapping *ltogcat)
1733: {
1734:   PetscInt i, cnt, m, *idx;

1736:   PetscFunctionBegin;
1737:   PetscCheck(n >= 0, comm, PETSC_ERR_ARG_OUTOFRANGE, "Must have a non-negative number of mappings, given %" PetscInt_FMT, n);
1738:   if (n > 0) PetscAssertPointer(ltogs, 3);
1740:   PetscAssertPointer(ltogcat, 4);
1741:   for (cnt = 0, i = 0; i < n; i++) {
1742:     PetscCall(ISLocalToGlobalMappingGetSize(ltogs[i], &m));
1743:     cnt += m;
1744:   }
1745:   PetscCall(PetscMalloc1(cnt, &idx));
1746:   for (cnt = 0, i = 0; i < n; i++) {
1747:     const PetscInt *subidx;
1748:     PetscCall(ISLocalToGlobalMappingGetSize(ltogs[i], &m));
1749:     PetscCall(ISLocalToGlobalMappingGetIndices(ltogs[i], &subidx));
1750:     PetscCall(PetscArraycpy(&idx[cnt], subidx, m));
1751:     PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogs[i], &subidx));
1752:     cnt += m;
1753:   }
1754:   PetscCall(ISLocalToGlobalMappingCreate(comm, 1, cnt, idx, PETSC_OWN_POINTER, ltogcat));
1755:   PetscFunctionReturn(PETSC_SUCCESS);
1756: }

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

1762:    Options Database Key:
1763: .   -islocaltoglobalmapping_type basic - select this method

1765:    Level: beginner

1767:    Developer Note:
1768:    This stores all the mapping information on each MPI rank.

1770: .seealso: [](sec_scatter), `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingSetType()`, `ISLOCALTOGLOBALMAPPINGHASH`
1771: M*/
1772: PETSC_EXTERN PetscErrorCode ISLocalToGlobalMappingCreate_Basic(ISLocalToGlobalMapping ltog)
1773: {
1774:   PetscFunctionBegin;
1775:   ltog->ops->globaltolocalmappingapply      = ISGlobalToLocalMappingApply_Basic;
1776:   ltog->ops->globaltolocalmappingsetup      = ISGlobalToLocalMappingSetUp_Basic;
1777:   ltog->ops->globaltolocalmappingapplyblock = ISGlobalToLocalMappingApplyBlock_Basic;
1778:   ltog->ops->destroy                        = ISLocalToGlobalMappingDestroy_Basic;
1779:   PetscFunctionReturn(PETSC_SUCCESS);
1780: }

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

1786:    Options Database Key:
1787: .   -islocaltoglobalmapping_type hash - select this method

1789:    Level: beginner

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

1794: .seealso: [](sec_scatter), `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingSetType()`, `ISLOCALTOGLOBALMAPPINGBASIC`
1795: M*/
1796: PETSC_EXTERN PetscErrorCode ISLocalToGlobalMappingCreate_Hash(ISLocalToGlobalMapping ltog)
1797: {
1798:   PetscFunctionBegin;
1799:   ltog->ops->globaltolocalmappingapply      = ISGlobalToLocalMappingApply_Hash;
1800:   ltog->ops->globaltolocalmappingsetup      = ISGlobalToLocalMappingSetUp_Hash;
1801:   ltog->ops->globaltolocalmappingapplyblock = ISGlobalToLocalMappingApplyBlock_Hash;
1802:   ltog->ops->destroy                        = ISLocalToGlobalMappingDestroy_Hash;
1803:   PetscFunctionReturn(PETSC_SUCCESS);
1804: }

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

1809:   Not Collective, No Fortran Support

1811:   Input Parameters:
1812: + sname    - name of a new method
1813: - function - routine to create method context

1815:   Example Usage:
1816: .vb
1817:    ISLocalToGlobalMappingRegister("my_mapper", MyCreate);
1818: .ve

1820:   Then, your mapping can be chosen with the procedural interface via
1821: .vb
1822:   ISLocalToGlobalMappingSetType(ltog, "my_mapper")
1823: .ve
1824:   or at runtime via the option
1825: .vb
1826:   -islocaltoglobalmapping_type my_mapper
1827: .ve

1829:   Level: advanced

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

1834: .seealso: [](sec_scatter), `ISLocalToGlobalMappingRegisterAll()`, `ISLocalToGlobalMappingRegisterDestroy()`, `ISLOCALTOGLOBALMAPPINGBASIC`,
1835:           `ISLOCALTOGLOBALMAPPINGHASH`, `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingApply()`
1836: @*/
1837: PetscErrorCode ISLocalToGlobalMappingRegister(const char sname[], PetscErrorCode (*function)(ISLocalToGlobalMapping))
1838: {
1839:   PetscFunctionBegin;
1840:   PetscCall(ISInitializePackage());
1841:   PetscCall(PetscFunctionListAdd(&ISLocalToGlobalMappingList, sname, function));
1842:   PetscFunctionReturn(PETSC_SUCCESS);
1843: }

1845: /*@
1846:   ISLocalToGlobalMappingSetType - Sets the implementation type `ISLocalToGlobalMapping` will use

1848:   Logically Collective

1850:   Input Parameters:
1851: + ltog - the `ISLocalToGlobalMapping` object
1852: - type - a known method

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

1857:   Level: intermediate

1859:   Notes:
1860:   See `ISLocalToGlobalMappingType` for available methods

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

1866:   Developer Notes:
1867:   `ISLocalToGlobalMappingRegister()` is used to add new types to `ISLocalToGlobalMappingList` from which they
1868:   are accessed by `ISLocalToGlobalMappingSetType()`.

1870: .seealso: [](sec_scatter), `ISLocalToGlobalMappingType`, `ISLocalToGlobalMappingRegister()`, `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingGetType()`
1871: @*/
1872: PetscErrorCode ISLocalToGlobalMappingSetType(ISLocalToGlobalMapping ltog, ISLocalToGlobalMappingType type)
1873: {
1874:   PetscBool match;
1875:   PetscErrorCode (*r)(ISLocalToGlobalMapping) = NULL;

1877:   PetscFunctionBegin;
1879:   if (type) PetscAssertPointer(type, 2);

1881:   PetscCall(PetscObjectTypeCompare((PetscObject)ltog, type, &match));
1882:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

1884:   /* L2G maps defer type setup at globaltolocal calls, allow passing NULL here */
1885:   if (type) {
1886:     PetscCall(PetscFunctionListFind(ISLocalToGlobalMappingList, type, &r));
1887:     PetscCheck(r, PetscObjectComm((PetscObject)ltog), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unable to find requested ISLocalToGlobalMapping type %s", type);
1888:   }
1889:   /* Destroy the previous private LTOG context */
1890:   PetscTryTypeMethod(ltog, destroy);
1891:   ltog->ops->destroy = NULL;

1893:   PetscCall(PetscObjectChangeTypeName((PetscObject)ltog, type));
1894:   if (r) PetscCall((*r)(ltog));
1895:   PetscFunctionReturn(PETSC_SUCCESS);
1896: }

1898: /*@
1899:   ISLocalToGlobalMappingGetType - Get the type of the `ISLocalToGlobalMapping`

1901:   Not Collective

1903:   Input Parameter:
1904: . ltog - the `ISLocalToGlobalMapping` object

1906:   Output Parameter:
1907: . type - the type

1909:   Level: intermediate

1911: .seealso: [](sec_scatter), `ISLocalToGlobalMappingType`, `ISLocalToGlobalMappingRegister()`, `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingSetType()`
1912: @*/
1913: PetscErrorCode ISLocalToGlobalMappingGetType(ISLocalToGlobalMapping ltog, ISLocalToGlobalMappingType *type)
1914: {
1915:   PetscFunctionBegin;
1917:   PetscAssertPointer(type, 2);
1918:   *type = ((PetscObject)ltog)->type_name;
1919:   PetscFunctionReturn(PETSC_SUCCESS);
1920: }

1922: PetscBool ISLocalToGlobalMappingRegisterAllCalled = PETSC_FALSE;

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

1927:   Not Collective

1929:   Level: advanced

1931: .seealso: [](sec_scatter), `ISRegister()`, `ISLocalToGlobalRegister()`
1932: @*/
1933: PetscErrorCode ISLocalToGlobalMappingRegisterAll(void)
1934: {
1935:   PetscFunctionBegin;
1936:   if (ISLocalToGlobalMappingRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
1937:   ISLocalToGlobalMappingRegisterAllCalled = PETSC_TRUE;
1938:   PetscCall(ISLocalToGlobalMappingRegister(ISLOCALTOGLOBALMAPPINGBASIC, ISLocalToGlobalMappingCreate_Basic));
1939:   PetscCall(ISLocalToGlobalMappingRegister(ISLOCALTOGLOBALMAPPINGHASH, ISLocalToGlobalMappingCreate_Hash));
1940:   PetscFunctionReturn(PETSC_SUCCESS);
1941: }