Actual source code: index.c

  1: /*
  2:    Defines the abstract operations on index sets, i.e. the public interface.
  3: */
  4: #include <petsc/private/isimpl.h>
  5: #include <petscviewer.h>
  6: #include <petscsf.h>

  8: /* Logging support */
  9: PetscClassId IS_CLASSID;
 10: /* TODO: Much more events are missing! */
 11: PetscLogEvent IS_View;
 12: PetscLogEvent IS_Load;

 14: /*@
 15:   ISRenumber - Renumbers the non-negative entries of an index set in a contiguous way, starting from 0.

 17:   Collective

 19:   Input Parameters:
 20: + subset      - the index set
 21: - subset_mult - the multiplicity of each entry in subset (optional, can be `NULL`)

 23:   Output Parameters:
 24: + N        - one past the largest entry of the new `IS`
 25: - subset_n - the new `IS`

 27:   Level: intermediate

 29:   Note:
 30:   All negative entries are mapped to -1. Indices with non positive multiplicities are skipped.

 32: .seealso: `IS`
 33: @*/
 34: PetscErrorCode ISRenumber(IS subset, IS subset_mult, PetscInt *N, IS *subset_n)
 35: {
 36:   PetscSF         sf;
 37:   PetscLayout     map;
 38:   const PetscInt *idxs, *idxs_mult = NULL;
 39:   PetscInt       *leaf_data, *root_data, *gidxs, *ilocal, *ilocalneg;
 40:   PetscInt        N_n, n, i, lbounds[2], gbounds[2], Nl, ibs;
 41:   PetscInt        n_n, nlocals, start, first_index, npos, nneg;
 42:   PetscMPIInt     commsize;
 43:   PetscBool       first_found, isblock;

 45:   PetscFunctionBegin;
 48:   if (N) PetscAssertPointer(N, 3);
 49:   else if (!subset_n) PetscFunctionReturn(PETSC_SUCCESS);
 50:   PetscCall(ISGetLocalSize(subset, &n));
 51:   if (subset_mult) {
 52:     PetscCall(ISGetLocalSize(subset_mult, &i));
 53:     PetscCheck(i == n, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Local subset and multiplicity sizes don't match! %" PetscInt_FMT " != %" PetscInt_FMT, n, i);
 54:   }
 55:   /* create workspace layout for computing global indices of subset */
 56:   PetscCall(PetscMalloc1(n, &ilocal));
 57:   PetscCall(PetscMalloc1(n, &ilocalneg));
 58:   PetscCall(ISGetIndices(subset, &idxs));
 59:   PetscCall(ISGetBlockSize(subset, &ibs));
 60:   PetscCall(PetscObjectTypeCompare((PetscObject)subset, ISBLOCK, &isblock));
 61:   if (subset_mult) PetscCall(ISGetIndices(subset_mult, &idxs_mult));
 62:   lbounds[0] = PETSC_MAX_INT;
 63:   lbounds[1] = PETSC_MIN_INT;
 64:   for (i = 0, npos = 0, nneg = 0; i < n; i++) {
 65:     if (idxs[i] < 0) {
 66:       ilocalneg[nneg++] = i;
 67:       continue;
 68:     }
 69:     if (idxs[i] < lbounds[0]) lbounds[0] = idxs[i];
 70:     if (idxs[i] > lbounds[1]) lbounds[1] = idxs[i];
 71:     ilocal[npos++] = i;
 72:   }
 73:   if (npos == n) {
 74:     PetscCall(PetscFree(ilocal));
 75:     PetscCall(PetscFree(ilocalneg));
 76:   }

 78:   /* create sf : leaf_data == multiplicity of indexes, root data == global index in layout */
 79:   PetscCall(PetscMalloc1(n, &leaf_data));
 80:   for (i = 0; i < n; i++) leaf_data[i] = idxs_mult ? PetscMax(idxs_mult[i], 0) : 1;

 82:   /* local size of new subset */
 83:   n_n = 0;
 84:   for (i = 0; i < n; i++) n_n += leaf_data[i];
 85:   if (ilocalneg)
 86:     for (i = 0; i < nneg; i++) leaf_data[ilocalneg[i]] = 0;
 87:   PetscCall(PetscFree(ilocalneg));
 88:   PetscCall(PetscMalloc1(PetscMax(n_n, n), &gidxs)); /* allocating extra space to reuse gidxs */
 89:   /* check for early termination (all negative) */
 90:   PetscCall(PetscGlobalMinMaxInt(PetscObjectComm((PetscObject)subset), lbounds, gbounds));
 91:   if (gbounds[1] < gbounds[0]) {
 92:     if (N) *N = 0;
 93:     if (subset_n) { /* all negative */
 94:       for (i = 0; i < n_n; i++) gidxs[i] = -1;
 95:       PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)subset), n_n, gidxs, PETSC_COPY_VALUES, subset_n));
 96:     }
 97:     PetscCall(PetscFree(leaf_data));
 98:     PetscCall(PetscFree(gidxs));
 99:     PetscCall(ISRestoreIndices(subset, &idxs));
100:     if (subset_mult) PetscCall(ISRestoreIndices(subset_mult, &idxs_mult));
101:     PetscCall(PetscFree(ilocal));
102:     PetscCall(PetscFree(ilocalneg));
103:     PetscFunctionReturn(PETSC_SUCCESS);
104:   }

106:   /* split work */
107:   N_n = gbounds[1] - gbounds[0] + 1;
108:   PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)subset), &map));
109:   PetscCall(PetscLayoutSetBlockSize(map, 1));
110:   PetscCall(PetscLayoutSetSize(map, N_n));
111:   PetscCall(PetscLayoutSetUp(map));
112:   PetscCall(PetscLayoutGetLocalSize(map, &Nl));

114:   /* global indexes in layout */
115:   for (i = 0; i < npos; i++) gidxs[i] = (ilocal ? idxs[ilocal[i]] : idxs[i]) - gbounds[0];
116:   PetscCall(ISRestoreIndices(subset, &idxs));
117:   PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)subset), &sf));
118:   PetscCall(PetscSFSetGraphLayout(sf, map, npos, ilocal, PETSC_USE_POINTER, gidxs));
119:   PetscCall(PetscLayoutDestroy(&map));

121:   /* reduce from leaves to roots */
122:   PetscCall(PetscCalloc1(Nl, &root_data));
123:   PetscCall(PetscSFReduceBegin(sf, MPIU_INT, leaf_data, root_data, MPI_MAX));
124:   PetscCall(PetscSFReduceEnd(sf, MPIU_INT, leaf_data, root_data, MPI_MAX));

126:   /* count indexes in local part of layout */
127:   nlocals     = 0;
128:   first_index = -1;
129:   first_found = PETSC_FALSE;
130:   for (i = 0; i < Nl; i++) {
131:     if (!first_found && root_data[i]) {
132:       first_found = PETSC_TRUE;
133:       first_index = i;
134:     }
135:     nlocals += root_data[i];
136:   }

138:   /* cumulative of number of indexes and size of subset without holes */
139: #if defined(PETSC_HAVE_MPI_EXSCAN)
140:   start = 0;
141:   PetscCallMPI(MPI_Exscan(&nlocals, &start, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)subset)));
142: #else
143:   PetscCallMPI(MPI_Scan(&nlocals, &start, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)subset)));
144:   start = start - nlocals;
145: #endif

147:   if (N) { /* compute total size of new subset if requested */
148:     *N = start + nlocals;
149:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)subset), &commsize));
150:     PetscCallMPI(MPI_Bcast(N, 1, MPIU_INT, commsize - 1, PetscObjectComm((PetscObject)subset)));
151:   }

153:   if (!subset_n) {
154:     PetscCall(PetscFree(gidxs));
155:     PetscCall(PetscSFDestroy(&sf));
156:     PetscCall(PetscFree(leaf_data));
157:     PetscCall(PetscFree(root_data));
158:     PetscCall(PetscFree(ilocal));
159:     if (subset_mult) PetscCall(ISRestoreIndices(subset_mult, &idxs_mult));
160:     PetscFunctionReturn(PETSC_SUCCESS);
161:   }

163:   /* adapt root data with cumulative */
164:   if (first_found) {
165:     PetscInt old_index;

167:     root_data[first_index] += start;
168:     old_index = first_index;
169:     for (i = first_index + 1; i < Nl; i++) {
170:       if (root_data[i]) {
171:         root_data[i] += root_data[old_index];
172:         old_index = i;
173:       }
174:     }
175:   }

177:   /* from roots to leaves */
178:   PetscCall(PetscSFBcastBegin(sf, MPIU_INT, root_data, leaf_data, MPI_REPLACE));
179:   PetscCall(PetscSFBcastEnd(sf, MPIU_INT, root_data, leaf_data, MPI_REPLACE));
180:   PetscCall(PetscSFDestroy(&sf));

182:   /* create new IS with global indexes without holes */
183:   for (i = 0; i < n_n; i++) gidxs[i] = -1;
184:   if (subset_mult) {
185:     PetscInt cum;

187:     isblock = PETSC_FALSE;
188:     for (i = 0, cum = 0; i < n; i++)
189:       for (PetscInt j = 0; j < idxs_mult[i]; j++) gidxs[cum++] = leaf_data[i] - idxs_mult[i] + j;
190:   } else
191:     for (i = 0; i < n; i++) gidxs[i] = leaf_data[i] - 1;

193:   if (isblock) {
194:     if (ibs > 1)
195:       for (i = 0; i < n_n / ibs; i++) gidxs[i] = gidxs[i * ibs] / ibs;
196:     PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)subset), ibs, n_n / ibs, gidxs, PETSC_COPY_VALUES, subset_n));
197:   } else {
198:     PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)subset), n_n, gidxs, PETSC_COPY_VALUES, subset_n));
199:   }
200:   if (subset_mult) PetscCall(ISRestoreIndices(subset_mult, &idxs_mult));
201:   PetscCall(PetscFree(gidxs));
202:   PetscCall(PetscFree(leaf_data));
203:   PetscCall(PetscFree(root_data));
204:   PetscCall(PetscFree(ilocal));
205:   PetscFunctionReturn(PETSC_SUCCESS);
206: }

