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_INT_MAX;
 63:   lbounds[1] = PETSC_INT_MIN;
 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 (PetscDefined(USE_DEBUG) && flg) {
330:       PetscInt        n;
331:       const PetscInt *indices;

333:       PetscCall(ISGetLocalSize(is, &n));
334:       PetscCall(ISGetIndices(is, &indices));
335:       PetscCall(PetscSortedInt(n, indices, &flg));
336:       if (type == IS_GLOBAL) PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &flg, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
337:       PetscCheck(flg, type == IS_GLOBAL ? PetscObjectComm((PetscObject)is) : PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "IS is not sorted");
338:       PetscCall(ISRestoreIndices(is, &indices));
339:     }
340:     if (flg && type == IS_GLOBAL) { /* an array that is globally sorted is also locally sorted */
341:       is->info[IS_LOCAL][(int)info] = IS_INFO_TRUE;
342:       /* global permanence implies local permanence */
343:       if (permanent_set && permanent) is->info_permanent[IS_LOCAL][(int)info] = PETSC_TRUE;
344:     }
345:     if (!flg) { /* if an array is not sorted, it cannot be an interval or the identity */
346:       is->info[itype][IS_INTERVAL] = IS_INFO_FALSE;
347:       is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
348:       if (permanent_set) {
349:         is->info_permanent[itype][IS_INTERVAL] = permanent;
350:         is->info_permanent[itype][IS_IDENTITY] = permanent;
351:       }
352:     }
353:     break;
354:   case IS_UNIQUE:
355:     if (flg && type == IS_GLOBAL) { /* an array that is globally unique is also locally unique */
356:       is->info[IS_LOCAL][(int)info] = IS_INFO_TRUE;
357:       /* global permanence implies local permanence */
358:       if (permanent_set && permanent) is->info_permanent[IS_LOCAL][(int)info] = PETSC_TRUE;
359:     }
360:     if (!flg) { /* if an array is not unique, it cannot be a permutation, and interval, or the identity */
361:       is->info[itype][IS_PERMUTATION] = IS_INFO_FALSE;
362:       is->info[itype][IS_INTERVAL]    = IS_INFO_FALSE;
363:       is->info[itype][IS_IDENTITY]    = IS_INFO_FALSE;
364:       if (permanent_set) {
365:         is->info_permanent[itype][IS_PERMUTATION] = permanent;
366:         is->info_permanent[itype][IS_INTERVAL]    = permanent;
367:         is->info_permanent[itype][IS_IDENTITY]    = permanent;
368:       }
369:     }
370:     break;
371:   case IS_PERMUTATION:
372:     if (flg) { /* an array that is a permutation is unique and is unique locally */
373:       is->info[itype][IS_UNIQUE]    = IS_INFO_TRUE;
374:       is->info[IS_LOCAL][IS_UNIQUE] = IS_INFO_TRUE;
375:       if (permanent_set && permanent) {
376:         is->info_permanent[itype][IS_UNIQUE]    = PETSC_TRUE;
377:         is->info_permanent[IS_LOCAL][IS_UNIQUE] = PETSC_TRUE;
378:       }
379:     } else { /* an array that is not a permutation cannot be the identity */
380:       is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
381:       if (permanent_set) is->info_permanent[itype][IS_IDENTITY] = permanent;
382:     }
383:     break;
384:   case IS_INTERVAL:
385:     if (flg) { /* an array that is an interval is sorted and unique */
386:       is->info[itype][IS_SORTED]    = IS_INFO_TRUE;
387:       is->info[IS_LOCAL][IS_SORTED] = IS_INFO_TRUE;
388:       is->info[itype][IS_UNIQUE]    = IS_INFO_TRUE;
389:       is->info[IS_LOCAL][IS_UNIQUE] = IS_INFO_TRUE;
390:       if (permanent_set && permanent) {
391:         is->info_permanent[itype][IS_SORTED]    = PETSC_TRUE;
392:         is->info_permanent[IS_LOCAL][IS_SORTED] = PETSC_TRUE;
393:         is->info_permanent[itype][IS_UNIQUE]    = PETSC_TRUE;
394:         is->info_permanent[IS_LOCAL][IS_UNIQUE] = PETSC_TRUE;
395:       }
396:     } else { /* an array that is not an interval cannot be the identity */
397:       is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
398:       if (permanent_set) is->info_permanent[itype][IS_IDENTITY] = permanent;
399:     }
400:     break;
401:   case IS_IDENTITY:
402:     if (flg) { /* an array that is the identity is sorted, unique, an interval, and a permutation */
403:       is->info[itype][IS_SORTED]      = IS_INFO_TRUE;
404:       is->info[IS_LOCAL][IS_SORTED]   = IS_INFO_TRUE;
405:       is->info[itype][IS_UNIQUE]      = IS_INFO_TRUE;
406:       is->info[IS_LOCAL][IS_UNIQUE]   = IS_INFO_TRUE;
407:       is->info[itype][IS_PERMUTATION] = IS_INFO_TRUE;
408:       is->info[itype][IS_INTERVAL]    = IS_INFO_TRUE;
409:       is->info[IS_LOCAL][IS_INTERVAL] = IS_INFO_TRUE;
410:       if (permanent_set && permanent) {
411:         is->info_permanent[itype][IS_SORTED]      = PETSC_TRUE;
412:         is->info_permanent[IS_LOCAL][IS_SORTED]   = PETSC_TRUE;
413:         is->info_permanent[itype][IS_UNIQUE]      = PETSC_TRUE;
414:         is->info_permanent[IS_LOCAL][IS_UNIQUE]   = PETSC_TRUE;
415:         is->info_permanent[itype][IS_PERMUTATION] = PETSC_TRUE;
416:         is->info_permanent[itype][IS_INTERVAL]    = PETSC_TRUE;
417:         is->info_permanent[IS_LOCAL][IS_INTERVAL] = PETSC_TRUE;
418:       }
419:     }
420:     break;
421:   default:
422:     PetscCheck(type != IS_LOCAL, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown IS property");
423:     SETERRQ(PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Unknown IS property");
424:   }
425:   PetscFunctionReturn(PETSC_SUCCESS);
426: }

428: // PetscClangLinter pragma disable: -fdoc-section-header-unknown
429: /*@
430:   ISSetInfo - Set known information about an index set.

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

434:   Input Parameters:
435: + is        - the index set
436: . info      - describing a property of the index set, one of those listed below,
437: . type      - `IS_LOCAL` if the information describes the local portion of the index set,
438:           `IS_GLOBAL` if it describes the whole index set
439: . permanent - `PETSC_TRUE` if it is known that the property will persist through changes to the index set, `PETSC_FALSE` otherwise
440:                If the user sets a property as permanently known, it will bypass computation of that property
441: - flg       - set the described property as true (`PETSC_TRUE`) or false (`PETSC_FALSE`)

443:   Values of `info` Describing `IS` Structure:
444: + `IS_SORTED`      - the [local part of the] index set is sorted in ascending order
445: . `IS_UNIQUE`      - each entry in the [local part of the] index set is unique
446: . `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
447: . `IS_INTERVAL`    - the [local part of the] index set is equal to a contiguous range of integers {f, f + 1, ..., f + N-1}
448: - `IS_IDENTITY`    - the [local part of the] index set is equal to the integers {0, 1, ..., N-1}

450:   Level: advanced

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

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

457: .seealso: `ISInfo`, `ISInfoType`, `IS`
458: @*/
459: PetscErrorCode ISSetInfo(IS is, ISInfo info, ISInfoType type, PetscBool permanent, PetscBool flg)
460: {
461:   MPI_Comm    comm, errcomm;
462:   PetscMPIInt size;

464:   PetscFunctionBegin;
467:   comm = PetscObjectComm((PetscObject)is);
468:   if (type == IS_GLOBAL) {
472:     errcomm = comm;
473:   } else {
474:     errcomm = PETSC_COMM_SELF;
475:   }

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

479:   PetscCallMPI(MPI_Comm_size(comm, &size));
480:   /* do not use global values if size == 1: it makes it easier to keep the implications straight */
481:   if (size == 1) type = IS_LOCAL;
482:   PetscCall(ISSetInfo_Internal(is, info, type, permanent ? IS_INFO_TRUE : IS_INFO_FALSE, flg));
483:   PetscFunctionReturn(PETSC_SUCCESS);
484: }

