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;
 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;

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()`, `ISClearInfoCache()`
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: /*@
418:    ISSetInfo - Set known information about an index set.

420:    Logically Collective on `is` if `ISInfoType` is `IS_GLOBAL`

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

431:   Info Describing IS Structure:
432: +    `IS_SORTED` - the [local part of the] index set is sorted in ascending order
433: .    `IS_UNIQUE` - each entry in the [local part of the] index set is unique
434: .    `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
435: .    `IS_INTERVAL` - the [local part of the] index set is equal to a contiguous range of integers {f, f + 1, ..., f + N-1}
436: -    `IS_IDENTITY` - the [local part of the] index set is equal to the integers {0, 1, ..., N-1}

438:    Level: advanced

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

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

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

452:   PetscFunctionBegin;
455:   comm = PetscObjectComm((PetscObject)is);
456:   if (type == IS_GLOBAL) {
460:     errcomm = comm;
461:   } else {
462:     errcomm = PETSC_COMM_SELF;
463:   }

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

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

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

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

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

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

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

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

525: PetscErrorCode ISGetIndicesCopy(IS is, PetscInt idx[]);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

751:    Input Parameters:
752: +  is - the index set
753: .  info - describing a property of the index set, one of those listed in the documentation of `ISSetInfo()`
754: .  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
755: -  type - whether the property is local (`IS_LOCAL`) or global (`IS_GLOBAL`)

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

760:    Level: advanced

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

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

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

777:   PetscFunctionBegin;
780:   comm = PetscObjectComm((PetscObject)is);
781:   if (type == IS_GLOBAL) {
783:     errcomm = comm;
784:   } else {
785:     errcomm = PETSC_COMM_SELF;
786:   }

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

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

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

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

844:    Collective

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

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

852:    Level: intermediate

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

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

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

875:    Logically Collective

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

880:    Level: intermediate

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

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

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

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

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

904:    Not Collective

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

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

915:    Level: developer

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

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

935:    Logically Collective

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

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

943:    Level: intermediate

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

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

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

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

969:    Logically Collective

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

974:    Level: intermediate

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

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

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

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

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

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

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

1016:    Collective

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

1021:    Level: beginner

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

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

1053:    Collective

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

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

1063:    Level: intermediate

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

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

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

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

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

1103:    Not Collective

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

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

1111:    Level: beginner

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

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

1127:    Not Collective

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

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

1135:    Level: beginner

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

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

1151:    Not Collective

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

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

1159:    Level: developer

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

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

1175:    Collective

1177:    Input Arguments:
1178: +  is - the index set
1179: -  map - the layout

1181:    Level: developer

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

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

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

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

1205:    Not Collective

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

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

1213:    Level: intermediate

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

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

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

1232:    Not Collective

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

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

1241:    Level: intermediate

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

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

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

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

1262:   Not Collective

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

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

1271:   Level: intermediate

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

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

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

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

1309:   Not Collective

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

1315:   Level: intermediate

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

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

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

1338:   PetscFunctionBegin;

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

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

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

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

1365:    Collective

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

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

1374:    Level: intermediate

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

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

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

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

1406:    Not Collective.

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

1412:    Level: intermediate

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

1420:   PetscFunctionBegin;
1423:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is), &size));
1424:   if (size == 1) {
1425:     PetscUseTypeMethod(is, restoreindices, indices);
1426:   } else {
1427:     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.");
1428:   }
1429:   PetscFunctionReturn(PETSC_SUCCESS);
1430: }

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

1436:    Collective

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

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

1447:    Level: intermediate

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

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

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

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

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

1481:    Not Collective.

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

1487:    Level: intermediate

1489: .seealso: `IS`, `ISGetTotalIndices()`, `ISGetNonlocalIndices()`, `ISRestoreTotalIndices()`
1490: @*/
1491: PetscErrorCode ISRestoreNonlocalIndices(IS is, const PetscInt *indices[])
1492: {
1493:   PetscFunctionBegin;
1496:   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.");
1497:   PetscFunctionReturn(PETSC_SUCCESS);
1498: }

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

1504:    Collective

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

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

1513:    Level: intermediate

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

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

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

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

1548:    Not collective.

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

1554:    Level: intermediate

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

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

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

1575:    Collective

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

1582:    Level: intermediate

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

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

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

1600:    Collective

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

1606:    Level: intermediate

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

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

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

1628:   Collective

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

1634:   Level: intermediate

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

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

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

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

1664:    Collective

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

1669:    Level: intermediate

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

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

1685:   Collective

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

1690:   Level: intermediate

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

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

1708:    Collective

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

1713:    Level: intermediate

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

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

