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, bounds[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: bounds[0] = PETSC_INT_MAX;
63: bounds[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] < bounds[0]) bounds[0] = idxs[i];
70: if (idxs[i] > bounds[1]) bounds[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), bounds, bounds));
91: if (bounds[1] < bounds[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 = bounds[1] - bounds[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]) - bounds[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, MPI_C_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_rank(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, MPI_C_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, MPI_C_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_rank(comm, &rank));
549: if (type == IS_GLOBAL && is->ops->uniqueglobal) PetscUseTypeMethod(is, uniqueglobal, flg);
550: else {
551: PetscBool uniqueLocal;
552: PetscInt n = -1;
553: PetscInt *idx = NULL;
555: /* determine if the array is locally unique */
556: if (type == IS_GLOBAL && size > 1) {
557: /* call ISGetInfo so that a cached value will be used if possible */
558: PetscCall(ISGetInfo(is, IS_UNIQUE, IS_LOCAL, PETSC_TRUE, &uniqueLocal));
559: } else if (is->ops->uniquelocal) {
560: PetscUseTypeMethod(is, uniquelocal, &uniqueLocal);
561: } else {
562: /* default: get the local indices and directly check */
563: uniqueLocal = PETSC_TRUE;
564: PetscCall(ISGetLocalSize(is, &n));
565: PetscCall(PetscMalloc1(n, &idx));
566: PetscCall(ISGetIndicesCopy_Private(is, idx));
567: PetscCall(PetscIntSortSemiOrdered(n, idx));
568: for (i = 1; i < n; i++)
569: if (idx[i] == idx[i - 1]) break;
570: if (i < n) uniqueLocal = PETSC_FALSE;
571: }
573: PetscCall(PetscFree(idx));
574: if (type == IS_LOCAL || size == 1) {
575: *flg = uniqueLocal;
576: } else {
577: PetscCallMPI(MPIU_Allreduce(&uniqueLocal, flg, 1, MPI_C_BOOL, MPI_LAND, comm));
578: if (*flg) {
579: PetscInt min = PETSC_INT_MAX, max = PETSC_INT_MIN, maxprev;
581: if (!idx) {
582: PetscCall(ISGetLocalSize(is, &n));
583: PetscCall(PetscMalloc1(n, &idx));
584: PetscCall(ISGetIndicesCopy_Private(is, idx));
585: }
586: PetscCall(PetscParallelSortInt(is->map, is->map, idx, idx));
587: if (n) {
588: min = idx[0];
589: max = idx[n - 1];
590: }
591: for (i = 1; i < n; i++)
592: if (idx[i] == idx[i - 1]) break;
593: if (i < n) uniqueLocal = PETSC_FALSE;
594: maxprev = PETSC_INT_MIN;
595: PetscCallMPI(MPI_Exscan(&max, &maxprev, 1, MPIU_INT, MPI_MAX, comm));
596: if (rank && (maxprev == min)) uniqueLocal = PETSC_FALSE;
597: PetscCallMPI(MPIU_Allreduce(&uniqueLocal, flg, 1, MPI_C_BOOL, MPI_LAND, comm));
598: }
599: }
600: PetscCall(PetscFree(idx));
601: }
602: PetscFunctionReturn(PETSC_SUCCESS);
603: }
605: static PetscErrorCode ISGetInfo_Permutation(IS is, ISInfoType type, PetscBool *flg)
606: {
607: MPI_Comm comm;
608: PetscMPIInt size;
610: PetscFunctionBegin;
611: comm = PetscObjectComm((PetscObject)is);
612: PetscCallMPI(MPI_Comm_size(comm, &size));
613: if (type == IS_GLOBAL && is->ops->permglobal) {
614: PetscUseTypeMethod(is, permglobal, flg);
615: } else if (type == IS_LOCAL && is->ops->permlocal) {
616: PetscUseTypeMethod(is, permlocal, flg);
617: } else {
618: PetscBool permLocal;
619: PetscInt n, i, rStart;
620: PetscInt *idx;
622: PetscCall(ISGetLocalSize(is, &n));
623: PetscCall(PetscMalloc1(n, &idx));
624: PetscCall(ISGetIndicesCopy_Private(is, idx));
625: if (type == IS_GLOBAL) {
626: PetscCall(PetscParallelSortInt(is->map, is->map, idx, idx));
627: PetscCall(PetscLayoutGetRange(is->map, &rStart, NULL));
628: } else {
629: PetscCall(PetscIntSortSemiOrdered(n, idx));
630: rStart = 0;
631: }
632: permLocal = PETSC_TRUE;
633: for (i = 0; i < n; i++) {
634: if (idx[i] != rStart + i) break;
635: }
636: if (i < n) permLocal = PETSC_FALSE;
637: if (type == IS_LOCAL || size == 1) {
638: *flg = permLocal;
639: } else {
640: PetscCallMPI(MPIU_Allreduce(&permLocal, flg, 1, MPI_C_BOOL, MPI_LAND, comm));
641: }
642: PetscCall(PetscFree(idx));
643: }
644: PetscFunctionReturn(PETSC_SUCCESS);
645: }
647: static PetscErrorCode ISGetInfo_Interval(IS is, ISInfoType type, PetscBool *flg)
648: {
649: MPI_Comm comm;
650: PetscMPIInt size, rank;
651: PetscInt i;
653: PetscFunctionBegin;
654: comm = PetscObjectComm((PetscObject)is);
655: PetscCallMPI(MPI_Comm_size(comm, &size));
656: PetscCallMPI(MPI_Comm_rank(comm, &rank));
657: if (type == IS_GLOBAL && is->ops->intervalglobal) {
658: PetscUseTypeMethod(is, intervalglobal, flg);
659: } else {
660: PetscBool intervalLocal;
662: /* determine if the array is locally an interval */
663: if (type == IS_GLOBAL && size > 1) {
664: /* call ISGetInfo so that a cached value will be used if possible */
665: PetscCall(ISGetInfo(is, IS_INTERVAL, IS_LOCAL, PETSC_TRUE, &intervalLocal));
666: } else if (is->ops->intervallocal) {
667: PetscUseTypeMethod(is, intervallocal, &intervalLocal);
668: } else {
669: PetscInt n;
670: const PetscInt *idx;
671: /* default: get the local indices and directly check */
672: intervalLocal = PETSC_TRUE;
673: PetscCall(ISGetLocalSize(is, &n));
674: PetscCall(ISGetIndices(is, &idx));
675: for (i = 1; i < n; i++)
676: if (idx[i] != idx[i - 1] + 1) break;
677: if (i < n) intervalLocal = PETSC_FALSE;
678: PetscCall(ISRestoreIndices(is, &idx));
679: }
681: if (type == IS_LOCAL || size == 1) {
682: *flg = intervalLocal;
683: } else {
684: PetscCallMPI(MPIU_Allreduce(&intervalLocal, flg, 1, MPI_C_BOOL, MPI_LAND, comm));
685: if (*flg) {
686: PetscInt n, min = PETSC_INT_MAX, max = PETSC_INT_MIN;
687: PetscInt maxprev;
689: PetscCall(ISGetLocalSize(is, &n));
690: if (n) PetscCall(ISGetMinMax(is, &min, &max));
691: maxprev = PETSC_INT_MIN;
692: PetscCallMPI(MPI_Exscan(&max, &maxprev, 1, MPIU_INT, MPI_MAX, comm));
693: if (rank && n && (maxprev != min - 1)) intervalLocal = PETSC_FALSE;
694: PetscCallMPI(MPIU_Allreduce(&intervalLocal, flg, 1, MPI_C_BOOL, MPI_LAND, comm));
695: }
696: }
697: }
698: PetscFunctionReturn(PETSC_SUCCESS);
699: }
701: static PetscErrorCode ISGetInfo_Identity(IS is, ISInfoType type, PetscBool *flg)
702: {
703: MPI_Comm comm;
704: PetscMPIInt size;
706: PetscFunctionBegin;
707: comm = PetscObjectComm((PetscObject)is);
708: PetscCallMPI(MPI_Comm_size(comm, &size));
709: if (type == IS_GLOBAL && is->ops->intervalglobal) {
710: PetscBool isinterval;
712: PetscUseTypeMethod(is, intervalglobal, &isinterval);
713: *flg = PETSC_FALSE;
714: if (isinterval) {
715: PetscInt min;
717: PetscCall(ISGetMinMax(is, &min, NULL));
718: PetscCallMPI(MPI_Bcast(&min, 1, MPIU_INT, 0, comm));
719: if (min == 0) *flg = PETSC_TRUE;
720: }
721: } else if (type == IS_LOCAL && is->ops->intervallocal) {
722: PetscBool isinterval;
724: PetscUseTypeMethod(is, intervallocal, &isinterval);
725: *flg = PETSC_FALSE;
726: if (isinterval) {
727: PetscInt min;
729: PetscCall(ISGetMinMax(is, &min, NULL));
730: if (min == 0) *flg = PETSC_TRUE;
731: }
732: } else {
733: PetscBool identLocal;
734: PetscInt n, i, rStart;
735: const PetscInt *idx;
737: PetscCall(ISGetLocalSize(is, &n));
738: PetscCall(ISGetIndices(is, &idx));
739: PetscCall(PetscLayoutGetRange(is->map, &rStart, NULL));
740: identLocal = PETSC_TRUE;
741: for (i = 0; i < n; i++) {
742: if (idx[i] != rStart + i) break;
743: }
744: if (i < n) identLocal = PETSC_FALSE;
745: if (type == IS_LOCAL || size == 1) {
746: *flg = identLocal;
747: } else {
748: PetscCallMPI(MPIU_Allreduce(&identLocal, flg, 1, MPI_C_BOOL, MPI_LAND, comm));
749: }
750: PetscCall(ISRestoreIndices(is, &idx));
751: }
752: PetscFunctionReturn(PETSC_SUCCESS);
753: }
755: /*@
756: ISGetInfo - Determine whether an index set satisfies a given property
758: Collective or Logically Collective if the type is `IS_GLOBAL` (logically collective if the value of the property has been permanently set with `ISSetInfo()`)
760: Input Parameters:
761: + is - the index set
762: . info - describing a property of the index set, one of those listed in the documentation of `ISSetInfo()`
763: . 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
764: - type - whether the property is local (`IS_LOCAL`) or global (`IS_GLOBAL`)
766: Output Parameter:
767: . flg - whether the property is true (`PETSC_TRUE`) or false (`PETSC_FALSE`)
769: Level: advanced
771: Notes:
772: `ISGetInfo()` uses cached values when possible, which will be incorrect if `ISSetInfo()` has been called with incorrect information.
774: To clear cached values, use `ISClearInfoCache()`.
776: .seealso: `IS`, `ISInfo`, `ISInfoType`, `ISSetInfo()`, `ISClearInfoCache()`
777: @*/
778: PetscErrorCode ISGetInfo(IS is, ISInfo info, ISInfoType type, PetscBool compute, PetscBool *flg)
779: {
780: MPI_Comm comm, errcomm;
781: PetscMPIInt rank, size;
782: PetscInt itype;
783: PetscBool hasprop;
784: PetscBool infer;
786: PetscFunctionBegin;
789: comm = PetscObjectComm((PetscObject)is);
790: if (type == IS_GLOBAL) {
792: errcomm = comm;
793: } else {
794: errcomm = PETSC_COMM_SELF;
795: }
797: PetscCallMPI(MPI_Comm_size(comm, &size));
798: PetscCallMPI(MPI_Comm_rank(comm, &rank));
800: PetscCheck((int)info > IS_INFO_MIN && (int)info < IS_INFO_MAX, errcomm, PETSC_ERR_ARG_OUTOFRANGE, "Option %d is out of range", (int)info);
801: if (size == 1) type = IS_LOCAL;
802: itype = (type == IS_LOCAL) ? 0 : 1;
803: hasprop = PETSC_FALSE;
804: infer = PETSC_FALSE;
805: if (is->info_permanent[itype][(int)info]) {
806: hasprop = (is->info[itype][(int)info] == IS_INFO_TRUE) ? PETSC_TRUE : PETSC_FALSE;
807: infer = PETSC_TRUE;
808: } else if ((itype == IS_LOCAL) && (is->info[IS_LOCAL][info] != IS_INFO_UNKNOWN)) {
809: /* we can cache local properties as long as we clear them when the IS changes */
810: /* NOTE: we only cache local values because there is no ISAssemblyBegin()/ISAssemblyEnd(),
811: so we have no way of knowing when a cached value has been invalidated by changes on a different process */
812: hasprop = (is->info[itype][(int)info] == IS_INFO_TRUE) ? PETSC_TRUE : PETSC_FALSE;
813: infer = PETSC_TRUE;
814: } else if (compute) {
815: switch (info) {
816: case IS_SORTED:
817: PetscCall(ISGetInfo_Sorted_Private(is, type, &hasprop));
818: break;
819: case IS_UNIQUE:
820: PetscCall(ISGetInfo_Unique_Private(is, type, &hasprop));
821: break;
822: case IS_PERMUTATION:
823: PetscCall(ISGetInfo_Permutation(is, type, &hasprop));
824: break;
825: case IS_INTERVAL:
826: PetscCall(ISGetInfo_Interval(is, type, &hasprop));
827: break;
828: case IS_IDENTITY:
829: PetscCall(ISGetInfo_Identity(is, type, &hasprop));
830: break;
831: default:
832: SETERRQ(errcomm, PETSC_ERR_ARG_OUTOFRANGE, "Unknown IS property");
833: }
834: infer = PETSC_TRUE;
835: }
836: /* call ISSetInfo_Internal to keep all of the implications straight */
837: if (infer) PetscCall(ISSetInfo_Internal(is, info, type, IS_INFO_UNKNOWN, hasprop));
838: *flg = hasprop;
839: PetscFunctionReturn(PETSC_SUCCESS);
840: }
842: static PetscErrorCode ISCopyInfo_Private(IS source, IS dest)
843: {
844: PetscFunctionBegin;
845: PetscCall(PetscArraycpy(&dest->info[0], &source->info[0], 2));
846: PetscCall(PetscArraycpy(&dest->info_permanent[0], &source->info_permanent[0], 2));
847: PetscFunctionReturn(PETSC_SUCCESS);
848: }
850: /*@
851: ISIdentity - Determines whether index set is the identity mapping.
853: Collective
855: Input Parameter:
856: . is - the index set
858: Output Parameter:
859: . ident - `PETSC_TRUE` if an identity, else `PETSC_FALSE`
861: Level: intermediate
863: Note:
864: If `ISSetIdentity()` (or `ISSetInfo()` for a permanent property) has been called,
865: `ISIdentity()` will return its answer without communication between processes, but
866: otherwise the output ident will be computed from `ISGetInfo()`,
867: which may require synchronization on the communicator of `is`. To avoid this computation,
868: call `ISGetInfo()` directly with the compute flag set to `PETSC_FALSE`, and ident will be assumed false.
870: .seealso: `IS`, `ISSetIdentity()`, `ISGetInfo()`
871: @*/
872: PetscErrorCode ISIdentity(IS is, PetscBool *ident)
873: {
874: PetscFunctionBegin;
876: PetscAssertPointer(ident, 2);
877: PetscCall(ISGetInfo(is, IS_IDENTITY, IS_GLOBAL, PETSC_TRUE, ident));
878: PetscFunctionReturn(PETSC_SUCCESS);
879: }
881: /*@
882: ISSetIdentity - Informs the index set that it is an identity.
884: Logically Collective
886: Input Parameter:
887: . is - the index set
889: Level: intermediate
891: Notes:
892: `is` will be considered the identity permanently, even if indices have been changes (for example, with
893: `ISGeneralSetIndices()`). It's a good idea to only set this property if `is` will not change in the future.
895: To clear this property, use `ISClearInfoCache()`.
897: Developer Notes:
898: Some of these info routines have statements about values changing in the `IS`, this seems to contradict the fact that `IS` cannot be changed?
900: .seealso: `IS`, `ISIdentity()`, `ISSetInfo()`, `ISClearInfoCache()`
901: @*/
902: PetscErrorCode ISSetIdentity(IS is)
903: {
904: PetscFunctionBegin;
906: PetscCall(ISSetInfo(is, IS_IDENTITY, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE));
907: PetscFunctionReturn(PETSC_SUCCESS);
908: }
910: /*@
911: ISContiguousLocal - Locates an index set with contiguous range within a global range, if possible
913: Not Collective
915: Input Parameters:
916: + is - the index set
917: . gstart - global start
918: - gend - global end
920: Output Parameters:
921: + start - start of contiguous block, as an offset from `gstart`
922: - contig - `PETSC_TRUE` if the index set refers to contiguous entries on this process, else `PETSC_FALSE`
924: Level: developer
926: .seealso: `IS`, `ISGetLocalSize()`, `VecGetOwnershipRange()`
927: @*/
928: PetscErrorCode ISContiguousLocal(IS is, PetscInt gstart, PetscInt gend, PetscInt *start, PetscBool *contig)
929: {
930: PetscFunctionBegin;
932: PetscAssertPointer(start, 4);
933: PetscAssertPointer(contig, 5);
934: 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);
935: *start = -1;
936: *contig = PETSC_FALSE;
937: PetscTryTypeMethod(is, contiguous, gstart, gend, start, contig);
938: PetscFunctionReturn(PETSC_SUCCESS);
939: }
941: /*@
942: ISPermutation - `PETSC_TRUE` or `PETSC_FALSE` depending on whether the
943: index set has been declared to be a permutation.
945: Logically Collective
947: Input Parameter:
948: . is - the index set
950: Output Parameter:
951: . perm - `PETSC_TRUE` if a permutation, else `PETSC_FALSE`
953: Level: intermediate
955: Note:
956: If it is not already known that `is` is a permutation (if `ISSetPermutation()`
957: or `ISSetInfo()` has not been called), this routine will not attempt to compute
958: whether the index set is a permutation and will assume `perm` is `PETSC_FALSE`.
959: To compute the value when it is not already known, use `ISGetInfo()` with
960: the compute flag set to `PETSC_TRUE`.
962: Developer Notes:
963: Perhaps some of these routines should use the `PetscBool3` enum to return appropriate values
965: .seealso: `IS`, `ISSetPermutation()`, `ISGetInfo()`
966: @*/
967: PetscErrorCode ISPermutation(IS is, PetscBool *perm)
968: {
969: PetscFunctionBegin;
971: PetscAssertPointer(perm, 2);
972: PetscCall(ISGetInfo(is, IS_PERMUTATION, IS_GLOBAL, PETSC_FALSE, perm));
973: PetscFunctionReturn(PETSC_SUCCESS);
974: }
976: /*@
977: ISSetPermutation - Informs the index set that it is a permutation.
979: Logically Collective
981: Input Parameter:
982: . is - the index set
984: Level: intermediate
986: Notes:
987: `is` will be considered a permutation permanently, even if indices have been changes (for example, with
988: `ISGeneralSetIndices()`). It's a good idea to only set this property if `is` will not change in the future.
990: To clear this property, use `ISClearInfoCache()`.
992: The debug version of the libraries (./configure --with-debugging=1) checks if the
993: index set is actually a permutation. The optimized version just believes you.
995: .seealso: `IS`, `ISPermutation()`, `ISSetInfo()`, `ISClearInfoCache().`
996: @*/
997: PetscErrorCode ISSetPermutation(IS is)
998: {
999: PetscFunctionBegin;
1001: if (PetscDefined(USE_DEBUG)) {
1002: PetscMPIInt size;
1004: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is), &size));
1005: if (size == 1) {
1006: PetscInt i, n, *idx;
1007: const PetscInt *iidx;
1009: PetscCall(ISGetSize(is, &n));
1010: PetscCall(PetscMalloc1(n, &idx));
1011: PetscCall(ISGetIndices(is, &iidx));
1012: PetscCall(PetscArraycpy(idx, iidx, n));
1013: PetscCall(PetscIntSortSemiOrdered(n, idx));
1014: for (i = 0; i < n; i++) PetscCheck(idx[i] == i, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Index set is not a permutation");
1015: PetscCall(PetscFree(idx));
1016: PetscCall(ISRestoreIndices(is, &iidx));
1017: }
1018: }
1019: PetscCall(ISSetInfo(is, IS_PERMUTATION, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE));
1020: PetscFunctionReturn(PETSC_SUCCESS);
1021: }
1023: /*@
1024: ISDestroy - Destroys an index set.
1026: Collective
1028: Input Parameter:
1029: . is - the index set
1031: Level: beginner
1033: .seealso: `IS`, `ISCreateGeneral()`, `ISCreateStride()`, `ISCreateBlock()`
1034: @*/
1035: PetscErrorCode ISDestroy(IS *is)
1036: {
1037: PetscFunctionBegin;
1038: if (!*is) PetscFunctionReturn(PETSC_SUCCESS);
1040: if (--((PetscObject)*is)->refct > 0) {
1041: *is = NULL;
1042: PetscFunctionReturn(PETSC_SUCCESS);
1043: }
1044: if ((*is)->complement) {
1045: PetscInt refcnt;
1046: PetscCall(PetscObjectGetReference((PetscObject)((*is)->complement), &refcnt));
1047: PetscCheck(refcnt <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Nonlocal IS has not been restored");
1048: PetscCall(ISDestroy(&(*is)->complement));
1049: }
1050: PetscTryTypeMethod(*is, destroy);
1051: PetscCall(PetscLayoutDestroy(&(*is)->map));
1052: /* Destroy local representations of offproc data. */
1053: PetscCall(PetscFree((*is)->total));
1054: PetscCall(PetscFree((*is)->nonlocal));
1055: PetscCall(PetscHeaderDestroy(is));
1056: PetscFunctionReturn(PETSC_SUCCESS);
1057: }
1059: /*@
1060: ISInvertPermutation - Creates a new permutation that is the inverse of
1061: a given permutation.
1063: Collective
1065: Input Parameters:
1066: + is - the index set
1067: - nlocal - number of indices on this processor in result (ignored for 1 processor) or
1068: use `PETSC_DECIDE`
1070: Output Parameter:
1071: . isout - the inverse permutation
1073: Level: intermediate
1075: Note:
1076: For parallel index sets this does the complete parallel permutation, but the
1077: code is not efficient for huge index sets (10,000,000 indices).
1079: .seealso: `IS`, `ISGetInfo()`, `ISSetPermutation()`, `ISGetPermutation()`
1080: @*/
1081: PetscErrorCode ISInvertPermutation(IS is, PetscInt nlocal, IS *isout)
1082: {
1083: PetscBool isperm, isidentity, issame;
1085: PetscFunctionBegin;
1087: PetscAssertPointer(isout, 3);
1088: PetscCall(ISGetInfo(is, IS_PERMUTATION, IS_GLOBAL, PETSC_TRUE, &isperm));
1089: PetscCheck(isperm, PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_WRONG, "Not a permutation");
1090: PetscCall(ISGetInfo(is, IS_IDENTITY, IS_GLOBAL, PETSC_TRUE, &isidentity));
1091: issame = PETSC_FALSE;
1092: if (isidentity) {
1093: PetscInt n;
1094: PetscBool isallsame;
1096: PetscCall(ISGetLocalSize(is, &n));
1097: issame = (PetscBool)(n == nlocal);
1098: PetscCallMPI(MPIU_Allreduce(&issame, &isallsame, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1099: issame = isallsame;
1100: }
1101: if (issame) PetscCall(ISDuplicate(is, isout));
1102: else {
1103: PetscUseTypeMethod(is, invertpermutation, nlocal, isout);
1104: PetscCall(ISSetPermutation(*isout));
1105: }
1106: PetscFunctionReturn(PETSC_SUCCESS);
1107: }
1109: /*@
1110: ISGetSize - Returns the global length of an index set.
1112: Not Collective
1114: Input Parameter:
1115: . is - the index set
1117: Output Parameter:
1118: . size - the global size
1120: Level: beginner
1122: .seealso: `IS`
1123: @*/
1124: PetscErrorCode ISGetSize(IS is, PetscInt *size)
1125: {
1126: PetscFunctionBegin;
1128: PetscAssertPointer(size, 2);
1129: *size = is->map->N;
1130: PetscFunctionReturn(PETSC_SUCCESS);
1131: }
1133: /*@
1134: ISGetLocalSize - Returns the local (processor) length of an index set.
1136: Not Collective
1138: Input Parameter:
1139: . is - the index set
1141: Output Parameter:
1142: . size - the local size
1144: Level: beginner
1146: .seealso: `IS`, `ISGetSize()`
1147: @*/
1148: PetscErrorCode ISGetLocalSize(IS is, PetscInt *size)
1149: {
1150: PetscFunctionBegin;
1152: PetscAssertPointer(size, 2);
1153: *size = is->map->n;
1154: PetscFunctionReturn(PETSC_SUCCESS);
1155: }
1157: /*@
1158: ISGetLayout - get `PetscLayout` describing index set layout
1160: Not Collective
1162: Input Parameter:
1163: . is - the index set
1165: Output Parameter:
1166: . map - the layout
1168: Level: developer
1170: .seealso: `IS`, `PetscLayout`, `ISSetLayout()`, `ISGetSize()`, `ISGetLocalSize()`
1171: @*/
1172: PetscErrorCode ISGetLayout(IS is, PetscLayout *map)
1173: {
1174: PetscFunctionBegin;
1176: PetscAssertPointer(map, 2);
1177: *map = is->map;
1178: PetscFunctionReturn(PETSC_SUCCESS);
1179: }
1181: /*@
1182: ISSetLayout - set `PetscLayout` describing index set layout
1184: Collective
1186: Input Parameters:
1187: + is - the index set
1188: - map - the layout
1190: Level: developer
1192: Notes:
1193: Users should typically use higher level functions such as `ISCreateGeneral()`.
1195: This function can be useful in some special cases of constructing a new `IS`, e.g. after `ISCreate()` and before `ISLoad()`.
1196: Otherwise, it is only valid to replace the layout with a layout known to be equivalent.
1198: .seealso: `IS`, `PetscLayout`, `ISCreate()`, `ISGetLayout()`, `ISGetSize()`, `ISGetLocalSize()`
1199: @*/
1200: PetscErrorCode ISSetLayout(IS is, PetscLayout map)
1201: {
1202: PetscFunctionBegin;
1204: PetscAssertPointer(map, 2);
1205: PetscCall(PetscLayoutReference(map, &is->map));
1206: PetscFunctionReturn(PETSC_SUCCESS);
1207: }
1209: /*@C
1210: ISGetIndices - Returns a pointer to the indices. The user should call
1211: `ISRestoreIndices()` after having looked at the indices. The user should
1212: NOT change the indices.
1214: Not Collective
1216: Input Parameter:
1217: . is - the index set
1219: Output Parameter:
1220: . ptr - the location to put the pointer to the indices
1222: Level: intermediate
1224: Fortran Note:
1225: .vb
1226: PetscInt, pointer :: ptr(:)
1227: .ve
1229: .seealso: `IS`, `ISRestoreIndices()`
1230: @*/
1231: PetscErrorCode ISGetIndices(IS is, const PetscInt *ptr[])
1232: {
1233: PetscFunctionBegin;
1235: PetscAssertPointer(ptr, 2);
1236: PetscUseTypeMethod(is, getindices, ptr);
1237: PetscFunctionReturn(PETSC_SUCCESS);
1238: }
1240: /*@
1241: ISGetMinMax - Gets the minimum and maximum values in an `IS`
1243: Not Collective
1245: Input Parameter:
1246: . is - the index set
1248: Output Parameters:
1249: + min - the minimum value, you may pass `NULL`
1250: - max - the maximum value, you may pass `NULL`
1252: Level: intermediate
1254: Notes:
1255: Empty index sets return min=`PETSC_INT_MAX` and max=`PETSC_INT_MIN`.
1257: In parallel, it returns the `min` and `max` of the local portion of `is`
1259: .seealso: `IS`, `ISGetIndices()`, `ISRestoreIndices()`
1260: @*/
1261: PetscErrorCode ISGetMinMax(IS is, PetscInt *min, PetscInt *max)
1262: {
1263: PetscFunctionBegin;
1265: if (min) *min = is->min;
1266: if (max) *max = is->max;
1267: PetscFunctionReturn(PETSC_SUCCESS);
1268: }
1270: /*@
1271: ISLocate - determine the location of an index within the local component of an index set
1273: Not Collective
1275: Input Parameters:
1276: + is - the index set
1277: - key - the search key
1279: Output Parameter:
1280: . location - if >= 0, a location within the index set that is equal to the key, otherwise the key is not in the index set
1282: Level: intermediate
1284: .seealso: `IS`
1285: @*/
1286: PetscErrorCode ISLocate(IS is, PetscInt key, PetscInt *location)
1287: {
1288: PetscFunctionBegin;
1289: if (is->ops->locate) PetscUseTypeMethod(is, locate, key, location);
1290: else {
1291: PetscInt numIdx;
1292: PetscBool sorted;
1293: const PetscInt *idx;
1295: PetscCall(ISGetLocalSize(is, &numIdx));
1296: PetscCall(ISGetIndices(is, &idx));
1297: PetscCall(ISSorted(is, &sorted));
1298: if (sorted) {
1299: PetscCall(PetscFindInt(key, numIdx, idx, location));
1300: } else {
1301: PetscInt i;
1303: *location = -1;
1304: for (i = 0; i < numIdx; i++) {
1305: if (idx[i] == key) {
1306: *location = i;
1307: break;
1308: }
1309: }
1310: }
1311: PetscCall(ISRestoreIndices(is, &idx));
1312: }
1313: PetscFunctionReturn(PETSC_SUCCESS);
1314: }
1316: /*@C
1317: ISRestoreIndices - Restores an index set to a usable state after a call to `ISGetIndices()`.
1319: Not Collective
1321: Input Parameters:
1322: + is - the index set
1323: - ptr - the pointer obtained by `ISGetIndices()`
1325: Level: intermediate
1327: Fortran Note:
1328: .vb
1329: PetscInt, pointer :: ptr(:)
1330: .ve
1332: .seealso: `IS`, `ISGetIndices()`
1333: @*/
1334: PetscErrorCode ISRestoreIndices(IS is, const PetscInt *ptr[])
1335: {
1336: PetscFunctionBegin;
1338: PetscAssertPointer(ptr, 2);
1339: PetscTryTypeMethod(is, restoreindices, ptr);
1340: PetscFunctionReturn(PETSC_SUCCESS);
1341: }
1343: static PetscErrorCode ISGatherTotal_Private(IS is)
1344: {
1345: PetscInt i, n, N;
1346: const PetscInt *lindices;
1347: MPI_Comm comm;
1348: PetscMPIInt rank, size, *sizes = NULL, *offsets = NULL, nn;
1350: PetscFunctionBegin;
1353: PetscCall(PetscObjectGetComm((PetscObject)is, &comm));
1354: PetscCallMPI(MPI_Comm_size(comm, &size));
1355: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1356: PetscCall(ISGetLocalSize(is, &n));
1357: PetscCall(PetscMalloc2(size, &sizes, size, &offsets));
1359: PetscCall(PetscMPIIntCast(n, &nn));
1360: PetscCallMPI(MPI_Allgather(&nn, 1, MPI_INT, sizes, 1, MPI_INT, comm));
1361: offsets[0] = 0;
1362: for (i = 1; i < size; ++i) offsets[i] = offsets[i - 1] + sizes[i - 1];
1363: N = offsets[size - 1] + sizes[size - 1];
1365: PetscCall(PetscMalloc1(N, &is->total));
1366: PetscCall(ISGetIndices(is, &lindices));
1367: PetscCallMPI(MPI_Allgatherv((void *)lindices, nn, MPIU_INT, is->total, sizes, offsets, MPIU_INT, comm));
1368: PetscCall(ISRestoreIndices(is, &lindices));
1369: is->local_offset = offsets[rank];
1370: PetscCall(PetscFree2(sizes, offsets));
1371: PetscFunctionReturn(PETSC_SUCCESS);
1372: }
1374: /*@C
1375: ISGetTotalIndices - Retrieve an array containing all indices across the communicator.
1377: Collective
1379: Input Parameter:
1380: . is - the index set
1382: Output Parameter:
1383: . indices - total indices with rank 0 indices first, and so on; total array size is
1384: the same as returned with `ISGetSize()`.
1386: Level: intermediate
1388: Notes:
1389: this is potentially nonscalable, but depends on the size of the total index set
1390: and the size of the communicator. This may be feasible for index sets defined on
1391: subcommunicators, such that the set size does not grow with `PETSC_WORLD_COMM`.
1392: Note also that there is no way to tell where the local part of the indices starts
1393: (use `ISGetIndices()` and `ISGetNonlocalIndices()` to retrieve just the local and just
1394: the nonlocal part (complement), respectively).
1396: .seealso: `IS`, `ISRestoreTotalIndices()`, `ISGetNonlocalIndices()`, `ISGetSize()`
1397: @*/
1398: PetscErrorCode ISGetTotalIndices(IS is, const PetscInt *indices[])
1399: {
1400: PetscMPIInt size;
1402: PetscFunctionBegin;
1404: PetscAssertPointer(indices, 2);
1405: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is), &size));
1406: if (size == 1) PetscUseTypeMethod(is, getindices, indices);
1407: else {
1408: if (!is->total) PetscCall(ISGatherTotal_Private(is));
1409: *indices = is->total;
1410: }
1411: PetscFunctionReturn(PETSC_SUCCESS);
1412: }
1414: /*@C
1415: ISRestoreTotalIndices - Restore the index array obtained with `ISGetTotalIndices()`.
1417: Not Collective.
1419: Input Parameters:
1420: + is - the index set
1421: - indices - index array; must be the array obtained with `ISGetTotalIndices()`
1423: Level: intermediate
1425: .seealso: `IS`, `ISGetNonlocalIndices()`
1426: @*/
1427: PetscErrorCode ISRestoreTotalIndices(IS is, const PetscInt *indices[])
1428: {
1429: PetscMPIInt size;
1431: PetscFunctionBegin;
1433: PetscAssertPointer(indices, 2);
1434: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is), &size));
1435: if (size == 1) {
1436: PetscUseTypeMethod(is, restoreindices, indices);
1437: } else {
1438: 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.");
1439: }
1440: PetscFunctionReturn(PETSC_SUCCESS);
1441: }
1443: /*@C
1444: ISGetNonlocalIndices - Retrieve an array of indices from remote processors
1445: in this communicator.
1447: Collective
1449: Input Parameter:
1450: . is - the index set
1452: Output Parameter:
1453: . indices - indices with rank 0 indices first, and so on, omitting
1454: the current rank. Total number of indices is the difference
1455: total and local, obtained with `ISGetSize()` and `ISGetLocalSize()`,
1456: respectively.
1458: Level: intermediate
1460: Notes:
1461: Restore the indices using `ISRestoreNonlocalIndices()`.
1463: The same scalability considerations as those for `ISGetTotalIndices()` apply here.
1465: .seealso: `IS`, `ISGetTotalIndices()`, `ISRestoreNonlocalIndices()`, `ISGetSize()`, `ISGetLocalSize().`
1466: @*/
1467: PetscErrorCode ISGetNonlocalIndices(IS is, const PetscInt *indices[])
1468: {
1469: PetscMPIInt size;
1470: PetscInt n, N;
1472: PetscFunctionBegin;
1474: PetscAssertPointer(indices, 2);
1475: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is), &size));
1476: if (size == 1) *indices = NULL;
1477: else {
1478: if (!is->total) PetscCall(ISGatherTotal_Private(is));
1479: PetscCall(ISGetLocalSize(is, &n));
1480: PetscCall(ISGetSize(is, &N));
1481: PetscCall(PetscMalloc1(N - n, &is->nonlocal));
1482: PetscCall(PetscArraycpy(is->nonlocal, is->total, is->local_offset));
1483: PetscCall(PetscArraycpy(is->nonlocal + is->local_offset, is->total + is->local_offset + n, N - is->local_offset - n));
1484: *indices = is->nonlocal;
1485: }
1486: PetscFunctionReturn(PETSC_SUCCESS);
1487: }
1489: /*@C
1490: ISRestoreNonlocalIndices - Restore the index array obtained with `ISGetNonlocalIndices()`.
1492: Not Collective.
1494: Input Parameters:
1495: + is - the index set
1496: - indices - index array; must be the array obtained with `ISGetNonlocalIndices()`
1498: Level: intermediate
1500: .seealso: `IS`, `ISGetTotalIndices()`, `ISGetNonlocalIndices()`, `ISRestoreTotalIndices()`
1501: @*/
1502: PetscErrorCode ISRestoreNonlocalIndices(IS is, const PetscInt *indices[])
1503: {
1504: PetscFunctionBegin;
1506: PetscAssertPointer(indices, 2);
1507: 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.");
1508: PetscFunctionReturn(PETSC_SUCCESS);
1509: }
1511: /*@
1512: ISGetNonlocalIS - Gather all nonlocal indices for this `IS` and present
1513: them as another sequential index set.
1515: Collective
1517: Input Parameter:
1518: . is - the index set
1520: Output Parameter:
1521: . complement - sequential `IS` with indices identical to the result of
1522: `ISGetNonlocalIndices()`
1524: Level: intermediate
1526: Notes:
1527: Complement represents the result of `ISGetNonlocalIndices()` as an `IS`.
1528: Therefore scalability issues similar to `ISGetNonlocalIndices()` apply.
1530: The resulting `IS` must be restored using `ISRestoreNonlocalIS()`.
1532: .seealso: `IS`, `ISGetNonlocalIndices()`, `ISRestoreNonlocalIndices()`, `ISAllGather()`, `ISGetSize()`
1533: @*/
1534: PetscErrorCode ISGetNonlocalIS(IS is, IS *complement)
1535: {
1536: PetscFunctionBegin;
1538: PetscAssertPointer(complement, 2);
1539: /* Check if the complement exists already. */
1540: if (is->complement) {
1541: *complement = is->complement;
1542: PetscCall(PetscObjectReference((PetscObject)is->complement));
1543: } else {
1544: PetscInt N, n;
1545: const PetscInt *idx;
1546: PetscCall(ISGetSize(is, &N));
1547: PetscCall(ISGetLocalSize(is, &n));
1548: PetscCall(ISGetNonlocalIndices(is, &idx));
1549: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, N - n, idx, PETSC_USE_POINTER, &is->complement));
1550: PetscCall(PetscObjectReference((PetscObject)is->complement));
1551: *complement = is->complement;
1552: }
1553: PetscFunctionReturn(PETSC_SUCCESS);
1554: }
1556: /*@
1557: ISRestoreNonlocalIS - Restore the `IS` obtained with `ISGetNonlocalIS()`.
1559: Not collective.
1561: Input Parameters:
1562: + is - the index set
1563: - complement - index set of `is`'s nonlocal indices
1565: Level: intermediate
1567: .seealso: `IS`, `ISGetNonlocalIS()`, `ISGetNonlocalIndices()`, `ISRestoreNonlocalIndices()`
1568: @*/
1569: PetscErrorCode ISRestoreNonlocalIS(IS is, IS *complement)
1570: {
1571: PetscInt refcnt;
1573: PetscFunctionBegin;
1575: PetscAssertPointer(complement, 2);
1576: PetscCheck(*complement == is->complement, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Complement IS being restored was not obtained with ISGetNonlocalIS()");
1577: PetscCall(PetscObjectGetReference((PetscObject)is->complement, &refcnt));
1578: PetscCheck(refcnt > 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Duplicate call to ISRestoreNonlocalIS() detected");
1579: PetscCall(PetscObjectDereference((PetscObject)is->complement));
1580: PetscFunctionReturn(PETSC_SUCCESS);
1581: }
1583: /*@
1584: ISViewFromOptions - View an `IS` based on options in the options database
1586: Collective
1588: Input Parameters:
1589: + A - the index set
1590: . obj - Optional object that provides the prefix for the options database
1591: - name - command line option
1593: Level: intermediate
1595: Note:
1596: See `PetscObjectViewFromOptions()` for possible `PetscViewer` and `PetscViewerFormat` values
1598: .seealso: `IS`, `ISView()`, `PetscObjectViewFromOptions()`, `ISCreate()`
1599: @*/
1600: PetscErrorCode ISViewFromOptions(IS A, PetscObject obj, const char name[])
1601: {
1602: PetscFunctionBegin;
1604: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
1605: PetscFunctionReturn(PETSC_SUCCESS);
1606: }
1608: /*@
1609: ISView - Displays an index set.
1611: Collective
1613: Input Parameters:
1614: + is - the index set
1615: - viewer - viewer used to display the set, for example `PETSC_VIEWER_STDOUT_SELF`.
1617: Level: intermediate
1619: .seealso: `IS`, `PetscViewer`, `PetscViewerASCIIOpen()`, `ISViewFromOptions()`
1620: @*/
1621: PetscErrorCode ISView(IS is, PetscViewer viewer)
1622: {
1623: PetscFunctionBegin;
1625: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)is), &viewer));
1627: PetscCheckSameComm(is, 1, viewer, 2);
1629: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)is, viewer));
1630: PetscCall(PetscLogEventBegin(IS_View, is, viewer, 0, 0));
1631: PetscUseTypeMethod(is, view, viewer);
1632: PetscCall(PetscLogEventEnd(IS_View, is, viewer, 0, 0));
1633: PetscFunctionReturn(PETSC_SUCCESS);
1634: }
1636: /*@
1637: ISLoad - Loads an index set that has been stored in binary or HDF5 format with `ISView()`.
1639: Collective
1641: Input Parameters:
1642: + is - the newly loaded index set, this needs to have been created with `ISCreate()` or some related function before a call to `ISLoad()`.
1643: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()` or HDF5 file viewer, obtained from `PetscViewerHDF5Open()`
1645: Level: intermediate
1647: Notes:
1648: IF using HDF5, you must assign the IS the same name as was used in `is`
1649: that was stored in the file using `PetscObjectSetName()`. Otherwise you will
1650: get the error message: "Cannot H5DOpen2() with Vec name NAMEOFOBJECT"
1652: .seealso: `IS`, `PetscViewerBinaryOpen()`, `ISView()`, `MatLoad()`, `VecLoad()`
1653: @*/
1654: PetscErrorCode ISLoad(IS is, PetscViewer viewer)
1655: {
1656: PetscBool isbinary, ishdf5;
1658: PetscFunctionBegin;
1661: PetscCheckSameComm(is, 1, viewer, 2);
1662: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1663: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1664: PetscCheck(isbinary || ishdf5, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1665: if (!((PetscObject)is)->type_name) PetscCall(ISSetType(is, ISGENERAL));
1666: PetscCall(PetscLogEventBegin(IS_Load, is, viewer, 0, 0));
1667: PetscUseTypeMethod(is, load, viewer);
1668: PetscCall(PetscLogEventEnd(IS_Load, is, viewer, 0, 0));
1669: PetscFunctionReturn(PETSC_SUCCESS);
1670: }
1672: /*@
1673: ISSort - Sorts the indices of an index set.
1675: Collective
1677: Input Parameter:
1678: . is - the index set
1680: Level: intermediate
1682: .seealso: `IS`, `ISSortRemoveDups()`, `ISSorted()`
1683: @*/
1684: PetscErrorCode ISSort(IS is)
1685: {
1686: PetscBool flg;
1688: PetscFunctionBegin;
1690: PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_FALSE, &flg));
1691: if (!flg) {
1692: PetscUseTypeMethod(is, sort);
1693: PetscCall(ISSetInfo(is, IS_SORTED, IS_LOCAL, is->info_permanent[IS_LOCAL][IS_SORTED], PETSC_TRUE));
1694: }
1695: PetscFunctionReturn(PETSC_SUCCESS);
1696: }
1698: /*@
1699: ISSortRemoveDups - Sorts the indices of an index set, removing duplicates.
1701: Collective
1703: Input Parameter:
1704: . is - the index set
1706: Level: intermediate
1708: .seealso: `IS`, `ISSort()`, `ISSorted()`
1709: @*/
1710: PetscErrorCode ISSortRemoveDups(IS is)
1711: {
1712: PetscFunctionBegin;
1714: PetscCall(ISClearInfoCache(is, PETSC_FALSE));
1715: PetscUseTypeMethod(is, sortremovedups);
1716: PetscCall(ISSetInfo(is, IS_SORTED, IS_LOCAL, is->info_permanent[IS_LOCAL][IS_SORTED], PETSC_TRUE));
1717: PetscCall(ISSetInfo(is, IS_UNIQUE, IS_LOCAL, is->info_permanent[IS_LOCAL][IS_UNIQUE], PETSC_TRUE));
1718: PetscFunctionReturn(PETSC_SUCCESS);
1719: }
1721: /*@
1722: ISToGeneral - Converts an IS object of any type to `ISGENERAL` type
1724: Collective
1726: Input Parameter:
1727: . is - the index set
1729: Level: intermediate
1731: .seealso: `IS`, `ISSorted()`
1732: @*/
1733: PetscErrorCode ISToGeneral(IS is)
1734: {
1735: PetscFunctionBegin;
1737: PetscUseTypeMethod(is, togeneral);
1738: PetscFunctionReturn(PETSC_SUCCESS);
1739: }
1741: /*@
1742: ISSorted - Checks the indices to determine whether they have been sorted.
1744: Not Collective
1746: Input Parameter:
1747: . is - the index set
1749: Output Parameter:
1750: . flg - output flag, either `PETSC_TRUE` if the index set is sorted,
1751: or `PETSC_FALSE` otherwise.
1753: Level: intermediate
1755: Note:
1756: For parallel IS objects this only indicates if the local part of `is`
1757: is sorted. So some processors may return `PETSC_TRUE` while others may
1758: return `PETSC_FALSE`.
1760: .seealso: `ISSort()`, `ISSortRemoveDups()`
1761: @*/
1762: PetscErrorCode ISSorted(IS is, PetscBool *flg)
1763: {
1764: PetscFunctionBegin;
1766: PetscAssertPointer(flg, 2);
1767: PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, flg));
1768: PetscFunctionReturn(PETSC_SUCCESS);
1769: }
1771: /*@
1772: ISDuplicate - Creates a duplicate copy of an index set.
1774: Collective
1776: Input Parameter:
1777: . is - the index set
1779: Output Parameter:
1780: . newIS - the copy of the index set
1782: Level: beginner
1784: .seealso: `IS`, `ISCreateGeneral()`, `ISCopy()`
1785: @*/
1786: PetscErrorCode ISDuplicate(IS is, IS *newIS)
1787: {
1788: PetscFunctionBegin;
1790: PetscAssertPointer(newIS, 2);
1791: PetscUseTypeMethod(is, duplicate, newIS);
1792: PetscCall(ISCopyInfo_Private(is, *newIS));
1793: PetscFunctionReturn(PETSC_SUCCESS);
1794: }
1796: /*@
1797: ISCopy - Copies an index set.
1799: Collective
1801: Input Parameter:
1802: . is - the index set
1804: Output Parameter:
1805: . isy - the copy of the index set
1807: Level: beginner
1809: .seealso: `IS`, `ISDuplicate()`, `ISShift()`
1810: @*/
1811: PetscErrorCode ISCopy(IS is, IS isy)
1812: {
1813: PetscInt bs, bsy;
1815: PetscFunctionBegin;
1818: PetscCheckSameComm(is, 1, isy, 2);
1819: if (is == isy) PetscFunctionReturn(PETSC_SUCCESS);
1820: PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
1821: PetscCall(PetscLayoutGetBlockSize(isy->map, &bsy));
1822: 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);
1823: 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);
1824: PetscCheck(bs == bsy, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Index sets have different block size %" PetscInt_FMT " != %" PetscInt_FMT, bs, bsy);
1825: PetscCall(ISCopyInfo_Private(is, isy));
1826: isy->max = is->max;
1827: isy->min = is->min;
1828: PetscUseTypeMethod(is, copy, isy);
1829: PetscFunctionReturn(PETSC_SUCCESS);
1830: }
1832: /*@
1833: ISShift - Shift all indices by given offset
1835: Collective
1837: Input Parameters:
1838: + is - the index set
1839: - offset - the offset
1841: Output Parameter:
1842: . isy - the shifted copy of the input index set
1844: Level: beginner
1846: Notes:
1847: The `offset` can be different across processes.
1849: `is` and `isy` can be the same.
1851: .seealso: `ISDuplicate()`, `ISCopy()`
1852: @*/
1853: PetscErrorCode ISShift(IS is, PetscInt offset, IS isy)
1854: {
1855: PetscFunctionBegin;
1858: PetscCheckSameComm(is, 1, isy, 3);
1859: if (!offset) {
1860: PetscCall(ISCopy(is, isy));
1861: PetscFunctionReturn(PETSC_SUCCESS);
1862: }
1863: 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);
1864: 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);
1865: 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);
1866: PetscCall(ISCopyInfo_Private(is, isy));
1867: if (offset > 0) {
1868: isy->max = is->max < PETSC_INT_MAX - offset ? is->max + offset : PETSC_INT_MAX;
1869: isy->min = is->min < PETSC_INT_MAX - offset ? is->min + offset : PETSC_INT_MAX;
1870: } else {
1871: isy->max = is->max > PETSC_INT_MIN - offset ? is->max + offset : PETSC_INT_MIN;
1872: isy->min = is->min > PETSC_INT_MIN - offset ? is->min + offset : PETSC_INT_MIN;
1873: }
1874: PetscUseMethod(is, "ISShift_C", (IS, PetscInt, IS), (is, offset, isy));
1875: PetscFunctionReturn(PETSC_SUCCESS);
1876: }
1878: /*@
1879: ISOnComm - Split a parallel `IS` on subcomms (usually self) or concatenate index sets on subcomms into a parallel index set
1881: Collective
1883: Input Parameters:
1884: + is - index set
1885: . comm - communicator for new index set
1886: - mode - copy semantics, `PETSC_USE_POINTER` for no-copy if possible, otherwise `PETSC_COPY_VALUES`
1888: Output Parameter:
1889: . newis - new `IS` on `comm`
1891: Level: advanced
1893: Notes:
1894: It is usually desirable to create a parallel `IS` and look at the local part when necessary.
1896: This function is useful if serial `IS`s must be created independently, or to view many
1897: logically independent serial `IS`s.
1899: The input `IS` must have the same type on every MPI process.
1901: .seealso: `IS`
1902: @*/
1903: PetscErrorCode ISOnComm(IS is, MPI_Comm comm, PetscCopyMode mode, IS *newis)
1904: {
1905: PetscMPIInt match;
1907: PetscFunctionBegin;
1909: PetscAssertPointer(newis, 4);
1910: PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)is), comm, &match));
1911: if (mode != PETSC_COPY_VALUES && (match == MPI_IDENT || match == MPI_CONGRUENT)) {
1912: PetscCall(PetscObjectReference((PetscObject)is));
1913: *newis = is;
1914: } else PetscUseTypeMethod(is, oncomm, comm, mode, newis);
1915: PetscFunctionReturn(PETSC_SUCCESS);
1916: }
1918: /*@
1919: ISSetBlockSize - informs an index set that it has a given block size
1921: Logicall Collective
1923: Input Parameters:
1924: + is - index set
1925: - bs - block size
1927: Level: intermediate
1929: Notes:
1930: This is much like the block size for `Vec`s. It indicates that one can think of the indices as
1931: being in a collection of equal size blocks. For `ISBLOCK` these collections of blocks are all contiguous
1932: 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
1933: have a block size of three set.
1935: `ISBlockGetIndices()` only works for `ISBLOCK`, not others.
1937: .seealso: `IS`, `ISGetBlockSize()`, `ISCreateBlock()`, `ISBlockGetIndices()`
1938: @*/
1939: PetscErrorCode ISSetBlockSize(IS is, PetscInt bs)
1940: {
1941: PetscFunctionBegin;
1944: PetscCheck(bs >= 1, PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Block size %" PetscInt_FMT ", must be positive", bs);
1945: if (PetscDefined(USE_DEBUG)) {
1946: const PetscInt *indices;
1947: PetscInt length, i, j;
1948: PetscCall(ISGetIndices(is, &indices));
1949: if (indices) {
1950: PetscCall(ISGetLocalSize(is, &length));
1951: 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);
1952: for (i = 1; i < length / bs; i += bs) {
1953: for (j = 1; j < bs - 1; j++) {
1954: 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);
1955: }
1956: }
1957: }
1958: PetscCall(ISRestoreIndices(is, &indices));
1959: }
1960: PetscUseTypeMethod(is, setblocksize, bs);
1961: PetscFunctionReturn(PETSC_SUCCESS);
1962: }
1964: /*@
1965: ISGetBlockSize - Returns the number of elements in a block.
1967: Not Collective
1969: Input Parameter:
1970: . is - the index set
1972: Output Parameter:
1973: . size - the number of elements in a block
1975: Level: intermediate
1977: Note:
1978: See `ISSetBlockSize()`
1980: .seealso: `IS`, `ISBlockGetSize()`, `ISGetSize()`, `ISCreateBlock()`, `ISSetBlockSize()`
1981: @*/
1982: PetscErrorCode ISGetBlockSize(IS is, PetscInt *size)
1983: {
1984: PetscFunctionBegin;
1985: PetscCall(PetscLayoutGetBlockSize(is->map, size));
1986: PetscFunctionReturn(PETSC_SUCCESS);
1987: }
1989: /*@
1990: ISSetCompressOutput - set the flag for output compression
1992: Logicall Collective
1994: Input Parameters:
1995: + is - index set
1996: - compress - flag for output compression
1998: Level: intermediate
2000: .seealso: `IS`, `ISGetCompressOutput()`, `ISView()`
2001: @*/
2002: PetscErrorCode ISSetCompressOutput(IS is, PetscBool compress)
2003: {
2004: PetscFunctionBegin;
2007: is->compressOutput = compress;
2008: PetscFunctionReturn(PETSC_SUCCESS);
2009: }
2011: /*@
2012: ISGetCompressOutput - Returns the flag for output compression
2014: Not Collective
2016: Input Parameter:
2017: . is - the index set
2019: Output Parameter:
2020: . compress - the flag to compress output
2022: Level: intermediate
2024: .seealso: `IS`, `ISSetCompressOutput()`, `ISView()`
2025: @*/
2026: PetscErrorCode ISGetCompressOutput(IS is, PetscBool *compress)
2027: {
2028: PetscFunctionBegin;
2030: PetscAssertPointer(compress, 2);
2031: *compress = is->compressOutput;
2032: PetscFunctionReturn(PETSC_SUCCESS);
2033: }
2035: static PetscErrorCode ISGetIndicesCopy_Private(IS is, PetscInt idx[])
2036: {
2037: PetscInt len, i;
2038: const PetscInt *ptr;
2040: PetscFunctionBegin;
2041: PetscCall(ISGetLocalSize(is, &len));
2042: PetscCall(ISGetIndices(is, &ptr));
2043: for (i = 0; i < len; i++) idx[i] = ptr[i];
2044: PetscCall(ISRestoreIndices(is, &ptr));
2045: PetscFunctionReturn(PETSC_SUCCESS);
2046: }