486: static PetscErrorCode ISGetInfo_Sorted_Private(IS is, ISInfoType type, PetscBool *flg)
487: {
488:   MPI_Comm    comm;
489:   PetscMPIInt size, rank;

491:   PetscFunctionBegin;
492:   comm = PetscObjectComm((PetscObject)is);
493:   PetscCallMPI(MPI_Comm_size(comm, &size));
494:   PetscCallMPI(MPI_Comm_size(comm, &rank));
495:   if (type == IS_GLOBAL && is->ops->sortedglobal) {
496:     PetscUseTypeMethod(is, sortedglobal, flg);
497:   } else {
498:     PetscBool sortedLocal = PETSC_FALSE;

500:     /* determine if the array is locally sorted */
501:     if (type == IS_GLOBAL && size > 1) {
502:       /* call ISGetInfo so that a cached value will be used if possible */
503:       PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, &sortedLocal));
504:     } else if (is->ops->sortedlocal) {
505:       PetscUseTypeMethod(is, sortedlocal, &sortedLocal);
506:     } else {
507:       /* default: get the local indices and directly check */
508:       const PetscInt *idx;
509:       PetscInt        n;

511:       PetscCall(ISGetIndices(is, &idx));
512:       PetscCall(ISGetLocalSize(is, &n));
513:       PetscCall(PetscSortedInt(n, idx, &sortedLocal));
514:       PetscCall(ISRestoreIndices(is, &idx));
515:     }

517:     if (type == IS_LOCAL || size == 1) {
518:       *flg = sortedLocal;
519:     } else {
520:       PetscCallMPI(MPIU_Allreduce(&sortedLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm));
521:       if (*flg) {
522:         PetscInt n, min = PETSC_INT_MAX, max = PETSC_INT_MIN;
523:         PetscInt maxprev;

525:         PetscCall(ISGetLocalSize(is, &n));
526:         if (n) PetscCall(ISGetMinMax(is, &min, &max));
527:         maxprev = PETSC_INT_MIN;
528:         PetscCallMPI(MPI_Exscan(&max, &maxprev, 1, MPIU_INT, MPI_MAX, comm));
529:         if (rank && (maxprev > min)) sortedLocal = PETSC_FALSE;
530:         PetscCallMPI(MPIU_Allreduce(&sortedLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm));
531:       }
532:     }
533:   }
534:   PetscFunctionReturn(PETSC_SUCCESS);
535: }

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

539: static PetscErrorCode ISGetInfo_Unique_Private(IS is, ISInfoType type, PetscBool *flg)
540: {
541:   MPI_Comm    comm;
542:   PetscMPIInt size, rank;
543:   PetscInt    i;

545:   PetscFunctionBegin;
546:   comm = PetscObjectComm((PetscObject)is);
547:   PetscCallMPI(MPI_Comm_size(comm, &size));
548:   PetscCallMPI(MPI_Comm_size(comm, &rank));
549:   if (type == IS_GLOBAL && is->ops->uniqueglobal) {
550:     PetscUseTypeMethod(is, uniqueglobal, flg);
551:   } else {
552:     PetscBool uniqueLocal;
553:     PetscInt  n   = -1;
554:     PetscInt *idx = NULL;

556:     /* determine if the array is locally unique */
557:     if (type == IS_GLOBAL && size > 1) {
558:       /* call ISGetInfo so that a cached value will be used if possible */
559:       PetscCall(ISGetInfo(is, IS_UNIQUE, IS_LOCAL, PETSC_TRUE, &uniqueLocal));
560:     } else if (is->ops->uniquelocal) {
561:       PetscUseTypeMethod(is, uniquelocal, &uniqueLocal);
562:     } else {
563:       /* default: get the local indices and directly check */
564:       uniqueLocal = PETSC_TRUE;
565:       PetscCall(ISGetLocalSize(is, &n));
566:       PetscCall(PetscMalloc1(n, &idx));
567:       PetscCall(ISGetIndicesCopy_Private(is, idx));
568:       PetscCall(PetscIntSortSemiOrdered(n, idx));
569:       for (i = 1; i < n; i++)
570:         if (idx[i] == idx[i - 1]) break;
571:       if (i < n) uniqueLocal = PETSC_FALSE;
572:     }

574:     PetscCall(PetscFree(idx));
575:     if (type == IS_LOCAL || size == 1) {
576:       *flg = uniqueLocal;
577:     } else {
578:       PetscCallMPI(MPIU_Allreduce(&uniqueLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm));
579:       if (*flg) {
580:         PetscInt min = PETSC_INT_MAX, max = PETSC_INT_MIN, maxprev;

582:         if (!idx) {
583:           PetscCall(ISGetLocalSize(is, &n));
584:           PetscCall(PetscMalloc1(n, &idx));
585:           PetscCall(ISGetIndicesCopy_Private(is, idx));
586:         }
587:         PetscCall(PetscParallelSortInt(is->map, is->map, idx, idx));
588:         if (n) {
589:           min = idx[0];
590:           max = idx[n - 1];
591:         }
592:         for (i = 1; i < n; i++)
593:           if (idx[i] == idx[i - 1]) break;
594:         if (i < n) uniqueLocal = PETSC_FALSE;
595:         maxprev = PETSC_INT_MIN;
596:         PetscCallMPI(MPI_Exscan(&max, &maxprev, 1, MPIU_INT, MPI_MAX, comm));
597:         if (rank && (maxprev == min)) uniqueLocal = PETSC_FALSE;
598:         PetscCallMPI(MPIU_Allreduce(&uniqueLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm));
599:       }
600:     }
601:     PetscCall(PetscFree(idx));
602:   }
603:   PetscFunctionReturn(PETSC_SUCCESS);
604: }

606: static PetscErrorCode ISGetInfo_Permutation(IS is, ISInfoType type, PetscBool *flg)
607: {
608:   MPI_Comm    comm;
609:   PetscMPIInt size, rank;

611:   PetscFunctionBegin;
612:   comm = PetscObjectComm((PetscObject)is);
613:   PetscCallMPI(MPI_Comm_size(comm, &size));
614:   PetscCallMPI(MPI_Comm_size(comm, &rank));
615:   if (type == IS_GLOBAL && is->ops->permglobal) {
616:     PetscUseTypeMethod(is, permglobal, flg);
617:   } else if (type == IS_LOCAL && is->ops->permlocal) {
618:     PetscUseTypeMethod(is, permlocal, flg);
619:   } else {
620:     PetscBool permLocal;
621:     PetscInt  n, i, rStart;
622:     PetscInt *idx;

624:     PetscCall(ISGetLocalSize(is, &n));
625:     PetscCall(PetscMalloc1(n, &idx));
626:     PetscCall(ISGetIndicesCopy_Private(is, idx));
627:     if (type == IS_GLOBAL) {
628:       PetscCall(PetscParallelSortInt(is->map, is->map, idx, idx));
629:       PetscCall(PetscLayoutGetRange(is->map, &rStart, NULL));
630:     } else {
631:       PetscCall(PetscIntSortSemiOrdered(n, idx));
632:       rStart = 0;
633:     }
634:     permLocal = PETSC_TRUE;
635:     for (i = 0; i < n; i++) {
636:       if (idx[i] != rStart + i) break;
637:     }
638:     if (i < n) permLocal = PETSC_FALSE;
639:     if (type == IS_LOCAL || size == 1) {
640:       *flg = permLocal;
641:     } else {
642:       PetscCallMPI(MPIU_Allreduce(&permLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm));
643:     }
644:     PetscCall(PetscFree(idx));
645:   }
646:   PetscFunctionReturn(PETSC_SUCCESS);
647: }