208: /*@
209:   ISCreateSubIS - Create a sub index set from a global index set selecting some components.

211:   Collective

213:   Input Parameters:
214: + is    - the index set
215: - comps - which components we will extract from `is`

217:   Output Parameters:
218: . subis - the new sub index set

220:   Example usage:
221:   We have an index set `is` living on 3 processes with the following values\:
222:   | 4 9 0 | 2 6 7 | 10 11 1|
223:   and another index set `comps` used to indicate which components of is  we want to take,
224:   | 7 5  | 1 2 | 0 4|
225:   The output index set `subis` should look like\:
226:   | 11 7 | 9 0 | 4 6|

228:   Level: intermediate

230: .seealso: `IS`, `VecGetSubVector()`, `MatCreateSubMatrix()`
231: @*/
232: PetscErrorCode ISCreateSubIS(IS is, IS comps, IS *subis)
233: {
234:   PetscSF         sf;
235:   const PetscInt *is_indices, *comps_indices;
236:   PetscInt       *subis_indices, nroots, nleaves, *mine, i, lidx;
237:   PetscMPIInt     owner;
238:   PetscSFNode    *remote;
239:   MPI_Comm        comm;

241:   PetscFunctionBegin;
244:   PetscAssertPointer(subis, 3);

246:   PetscCall(PetscObjectGetComm((PetscObject)is, &comm));
247:   PetscCall(ISGetLocalSize(comps, &nleaves));
248:   PetscCall(ISGetLocalSize(is, &nroots));
249:   PetscCall(PetscMalloc1(nleaves, &remote));
250:   PetscCall(PetscMalloc1(nleaves, &mine));
251:   PetscCall(ISGetIndices(comps, &comps_indices));
252:   /*
253:    * Construct a PetscSF in which "is" data serves as roots and "subis" is leaves.
254:    * Root data are sent to leaves using PetscSFBcast().
255:    * */
256:   for (i = 0; i < nleaves; i++) {
257:     mine[i] = i;
258:     /* Connect a remote root with the current leaf. The value on the remote root
259:      * will be received by the current local leaf.
260:      * */
261:     owner = -1;
262:     lidx  = -1;
263:     PetscCall(PetscLayoutFindOwnerIndex(is->map, comps_indices[i], &owner, &lidx));
264:     remote[i].rank  = owner;
265:     remote[i].index = lidx;
266:   }
267:   PetscCall(ISRestoreIndices(comps, &comps_indices));
268:   PetscCall(PetscSFCreate(comm, &sf));
269:   PetscCall(PetscSFSetFromOptions(sf));
270:   PetscCall(PetscSFSetGraph(sf, nroots, nleaves, mine, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER));

272:   PetscCall(PetscMalloc1(nleaves, &subis_indices));
273:   PetscCall(ISGetIndices(is, &is_indices));
274:   PetscCall(PetscSFBcastBegin(sf, MPIU_INT, is_indices, subis_indices, MPI_REPLACE));
275:   PetscCall(PetscSFBcastEnd(sf, MPIU_INT, is_indices, subis_indices, MPI_REPLACE));
276:   PetscCall(ISRestoreIndices(is, &is_indices));
277:   PetscCall(PetscSFDestroy(&sf));
278:   PetscCall(ISCreateGeneral(comm, nleaves, subis_indices, PETSC_OWN_POINTER, subis));
279:   PetscFunctionReturn(PETSC_SUCCESS);
280: }

282: /*@
283:   ISClearInfoCache - clear the cache of computed index set properties

285:   Not Collective

287:   Input Parameters:
288: + is                    - the index set
289: - clear_permanent_local - whether to remove the permanent status of local properties

291:   Level: developer

293:   Note:
294:   Because all processes must agree on the global permanent status of a property,
295:   the permanent status can only be changed with `ISSetInfo()`, because this routine is not collective

297: .seealso: `IS`, `ISInfo`, `ISInfoType`, `ISSetInfo()`
298: @*/
299: PetscErrorCode ISClearInfoCache(IS is, PetscBool clear_permanent_local)
300: {
301:   PetscInt i, j;

303:   PetscFunctionBegin;
306:   for (i = 0; i < IS_INFO_MAX; i++) {
307:     if (clear_permanent_local) is->info_permanent[IS_LOCAL][i] = PETSC_FALSE;
308:     for (j = 0; j < 2; j++) {
309:       if (!is->info_permanent[j][i]) is->info[j][i] = IS_INFO_UNKNOWN;
310:     }
311:   }
312:   PetscFunctionReturn(PETSC_SUCCESS);
313: }

315: static PetscErrorCode ISSetInfo_Internal(IS is, ISInfo info, ISInfoType type, ISInfoBool ipermanent, PetscBool flg)
316: {
317:   ISInfoBool iflg          = flg ? IS_INFO_TRUE : IS_INFO_FALSE;
318:   PetscInt   itype         = (type == IS_LOCAL) ? 0 : 1;
319:   PetscBool  permanent_set = (ipermanent == IS_INFO_UNKNOWN) ? PETSC_FALSE : PETSC_TRUE;
320:   PetscBool  permanent     = (ipermanent == IS_INFO_TRUE) ? PETSC_TRUE : PETSC_FALSE;

322:   PetscFunctionBegin;
323:   /* set this property */
324:   is->info[itype][(int)info] = iflg;
325:   if (permanent_set) is->info_permanent[itype][(int)info] = permanent;
326:   /* set implications */
327:   switch (info) {
328:   case IS_SORTED:
329:     if (flg && type == IS_GLOBAL) { /* an array that is globally sorted is also locally sorted */
330:       is->info[IS_LOCAL][(int)info] = IS_INFO_TRUE;
331:       /* global permanence implies local permanence */
332:       if (permanent_set && permanent) is->info_permanent[IS_LOCAL][(int)info] = PETSC_TRUE;
333:     }
334:     if (!flg) { /* if an array is not sorted, it cannot be an interval or the identity */
335:       is->info[itype][IS_INTERVAL] = IS_INFO_FALSE;
336:       is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
337:       if (permanent_set) {
338:         is->info_permanent[itype][IS_INTERVAL] = permanent;
339:         is->info_permanent[itype][IS_IDENTITY] = permanent;
340:       }
341:     }
342:     break;
343:   case IS_UNIQUE:
344:     if (flg && type == IS_GLOBAL) { /* an array that is globally unique is also locally unique */
345:       is->info[IS_LOCAL][(int)info] = IS_INFO_TRUE;
346:       /* global permanence implies local permanence */
347:       if (permanent_set && permanent) is->info_permanent[IS_LOCAL][(int)info] = PETSC_TRUE;
348:     }
349:     if (!flg) { /* if an array is not unique, it cannot be a permutation, and interval, or the identity */
350:       is->info[itype][IS_PERMUTATION] = IS_INFO_FALSE;
351:       is->info[itype][IS_INTERVAL]    = IS_INFO_FALSE;
352:       is->info[itype][IS_IDENTITY]    = IS_INFO_FALSE;
353:       if (permanent_set) {
354:         is->info_permanent[itype][IS_PERMUTATION] = permanent;
355:         is->info_permanent[itype][IS_INTERVAL]    = permanent;
356:         is->info_permanent[itype][IS_IDENTITY]    = permanent;
357:       }
358:     }
359:     break;
360:   case IS_PERMUTATION:
361:     if (flg) { /* an array that is a permutation is unique and is unique locally */
362:       is->info[itype][IS_UNIQUE]    = IS_INFO_TRUE;
363:       is->info[IS_LOCAL][IS_UNIQUE] = IS_INFO_TRUE;
364:       if (permanent_set && permanent) {
365:         is->info_permanent[itype][IS_UNIQUE]    = PETSC_TRUE;
366:         is->info_permanent[IS_LOCAL][IS_UNIQUE] = PETSC_TRUE;
367:       }
368:     } else { /* an array that is not a permutation cannot be the identity */
369:       is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
370:       if (permanent_set) is->info_permanent[itype][IS_IDENTITY] = permanent;
371:     }
372:     break;
373:   case IS_INTERVAL:
374:     if (flg) { /* an array that is an interval is sorted and unique */
375:       is->info[itype][IS_SORTED]    = IS_INFO_TRUE;
376:       is->info[IS_LOCAL][IS_SORTED] = IS_INFO_TRUE;
377:       is->info[itype][IS_UNIQUE]    = IS_INFO_TRUE;
378:       is->info[IS_LOCAL][IS_UNIQUE] = IS_INFO_TRUE;
379:       if (permanent_set && permanent) {
380:         is->info_permanent[itype][IS_SORTED]    = PETSC_TRUE;
381:         is->info_permanent[IS_LOCAL][IS_SORTED] = PETSC_TRUE;
382:         is->info_permanent[itype][IS_UNIQUE]    = PETSC_TRUE;
383:         is->info_permanent[IS_LOCAL][IS_UNIQUE] = PETSC_TRUE;
384:       }
385:     } else { /* an array that is not an interval cannot be the identity */
386:       is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
387:       if (permanent_set) is->info_permanent[itype][IS_IDENTITY] = permanent;
388:     }
389:     break;
390:   case IS_IDENTITY:
391:     if (flg) { /* an array that is the identity is sorted, unique, an interval, and a permutation */
392:       is->info[itype][IS_SORTED]      = IS_INFO_TRUE;
393:       is->info[IS_LOCAL][IS_SORTED]   = IS_INFO_TRUE;
394:       is->info[itype][IS_UNIQUE]      = IS_INFO_TRUE;
395:       is->info[IS_LOCAL][IS_UNIQUE]   = IS_INFO_TRUE;
396:       is->info[itype][IS_PERMUTATION] = IS_INFO_TRUE;
397:       is->info[itype][IS_INTERVAL]    = IS_INFO_TRUE;
398:       is->info[IS_LOCAL][IS_INTERVAL] = IS_INFO_TRUE;
399:       if (permanent_set && permanent) {
400:         is->info_permanent[itype][IS_SORTED]      = PETSC_TRUE;
401:         is->info_permanent[IS_LOCAL][IS_SORTED]   = PETSC_TRUE;
402:         is->info_permanent[itype][IS_UNIQUE]      = PETSC_TRUE;
403:         is->info_permanent[IS_LOCAL][IS_UNIQUE]   = PETSC_TRUE;
404:         is->info_permanent[itype][IS_PERMUTATION] = PETSC_TRUE;
405:         is->info_permanent[itype][IS_INTERVAL]    = PETSC_TRUE;
406:         is->info_permanent[IS_LOCAL][IS_INTERVAL] = PETSC_TRUE;
407:       }
408:     }
409:     break;
410:   default:
411:     PetscCheck(type != IS_LOCAL, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown IS property");
412:     SETERRQ(PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Unknown IS property");
413:   }
414:   PetscFunctionReturn(PETSC_SUCCESS);
415: }

417: // PetscClangLinter pragma disable: -fdoc-section-header-unknown
418: /*@
419:   ISSetInfo - Set known information about an index set.

421:   Logically Collective if `ISInfoType` is `IS_GLOBAL`

423:   Input Parameters:
424: + is        - the index set
425: . info      - describing a property of the index set, one of those listed below,
426: . type      - `IS_LOCAL` if the information describes the local portion of the index set,
427:           `IS_GLOBAL` if it describes the whole index set
428: . permanent - `PETSC_TRUE` if it is known that the property will persist through changes to the index set, `PETSC_FALSE` otherwise
429:                If the user sets a property as permanently known, it will bypass computation of that property
430: - flg       - set the described property as true (`PETSC_TRUE`) or false (`PETSC_FALSE`)

432:   Values of `info` Describing `IS` Structure:
433: + `IS_SORTED`      - the [local part of the] index set is sorted in ascending order
434: . `IS_UNIQUE`      - each entry in the [local part of the] index set is unique
435: . `IS_PERMUTATION` - the [local part of the] index set is a permutation of the integers {0, 1, ..., N-1}, where N is the size of the [local part of the] index set
436: . `IS_INTERVAL`    - the [local part of the] index set is equal to a contiguous range of integers {f, f + 1, ..., f + N-1}
437: - `IS_IDENTITY`    - the [local part of the] index set is equal to the integers {0, 1, ..., N-1}

439:   Level: advanced

441:   Notes:
442:   If type is `IS_GLOBAL`, all processes that share the index set must pass the same value in flg

444:   It is possible to set a property with `ISSetInfo()` that contradicts what would be previously computed with `ISGetInfo()`

446: .seealso: `ISInfo`, `ISInfoType`, `IS`
447: @*/
448: PetscErrorCode ISSetInfo(IS is, ISInfo info, ISInfoType type, PetscBool permanent, PetscBool flg)
449: {
450:   MPI_Comm    comm, errcomm;
451:   PetscMPIInt size;

453:   PetscFunctionBegin;
456:   comm = PetscObjectComm((PetscObject)is);
457:   if (type == IS_GLOBAL) {
461:     errcomm = comm;
462:   } else {
463:     errcomm = PETSC_COMM_SELF;
464:   }

466:   PetscCheck(((int)info) > IS_INFO_MIN && ((int)info) < IS_INFO_MAX, errcomm, PETSC_ERR_ARG_OUTOFRANGE, "Options %d is out of range", (int)info);

468:   PetscCallMPI(MPI_Comm_size(comm, &size));
469:   /* do not use global values if size == 1: it makes it easier to keep the implications straight */
470:   if (size == 1) type = IS_LOCAL;
471:   PetscCall(ISSetInfo_Internal(is, info, type, permanent ? IS_INFO_TRUE : IS_INFO_FALSE, flg));
472:   PetscFunctionReturn(PETSC_SUCCESS);
473: }

475: static PetscErrorCode ISGetInfo_Sorted_Private(IS is, ISInfoType type, PetscBool *flg)
476: {
477:   MPI_Comm    comm;
478:   PetscMPIInt size, rank;

480:   PetscFunctionBegin;
481:   comm = PetscObjectComm((PetscObject)is);
482:   PetscCallMPI(MPI_Comm_size(comm, &size));
483:   PetscCallMPI(MPI_Comm_size(comm, &rank));
484:   if (type == IS_GLOBAL && is->ops->sortedglobal) {
485:     PetscUseTypeMethod(is, sortedglobal, flg);
486:   } else {
487:     PetscBool sortedLocal = PETSC_FALSE;

489:     /* determine if the array is locally sorted */
490:     if (type == IS_GLOBAL && size > 1) {
491:       /* call ISGetInfo so that a cached value will be used if possible */
492:       PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, &sortedLocal));
493:     } else if (is->ops->sortedlocal) {
494:       PetscUseTypeMethod(is, sortedlocal, &sortedLocal);
495:     } else {
496:       /* default: get the local indices and directly check */
497:       const PetscInt *idx;
498:       PetscInt        n;

500:       PetscCall(ISGetIndices(is, &idx));
501:       PetscCall(ISGetLocalSize(is, &n));
502:       PetscCall(PetscSortedInt(n, idx, &sortedLocal));
503:       PetscCall(ISRestoreIndices(is, &idx));
504:     }

506:     if (type == IS_LOCAL || size == 1) {
507:       *flg = sortedLocal;
508:     } else {
509:       PetscCall(MPIU_Allreduce(&sortedLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm));
510:       if (*flg) {
511:         PetscInt n, min = PETSC_MAX_INT, max = PETSC_MIN_INT;
512:         PetscInt maxprev;

514:         PetscCall(ISGetLocalSize(is, &n));
515:         if (n) PetscCall(ISGetMinMax(is, &min, &max));
516:         maxprev = PETSC_MIN_INT;
517:         PetscCallMPI(MPI_Exscan(&max, &maxprev, 1, MPIU_INT, MPI_MAX, comm));
518:         if (rank && (maxprev > min)) sortedLocal = PETSC_FALSE;
519:         PetscCall(MPIU_Allreduce(&sortedLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm));
520:       }
521:     }
522:   }
523:   PetscFunctionReturn(PETSC_SUCCESS);
524: }

