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: Notes: All negative entries are mapped to -1. Indices with non positive multiplicities are skipped.
29: Level: intermediate
31: .seealso:
32: @*/
33: PetscErrorCode ISRenumber(IS subset, IS subset_mult, PetscInt *N, IS *subset_n)
34: {
35: PetscSF sf;
36: PetscLayout map;
37: const PetscInt *idxs, *idxs_mult = NULL;
38: PetscInt *leaf_data, *root_data, *gidxs, *ilocal, *ilocalneg;
39: PetscInt N_n, n, i, lbounds[2], gbounds[2], Nl, ibs;
40: PetscInt n_n, nlocals, start, first_index, npos, nneg;
41: PetscMPIInt commsize;
42: PetscBool first_found, isblock;
47: else if (!subset_n) return 0;
48: ISGetLocalSize(subset, &n);
49: if (subset_mult) {
50: ISGetLocalSize(subset_mult, &i);
52: }
53: /* create workspace layout for computing global indices of subset */
54: PetscMalloc1(n, &ilocal);
55: PetscMalloc1(n, &ilocalneg);
56: ISGetIndices(subset, &idxs);
57: ISGetBlockSize(subset, &ibs);
58: PetscObjectTypeCompare((PetscObject)subset, ISBLOCK, &isblock);
59: if (subset_mult) ISGetIndices(subset_mult, &idxs_mult);
60: lbounds[0] = PETSC_MAX_INT;
61: lbounds[1] = PETSC_MIN_INT;
62: for (i = 0, npos = 0, nneg = 0; i < n; i++) {
63: if (idxs[i] < 0) {
64: ilocalneg[nneg++] = i;
65: continue;
66: }
67: if (idxs[i] < lbounds[0]) lbounds[0] = idxs[i];
68: if (idxs[i] > lbounds[1]) lbounds[1] = idxs[i];
69: ilocal[npos++] = i;
70: }
71: if (npos == n) {
72: PetscFree(ilocal);
73: PetscFree(ilocalneg);
74: }
76: /* create sf : leaf_data == multiplicity of indexes, root data == global index in layout */
77: PetscMalloc1(n, &leaf_data);
78: for (i = 0; i < n; i++) leaf_data[i] = idxs_mult ? PetscMax(idxs_mult[i], 0) : 1;
80: /* local size of new subset */
81: n_n = 0;
82: for (i = 0; i < n; i++) n_n += leaf_data[i];
83: if (ilocalneg)
84: for (i = 0; i < nneg; i++) leaf_data[ilocalneg[i]] = 0;
85: PetscFree(ilocalneg);
86: PetscMalloc1(PetscMax(n_n, n), &gidxs); /* allocating extra space to reuse gidxs */
87: /* check for early termination (all negative) */
88: PetscGlobalMinMaxInt(PetscObjectComm((PetscObject)subset), lbounds, gbounds);
89: if (gbounds[1] < gbounds[0]) {
90: if (N) *N = 0;
91: if (subset_n) { /* all negative */
92: for (i = 0; i < n_n; i++) gidxs[i] = -1;
93: ISCreateGeneral(PetscObjectComm((PetscObject)subset), n_n, gidxs, PETSC_COPY_VALUES, subset_n);
94: }
95: PetscFree(leaf_data);
96: PetscFree(gidxs);
97: ISRestoreIndices(subset, &idxs);
98: if (subset_mult) ISRestoreIndices(subset_mult, &idxs_mult);
99: PetscFree(ilocal);
100: PetscFree(ilocalneg);
101: return 0;
102: }
104: /* split work */
105: N_n = gbounds[1] - gbounds[0] + 1;
106: PetscLayoutCreate(PetscObjectComm((PetscObject)subset), &map);
107: PetscLayoutSetBlockSize(map, 1);
108: PetscLayoutSetSize(map, N_n);
109: PetscLayoutSetUp(map);
110: PetscLayoutGetLocalSize(map, &Nl);
112: /* global indexes in layout */
113: for (i = 0; i < npos; i++) gidxs[i] = (ilocal ? idxs[ilocal[i]] : idxs[i]) - gbounds[0];
114: ISRestoreIndices(subset, &idxs);
115: PetscSFCreate(PetscObjectComm((PetscObject)subset), &sf);
116: PetscSFSetGraphLayout(sf, map, npos, ilocal, PETSC_USE_POINTER, gidxs);
117: PetscLayoutDestroy(&map);
119: /* reduce from leaves to roots */
120: PetscCalloc1(Nl, &root_data);
121: PetscSFReduceBegin(sf, MPIU_INT, leaf_data, root_data, MPI_MAX);
122: PetscSFReduceEnd(sf, MPIU_INT, leaf_data, root_data, MPI_MAX);
124: /* count indexes in local part of layout */
125: nlocals = 0;
126: first_index = -1;
127: first_found = PETSC_FALSE;
128: for (i = 0; i < Nl; i++) {
129: if (!first_found && root_data[i]) {
130: first_found = PETSC_TRUE;
131: first_index = i;
132: }
133: nlocals += root_data[i];
134: }
136: /* cumulative of number of indexes and size of subset without holes */
137: #if defined(PETSC_HAVE_MPI_EXSCAN)
138: start = 0;
139: MPI_Exscan(&nlocals, &start, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)subset));
140: #else
141: MPI_Scan(&nlocals, &start, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)subset));
142: start = start - nlocals;
143: #endif
145: if (N) { /* compute total size of new subset if requested */
146: *N = start + nlocals;
147: MPI_Comm_size(PetscObjectComm((PetscObject)subset), &commsize);
148: MPI_Bcast(N, 1, MPIU_INT, commsize - 1, PetscObjectComm((PetscObject)subset));
149: }
151: if (!subset_n) {
152: PetscFree(gidxs);
153: PetscSFDestroy(&sf);
154: PetscFree(leaf_data);
155: PetscFree(root_data);
156: PetscFree(ilocal);
157: if (subset_mult) ISRestoreIndices(subset_mult, &idxs_mult);
158: return 0;
159: }
161: /* adapt root data with cumulative */
162: if (first_found) {
163: PetscInt old_index;
165: root_data[first_index] += start;
166: old_index = first_index;
167: for (i = first_index + 1; i < Nl; i++) {
168: if (root_data[i]) {
169: root_data[i] += root_data[old_index];
170: old_index = i;
171: }
172: }
173: }
175: /* from roots to leaves */
176: PetscSFBcastBegin(sf, MPIU_INT, root_data, leaf_data, MPI_REPLACE);
177: PetscSFBcastEnd(sf, MPIU_INT, root_data, leaf_data, MPI_REPLACE);
178: PetscSFDestroy(&sf);
180: /* create new IS with global indexes without holes */
181: for (i = 0; i < n_n; i++) gidxs[i] = -1;
182: if (subset_mult) {
183: PetscInt cum;
185: isblock = PETSC_FALSE;
186: for (i = 0, cum = 0; i < n; i++)
187: for (PetscInt j = 0; j < idxs_mult[i]; j++) gidxs[cum++] = leaf_data[i] - idxs_mult[i] + j;
188: } else
189: for (i = 0; i < n; i++) gidxs[i] = leaf_data[i] - 1;
191: if (isblock) {
192: if (ibs > 1)
193: for (i = 0; i < n_n / ibs; i++) gidxs[i] = gidxs[i * ibs] / ibs;
194: ISCreateBlock(PetscObjectComm((PetscObject)subset), ibs, n_n / ibs, gidxs, PETSC_COPY_VALUES, subset_n);
195: } else {
196: ISCreateGeneral(PetscObjectComm((PetscObject)subset), n_n, gidxs, PETSC_COPY_VALUES, subset_n);
197: }
198: if (subset_mult) ISRestoreIndices(subset_mult, &idxs_mult);
199: PetscFree(gidxs);
200: PetscFree(leaf_data);
201: PetscFree(root_data);
202: PetscFree(ilocal);
203: return 0;
204: }
206: /*@
207: ISCreateSubIS - Create a sub index set from a global index set selecting some components.
209: Collective
211: Input Parameters:
212: + is - the index set
213: - comps - which components we will extract from is
215: Output Parameters:
216: . subis - the new sub index set
218: Level: intermediate
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: .seealso: `VecGetSubVector()`, `MatCreateSubMatrix()`
229: @*/
230: PetscErrorCode ISCreateSubIS(IS is, IS comps, IS *subis)
231: {
232: PetscSF sf;
233: const PetscInt *is_indices, *comps_indices;
234: PetscInt *subis_indices, nroots, nleaves, *mine, i, lidx;
235: PetscMPIInt owner;
236: PetscSFNode *remote;
237: MPI_Comm comm;
243: PetscObjectGetComm((PetscObject)is, &comm);
244: ISGetLocalSize(comps, &nleaves);
245: ISGetLocalSize(is, &nroots);
246: PetscMalloc1(nleaves, &remote);
247: PetscMalloc1(nleaves, &mine);
248: ISGetIndices(comps, &comps_indices);
249: /*
250: * Construct a PetscSF in which "is" data serves as roots and "subis" is leaves.
251: * Root data are sent to leaves using PetscSFBcast().
252: * */
253: for (i = 0; i < nleaves; i++) {
254: mine[i] = i;
255: /* Connect a remote root with the current leaf. The value on the remote root
256: * will be received by the current local leaf.
257: * */
258: owner = -1;
259: lidx = -1;
260: PetscLayoutFindOwnerIndex(is->map, comps_indices[i], &owner, &lidx);
261: remote[i].rank = owner;
262: remote[i].index = lidx;
263: }
264: ISRestoreIndices(comps, &comps_indices);
265: PetscSFCreate(comm, &sf);
266: PetscSFSetFromOptions(sf);
267: PetscSFSetGraph(sf, nroots, nleaves, mine, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
269: PetscMalloc1(nleaves, &subis_indices);
270: ISGetIndices(is, &is_indices);
271: PetscSFBcastBegin(sf, MPIU_INT, is_indices, subis_indices, MPI_REPLACE);
272: PetscSFBcastEnd(sf, MPIU_INT, is_indices, subis_indices, MPI_REPLACE);
273: ISRestoreIndices(is, &is_indices);
274: PetscSFDestroy(&sf);
275: ISCreateGeneral(comm, nleaves, subis_indices, PETSC_OWN_POINTER, subis);
276: return 0;
277: }
279: /*@
280: ISClearInfoCache - clear the cache of computed index set properties
282: Not collective
284: Input Parameters:
285: + is - the index set
286: - clear_permanent_local - whether to remove the permanent status of local properties
288: NOTE: because all processes must agree on the global permanent status of a property,
289: the permanent status can only be changed with ISSetInfo(), because this routine is not collective
291: Level: developer
293: .seealso: `ISInfo`, `ISInfoType`, `ISSetInfo()`, `ISClearInfoCache()`
295: @*/
296: PetscErrorCode ISClearInfoCache(IS is, PetscBool clear_permanent_local)
297: {
298: PetscInt i, j;
302: for (i = 0; i < IS_INFO_MAX; i++) {
303: if (clear_permanent_local) is->info_permanent[IS_LOCAL][i] = PETSC_FALSE;
304: for (j = 0; j < 2; j++) {
305: if (!is->info_permanent[j][i]) is->info[j][i] = IS_INFO_UNKNOWN;
306: }
307: }
308: return 0;
309: }
311: static PetscErrorCode ISSetInfo_Internal(IS is, ISInfo info, ISInfoType type, ISInfoBool ipermanent, PetscBool flg)
312: {
313: ISInfoBool iflg = flg ? IS_INFO_TRUE : IS_INFO_FALSE;
314: PetscInt itype = (type == IS_LOCAL) ? 0 : 1;
315: PetscBool permanent_set = (ipermanent == IS_INFO_UNKNOWN) ? PETSC_FALSE : PETSC_TRUE;
316: PetscBool permanent = (ipermanent == IS_INFO_TRUE) ? PETSC_TRUE : PETSC_FALSE;
318: /* set this property */
319: is->info[itype][(int)info] = iflg;
320: if (permanent_set) is->info_permanent[itype][(int)info] = permanent;
321: /* set implications */
322: switch (info) {
323: case IS_SORTED:
324: if (flg && type == IS_GLOBAL) { /* an array that is globally sorted is also locally sorted */
325: is->info[IS_LOCAL][(int)info] = IS_INFO_TRUE;
326: /* global permanence implies local permanence */
327: if (permanent_set && permanent) is->info_permanent[IS_LOCAL][(int)info] = PETSC_TRUE;
328: }
329: if (!flg) { /* if an array is not sorted, it cannot be an interval or the identity */
330: is->info[itype][IS_INTERVAL] = IS_INFO_FALSE;
331: is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
332: if (permanent_set) {
333: is->info_permanent[itype][IS_INTERVAL] = permanent;
334: is->info_permanent[itype][IS_IDENTITY] = permanent;
335: }
336: }
337: break;
338: case IS_UNIQUE:
339: if (flg && type == IS_GLOBAL) { /* an array that is globally unique is also locally unique */
340: is->info[IS_LOCAL][(int)info] = IS_INFO_TRUE;
341: /* global permanence implies local permanence */
342: if (permanent_set && permanent) is->info_permanent[IS_LOCAL][(int)info] = PETSC_TRUE;
343: }
344: if (!flg) { /* if an array is not unique, it cannot be a permutation, and interval, or the identity */
345: is->info[itype][IS_PERMUTATION] = IS_INFO_FALSE;
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_PERMUTATION] = permanent;
350: is->info_permanent[itype][IS_INTERVAL] = permanent;
351: is->info_permanent[itype][IS_IDENTITY] = permanent;
352: }
353: }
354: break;
355: case IS_PERMUTATION:
356: if (flg) { /* an array that is a permutation is unique and is unique locally */
357: is->info[itype][IS_UNIQUE] = IS_INFO_TRUE;
358: is->info[IS_LOCAL][IS_UNIQUE] = IS_INFO_TRUE;
359: if (permanent_set && permanent) {
360: is->info_permanent[itype][IS_UNIQUE] = PETSC_TRUE;
361: is->info_permanent[IS_LOCAL][IS_UNIQUE] = PETSC_TRUE;
362: }
363: } else { /* an array that is not a permutation cannot be the identity */
364: is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
365: if (permanent_set) is->info_permanent[itype][IS_IDENTITY] = permanent;
366: }
367: break;
368: case IS_INTERVAL:
369: if (flg) { /* an array that is an interval is sorted and unique */
370: is->info[itype][IS_SORTED] = IS_INFO_TRUE;
371: is->info[IS_LOCAL][IS_SORTED] = IS_INFO_TRUE;
372: is->info[itype][IS_UNIQUE] = IS_INFO_TRUE;
373: is->info[IS_LOCAL][IS_UNIQUE] = IS_INFO_TRUE;
374: if (permanent_set && permanent) {
375: is->info_permanent[itype][IS_SORTED] = PETSC_TRUE;
376: is->info_permanent[IS_LOCAL][IS_SORTED] = PETSC_TRUE;
377: is->info_permanent[itype][IS_UNIQUE] = PETSC_TRUE;
378: is->info_permanent[IS_LOCAL][IS_UNIQUE] = PETSC_TRUE;
379: }
380: } else { /* an array that is not an interval cannot be the identity */
381: is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
382: if (permanent_set) is->info_permanent[itype][IS_IDENTITY] = permanent;
383: }
384: break;
385: case IS_IDENTITY:
386: if (flg) { /* an array that is the identity is sorted, unique, an interval, and a permutation */
387: is->info[itype][IS_SORTED] = IS_INFO_TRUE;
388: is->info[IS_LOCAL][IS_SORTED] = IS_INFO_TRUE;
389: is->info[itype][IS_UNIQUE] = IS_INFO_TRUE;
390: is->info[IS_LOCAL][IS_UNIQUE] = IS_INFO_TRUE;
391: is->info[itype][IS_PERMUTATION] = IS_INFO_TRUE;
392: is->info[itype][IS_INTERVAL] = IS_INFO_TRUE;
393: is->info[IS_LOCAL][IS_INTERVAL] = IS_INFO_TRUE;
394: if (permanent_set && permanent) {
395: is->info_permanent[itype][IS_SORTED] = PETSC_TRUE;
396: is->info_permanent[IS_LOCAL][IS_SORTED] = PETSC_TRUE;
397: is->info_permanent[itype][IS_UNIQUE] = PETSC_TRUE;
398: is->info_permanent[IS_LOCAL][IS_UNIQUE] = PETSC_TRUE;
399: is->info_permanent[itype][IS_PERMUTATION] = PETSC_TRUE;
400: is->info_permanent[itype][IS_INTERVAL] = PETSC_TRUE;
401: is->info_permanent[IS_LOCAL][IS_INTERVAL] = PETSC_TRUE;
402: }
403: }
404: break;
405: default:
407: SETERRQ(PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Unknown IS property");
408: }
409: return 0;
410: }
412: /*@
413: ISSetInfo - Set known information about an index set.
415: Logically Collective on IS if type is IS_GLOBAL
417: Input Parameters:
418: + is - the index set
419: . info - describing a property of the index set, one of those listed below,
420: . type - IS_LOCAL if the information describes the local portion of the index set,
421: IS_GLOBAL if it describes the whole index set
422: . permanent - PETSC_TRUE if it is known that the property will persist through changes to the index set, PETSC_FALSE otherwise
423: If the user sets a property as permanently known, it will bypass computation of that property
424: - flg - set the described property as true (PETSC_TRUE) or false (PETSC_FALSE)
426: Info Describing IS Structure:
427: + IS_SORTED - the [local part of the] index set is sorted in ascending order
428: . IS_UNIQUE - each entry in the [local part of the] index set is unique
429: . 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
430: . IS_INTERVAL - the [local part of the] index set is equal to a contiguous range of integers {f, f + 1, ..., f + N-1}
431: - IS_IDENTITY - the [local part of the] index set is equal to the integers {0, 1, ..., N-1}
433: Notes:
434: If type is IS_GLOBAL, all processes that share the index set must pass the same value in flg
436: It is possible to set a property with ISSetInfo() that contradicts what would be previously computed with ISGetInfo()
438: Level: advanced
440: .seealso: `ISInfo`, `ISInfoType`, `IS`
442: @*/
443: PetscErrorCode ISSetInfo(IS is, ISInfo info, ISInfoType type, PetscBool permanent, PetscBool flg)
444: {
445: MPI_Comm comm, errcomm;
446: PetscMPIInt size;
450: comm = PetscObjectComm((PetscObject)is);
451: if (type == IS_GLOBAL) {
455: errcomm = comm;
456: } else {
457: errcomm = PETSC_COMM_SELF;
458: }
462: MPI_Comm_size(comm, &size);
463: /* do not use global values if size == 1: it makes it easier to keep the implications straight */
464: if (size == 1) type = IS_LOCAL;
465: ISSetInfo_Internal(is, info, type, permanent ? IS_INFO_TRUE : IS_INFO_FALSE, flg);
466: return 0;
467: }
469: static PetscErrorCode ISGetInfo_Sorted(IS is, ISInfoType type, PetscBool *flg)
470: {
471: MPI_Comm comm;
472: PetscMPIInt size, rank;
474: comm = PetscObjectComm((PetscObject)is);
475: MPI_Comm_size(comm, &size);
476: MPI_Comm_size(comm, &rank);
477: if (type == IS_GLOBAL && is->ops->sortedglobal) {
478: PetscUseTypeMethod(is, sortedglobal, flg);
479: } else {
480: PetscBool sortedLocal = PETSC_FALSE;
482: /* determine if the array is locally sorted */
483: if (type == IS_GLOBAL && size > 1) {
484: /* call ISGetInfo so that a cached value will be used if possible */
485: ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, &sortedLocal);
486: } else if (is->ops->sortedlocal) {
487: PetscUseTypeMethod(is, sortedlocal, &sortedLocal);
488: } else {
489: /* default: get the local indices and directly check */
490: const PetscInt *idx;
491: PetscInt n;
493: ISGetIndices(is, &idx);
494: ISGetLocalSize(is, &n);
495: PetscSortedInt(n, idx, &sortedLocal);
496: ISRestoreIndices(is, &idx);
497: }
499: if (type == IS_LOCAL || size == 1) {
500: *flg = sortedLocal;
501: } else {
502: MPI_Allreduce(&sortedLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
503: if (*flg) {
504: PetscInt n, min = PETSC_MAX_INT, max = PETSC_MIN_INT;
505: PetscInt maxprev;
507: ISGetLocalSize(is, &n);
508: if (n) ISGetMinMax(is, &min, &max);
509: maxprev = PETSC_MIN_INT;
510: MPI_Exscan(&max, &maxprev, 1, MPIU_INT, MPI_MAX, comm);
511: if (rank && (maxprev > min)) sortedLocal = PETSC_FALSE;
512: MPI_Allreduce(&sortedLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
513: }
514: }
515: }
516: return 0;
517: }
519: PetscErrorCode ISGetIndicesCopy(IS is, PetscInt idx[]);
521: static PetscErrorCode ISGetInfo_Unique(IS is, ISInfoType type, PetscBool *flg)
522: {
523: MPI_Comm comm;
524: PetscMPIInt size, rank;
525: PetscInt i;
527: comm = PetscObjectComm((PetscObject)is);
528: MPI_Comm_size(comm, &size);
529: MPI_Comm_size(comm, &rank);
530: if (type == IS_GLOBAL && is->ops->uniqueglobal) {
531: PetscUseTypeMethod(is, uniqueglobal, flg);
532: } else {
533: PetscBool uniqueLocal;
534: PetscInt n = -1;
535: PetscInt *idx = NULL;
537: /* determine if the array is locally unique */
538: if (type == IS_GLOBAL && size > 1) {
539: /* call ISGetInfo so that a cached value will be used if possible */
540: ISGetInfo(is, IS_UNIQUE, IS_LOCAL, PETSC_TRUE, &uniqueLocal);
541: } else if (is->ops->uniquelocal) {
542: PetscUseTypeMethod(is, uniquelocal, &uniqueLocal);
543: } else {
544: /* default: get the local indices and directly check */
545: uniqueLocal = PETSC_TRUE;
546: ISGetLocalSize(is, &n);
547: PetscMalloc1(n, &idx);
548: ISGetIndicesCopy(is, idx);
549: PetscIntSortSemiOrdered(n, idx);
550: for (i = 1; i < n; i++)
551: if (idx[i] == idx[i - 1]) break;
552: if (i < n) uniqueLocal = PETSC_FALSE;
553: }
555: PetscFree(idx);
556: if (type == IS_LOCAL || size == 1) {
557: *flg = uniqueLocal;
558: } else {
559: MPI_Allreduce(&uniqueLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
560: if (*flg) {
561: PetscInt min = PETSC_MAX_INT, max = PETSC_MIN_INT, maxprev;
563: if (!idx) {
564: ISGetLocalSize(is, &n);
565: PetscMalloc1(n, &idx);
566: ISGetIndicesCopy(is, idx);
567: }
568: PetscParallelSortInt(is->map, is->map, idx, idx);
569: if (n) {
570: min = idx[0];
571: max = idx[n - 1];
572: }
573: for (i = 1; i < n; i++)
574: if (idx[i] == idx[i - 1]) break;
575: if (i < n) uniqueLocal = PETSC_FALSE;
576: maxprev = PETSC_MIN_INT;
577: MPI_Exscan(&max, &maxprev, 1, MPIU_INT, MPI_MAX, comm);
578: if (rank && (maxprev == min)) uniqueLocal = PETSC_FALSE;
579: MPI_Allreduce(&uniqueLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
580: }
581: }
582: PetscFree(idx);
583: }
584: return 0;
585: }
587: static PetscErrorCode ISGetInfo_Permutation(IS is, ISInfoType type, PetscBool *flg)
588: {
589: MPI_Comm comm;
590: PetscMPIInt size, rank;
592: comm = PetscObjectComm((PetscObject)is);
593: MPI_Comm_size(comm, &size);
594: MPI_Comm_size(comm, &rank);
595: if (type == IS_GLOBAL && is->ops->permglobal) {
596: PetscUseTypeMethod(is, permglobal, flg);
597: } else if (type == IS_LOCAL && is->ops->permlocal) {
598: PetscUseTypeMethod(is, permlocal, flg);
599: } else {
600: PetscBool permLocal;
601: PetscInt n, i, rStart;
602: PetscInt *idx;
604: ISGetLocalSize(is, &n);
605: PetscMalloc1(n, &idx);
606: ISGetIndicesCopy(is, idx);
607: if (type == IS_GLOBAL) {
608: PetscParallelSortInt(is->map, is->map, idx, idx);
609: PetscLayoutGetRange(is->map, &rStart, NULL);
610: } else {
611: PetscIntSortSemiOrdered(n, idx);
612: rStart = 0;
613: }
614: permLocal = PETSC_TRUE;
615: for (i = 0; i < n; i++) {
616: if (idx[i] != rStart + i) break;
617: }
618: if (i < n) permLocal = PETSC_FALSE;
619: if (type == IS_LOCAL || size == 1) {
620: *flg = permLocal;
621: } else {
622: MPI_Allreduce(&permLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
623: }
624: PetscFree(idx);
625: }
626: return 0;
627: }
629: static PetscErrorCode ISGetInfo_Interval(IS is, ISInfoType type, PetscBool *flg)
630: {
631: MPI_Comm comm;
632: PetscMPIInt size, rank;
633: PetscInt i;
635: comm = PetscObjectComm((PetscObject)is);
636: MPI_Comm_size(comm, &size);
637: MPI_Comm_size(comm, &rank);
638: if (type == IS_GLOBAL && is->ops->intervalglobal) {
639: PetscUseTypeMethod(is, intervalglobal, flg);
640: } else {
641: PetscBool intervalLocal;
643: /* determine if the array is locally an interval */
644: if (type == IS_GLOBAL && size > 1) {
645: /* call ISGetInfo so that a cached value will be used if possible */
646: ISGetInfo(is, IS_INTERVAL, IS_LOCAL, PETSC_TRUE, &intervalLocal);
647: } else if (is->ops->intervallocal) {
648: PetscUseTypeMethod(is, intervallocal, &intervalLocal);
649: } else {
650: PetscInt n;
651: const PetscInt *idx;
652: /* default: get the local indices and directly check */
653: intervalLocal = PETSC_TRUE;
654: ISGetLocalSize(is, &n);
655: ISGetIndices(is, &idx);
656: for (i = 1; i < n; i++)
657: if (idx[i] != idx[i - 1] + 1) break;
658: if (i < n) intervalLocal = PETSC_FALSE;
659: ISRestoreIndices(is, &idx);
660: }
662: if (type == IS_LOCAL || size == 1) {
663: *flg = intervalLocal;
664: } else {
665: MPI_Allreduce(&intervalLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
666: if (*flg) {
667: PetscInt n, min = PETSC_MAX_INT, max = PETSC_MIN_INT;
668: PetscInt maxprev;
670: ISGetLocalSize(is, &n);
671: if (n) ISGetMinMax(is, &min, &max);
672: maxprev = PETSC_MIN_INT;
673: MPI_Exscan(&max, &maxprev, 1, MPIU_INT, MPI_MAX, comm);
674: if (rank && n && (maxprev != min - 1)) intervalLocal = PETSC_FALSE;
675: MPI_Allreduce(&intervalLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
676: }
677: }
678: }
679: return 0;
680: }
682: static PetscErrorCode ISGetInfo_Identity(IS is, ISInfoType type, PetscBool *flg)
683: {
684: MPI_Comm comm;
685: PetscMPIInt size, rank;
687: comm = PetscObjectComm((PetscObject)is);
688: MPI_Comm_size(comm, &size);
689: MPI_Comm_size(comm, &rank);
690: if (type == IS_GLOBAL && is->ops->intervalglobal) {
691: PetscBool isinterval;
693: PetscUseTypeMethod(is, intervalglobal, &isinterval);
694: *flg = PETSC_FALSE;
695: if (isinterval) {
696: PetscInt min;
698: ISGetMinMax(is, &min, NULL);
699: MPI_Bcast(&min, 1, MPIU_INT, 0, comm);
700: if (min == 0) *flg = PETSC_TRUE;
701: }
702: } else if (type == IS_LOCAL && is->ops->intervallocal) {
703: PetscBool isinterval;
705: PetscUseTypeMethod(is, intervallocal, &isinterval);
706: *flg = PETSC_FALSE;
707: if (isinterval) {
708: PetscInt min;
710: ISGetMinMax(is, &min, NULL);
711: if (min == 0) *flg = PETSC_TRUE;
712: }
713: } else {
714: PetscBool identLocal;
715: PetscInt n, i, rStart;
716: const PetscInt *idx;
718: ISGetLocalSize(is, &n);
719: ISGetIndices(is, &idx);
720: PetscLayoutGetRange(is->map, &rStart, NULL);
721: identLocal = PETSC_TRUE;
722: for (i = 0; i < n; i++) {
723: if (idx[i] != rStart + i) break;
724: }
725: if (i < n) identLocal = PETSC_FALSE;
726: if (type == IS_LOCAL || size == 1) {
727: *flg = identLocal;
728: } else {
729: MPI_Allreduce(&identLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
730: }
731: ISRestoreIndices(is, &idx);
732: }
733: return 0;
734: }
736: /*@
737: ISGetInfo - Determine whether an index set satisfies a given property
739: Collective or logically collective on IS if the type is IS_GLOBAL (logically collective if the value of the property has been permanently set with ISSetInfo())
741: Input Parameters:
742: + is - the index set
743: . info - describing a property of the index set, one of those listed in the documentation of ISSetInfo()
744: . 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
745: - type - whether the property is local (IS_LOCAL) or global (IS_GLOBAL)
747: Output Parameter:
748: . flg - wheter the property is true (PETSC_TRUE) or false (PETSC_FALSE)
750: Note: ISGetInfo uses cached values when possible, which will be incorrect if ISSetInfo() has been called with incorrect information. To clear cached values, use ISClearInfoCache().
752: Level: advanced
754: .seealso: `ISInfo`, `ISInfoType`, `ISSetInfo()`, `ISClearInfoCache()`
756: @*/
757: PetscErrorCode ISGetInfo(IS is, ISInfo info, ISInfoType type, PetscBool compute, PetscBool *flg)
758: {
759: MPI_Comm comm, errcomm;
760: PetscMPIInt rank, size;
761: PetscInt itype;
762: PetscBool hasprop;
763: PetscBool infer;
767: comm = PetscObjectComm((PetscObject)is);
768: if (type == IS_GLOBAL) {
770: errcomm = comm;
771: } else {
772: errcomm = PETSC_COMM_SELF;
773: }
775: MPI_Comm_size(comm, &size);
776: MPI_Comm_rank(comm, &rank);
779: if (size == 1) type = IS_LOCAL;
780: itype = (type == IS_LOCAL) ? 0 : 1;
781: hasprop = PETSC_FALSE;
782: infer = PETSC_FALSE;
783: if (is->info_permanent[itype][(int)info]) {
784: hasprop = (is->info[itype][(int)info] == IS_INFO_TRUE) ? PETSC_TRUE : PETSC_FALSE;
785: infer = PETSC_TRUE;
786: } else if ((itype == IS_LOCAL) && (is->info[IS_LOCAL][info] != IS_INFO_UNKNOWN)) {
787: /* we can cache local properties as long as we clear them when the IS changes */
788: /* NOTE: we only cache local values because there is no ISAssemblyBegin()/ISAssemblyEnd(),
789: so we have no way of knowing when a cached value has been invalidated by changes on a different process */
790: hasprop = (is->info[itype][(int)info] == IS_INFO_TRUE) ? PETSC_TRUE : PETSC_FALSE;
791: infer = PETSC_TRUE;
792: } else if (compute) {
793: switch (info) {
794: case IS_SORTED:
795: ISGetInfo_Sorted(is, type, &hasprop);
796: break;
797: case IS_UNIQUE:
798: ISGetInfo_Unique(is, type, &hasprop);
799: break;
800: case IS_PERMUTATION:
801: ISGetInfo_Permutation(is, type, &hasprop);
802: break;
803: case IS_INTERVAL:
804: ISGetInfo_Interval(is, type, &hasprop);
805: break;
806: case IS_IDENTITY:
807: ISGetInfo_Identity(is, type, &hasprop);
808: break;
809: default:
810: SETERRQ(errcomm, PETSC_ERR_ARG_OUTOFRANGE, "Unknown IS property");
811: }
812: infer = PETSC_TRUE;
813: }
814: /* call ISSetInfo_Internal to keep all of the implications straight */
815: if (infer) ISSetInfo_Internal(is, info, type, IS_INFO_UNKNOWN, hasprop);
816: *flg = hasprop;
817: return 0;
818: }
820: static PetscErrorCode ISCopyInfo(IS source, IS dest)
821: {
822: PetscArraycpy(&dest->info[0], &source->info[0], 2);
823: PetscArraycpy(&dest->info_permanent[0], &source->info_permanent[0], 2);
824: return 0;
825: }
827: /*@
828: ISIdentity - Determines whether index set is the identity mapping.
830: Collective
832: Input Parameters:
833: . is - the index set
835: Output Parameters:
836: . ident - PETSC_TRUE if an identity, else PETSC_FALSE
838: Level: intermediate
840: Note: If ISSetIdentity() (or ISSetInfo() for a permanent property) has been called,
841: ISIdentity() will return its answer without communication between processes, but
842: otherwise the output ident will be computed from ISGetInfo(),
843: which may require synchronization on the communicator of IS. To avoid this computation,
844: call ISGetInfo() directly with the compute flag set to PETSC_FALSE, and ident will be assumed false.
846: .seealso: `ISSetIdentity()`, `ISGetInfo()`
847: @*/
848: PetscErrorCode ISIdentity(IS is, PetscBool *ident)
849: {
852: ISGetInfo(is, IS_IDENTITY, IS_GLOBAL, PETSC_TRUE, ident);
853: return 0;
854: }
856: /*@
857: ISSetIdentity - Informs the index set that it is an identity.
859: Logically Collective
861: Input Parameter:
862: . is - the index set
864: Level: intermediate
866: Note: The IS will be considered the identity permanently, even if indices have been changes (for example, with
867: ISGeneralSetIndices()). It's a good idea to only set this property if the IS will not change in the future.
868: To clear this property, use ISClearInfoCache().
870: .seealso: `ISIdentity()`, `ISSetInfo()`, `ISClearInfoCache()`
871: @*/
872: PetscErrorCode ISSetIdentity(IS is)
873: {
875: ISSetInfo(is, IS_IDENTITY, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE);
876: return 0;
877: }
879: /*@
880: ISContiguousLocal - Locates an index set with contiguous range within a global range, if possible
882: Not Collective
884: Input Parameters:
885: + is - the index set
886: . gstart - global start
887: - gend - global end
889: Output Parameters:
890: + start - start of contiguous block, as an offset from gstart
891: - contig - PETSC_TRUE if the index set refers to contiguous entries on this process, else PETSC_FALSE
893: Level: developer
895: .seealso: `ISGetLocalSize()`, `VecGetOwnershipRange()`
896: @*/
897: PetscErrorCode ISContiguousLocal(IS is, PetscInt gstart, PetscInt gend, PetscInt *start, PetscBool *contig)
898: {
902: *start = -1;
903: *contig = PETSC_FALSE;
904: PetscTryTypeMethod(is, contiguous, gstart, gend, start, contig);
905: return 0;
906: }
908: /*@
909: ISPermutation - PETSC_TRUE or PETSC_FALSE depending on whether the
910: index set has been declared to be a permutation.
912: Logically Collective
914: Input Parameter:
915: . is - the index set
917: Output Parameter:
918: . perm - PETSC_TRUE if a permutation, else PETSC_FALSE
920: Level: intermediate
922: Note: If it is not already known that the IS is a permutation (if ISSetPermutation()
923: or ISSetInfo() has not been called), this routine will not attempt to compute
924: whether the index set is a permutation and will assume perm is PETSC_FALSE.
925: To compute the value when it is not already known, use ISGetInfo() with
926: the compute flag set to PETSC_TRUE.
928: .seealso: `ISSetPermutation()`, `ISGetInfo()`
929: @*/
930: PetscErrorCode ISPermutation(IS is, PetscBool *perm)
931: {
934: ISGetInfo(is, IS_PERMUTATION, IS_GLOBAL, PETSC_FALSE, perm);
935: return 0;
936: }
938: /*@
939: ISSetPermutation - Informs the index set that it is a permutation.
941: Logically Collective
943: Input Parameter:
944: . is - the index set
946: Level: intermediate
948: The debug version of the libraries (./configure --with-debugging=1) checks if the
949: index set is actually a permutation. The optimized version just believes you.
951: Note: The IS will be considered a permutation permanently, even if indices have been changes (for example, with
952: ISGeneralSetIndices()). It's a good idea to only set this property if the IS will not change in the future.
953: To clear this property, use ISClearInfoCache().
955: .seealso: `ISPermutation()`, `ISSetInfo()`, `ISClearInfoCache().`
956: @*/
957: PetscErrorCode ISSetPermutation(IS is)
958: {
960: if (PetscDefined(USE_DEBUG)) {
961: PetscMPIInt size;
963: MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
964: if (size == 1) {
965: PetscInt i, n, *idx;
966: const PetscInt *iidx;
968: ISGetSize(is, &n);
969: PetscMalloc1(n, &idx);
970: ISGetIndices(is, &iidx);
971: PetscArraycpy(idx, iidx, n);
972: PetscIntSortSemiOrdered(n, idx);
974: PetscFree(idx);
975: ISRestoreIndices(is, &iidx);
976: }
977: }
978: ISSetInfo(is, IS_PERMUTATION, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE);
979: return 0;
980: }
982: /*@C
983: ISDestroy - Destroys an index set.
985: Collective
987: Input Parameters:
988: . is - the index set
990: Level: beginner
992: .seealso: `ISCreateGeneral()`, `ISCreateStride()`, `ISCreateBlocked()`
993: @*/
994: PetscErrorCode ISDestroy(IS *is)
995: {
996: if (!*is) return 0;
998: if (--((PetscObject)(*is))->refct > 0) {
999: *is = NULL;
1000: return 0;
1001: }
1002: if ((*is)->complement) {
1003: PetscInt refcnt;
1004: PetscObjectGetReference((PetscObject)((*is)->complement), &refcnt);
1006: ISDestroy(&(*is)->complement);
1007: }
1008: if ((*is)->ops->destroy) (*(*is)->ops->destroy)(*is);
1009: PetscLayoutDestroy(&(*is)->map);
1010: /* Destroy local representations of offproc data. */
1011: PetscFree((*is)->total);
1012: PetscFree((*is)->nonlocal);
1013: PetscHeaderDestroy(is);
1014: return 0;
1015: }
1017: /*@
1018: ISInvertPermutation - Creates a new permutation that is the inverse of
1019: a given permutation.
1021: Collective
1023: Input Parameters:
1024: + is - the index set
1025: - nlocal - number of indices on this processor in result (ignored for 1 proccessor) or
1026: use PETSC_DECIDE
1028: Output Parameter:
1029: . isout - the inverse permutation
1031: Level: intermediate
1033: Notes:
1034: For parallel index sets this does the complete parallel permutation, but the
1035: code is not efficient for huge index sets (10,000,000 indices).
1037: @*/
1038: PetscErrorCode ISInvertPermutation(IS is, PetscInt nlocal, IS *isout)
1039: {
1040: PetscBool isperm, isidentity, issame;
1044: ISGetInfo(is, IS_PERMUTATION, IS_GLOBAL, PETSC_TRUE, &isperm);
1046: ISGetInfo(is, IS_IDENTITY, IS_GLOBAL, PETSC_TRUE, &isidentity);
1047: issame = PETSC_FALSE;
1048: if (isidentity) {
1049: PetscInt n;
1050: PetscBool isallsame;
1052: ISGetLocalSize(is, &n);
1053: issame = (PetscBool)(n == nlocal);
1054: MPI_Allreduce(&issame, &isallsame, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is));
1055: issame = isallsame;
1056: }
1057: if (issame) {
1058: ISDuplicate(is, isout);
1059: } else {
1060: PetscUseTypeMethod(is, invertpermutation, nlocal, isout);
1061: ISSetPermutation(*isout);
1062: }
1063: return 0;
1064: }
1066: /*@
1067: ISGetSize - Returns the global length of an index set.
1069: Not Collective
1071: Input Parameter:
1072: . is - the index set
1074: Output Parameter:
1075: . size - the global size
1077: Level: beginner
1079: @*/
1080: PetscErrorCode ISGetSize(IS is, PetscInt *size)
1081: {
1084: *size = is->map->N;
1085: return 0;
1086: }
1088: /*@
1089: ISGetLocalSize - Returns the local (processor) length of an index set.
1091: Not Collective
1093: Input Parameter:
1094: . is - the index set
1096: Output Parameter:
1097: . size - the local size
1099: Level: beginner
1101: @*/
1102: PetscErrorCode ISGetLocalSize(IS is, PetscInt *size)
1103: {
1106: *size = is->map->n;
1107: return 0;
1108: }
1110: /*@
1111: ISGetLayout - get PetscLayout describing index set layout
1113: Not Collective
1115: Input Parameter:
1116: . is - the index set
1118: Output Parameter:
1119: . map - the layout
1121: Level: developer
1123: .seealso: `ISSetLayout()`, `ISGetSize()`, `ISGetLocalSize()`
1124: @*/
1125: PetscErrorCode ISGetLayout(IS is, PetscLayout *map)
1126: {
1129: *map = is->map;
1130: return 0;
1131: }
1133: /*@
1134: ISSetLayout - set PetscLayout describing index set layout
1136: Collective
1138: Input Arguments:
1139: + is - the index set
1140: - map - the layout
1142: Level: developer
1144: Notes:
1145: Users should typically use higher level functions such as ISCreateGeneral().
1147: This function can be useful in some special cases of constructing a new IS, e.g. after ISCreate() and before ISLoad().
1148: Otherwise, it is only valid to replace the layout with a layout known to be equivalent.
1150: .seealso: `ISCreate()`, `ISGetLayout()`, `ISGetSize()`, `ISGetLocalSize()`
1151: @*/
1152: PetscErrorCode ISSetLayout(IS is, PetscLayout map)
1153: {
1156: PetscLayoutReference(map, &is->map);
1157: return 0;
1158: }
1160: /*@C
1161: ISGetIndices - Returns a pointer to the indices. The user should call
1162: ISRestoreIndices() after having looked at the indices. The user should
1163: NOT change the indices.
1165: Not Collective
1167: Input Parameter:
1168: . is - the index set
1170: Output Parameter:
1171: . ptr - the location to put the pointer to the indices
1173: Fortran Note:
1174: This routine has two different interfaces from Fortran; the first is not recommend, it does not require Fortran 90
1175: $ IS is
1176: $ integer is_array(1)
1177: $ PetscOffset i_is
1178: $ int ierr
1179: $ call ISGetIndices(is,is_array,i_is,ierr)
1180: $
1181: $ Access first local entry in list
1182: $ value = is_array(i_is + 1)
1183: $
1184: $ ...... other code
1185: $ call ISRestoreIndices(is,is_array,i_is,ierr)
1186: The second Fortran interface is recommended.
1187: $ use petscisdef
1188: $ PetscInt, pointer :: array(:)
1189: $ PetscErrorCode ierr
1190: $ IS i
1191: $ call ISGetIndicesF90(i,array,ierr)
1193: See the Fortran chapter of the users manual and
1194: petsc/src/is/[tutorials,tests] for details.
1196: Level: intermediate
1198: .seealso: `ISRestoreIndices()`, `ISGetIndicesF90()`
1199: @*/
1200: PetscErrorCode ISGetIndices(IS is, const PetscInt *ptr[])
1201: {
1204: PetscUseTypeMethod(is, getindices, ptr);
1205: return 0;
1206: }
1208: /*@C
1209: ISGetMinMax - Gets the minimum and maximum values in an IS
1211: Not Collective
1213: Input Parameter:
1214: . is - the index set
1216: Output Parameters:
1217: + min - the minimum value
1218: - max - the maximum value
1220: Level: intermediate
1222: Notes:
1223: Empty index sets return min=PETSC_MAX_INT and max=PETSC_MIN_INT.
1224: In parallel, it returns the min and max of the local portion of the IS
1226: .seealso: `ISGetIndices()`, `ISRestoreIndices()`, `ISGetIndicesF90()`
1227: @*/
1228: PetscErrorCode ISGetMinMax(IS is, PetscInt *min, PetscInt *max)
1229: {
1231: if (min) *min = is->min;
1232: if (max) *max = is->max;
1233: return 0;
1234: }
1236: /*@
1237: ISLocate - determine the location of an index within the local component of an index set
1239: Not Collective
1241: Input Parameters:
1242: + is - the index set
1243: - key - the search key
1245: Output Parameter:
1246: . location - if >= 0, a location within the index set that is equal to the key, otherwise the key is not in the index set
1248: Level: intermediate
1249: @*/
1250: PetscErrorCode ISLocate(IS is, PetscInt key, PetscInt *location)
1251: {
1252: if (is->ops->locate) {
1253: PetscUseTypeMethod(is, locate, key, location);
1254: } else {
1255: PetscInt numIdx;
1256: PetscBool sorted;
1257: const PetscInt *idx;
1259: ISGetLocalSize(is, &numIdx);
1260: ISGetIndices(is, &idx);
1261: ISSorted(is, &sorted);
1262: if (sorted) {
1263: PetscFindInt(key, numIdx, idx, location);
1264: } else {
1265: PetscInt i;
1267: *location = -1;
1268: for (i = 0; i < numIdx; i++) {
1269: if (idx[i] == key) {
1270: *location = i;
1271: break;
1272: }
1273: }
1274: }
1275: ISRestoreIndices(is, &idx);
1276: }
1277: return 0;
1278: }
1280: /*@C
1281: ISRestoreIndices - Restores an index set to a usable state after a call
1282: to ISGetIndices().
1284: Not Collective
1286: Input Parameters:
1287: + is - the index set
1288: - ptr - the pointer obtained by ISGetIndices()
1290: Fortran Note:
1291: This routine is used differently from Fortran
1292: $ IS is
1293: $ integer is_array(1)
1294: $ PetscOffset i_is
1295: $ int ierr
1296: $ call ISGetIndices(is,is_array,i_is,ierr)
1297: $
1298: $ Access first local entry in list
1299: $ value = is_array(i_is + 1)
1300: $
1301: $ ...... other code
1302: $ call ISRestoreIndices(is,is_array,i_is,ierr)
1304: See the Fortran chapter of the users manual and
1305: petsc/src/vec/is/tests for details.
1307: Level: intermediate
1309: Note:
1310: This routine zeros out ptr. This is to prevent accidental us of the array after it has been restored.
1312: .seealso: `ISGetIndices()`, `ISRestoreIndicesF90()`
1313: @*/
1314: PetscErrorCode ISRestoreIndices(IS is, const PetscInt *ptr[])
1315: {
1318: PetscTryTypeMethod(is, restoreindices, ptr);
1319: return 0;
1320: }
1322: static PetscErrorCode ISGatherTotal_Private(IS is)
1323: {
1324: PetscInt i, n, N;
1325: const PetscInt *lindices;
1326: MPI_Comm comm;
1327: PetscMPIInt rank, size, *sizes = NULL, *offsets = NULL, nn;
1331: PetscObjectGetComm((PetscObject)is, &comm);
1332: MPI_Comm_size(comm, &size);
1333: MPI_Comm_rank(comm, &rank);
1334: ISGetLocalSize(is, &n);
1335: PetscMalloc2(size, &sizes, size, &offsets);
1337: PetscMPIIntCast(n, &nn);
1338: MPI_Allgather(&nn, 1, MPI_INT, sizes, 1, MPI_INT, comm);
1339: offsets[0] = 0;
1340: for (i = 1; i < size; ++i) offsets[i] = offsets[i - 1] + sizes[i - 1];
1341: N = offsets[size - 1] + sizes[size - 1];
1343: PetscMalloc1(N, &(is->total));
1344: ISGetIndices(is, &lindices);
1345: MPI_Allgatherv((void *)lindices, nn, MPIU_INT, is->total, sizes, offsets, MPIU_INT, comm);
1346: ISRestoreIndices(is, &lindices);
1347: is->local_offset = offsets[rank];
1348: PetscFree2(sizes, offsets);
1349: return 0;
1350: }
1352: /*@C
1353: ISGetTotalIndices - Retrieve an array containing all indices across the communicator.
1355: Collective
1357: Input Parameter:
1358: . is - the index set
1360: Output Parameter:
1361: . indices - total indices with rank 0 indices first, and so on; total array size is
1362: the same as returned with ISGetSize().
1364: Level: intermediate
1366: Notes:
1367: this is potentially nonscalable, but depends on the size of the total index set
1368: and the size of the communicator. This may be feasible for index sets defined on
1369: subcommunicators, such that the set size does not grow with PETSC_WORLD_COMM.
1370: Note also that there is no way to tell where the local part of the indices starts
1371: (use ISGetIndices() and ISGetNonlocalIndices() to retrieve just the local and just
1372: the nonlocal part (complement), respectively).
1374: .seealso: `ISRestoreTotalIndices()`, `ISGetNonlocalIndices()`, `ISGetSize()`
1375: @*/
1376: PetscErrorCode ISGetTotalIndices(IS is, const PetscInt *indices[])
1377: {
1378: PetscMPIInt size;
1382: MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
1383: if (size == 1) {
1384: PetscUseTypeMethod(is, getindices, indices);
1385: } else {
1386: if (!is->total) ISGatherTotal_Private(is);
1387: *indices = is->total;
1388: }
1389: return 0;
1390: }
1392: /*@C
1393: ISRestoreTotalIndices - Restore the index array obtained with ISGetTotalIndices().
1395: Not Collective.
1397: Input Parameters:
1398: + is - the index set
1399: - indices - index array; must be the array obtained with ISGetTotalIndices()
1401: Level: intermediate
1403: .seealso: `ISRestoreTotalIndices()`, `ISGetNonlocalIndices()`
1404: @*/
1405: PetscErrorCode ISRestoreTotalIndices(IS is, const PetscInt *indices[])
1406: {
1407: PetscMPIInt size;
1411: MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
1412: if (size == 1) {
1413: PetscUseTypeMethod(is, restoreindices, indices);
1414: } else {
1416: }
1417: return 0;
1418: }
1420: /*@C
1421: ISGetNonlocalIndices - Retrieve an array of indices from remote processors
1422: in this communicator.
1424: Collective
1426: Input Parameter:
1427: . is - the index set
1429: Output Parameter:
1430: . indices - indices with rank 0 indices first, and so on, omitting
1431: the current rank. Total number of indices is the difference
1432: total and local, obtained with ISGetSize() and ISGetLocalSize(),
1433: respectively.
1435: Level: intermediate
1437: Notes:
1438: restore the indices using ISRestoreNonlocalIndices().
1439: The same scalability considerations as those for ISGetTotalIndices
1440: apply here.
1442: .seealso: `ISGetTotalIndices()`, `ISRestoreNonlocalIndices()`, `ISGetSize()`, `ISGetLocalSize().`
1443: @*/
1444: PetscErrorCode ISGetNonlocalIndices(IS is, const PetscInt *indices[])
1445: {
1446: PetscMPIInt size;
1447: PetscInt n, N;
1451: MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
1452: if (size == 1) *indices = NULL;
1453: else {
1454: if (!is->total) ISGatherTotal_Private(is);
1455: ISGetLocalSize(is, &n);
1456: ISGetSize(is, &N);
1457: PetscMalloc1(N - n, &(is->nonlocal));
1458: PetscArraycpy(is->nonlocal, is->total, is->local_offset);
1459: PetscArraycpy(is->nonlocal + is->local_offset, is->total + is->local_offset + n, N - is->local_offset - n);
1460: *indices = is->nonlocal;
1461: }
1462: return 0;
1463: }
1465: /*@C
1466: ISRestoreNonlocalIndices - Restore the index array obtained with ISGetNonlocalIndices().
1468: Not Collective.
1470: Input Parameters:
1471: + is - the index set
1472: - indices - index array; must be the array obtained with ISGetNonlocalIndices()
1474: Level: intermediate
1476: .seealso: `ISGetTotalIndices()`, `ISGetNonlocalIndices()`, `ISRestoreTotalIndices()`
1477: @*/
1478: PetscErrorCode ISRestoreNonlocalIndices(IS is, const PetscInt *indices[])
1479: {
1483: return 0;
1484: }
1486: /*@
1487: ISGetNonlocalIS - Gather all nonlocal indices for this IS and present
1488: them as another sequential index set.
1490: Collective
1492: Input Parameter:
1493: . is - the index set
1495: Output Parameter:
1496: . complement - sequential IS with indices identical to the result of
1497: ISGetNonlocalIndices()
1499: Level: intermediate
1501: Notes:
1502: complement represents the result of ISGetNonlocalIndices as an IS.
1503: Therefore scalability issues similar to ISGetNonlocalIndices apply.
1504: The resulting IS must be restored using ISRestoreNonlocalIS().
1506: .seealso: `ISGetNonlocalIndices()`, `ISRestoreNonlocalIndices()`, `ISAllGather()`, `ISGetSize()`
1507: @*/
1508: PetscErrorCode ISGetNonlocalIS(IS is, IS *complement)
1509: {
1512: /* Check if the complement exists already. */
1513: if (is->complement) {
1514: *complement = is->complement;
1515: PetscObjectReference((PetscObject)(is->complement));
1516: } else {
1517: PetscInt N, n;
1518: const PetscInt *idx;
1519: ISGetSize(is, &N);
1520: ISGetLocalSize(is, &n);
1521: ISGetNonlocalIndices(is, &idx);
1522: ISCreateGeneral(PETSC_COMM_SELF, N - n, idx, PETSC_USE_POINTER, &(is->complement));
1523: PetscObjectReference((PetscObject)is->complement);
1524: *complement = is->complement;
1525: }
1526: return 0;
1527: }
1529: /*@
1530: ISRestoreNonlocalIS - Restore the IS obtained with ISGetNonlocalIS().
1532: Not collective.
1534: Input Parameters:
1535: + is - the index set
1536: - complement - index set of is's nonlocal indices
1538: Level: intermediate
1540: .seealso: `ISGetNonlocalIS()`, `ISGetNonlocalIndices()`, `ISRestoreNonlocalIndices()`
1541: @*/
1542: PetscErrorCode ISRestoreNonlocalIS(IS is, IS *complement)
1543: {
1544: PetscInt refcnt;
1549: PetscObjectGetReference((PetscObject)(is->complement), &refcnt);
1551: PetscObjectDereference((PetscObject)(is->complement));
1552: return 0;
1553: }
1555: /*@C
1556: ISViewFromOptions - View from Options
1558: Collective
1560: Input Parameters:
1561: + A - the index set
1562: . obj - Optional object
1563: - name - command line option
1565: Level: intermediate
1566: .seealso: `IS`, `ISView`, `PetscObjectViewFromOptions()`, `ISCreate()`
1567: @*/
1568: PetscErrorCode ISViewFromOptions(IS A, PetscObject obj, const char name[])
1569: {
1571: PetscObjectViewFromOptions((PetscObject)A, obj, name);
1572: return 0;
1573: }
1575: /*@C
1576: ISView - Displays an index set.
1578: Collective
1580: Input Parameters:
1581: + is - the index set
1582: - viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.
1584: Level: intermediate
1586: .seealso: `PetscViewerASCIIOpen()`
1587: @*/
1588: PetscErrorCode ISView(IS is, PetscViewer viewer)
1589: {
1591: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)is), &viewer);
1595: PetscObjectPrintClassNamePrefixType((PetscObject)is, viewer);
1596: PetscLogEventBegin(IS_View, is, viewer, 0, 0);
1597: PetscUseTypeMethod(is, view, viewer);
1598: PetscLogEventEnd(IS_View, is, viewer, 0, 0);
1599: return 0;
1600: }
1602: /*@
1603: ISLoad - Loads a vector that has been stored in binary or HDF5 format with ISView().
1605: Collective
1607: Input Parameters:
1608: + is - the newly loaded vector, this needs to have been created with ISCreate() or some related function before a call to ISLoad().
1609: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or HDF5 file viewer, obtained from PetscViewerHDF5Open()
1611: Level: intermediate
1613: Notes:
1614: IF using HDF5, you must assign the IS the same name as was used in the IS
1615: that was stored in the file using PetscObjectSetName(). Otherwise you will
1616: get the error message: "Cannot H5DOpen2() with Vec name NAMEOFOBJECT"
1618: .seealso: `PetscViewerBinaryOpen()`, `ISView()`, `MatLoad()`, `VecLoad()`
1619: @*/
1620: PetscErrorCode ISLoad(IS is, PetscViewer viewer)
1621: {
1622: PetscBool isbinary, ishdf5;
1627: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary);
1628: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
1630: if (!((PetscObject)is)->type_name) ISSetType(is, ISGENERAL);
1631: PetscLogEventBegin(IS_Load, is, viewer, 0, 0);
1632: PetscUseTypeMethod(is, load, viewer);
1633: PetscLogEventEnd(IS_Load, is, viewer, 0, 0);
1634: return 0;
1635: }
1637: /*@
1638: ISSort - Sorts the indices of an index set.
1640: Collective
1642: Input Parameters:
1643: . is - the index set
1645: Level: intermediate
1647: .seealso: `ISSortRemoveDups()`, `ISSorted()`
1648: @*/
1649: PetscErrorCode ISSort(IS is)
1650: {
1652: PetscUseTypeMethod(is, sort);
1653: ISSetInfo(is, IS_SORTED, IS_LOCAL, is->info_permanent[IS_LOCAL][IS_SORTED], PETSC_TRUE);
1654: return 0;
1655: }
1657: /*@
1658: ISSortRemoveDups - Sorts the indices of an index set, removing duplicates.
1660: Collective
1662: Input Parameters:
1663: . is - the index set
1665: Level: intermediate
1667: .seealso: `ISSort()`, `ISSorted()`
1668: @*/
1669: PetscErrorCode ISSortRemoveDups(IS is)
1670: {
1672: ISClearInfoCache(is, PETSC_FALSE);
1673: PetscUseTypeMethod(is, sortremovedups);
1674: ISSetInfo(is, IS_SORTED, IS_LOCAL, is->info_permanent[IS_LOCAL][IS_SORTED], PETSC_TRUE);
1675: ISSetInfo(is, IS_UNIQUE, IS_LOCAL, is->info_permanent[IS_LOCAL][IS_UNIQUE], PETSC_TRUE);
1676: return 0;
1677: }
1679: /*@
1680: ISToGeneral - Converts an IS object of any type to ISGENERAL type
1682: Collective
1684: Input Parameters:
1685: . is - the index set
1687: Level: intermediate
1689: .seealso: `ISSorted()`
1690: @*/
1691: PetscErrorCode ISToGeneral(IS is)
1692: {
1694: PetscUseTypeMethod(is, togeneral);
1695: return 0;
1696: }
1698: /*@
1699: ISSorted - Checks the indices to determine whether they have been sorted.
1701: Collective
1703: Input Parameter:
1704: . is - the index set
1706: Output Parameter:
1707: . flg - output flag, either PETSC_TRUE if the index set is sorted,
1708: or PETSC_FALSE otherwise.
1710: Notes:
1711: For parallel IS objects this only indicates if the local part of the IS
1712: is sorted. So some processors may return PETSC_TRUE while others may
1713: return PETSC_FALSE.
1715: Level: intermediate
1717: .seealso: `ISSort()`, `ISSortRemoveDups()`
1718: @*/
1719: PetscErrorCode ISSorted(IS is, PetscBool *flg)
1720: {
1723: ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, flg);
1724: return 0;
1725: }
1727: /*@
1728: ISDuplicate - Creates a duplicate copy of an index set.
1730: Collective
1732: Input Parameter:
1733: . is - the index set
1735: Output Parameter:
1736: . isnew - the copy of the index set
1738: Level: beginner
1740: .seealso: `ISCreateGeneral()`, `ISCopy()`
1741: @*/
1742: PetscErrorCode ISDuplicate(IS is, IS *newIS)
1743: {
1746: PetscUseTypeMethod(is, duplicate, newIS);
1747: ISCopyInfo(is, *newIS);
1748: return 0;
1749: }
1751: /*@
1752: ISCopy - Copies an index set.
1754: Collective
1756: Input Parameter:
1757: . is - the index set
1759: Output Parameter:
1760: . isy - the copy of the index set
1762: Level: beginner
1764: .seealso: `ISDuplicate()`, `ISShift()`
1765: @*/
1766: PetscErrorCode ISCopy(IS is, IS isy)
1767: {
1768: PetscInt bs, bsy;
1773: if (is == isy) return 0;
1774: PetscLayoutGetBlockSize(is->map, &bs);
1775: PetscLayoutGetBlockSize(isy->map, &bsy);
1779: ISCopyInfo(is, isy);
1780: isy->max = is->max;
1781: isy->min = is->min;
1782: PetscUseTypeMethod(is, copy, isy);
1783: return 0;
1784: }
1786: /*@
1787: ISShift - Shift all indices by given offset
1789: Collective
1791: Input Parameters:
1792: + is - the index set
1793: - offset - the offset
1795: Output Parameter:
1796: . isy - the shifted copy of the input index set
1798: Notes:
1799: The offset can be different across processes.
1800: IS is and isy can be the same.
1802: Level: beginner
1804: .seealso: `ISDuplicate()`, `ISCopy()`
1805: @*/
1806: PetscErrorCode ISShift(IS is, PetscInt offset, IS isy)
1807: {
1811: if (!offset) {
1812: ISCopy(is, isy);
1813: return 0;
1814: }
1818: ISCopyInfo(is, isy);
1819: isy->max = is->max + offset;
1820: isy->min = is->min + offset;
1821: PetscUseMethod(is, "ISShift_C", (IS, PetscInt, IS), (is, offset, isy));
1822: return 0;
1823: }
1825: /*@
1826: ISOnComm - Split a parallel IS on subcomms (usually self) or concatenate index sets on subcomms into a parallel index set
1828: Collective
1830: Input Parameters:
1831: + is - index set
1832: . comm - communicator for new index set
1833: - mode - copy semantics, PETSC_USE_POINTER for no-copy if possible, otherwise PETSC_COPY_VALUES
1835: Output Parameter:
1836: . newis - new IS on comm
1838: Level: advanced
1840: Notes:
1841: It is usually desirable to create a parallel IS and look at the local part when necessary.
1843: This function is useful if serial ISs must be created independently, or to view many
1844: logically independent serial ISs.
1846: The input IS must have the same type on every process.
1847: @*/
1848: PetscErrorCode ISOnComm(IS is, MPI_Comm comm, PetscCopyMode mode, IS *newis)
1849: {
1850: PetscMPIInt match;
1854: MPI_Comm_compare(PetscObjectComm((PetscObject)is), comm, &match);
1855: if (mode != PETSC_COPY_VALUES && (match == MPI_IDENT || match == MPI_CONGRUENT)) {
1856: PetscObjectReference((PetscObject)is);
1857: *newis = is;
1858: } else PetscUseTypeMethod(is, oncomm, comm, mode, newis);
1859: return 0;
1860: }
1862: /*@
1863: ISSetBlockSize - informs an index set that it has a given block size
1865: Logicall Collective
1867: Input Parameters:
1868: + is - index set
1869: - bs - block size
1871: Level: intermediate
1873: Notes:
1874: This is much like the block size for Vecs. It indicates that one can think of the indices as
1875: being in a collection of equal size blocks. For ISBlock() these collections of blocks are all contiquous
1876: within a block but this is not the case for other IS.
1877: ISBlockGetIndices() only works for ISBlock IS, not others.
1879: .seealso: `ISGetBlockSize()`, `ISCreateBlock()`, `ISBlockGetIndices()`,
1880: @*/
1881: PetscErrorCode ISSetBlockSize(IS is, PetscInt bs)
1882: {
1886: if (PetscDefined(USE_DEBUG)) {
1887: const PetscInt *indices;
1888: PetscInt length, i, j;
1889: ISGetIndices(is, &indices);
1890: if (indices) {
1891: ISGetLocalSize(is, &length);
1893: for (i = 0; i < length / bs; i += bs) {
1894: for (j = 0; j < bs - 1; j++) {
1896: }
1897: }
1898: }
1899: ISRestoreIndices(is, &indices);
1900: }
1901: PetscUseTypeMethod(is, setblocksize, bs);
1902: return 0;
1903: }
1905: /*@
1906: ISGetBlockSize - Returns the number of elements in a block.
1908: Not Collective
1910: Input Parameter:
1911: . is - the index set
1913: Output Parameter:
1914: . size - the number of elements in a block
1916: Level: intermediate
1918: Notes:
1919: This is much like the block size for Vecs. It indicates that one can think of the indices as
1920: being in a collection of equal size blocks. For ISBlock() these collections of blocks are all contiquous
1921: within a block but this is not the case for other IS.
1922: ISBlockGetIndices() only works for ISBlock IS, not others.
1924: .seealso: `ISBlockGetSize()`, `ISGetSize()`, `ISCreateBlock()`, `ISSetBlockSize()`
1925: @*/
1926: PetscErrorCode ISGetBlockSize(IS is, PetscInt *size)
1927: {
1928: PetscLayoutGetBlockSize(is->map, size);
1929: return 0;
1930: }
1932: PetscErrorCode ISGetIndicesCopy(IS is, PetscInt idx[])
1933: {
1934: PetscInt len, i;
1935: const PetscInt *ptr;
1937: ISGetLocalSize(is, &len);
1938: ISGetIndices(is, &ptr);
1939: for (i = 0; i < len; i++) idx[i] = ptr[i];
1940: ISRestoreIndices(is, &ptr);
1941: return 0;
1942: }
1944: /*MC
1945: ISGetIndicesF90 - Accesses the elements of an index set from Fortran90.
1946: The users should call ISRestoreIndicesF90() after having looked at the
1947: indices. The user should NOT change the indices.
1949: Synopsis:
1950: ISGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
1952: Not collective
1954: Input Parameter:
1955: . x - index set
1957: Output Parameters:
1958: + xx_v - the Fortran90 pointer to the array
1959: - ierr - error code
1961: Example of Usage:
1962: .vb
1963: PetscInt, pointer xx_v(:)
1964: ....
1965: call ISGetIndicesF90(x,xx_v,ierr)
1966: a = xx_v(3)
1967: call ISRestoreIndicesF90(x,xx_v,ierr)
1968: .ve
1970: Level: intermediate
1972: .seealso: `ISRestoreIndicesF90()`, `ISGetIndices()`, `ISRestoreIndices()`
1974: M*/
1976: /*MC
1977: ISRestoreIndicesF90 - Restores an index set to a usable state after
1978: a call to ISGetIndicesF90().
1980: Synopsis:
1981: ISRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
1983: Not collective
1985: Input Parameters:
1986: + x - index set
1987: - xx_v - the Fortran90 pointer to the array
1989: Output Parameter:
1990: . ierr - error code
1992: Example of Usage:
1993: .vb
1994: PetscInt, pointer xx_v(:)
1995: ....
1996: call ISGetIndicesF90(x,xx_v,ierr)
1997: a = xx_v(3)
1998: call ISRestoreIndicesF90(x,xx_v,ierr)
1999: .ve
2001: Level: intermediate
2003: .seealso: `ISGetIndicesF90()`, `ISGetIndices()`, `ISRestoreIndices()`
2005: M*/
2007: /*MC
2008: ISBlockGetIndicesF90 - Accesses the elements of an index set from Fortran90.
2009: The users should call ISBlockRestoreIndicesF90() after having looked at the
2010: indices. The user should NOT change the indices.
2012: Synopsis:
2013: ISBlockGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
2015: Not collective
2017: Input Parameter:
2018: . x - index set
2020: Output Parameters:
2021: + xx_v - the Fortran90 pointer to the array
2022: - ierr - error code
2023: Example of Usage:
2024: .vb
2025: PetscInt, pointer xx_v(:)
2026: ....
2027: call ISBlockGetIndicesF90(x,xx_v,ierr)
2028: a = xx_v(3)
2029: call ISBlockRestoreIndicesF90(x,xx_v,ierr)
2030: .ve
2032: Level: intermediate
2034: .seealso: `ISBlockRestoreIndicesF90()`, `ISGetIndices()`, `ISRestoreIndices()`,
2035: `ISRestoreIndices()`
2037: M*/
2039: /*MC
2040: ISBlockRestoreIndicesF90 - Restores an index set to a usable state after
2041: a call to ISBlockGetIndicesF90().
2043: Synopsis:
2044: ISBlockRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
2046: Not Collective
2048: Input Parameters:
2049: + x - index set
2050: - xx_v - the Fortran90 pointer to the array
2052: Output Parameter:
2053: . ierr - error code
2055: Example of Usage:
2056: .vb
2057: PetscInt, pointer xx_v(:)
2058: ....
2059: call ISBlockGetIndicesF90(x,xx_v,ierr)
2060: a = xx_v(3)
2061: call ISBlockRestoreIndicesF90(x,xx_v,ierr)
2062: .ve
2064: Notes:
2065: Not yet supported for all F90 compilers
2067: Level: intermediate
2069: .seealso: `ISBlockGetIndicesF90()`, `ISGetIndices()`, `ISRestoreIndices()`, `ISRestoreIndicesF90()`
2071: M*/