649: static PetscErrorCode ISGetInfo_Interval(IS is, ISInfoType type, PetscBool *flg)
650: {
651:   MPI_Comm    comm;
652:   PetscMPIInt size, rank;
653:   PetscInt    i;

655:   PetscFunctionBegin;
656:   comm = PetscObjectComm((PetscObject)is);
657:   PetscCallMPI(MPI_Comm_size(comm, &size));
658:   PetscCallMPI(MPI_Comm_size(comm, &rank));
659:   if (type == IS_GLOBAL && is->ops->intervalglobal) {
660:     PetscUseTypeMethod(is, intervalglobal, flg);
661:   } else {
662:     PetscBool intervalLocal;

664:     /* determine if the array is locally an interval */
665:     if (type == IS_GLOBAL && size > 1) {
666:       /* call ISGetInfo so that a cached value will be used if possible */
667:       PetscCall(ISGetInfo(is, IS_INTERVAL, IS_LOCAL, PETSC_TRUE, &intervalLocal));
668:     } else if (is->ops->intervallocal) {
669:       PetscUseTypeMethod(is, intervallocal, &intervalLocal);
670:     } else {
671:       PetscInt        n;
672:       const PetscInt *idx;
673:       /* default: get the local indices and directly check */
674:       intervalLocal = PETSC_TRUE;
675:       PetscCall(ISGetLocalSize(is, &n));
676:       PetscCall(ISGetIndices(is, &idx));
677:       for (i = 1; i < n; i++)
678:         if (idx[i] != idx[i - 1] + 1) break;
679:       if (i < n) intervalLocal = PETSC_FALSE;
680:       PetscCall(ISRestoreIndices(is, &idx));
681:     }

683:     if (type == IS_LOCAL || size == 1) {
684:       *flg = intervalLocal;
685:     } else {
686:       PetscCallMPI(MPIU_Allreduce(&intervalLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm));
687:       if (*flg) {
688:         PetscInt n, min = PETSC_INT_MAX, max = PETSC_INT_MIN;
689:         PetscInt maxprev;

691:         PetscCall(ISGetLocalSize(is, &n));
692:         if (n) PetscCall(ISGetMinMax(is, &min, &max));
693:         maxprev = PETSC_INT_MIN;
694:         PetscCallMPI(MPI_Exscan(&max, &maxprev, 1, MPIU_INT, MPI_MAX, comm));
695:         if (rank && n && (maxprev != min - 1)) intervalLocal = PETSC_FALSE;
696:         PetscCallMPI(MPIU_Allreduce(&intervalLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm));
697:       }
698:     }
699:   }
700:   PetscFunctionReturn(PETSC_SUCCESS);
701: }

703: static PetscErrorCode ISGetInfo_Identity(IS is, ISInfoType type, PetscBool *flg)
704: {
705:   MPI_Comm    comm;
706:   PetscMPIInt size, rank;

708:   PetscFunctionBegin;
709:   comm = PetscObjectComm((PetscObject)is);
710:   PetscCallMPI(MPI_Comm_size(comm, &size));
711:   PetscCallMPI(MPI_Comm_size(comm, &rank));
712:   if (type == IS_GLOBAL && is->ops->intervalglobal) {
713:     PetscBool isinterval;

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

720:       PetscCall(ISGetMinMax(is, &min, NULL));
721:       PetscCallMPI(MPI_Bcast(&min, 1, MPIU_INT, 0, comm));
722:       if (min == 0) *flg = PETSC_TRUE;
723:     }
724:   } else if (type == IS_LOCAL && is->ops->intervallocal) {
725:     PetscBool isinterval;

727:     PetscUseTypeMethod(is, intervallocal, &isinterval);
728:     *flg = PETSC_FALSE;
729:     if (isinterval) {
730:       PetscInt min;

732:       PetscCall(ISGetMinMax(is, &min, NULL));
733:       if (min == 0) *flg = PETSC_TRUE;
734:     }
735:   } else {
736:     PetscBool       identLocal;
737:     PetscInt        n, i, rStart;
738:     const PetscInt *idx;

740:     PetscCall(ISGetLocalSize(is, &n));
741:     PetscCall(ISGetIndices(is, &idx));
742:     PetscCall(PetscLayoutGetRange(is->map, &rStart, NULL));
743:     identLocal = PETSC_TRUE;
744:     for (i = 0; i < n; i++) {
745:       if (idx[i] != rStart + i) break;
746:     }
747:     if (i < n) identLocal = PETSC_FALSE;
748:     if (type == IS_LOCAL || size == 1) {
749:       *flg = identLocal;
750:     } else {
751:       PetscCallMPI(MPIU_Allreduce(&identLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm));
752:     }
753:     PetscCall(ISRestoreIndices(is, &idx));
754:   }
755:   PetscFunctionReturn(PETSC_SUCCESS);
756: }

758: /*@
759:   ISGetInfo - Determine whether an index set satisfies a given property

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

763:   Input Parameters:
764: + is      - the index set
765: . info    - describing a property of the index set, one of those listed in the documentation of `ISSetInfo()`
766: . 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
767: - type    - whether the property is local (`IS_LOCAL`) or global (`IS_GLOBAL`)

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

772:   Level: advanced

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

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

779: .seealso: `IS`, `ISInfo`, `ISInfoType`, `ISSetInfo()`, `ISClearInfoCache()`
780: @*/
781: PetscErrorCode ISGetInfo(IS is, ISInfo info, ISInfoType type, PetscBool compute, PetscBool *flg)
782: {
783:   MPI_Comm    comm, errcomm;
784:   PetscMPIInt rank, size;
785:   PetscInt    itype;
786:   PetscBool   hasprop;
787:   PetscBool   infer;

789:   PetscFunctionBegin;
792:   comm = PetscObjectComm((PetscObject)is);
793:   if (type == IS_GLOBAL) {
795:     errcomm = comm;
796:   } else {
797:     errcomm = PETSC_COMM_SELF;
798:   }

800:   PetscCallMPI(MPI_Comm_size(comm, &size));
801:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

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

845: static PetscErrorCode ISCopyInfo_Private(IS source, IS dest)
846: {
847:   PetscFunctionBegin;
848:   PetscCall(PetscArraycpy(&dest->info[0], &source->info[0], 2));
849:   PetscCall(PetscArraycpy(&dest->info_permanent[0], &source->info_permanent[0], 2));
850:   PetscFunctionReturn(PETSC_SUCCESS);
851: }

853: /*@
854:   ISIdentity - Determines whether index set is the identity mapping.

856:   Collective

858:   Input Parameter:
859: . is - the index set

861:   Output Parameter:
862: . ident - `PETSC_TRUE` if an identity, else `PETSC_FALSE`

864:   Level: intermediate

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

873: .seealso: `IS`, `ISSetIdentity()`, `ISGetInfo()`
874: @*/
875: PetscErrorCode ISIdentity(IS is, PetscBool *ident)
876: {
877:   PetscFunctionBegin;
879:   PetscAssertPointer(ident, 2);
880:   PetscCall(ISGetInfo(is, IS_IDENTITY, IS_GLOBAL, PETSC_TRUE, ident));
881:   PetscFunctionReturn(PETSC_SUCCESS);
882: }

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

887:   Logically Collective

889:   Input Parameter:
890: . is - the index set

892:   Level: intermediate

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

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

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

903: .seealso: `IS`, `ISIdentity()`, `ISSetInfo()`, `ISClearInfoCache()`
904: @*/
905: PetscErrorCode ISSetIdentity(IS is)
906: {
907:   PetscFunctionBegin;
909:   PetscCall(ISSetInfo(is, IS_IDENTITY, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE));
910:   PetscFunctionReturn(PETSC_SUCCESS);
911: }

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

916:   Not Collective

918:   Input Parameters:
919: + is     - the index set
920: . gstart - global start
921: - gend   - global end

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

927:   Level: developer

929: .seealso: `IS`, `ISGetLocalSize()`, `VecGetOwnershipRange()`
930: @*/
931: PetscErrorCode ISContiguousLocal(IS is, PetscInt gstart, PetscInt gend, PetscInt *start, PetscBool *contig)
932: {
933:   PetscFunctionBegin;
935:   PetscAssertPointer(start, 4);
936:   PetscAssertPointer(contig, 5);
937:   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);
938:   *start  = -1;
939:   *contig = PETSC_FALSE;
940:   PetscTryTypeMethod(is, contiguous, gstart, gend, start, contig);
941:   PetscFunctionReturn(PETSC_SUCCESS);
942: }

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