526: static PetscErrorCode ISGetIndicesCopy_Private(IS is, PetscInt idx[]);

528: static PetscErrorCode ISGetInfo_Unique_Private(IS is, ISInfoType type, PetscBool *flg)
529: {
530:   MPI_Comm    comm;
531:   PetscMPIInt size, rank;
532:   PetscInt    i;

534:   PetscFunctionBegin;
535:   comm = PetscObjectComm((PetscObject)is);
536:   PetscCallMPI(MPI_Comm_size(comm, &size));
537:   PetscCallMPI(MPI_Comm_size(comm, &rank));
538:   if (type == IS_GLOBAL && is->ops->uniqueglobal) {
539:     PetscUseTypeMethod(is, uniqueglobal, flg);
540:   } else {
541:     PetscBool uniqueLocal;
542:     PetscInt  n   = -1;
543:     PetscInt *idx = NULL;

545:     /* determine if the array is locally unique */
546:     if (type == IS_GLOBAL && size > 1) {
547:       /* call ISGetInfo so that a cached value will be used if possible */
548:       PetscCall(ISGetInfo(is, IS_UNIQUE, IS_LOCAL, PETSC_TRUE, &uniqueLocal));
549:     } else if (is->ops->uniquelocal) {
550:       PetscUseTypeMethod(is, uniquelocal, &uniqueLocal);
551:     } else {
552:       /* default: get the local indices and directly check */
553:       uniqueLocal = PETSC_TRUE;
554:       PetscCall(ISGetLocalSize(is, &n));
555:       PetscCall(PetscMalloc1(n, &idx));
556:       PetscCall(ISGetIndicesCopy_Private(is, idx));
557:       PetscCall(PetscIntSortSemiOrdered(n, idx));
558:       for (i = 1; i < n; i++)
559:         if (idx[i] == idx[i - 1]) break;
560:       if (i < n) uniqueLocal = PETSC_FALSE;
561:     }

563:     PetscCall(PetscFree(idx));
564:     if (type == IS_LOCAL || size == 1) {
565:       *flg = uniqueLocal;
566:     } else {
567:       PetscCall(MPIU_Allreduce(&uniqueLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm));
568:       if (*flg) {
569:         PetscInt min = PETSC_MAX_INT, max = PETSC_MIN_INT, maxprev;

571:         if (!idx) {
572:           PetscCall(ISGetLocalSize(is, &n));
573:           PetscCall(PetscMalloc1(n, &idx));
574:           PetscCall(ISGetIndicesCopy_Private(is, idx));
575:         }
576:         PetscCall(PetscParallelSortInt(is->map, is->map, idx, idx));
577:         if (n) {
578:           min = idx[0];
579:           max = idx[n - 1];
580:         }
581:         for (i = 1; i < n; i++)
582:           if (idx[i] == idx[i - 1]) break;
583:         if (i < n) uniqueLocal = PETSC_FALSE;
584:         maxprev = PETSC_MIN_INT;
585:         PetscCallMPI(MPI_Exscan(&max, &maxprev, 1, MPIU_INT, MPI_MAX, comm));
586:         if (rank && (maxprev == min)) uniqueLocal = PETSC_FALSE;
587:         PetscCall(MPIU_Allreduce(&uniqueLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm));
588:       }
589:     }
590:     PetscCall(PetscFree(idx));
591:   }
592:   PetscFunctionReturn(PETSC_SUCCESS);
593: }

595: static PetscErrorCode ISGetInfo_Permutation(IS is, ISInfoType type, PetscBool *flg)
596: {
597:   MPI_Comm    comm;
598:   PetscMPIInt size, rank;

600:   PetscFunctionBegin;
601:   comm = PetscObjectComm((PetscObject)is);
602:   PetscCallMPI(MPI_Comm_size(comm, &size));
603:   PetscCallMPI(MPI_Comm_size(comm, &rank));
604:   if (type == IS_GLOBAL && is->ops->permglobal) {
605:     PetscUseTypeMethod(is, permglobal, flg);
606:   } else if (type == IS_LOCAL && is->ops->permlocal) {
607:     PetscUseTypeMethod(is, permlocal, flg);
608:   } else {
609:     PetscBool permLocal;
610:     PetscInt  n, i, rStart;
611:     PetscInt *idx;

613:     PetscCall(ISGetLocalSize(is, &n));
614:     PetscCall(PetscMalloc1(n, &idx));
615:     PetscCall(ISGetIndicesCopy_Private(is, idx));
616:     if (type == IS_GLOBAL) {
617:       PetscCall(PetscParallelSortInt(is->map, is->map, idx, idx));
618:       PetscCall(PetscLayoutGetRange(is->map, &rStart, NULL));
619:     } else {
620:       PetscCall(PetscIntSortSemiOrdered(n, idx));
621:       rStart = 0;
622:     }
623:     permLocal = PETSC_TRUE;
624:     for (i = 0; i < n; i++) {
625:       if (idx[i] != rStart + i) break;
626:     }
627:     if (i < n) permLocal = PETSC_FALSE;
628:     if (type == IS_LOCAL || size == 1) {
629:       *flg = permLocal;
630:     } else {
631:       PetscCall(MPIU_Allreduce(&permLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm));
632:     }
633:     PetscCall(PetscFree(idx));
634:   }
635:   PetscFunctionReturn(PETSC_SUCCESS);
636: }

638: static PetscErrorCode ISGetInfo_Interval(IS is, ISInfoType type, PetscBool *flg)
639: {
640:   MPI_Comm    comm;
641:   PetscMPIInt size, rank;
642:   PetscInt    i;

644:   PetscFunctionBegin;
645:   comm = PetscObjectComm((PetscObject)is);
646:   PetscCallMPI(MPI_Comm_size(comm, &size));
647:   PetscCallMPI(MPI_Comm_size(comm, &rank));
648:   if (type == IS_GLOBAL && is->ops->intervalglobal) {
649:     PetscUseTypeMethod(is, intervalglobal, flg);
650:   } else {
651:     PetscBool intervalLocal;

653:     /* determine if the array is locally an interval */
654:     if (type == IS_GLOBAL && size > 1) {
655:       /* call ISGetInfo so that a cached value will be used if possible */
656:       PetscCall(ISGetInfo(is, IS_INTERVAL, IS_LOCAL, PETSC_TRUE, &intervalLocal));
657:     } else if (is->ops->intervallocal) {
658:       PetscUseTypeMethod(is, intervallocal, &intervalLocal);
659:     } else {
660:       PetscInt        n;
661:       const PetscInt *idx;
662:       /* default: get the local indices and directly check */
663:       intervalLocal = PETSC_TRUE;
664:       PetscCall(ISGetLocalSize(is, &n));
665:       PetscCall(ISGetIndices(is, &idx));
666:       for (i = 1; i < n; i++)
667:         if (idx[i] != idx[i - 1] + 1) break;
668:       if (i < n) intervalLocal = PETSC_FALSE;
669:       PetscCall(ISRestoreIndices(is, &idx));
670:     }

672:     if (type == IS_LOCAL || size == 1) {
673:       *flg = intervalLocal;
674:     } else {
675:       PetscCall(MPIU_Allreduce(&intervalLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm));
676:       if (*flg) {
677:         PetscInt n, min = PETSC_MAX_INT, max = PETSC_MIN_INT;
678:         PetscInt maxprev;

680:         PetscCall(ISGetLocalSize(is, &n));
681:         if (n) PetscCall(ISGetMinMax(is, &min, &max));
682:         maxprev = PETSC_MIN_INT;
683:         PetscCallMPI(MPI_Exscan(&max, &maxprev, 1, MPIU_INT, MPI_MAX, comm));
684:         if (rank && n && (maxprev != min - 1)) intervalLocal = PETSC_FALSE;
685:         PetscCall(MPIU_Allreduce(&intervalLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm));
686:       }
687:     }
688:   }
689:   PetscFunctionReturn(PETSC_SUCCESS);
690: }