1728:    Not Collective

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

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

1737:    Level: intermediate

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

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

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

1758:    Collective

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

1763:    Output Parameter:
1764: .  isnew - the copy of the index set

1766:    Level: beginner

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

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

1783:    Collective

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

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

1791:    Level: beginner

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

1799:   PetscFunctionBegin;
1802:   PetscCheckSameComm(is, 1, isy, 2);
1803:   if (is == isy) PetscFunctionReturn(PETSC_SUCCESS);
1804:   PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
1805:   PetscCall(PetscLayoutGetBlockSize(isy->map, &bsy));
1806:   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);
1807:   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);
1808:   PetscCheck(bs == bsy, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Index sets have different block size %" PetscInt_FMT " != %" PetscInt_FMT, bs, bsy);
1809:   PetscCall(ISCopyInfo(is, isy));
1810:   isy->max = is->max;
1811:   isy->min = is->min;
1812:   PetscUseTypeMethod(is, copy, isy);
1813:   PetscFunctionReturn(PETSC_SUCCESS);
1814: }

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

1819:    Collective

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

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

1828:    Level: beginner

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

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

1835: .seealso: `ISDuplicate()`, `ISCopy()`
1836: @*/
1837: PetscErrorCode ISShift(IS is, PetscInt offset, IS isy)
1838: {
1839:   PetscFunctionBegin;
1842:   PetscCheckSameComm(is, 1, isy, 3);
1843:   if (!offset) {
1844:     PetscCall(ISCopy(is, isy));
1845:     PetscFunctionReturn(PETSC_SUCCESS);
1846:   }
1847:   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);
1848:   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);
1849:   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);
1850:   PetscCall(ISCopyInfo(is, isy));
1851:   isy->max = is->max + offset;
1852:   isy->min = is->min + offset;
1853:   PetscUseMethod(is, "ISShift_C", (IS, PetscInt, IS), (is, offset, isy));
1854:   PetscFunctionReturn(PETSC_SUCCESS);
1855: }

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

1860:    Collective

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

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

1870:    Level: advanced

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

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

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

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

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

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

1900:    Logicall Collective

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

1906:    Level: intermediate

1908:    Notes:
1909:    This is much like the block size for `Vec`s. It indicates that one can think of the indices as
1910:    being in a collection of equal size blocks. For `ISBLOCK` these collections of blocks are all contiquous
1911:    within a block but this is not the case for other `IS`.

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

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

1942: /*@
1943:    ISGetBlockSize - Returns the number of elements in a block.

1945:    Not Collective

1947:    Input Parameter:
1948: .  is - the index set

1950:    Output Parameter:
1951: .  size - the number of elements in a block

1953:    Level: intermediate

1955:    Note:
1956:    See `ISSetBlockSize()`

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

1967: PetscErrorCode ISGetIndicesCopy(IS is, PetscInt idx[])
1968: {
1969:   PetscInt        len, i;
1970:   const PetscInt *ptr;

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

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

1985:     Synopsis:
1986:     ISGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

1988:     Not Collective

1990:     Input Parameter:
1991: .   x - index set

1993:     Output Parameters:
1994: +   xx_v - the Fortran pointer to the array
1995: -   ierr - error code

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

2006:     Level: intermediate

2008: .seealso: `ISRestoreIndicesF90()`, `ISGetIndices()`, `ISRestoreIndices()`
2009: M*/

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

2015:     Synopsis:
2016:     ISRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

2018:     Not Collective

2020:     Input Parameters:
2021: +   x - index set
2022: -   xx_v - the Fortran pointer to the array

2024:     Output Parameter:
2025: .   ierr - error code

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

2036:     Level: intermediate

2038: .seealso: `ISGetIndicesF90()`, `ISGetIndices()`, `ISRestoreIndices()`
2039: M*/

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

2046:     Synopsis:
2047:     ISBlockGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

2049:     Not Collective

2051:     Input Parameter:
2052: .   x - index set

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

2066:     Level: intermediate

2068: .seealso: `ISBlockRestoreIndicesF90()`, `ISGetIndices()`, `ISRestoreIndices()`,
2069:           `ISRestoreIndices()`
2070: M*/

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

2076:     Synopsis:
2077:     ISBlockRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

2079:     Not Collective

2081:     Input Parameters:
2082: +   x - index set
2083: -   xx_v - the Fortran pointer to the array

2085:     Output Parameter:
2086: .   ierr - error code

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

2097:     Level: intermediate

2099: .seealso: `ISBlockGetIndicesF90()`, `ISGetIndices()`, `ISRestoreIndices()`, `ISRestoreIndicesF90()`
2100: M*/