948:   Logically Collective

950:   Input Parameter:
951: . is - the index set

953:   Output Parameter:
954: . perm - `PETSC_TRUE` if a permutation, else `PETSC_FALSE`

956:   Level: intermediate

958:   Note:
959:   If it is not already known that `is` is a permutation (if `ISSetPermutation()`
960:   or `ISSetInfo()` has not been called), this routine will not attempt to compute
961:   whether the index set is a permutation and will assume `perm` is `PETSC_FALSE`.
962:   To compute the value when it is not already known, use `ISGetInfo()` with
963:   the compute flag set to `PETSC_TRUE`.

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

968: .seealso: `IS`, `ISSetPermutation()`, `ISGetInfo()`
969: @*/
970: PetscErrorCode ISPermutation(IS is, PetscBool *perm)
971: {
972:   PetscFunctionBegin;
974:   PetscAssertPointer(perm, 2);
975:   PetscCall(ISGetInfo(is, IS_PERMUTATION, IS_GLOBAL, PETSC_FALSE, perm));
976:   PetscFunctionReturn(PETSC_SUCCESS);
977: }

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

982:   Logically Collective

984:   Input Parameter:
985: . is - the index set

987:   Level: intermediate

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

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

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

998: .seealso: `IS`, `ISPermutation()`, `ISSetInfo()`, `ISClearInfoCache().`
999: @*/
1000: PetscErrorCode ISSetPermutation(IS is)
1001: {
1002:   PetscFunctionBegin;
1004:   if (PetscDefined(USE_DEBUG)) {
1005:     PetscMPIInt size;

1007:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is), &size));
1008:     if (size == 1) {
1009:       PetscInt        i, n, *idx;
1010:       const PetscInt *iidx;

1012:       PetscCall(ISGetSize(is, &n));
1013:       PetscCall(PetscMalloc1(n, &idx));
1014:       PetscCall(ISGetIndices(is, &iidx));
1015:       PetscCall(PetscArraycpy(idx, iidx, n));
1016:       PetscCall(PetscIntSortSemiOrdered(n, idx));
1017:       for (i = 0; i < n; i++) PetscCheck(idx[i] == i, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Index set is not a permutation");
1018:       PetscCall(PetscFree(idx));
1019:       PetscCall(ISRestoreIndices(is, &iidx));
1020:     }
1021:   }
1022:   PetscCall(ISSetInfo(is, IS_PERMUTATION, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE));
1023:   PetscFunctionReturn(PETSC_SUCCESS);
1024: }