692: static PetscErrorCode ISGetInfo_Identity(IS is, ISInfoType type, PetscBool *flg)
693: {
694:   MPI_Comm    comm;
695:   PetscMPIInt size, rank;

697:   PetscFunctionBegin;
698:   comm = PetscObjectComm((PetscObject)is);
699:   PetscCallMPI(MPI_Comm_size(comm, &size));
700:   PetscCallMPI(MPI_Comm_size(comm, &rank));
701:   if (type == IS_GLOBAL && is->ops->intervalglobal) {
702:     PetscBool isinterval;

704:     PetscUseTypeMethod(is, intervalglobal, &isinterval);
705:     *flg = PETSC_FALSE;
706:     if (isinterval) {
707:       PetscInt min;

709:       PetscCall(ISGetMinMax(is, &min, NULL));
710:       PetscCallMPI(MPI_Bcast(&min, 1, MPIU_INT, 0, comm));
711:       if (min == 0) *flg = PETSC_TRUE;
712:     }
713:   } else if (type == IS_LOCAL && is->ops->intervallocal) {
714:     PetscBool isinterval;

716:     PetscUseTypeMethod(is, intervallocal, &isinterval);
717:     *flg = PETSC_FALSE;
718:     if (isinterval) {
719:       PetscInt min;

721:       PetscCall(ISGetMinMax(is, &min, NULL));
722:       if (min == 0) *flg = PETSC_TRUE;
723:     }
724:   } else {
725:     PetscBool       identLocal;
726:     PetscInt        n, i, rStart;
727:     const PetscInt *idx;

729:     PetscCall(ISGetLocalSize(is, &n));
730:     PetscCall(ISGetIndices(is, &idx));
731:     PetscCall(PetscLayoutGetRange(is->map, &rStart, NULL));
732:     identLocal = PETSC_TRUE;
733:     for (i = 0; i < n; i++) {
734:       if (idx[i] != rStart + i) break;
735:     }
736:     if (i < n) identLocal = PETSC_FALSE;
737:     if (type == IS_LOCAL || size == 1) {
738:       *flg = identLocal;
739:     } else {
740:       PetscCall(MPIU_Allreduce(&identLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm));
741:     }
742:     PetscCall(ISRestoreIndices(is, &idx));
743:   }
744:   PetscFunctionReturn(PETSC_SUCCESS);
745: }

747: /*@
748:   ISGetInfo - Determine whether an index set satisfies a given property

750:   Collective or Logically Collective if the type is `IS_GLOBAL` (logically collective if the value of the property has been permanently set with `ISSetInfo()`)

752:   Input Parameters:
753: + is      - the index set
754: . info    - describing a property of the index set, one of those listed in the documentation of `ISSetInfo()`
755: . compute - if `PETSC_FALSE`, the property will not be computed if it is not already known and the property will be assumed to be false
756: - type    - whether the property is local (`IS_LOCAL`) or global (`IS_GLOBAL`)

758:   Output Parameter:
759: . flg - whether the property is true (`PETSC_TRUE`) or false (`PETSC_FALSE`)

761:   Level: advanced

763:   Notes:
764:   `ISGetInfo()` uses cached values when possible, which will be incorrect if `ISSetInfo()` has been called with incorrect information.

766:   To clear cached values, use `ISClearInfoCache()`.

768: .seealso: `IS`, `ISInfo`, `ISInfoType`, `ISSetInfo()`, `ISClearInfoCache()`
769: @*/
770: PetscErrorCode ISGetInfo(IS is, ISInfo info, ISInfoType type, PetscBool compute, PetscBool *flg)
771: {
772:   MPI_Comm    comm, errcomm;
773:   PetscMPIInt rank, size;
774:   PetscInt    itype;
775:   PetscBool   hasprop;
776:   PetscBool   infer;

778:   PetscFunctionBegin;
781:   comm = PetscObjectComm((PetscObject)is);
782:   if (type == IS_GLOBAL) {
784:     errcomm = comm;
785:   } else {
786:     errcomm = PETSC_COMM_SELF;
787:   }

789:   PetscCallMPI(MPI_Comm_size(comm, &size));
790:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

792:   PetscCheck(((int)info) > IS_INFO_MIN && ((int)info) < IS_INFO_MAX, errcomm, PETSC_ERR_ARG_OUTOFRANGE, "Options %d is out of range", (int)info);
793:   if (size == 1) type = IS_LOCAL;
794:   itype   = (type == IS_LOCAL) ? 0 : 1;
795:   hasprop = PETSC_FALSE;
796:   infer   = PETSC_FALSE;
797:   if (is->info_permanent[itype][(int)info]) {
798:     hasprop = (is->info[itype][(int)info] == IS_INFO_TRUE) ? PETSC_TRUE : PETSC_FALSE;
799:     infer   = PETSC_TRUE;
800:   } else if ((itype == IS_LOCAL) && (is->info[IS_LOCAL][info] != IS_INFO_UNKNOWN)) {
801:     /* we can cache local properties as long as we clear them when the IS changes */
802:     /* NOTE: we only cache local values because there is no ISAssemblyBegin()/ISAssemblyEnd(),
803:      so we have no way of knowing when a cached value has been invalidated by changes on a different process */
804:     hasprop = (is->info[itype][(int)info] == IS_INFO_TRUE) ? PETSC_TRUE : PETSC_FALSE;
805:     infer   = PETSC_TRUE;
806:   } else if (compute) {
807:     switch (info) {
808:     case IS_SORTED:
809:       PetscCall(ISGetInfo_Sorted_Private(is, type, &hasprop));
810:       break;
811:     case IS_UNIQUE:
812:       PetscCall(ISGetInfo_Unique_Private(is, type, &hasprop));
813:       break;
814:     case IS_PERMUTATION:
815:       PetscCall(ISGetInfo_Permutation(is, type, &hasprop));
816:       break;
817:     case IS_INTERVAL:
818:       PetscCall(ISGetInfo_Interval(is, type, &hasprop));
819:       break;
820:     case IS_IDENTITY:
821:       PetscCall(ISGetInfo_Identity(is, type, &hasprop));
822:       break;
823:     default:
824:       SETERRQ(errcomm, PETSC_ERR_ARG_OUTOFRANGE, "Unknown IS property");
825:     }
826:     infer = PETSC_TRUE;
827:   }
828:   /* call ISSetInfo_Internal to keep all of the implications straight */
829:   if (infer) PetscCall(ISSetInfo_Internal(is, info, type, IS_INFO_UNKNOWN, hasprop));
830:   *flg = hasprop;
831:   PetscFunctionReturn(PETSC_SUCCESS);
832: }

834: static PetscErrorCode ISCopyInfo_Private(IS source, IS dest)
835: {
836:   PetscFunctionBegin;
837:   PetscCall(PetscArraycpy(&dest->info[0], &source->info[0], 2));
838:   PetscCall(PetscArraycpy(&dest->info_permanent[0], &source->info_permanent[0], 2));
839:   PetscFunctionReturn(PETSC_SUCCESS);
840: }

842: /*@
843:   ISIdentity - Determines whether index set is the identity mapping.

845:   Collective

847:   Input Parameter:
848: . is - the index set

850:   Output Parameter:
851: . ident - `PETSC_TRUE` if an identity, else `PETSC_FALSE`

853:   Level: intermediate

855:   Note:
856:   If `ISSetIdentity()` (or `ISSetInfo()` for a permanent property) has been called,
857:   `ISIdentity()` will return its answer without communication between processes, but
858:   otherwise the output ident will be computed from `ISGetInfo()`,
859:   which may require synchronization on the communicator of `is`.  To avoid this computation,
860:   call `ISGetInfo()` directly with the compute flag set to `PETSC_FALSE`, and ident will be assumed false.

862: .seealso: `IS`, `ISSetIdentity()`, `ISGetInfo()`
863: @*/
864: PetscErrorCode ISIdentity(IS is, PetscBool *ident)
865: {
866:   PetscFunctionBegin;
868:   PetscAssertPointer(ident, 2);
869:   PetscCall(ISGetInfo(is, IS_IDENTITY, IS_GLOBAL, PETSC_TRUE, ident));
870:   PetscFunctionReturn(PETSC_SUCCESS);
871: }

873: /*@
874:   ISSetIdentity - Informs the index set that it is an identity.

876:   Logically Collective

878:   Input Parameter:
879: . is - the index set

881:   Level: intermediate

883:   Notes:
884:   `is` will be considered the identity permanently, even if indices have been changes (for example, with
885:   `ISGeneralSetIndices()`).  It's a good idea to only set this property if `is` will not change in the future.

887:   To clear this property, use `ISClearInfoCache()`.

889:   Developer Notes:
890:   Some of these info routines have statements about values changing in the `IS`, this seems to contradict the fact that `IS` cannot be changed?

892: .seealso: `IS`, `ISIdentity()`, `ISSetInfo()`, `ISClearInfoCache()`
893: @*/
894: PetscErrorCode ISSetIdentity(IS is)
895: {
896:   PetscFunctionBegin;
898:   PetscCall(ISSetInfo(is, IS_IDENTITY, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE));
899:   PetscFunctionReturn(PETSC_SUCCESS);
900: }

902: /*@
903:   ISContiguousLocal - Locates an index set with contiguous range within a global range, if possible

905:   Not Collective

907:   Input Parameters:
908: + is     - the index set
909: . gstart - global start
910: - gend   - global end

912:   Output Parameters:
913: + start  - start of contiguous block, as an offset from `gstart`
914: - contig - `PETSC_TRUE` if the index set refers to contiguous entries on this process, else `PETSC_FALSE`

916:   Level: developer

918: .seealso: `IS`, `ISGetLocalSize()`, `VecGetOwnershipRange()`
919: @*/
920: PetscErrorCode ISContiguousLocal(IS is, PetscInt gstart, PetscInt gend, PetscInt *start, PetscBool *contig)
921: {
922:   PetscFunctionBegin;
924:   PetscAssertPointer(start, 4);
925:   PetscAssertPointer(contig, 5);
926:   *start  = -1;
927:   *contig = PETSC_FALSE;
928:   PetscTryTypeMethod(is, contiguous, gstart, gend, start, contig);
929:   PetscFunctionReturn(PETSC_SUCCESS);
930: }

