Actual source code: index.c
1: /*
2: Defines the abstract operations on index sets, i.e. the public interface.
3: */
4: #include <petsc/private/isimpl.h>
5: #include <petscviewer.h>
6: #include <petscsf.h>
8: /* Logging support */
9: PetscClassId IS_CLASSID;
10: /* TODO: Much more events are missing! */
11: PetscLogEvent IS_View;
12: PetscLogEvent IS_Load;
14: /*@
15: ISRenumber - Renumbers the non-negative entries of an index set in a contiguous way, starting from 0.
17: Collective
19: Input Parameters:
20: + subset - the index set
21: - subset_mult - the multiplicity of each entry in subset (optional, can be `NULL`)
23: Output Parameters:
24: + N - one past the largest entry of the new `IS`
25: - subset_n - the new `IS`
27: Level: intermediate
29: Note:
30: All negative entries are mapped to -1. Indices with non positive multiplicities are skipped.
32: .seealso: `IS`
33: @*/
34: PetscErrorCode ISRenumber(IS subset, IS subset_mult, PetscInt *N, IS *subset_n)
35: {
36: PetscSF sf;
37: PetscLayout map;
38: const PetscInt *idxs, *idxs_mult = NULL;
39: PetscInt *leaf_data, *root_data, *gidxs, *ilocal, *ilocalneg;
40: PetscInt N_n, n, i, lbounds[2], gbounds[2], Nl, ibs;
41: PetscInt n_n, nlocals, start, first_index, npos, nneg;
42: PetscMPIInt commsize;
43: PetscBool first_found, isblock;
45: PetscFunctionBegin;
48: if (N) PetscAssertPointer(N, 3);
49: else if (!subset_n) PetscFunctionReturn(PETSC_SUCCESS);
50: PetscCall(ISGetLocalSize(subset, &n));
51: if (subset_mult) {
52: PetscCall(ISGetLocalSize(subset_mult, &i));
53: PetscCheck(i == n, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Local subset and multiplicity sizes don't match! %" PetscInt_FMT " != %" PetscInt_FMT, n, i);
54: }
55: /* create workspace layout for computing global indices of subset */
56: PetscCall(PetscMalloc1(n, &ilocal));
57: PetscCall(PetscMalloc1(n, &ilocalneg));
58: PetscCall(ISGetIndices(subset, &idxs));
59: PetscCall(ISGetBlockSize(subset, &ibs));
60: PetscCall(PetscObjectTypeCompare((PetscObject)subset, ISBLOCK, &isblock));
61: if (subset_mult) PetscCall(ISGetIndices(subset_mult, &idxs_mult));
62: lbounds[0] = PETSC_MAX_INT;
63: lbounds[1] = PETSC_MIN_INT;
64: for (i = 0, npos = 0, nneg = 0; i < n; i++) {
65: if (idxs[i] < 0) {
66: ilocalneg[nneg++] = i;
67: continue;
68: }
69: if (idxs[i] < lbounds[0]) lbounds[0] = idxs[i];
70: if (idxs[i] > lbounds[1]) lbounds[1] = idxs[i];
71: ilocal[npos++] = i;
72: }
73: if (npos == n) {
74: PetscCall(PetscFree(ilocal));
75: PetscCall(PetscFree(ilocalneg));
76: }
78: /* create sf : leaf_data == multiplicity of indexes, root data == global index in layout */
79: PetscCall(PetscMalloc1(n, &leaf_data));
80: for (i = 0; i < n; i++) leaf_data[i] = idxs_mult ? PetscMax(idxs_mult[i], 0) : 1;
82: /* local size of new subset */
83: n_n = 0;
84: for (i = 0; i < n; i++) n_n += leaf_data[i];
85: if (ilocalneg)
86: for (i = 0; i < nneg; i++) leaf_data[ilocalneg[i]] = 0;
87: PetscCall(PetscFree(ilocalneg));
88: PetscCall(PetscMalloc1(PetscMax(n_n, n), &gidxs)); /* allocating extra space to reuse gidxs */
89: /* check for early termination (all negative) */
90: PetscCall(PetscGlobalMinMaxInt(PetscObjectComm((PetscObject)subset), lbounds, gbounds));
91: if (gbounds[1] < gbounds[0]) {
92: if (N) *N = 0;
93: if (subset_n) { /* all negative */
94: for (i = 0; i < n_n; i++) gidxs[i] = -1;
95: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)subset), n_n, gidxs, PETSC_COPY_VALUES, subset_n));
96: }
97: PetscCall(PetscFree(leaf_data));
98: PetscCall(PetscFree(gidxs));
99: PetscCall(ISRestoreIndices(subset, &idxs));
100: if (subset_mult) PetscCall(ISRestoreIndices(subset_mult, &idxs_mult));
101: PetscCall(PetscFree(ilocal));
102: PetscCall(PetscFree(ilocalneg));
103: PetscFunctionReturn(PETSC_SUCCESS);
104: }
106: /* split work */
107: N_n = gbounds[1] - gbounds[0] + 1;
108: PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)subset), &map));
109: PetscCall(PetscLayoutSetBlockSize(map, 1));
110: PetscCall(PetscLayoutSetSize(map, N_n));
111: PetscCall(PetscLayoutSetUp(map));
112: PetscCall(PetscLayoutGetLocalSize(map, &Nl));
114: /* global indexes in layout */
115: for (i = 0; i < npos; i++) gidxs[i] = (ilocal ? idxs[ilocal[i]] : idxs[i]) - gbounds[0];
116: PetscCall(ISRestoreIndices(subset, &idxs));
117: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)subset), &sf));
118: PetscCall(PetscSFSetGraphLayout(sf, map, npos, ilocal, PETSC_USE_POINTER, gidxs));
119: PetscCall(PetscLayoutDestroy(&map));
121: /* reduce from leaves to roots */
122: PetscCall(PetscCalloc1(Nl, &root_data));
123: PetscCall(PetscSFReduceBegin(sf, MPIU_INT, leaf_data, root_data, MPI_MAX));
124: PetscCall(PetscSFReduceEnd(sf, MPIU_INT, leaf_data, root_data, MPI_MAX));
126: /* count indexes in local part of layout */
127: nlocals = 0;
128: first_index = -1;
129: first_found = PETSC_FALSE;
130: for (i = 0; i < Nl; i++) {
131: if (!first_found && root_data[i]) {
132: first_found = PETSC_TRUE;
133: first_index = i;
134: }
135: nlocals += root_data[i];
136: }
138: /* cumulative of number of indexes and size of subset without holes */
139: #if defined(PETSC_HAVE_MPI_EXSCAN)
140: start = 0;
141: PetscCallMPI(MPI_Exscan(&nlocals, &start, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)subset)));
142: #else
143: PetscCallMPI(MPI_Scan(&nlocals, &start, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)subset)));
144: start = start - nlocals;
145: #endif
147: if (N) { /* compute total size of new subset if requested */
148: *N = start + nlocals;
149: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)subset), &commsize));
150: PetscCallMPI(MPI_Bcast(N, 1, MPIU_INT, commsize - 1, PetscObjectComm((PetscObject)subset)));
151: }
153: if (!subset_n) {
154: PetscCall(PetscFree(gidxs));
155: PetscCall(PetscSFDestroy(&sf));
156: PetscCall(PetscFree(leaf_data));
157: PetscCall(PetscFree(root_data));
158: PetscCall(PetscFree(ilocal));
159: if (subset_mult) PetscCall(ISRestoreIndices(subset_mult, &idxs_mult));
160: PetscFunctionReturn(PETSC_SUCCESS);
161: }
163: /* adapt root data with cumulative */
164: if (first_found) {
165: PetscInt old_index;
167: root_data[first_index] += start;
168: old_index = first_index;
169: for (i = first_index + 1; i < Nl; i++) {
170: if (root_data[i]) {
171: root_data[i] += root_data[old_index];
172: old_index = i;
173: }
174: }
175: }
177: /* from roots to leaves */
178: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, root_data, leaf_data, MPI_REPLACE));
179: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, root_data, leaf_data, MPI_REPLACE));
180: PetscCall(PetscSFDestroy(&sf));
182: /* create new IS with global indexes without holes */
183: for (i = 0; i < n_n; i++) gidxs[i] = -1;
184: if (subset_mult) {
185: PetscInt cum;
187: isblock = PETSC_FALSE;
188: for (i = 0, cum = 0; i < n; i++)
189: for (PetscInt j = 0; j < idxs_mult[i]; j++) gidxs[cum++] = leaf_data[i] - idxs_mult[i] + j;
190: } else
191: for (i = 0; i < n; i++) gidxs[i] = leaf_data[i] - 1;
193: if (isblock) {
194: if (ibs > 1)
195: for (i = 0; i < n_n / ibs; i++) gidxs[i] = gidxs[i * ibs] / ibs;
196: PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)subset), ibs, n_n / ibs, gidxs, PETSC_COPY_VALUES, subset_n));
197: } else {
198: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)subset), n_n, gidxs, PETSC_COPY_VALUES, subset_n));
199: }
200: if (subset_mult) PetscCall(ISRestoreIndices(subset_mult, &idxs_mult));
201: PetscCall(PetscFree(gidxs));
202: PetscCall(PetscFree(leaf_data));
203: PetscCall(PetscFree(root_data));
204: PetscCall(PetscFree(ilocal));
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: /*@
209: ISCreateSubIS - Create a sub index set from a global index set selecting some components.
211: Collective
213: Input Parameters:
214: + is - the index set
215: - comps - which components we will extract from `is`
217: Output Parameters:
218: . subis - the new sub index set
220: Example usage:
221: We have an index set `is` living on 3 processes with the following values\:
222: | 4 9 0 | 2 6 7 | 10 11 1|
223: and another index set `comps` used to indicate which components of is we want to take,
224: | 7 5 | 1 2 | 0 4|
225: The output index set `subis` should look like\:
226: | 11 7 | 9 0 | 4 6|
228: Level: intermediate
230: .seealso: `IS`, `VecGetSubVector()`, `MatCreateSubMatrix()`
231: @*/
232: PetscErrorCode ISCreateSubIS(IS is, IS comps, IS *subis)
233: {
234: PetscSF sf;
235: const PetscInt *is_indices, *comps_indices;
236: PetscInt *subis_indices, nroots, nleaves, *mine, i, lidx;
237: PetscMPIInt owner;
238: PetscSFNode *remote;
239: MPI_Comm comm;
241: PetscFunctionBegin;
244: PetscAssertPointer(subis, 3);
246: PetscCall(PetscObjectGetComm((PetscObject)is, &comm));
247: PetscCall(ISGetLocalSize(comps, &nleaves));
248: PetscCall(ISGetLocalSize(is, &nroots));
249: PetscCall(PetscMalloc1(nleaves, &remote));
250: PetscCall(PetscMalloc1(nleaves, &mine));
251: PetscCall(ISGetIndices(comps, &comps_indices));
252: /*
253: * Construct a PetscSF in which "is" data serves as roots and "subis" is leaves.
254: * Root data are sent to leaves using PetscSFBcast().
255: * */
256: for (i = 0; i < nleaves; i++) {
257: mine[i] = i;
258: /* Connect a remote root with the current leaf. The value on the remote root
259: * will be received by the current local leaf.
260: * */
261: owner = -1;
262: lidx = -1;
263: PetscCall(PetscLayoutFindOwnerIndex(is->map, comps_indices[i], &owner, &lidx));
264: remote[i].rank = owner;
265: remote[i].index = lidx;
266: }
267: PetscCall(ISRestoreIndices(comps, &comps_indices));
268: PetscCall(PetscSFCreate(comm, &sf));
269: PetscCall(PetscSFSetFromOptions(sf));
270: PetscCall(PetscSFSetGraph(sf, nroots, nleaves, mine, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER));
272: PetscCall(PetscMalloc1(nleaves, &subis_indices));
273: PetscCall(ISGetIndices(is, &is_indices));
274: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, is_indices, subis_indices, MPI_REPLACE));
275: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, is_indices, subis_indices, MPI_REPLACE));
276: PetscCall(ISRestoreIndices(is, &is_indices));
277: PetscCall(PetscSFDestroy(&sf));
278: PetscCall(ISCreateGeneral(comm, nleaves, subis_indices, PETSC_OWN_POINTER, subis));
279: PetscFunctionReturn(PETSC_SUCCESS);
280: }
282: /*@
283: ISClearInfoCache - clear the cache of computed index set properties
285: Not Collective
287: Input Parameters:
288: + is - the index set
289: - clear_permanent_local - whether to remove the permanent status of local properties
291: Level: developer
293: Note:
294: Because all processes must agree on the global permanent status of a property,
295: the permanent status can only be changed with `ISSetInfo()`, because this routine is not collective
297: .seealso: `IS`, `ISInfo`, `ISInfoType`, `ISSetInfo()`
298: @*/
299: PetscErrorCode ISClearInfoCache(IS is, PetscBool clear_permanent_local)
300: {
301: PetscInt i, j;
303: PetscFunctionBegin;
306: for (i = 0; i < IS_INFO_MAX; i++) {
307: if (clear_permanent_local) is->info_permanent[IS_LOCAL][i] = PETSC_FALSE;
308: for (j = 0; j < 2; j++) {
309: if (!is->info_permanent[j][i]) is->info[j][i] = IS_INFO_UNKNOWN;
310: }
311: }
312: PetscFunctionReturn(PETSC_SUCCESS);
313: }
315: static PetscErrorCode ISSetInfo_Internal(IS is, ISInfo info, ISInfoType type, ISInfoBool ipermanent, PetscBool flg)
316: {
317: ISInfoBool iflg = flg ? IS_INFO_TRUE : IS_INFO_FALSE;
318: PetscInt itype = (type == IS_LOCAL) ? 0 : 1;
319: PetscBool permanent_set = (ipermanent == IS_INFO_UNKNOWN) ? PETSC_FALSE : PETSC_TRUE;
320: PetscBool permanent = (ipermanent == IS_INFO_TRUE) ? PETSC_TRUE : PETSC_FALSE;
322: PetscFunctionBegin;
323: /* set this property */
324: is->info[itype][(int)info] = iflg;
325: if (permanent_set) is->info_permanent[itype][(int)info] = permanent;
326: /* set implications */
327: switch (info) {
328: case IS_SORTED:
329: if (flg && type == IS_GLOBAL) { /* an array that is globally sorted is also locally sorted */
330: is->info[IS_LOCAL][(int)info] = IS_INFO_TRUE;
331: /* global permanence implies local permanence */
332: if (permanent_set && permanent) is->info_permanent[IS_LOCAL][(int)info] = PETSC_TRUE;
333: }
334: if (!flg) { /* if an array is not sorted, it cannot be an interval or the identity */
335: is->info[itype][IS_INTERVAL] = IS_INFO_FALSE;
336: is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
337: if (permanent_set) {
338: is->info_permanent[itype][IS_INTERVAL] = permanent;
339: is->info_permanent[itype][IS_IDENTITY] = permanent;
340: }
341: }
342: break;
343: case IS_UNIQUE:
344: if (flg && type == IS_GLOBAL) { /* an array that is globally unique is also locally unique */
345: is->info[IS_LOCAL][(int)info] = IS_INFO_TRUE;
346: /* global permanence implies local permanence */
347: if (permanent_set && permanent) is->info_permanent[IS_LOCAL][(int)info] = PETSC_TRUE;
348: }
349: if (!flg) { /* if an array is not unique, it cannot be a permutation, and interval, or the identity */
350: is->info[itype][IS_PERMUTATION] = IS_INFO_FALSE;
351: is->info[itype][IS_INTERVAL] = IS_INFO_FALSE;
352: is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
353: if (permanent_set) {
354: is->info_permanent[itype][IS_PERMUTATION] = permanent;
355: is->info_permanent[itype][IS_INTERVAL] = permanent;
356: is->info_permanent[itype][IS_IDENTITY] = permanent;
357: }
358: }
359: break;
360: case IS_PERMUTATION:
361: if (flg) { /* an array that is a permutation is unique and is unique locally */
362: is->info[itype][IS_UNIQUE] = IS_INFO_TRUE;
363: is->info[IS_LOCAL][IS_UNIQUE] = IS_INFO_TRUE;
364: if (permanent_set && permanent) {
365: is->info_permanent[itype][IS_UNIQUE] = PETSC_TRUE;
366: is->info_permanent[IS_LOCAL][IS_UNIQUE] = PETSC_TRUE;
367: }
368: } else { /* an array that is not a permutation cannot be the identity */
369: is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
370: if (permanent_set) is->info_permanent[itype][IS_IDENTITY] = permanent;
371: }
372: break;
373: case IS_INTERVAL:
374: if (flg) { /* an array that is an interval is sorted and unique */
375: is->info[itype][IS_SORTED] = IS_INFO_TRUE;
376: is->info[IS_LOCAL][IS_SORTED] = IS_INFO_TRUE;
377: is->info[itype][IS_UNIQUE] = IS_INFO_TRUE;
378: is->info[IS_LOCAL][IS_UNIQUE] = IS_INFO_TRUE;
379: if (permanent_set && permanent) {
380: is->info_permanent[itype][IS_SORTED] = PETSC_TRUE;
381: is->info_permanent[IS_LOCAL][IS_SORTED] = PETSC_TRUE;
382: is->info_permanent[itype][IS_UNIQUE] = PETSC_TRUE;
383: is->info_permanent[IS_LOCAL][IS_UNIQUE] = PETSC_TRUE;
384: }
385: } else { /* an array that is not an interval cannot be the identity */
386: is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
387: if (permanent_set) is->info_permanent[itype][IS_IDENTITY] = permanent;
388: }
389: break;
390: case IS_IDENTITY:
391: if (flg) { /* an array that is the identity is sorted, unique, an interval, and a permutation */
392: is->info[itype][IS_SORTED] = IS_INFO_TRUE;
393: is->info[IS_LOCAL][IS_SORTED] = IS_INFO_TRUE;
394: is->info[itype][IS_UNIQUE] = IS_INFO_TRUE;
395: is->info[IS_LOCAL][IS_UNIQUE] = IS_INFO_TRUE;
396: is->info[itype][IS_PERMUTATION] = IS_INFO_TRUE;
397: is->info[itype][IS_INTERVAL] = IS_INFO_TRUE;
398: is->info[IS_LOCAL][IS_INTERVAL] = IS_INFO_TRUE;
399: if (permanent_set && permanent) {
400: is->info_permanent[itype][IS_SORTED] = PETSC_TRUE;
401: is->info_permanent[IS_LOCAL][IS_SORTED] = PETSC_TRUE;
402: is->info_permanent[itype][IS_UNIQUE] = PETSC_TRUE;
403: is->info_permanent[IS_LOCAL][IS_UNIQUE] = PETSC_TRUE;
404: is->info_permanent[itype][IS_PERMUTATION] = PETSC_TRUE;
405: is->info_permanent[itype][IS_INTERVAL] = PETSC_TRUE;
406: is->info_permanent[IS_LOCAL][IS_INTERVAL] = PETSC_TRUE;
407: }
408: }
409: break;
410: default:
411: PetscCheck(type != IS_LOCAL, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown IS property");
412: SETERRQ(PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Unknown IS property");
413: }
414: PetscFunctionReturn(PETSC_SUCCESS);
415: }
417: // PetscClangLinter pragma disable: -fdoc-section-header-unknown
418: /*@
419: ISSetInfo - Set known information about an index set.
421: Logically Collective if `ISInfoType` is `IS_GLOBAL`
423: Input Parameters:
424: + is - the index set
425: . info - describing a property of the index set, one of those listed below,
426: . type - `IS_LOCAL` if the information describes the local portion of the index set,
427: `IS_GLOBAL` if it describes the whole index set
428: . permanent - `PETSC_TRUE` if it is known that the property will persist through changes to the index set, `PETSC_FALSE` otherwise
429: If the user sets a property as permanently known, it will bypass computation of that property
430: - flg - set the described property as true (`PETSC_TRUE`) or false (`PETSC_FALSE`)
432: Values of `info` Describing `IS` Structure:
433: + `IS_SORTED` - the [local part of the] index set is sorted in ascending order
434: . `IS_UNIQUE` - each entry in the [local part of the] index set is unique
435: . `IS_PERMUTATION` - the [local part of the] index set is a permutation of the integers {0, 1, ..., N-1}, where N is the size of the [local part of the] index set
436: . `IS_INTERVAL` - the [local part of the] index set is equal to a contiguous range of integers {f, f + 1, ..., f + N-1}
437: - `IS_IDENTITY` - the [local part of the] index set is equal to the integers {0, 1, ..., N-1}
439: Level: advanced
441: Notes:
442: If type is `IS_GLOBAL`, all processes that share the index set must pass the same value in flg
444: It is possible to set a property with `ISSetInfo()` that contradicts what would be previously computed with `ISGetInfo()`
446: .seealso: `ISInfo`, `ISInfoType`, `IS`
447: @*/
448: PetscErrorCode ISSetInfo(IS is, ISInfo info, ISInfoType type, PetscBool permanent, PetscBool flg)
449: {
450: MPI_Comm comm, errcomm;
451: PetscMPIInt size;
453: PetscFunctionBegin;
456: comm = PetscObjectComm((PetscObject)is);
457: if (type == IS_GLOBAL) {
461: errcomm = comm;
462: } else {
463: errcomm = PETSC_COMM_SELF;
464: }
466: PetscCheck(((int)info) > IS_INFO_MIN && ((int)info) < IS_INFO_MAX, errcomm, PETSC_ERR_ARG_OUTOFRANGE, "Options %d is out of range", (int)info);
468: PetscCallMPI(MPI_Comm_size(comm, &size));
469: /* do not use global values if size == 1: it makes it easier to keep the implications straight */
470: if (size == 1) type = IS_LOCAL;
471: PetscCall(ISSetInfo_Internal(is, info, type, permanent ? IS_INFO_TRUE : IS_INFO_FALSE, flg));
472: PetscFunctionReturn(PETSC_SUCCESS);
473: }
475: static PetscErrorCode ISGetInfo_Sorted_Private(IS is, ISInfoType type, PetscBool *flg)
476: {
477: MPI_Comm comm;
478: PetscMPIInt size, rank;
480: PetscFunctionBegin;
481: comm = PetscObjectComm((PetscObject)is);
482: PetscCallMPI(MPI_Comm_size(comm, &size));
483: PetscCallMPI(MPI_Comm_size(comm, &rank));
484: if (type == IS_GLOBAL && is->ops->sortedglobal) {
485: PetscUseTypeMethod(is, sortedglobal, flg);
486: } else {
487: PetscBool sortedLocal = PETSC_FALSE;
489: /* determine if the array is locally sorted */
490: if (type == IS_GLOBAL && size > 1) {
491: /* call ISGetInfo so that a cached value will be used if possible */
492: PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, &sortedLocal));
493: } else if (is->ops->sortedlocal) {
494: PetscUseTypeMethod(is, sortedlocal, &sortedLocal);
495: } else {
496: /* default: get the local indices and directly check */
497: const PetscInt *idx;
498: PetscInt n;
500: PetscCall(ISGetIndices(is, &idx));
501: PetscCall(ISGetLocalSize(is, &n));
502: PetscCall(PetscSortedInt(n, idx, &sortedLocal));
503: PetscCall(ISRestoreIndices(is, &idx));
504: }
506: if (type == IS_LOCAL || size == 1) {
507: *flg = sortedLocal;
508: } else {
509: PetscCall(MPIU_Allreduce(&sortedLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm));
510: if (*flg) {
511: PetscInt n, min = PETSC_MAX_INT, max = PETSC_MIN_INT;
512: PetscInt maxprev;
514: PetscCall(ISGetLocalSize(is, &n));
515: if (n) PetscCall(ISGetMinMax(is, &min, &max));
516: maxprev = PETSC_MIN_INT;
517: PetscCallMPI(MPI_Exscan(&max, &maxprev, 1, MPIU_INT, MPI_MAX, comm));
518: if (rank && (maxprev > min)) sortedLocal = PETSC_FALSE;
519: PetscCall(MPIU_Allreduce(&sortedLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm));
520: }
521: }
522: }
523: PetscFunctionReturn(PETSC_SUCCESS);
524: }
526: static PetscErrorCode ISGetIndicesCopy_Private(IS is, PetscInt idx[]);
528: static PetscErrorCode ISGetInfo_Unique_Private(IS is, ISInfoType type, PetscBool *flg)
529: {
530: MPI_Comm comm;
531: PetscMPIInt size, rank;
532: PetscInt i;
534: PetscFunctionBegin;
535: comm = PetscObjectComm((PetscObject)is);
536: PetscCallMPI(MPI_Comm_size(comm, &size));
537: PetscCallMPI(MPI_Comm_size(comm, &rank));
538: if (type == IS_GLOBAL && is->ops->uniqueglobal) {
539: PetscUseTypeMethod(is, uniqueglobal, flg);
540: } else {
541: PetscBool uniqueLocal;
542: PetscInt n = -1;
543: PetscInt *idx = NULL;
545: /* determine if the array is locally unique */
546: if (type == IS_GLOBAL && size > 1) {
547: /* call ISGetInfo so that a cached value will be used if possible */
548: PetscCall(ISGetInfo(is, IS_UNIQUE, IS_LOCAL, PETSC_TRUE, &uniqueLocal));
549: } else if (is->ops->uniquelocal) {
550: PetscUseTypeMethod(is, uniquelocal, &uniqueLocal);
551: } else {
552: /* default: get the local indices and directly check */
553: uniqueLocal = PETSC_TRUE;
554: PetscCall(ISGetLocalSize(is, &n));
555: PetscCall(PetscMalloc1(n, &idx));
556: PetscCall(ISGetIndicesCopy_Private(is, idx));
557: PetscCall(PetscIntSortSemiOrdered(n, idx));
558: for (i = 1; i < n; i++)
559: if (idx[i] == idx[i - 1]) break;
560: if (i < n) uniqueLocal = PETSC_FALSE;
561: }
563: PetscCall(PetscFree(idx));
564: if (type == IS_LOCAL || size == 1) {
565: *flg = uniqueLocal;
566: } else {
567: PetscCall(MPIU_Allreduce(&uniqueLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm));
568: if (*flg) {
569: PetscInt min = PETSC_MAX_INT, max = PETSC_MIN_INT, maxprev;
571: if (!idx) {
572: PetscCall(ISGetLocalSize(is, &n));
573: PetscCall(PetscMalloc1(n, &idx));
574: PetscCall(ISGetIndicesCopy_Private(is, idx));
575: }
576: PetscCall(PetscParallelSortInt(is->map, is->map, idx, idx));
577: if (n) {
578: min = idx[0];
579: max = idx[n - 1];
580: }
581: for (i = 1; i < n; i++)
582: if (idx[i] == idx[i - 1]) break;
583: if (i < n) uniqueLocal = PETSC_FALSE;
584: maxprev = PETSC_MIN_INT;
585: PetscCallMPI(MPI_Exscan(&max, &maxprev, 1, MPIU_INT, MPI_MAX, comm));
586: if (rank && (maxprev == min)) uniqueLocal = PETSC_FALSE;
587: PetscCall(MPIU_Allreduce(&uniqueLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm));
588: }
589: }
590: PetscCall(PetscFree(idx));
591: }
592: PetscFunctionReturn(PETSC_SUCCESS);
593: }
595: static PetscErrorCode ISGetInfo_Permutation(IS is, ISInfoType type, PetscBool *flg)
596: {
597: MPI_Comm comm;
598: PetscMPIInt size, rank;
600: PetscFunctionBegin;
601: comm = PetscObjectComm((PetscObject)is);
602: PetscCallMPI(MPI_Comm_size(comm, &size));
603: PetscCallMPI(MPI_Comm_size(comm, &rank));
604: if (type == IS_GLOBAL && is->ops->permglobal) {
605: PetscUseTypeMethod(is, permglobal, flg);
606: } else if (type == IS_LOCAL && is->ops->permlocal) {
607: PetscUseTypeMethod(is, permlocal, flg);
608: } else {
609: PetscBool permLocal;
610: PetscInt n, i, rStart;
611: PetscInt *idx;
613: PetscCall(ISGetLocalSize(is, &n));
614: PetscCall(PetscMalloc1(n, &idx));
615: PetscCall(ISGetIndicesCopy_Private(is, idx));
616: if (type == IS_GLOBAL) {
617: PetscCall(PetscParallelSortInt(is->map, is->map, idx, idx));
618: PetscCall(PetscLayoutGetRange(is->map, &rStart, NULL));
619: } else {
620: PetscCall(PetscIntSortSemiOrdered(n, idx));
621: rStart = 0;
622: }
623: permLocal = PETSC_TRUE;
624: for (i = 0; i < n; i++) {
625: if (idx[i] != rStart + i) break;
626: }
627: if (i < n) permLocal = PETSC_FALSE;
628: if (type == IS_LOCAL || size == 1) {
629: *flg = permLocal;
630: } else {
631: PetscCall(MPIU_Allreduce(&permLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm));
632: }
633: PetscCall(PetscFree(idx));
634: }
635: PetscFunctionReturn(PETSC_SUCCESS);
636: }
638: static PetscErrorCode ISGetInfo_Interval(IS is, ISInfoType type, PetscBool *flg)
639: {
640: MPI_Comm comm;
641: PetscMPIInt size, rank;
642: PetscInt i;
644: PetscFunctionBegin;
645: comm = PetscObjectComm((PetscObject)is);
646: PetscCallMPI(MPI_Comm_size(comm, &size));
647: PetscCallMPI(MPI_Comm_size(comm, &rank));
648: if (type == IS_GLOBAL && is->ops->intervalglobal) {
649: PetscUseTypeMethod(is, intervalglobal, flg);
650: } else {
651: PetscBool intervalLocal;
653: /* determine if the array is locally an interval */
654: if (type == IS_GLOBAL && size > 1) {
655: /* call ISGetInfo so that a cached value will be used if possible */
656: PetscCall(ISGetInfo(is, IS_INTERVAL, IS_LOCAL, PETSC_TRUE, &intervalLocal));
657: } else if (is->ops->intervallocal) {
658: PetscUseTypeMethod(is, intervallocal, &intervalLocal);
659: } else {
660: PetscInt n;
661: const PetscInt *idx;
662: /* default: get the local indices and directly check */
663: intervalLocal = PETSC_TRUE;
664: PetscCall(ISGetLocalSize(is, &n));
665: PetscCall(ISGetIndices(is, &idx));
666: for (i = 1; i < n; i++)
667: if (idx[i] != idx[i - 1] + 1) break;
668: if (i < n) intervalLocal = PETSC_FALSE;
669: PetscCall(ISRestoreIndices(is, &idx));
670: }
672: if (type == IS_LOCAL || size == 1) {
673: *flg = intervalLocal;
674: } else {
675: PetscCall(MPIU_Allreduce(&intervalLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm));
676: if (*flg) {
677: PetscInt n, min = PETSC_MAX_INT, max = PETSC_MIN_INT;
678: PetscInt maxprev;
680: PetscCall(ISGetLocalSize(is, &n));
681: if (n) PetscCall(ISGetMinMax(is, &min, &max));
682: maxprev = PETSC_MIN_INT;
683: PetscCallMPI(MPI_Exscan(&max, &maxprev, 1, MPIU_INT, MPI_MAX, comm));
684: if (rank && n && (maxprev != min - 1)) intervalLocal = PETSC_FALSE;
685: PetscCall(MPIU_Allreduce(&intervalLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm));
686: }
687: }
688: }
689: PetscFunctionReturn(PETSC_SUCCESS);
690: }
692: static PetscErrorCode ISGetInfo_Identity(IS is, ISInfoType type, PetscBool *flg)
693: {
694: MPI_Comm comm;
695: PetscMPIInt size, rank;
697: PetscFunctionBegin;
698: comm = PetscObjectComm((PetscObject)is);
699: PetscCallMPI(MPI_Comm_size(comm, &size));
700: PetscCallMPI(MPI_Comm_size(comm, &rank));
701: if (type == IS_GLOBAL && is->ops->intervalglobal) {
702: PetscBool isinterval;
704: PetscUseTypeMethod(is, intervalglobal, &isinterval);
705: *flg = PETSC_FALSE;
706: if (isinterval) {
707: PetscInt min;
709: PetscCall(ISGetMinMax(is, &min, NULL));
710: PetscCallMPI(MPI_Bcast(&min, 1, MPIU_INT, 0, comm));
711: if (min == 0) *flg = PETSC_TRUE;
712: }
713: } else if (type == IS_LOCAL && is->ops->intervallocal) {
714: PetscBool isinterval;
716: PetscUseTypeMethod(is, intervallocal, &isinterval);
717: *flg = PETSC_FALSE;
718: if (isinterval) {
719: PetscInt min;
721: PetscCall(ISGetMinMax(is, &min, NULL));
722: if (min == 0) *flg = PETSC_TRUE;
723: }
724: } else {
725: PetscBool identLocal;
726: PetscInt n, i, rStart;
727: const PetscInt *idx;
729: PetscCall(ISGetLocalSize(is, &n));
730: PetscCall(ISGetIndices(is, &idx));
731: PetscCall(PetscLayoutGetRange(is->map, &rStart, NULL));
732: identLocal = PETSC_TRUE;
733: for (i = 0; i < n; i++) {
734: if (idx[i] != rStart + i) break;
735: }
736: if (i < n) identLocal = PETSC_FALSE;
737: if (type == IS_LOCAL || size == 1) {
738: *flg = identLocal;
739: } else {
740: PetscCall(MPIU_Allreduce(&identLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm));
741: }
742: PetscCall(ISRestoreIndices(is, &idx));
743: }
744: PetscFunctionReturn(PETSC_SUCCESS);
745: }
747: /*@
748: ISGetInfo - Determine whether an index set satisfies a given property
750: Collective or Logically Collective if the type is `IS_GLOBAL` (logically collective if the value of the property has been permanently set with `ISSetInfo()`)
752: Input Parameters:
753: + is - the index set
754: . info - describing a property of the index set, one of those listed in the documentation of `ISSetInfo()`
755: . compute - if `PETSC_FALSE`, the property will not be computed if it is not already known and the property will be assumed to be false
756: - type - whether the property is local (`IS_LOCAL`) or global (`IS_GLOBAL`)
758: Output Parameter:
759: . flg - whether the property is true (`PETSC_TRUE`) or false (`PETSC_FALSE`)
761: Level: advanced
763: Notes:
764: `ISGetInfo()` uses cached values when possible, which will be incorrect if `ISSetInfo()` has been called with incorrect information.
766: To clear cached values, use `ISClearInfoCache()`.
768: .seealso: `IS`, `ISInfo`, `ISInfoType`, `ISSetInfo()`, `ISClearInfoCache()`
769: @*/
770: PetscErrorCode ISGetInfo(IS is, ISInfo info, ISInfoType type, PetscBool compute, PetscBool *flg)
771: {
772: MPI_Comm comm, errcomm;
773: PetscMPIInt rank, size;
774: PetscInt itype;
775: PetscBool hasprop;
776: PetscBool infer;
778: PetscFunctionBegin;
781: comm = PetscObjectComm((PetscObject)is);
782: if (type == IS_GLOBAL) {
784: errcomm = comm;
785: } else {
786: errcomm = PETSC_COMM_SELF;
787: }
789: PetscCallMPI(MPI_Comm_size(comm, &size));
790: PetscCallMPI(MPI_Comm_rank(comm, &rank));
792: PetscCheck(((int)info) > IS_INFO_MIN && ((int)info) < IS_INFO_MAX, errcomm, PETSC_ERR_ARG_OUTOFRANGE, "Options %d is out of range", (int)info);
793: if (size == 1) type = IS_LOCAL;
794: itype = (type == IS_LOCAL) ? 0 : 1;
795: hasprop = PETSC_FALSE;
796: infer = PETSC_FALSE;
797: if (is->info_permanent[itype][(int)info]) {
798: hasprop = (is->info[itype][(int)info] == IS_INFO_TRUE) ? PETSC_TRUE : PETSC_FALSE;
799: infer = PETSC_TRUE;
800: } else if ((itype == IS_LOCAL) && (is->info[IS_LOCAL][info] != IS_INFO_UNKNOWN)) {
801: /* we can cache local properties as long as we clear them when the IS changes */
802: /* NOTE: we only cache local values because there is no ISAssemblyBegin()/ISAssemblyEnd(),
803: so we have no way of knowing when a cached value has been invalidated by changes on a different process */
804: hasprop = (is->info[itype][(int)info] == IS_INFO_TRUE) ? PETSC_TRUE : PETSC_FALSE;
805: infer = PETSC_TRUE;
806: } else if (compute) {
807: switch (info) {
808: case IS_SORTED:
809: PetscCall(ISGetInfo_Sorted_Private(is, type, &hasprop));
810: break;
811: case IS_UNIQUE:
812: PetscCall(ISGetInfo_Unique_Private(is, type, &hasprop));
813: break;
814: case IS_PERMUTATION:
815: PetscCall(ISGetInfo_Permutation(is, type, &hasprop));
816: break;
817: case IS_INTERVAL:
818: PetscCall(ISGetInfo_Interval(is, type, &hasprop));
819: break;
820: case IS_IDENTITY:
821: PetscCall(ISGetInfo_Identity(is, type, &hasprop));
822: break;
823: default:
824: SETERRQ(errcomm, PETSC_ERR_ARG_OUTOFRANGE, "Unknown IS property");
825: }
826: infer = PETSC_TRUE;
827: }
828: /* call ISSetInfo_Internal to keep all of the implications straight */
829: if (infer) PetscCall(ISSetInfo_Internal(is, info, type, IS_INFO_UNKNOWN, hasprop));
830: *flg = hasprop;
831: PetscFunctionReturn(PETSC_SUCCESS);
832: }
834: static PetscErrorCode ISCopyInfo_Private(IS source, IS dest)
835: {
836: PetscFunctionBegin;
837: PetscCall(PetscArraycpy(&dest->info[0], &source->info[0], 2));
838: PetscCall(PetscArraycpy(&dest->info_permanent[0], &source->info_permanent[0], 2));
839: PetscFunctionReturn(PETSC_SUCCESS);
840: }
842: /*@
843: ISIdentity - Determines whether index set is the identity mapping.
845: Collective
847: Input Parameter:
848: . is - the index set
850: Output Parameter:
851: . ident - `PETSC_TRUE` if an identity, else `PETSC_FALSE`
853: Level: intermediate
855: Note:
856: If `ISSetIdentity()` (or `ISSetInfo()` for a permanent property) has been called,
857: `ISIdentity()` will return its answer without communication between processes, but
858: otherwise the output ident will be computed from `ISGetInfo()`,
859: which may require synchronization on the communicator of `is`. To avoid this computation,
860: call `ISGetInfo()` directly with the compute flag set to `PETSC_FALSE`, and ident will be assumed false.
862: .seealso: `IS`, `ISSetIdentity()`, `ISGetInfo()`
863: @*/
864: PetscErrorCode ISIdentity(IS is, PetscBool *ident)
865: {
866: PetscFunctionBegin;
868: PetscAssertPointer(ident, 2);
869: PetscCall(ISGetInfo(is, IS_IDENTITY, IS_GLOBAL, PETSC_TRUE, ident));
870: PetscFunctionReturn(PETSC_SUCCESS);
871: }
873: /*@
874: ISSetIdentity - Informs the index set that it is an identity.
876: Logically Collective
878: Input Parameter:
879: . is - the index set
881: Level: intermediate
883: Notes:
884: `is` will be considered the identity permanently, even if indices have been changes (for example, with
885: `ISGeneralSetIndices()`). It's a good idea to only set this property if `is` will not change in the future.
887: To clear this property, use `ISClearInfoCache()`.
889: Developer Notes:
890: Some of these info routines have statements about values changing in the `IS`, this seems to contradict the fact that `IS` cannot be changed?
892: .seealso: `IS`, `ISIdentity()`, `ISSetInfo()`, `ISClearInfoCache()`
893: @*/
894: PetscErrorCode ISSetIdentity(IS is)
895: {
896: PetscFunctionBegin;
898: PetscCall(ISSetInfo(is, IS_IDENTITY, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE));
899: PetscFunctionReturn(PETSC_SUCCESS);
900: }
902: /*@
903: ISContiguousLocal - Locates an index set with contiguous range within a global range, if possible
905: Not Collective
907: Input Parameters:
908: + is - the index set
909: . gstart - global start
910: - gend - global end
912: Output Parameters:
913: + start - start of contiguous block, as an offset from `gstart`
914: - contig - `PETSC_TRUE` if the index set refers to contiguous entries on this process, else `PETSC_FALSE`
916: Level: developer
918: .seealso: `IS`, `ISGetLocalSize()`, `VecGetOwnershipRange()`
919: @*/
920: PetscErrorCode ISContiguousLocal(IS is, PetscInt gstart, PetscInt gend, PetscInt *start, PetscBool *contig)
921: {
922: PetscFunctionBegin;
924: PetscAssertPointer(start, 4);
925: PetscAssertPointer(contig, 5);
926: PetscCheck(gstart <= gend, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "gstart %" PetscInt_FMT " must be less than or equal to gend %" PetscInt_FMT, gstart, gend);
927: *start = -1;
928: *contig = PETSC_FALSE;
929: PetscTryTypeMethod(is, contiguous, gstart, gend, start, contig);
930: PetscFunctionReturn(PETSC_SUCCESS);
931: }
933: /*@
934: ISPermutation - `PETSC_TRUE` or `PETSC_FALSE` depending on whether the
935: index set has been declared to be a permutation.
937: Logically Collective
939: Input Parameter:
940: . is - the index set
942: Output Parameter:
943: . perm - `PETSC_TRUE` if a permutation, else `PETSC_FALSE`
945: Level: intermediate
947: Note:
948: If it is not already known that `is` is a permutation (if `ISSetPermutation()`
949: or `ISSetInfo()` has not been called), this routine will not attempt to compute
950: whether the index set is a permutation and will assume `perm` is `PETSC_FALSE`.
951: To compute the value when it is not already known, use `ISGetInfo()` with
952: the compute flag set to `PETSC_TRUE`.
954: Developer Notes:
955: Perhaps some of these routines should use the `PetscBool3` enum to return appropriate values
957: .seealso: `IS`, `ISSetPermutation()`, `ISGetInfo()`
958: @*/
959: PetscErrorCode ISPermutation(IS is, PetscBool *perm)
960: {
961: PetscFunctionBegin;
963: PetscAssertPointer(perm, 2);
964: PetscCall(ISGetInfo(is, IS_PERMUTATION, IS_GLOBAL, PETSC_FALSE, perm));
965: PetscFunctionReturn(PETSC_SUCCESS);
966: }
968: /*@
969: ISSetPermutation - Informs the index set that it is a permutation.
971: Logically Collective
973: Input Parameter:
974: . is - the index set
976: Level: intermediate
978: Notes:
979: `is` will be considered a permutation permanently, even if indices have been changes (for example, with
980: `ISGeneralSetIndices()`). It's a good idea to only set this property if `is` will not change in the future.
982: To clear this property, use `ISClearInfoCache()`.
984: The debug version of the libraries (./configure --with-debugging=1) checks if the
985: index set is actually a permutation. The optimized version just believes you.
987: .seealso: `IS`, `ISPermutation()`, `ISSetInfo()`, `ISClearInfoCache().`
988: @*/
989: PetscErrorCode ISSetPermutation(IS is)
990: {
991: PetscFunctionBegin;
993: if (PetscDefined(USE_DEBUG)) {
994: PetscMPIInt size;
996: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is), &size));
997: if (size == 1) {
998: PetscInt i, n, *idx;
999: const PetscInt *iidx;
1001: PetscCall(ISGetSize(is, &n));
1002: PetscCall(PetscMalloc1(n, &idx));
1003: PetscCall(ISGetIndices(is, &iidx));
1004: PetscCall(PetscArraycpy(idx, iidx, n));
1005: PetscCall(PetscIntSortSemiOrdered(n, idx));
1006: for (i = 0; i < n; i++) PetscCheck(idx[i] == i, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Index set is not a permutation");
1007: PetscCall(PetscFree(idx));
1008: PetscCall(ISRestoreIndices(is, &iidx));
1009: }
1010: }
1011: PetscCall(ISSetInfo(is, IS_PERMUTATION, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE));
1012: PetscFunctionReturn(PETSC_SUCCESS);
1013: }
1015: /*@C
1016: ISDestroy - Destroys an index set.
1018: Collective
1020: Input Parameter:
1021: . is - the index set
1023: Level: beginner
1025: .seealso: `IS`, `ISCreateGeneral()`, `ISCreateStride()`, `ISCreateBlocked()`
1026: @*/
1027: PetscErrorCode ISDestroy(IS *is)
1028: {
1029: PetscFunctionBegin;
1030: if (!*is) PetscFunctionReturn(PETSC_SUCCESS);
1032: if (--((PetscObject)*is)->refct > 0) {
1033: *is = NULL;
1034: PetscFunctionReturn(PETSC_SUCCESS);
1035: }
1036: if ((*is)->complement) {
1037: PetscInt refcnt;
1038: PetscCall(PetscObjectGetReference((PetscObject)((*is)->complement), &refcnt));
1039: PetscCheck(refcnt <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Nonlocal IS has not been restored");
1040: PetscCall(ISDestroy(&(*is)->complement));
1041: }
1042: PetscTryTypeMethod(*is, destroy);
1043: PetscCall(PetscLayoutDestroy(&(*is)->map));
1044: /* Destroy local representations of offproc data. */
1045: PetscCall(PetscFree((*is)->total));
1046: PetscCall(PetscFree((*is)->nonlocal));
1047: PetscCall(PetscHeaderDestroy(is));
1048: PetscFunctionReturn(PETSC_SUCCESS);
1049: }
1051: /*@
1052: ISInvertPermutation - Creates a new permutation that is the inverse of
1053: a given permutation.
1055: Collective
1057: Input Parameters:
1058: + is - the index set
1059: - nlocal - number of indices on this processor in result (ignored for 1 processor) or
1060: use `PETSC_DECIDE`
1062: Output Parameter:
1063: . isout - the inverse permutation
1065: Level: intermediate
1067: Note:
1068: For parallel index sets this does the complete parallel permutation, but the
1069: code is not efficient for huge index sets (10,000,000 indices).
1071: .seealso: `IS`, `ISGetInfo()`, `ISSetPermutation()`, `ISGetPermutation()`
1072: @*/
1073: PetscErrorCode ISInvertPermutation(IS is, PetscInt nlocal, IS *isout)
1074: {
1075: PetscBool isperm, isidentity, issame;
1077: PetscFunctionBegin;
1079: PetscAssertPointer(isout, 3);
1080: PetscCall(ISGetInfo(is, IS_PERMUTATION, IS_GLOBAL, PETSC_TRUE, &isperm));
1081: PetscCheck(isperm, PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_WRONG, "Not a permutation");
1082: PetscCall(ISGetInfo(is, IS_IDENTITY, IS_GLOBAL, PETSC_TRUE, &isidentity));
1083: issame = PETSC_FALSE;
1084: if (isidentity) {
1085: PetscInt n;
1086: PetscBool isallsame;
1088: PetscCall(ISGetLocalSize(is, &n));
1089: issame = (PetscBool)(n == nlocal);
1090: PetscCall(MPIU_Allreduce(&issame, &isallsame, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1091: issame = isallsame;
1092: }
1093: if (issame) {
1094: PetscCall(ISDuplicate(is, isout));
1095: } else {
1096: PetscUseTypeMethod(is, invertpermutation, nlocal, isout);
1097: PetscCall(ISSetPermutation(*isout));
1098: }
1099: PetscFunctionReturn(PETSC_SUCCESS);
1100: }
1102: /*@
1103: ISGetSize - Returns the global length of an index set.
1105: Not Collective
1107: Input Parameter:
1108: . is - the index set
1110: Output Parameter:
1111: . size - the global size
1113: Level: beginner
1115: .seealso: `IS`
1116: @*/
1117: PetscErrorCode ISGetSize(IS is, PetscInt *size)
1118: {
1119: PetscFunctionBegin;
1121: PetscAssertPointer(size, 2);
1122: *size = is->map->N;
1123: PetscFunctionReturn(PETSC_SUCCESS);
1124: }
1126: /*@
1127: ISGetLocalSize - Returns the local (processor) length of an index set.
1129: Not Collective
1131: Input Parameter:
1132: . is - the index set
1134: Output Parameter:
1135: . size - the local size
1137: Level: beginner
1139: .seealso: `IS`, `ISGetSize()`
1140: @*/
1141: PetscErrorCode ISGetLocalSize(IS is, PetscInt *size)
1142: {
1143: PetscFunctionBegin;
1145: PetscAssertPointer(size, 2);
1146: *size = is->map->n;
1147: PetscFunctionReturn(PETSC_SUCCESS);
1148: }
1150: /*@
1151: ISGetLayout - get `PetscLayout` describing index set layout
1153: Not Collective
1155: Input Parameter:
1156: . is - the index set
1158: Output Parameter:
1159: . map - the layout
1161: Level: developer
1163: .seealso: `IS`, `PetscLayout`, `ISSetLayout()`, `ISGetSize()`, `ISGetLocalSize()`
1164: @*/
1165: PetscErrorCode ISGetLayout(IS is, PetscLayout *map)
1166: {
1167: PetscFunctionBegin;
1169: PetscAssertPointer(map, 2);
1170: *map = is->map;
1171: PetscFunctionReturn(PETSC_SUCCESS);
1172: }
1174: /*@
1175: ISSetLayout - set `PetscLayout` describing index set layout
1177: Collective
1179: Input Parameters:
1180: + is - the index set
1181: - map - the layout
1183: Level: developer
1185: Notes:
1186: Users should typically use higher level functions such as `ISCreateGeneral()`.
1188: This function can be useful in some special cases of constructing a new `IS`, e.g. after `ISCreate()` and before `ISLoad()`.
1189: Otherwise, it is only valid to replace the layout with a layout known to be equivalent.
1191: .seealso: `IS`, `PetscLayout`, `ISCreate()`, `ISGetLayout()`, `ISGetSize()`, `ISGetLocalSize()`
1192: @*/
1193: PetscErrorCode ISSetLayout(IS is, PetscLayout map)
1194: {
1195: PetscFunctionBegin;
1197: PetscAssertPointer(map, 2);
1198: PetscCall(PetscLayoutReference(map, &is->map));
1199: PetscFunctionReturn(PETSC_SUCCESS);
1200: }
1202: /*@C
1203: ISGetIndices - Returns a pointer to the indices. The user should call
1204: `ISRestoreIndices()` after having looked at the indices. The user should
1205: NOT change the indices.
1207: Not Collective
1209: Input Parameter:
1210: . is - the index set
1212: Output Parameter:
1213: . ptr - the location to put the pointer to the indices
1215: Level: intermediate
1217: Fortran Notes:
1218: `ISGetIndices()` Fortran binding is deprecated (since PETSc 3.19), use `ISGetIndicesF90()`
1220: .seealso: `IS`, `ISRestoreIndices()`, `ISGetIndicesF90()`
1221: @*/
1222: PetscErrorCode ISGetIndices(IS is, const PetscInt *ptr[])
1223: {
1224: PetscFunctionBegin;
1226: PetscAssertPointer(ptr, 2);
1227: PetscUseTypeMethod(is, getindices, ptr);
1228: PetscFunctionReturn(PETSC_SUCCESS);
1229: }
1231: /*@C
1232: ISGetMinMax - Gets the minimum and maximum values in an `IS`
1234: Not Collective
1236: Input Parameter:
1237: . is - the index set
1239: Output Parameters:
1240: + min - the minimum value
1241: - max - the maximum value
1243: Level: intermediate
1245: Notes:
1246: Empty index sets return min=`PETSC_MAX_INT` and max=`PETSC_MIN_INT`.
1248: In parallel, it returns the `min` and `max` of the local portion of `is`
1250: .seealso: `IS`, `ISGetIndices()`, `ISRestoreIndices()`, `ISGetIndicesF90()`
1251: @*/
1252: PetscErrorCode ISGetMinMax(IS is, PetscInt *min, PetscInt *max)
1253: {
1254: PetscFunctionBegin;
1256: if (min) *min = is->min;
1257: if (max) *max = is->max;
1258: PetscFunctionReturn(PETSC_SUCCESS);
1259: }
1261: /*@
1262: ISLocate - determine the location of an index within the local component of an index set
1264: Not Collective
1266: Input Parameters:
1267: + is - the index set
1268: - key - the search key
1270: Output Parameter:
1271: . location - if >= 0, a location within the index set that is equal to the key, otherwise the key is not in the index set
1273: Level: intermediate
1275: .seealso: `IS`
1276: @*/
1277: PetscErrorCode ISLocate(IS is, PetscInt key, PetscInt *location)
1278: {
1279: PetscFunctionBegin;
1280: if (is->ops->locate) {
1281: PetscUseTypeMethod(is, locate, key, location);
1282: } else {
1283: PetscInt numIdx;
1284: PetscBool sorted;
1285: const PetscInt *idx;
1287: PetscCall(ISGetLocalSize(is, &numIdx));
1288: PetscCall(ISGetIndices(is, &idx));
1289: PetscCall(ISSorted(is, &sorted));
1290: if (sorted) {
1291: PetscCall(PetscFindInt(key, numIdx, idx, location));
1292: } else {
1293: PetscInt i;
1295: *location = -1;
1296: for (i = 0; i < numIdx; i++) {
1297: if (idx[i] == key) {
1298: *location = i;
1299: break;
1300: }
1301: }
1302: }
1303: PetscCall(ISRestoreIndices(is, &idx));
1304: }
1305: PetscFunctionReturn(PETSC_SUCCESS);
1306: }
1308: /*@C
1309: ISRestoreIndices - Restores an index set to a usable state after a call to `ISGetIndices()`.
1311: Not Collective
1313: Input Parameters:
1314: + is - the index set
1315: - ptr - the pointer obtained by `ISGetIndices()`
1317: Level: intermediate
1319: Fortran Notes:
1320: `ISRestoreIndices()` Fortran binding is deprecated (since PETSc 3.19), use `ISRestoreIndicesF90()`
1322: .seealso: `IS`, `ISGetIndices()`, `ISRestoreIndicesF90()`
1323: @*/
1324: PetscErrorCode ISRestoreIndices(IS is, const PetscInt *ptr[])
1325: {
1326: PetscFunctionBegin;
1328: PetscAssertPointer(ptr, 2);
1329: PetscTryTypeMethod(is, restoreindices, ptr);
1330: PetscFunctionReturn(PETSC_SUCCESS);
1331: }
1333: static PetscErrorCode ISGatherTotal_Private(IS is)
1334: {
1335: PetscInt i, n, N;
1336: const PetscInt *lindices;
1337: MPI_Comm comm;
1338: PetscMPIInt rank, size, *sizes = NULL, *offsets = NULL, nn;
1340: PetscFunctionBegin;
1343: PetscCall(PetscObjectGetComm((PetscObject)is, &comm));
1344: PetscCallMPI(MPI_Comm_size(comm, &size));
1345: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1346: PetscCall(ISGetLocalSize(is, &n));
1347: PetscCall(PetscMalloc2(size, &sizes, size, &offsets));
1349: PetscCall(PetscMPIIntCast(n, &nn));
1350: PetscCallMPI(MPI_Allgather(&nn, 1, MPI_INT, sizes, 1, MPI_INT, comm));
1351: offsets[0] = 0;
1352: for (i = 1; i < size; ++i) offsets[i] = offsets[i - 1] + sizes[i - 1];
1353: N = offsets[size - 1] + sizes[size - 1];
1355: PetscCall(PetscMalloc1(N, &is->total));
1356: PetscCall(ISGetIndices(is, &lindices));
1357: PetscCallMPI(MPI_Allgatherv((void *)lindices, nn, MPIU_INT, is->total, sizes, offsets, MPIU_INT, comm));
1358: PetscCall(ISRestoreIndices(is, &lindices));
1359: is->local_offset = offsets[rank];
1360: PetscCall(PetscFree2(sizes, offsets));
1361: PetscFunctionReturn(PETSC_SUCCESS);
1362: }
1364: /*@C
1365: ISGetTotalIndices - Retrieve an array containing all indices across the communicator.
1367: Collective
1369: Input Parameter:
1370: . is - the index set
1372: Output Parameter:
1373: . indices - total indices with rank 0 indices first, and so on; total array size is
1374: the same as returned with `ISGetSize()`.
1376: Level: intermediate
1378: Notes:
1379: this is potentially nonscalable, but depends on the size of the total index set
1380: and the size of the communicator. This may be feasible for index sets defined on
1381: subcommunicators, such that the set size does not grow with `PETSC_WORLD_COMM`.
1382: Note also that there is no way to tell where the local part of the indices starts
1383: (use `ISGetIndices()` and `ISGetNonlocalIndices()` to retrieve just the local and just
1384: the nonlocal part (complement), respectively).
1386: .seealso: `IS`, `ISRestoreTotalIndices()`, `ISGetNonlocalIndices()`, `ISGetSize()`
1387: @*/
1388: PetscErrorCode ISGetTotalIndices(IS is, const PetscInt *indices[])
1389: {
1390: PetscMPIInt size;
1392: PetscFunctionBegin;
1394: PetscAssertPointer(indices, 2);
1395: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is), &size));
1396: if (size == 1) {
1397: PetscUseTypeMethod(is, getindices, indices);
1398: } else {
1399: if (!is->total) PetscCall(ISGatherTotal_Private(is));
1400: *indices = is->total;
1401: }
1402: PetscFunctionReturn(PETSC_SUCCESS);
1403: }
1405: /*@C
1406: ISRestoreTotalIndices - Restore the index array obtained with `ISGetTotalIndices()`.
1408: Not Collective.
1410: Input Parameters:
1411: + is - the index set
1412: - indices - index array; must be the array obtained with `ISGetTotalIndices()`
1414: Level: intermediate
1416: .seealso: `IS`, `ISGetNonlocalIndices()`
1417: @*/
1418: PetscErrorCode ISRestoreTotalIndices(IS is, const PetscInt *indices[])
1419: {
1420: PetscMPIInt size;
1422: PetscFunctionBegin;
1424: PetscAssertPointer(indices, 2);
1425: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is), &size));
1426: if (size == 1) {
1427: PetscUseTypeMethod(is, restoreindices, indices);
1428: } else {
1429: PetscCheck(is->total == *indices, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Index array pointer being restored does not point to the array obtained from the IS.");
1430: }
1431: PetscFunctionReturn(PETSC_SUCCESS);
1432: }
1434: /*@C
1435: ISGetNonlocalIndices - Retrieve an array of indices from remote processors
1436: in this communicator.
1438: Collective
1440: Input Parameter:
1441: . is - the index set
1443: Output Parameter:
1444: . indices - indices with rank 0 indices first, and so on, omitting
1445: the current rank. Total number of indices is the difference
1446: total and local, obtained with `ISGetSize()` and `ISGetLocalSize()`,
1447: respectively.
1449: Level: intermediate
1451: Notes:
1452: Restore the indices using `ISRestoreNonlocalIndices()`.
1454: The same scalability considerations as those for `ISGetTotalIndices()` apply here.
1456: .seealso: `IS`, `ISGetTotalIndices()`, `ISRestoreNonlocalIndices()`, `ISGetSize()`, `ISGetLocalSize().`
1457: @*/
1458: PetscErrorCode ISGetNonlocalIndices(IS is, const PetscInt *indices[])
1459: {
1460: PetscMPIInt size;
1461: PetscInt n, N;
1463: PetscFunctionBegin;
1465: PetscAssertPointer(indices, 2);
1466: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is), &size));
1467: if (size == 1) *indices = NULL;
1468: else {
1469: if (!is->total) PetscCall(ISGatherTotal_Private(is));
1470: PetscCall(ISGetLocalSize(is, &n));
1471: PetscCall(ISGetSize(is, &N));
1472: PetscCall(PetscMalloc1(N - n, &is->nonlocal));
1473: PetscCall(PetscArraycpy(is->nonlocal, is->total, is->local_offset));
1474: PetscCall(PetscArraycpy(is->nonlocal + is->local_offset, is->total + is->local_offset + n, N - is->local_offset - n));
1475: *indices = is->nonlocal;
1476: }
1477: PetscFunctionReturn(PETSC_SUCCESS);
1478: }
1480: /*@C
1481: ISRestoreNonlocalIndices - Restore the index array obtained with `ISGetNonlocalIndices()`.
1483: Not Collective.
1485: Input Parameters:
1486: + is - the index set
1487: - indices - index array; must be the array obtained with `ISGetNonlocalIndices()`
1489: Level: intermediate
1491: .seealso: `IS`, `ISGetTotalIndices()`, `ISGetNonlocalIndices()`, `ISRestoreTotalIndices()`
1492: @*/
1493: PetscErrorCode ISRestoreNonlocalIndices(IS is, const PetscInt *indices[])
1494: {
1495: PetscFunctionBegin;
1497: PetscAssertPointer(indices, 2);
1498: PetscCheck(is->nonlocal == *indices, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Index array pointer being restored does not point to the array obtained from the IS.");
1499: PetscFunctionReturn(PETSC_SUCCESS);
1500: }
1502: /*@
1503: ISGetNonlocalIS - Gather all nonlocal indices for this `IS` and present
1504: them as another sequential index set.
1506: Collective
1508: Input Parameter:
1509: . is - the index set
1511: Output Parameter:
1512: . complement - sequential `IS` with indices identical to the result of
1513: `ISGetNonlocalIndices()`
1515: Level: intermediate
1517: Notes:
1518: Complement represents the result of `ISGetNonlocalIndices()` as an `IS`.
1519: Therefore scalability issues similar to `ISGetNonlocalIndices()` apply.
1521: The resulting `IS` must be restored using `ISRestoreNonlocalIS()`.
1523: .seealso: `IS`, `ISGetNonlocalIndices()`, `ISRestoreNonlocalIndices()`, `ISAllGather()`, `ISGetSize()`
1524: @*/
1525: PetscErrorCode ISGetNonlocalIS(IS is, IS *complement)
1526: {
1527: PetscFunctionBegin;
1529: PetscAssertPointer(complement, 2);
1530: /* Check if the complement exists already. */
1531: if (is->complement) {
1532: *complement = is->complement;
1533: PetscCall(PetscObjectReference((PetscObject)is->complement));
1534: } else {
1535: PetscInt N, n;
1536: const PetscInt *idx;
1537: PetscCall(ISGetSize(is, &N));
1538: PetscCall(ISGetLocalSize(is, &n));
1539: PetscCall(ISGetNonlocalIndices(is, &idx));
1540: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, N - n, idx, PETSC_USE_POINTER, &is->complement));
1541: PetscCall(PetscObjectReference((PetscObject)is->complement));
1542: *complement = is->complement;
1543: }
1544: PetscFunctionReturn(PETSC_SUCCESS);
1545: }
1547: /*@
1548: ISRestoreNonlocalIS - Restore the `IS` obtained with `ISGetNonlocalIS()`.
1550: Not collective.
1552: Input Parameters:
1553: + is - the index set
1554: - complement - index set of `is`'s nonlocal indices
1556: Level: intermediate
1558: .seealso: `IS`, `ISGetNonlocalIS()`, `ISGetNonlocalIndices()`, `ISRestoreNonlocalIndices()`
1559: @*/
1560: PetscErrorCode ISRestoreNonlocalIS(IS is, IS *complement)
1561: {
1562: PetscInt refcnt;
1564: PetscFunctionBegin;
1566: PetscAssertPointer(complement, 2);
1567: PetscCheck(*complement == is->complement, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Complement IS being restored was not obtained with ISGetNonlocalIS()");
1568: PetscCall(PetscObjectGetReference((PetscObject)is->complement, &refcnt));
1569: PetscCheck(refcnt > 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Duplicate call to ISRestoreNonlocalIS() detected");
1570: PetscCall(PetscObjectDereference((PetscObject)is->complement));
1571: PetscFunctionReturn(PETSC_SUCCESS);
1572: }
1574: /*@C
1575: ISViewFromOptions - View an `IS` based on options in the options database
1577: Collective
1579: Input Parameters:
1580: + A - the index set
1581: . obj - Optional object that provides the prefix for the options database
1582: - name - command line option
1584: Level: intermediate
1586: Note:
1587: See `PetscObjectViewFromOptions()` for possible `PetscViewer` and `PetscViewerFormat` values
1589: .seealso: `IS`, `ISView()`, `PetscObjectViewFromOptions()`, `ISCreate()`
1590: @*/
1591: PetscErrorCode ISViewFromOptions(IS A, PetscObject obj, const char name[])
1592: {
1593: PetscFunctionBegin;
1595: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
1596: PetscFunctionReturn(PETSC_SUCCESS);
1597: }
1599: /*@C
1600: ISView - Displays an index set.
1602: Collective
1604: Input Parameters:
1605: + is - the index set
1606: - viewer - viewer used to display the set, for example `PETSC_VIEWER_STDOUT_SELF`.
1608: Level: intermediate
1610: .seealso: `IS`, `PetscViewer`, `PetscViewerASCIIOpen()`, `ISViewFromOptions()`
1611: @*/
1612: PetscErrorCode ISView(IS is, PetscViewer viewer)
1613: {
1614: PetscFunctionBegin;
1616: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)is), &viewer));
1618: PetscCheckSameComm(is, 1, viewer, 2);
1620: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)is, viewer));
1621: PetscCall(PetscLogEventBegin(IS_View, is, viewer, 0, 0));
1622: PetscUseTypeMethod(is, view, viewer);
1623: PetscCall(PetscLogEventEnd(IS_View, is, viewer, 0, 0));
1624: PetscFunctionReturn(PETSC_SUCCESS);
1625: }
1627: /*@
1628: ISLoad - Loads an index set that has been stored in binary or HDF5 format with `ISView()`.
1630: Collective
1632: Input Parameters:
1633: + is - the newly loaded index set, this needs to have been created with `ISCreate()` or some related function before a call to `ISLoad()`.
1634: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()` or HDF5 file viewer, obtained from `PetscViewerHDF5Open()`
1636: Level: intermediate
1638: Notes:
1639: IF using HDF5, you must assign the IS the same name as was used in `is`
1640: that was stored in the file using `PetscObjectSetName()`. Otherwise you will
1641: get the error message: "Cannot H5DOpen2() with Vec name NAMEOFOBJECT"
1643: .seealso: `IS`, `PetscViewerBinaryOpen()`, `ISView()`, `MatLoad()`, `VecLoad()`
1644: @*/
1645: PetscErrorCode ISLoad(IS is, PetscViewer viewer)
1646: {
1647: PetscBool isbinary, ishdf5;
1649: PetscFunctionBegin;
1652: PetscCheckSameComm(is, 1, viewer, 2);
1653: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1654: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1655: PetscCheck(isbinary || ishdf5, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1656: if (!((PetscObject)is)->type_name) PetscCall(ISSetType(is, ISGENERAL));
1657: PetscCall(PetscLogEventBegin(IS_Load, is, viewer, 0, 0));
1658: PetscUseTypeMethod(is, load, viewer);
1659: PetscCall(PetscLogEventEnd(IS_Load, is, viewer, 0, 0));
1660: PetscFunctionReturn(PETSC_SUCCESS);
1661: }
1663: /*@
1664: ISSort - Sorts the indices of an index set.
1666: Collective
1668: Input Parameter:
1669: . is - the index set
1671: Level: intermediate
1673: .seealso: `IS`, `ISSortRemoveDups()`, `ISSorted()`
1674: @*/
1675: PetscErrorCode ISSort(IS is)
1676: {
1677: PetscFunctionBegin;
1679: PetscUseTypeMethod(is, sort);
1680: PetscCall(ISSetInfo(is, IS_SORTED, IS_LOCAL, is->info_permanent[IS_LOCAL][IS_SORTED], PETSC_TRUE));
1681: PetscFunctionReturn(PETSC_SUCCESS);
1682: }
1684: /*@
1685: ISSortRemoveDups - Sorts the indices of an index set, removing duplicates.
1687: Collective
1689: Input Parameter:
1690: . is - the index set
1692: Level: intermediate
1694: .seealso: `IS`, `ISSort()`, `ISSorted()`
1695: @*/
1696: PetscErrorCode ISSortRemoveDups(IS is)
1697: {
1698: PetscFunctionBegin;
1700: PetscCall(ISClearInfoCache(is, PETSC_FALSE));
1701: PetscUseTypeMethod(is, sortremovedups);
1702: PetscCall(ISSetInfo(is, IS_SORTED, IS_LOCAL, is->info_permanent[IS_LOCAL][IS_SORTED], PETSC_TRUE));
1703: PetscCall(ISSetInfo(is, IS_UNIQUE, IS_LOCAL, is->info_permanent[IS_LOCAL][IS_UNIQUE], PETSC_TRUE));
1704: PetscFunctionReturn(PETSC_SUCCESS);
1705: }
1707: /*@
1708: ISToGeneral - Converts an IS object of any type to `ISGENERAL` type
1710: Collective
1712: Input Parameter:
1713: . is - the index set
1715: Level: intermediate
1717: .seealso: `IS`, `ISSorted()`
1718: @*/
1719: PetscErrorCode ISToGeneral(IS is)
1720: {
1721: PetscFunctionBegin;
1723: PetscUseTypeMethod(is, togeneral);
1724: PetscFunctionReturn(PETSC_SUCCESS);
1725: }
1727: /*@
1728: ISSorted - Checks the indices to determine whether they have been sorted.
1730: Not Collective
1732: Input Parameter:
1733: . is - the index set
1735: Output Parameter:
1736: . flg - output flag, either `PETSC_TRUE` if the index set is sorted,
1737: or `PETSC_FALSE` otherwise.
1739: Level: intermediate
1741: Note:
1742: For parallel IS objects this only indicates if the local part of `is`
1743: is sorted. So some processors may return `PETSC_TRUE` while others may
1744: return `PETSC_FALSE`.
1746: .seealso: `ISSort()`, `ISSortRemoveDups()`
1747: @*/
1748: PetscErrorCode ISSorted(IS is, PetscBool *flg)
1749: {
1750: PetscFunctionBegin;
1752: PetscAssertPointer(flg, 2);
1753: PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, flg));
1754: PetscFunctionReturn(PETSC_SUCCESS);
1755: }
1757: /*@
1758: ISDuplicate - Creates a duplicate copy of an index set.
1760: Collective
1762: Input Parameter:
1763: . is - the index set
1765: Output Parameter:
1766: . newIS - the copy of the index set
1768: Level: beginner
1770: .seealso: `IS`, `ISCreateGeneral()`, `ISCopy()`
1771: @*/
1772: PetscErrorCode ISDuplicate(IS is, IS *newIS)
1773: {
1774: PetscFunctionBegin;
1776: PetscAssertPointer(newIS, 2);
1777: PetscUseTypeMethod(is, duplicate, newIS);
1778: PetscCall(ISCopyInfo_Private(is, *newIS));
1779: PetscFunctionReturn(PETSC_SUCCESS);
1780: }
1782: /*@
1783: ISCopy - Copies an index set.
1785: Collective
1787: Input Parameter:
1788: . is - the index set
1790: Output Parameter:
1791: . isy - the copy of the index set
1793: Level: beginner
1795: .seealso: `IS`, `ISDuplicate()`, `ISShift()`
1796: @*/
1797: PetscErrorCode ISCopy(IS is, IS isy)
1798: {
1799: PetscInt bs, bsy;
1801: PetscFunctionBegin;
1804: PetscCheckSameComm(is, 1, isy, 2);
1805: if (is == isy) PetscFunctionReturn(PETSC_SUCCESS);
1806: PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
1807: PetscCall(PetscLayoutGetBlockSize(isy->map, &bsy));
1808: PetscCheck(is->map->N == isy->map->N, PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_INCOMP, "Index sets have different global size %" PetscInt_FMT " != %" PetscInt_FMT, is->map->N, isy->map->N);
1809: PetscCheck(is->map->n == isy->map->n, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Index sets have different local size %" PetscInt_FMT " != %" PetscInt_FMT, is->map->n, isy->map->n);
1810: PetscCheck(bs == bsy, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Index sets have different block size %" PetscInt_FMT " != %" PetscInt_FMT, bs, bsy);
1811: PetscCall(ISCopyInfo_Private(is, isy));
1812: isy->max = is->max;
1813: isy->min = is->min;
1814: PetscUseTypeMethod(is, copy, isy);
1815: PetscFunctionReturn(PETSC_SUCCESS);
1816: }
1818: /*@
1819: ISShift - Shift all indices by given offset
1821: Collective
1823: Input Parameters:
1824: + is - the index set
1825: - offset - the offset
1827: Output Parameter:
1828: . isy - the shifted copy of the input index set
1830: Level: beginner
1832: Notes:
1833: The `offset` can be different across processes.
1835: `is` and `isy` can be the same.
1837: .seealso: `ISDuplicate()`, `ISCopy()`
1838: @*/
1839: PetscErrorCode ISShift(IS is, PetscInt offset, IS isy)
1840: {
1841: PetscFunctionBegin;
1844: PetscCheckSameComm(is, 1, isy, 3);
1845: if (!offset) {
1846: PetscCall(ISCopy(is, isy));
1847: PetscFunctionReturn(PETSC_SUCCESS);
1848: }
1849: PetscCheck(is->map->N == isy->map->N, PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_INCOMP, "Index sets have different global size %" PetscInt_FMT " != %" PetscInt_FMT, is->map->N, isy->map->N);
1850: PetscCheck(is->map->n == isy->map->n, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Index sets have different local size %" PetscInt_FMT " != %" PetscInt_FMT, is->map->n, isy->map->n);
1851: PetscCheck(is->map->bs == isy->map->bs, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Index sets have different block size %" PetscInt_FMT " != %" PetscInt_FMT, is->map->bs, isy->map->bs);
1852: PetscCall(ISCopyInfo_Private(is, isy));
1853: isy->max = is->max + offset;
1854: isy->min = is->min + offset;
1855: PetscUseMethod(is, "ISShift_C", (IS, PetscInt, IS), (is, offset, isy));
1856: PetscFunctionReturn(PETSC_SUCCESS);
1857: }
1859: /*@
1860: ISOnComm - Split a parallel `IS` on subcomms (usually self) or concatenate index sets on subcomms into a parallel index set
1862: Collective
1864: Input Parameters:
1865: + is - index set
1866: . comm - communicator for new index set
1867: - mode - copy semantics, `PETSC_USE_POINTER` for no-copy if possible, otherwise `PETSC_COPY_VALUES`
1869: Output Parameter:
1870: . newis - new `IS` on `comm`
1872: Level: advanced
1874: Notes:
1875: It is usually desirable to create a parallel `IS` and look at the local part when necessary.
1877: This function is useful if serial `IS`s must be created independently, or to view many
1878: logically independent serial `IS`s.
1880: The input `IS` must have the same type on every MPI process.
1882: .seealso: `IS`
1883: @*/
1884: PetscErrorCode ISOnComm(IS is, MPI_Comm comm, PetscCopyMode mode, IS *newis)
1885: {
1886: PetscMPIInt match;
1888: PetscFunctionBegin;
1890: PetscAssertPointer(newis, 4);
1891: PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)is), comm, &match));
1892: if (mode != PETSC_COPY_VALUES && (match == MPI_IDENT || match == MPI_CONGRUENT)) {
1893: PetscCall(PetscObjectReference((PetscObject)is));
1894: *newis = is;
1895: } else PetscUseTypeMethod(is, oncomm, comm, mode, newis);
1896: PetscFunctionReturn(PETSC_SUCCESS);
1897: }
1899: /*@
1900: ISSetBlockSize - informs an index set that it has a given block size
1902: Logicall Collective
1904: Input Parameters:
1905: + is - index set
1906: - bs - block size
1908: Level: intermediate
1910: Notes:
1911: This is much like the block size for `Vec`s. It indicates that one can think of the indices as
1912: being in a collection of equal size blocks. For `ISBLOCK` these collections of blocks are all contiguous
1913: within a block but this is not the case for other `IS`. For example, an `IS` with entries {0, 2, 3, 4, 6, 7} could
1914: have a block size of three set.
1916: `ISBlockGetIndices()` only works for `ISBLOCK`, not others.
1918: .seealso: `IS`, `ISGetBlockSize()`, `ISCreateBlock()`, `ISBlockGetIndices()`,
1919: @*/
1920: PetscErrorCode ISSetBlockSize(IS is, PetscInt bs)
1921: {
1922: PetscFunctionBegin;
1925: PetscCheck(bs >= 1, PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Block size %" PetscInt_FMT ", must be positive", bs);
1926: if (PetscDefined(USE_DEBUG)) {
1927: const PetscInt *indices;
1928: PetscInt length, i, j;
1929: PetscCall(ISGetIndices(is, &indices));
1930: if (indices) {
1931: PetscCall(ISGetLocalSize(is, &length));
1932: PetscCheck(length % bs == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local size %" PetscInt_FMT " not compatible with proposed block size %" PetscInt_FMT, length, bs);
1933: for (i = 1; i < length / bs; i += bs) {
1934: for (j = 1; j < bs - 1; j++) {
1935: PetscCheck(indices[i * bs + j] == indices[(i - 1) * bs + j] + indices[i * bs] - indices[(i - 1) * bs], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Proposed block size %" PetscInt_FMT " is incompatible with the indices", bs);
1936: }
1937: }
1938: }
1939: PetscCall(ISRestoreIndices(is, &indices));
1940: }
1941: PetscUseTypeMethod(is, setblocksize, bs);
1942: PetscFunctionReturn(PETSC_SUCCESS);
1943: }
1945: /*@
1946: ISGetBlockSize - Returns the number of elements in a block.
1948: Not Collective
1950: Input Parameter:
1951: . is - the index set
1953: Output Parameter:
1954: . size - the number of elements in a block
1956: Level: intermediate
1958: Note:
1959: See `ISSetBlockSize()`
1961: .seealso: `IS`, `ISBlockGetSize()`, `ISGetSize()`, `ISCreateBlock()`, `ISSetBlockSize()`
1962: @*/
1963: PetscErrorCode ISGetBlockSize(IS is, PetscInt *size)
1964: {
1965: PetscFunctionBegin;
1966: PetscCall(PetscLayoutGetBlockSize(is->map, size));
1967: PetscFunctionReturn(PETSC_SUCCESS);
1968: }
1970: static PetscErrorCode ISGetIndicesCopy_Private(IS is, PetscInt idx[])
1971: {
1972: PetscInt len, i;
1973: const PetscInt *ptr;
1975: PetscFunctionBegin;
1976: PetscCall(ISGetLocalSize(is, &len));
1977: PetscCall(ISGetIndices(is, &ptr));
1978: for (i = 0; i < len; i++) idx[i] = ptr[i];
1979: PetscCall(ISRestoreIndices(is, &ptr));
1980: PetscFunctionReturn(PETSC_SUCCESS);
1981: }
1983: /*MC
1984: ISGetIndicesF90 - Accesses the elements of an index set from Fortran.
1985: The users should call `ISRestoreIndicesF90()` after having looked at the
1986: indices. The user should NOT change the indices.
1988: Synopsis:
1989: ISGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
1991: Not Collective
1993: Input Parameter:
1994: . x - index set
1996: Output Parameters:
1997: + xx_v - the Fortran pointer to the array
1998: - ierr - error code
2000: Example of Usage:
2001: .vb
2002: PetscInt, pointer xx_v(:)
2003: ....
2004: call ISGetIndicesF90(x,xx_v,ierr)
2005: a = xx_v(3)
2006: call ISRestoreIndicesF90(x,xx_v,ierr)
2007: .ve
2009: Level: intermediate
2011: .seealso: `ISRestoreIndicesF90()`, `ISGetIndices()`, `ISRestoreIndices()`
2012: M*/
2014: /*MC
2015: ISRestoreIndicesF90 - Restores an index set to a usable state after
2016: a call to `ISGetIndicesF90()`.
2018: Synopsis:
2019: ISRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
2021: Not Collective
2023: Input Parameters:
2024: + x - index set
2025: - xx_v - the Fortran pointer to the array
2027: Output Parameter:
2028: . ierr - error code
2030: Example of Usage:
2031: .vb
2032: PetscInt, pointer xx_v(:)
2033: ....
2034: call ISGetIndicesF90(x,xx_v,ierr)
2035: a = xx_v(3)
2036: call ISRestoreIndicesF90(x,xx_v,ierr)
2037: .ve
2039: Level: intermediate
2041: .seealso: `ISGetIndicesF90()`, `ISGetIndices()`, `ISRestoreIndices()`
2042: M*/
2044: /*MC
2045: ISBlockGetIndicesF90 - Accesses the elements of an index set from Fortran.
2046: The users should call `ISBlockRestoreIndicesF90()` after having looked at the
2047: indices. The user should NOT change the indices.
2049: Synopsis:
2050: ISBlockGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
2052: Not Collective
2054: Input Parameter:
2055: . x - index set
2057: Output Parameters:
2058: + xx_v - the Fortran pointer to the array
2059: - ierr - error code
2060: Example of Usage:
2061: .vb
2062: PetscInt, pointer xx_v(:)
2063: ....
2064: call ISBlockGetIndicesF90(x,xx_v,ierr)
2065: a = xx_v(3)
2066: call ISBlockRestoreIndicesF90(x,xx_v,ierr)
2067: .ve
2069: Level: intermediate
2071: .seealso: `ISBlockRestoreIndicesF90()`, `ISGetIndices()`, `ISRestoreIndices()`,
2072: `ISRestoreIndices()`
2073: M*/
2075: /*MC
2076: ISBlockRestoreIndicesF90 - Restores an index set to a usable state after
2077: a call to `ISBlockGetIndicesF90()`.
2079: Synopsis:
2080: ISBlockRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
2082: Not Collective
2084: Input Parameters:
2085: + x - index set
2086: - xx_v - the Fortran pointer to the array
2088: Output Parameter:
2089: . ierr - error code
2091: Example of Usage:
2092: .vb
2093: PetscInt, pointer xx_v(:)
2094: ....
2095: call ISBlockGetIndicesF90(x,xx_v,ierr)
2096: a = xx_v(3)
2097: call ISBlockRestoreIndicesF90(x,xx_v,ierr)
2098: .ve
2100: Level: intermediate
2102: .seealso: `ISBlockGetIndicesF90()`, `ISGetIndices()`, `ISRestoreIndices()`, `ISRestoreIndicesF90()`
2103: M*/