1026: /*@
1027:   ISDestroy - Destroys an index set.

1029:   Collective

1031:   Input Parameter:
1032: . is - the index set

1034:   Level: beginner

1036: .seealso: `IS`, `ISCreateGeneral()`, `ISCreateStride()`, `ISCreateBlock()`
1037: @*/
1038: PetscErrorCode ISDestroy(IS *is)
1039: {
1040:   PetscFunctionBegin;
1041:   if (!*is) PetscFunctionReturn(PETSC_SUCCESS);
1043:   if (--((PetscObject)*is)->refct > 0) {
1044:     *is = NULL;
1045:     PetscFunctionReturn(PETSC_SUCCESS);
1046:   }
1047:   if ((*is)->complement) {
1048:     PetscInt refcnt;
1049:     PetscCall(PetscObjectGetReference((PetscObject)((*is)->complement), &refcnt));
1050:     PetscCheck(refcnt <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Nonlocal IS has not been restored");
1051:     PetscCall(ISDestroy(&(*is)->complement));
1052:   }
1053:   PetscTryTypeMethod(*is, destroy);
1054:   PetscCall(PetscLayoutDestroy(&(*is)->map));
1055:   /* Destroy local representations of offproc data. */
1056:   PetscCall(PetscFree((*is)->total));
1057:   PetscCall(PetscFree((*is)->nonlocal));
1058:   PetscCall(PetscHeaderDestroy(is));
1059:   PetscFunctionReturn(PETSC_SUCCESS);
1060: }

1062: /*@
1063:   ISInvertPermutation - Creates a new permutation that is the inverse of
1064:   a given permutation.

1066:   Collective

1068:   Input Parameters:
1069: + is     - the index set
1070: - nlocal - number of indices on this processor in result (ignored for 1 processor) or
1071:             use `PETSC_DECIDE`

1073:   Output Parameter:
1074: . isout - the inverse permutation

1076:   Level: intermediate

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

1082: .seealso: `IS`, `ISGetInfo()`, `ISSetPermutation()`, `ISGetPermutation()`
1083: @*/
1084: PetscErrorCode ISInvertPermutation(IS is, PetscInt nlocal, IS *isout)
1085: {
1086:   PetscBool isperm, isidentity, issame;

1088:   PetscFunctionBegin;
1090:   PetscAssertPointer(isout, 3);
1091:   PetscCall(ISGetInfo(is, IS_PERMUTATION, IS_GLOBAL, PETSC_TRUE, &isperm));
1092:   PetscCheck(isperm, PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_WRONG, "Not a permutation");
1093:   PetscCall(ISGetInfo(is, IS_IDENTITY, IS_GLOBAL, PETSC_TRUE, &isidentity));
1094:   issame = PETSC_FALSE;
1095:   if (isidentity) {
1096:     PetscInt  n;
1097:     PetscBool isallsame;

1099:     PetscCall(ISGetLocalSize(is, &n));
1100:     issame = (PetscBool)(n == nlocal);
1101:     PetscCallMPI(MPIU_Allreduce(&issame, &isallsame, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1102:     issame = isallsame;
1103:   }
1104:   if (issame) {
1105:     PetscCall(ISDuplicate(is, isout));
1106:   } else {
1107:     PetscUseTypeMethod(is, invertpermutation, nlocal, isout);
1108:     PetscCall(ISSetPermutation(*isout));
1109:   }
1110:   PetscFunctionReturn(PETSC_SUCCESS);
1111: }

1113: /*@
1114:   ISGetSize - Returns the global length of an index set.

1116:   Not Collective

1118:   Input Parameter:
1119: . is - the index set

1121:   Output Parameter:
1122: . size - the global size

1124:   Level: beginner

1126: .seealso: `IS`
1127: @*/
1128: PetscErrorCode ISGetSize(IS is, PetscInt *size)
1129: {
1130:   PetscFunctionBegin;
1132:   PetscAssertPointer(size, 2);
1133:   *size = is->map->N;
1134:   PetscFunctionReturn(PETSC_SUCCESS);
1135: }

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

1140:   Not Collective

1142:   Input Parameter:
1143: . is - the index set

1145:   Output Parameter:
1146: . size - the local size

1148:   Level: beginner

1150: .seealso: `IS`, `ISGetSize()`
1151: @*/
1152: PetscErrorCode ISGetLocalSize(IS is, PetscInt *size)
1153: {
1154:   PetscFunctionBegin;
1156:   PetscAssertPointer(size, 2);
1157:   *size = is->map->n;
1158:   PetscFunctionReturn(PETSC_SUCCESS);
1159: }

1161: /*@
1162:   ISGetLayout - get `PetscLayout` describing index set layout

1164:   Not Collective

1166:   Input Parameter:
1167: . is - the index set

1169:   Output Parameter:
1170: . map - the layout

1172:   Level: developer

1174: .seealso: `IS`, `PetscLayout`, `ISSetLayout()`, `ISGetSize()`, `ISGetLocalSize()`
1175: @*/
1176: PetscErrorCode ISGetLayout(IS is, PetscLayout *map)
1177: {
1178:   PetscFunctionBegin;
1180:   PetscAssertPointer(map, 2);
1181:   *map = is->map;
1182:   PetscFunctionReturn(PETSC_SUCCESS);
1183: }

1185: /*@
1186:   ISSetLayout - set `PetscLayout` describing index set layout

1188:   Collective

1190:   Input Parameters:
1191: + is  - the index set
1192: - map - the layout

1194:   Level: developer

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

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

1202: .seealso: `IS`, `PetscLayout`, `ISCreate()`, `ISGetLayout()`, `ISGetSize()`, `ISGetLocalSize()`
1203: @*/
1204: PetscErrorCode ISSetLayout(IS is, PetscLayout map)
1205: {
1206:   PetscFunctionBegin;
1208:   PetscAssertPointer(map, 2);
1209:   PetscCall(PetscLayoutReference(map, &is->map));
1210:   PetscFunctionReturn(PETSC_SUCCESS);
1211: }

1213: /*@C
1214:   ISGetIndices - Returns a pointer to the indices.  The user should call
1215:   `ISRestoreIndices()` after having looked at the indices.  The user should
1216:   NOT change the indices.

1218:   Not Collective

1220:   Input Parameter:
1221: . is - the index set

1223:   Output Parameter:
1224: . ptr - the location to put the pointer to the indices

1226:   Level: intermediate

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

1231: .seealso: `IS`, `ISRestoreIndices()`, `ISGetIndicesF90()`
1232: @*/
1233: PetscErrorCode ISGetIndices(IS is, const PetscInt *ptr[])
1234: {
1235:   PetscFunctionBegin;
1237:   PetscAssertPointer(ptr, 2);
1238:   PetscUseTypeMethod(is, getindices, ptr);
1239:   PetscFunctionReturn(PETSC_SUCCESS);
1240: }

1242: /*@
1243:   ISGetMinMax - Gets the minimum and maximum values in an `IS`

1245:   Not Collective

1247:   Input Parameter:
1248: . is - the index set

1250:   Output Parameters:
1251: + min - the minimum value, you may pass `NULL`
1252: - max - the maximum value, you may pass `NULL`

1254:   Level: intermediate

1256:   Notes:
1257:   Empty index sets return min=`PETSC_INT_MAX` and max=`PETSC_INT_MIN`.

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

1261: .seealso: `IS`, `ISGetIndices()`, `ISRestoreIndices()`, `ISGetIndicesF90()`
1262: @*/
1263: PetscErrorCode ISGetMinMax(IS is, PetscInt *min, PetscInt *max)
1264: {
1265:   PetscFunctionBegin;
1267:   if (min) *min = is->min;
1268:   if (max) *max = is->max;
1269:   PetscFunctionReturn(PETSC_SUCCESS);
1270: }

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

1275:   Not Collective

1277:   Input Parameters:
1278: + is  - the index set
1279: - key - the search key

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

1284:   Level: intermediate

1286: .seealso: `IS`
1287:  @*/
1288: PetscErrorCode ISLocate(IS is, PetscInt key, PetscInt *location)
1289: {
1290:   PetscFunctionBegin;
1291:   if (is->ops->locate) {
1292:     PetscUseTypeMethod(is, locate, key, location);
1293:   } else {
1294:     PetscInt        numIdx;
1295:     PetscBool       sorted;
1296:     const PetscInt *idx;

1298:     PetscCall(ISGetLocalSize(is, &numIdx));
1299:     PetscCall(ISGetIndices(is, &idx));
1300:     PetscCall(ISSorted(is, &sorted));
1301:     if (sorted) {
1302:       PetscCall(PetscFindInt(key, numIdx, idx, location));
1303:     } else {
1304:       PetscInt i;

1306:       *location = -1;
1307:       for (i = 0; i < numIdx; i++) {
1308:         if (idx[i] == key) {
1309:           *location = i;
1310:           break;
1311:         }
1312:       }
1313:     }
1314:     PetscCall(ISRestoreIndices(is, &idx));
1315:   }
1316:   PetscFunctionReturn(PETSC_SUCCESS);
1317: }

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

1322:   Not Collective

1324:   Input Parameters:
1325: + is  - the index set
1326: - ptr - the pointer obtained by `ISGetIndices()`

1328:   Level: intermediate

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

1333: .seealso: `IS`, `ISGetIndices()`, `ISRestoreIndicesF90()`
1334: @*/
1335: PetscErrorCode ISRestoreIndices(IS is, const PetscInt *ptr[])
1336: {
1337:   PetscFunctionBegin;
1339:   PetscAssertPointer(ptr, 2);
1340:   PetscTryTypeMethod(is, restoreindices, ptr);
1341:   PetscFunctionReturn(PETSC_SUCCESS);
1342: }

1344: static PetscErrorCode ISGatherTotal_Private(IS is)
1345: {
1346:   PetscInt        i, n, N;
1347:   const PetscInt *lindices;
1348:   MPI_Comm        comm;
1349:   PetscMPIInt     rank, size, *sizes = NULL, *offsets = NULL, nn;

1351:   PetscFunctionBegin;

1354:   PetscCall(PetscObjectGetComm((PetscObject)is, &comm));
1355:   PetscCallMPI(MPI_Comm_size(comm, &size));
1356:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
1357:   PetscCall(ISGetLocalSize(is, &n));
1358:   PetscCall(PetscMalloc2(size, &sizes, size, &offsets));

1360:   PetscCall(PetscMPIIntCast(n, &nn));
1361:   PetscCallMPI(MPI_Allgather(&nn, 1, MPI_INT, sizes, 1, MPI_INT, comm));
1362:   offsets[0] = 0;
1363:   for (i = 1; i < size; ++i) offsets[i] = offsets[i - 1] + sizes[i - 1];
1364:   N = offsets[size - 1] + sizes[size - 1];

1366:   PetscCall(PetscMalloc1(N, &is->total));
1367:   PetscCall(ISGetIndices(is, &lindices));
1368:   PetscCallMPI(MPI_Allgatherv((void *)lindices, nn, MPIU_INT, is->total, sizes, offsets, MPIU_INT, comm));
1369:   PetscCall(ISRestoreIndices(is, &lindices));
1370:   is->local_offset = offsets[rank];
1371:   PetscCall(PetscFree2(sizes, offsets));
1372:   PetscFunctionReturn(PETSC_SUCCESS);
1373: }

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

1378:   Collective

1380:   Input Parameter:
1381: . is - the index set

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

1387:   Level: intermediate

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

1397: .seealso: `IS`, `ISRestoreTotalIndices()`, `ISGetNonlocalIndices()`, `ISGetSize()`
1398: @*/
1399: PetscErrorCode ISGetTotalIndices(IS is, const PetscInt *indices[])
1400: {
1401:   PetscMPIInt size;

1403:   PetscFunctionBegin;
1405:   PetscAssertPointer(indices, 2);
1406:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is), &size));
1407:   if (size == 1) {
1408:     PetscUseTypeMethod(is, getindices, indices);
1409:   } else {
1410:     if (!is->total) PetscCall(ISGatherTotal_Private(is));
1411:     *indices = is->total;
1412:   }
1413:   PetscFunctionReturn(PETSC_SUCCESS);
1414: }

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