932: /*@
933:   ISPermutation - `PETSC_TRUE` or `PETSC_FALSE` depending on whether the
934:   index set has been declared to be a permutation.

936:   Logically Collective

938:   Input Parameter:
939: . is - the index set

941:   Output Parameter:
942: . perm - `PETSC_TRUE` if a permutation, else `PETSC_FALSE`

944:   Level: intermediate

946:   Note:
947:   If it is not already known that `is` is a permutation (if `ISSetPermutation()`
948:   or `ISSetInfo()` has not been called), this routine will not attempt to compute
949:   whether the index set is a permutation and will assume `perm` is `PETSC_FALSE`.
950:   To compute the value when it is not already known, use `ISGetInfo()` with
951:   the compute flag set to `PETSC_TRUE`.

953:   Developer Notes:
954:   Perhaps some of these routines should use the `PetscBool3` enum to return appropriate values

956: .seealso: `IS`, `ISSetPermutation()`, `ISGetInfo()`
957: @*/
958: PetscErrorCode ISPermutation(IS is, PetscBool *perm)
959: {
960:   PetscFunctionBegin;
962:   PetscAssertPointer(perm, 2);
963:   PetscCall(ISGetInfo(is, IS_PERMUTATION, IS_GLOBAL, PETSC_FALSE, perm));
964:   PetscFunctionReturn(PETSC_SUCCESS);
965: }

967: /*@
968:   ISSetPermutation - Informs the index set that it is a permutation.

970:   Logically Collective

972:   Input Parameter:
973: . is - the index set

975:   Level: intermediate

977:   Notes:
978:   `is` will be considered a permutation permanently, even if indices have been changes (for example, with
979:   `ISGeneralSetIndices()`).  It's a good idea to only set this property if `is` will not change in the future.

981:   To clear this property, use `ISClearInfoCache()`.

983:   The debug version of the libraries (./configure --with-debugging=1) checks if the
984:   index set is actually a permutation. The optimized version just believes you.

986: .seealso: `IS`, `ISPermutation()`, `ISSetInfo()`, `ISClearInfoCache().`
987: @*/
988: PetscErrorCode ISSetPermutation(IS is)
989: {
990:   PetscFunctionBegin;
992:   if (PetscDefined(USE_DEBUG)) {
993:     PetscMPIInt size;

995:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is), &size));
996:     if (size == 1) {
997:       PetscInt        i, n, *idx;
998:       const PetscInt *iidx;

1000:       PetscCall(ISGetSize(is, &n));
1001:       PetscCall(PetscMalloc1(n, &idx));
1002:       PetscCall(ISGetIndices(is, &iidx));
1003:       PetscCall(PetscArraycpy(idx, iidx, n));
1004:       PetscCall(PetscIntSortSemiOrdered(n, idx));
1005:       for (i = 0; i < n; i++) PetscCheck(idx[i] == i, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Index set is not a permutation");
1006:       PetscCall(PetscFree(idx));
1007:       PetscCall(ISRestoreIndices(is, &iidx));
1008:     }
1009:   }
1010:   PetscCall(ISSetInfo(is, IS_PERMUTATION, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE));
1011:   PetscFunctionReturn(PETSC_SUCCESS);
1012: }

1014: /*@C
1015:   ISDestroy - Destroys an index set.

1017:   Collective

1019:   Input Parameter:
1020: . is - the index set

1022:   Level: beginner

1024: .seealso: `IS`, `ISCreateGeneral()`, `ISCreateStride()`, `ISCreateBlocked()`
1025: @*/
1026: PetscErrorCode ISDestroy(IS *is)
1027: {
1028:   PetscFunctionBegin;
1029:   if (!*is) PetscFunctionReturn(PETSC_SUCCESS);
1031:   if (--((PetscObject)(*is))->refct > 0) {
1032:     *is = NULL;
1033:     PetscFunctionReturn(PETSC_SUCCESS);
1034:   }
1035:   if ((*is)->complement) {
1036:     PetscInt refcnt;
1037:     PetscCall(PetscObjectGetReference((PetscObject)((*is)->complement), &refcnt));
1038:     PetscCheck(refcnt <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Nonlocal IS has not been restored");
1039:     PetscCall(ISDestroy(&(*is)->complement));
1040:   }
1041:   if ((*is)->ops->destroy) PetscCall((*(*is)->ops->destroy)(*is));
1042:   PetscCall(PetscLayoutDestroy(&(*is)->map));
1043:   /* Destroy local representations of offproc data. */
1044:   PetscCall(PetscFree((*is)->total));
1045:   PetscCall(PetscFree((*is)->nonlocal));
1046:   PetscCall(PetscHeaderDestroy(is));
1047:   PetscFunctionReturn(PETSC_SUCCESS);
1048: }

1050: /*@
1051:   ISInvertPermutation - Creates a new permutation that is the inverse of
1052:   a given permutation.

1054:   Collective

1056:   Input Parameters:
1057: + is     - the index set
1058: - nlocal - number of indices on this processor in result (ignored for 1 processor) or
1059:             use `PETSC_DECIDE`

1061:   Output Parameter:
1062: . isout - the inverse permutation

1064:   Level: intermediate

1066:   Note:
1067:   For parallel index sets this does the complete parallel permutation, but the
1068:   code is not efficient for huge index sets (10,000,000 indices).

1070: .seealso: `IS`, `ISGetInfo()`, `ISSetPermutation()`, `ISGetPermutation()`
1071: @*/
1072: PetscErrorCode ISInvertPermutation(IS is, PetscInt nlocal, IS *isout)
1073: {
1074:   PetscBool isperm, isidentity, issame;

1076:   PetscFunctionBegin;
1078:   PetscAssertPointer(isout, 3);
1079:   PetscCall(ISGetInfo(is, IS_PERMUTATION, IS_GLOBAL, PETSC_TRUE, &isperm));
1080:   PetscCheck(isperm, PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_WRONG, "Not a permutation");
1081:   PetscCall(ISGetInfo(is, IS_IDENTITY, IS_GLOBAL, PETSC_TRUE, &isidentity));
1082:   issame = PETSC_FALSE;
1083:   if (isidentity) {
1084:     PetscInt  n;
1085:     PetscBool isallsame;

1087:     PetscCall(ISGetLocalSize(is, &n));
1088:     issame = (PetscBool)(n == nlocal);
1089:     PetscCall(MPIU_Allreduce(&issame, &isallsame, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1090:     issame = isallsame;
1091:   }
1092:   if (issame) {
1093:     PetscCall(ISDuplicate(is, isout));
1094:   } else {
1095:     PetscUseTypeMethod(is, invertpermutation, nlocal, isout);
1096:     PetscCall(ISSetPermutation(*isout));
1097:   }
1098:   PetscFunctionReturn(PETSC_SUCCESS);
1099: }

1101: /*@
1102:   ISGetSize - Returns the global length of an index set.

1104:   Not Collective

1106:   Input Parameter:
1107: . is - the index set

1109:   Output Parameter:
1110: . size - the global size

1112:   Level: beginner

1114: .seealso: `IS`
1115: @*/
1116: PetscErrorCode ISGetSize(IS is, PetscInt *size)
1117: {
1118:   PetscFunctionBegin;
1120:   PetscAssertPointer(size, 2);
1121:   *size = is->map->N;
1122:   PetscFunctionReturn(PETSC_SUCCESS);
1123: }

1125: /*@
1126:   ISGetLocalSize - Returns the local (processor) length of an index set.

1128:   Not Collective

1130:   Input Parameter:
1131: . is - the index set

1133:   Output Parameter:
1134: . size - the local size

1136:   Level: beginner

1138: .seealso: `IS`, `ISGetSize()`
1139: @*/
1140: PetscErrorCode ISGetLocalSize(IS is, PetscInt *size)
1141: {
1142:   PetscFunctionBegin;
1144:   PetscAssertPointer(size, 2);
1145:   *size = is->map->n;
1146:   PetscFunctionReturn(PETSC_SUCCESS);
1147: }

1149: /*@
1150:   ISGetLayout - get `PetscLayout` describing index set layout

1152:   Not Collective

1154:   Input Parameter:
1155: . is - the index set

1157:   Output Parameter:
1158: . map - the layout

1160:   Level: developer

1162: .seealso: `IS`, `PetscLayout`, `ISSetLayout()`, `ISGetSize()`, `ISGetLocalSize()`
1163: @*/
1164: PetscErrorCode ISGetLayout(IS is, PetscLayout *map)
1165: {
1166:   PetscFunctionBegin;
1168:   PetscAssertPointer(map, 2);
1169:   *map = is->map;
1170:   PetscFunctionReturn(PETSC_SUCCESS);
1171: }

1173: /*@
1174:   ISSetLayout - set `PetscLayout` describing index set layout

1176:   Collective

1178:   Input Parameters:
1179: + is  - the index set
1180: - map - the layout

1182:   Level: developer

1184:   Notes:
1185:   Users should typically use higher level functions such as `ISCreateGeneral()`.

1187:   This function can be useful in some special cases of constructing a new `IS`, e.g. after `ISCreate()` and before `ISLoad()`.
1188:   Otherwise, it is only valid to replace the layout with a layout known to be equivalent.

1190: .seealso: `IS`, `PetscLayout`, `ISCreate()`, `ISGetLayout()`, `ISGetSize()`, `ISGetLocalSize()`
1191: @*/
1192: PetscErrorCode ISSetLayout(IS is, PetscLayout map)
1193: {
1194:   PetscFunctionBegin;
1196:   PetscAssertPointer(map, 2);
1197:   PetscCall(PetscLayoutReference(map, &is->map));
1198:   PetscFunctionReturn(PETSC_SUCCESS);
1199: }

1201: /*@C
1202:   ISGetIndices - Returns a pointer to the indices.  The user should call
1203:   `ISRestoreIndices()` after having looked at the indices.  The user should
1204:   NOT change the indices.

1206:   Not Collective

1208:   Input Parameter:
1209: . is - the index set

1211:   Output Parameter:
1212: . ptr - the location to put the pointer to the indices

1214:   Level: intermediate

1216:   Fortran Notes:
1217:   `ISGetIndices()` Fortran binding is deprecated (since PETSc 3.19), use `ISGetIndicesF90()`

1219: .seealso: `IS`, `ISRestoreIndices()`, `ISGetIndicesF90()`
1220: @*/
1221: PetscErrorCode ISGetIndices(IS is, const PetscInt *ptr[])
1222: {
1223:   PetscFunctionBegin;
1225:   PetscAssertPointer(ptr, 2);
1226:   PetscUseTypeMethod(is, getindices, ptr);
1227:   PetscFunctionReturn(PETSC_SUCCESS);
1228: }

1230: /*@C
1231:   ISGetMinMax - Gets the minimum and maximum values in an `IS`

1233:   Not Collective

1235:   Input Parameter:
1236: . is - the index set

1238:   Output Parameters:
1239: + min - the minimum value
1240: - max - the maximum value

1242:   Level: intermediate

1244:   Notes:
1245:   Empty index sets return min=`PETSC_MAX_INT` and max=`PETSC_MIN_INT`.

1247:   In parallel, it returns the `min` and `max` of the local portion of `is`

1249: .seealso: `IS`, `ISGetIndices()`, `ISRestoreIndices()`, `ISGetIndicesF90()`
1250: @*/
1251: PetscErrorCode ISGetMinMax(IS is, PetscInt *min, PetscInt *max)
1252: {
1253:   PetscFunctionBegin;
1255:   if (min) *min = is->min;
1256:   if (max) *max = is->max;
1257:   PetscFunctionReturn(PETSC_SUCCESS);
1258: }

1260: /*@
1261:   ISLocate - determine the location of an index within the local component of an index set

1263:   Not Collective

1265:   Input Parameters:
1266: + is  - the index set
1267: - key - the search key

1269:   Output Parameter:
1270: . location - if >= 0, a location within the index set that is equal to the key, otherwise the key is not in the index set

1272:   Level: intermediate

1274: .seealso: `IS`
1275:  @*/
1276: PetscErrorCode ISLocate(IS is, PetscInt key, PetscInt *location)
1277: {
1278:   PetscFunctionBegin;
1279:   if (is->ops->locate) {
1280:     PetscUseTypeMethod(is, locate, key, location);
1281:   } else {
1282:     PetscInt        numIdx;
1283:     PetscBool       sorted;
1284:     const PetscInt *idx;

1286:     PetscCall(ISGetLocalSize(is, &numIdx));
1287:     PetscCall(ISGetIndices(is, &idx));
1288:     PetscCall(ISSorted(is, &sorted));
1289:     if (sorted) {
1290:       PetscCall(PetscFindInt(key, numIdx, idx, location));
1291:     } else {
1292:       PetscInt i;

1294:       *location = -1;
1295:       for (i = 0; i < numIdx; i++) {
1296:         if (idx[i] == key) {
1297:           *location = i;
1298:           break;
1299:         }
1300:       }
1301:     }
1302:     PetscCall(ISRestoreIndices(is, &idx));
1303:   }
1304:   PetscFunctionReturn(PETSC_SUCCESS);
1305: }

1307: /*@C
1308:   ISRestoreIndices - Restores an index set to a usable state after a call to `ISGetIndices()`.

1310:   Not Collective

1312:   Input Parameters:
1313: + is  - the index set
1314: - ptr - the pointer obtained by `ISGetIndices()`

1316:   Level: intermediate

1318:   Fortran Notes:
1319:   `ISRestoreIndices()` Fortran binding is deprecated (since PETSc 3.19), use `ISRestoreIndicesF90()`

1321: .seealso: `IS`, `ISGetIndices()`, `ISRestoreIndicesF90()`
1322: @*/
1323: PetscErrorCode ISRestoreIndices(IS is, const PetscInt *ptr[])
1324: {
1325:   PetscFunctionBegin;
1327:   PetscAssertPointer(ptr, 2);
1328:   PetscTryTypeMethod(is, restoreindices, ptr);
1329:   PetscFunctionReturn(PETSC_SUCCESS);
1330: }

1332: static PetscErrorCode ISGatherTotal_Private(IS is)
1333: {
1334:   PetscInt        i, n, N;
1335:   const PetscInt *lindices;
1336:   MPI_Comm        comm;
1337:   PetscMPIInt     rank, size, *sizes = NULL, *offsets = NULL, nn;

1339:   PetscFunctionBegin;

1342:   PetscCall(PetscObjectGetComm((PetscObject)is, &comm));
1343:   PetscCallMPI(MPI_Comm_size(comm, &size));
1344:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
1345:   PetscCall(ISGetLocalSize(is, &n));
1346:   PetscCall(PetscMalloc2(size, &sizes, size, &offsets));

1348:   PetscCall(PetscMPIIntCast(n, &nn));
1349:   PetscCallMPI(MPI_Allgather(&nn, 1, MPI_INT, sizes, 1, MPI_INT, comm));
1350:   offsets[0] = 0;
1351:   for (i = 1; i < size; ++i) offsets[i] = offsets[i - 1] + sizes[i - 1];
1352:   N = offsets[size - 1] + sizes[size - 1];

1354:   PetscCall(PetscMalloc1(N, &(is->total)));
1355:   PetscCall(ISGetIndices(is, &lindices));
1356:   PetscCallMPI(MPI_Allgatherv((void *)lindices, nn, MPIU_INT, is->total, sizes, offsets, MPIU_INT, comm));
1357:   PetscCall(ISRestoreIndices(is, &lindices));
1358:   is->local_offset = offsets[rank];
1359:   PetscCall(PetscFree2(sizes, offsets));
1360:   PetscFunctionReturn(PETSC_SUCCESS);
1361: }

1363: /*@C
1364:   ISGetTotalIndices - Retrieve an array containing all indices across the communicator.

1366:   Collective

1368:   Input Parameter:
1369: . is - the index set

1371:   Output Parameter:
1372: . indices - total indices with rank 0 indices first, and so on; total array size is
1373:              the same as returned with `ISGetSize()`.

1375:   Level: intermediate

1377:   Notes:
1378:   this is potentially nonscalable, but depends on the size of the total index set
1379:   and the size of the communicator. This may be feasible for index sets defined on
1380:   subcommunicators, such that the set size does not grow with `PETSC_WORLD_COMM`.
1381:   Note also that there is no way to tell where the local part of the indices starts
1382:   (use `ISGetIndices()` and `ISGetNonlocalIndices()` to retrieve just the local and just
1383:   the nonlocal part (complement), respectively).

1385: .seealso: `IS`, `ISRestoreTotalIndices()`, `ISGetNonlocalIndices()`, `ISGetSize()`
1386: @*/
1387: PetscErrorCode ISGetTotalIndices(IS is, const PetscInt *indices[])
1388: {
1389:   PetscMPIInt size;

1391:   PetscFunctionBegin;
1393:   PetscAssertPointer(indices, 2);
1394:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is), &size));
1395:   if (size == 1) {
1396:     PetscUseTypeMethod(is, getindices, indices);
1397:   } else {
1398:     if (!is->total) PetscCall(ISGatherTotal_Private(is));
1399:     *indices = is->total;
1400:   }
1401:   PetscFunctionReturn(PETSC_SUCCESS);
1402: }

