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:   PetscCheck(gstart <= gend, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "gstart %" PetscInt_FMT " must be less than or equal to gend %" PetscInt_FMT, gstart, gend);
927:   *start  = -1;
928:   *contig = PETSC_FALSE;
929:   PetscTryTypeMethod(is, contiguous, gstart, gend, start, contig);
930:   PetscFunctionReturn(PETSC_SUCCESS);
931: }

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

937:   Logically Collective

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

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

945:   Level: intermediate

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

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

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

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

971:   Logically Collective

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

976:   Level: intermediate

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

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

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

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

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

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

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

1018:   Collective

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

1023:   Level: beginner

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

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

1055:   Collective

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

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

1065:   Level: intermediate

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

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

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

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

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

1105:   Not Collective

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

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

1113:   Level: beginner

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

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

1129:   Not Collective

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

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

1137:   Level: beginner

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

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

1153:   Not Collective

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

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

1161:   Level: developer

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

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

1177:   Collective

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

1183:   Level: developer

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

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

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

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

1207:   Not Collective

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

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

1215:   Level: intermediate

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

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

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

1234:   Not Collective

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

1239:   Output Parameters:
1240: + min - the minimum value, you may pass `NULL`
1241: - max - the maximum value, you may pass `NULL`

1243:   Level: intermediate

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

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

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

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

1264:   Not Collective

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

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

1273:   Level: intermediate

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

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

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

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

1311:   Not Collective

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

1317:   Level: intermediate

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

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

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

1340:   PetscFunctionBegin;

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

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

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

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

1367:   Collective

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

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

1376:   Level: intermediate

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

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

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

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

1408:   Not Collective.

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

1414:   Level: intermediate

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

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

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

1438:   Collective

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

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

1449:   Level: intermediate

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

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

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

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

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

1483:   Not Collective.

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

1489:   Level: intermediate

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

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

1506:   Collective

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

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

1515:   Level: intermediate

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

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

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

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

1550:   Not collective.

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

1556:   Level: intermediate

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

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

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

1577:   Collective

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

1584:   Level: intermediate

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

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

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

1602:   Collective

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

1608:   Level: intermediate

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

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

1627: /*@
1628:   ISLoad - Loads an index set that has been stored in binary or HDF5 format with `ISView()`.

1630:   Collective

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

1636:   Level: intermediate

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

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

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

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

1666:   Collective

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

1671:   Level: intermediate

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

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

1687:   Collective

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

1692:   Level: intermediate

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

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

1710:   Collective

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

1715:   Level: intermediate

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

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

1730:   Not Collective

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

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

1739:   Level: intermediate

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

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

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

1760:   Collective

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

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

1768:   Level: beginner

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

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

1785:   Collective

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

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

1793:   Level: beginner

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

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

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

1821:   Collective

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

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

1830:   Level: beginner

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

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

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

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

1862:   Collective

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

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

1872:   Level: advanced

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

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

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

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

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

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

1902:   Logicall Collective

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

1908:   Level: intermediate

1910:   Notes:
1911:   This is much like the block size for `Vec`s. It indicates that one can think of the indices as
1912:   being in a collection of equal size blocks. For `ISBLOCK` these collections of blocks are all contiguous
1913:   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
1914:   have a block size of three set.

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

1918: .seealso: `IS`, `ISGetBlockSize()`, `ISCreateBlock()`, `ISBlockGetIndices()`,
1919: @*/
1920: PetscErrorCode ISSetBlockSize(IS is, PetscInt bs)
1921: {
1922:   PetscFunctionBegin;
1925:   PetscCheck(bs >= 1, PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Block size %" PetscInt_FMT ", must be positive", bs);
1926:   if (PetscDefined(USE_DEBUG)) {
1927:     const PetscInt *indices;
1928:     PetscInt        length, i, j;
1929:     PetscCall(ISGetIndices(is, &indices));
1930:     if (indices) {
1931:       PetscCall(ISGetLocalSize(is, &length));
1932:       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);
1933:       for (i = 1; i < length / bs; i += bs) {
1934:         for (j = 1; j < bs - 1; j++) {
1935:           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);
1936:         }
1937:       }
1938:     }
1939:     PetscCall(ISRestoreIndices(is, &indices));
1940:   }
1941:   PetscUseTypeMethod(is, setblocksize, bs);
1942:   PetscFunctionReturn(PETSC_SUCCESS);
1943: }

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

1948:   Not Collective

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

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

1956:   Level: intermediate

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

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

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

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

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

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

1991:     Not Collective

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

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

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

2009:     Level: intermediate

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

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

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

2021:     Not Collective

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

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

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

2039:     Level: intermediate

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

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

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

2052:     Not Collective

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

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

2069:     Level: intermediate

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

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

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

2082:     Not Collective

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

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

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

2100:     Level: intermediate

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