1419:   Not Collective.

1421:   Input Parameters:
1422: + is      - the index set
1423: - indices - index array; must be the array obtained with `ISGetTotalIndices()`

1425:   Level: intermediate

1427: .seealso: `IS`, `ISGetNonlocalIndices()`
1428: @*/
1429: PetscErrorCode ISRestoreTotalIndices(IS is, const PetscInt *indices[])
1430: {
1431:   PetscMPIInt size;

1433:   PetscFunctionBegin;
1435:   PetscAssertPointer(indices, 2);
1436:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is), &size));
1437:   if (size == 1) {
1438:     PetscUseTypeMethod(is, restoreindices, indices);
1439:   } else {
1440:     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.");
1441:   }
1442:   PetscFunctionReturn(PETSC_SUCCESS);
1443: }

1445: /*@C
1446:   ISGetNonlocalIndices - Retrieve an array of indices from remote processors
1447:   in this communicator.

1449:   Collective

1451:   Input Parameter:
1452: . is - the index set

1454:   Output Parameter:
1455: . indices - indices with rank 0 indices first, and so on,  omitting
1456:              the current rank.  Total number of indices is the difference
1457:              total and local, obtained with `ISGetSize()` and `ISGetLocalSize()`,
1458:              respectively.

1460:   Level: intermediate

1462:   Notes:
1463:   Restore the indices using `ISRestoreNonlocalIndices()`.

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

1467: .seealso: `IS`, `ISGetTotalIndices()`, `ISRestoreNonlocalIndices()`, `ISGetSize()`, `ISGetLocalSize().`
1468: @*/
1469: PetscErrorCode ISGetNonlocalIndices(IS is, const PetscInt *indices[])
1470: {
1471:   PetscMPIInt size;
1472:   PetscInt    n, N;

1474:   PetscFunctionBegin;
1476:   PetscAssertPointer(indices, 2);
1477:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is), &size));
1478:   if (size == 1) *indices = NULL;
1479:   else {
1480:     if (!is->total) PetscCall(ISGatherTotal_Private(is));
1481:     PetscCall(ISGetLocalSize(is, &n));
1482:     PetscCall(ISGetSize(is, &N));
1483:     PetscCall(PetscMalloc1(N - n, &is->nonlocal));
1484:     PetscCall(PetscArraycpy(is->nonlocal, is->total, is->local_offset));
1485:     PetscCall(PetscArraycpy(is->nonlocal + is->local_offset, is->total + is->local_offset + n, N - is->local_offset - n));
1486:     *indices = is->nonlocal;
1487:   }
1488:   PetscFunctionReturn(PETSC_SUCCESS);
1489: }

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

1494:   Not Collective.

1496:   Input Parameters:
1497: + is      - the index set
1498: - indices - index array; must be the array obtained with `ISGetNonlocalIndices()`

1500:   Level: intermediate

1502: .seealso: `IS`, `ISGetTotalIndices()`, `ISGetNonlocalIndices()`, `ISRestoreTotalIndices()`
1503: @*/
1504: PetscErrorCode ISRestoreNonlocalIndices(IS is, const PetscInt *indices[])
1505: {
1506:   PetscFunctionBegin;
1508:   PetscAssertPointer(indices, 2);
1509:   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.");
1510:   PetscFunctionReturn(PETSC_SUCCESS);
1511: }

1513: /*@
1514:   ISGetNonlocalIS - Gather all nonlocal indices for this `IS` and present
1515:   them as another sequential index set.

1517:   Collective

1519:   Input Parameter:
1520: . is - the index set

1522:   Output Parameter:
1523: . complement - sequential `IS` with indices identical to the result of
1524:                 `ISGetNonlocalIndices()`

1526:   Level: intermediate

1528:   Notes:
1529:   Complement represents the result of `ISGetNonlocalIndices()` as an `IS`.
1530:   Therefore scalability issues similar to `ISGetNonlocalIndices()` apply.

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

1534: .seealso: `IS`, `ISGetNonlocalIndices()`, `ISRestoreNonlocalIndices()`, `ISAllGather()`, `ISGetSize()`
1535: @*/
1536: PetscErrorCode ISGetNonlocalIS(IS is, IS *complement)
1537: {
1538:   PetscFunctionBegin;
1540:   PetscAssertPointer(complement, 2);
1541:   /* Check if the complement exists already. */
1542:   if (is->complement) {
1543:     *complement = is->complement;
1544:     PetscCall(PetscObjectReference((PetscObject)is->complement));
1545:   } else {
1546:     PetscInt        N, n;
1547:     const PetscInt *idx;
1548:     PetscCall(ISGetSize(is, &N));
1549:     PetscCall(ISGetLocalSize(is, &n));
1550:     PetscCall(ISGetNonlocalIndices(is, &idx));
1551:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, N - n, idx, PETSC_USE_POINTER, &is->complement));
1552:     PetscCall(PetscObjectReference((PetscObject)is->complement));
1553:     *complement = is->complement;
1554:   }
1555:   PetscFunctionReturn(PETSC_SUCCESS);
1556: }

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

1561:   Not collective.

1563:   Input Parameters:
1564: + is         - the index set
1565: - complement - index set of `is`'s nonlocal indices

1567:   Level: intermediate

1569: .seealso: `IS`, `ISGetNonlocalIS()`, `ISGetNonlocalIndices()`, `ISRestoreNonlocalIndices()`
1570: @*/
1571: PetscErrorCode ISRestoreNonlocalIS(IS is, IS *complement)
1572: {
1573:   PetscInt refcnt;

1575:   PetscFunctionBegin;
1577:   PetscAssertPointer(complement, 2);
1578:   PetscCheck(*complement == is->complement, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Complement IS being restored was not obtained with ISGetNonlocalIS()");
1579:   PetscCall(PetscObjectGetReference((PetscObject)is->complement, &refcnt));
1580:   PetscCheck(refcnt > 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Duplicate call to ISRestoreNonlocalIS() detected");
1581:   PetscCall(PetscObjectDereference((PetscObject)is->complement));
1582:   PetscFunctionReturn(PETSC_SUCCESS);
1583: }

1585: /*@
1586:   ISViewFromOptions - View an `IS` based on options in the options database

1588:   Collective

1590:   Input Parameters:
1591: + A    - the index set
1592: . obj  - Optional object that provides the prefix for the options database
1593: - name - command line option

1595:   Level: intermediate

1597:   Note:
1598:   See `PetscObjectViewFromOptions()` for possible `PetscViewer` and `PetscViewerFormat` values

1600: .seealso: `IS`, `ISView()`, `PetscObjectViewFromOptions()`, `ISCreate()`
1601: @*/
1602: PetscErrorCode ISViewFromOptions(IS A, PetscObject obj, const char name[])
1603: {
1604:   PetscFunctionBegin;
1606:   PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
1607:   PetscFunctionReturn(PETSC_SUCCESS);
1608: }

1610: /*@
1611:   ISView - Displays an index set.

1613:   Collective

1615:   Input Parameters:
1616: + is     - the index set
1617: - viewer - viewer used to display the set, for example `PETSC_VIEWER_STDOUT_SELF`.

1619:   Level: intermediate

1621: .seealso: `IS`, `PetscViewer`, `PetscViewerASCIIOpen()`, `ISViewFromOptions()`
1622: @*/
1623: PetscErrorCode ISView(IS is, PetscViewer viewer)
1624: {
1625:   PetscFunctionBegin;
1627:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)is), &viewer));
1629:   PetscCheckSameComm(is, 1, viewer, 2);