1404: /*@C
1405:   ISRestoreTotalIndices - Restore the index array obtained with `ISGetTotalIndices()`.

1407:   Not Collective.

1409:   Input Parameters:
1410: + is      - the index set
1411: - indices - index array; must be the array obtained with `ISGetTotalIndices()`

1413:   Level: intermediate

1415: .seealso: `IS`, `ISGetNonlocalIndices()`
1416: @*/
1417: PetscErrorCode ISRestoreTotalIndices(IS is, const PetscInt *indices[])
1418: {
1419:   PetscMPIInt size;

1421:   PetscFunctionBegin;
1423:   PetscAssertPointer(indices, 2);
1424:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is), &size));
1425:   if (size == 1) {
1426:     PetscUseTypeMethod(is, restoreindices, indices);
1427:   } else {
1428:     PetscCheck(is->total == *indices, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Index array pointer being restored does not point to the array obtained from the IS.");
1429:   }
1430:   PetscFunctionReturn(PETSC_SUCCESS);
1431: }

1433: /*@C
1434:   ISGetNonlocalIndices - Retrieve an array of indices from remote processors
1435:   in this communicator.

1437:   Collective

1439:   Input Parameter:
1440: . is - the index set

1442:   Output Parameter:
1443: . indices - indices with rank 0 indices first, and so on,  omitting
1444:              the current rank.  Total number of indices is the difference
1445:              total and local, obtained with `ISGetSize()` and `ISGetLocalSize()`,
1446:              respectively.

1448:   Level: intermediate

1450:   Notes:
1451:   Restore the indices using `ISRestoreNonlocalIndices()`.

1453:   The same scalability considerations as those for `ISGetTotalIndices()` apply here.

1455: .seealso: `IS`, `ISGetTotalIndices()`, `ISRestoreNonlocalIndices()`, `ISGetSize()`, `ISGetLocalSize().`
1456: @*/
1457: PetscErrorCode ISGetNonlocalIndices(IS is, const PetscInt *indices[])
1458: {
1459:   PetscMPIInt size;
1460:   PetscInt    n, N;

1462:   PetscFunctionBegin;
1464:   PetscAssertPointer(indices, 2);
1465:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is), &size));
1466:   if (size == 1) *indices = NULL;
1467:   else {
1468:     if (!is->total) PetscCall(ISGatherTotal_Private(is));
1469:     PetscCall(ISGetLocalSize(is, &n));
1470:     PetscCall(ISGetSize(is, &N));
1471:     PetscCall(PetscMalloc1(N - n, &(is->nonlocal)));
1472:     PetscCall(PetscArraycpy(is->nonlocal, is->total, is->local_offset));
1473:     PetscCall(PetscArraycpy(is->nonlocal + is->local_offset, is->total + is->local_offset + n, N - is->local_offset - n));
1474:     *indices = is->nonlocal;
1475:   }
1476:   PetscFunctionReturn(PETSC_SUCCESS);
1477: }

1479: /*@C
1480:   ISRestoreNonlocalIndices - Restore the index array obtained with `ISGetNonlocalIndices()`.

1482:   Not Collective.

1484:   Input Parameters:
1485: + is      - the index set
1486: - indices - index array; must be the array obtained with `ISGetNonlocalIndices()`

1488:   Level: intermediate

1490: .seealso: `IS`, `ISGetTotalIndices()`, `ISGetNonlocalIndices()`, `ISRestoreTotalIndices()`
1491: @*/
1492: PetscErrorCode ISRestoreNonlocalIndices(IS is, const PetscInt *indices[])
1493: {
1494:   PetscFunctionBegin;
1496:   PetscAssertPointer(indices, 2);
1497:   PetscCheck(is->nonlocal == *indices, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Index array pointer being restored does not point to the array obtained from the IS.");
1498:   PetscFunctionReturn(PETSC_SUCCESS);
1499: }

1501: /*@
1502:   ISGetNonlocalIS - Gather all nonlocal indices for this `IS` and present
1503:   them as another sequential index set.

1505:   Collective

1507:   Input Parameter:
1508: . is - the index set

1510:   Output Parameter:
1511: . complement - sequential `IS` with indices identical to the result of
1512:                 `ISGetNonlocalIndices()`

1514:   Level: intermediate

1516:   Notes:
1517:   Complement represents the result of `ISGetNonlocalIndices()` as an `IS`.
1518:   Therefore scalability issues similar to `ISGetNonlocalIndices()` apply.

1520:   The resulting `IS` must be restored using `ISRestoreNonlocalIS()`.

1522: .seealso: `IS`, `ISGetNonlocalIndices()`, `ISRestoreNonlocalIndices()`, `ISAllGather()`, `ISGetSize()`
1523: @*/
1524: PetscErrorCode ISGetNonlocalIS(IS is, IS *complement)
1525: {
1526:   PetscFunctionBegin;
1528:   PetscAssertPointer(complement, 2);
1529:   /* Check if the complement exists already. */
1530:   if (is->complement) {
1531:     *complement = is->complement;
1532:     PetscCall(PetscObjectReference((PetscObject)(is->complement)));
1533:   } else {
1534:     PetscInt        N, n;
1535:     const PetscInt *idx;
1536:     PetscCall(ISGetSize(is, &N));
1537:     PetscCall(ISGetLocalSize(is, &n));
1538:     PetscCall(ISGetNonlocalIndices(is, &idx));
1539:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, N - n, idx, PETSC_USE_POINTER, &(is->complement)));
1540:     PetscCall(PetscObjectReference((PetscObject)is->complement));
1541:     *complement = is->complement;
1542:   }
1543:   PetscFunctionReturn(PETSC_SUCCESS);
1544: }

1546: /*@
1547:   ISRestoreNonlocalIS - Restore the `IS` obtained with `ISGetNonlocalIS()`.

1549:   Not collective.

1551:   Input Parameters:
1552: + is         - the index set
1553: - complement - index set of `is`'s nonlocal indices

1555:   Level: intermediate

1557: .seealso: `IS`, `ISGetNonlocalIS()`, `ISGetNonlocalIndices()`, `ISRestoreNonlocalIndices()`
1558: @*/
1559: PetscErrorCode ISRestoreNonlocalIS(IS is, IS *complement)
1560: {
1561:   PetscInt refcnt;

1563:   PetscFunctionBegin;
1565:   PetscAssertPointer(complement, 2);
1566:   PetscCheck(*complement == is->complement, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Complement IS being restored was not obtained with ISGetNonlocalIS()");
1567:   PetscCall(PetscObjectGetReference((PetscObject)(is->complement), &refcnt));
1568:   PetscCheck(refcnt > 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Duplicate call to ISRestoreNonlocalIS() detected");
1569:   PetscCall(PetscObjectDereference((PetscObject)(is->complement)));
1570:   PetscFunctionReturn(PETSC_SUCCESS);
1571: }

1573: /*@C
1574:   ISViewFromOptions - View an `IS` based on options in the options database

1576:   Collective

1578:   Input Parameters:
1579: + A    - the index set
1580: . obj  - Optional object that provides the prefix for the options database
1581: - name - command line option

1583:   Level: intermediate

1585:   Note:
1586:   See `PetscObjectViewFromOptions()` for possible `PetscViewer` and `PetscViewerFormat` values

1588: .seealso: `IS`, `ISView()`, `PetscObjectViewFromOptions()`, `ISCreate()`
1589: @*/
1590: PetscErrorCode ISViewFromOptions(IS A, PetscObject obj, const char name[])
1591: {
1592:   PetscFunctionBegin;
1594:   PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
1595:   PetscFunctionReturn(PETSC_SUCCESS);
1596: }

1598: /*@C
1599:   ISView - Displays an index set.

1601:   Collective

1603:   Input Parameters:
1604: + is     - the index set
1605: - viewer - viewer used to display the set, for example `PETSC_VIEWER_STDOUT_SELF`.

1607:   Level: intermediate

1609: .seealso: `IS`, `PetscViewer`, `PetscViewerASCIIOpen()`, `ISViewFromOptions()`
1610: @*/
1611: PetscErrorCode ISView(IS is, PetscViewer viewer)
1612: {
1613:   PetscFunctionBegin;
1615:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)is), &viewer));
1617:   PetscCheckSameComm(is, 1, viewer, 2);

1619:   PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)is, viewer));
1620:   PetscCall(PetscLogEventBegin(IS_View, is, viewer, 0, 0));
1621:   PetscUseTypeMethod(is, view, viewer);
1622:   PetscCall(PetscLogEventEnd(IS_View, is, viewer, 0, 0));
1623:   PetscFunctionReturn(PETSC_SUCCESS);
1624: }

1626: /*@
1627:   ISLoad - Loads a vector that has been stored in binary or HDF5 format with `ISView()`.

1629:   Collective

1631:   Input Parameters:
1632: + is     - the newly loaded index set, this needs to have been created with `ISCreate()` or some related function before a call to `ISLoad()`.
1633: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()` or HDF5 file viewer, obtained from `PetscViewerHDF5Open()`

1635:   Level: intermediate

1637:   Notes:
1638:   IF using HDF5, you must assign the IS the same name as was used in `is`
1639:   that was stored in the file using `PetscObjectSetName()`. Otherwise you will
1640:   get the error message: "Cannot H5DOpen2() with Vec name NAMEOFOBJECT"

1642: .seealso: `IS`, `PetscViewerBinaryOpen()`, `ISView()`, `MatLoad()`, `VecLoad()`
1643: @*/
1644: PetscErrorCode ISLoad(IS is, PetscViewer viewer)
1645: {
1646:   PetscBool isbinary, ishdf5;

1648:   PetscFunctionBegin;
1651:   PetscCheckSameComm(is, 1, viewer, 2);
1652:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1653:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1654:   PetscCheck(isbinary || ishdf5, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1655:   if (!((PetscObject)is)->type_name) PetscCall(ISSetType(is, ISGENERAL));
1656:   PetscCall(PetscLogEventBegin(IS_Load, is, viewer, 0, 0));
1657:   PetscUseTypeMethod(is, load, viewer);
1658:   PetscCall(PetscLogEventEnd(IS_Load, is, viewer, 0, 0));
1659:   PetscFunctionReturn(PETSC_SUCCESS);
1660: }

1662: /*@
1663:   ISSort - Sorts the indices of an index set.

1665:   Collective

1667:   Input Parameter:
1668: . is - the index set

1670:   Level: intermediate

1672: .seealso: `IS`, `ISSortRemoveDups()`, `ISSorted()`
1673: @*/
1674: PetscErrorCode ISSort(IS is)
1675: {
1676:   PetscFunctionBegin;
1678:   PetscUseTypeMethod(is, sort);
1679:   PetscCall(ISSetInfo(is, IS_SORTED, IS_LOCAL, is->info_permanent[IS_LOCAL][IS_SORTED], PETSC_TRUE));
1680:   PetscFunctionReturn(PETSC_SUCCESS);
1681: }

1683: /*@
1684:   ISSortRemoveDups - Sorts the indices of an index set, removing duplicates.

1686:   Collective

1688:   Input Parameter:
1689: . is - the index set

1691:   Level: intermediate

1693: .seealso: `IS`, `ISSort()`, `ISSorted()`
1694: @*/
1695: PetscErrorCode ISSortRemoveDups(IS is)
1696: {
1697:   PetscFunctionBegin;
1699:   PetscCall(ISClearInfoCache(is, PETSC_FALSE));
1700:   PetscUseTypeMethod(is, sortremovedups);
1701:   PetscCall(ISSetInfo(is, IS_SORTED, IS_LOCAL, is->info_permanent[IS_LOCAL][IS_SORTED], PETSC_TRUE));
1702:   PetscCall(ISSetInfo(is, IS_UNIQUE, IS_LOCAL, is->info_permanent[IS_LOCAL][IS_UNIQUE], PETSC_TRUE));
1703:   PetscFunctionReturn(PETSC_SUCCESS);
1704: }

1706: /*@
1707:   ISToGeneral - Converts an IS object of any type to `ISGENERAL` type

1709:   Collective

1711:   Input Parameter:
1712: . is - the index set

1714:   Level: intermediate

1716: .seealso: `IS`, `ISSorted()`
1717: @*/
1718: PetscErrorCode ISToGeneral(IS is)
1719: {
1720:   PetscFunctionBegin;
1722:   PetscUseTypeMethod(is, togeneral);
1723:   PetscFunctionReturn(PETSC_SUCCESS);
1724: }

1726: /*@
1727:   ISSorted - Checks the indices to determine whether they have been sorted.

1729:   Not Collective

1731:   Input Parameter:
1732: . is - the index set

1734:   Output Parameter:
1735: . flg - output flag, either `PETSC_TRUE` if the index set is sorted,
1736:          or `PETSC_FALSE` otherwise.

1738:   Level: intermediate

1740:   Note:
1741:   For parallel IS objects this only indicates if the local part of `is`
1742:   is sorted. So some processors may return `PETSC_TRUE` while others may
1743:   return `PETSC_FALSE`.

1745: .seealso: `ISSort()`, `ISSortRemoveDups()`
1746: @*/
1747: PetscErrorCode ISSorted(IS is, PetscBool *flg)
1748: {
1749:   PetscFunctionBegin;
1751:   PetscAssertPointer(flg, 2);
1752:   PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, flg));
1753:   PetscFunctionReturn(PETSC_SUCCESS);
1754: }

1756: /*@
1757:   ISDuplicate - Creates a duplicate copy of an index set.

1759:   Collective

1761:   Input Parameter:
1762: . is - the index set

1764:   Output Parameter:
1765: . newIS - the copy of the index set

1767:   Level: beginner

1769: .seealso: `IS`, `ISCreateGeneral()`, `ISCopy()`
1770: @*/
1771: PetscErrorCode ISDuplicate(IS is, IS *newIS)
1772: {
1773:   PetscFunctionBegin;
1775:   PetscAssertPointer(newIS, 2);
1776:   PetscUseTypeMethod(is, duplicate, newIS);
1777:   PetscCall(ISCopyInfo_Private(is, *newIS));
1778:   PetscFunctionReturn(PETSC_SUCCESS);
1779: }