1631:   PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)is, viewer));
1632:   PetscCall(PetscLogEventBegin(IS_View, is, viewer, 0, 0));
1633:   PetscUseTypeMethod(is, view, viewer);
1634:   PetscCall(PetscLogEventEnd(IS_View, is, viewer, 0, 0));
1635:   PetscFunctionReturn(PETSC_SUCCESS);
1636: }

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

1641:   Collective

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

1647:   Level: intermediate

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

1654: .seealso: `IS`, `PetscViewerBinaryOpen()`, `ISView()`, `MatLoad()`, `VecLoad()`
1655: @*/
1656: PetscErrorCode ISLoad(IS is, PetscViewer viewer)
1657: {
1658:   PetscBool isbinary, ishdf5;

1660:   PetscFunctionBegin;
1663:   PetscCheckSameComm(is, 1, viewer, 2);
1664:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1665:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1666:   PetscCheck(isbinary || ishdf5, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1667:   if (!((PetscObject)is)->type_name) PetscCall(ISSetType(is, ISGENERAL));
1668:   PetscCall(PetscLogEventBegin(IS_Load, is, viewer, 0, 0));
1669:   PetscUseTypeMethod(is, load, viewer);
1670:   PetscCall(PetscLogEventEnd(IS_Load, is, viewer, 0, 0));
1671:   PetscFunctionReturn(PETSC_SUCCESS);
1672: }

1674: /*@
1675:   ISSort - Sorts the indices of an index set.

1677:   Collective

1679:   Input Parameter:
1680: . is - the index set

1682:   Level: intermediate

1684: .seealso: `IS`, `ISSortRemoveDups()`, `ISSorted()`
1685: @*/
1686: PetscErrorCode ISSort(IS is)
1687: {
1688:   PetscBool flg;

1690:   PetscFunctionBegin;
1692:   PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_FALSE, &flg));
1693:   if (!flg) {
1694:     PetscUseTypeMethod(is, sort);
1695:     PetscCall(ISSetInfo(is, IS_SORTED, IS_LOCAL, is->info_permanent[IS_LOCAL][IS_SORTED], PETSC_TRUE));
1696:   }
1697:   PetscFunctionReturn(PETSC_SUCCESS);
1698: }

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

1703:   Collective

1705:   Input Parameter:
1706: . is - the index set

1708:   Level: intermediate

1710: .seealso: `IS`, `ISSort()`, `ISSorted()`
1711: @*/
1712: PetscErrorCode ISSortRemoveDups(IS is)
1713: {
1714:   PetscFunctionBegin;
1716:   PetscCall(ISClearInfoCache(is, PETSC_FALSE));
1717:   PetscUseTypeMethod(is, sortremovedups);
1718:   PetscCall(ISSetInfo(is, IS_SORTED, IS_LOCAL, is->info_permanent[IS_LOCAL][IS_SORTED], PETSC_TRUE));
1719:   PetscCall(ISSetInfo(is, IS_UNIQUE, IS_LOCAL, is->info_permanent[IS_LOCAL][IS_UNIQUE], PETSC_TRUE));
1720:   PetscFunctionReturn(PETSC_SUCCESS);
1721: }

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

1726:   Collective

1728:   Input Parameter:
1729: . is - the index set

1731:   Level: intermediate

1733: .seealso: `IS`, `ISSorted()`
1734: @*/
1735: PetscErrorCode ISToGeneral(IS is)
1736: {
1737:   PetscFunctionBegin;
1739:   PetscUseTypeMethod(is, togeneral);
1740:   PetscFunctionReturn(PETSC_SUCCESS);
1741: }

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

1746:   Not Collective

1748:   Input Parameter:
1749: . is - the index set

1751:   Output Parameter:
1752: . flg - output flag, either `PETSC_TRUE` if the index set is sorted,
1753:          or `PETSC_FALSE` otherwise.

1755:   Level: intermediate

1757:   Note:
1758:   For parallel IS objects this only indicates if the local part of `is`
1759:   is sorted. So some processors may return `PETSC_TRUE` while others may
1760:   return `PETSC_FALSE`.

1762: .seealso: `ISSort()`, `ISSortRemoveDups()`
1763: @*/
1764: PetscErrorCode ISSorted(IS is, PetscBool *flg)
1765: {
1766:   PetscFunctionBegin;
1768:   PetscAssertPointer(flg, 2);
1769:   PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, flg));
1770:   PetscFunctionReturn(PETSC_SUCCESS);
1771: }

1773: /*@
1774:   ISDuplicate - Creates a duplicate copy of an index set.

1776:   Collective

1778:   Input Parameter:
1779: . is - the index set

1781:   Output Parameter:
1782: . newIS - the copy of the index set

1784:   Level: beginner

1786: .seealso: `IS`, `ISCreateGeneral()`, `ISCopy()`
1787: @*/
1788: PetscErrorCode ISDuplicate(IS is, IS *newIS)
1789: {
1790:   PetscFunctionBegin;
1792:   PetscAssertPointer(newIS, 2);
1793:   PetscUseTypeMethod(is, duplicate, newIS);
1794:   PetscCall(ISCopyInfo_Private(is, *newIS));
1795:   PetscFunctionReturn(PETSC_SUCCESS);
1796: }

1798: /*@
1799:   ISCopy - Copies an index set.

1801:   Collective

1803:   Input Parameter:
1804: . is - the index set

1806:   Output Parameter:
1807: . isy - the copy of the index set

1809:   Level: beginner

1811: .seealso: `IS`, `ISDuplicate()`, `ISShift()`
1812: @*/
1813: PetscErrorCode ISCopy(IS is, IS isy)
1814: {
1815:   PetscInt bs, bsy;

1817:   PetscFunctionBegin;
1820:   PetscCheckSameComm(is, 1, isy, 2);
1821:   if (is == isy) PetscFunctionReturn(PETSC_SUCCESS);
1822:   PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
1823:   PetscCall(PetscLayoutGetBlockSize(isy->map, &bsy));
1824:   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);
1825:   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);
1826:   PetscCheck(bs == bsy, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Index sets have different block size %" PetscInt_FMT " != %" PetscInt_FMT, bs, bsy);
1827:   PetscCall(ISCopyInfo_Private(is, isy));
1828:   isy->max = is->max;
1829:   isy->min = is->min;
1830:   PetscUseTypeMethod(is, copy, isy);
1831:   PetscFunctionReturn(PETSC_SUCCESS);
1832: }

1834: /*@
1835:   ISShift - Shift all indices by given offset

1837:   Collective

1839:   Input Parameters:
1840: + is     - the index set
1841: - offset - the offset

1843:   Output Parameter:
1844: . isy - the shifted copy of the input index set

1846:   Level: beginner

1848:   Notes:
1849:   The `offset` can be different across processes.

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

1853: .seealso: `ISDuplicate()`, `ISCopy()`
1854: @*/
1855: PetscErrorCode ISShift(IS is, PetscInt offset, IS isy)
1856: {
1857:   PetscFunctionBegin;
1860:   PetscCheckSameComm(is, 1, isy, 3);
1861:   if (!offset) {
1862:     PetscCall(ISCopy(is, isy));
1863:     PetscFunctionReturn(PETSC_SUCCESS);
1864:   }
1865:   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);
1866:   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);
1867:   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);
1868:   PetscCall(ISCopyInfo_Private(is, isy));
1869:   isy->max = is->max + offset;
1870:   isy->min = is->min + offset;
1871:   PetscUseMethod(is, "ISShift_C", (IS, PetscInt, IS), (is, offset, isy));
1872:   PetscFunctionReturn(PETSC_SUCCESS);
1873: }

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