1781: /*@
1782:   ISCopy - Copies an index set.

1784:   Collective

1786:   Input Parameter:
1787: . is - the index set

1789:   Output Parameter:
1790: . isy - the copy of the index set

1792:   Level: beginner

1794: .seealso: `IS`, `ISDuplicate()`, `ISShift()`
1795: @*/
1796: PetscErrorCode ISCopy(IS is, IS isy)
1797: {
1798:   PetscInt bs, bsy;

1800:   PetscFunctionBegin;
1803:   PetscCheckSameComm(is, 1, isy, 2);
1804:   if (is == isy) PetscFunctionReturn(PETSC_SUCCESS);
1805:   PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
1806:   PetscCall(PetscLayoutGetBlockSize(isy->map, &bsy));
1807:   PetscCheck(is->map->N == isy->map->N, PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_INCOMP, "Index sets have different global size %" PetscInt_FMT " != %" PetscInt_FMT, is->map->N, isy->map->N);
1808:   PetscCheck(is->map->n == isy->map->n, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Index sets have different local size %" PetscInt_FMT " != %" PetscInt_FMT, is->map->n, isy->map->n);
1809:   PetscCheck(bs == bsy, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Index sets have different block size %" PetscInt_FMT " != %" PetscInt_FMT, bs, bsy);
1810:   PetscCall(ISCopyInfo_Private(is, isy));
1811:   isy->max = is->max;
1812:   isy->min = is->min;
1813:   PetscUseTypeMethod(is, copy, isy);
1814:   PetscFunctionReturn(PETSC_SUCCESS);
1815: }

1817: /*@
1818:   ISShift - Shift all indices by given offset

1820:   Collective

1822:   Input Parameters:
1823: + is     - the index set
1824: - offset - the offset

1826:   Output Parameter:
1827: . isy - the shifted copy of the input index set

1829:   Level: beginner

1831:   Notes:
1832:   The `offset` can be different across processes.

1834:   `is` and `isy` can be the same.

1836: .seealso: `ISDuplicate()`, `ISCopy()`
1837: @*/
1838: PetscErrorCode ISShift(IS is, PetscInt offset, IS isy)
1839: {
1840:   PetscFunctionBegin;
1843:   PetscCheckSameComm(is, 1, isy, 3);
1844:   if (!offset) {
1845:     PetscCall(ISCopy(is, isy));
1846:     PetscFunctionReturn(PETSC_SUCCESS);
1847:   }
1848:   PetscCheck(is->map->N == isy->map->N, PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_INCOMP, "Index sets have different global size %" PetscInt_FMT " != %" PetscInt_FMT, is->map->N, isy->map->N);
1849:   PetscCheck(is->map->n == isy->map->n, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Index sets have different local size %" PetscInt_FMT " != %" PetscInt_FMT, is->map->n, isy->map->n);
1850:   PetscCheck(is->map->bs == isy->map->bs, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Index sets have different block size %" PetscInt_FMT " != %" PetscInt_FMT, is->map->bs, isy->map->bs);
1851:   PetscCall(ISCopyInfo_Private(is, isy));
1852:   isy->max = is->max + offset;
1853:   isy->min = is->min + offset;
1854:   PetscUseMethod(is, "ISShift_C", (IS, PetscInt, IS), (is, offset, isy));
1855:   PetscFunctionReturn(PETSC_SUCCESS);
1856: }

1858: /*@
1859:   ISOnComm - Split a parallel `IS` on subcomms (usually self) or concatenate index sets on subcomms into a parallel index set

1861:   Collective

1863:   Input Parameters:
1864: + is   - index set
1865: . comm - communicator for new index set
1866: - mode - copy semantics, `PETSC_USE_POINTER` for no-copy if possible, otherwise `PETSC_COPY_VALUES`

1868:   Output Parameter:
1869: . newis - new `IS` on `comm`

1871:   Level: advanced

1873:   Notes:
1874:   It is usually desirable to create a parallel `IS` and look at the local part when necessary.

1876:   This function is useful if serial `IS`s must be created independently, or to view many
1877:   logically independent serial `IS`s.

1879:   The input `IS` must have the same type on every MPI process.

1881: .seealso: `IS`
1882: @*/
1883: PetscErrorCode ISOnComm(IS is, MPI_Comm comm, PetscCopyMode mode, IS *newis)
1884: {
1885:   PetscMPIInt match;

1887:   PetscFunctionBegin;
1889:   PetscAssertPointer(newis, 4);
1890:   PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)is), comm, &match));
1891:   if (mode != PETSC_COPY_VALUES && (match == MPI_IDENT || match == MPI_CONGRUENT)) {
1892:     PetscCall(PetscObjectReference((PetscObject)is));
1893:     *newis = is;
1894:   } else PetscUseTypeMethod(is, oncomm, comm, mode, newis);
1895:   PetscFunctionReturn(PETSC_SUCCESS);
1896: }

1898: /*@
1899:   ISSetBlockSize - informs an index set that it has a given block size

1901:   Logicall Collective

1903:   Input Parameters:
1904: + is - index set
1905: - bs - block size

1907:   Level: intermediate

1909:   Notes:
1910:   This is much like the block size for `Vec`s. It indicates that one can think of the indices as
1911:   being in a collection of equal size blocks. For `ISBLOCK` these collections of blocks are all contiguous
1912:   within a block but this is not the case for other `IS`. For example, an `IS` with entries {0, 2, 3, 4, 6, 7} could
1913:   have a block size of three set.

1915:   `ISBlockGetIndices()` only works for `ISBLOCK`, not others.

1917: .seealso: `IS`, `ISGetBlockSize()`, `ISCreateBlock()`, `ISBlockGetIndices()`,
1918: @*/
1919: PetscErrorCode ISSetBlockSize(IS is, PetscInt bs)
1920: {
1921:   PetscFunctionBegin;
1924:   PetscCheck(bs >= 1, PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Block size %" PetscInt_FMT ", must be positive", bs);
1925:   if (PetscDefined(USE_DEBUG)) {
1926:     const PetscInt *indices;
1927:     PetscInt        length, i, j;
1928:     PetscCall(ISGetIndices(is, &indices));
1929:     if (indices) {
1930:       PetscCall(ISGetLocalSize(is, &length));
1931:       PetscCheck(length % bs == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local size %" PetscInt_FMT " not compatible with proposed block size %" PetscInt_FMT, length, bs);
1932:       for (i = 1; i < length / bs; i += bs) {
1933:         for (j = 1; j < bs - 1; j++) {
1934:           PetscCheck(indices[i * bs + j] == indices[(i - 1) * bs + j] + indices[i * bs] - indices[(i - 1) * bs], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Proposed block size %" PetscInt_FMT " is incompatible with the indices", bs);
1935:         }
1936:       }
1937:     }
1938:     PetscCall(ISRestoreIndices(is, &indices));
1939:   }
1940:   PetscUseTypeMethod(is, setblocksize, bs);
1941:   PetscFunctionReturn(PETSC_SUCCESS);
1942: }

1944: /*@
1945:   ISGetBlockSize - Returns the number of elements in a block.

1947:   Not Collective

1949:   Input Parameter:
1950: . is - the index set

1952:   Output Parameter:
1953: . size - the number of elements in a block

1955:   Level: intermediate

1957:   Note:
1958:   See `ISSetBlockSize()`

1960: .seealso: `IS`, `ISBlockGetSize()`, `ISGetSize()`, `ISCreateBlock()`, `ISSetBlockSize()`
1961: @*/
1962: PetscErrorCode ISGetBlockSize(IS is, PetscInt *size)
1963: {
1964:   PetscFunctionBegin;
1965:   PetscCall(PetscLayoutGetBlockSize(is->map, size));
1966:   PetscFunctionReturn(PETSC_SUCCESS);
1967: }

1969: static PetscErrorCode ISGetIndicesCopy_Private(IS is, PetscInt idx[])
1970: {
1971:   PetscInt        len, i;
1972:   const PetscInt *ptr;

1974:   PetscFunctionBegin;
1975:   PetscCall(ISGetLocalSize(is, &len));
1976:   PetscCall(ISGetIndices(is, &ptr));
1977:   for (i = 0; i < len; i++) idx[i] = ptr[i];
1978:   PetscCall(ISRestoreIndices(is, &ptr));
1979:   PetscFunctionReturn(PETSC_SUCCESS);
1980: }

1982: /*MC
1983:     ISGetIndicesF90 - Accesses the elements of an index set from Fortran.
1984:     The users should call `ISRestoreIndicesF90()` after having looked at the
1985:     indices.  The user should NOT change the indices.

1987:     Synopsis:
1988:     ISGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

1990:     Not Collective

1992:     Input Parameter:
1993: .   x - index set

1995:     Output Parameters:
1996: +   xx_v - the Fortran pointer to the array
1997: -   ierr - error code

1999:     Example of Usage:
2000: .vb
2001:     PetscInt, pointer xx_v(:)
2002:     ....
2003:     call ISGetIndicesF90(x,xx_v,ierr)
2004:     a = xx_v(3)
2005:     call ISRestoreIndicesF90(x,xx_v,ierr)
2006: .ve

2008:     Level: intermediate

2010: .seealso: `ISRestoreIndicesF90()`, `ISGetIndices()`, `ISRestoreIndices()`
2011: M*/

2013: /*MC
2014:     ISRestoreIndicesF90 - Restores an index set to a usable state after
2015:     a call to `ISGetIndicesF90()`.

2017:     Synopsis:
2018:     ISRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

2020:     Not Collective

2022:     Input Parameters:
2023: +   x - index set
2024: -   xx_v - the Fortran pointer to the array

2026:     Output Parameter:
2027: .   ierr - error code

2029:     Example of Usage:
2030: .vb
2031:     PetscInt, pointer xx_v(:)
2032:     ....
2033:     call ISGetIndicesF90(x,xx_v,ierr)
2034:     a = xx_v(3)
2035:     call ISRestoreIndicesF90(x,xx_v,ierr)
2036: .ve

2038:     Level: intermediate

2040: .seealso: `ISGetIndicesF90()`, `ISGetIndices()`, `ISRestoreIndices()`
2041: M*/

2043: /*MC
2044:     ISBlockGetIndicesF90 - Accesses the elements of an index set from Fortran.
2045:     The users should call `ISBlockRestoreIndicesF90()` after having looked at the
2046:     indices.  The user should NOT change the indices.

2048:     Synopsis:
2049:     ISBlockGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

2051:     Not Collective

2053:     Input Parameter:
2054: .   x - index set

2056:     Output Parameters:
2057: +   xx_v - the Fortran pointer to the array
2058: -   ierr - error code
2059:     Example of Usage:
2060: .vb
2061:     PetscInt, pointer xx_v(:)
2062:     ....
2063:     call ISBlockGetIndicesF90(x,xx_v,ierr)
2064:     a = xx_v(3)
2065:     call ISBlockRestoreIndicesF90(x,xx_v,ierr)
2066: .ve

2068:     Level: intermediate

2070: .seealso: `ISBlockRestoreIndicesF90()`, `ISGetIndices()`, `ISRestoreIndices()`,
2071:           `ISRestoreIndices()`
2072: M*/

2074: /*MC
2075:     ISBlockRestoreIndicesF90 - Restores an index set to a usable state after
2076:     a call to `ISBlockGetIndicesF90()`.

2078:     Synopsis:
2079:     ISBlockRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

2081:     Not Collective

2083:     Input Parameters:
2084: +   x - index set
2085: -   xx_v - the Fortran pointer to the array

2087:     Output Parameter:
2088: .   ierr - error code

2090:     Example of Usage:
2091: .vb
2092:     PetscInt, pointer xx_v(:)
2093:     ....
2094:     call ISBlockGetIndicesF90(x,xx_v,ierr)
2095:     a = xx_v(3)
2096:     call ISBlockRestoreIndicesF90(x,xx_v,ierr)
2097: .ve

2099:     Level: intermediate

2101: .seealso: `ISBlockGetIndicesF90()`, `ISGetIndices()`, `ISRestoreIndices()`, `ISRestoreIndicesF90()`
2102: M*/