1878:   Collective

1880:   Input Parameters:
1881: + is   - index set
1882: . comm - communicator for new index set
1883: - mode - copy semantics, `PETSC_USE_POINTER` for no-copy if possible, otherwise `PETSC_COPY_VALUES`

1885:   Output Parameter:
1886: . newis - new `IS` on `comm`

1888:   Level: advanced

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

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

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

1898: .seealso: `IS`
1899: @*/
1900: PetscErrorCode ISOnComm(IS is, MPI_Comm comm, PetscCopyMode mode, IS *newis)
1901: {
1902:   PetscMPIInt match;

1904:   PetscFunctionBegin;
1906:   PetscAssertPointer(newis, 4);
1907:   PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)is), comm, &match));
1908:   if (mode != PETSC_COPY_VALUES && (match == MPI_IDENT || match == MPI_CONGRUENT)) {
1909:     PetscCall(PetscObjectReference((PetscObject)is));
1910:     *newis = is;
1911:   } else PetscUseTypeMethod(is, oncomm, comm, mode, newis);
1912:   PetscFunctionReturn(PETSC_SUCCESS);
1913: }

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

1918:   Logicall Collective

1920:   Input Parameters:
1921: + is - index set
1922: - bs - block size

1924:   Level: intermediate

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

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

1934: .seealso: `IS`, `ISGetBlockSize()`, `ISCreateBlock()`, `ISBlockGetIndices()`,
1935: @*/
1936: PetscErrorCode ISSetBlockSize(IS is, PetscInt bs)
1937: {
1938:   PetscFunctionBegin;
1941:   PetscCheck(bs >= 1, PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Block size %" PetscInt_FMT ", must be positive", bs);
1942:   if (PetscDefined(USE_DEBUG)) {
1943:     const PetscInt *indices;
1944:     PetscInt        length, i, j;
1945:     PetscCall(ISGetIndices(is, &indices));
1946:     if (indices) {
1947:       PetscCall(ISGetLocalSize(is, &length));
1948:       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);
1949:       for (i = 1; i < length / bs; i += bs) {
1950:         for (j = 1; j < bs - 1; j++) {
1951:           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);
1952:         }
1953:       }
1954:     }
1955:     PetscCall(ISRestoreIndices(is, &indices));
1956:   }
1957:   PetscUseTypeMethod(is, setblocksize, bs);
1958:   PetscFunctionReturn(PETSC_SUCCESS);
1959: }

1961: /*@
1962:   ISGetBlockSize - Returns the number of elements in a block.

1964:   Not Collective

1966:   Input Parameter:
1967: . is - the index set

1969:   Output Parameter:
1970: . size - the number of elements in a block

1972:   Level: intermediate

1974:   Note:
1975:   See `ISSetBlockSize()`

1977: .seealso: `IS`, `ISBlockGetSize()`, `ISGetSize()`, `ISCreateBlock()`, `ISSetBlockSize()`
1978: @*/
1979: PetscErrorCode ISGetBlockSize(IS is, PetscInt *size)
1980: {
1981:   PetscFunctionBegin;
1982:   PetscCall(PetscLayoutGetBlockSize(is->map, size));
1983:   PetscFunctionReturn(PETSC_SUCCESS);
1984: }

1986: static PetscErrorCode ISGetIndicesCopy_Private(IS is, PetscInt idx[])
1987: {
1988:   PetscInt        len, i;
1989:   const PetscInt *ptr;

1991:   PetscFunctionBegin;
1992:   PetscCall(ISGetLocalSize(is, &len));
1993:   PetscCall(ISGetIndices(is, &ptr));
1994:   for (i = 0; i < len; i++) idx[i] = ptr[i];
1995:   PetscCall(ISRestoreIndices(is, &ptr));
1996:   PetscFunctionReturn(PETSC_SUCCESS);
1997: }

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

2004:     Synopsis:
2005:     ISGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

2007:     Not Collective

2009:     Input Parameter:
2010: .   x - index set

2012:     Output Parameters:
2013: +   xx_v - the Fortran pointer to the array
2014: -   ierr - error code

2016:     Example of Usage:
2017: .vb
2018:     PetscInt, pointer xx_v(:)
2019:     ....
2020:     call ISGetIndicesF90(x,xx_v,ierr)
2021:     a = xx_v(3)
2022:     call ISRestoreIndicesF90(x,xx_v,ierr)
2023: .ve

2025:     Level: intermediate

2027: .seealso: `ISRestoreIndicesF90()`, `ISGetIndices()`, `ISRestoreIndices()`
2028: M*/

2030: /*MC
2031:     ISRestoreIndicesF90 - Restores an index set to a usable state after
2032:     a call to `ISGetIndicesF90()`.

2034:     Synopsis:
2035:     ISRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

2037:     Not Collective

2039:     Input Parameters:
2040: +   x - index set
2041: -   xx_v - the Fortran pointer to the array

2043:     Output Parameter:
2044: .   ierr - error code

2046:     Example of Usage:
2047: .vb
2048:     PetscInt, pointer xx_v(:)
2049:     ....
2050:     call ISGetIndicesF90(x,xx_v,ierr)
2051:     a = xx_v(3)
2052:     call ISRestoreIndicesF90(x,xx_v,ierr)
2053: .ve

2055:     Level: intermediate

2057: .seealso: `ISGetIndicesF90()`, `ISGetIndices()`, `ISRestoreIndices()`
2058: M*/

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

2065:     Synopsis:
2066:     ISBlockGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

2068:     Not Collective

2070:     Input Parameter:
2071: .   x - index set

2073:     Output Parameters:
2074: +   xx_v - the Fortran pointer to the array
2075: -   ierr - error code
2076:     Example of Usage:
2077: .vb
2078:     PetscInt, pointer xx_v(:)
2079:     ....
2080:     call ISBlockGetIndicesF90(x,xx_v,ierr)
2081:     a = xx_v(3)
2082:     call ISBlockRestoreIndicesF90(x,xx_v,ierr)
2083: .ve

2085:     Level: intermediate

2087: .seealso: `ISBlockRestoreIndicesF90()`, `ISGetIndices()`, `ISRestoreIndices()`,
2088:           `ISRestoreIndices()`
2089: M*/

2091: /*MC
2092:     ISBlockRestoreIndicesF90 - Restores an index set to a usable state after
2093:     a call to `ISBlockGetIndicesF90()`.

2095:     Synopsis:
2096:     ISBlockRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

2098:     Not Collective

2100:     Input Parameters:
2101: +   x - index set
2102: -   xx_v - the Fortran pointer to the array

2104:     Output Parameter:
2105: .   ierr - error code

2107:     Example of Usage:
2108: .vb
2109:     PetscInt, pointer xx_v(:)
2110:     ....
2111:     call ISBlockGetIndicesF90(x,xx_v,ierr)
2112:     a = xx_v(3)
2113:     call ISBlockRestoreIndicesF90(x,xx_v,ierr)
2114: .ve

2116:     Level: intermediate

2118: .seealso: `ISBlockGetIndicesF90()`, `ISGetIndices()`, `ISRestoreIndices()`, `ISRestoreIndicesF90()`
2119: M*/