Actual source code: isdiff.c
2: #include <petsc/private/isimpl.h>
3: #include <petsc/private/sectionimpl.h>
4: #include <petscbt.h>
6: /*@
7: ISDifference - Computes the difference between two index sets.
9: Collective on is1
11: Input Parameters:
12: + is1 - first index, to have items removed from it
13: - is2 - index values to be removed
15: Output Parameters:
16: . isout - is1 - is2
18: Level: intermediate
20: Notes:
21: Negative values are removed from the lists. is2 may have values
22: that are not in is1. This requires O(imax-imin) memory and O(imax-imin)
23: work, where imin and imax are the bounds on the indices in is1.
25: If is2 is NULL, the result is the same as for an empty IS, i.e., a duplicate of is1.
27: The difference is computed separately on each MPI rank
29: .seealso: [](sec_scatter), `IS`, `ISDestroy()`, `ISView()`, `ISSum()`, `ISExpand()`
30: @*/
31: PetscErrorCode ISDifference(IS is1, IS is2, IS *isout)
32: {
33: PetscInt i, n1, n2, imin, imax, nout, *iout;
34: const PetscInt *i1, *i2;
35: PetscBT mask;
36: MPI_Comm comm;
40: if (!is2) {
41: ISDuplicate(is1, isout);
42: return 0;
43: }
46: ISGetIndices(is1, &i1);
47: ISGetLocalSize(is1, &n1);
49: /* Create a bit mask array to contain required values */
50: if (n1) {
51: imin = PETSC_MAX_INT;
52: imax = 0;
53: for (i = 0; i < n1; i++) {
54: if (i1[i] < 0) continue;
55: imin = PetscMin(imin, i1[i]);
56: imax = PetscMax(imax, i1[i]);
57: }
58: } else imin = imax = 0;
60: PetscBTCreate(imax - imin, &mask);
61: /* Put the values from is1 */
62: for (i = 0; i < n1; i++) {
63: if (i1[i] < 0) continue;
64: PetscBTSet(mask, i1[i] - imin);
65: }
66: ISRestoreIndices(is1, &i1);
67: /* Remove the values from is2 */
68: ISGetIndices(is2, &i2);
69: ISGetLocalSize(is2, &n2);
70: for (i = 0; i < n2; i++) {
71: if (i2[i] < imin || i2[i] > imax) continue;
72: PetscBTClear(mask, i2[i] - imin);
73: }
74: ISRestoreIndices(is2, &i2);
76: /* Count the number in the difference */
77: nout = 0;
78: for (i = 0; i < imax - imin + 1; i++) {
79: if (PetscBTLookup(mask, i)) nout++;
80: }
82: /* create the new IS containing the difference */
83: PetscMalloc1(nout, &iout);
84: nout = 0;
85: for (i = 0; i < imax - imin + 1; i++) {
86: if (PetscBTLookup(mask, i)) iout[nout++] = i + imin;
87: }
88: PetscObjectGetComm((PetscObject)is1, &comm);
89: ISCreateGeneral(comm, nout, iout, PETSC_OWN_POINTER, isout);
91: PetscBTDestroy(&mask);
92: return 0;
93: }
95: /*@
96: ISSum - Computes the sum (union) of two index sets.
98: Only sequential version (at the moment)
100: Input Parameters:
101: + is1 - index set to be extended
102: - is2 - index values to be added
104: Output Parameter:
105: . is3 - the sum; this can not be is1 or is2
107: Level: intermediate
109: Notes:
110: If n1 and n2 are the sizes of the sets, this takes O(n1+n2) time;
112: Both index sets need to be sorted on input.
114: The sum is computed separately on each MPI rank
116: .seealso: [](sec_scatter), `IS`, `ISDestroy()`, `ISView()`, `ISDifference()`, `ISExpand()`
117: @*/
118: PetscErrorCode ISSum(IS is1, IS is2, IS *is3)
119: {
120: MPI_Comm comm;
121: PetscBool f;
122: PetscMPIInt size;
123: const PetscInt *i1, *i2;
124: PetscInt n1, n2, n3, p1, p2, *iout;
128: PetscObjectGetComm((PetscObject)(is1), &comm);
129: MPI_Comm_size(comm, &size);
132: ISSorted(is1, &f);
134: ISSorted(is2, &f);
137: ISGetLocalSize(is1, &n1);
138: ISGetLocalSize(is2, &n2);
139: if (!n2) {
140: ISDuplicate(is1, is3);
141: return 0;
142: }
143: ISGetIndices(is1, &i1);
144: ISGetIndices(is2, &i2);
146: p1 = 0;
147: p2 = 0;
148: n3 = 0;
149: do {
150: if (p1 == n1) { /* cleanup for is2 */
151: n3 += n2 - p2;
152: break;
153: } else {
154: while (p2 < n2 && i2[p2] < i1[p1]) {
155: n3++;
156: p2++;
157: }
158: if (p2 == n2) {
159: /* cleanup for is1 */
160: n3 += n1 - p1;
161: break;
162: } else {
163: if (i2[p2] == i1[p1]) {
164: n3++;
165: p1++;
166: p2++;
167: }
168: }
169: }
170: if (p2 == n2) {
171: /* cleanup for is1 */
172: n3 += n1 - p1;
173: break;
174: } else {
175: while (p1 < n1 && i1[p1] < i2[p2]) {
176: n3++;
177: p1++;
178: }
179: if (p1 == n1) {
180: /* clean up for is2 */
181: n3 += n2 - p2;
182: break;
183: } else {
184: if (i1[p1] == i2[p2]) {
185: n3++;
186: p1++;
187: p2++;
188: }
189: }
190: }
191: } while (p1 < n1 || p2 < n2);
193: if (n3 == n1) { /* no new elements to be added */
194: ISRestoreIndices(is1, &i1);
195: ISRestoreIndices(is2, &i2);
196: ISDuplicate(is1, is3);
197: return 0;
198: }
199: PetscMalloc1(n3, &iout);
201: p1 = 0;
202: p2 = 0;
203: n3 = 0;
204: do {
205: if (p1 == n1) { /* cleanup for is2 */
206: while (p2 < n2) iout[n3++] = i2[p2++];
207: break;
208: } else {
209: while (p2 < n2 && i2[p2] < i1[p1]) iout[n3++] = i2[p2++];
210: if (p2 == n2) { /* cleanup for is1 */
211: while (p1 < n1) iout[n3++] = i1[p1++];
212: break;
213: } else {
214: if (i2[p2] == i1[p1]) {
215: iout[n3++] = i1[p1++];
216: p2++;
217: }
218: }
219: }
220: if (p2 == n2) { /* cleanup for is1 */
221: while (p1 < n1) iout[n3++] = i1[p1++];
222: break;
223: } else {
224: while (p1 < n1 && i1[p1] < i2[p2]) iout[n3++] = i1[p1++];
225: if (p1 == n1) { /* clean up for is2 */
226: while (p2 < n2) iout[n3++] = i2[p2++];
227: break;
228: } else {
229: if (i1[p1] == i2[p2]) {
230: iout[n3++] = i1[p1++];
231: p2++;
232: }
233: }
234: }
235: } while (p1 < n1 || p2 < n2);
237: ISRestoreIndices(is1, &i1);
238: ISRestoreIndices(is2, &i2);
239: ISCreateGeneral(comm, n3, iout, PETSC_OWN_POINTER, is3);
240: return 0;
241: }
243: /*@
244: ISExpand - Computes the union of two index sets, by concatenating 2 lists and
245: removing duplicates.
247: Collective on is1
249: Input Parameters:
250: + is1 - first index set
251: - is2 - index values to be added
253: Output Parameters:
254: . isout - is1 + is2 The index set is2 is appended to is1 removing duplicates
256: Level: intermediate
258: Notes:
259: Negative values are removed from the lists. This requires O(imax-imin)
260: memory and O(imax-imin) work, where imin and imax are the bounds on the
261: indices in is1 and is2.
263: The IS's do not need to be sorted.
265: The operations are performed separately on each MPI rank
267: .seealso: [](sec_scatter), `IS`, `ISDestroy()`, `ISView()`, `ISDifference()`, `ISSum()`, `ISIntersect()`
268: @*/
269: PetscErrorCode ISExpand(IS is1, IS is2, IS *isout)
270: {
271: PetscInt i, n1, n2, imin, imax, nout, *iout;
272: const PetscInt *i1, *i2;
273: PetscBT mask;
274: MPI_Comm comm;
281: if (!is1) {
282: ISDuplicate(is2, isout);
283: return 0;
284: }
285: if (!is2) {
286: ISDuplicate(is1, isout);
287: return 0;
288: }
289: ISGetIndices(is1, &i1);
290: ISGetLocalSize(is1, &n1);
291: ISGetIndices(is2, &i2);
292: ISGetLocalSize(is2, &n2);
294: /* Create a bit mask array to contain required values */
295: if (n1 || n2) {
296: imin = PETSC_MAX_INT;
297: imax = 0;
298: for (i = 0; i < n1; i++) {
299: if (i1[i] < 0) continue;
300: imin = PetscMin(imin, i1[i]);
301: imax = PetscMax(imax, i1[i]);
302: }
303: for (i = 0; i < n2; i++) {
304: if (i2[i] < 0) continue;
305: imin = PetscMin(imin, i2[i]);
306: imax = PetscMax(imax, i2[i]);
307: }
308: } else imin = imax = 0;
310: PetscMalloc1(n1 + n2, &iout);
311: nout = 0;
312: PetscBTCreate(imax - imin, &mask);
313: /* Put the values from is1 */
314: for (i = 0; i < n1; i++) {
315: if (i1[i] < 0) continue;
316: if (!PetscBTLookupSet(mask, i1[i] - imin)) iout[nout++] = i1[i];
317: }
318: ISRestoreIndices(is1, &i1);
319: /* Put the values from is2 */
320: for (i = 0; i < n2; i++) {
321: if (i2[i] < 0) continue;
322: if (!PetscBTLookupSet(mask, i2[i] - imin)) iout[nout++] = i2[i];
323: }
324: ISRestoreIndices(is2, &i2);
326: /* create the new IS containing the sum */
327: PetscObjectGetComm((PetscObject)is1, &comm);
328: ISCreateGeneral(comm, nout, iout, PETSC_OWN_POINTER, isout);
330: PetscBTDestroy(&mask);
331: return 0;
332: }
334: /*@
335: ISIntersect - Computes the intersection of two index sets, by sorting and comparing.
337: Collective on is1
339: Input Parameters:
340: + is1 - first index set
341: - is2 - second index set
343: Output Parameters:
344: . isout - the sorted intersection of is1 and is2
346: Level: intermediate
348: Notes:
349: Negative values are removed from the lists. This requires O(min(is1,is2))
350: memory and O(max(is1,is2)log(max(is1,is2))) work
352: The IS's do not need to be sorted.
354: The operations are performed separately on each MPI rank
356: .seealso: [](sec_scatter), `IS`, `ISDestroy()`, `ISView()`, `ISDifference()`, `ISSum()`, `ISExpand()`, `ISConcatenate()`
357: @*/
358: PetscErrorCode ISIntersect(IS is1, IS is2, IS *isout)
359: {
360: PetscInt i, n1, n2, nout, *iout;
361: const PetscInt *i1, *i2;
362: IS is1sorted = NULL, is2sorted = NULL;
363: PetscBool sorted, lsorted;
364: MPI_Comm comm;
370: PetscObjectGetComm((PetscObject)is1, &comm);
372: ISGetLocalSize(is1, &n1);
373: ISGetLocalSize(is2, &n2);
374: if (n1 < n2) {
375: IS tempis = is1;
376: PetscInt ntemp = n1;
378: is1 = is2;
379: is2 = tempis;
380: n1 = n2;
381: n2 = ntemp;
382: }
383: ISSorted(is1, &lsorted);
384: MPIU_Allreduce(&lsorted, &sorted, 1, MPIU_BOOL, MPI_LAND, comm);
385: if (!sorted) {
386: ISDuplicate(is1, &is1sorted);
387: ISSort(is1sorted);
388: ISGetIndices(is1sorted, &i1);
389: } else {
390: is1sorted = is1;
391: PetscObjectReference((PetscObject)is1);
392: ISGetIndices(is1, &i1);
393: }
394: ISSorted(is2, &lsorted);
395: MPIU_Allreduce(&lsorted, &sorted, 1, MPIU_BOOL, MPI_LAND, comm);
396: if (!sorted) {
397: ISDuplicate(is2, &is2sorted);
398: ISSort(is2sorted);
399: ISGetIndices(is2sorted, &i2);
400: } else {
401: is2sorted = is2;
402: PetscObjectReference((PetscObject)is2);
403: ISGetIndices(is2, &i2);
404: }
406: PetscMalloc1(n2, &iout);
408: for (i = 0, nout = 0; i < n2; i++) {
409: PetscInt key = i2[i];
410: PetscInt loc;
412: ISLocate(is1sorted, key, &loc);
413: if (loc >= 0) {
414: if (!nout || iout[nout - 1] < key) iout[nout++] = key;
415: }
416: }
417: PetscRealloc(nout * sizeof(PetscInt), &iout);
419: /* create the new IS containing the sum */
420: ISCreateGeneral(comm, nout, iout, PETSC_OWN_POINTER, isout);
422: ISRestoreIndices(is2sorted, &i2);
423: ISDestroy(&is2sorted);
424: ISRestoreIndices(is1sorted, &i1);
425: ISDestroy(&is1sorted);
426: return 0;
427: }
429: PetscErrorCode ISIntersect_Caching_Internal(IS is1, IS is2, IS *isect)
430: {
431: *isect = NULL;
432: if (is2 && is1) {
433: char composeStr[33] = {0};
434: PetscObjectId is2id;
436: PetscObjectGetId((PetscObject)is2, &is2id);
437: PetscSNPrintf(composeStr, 32, "ISIntersect_Caching_%" PetscInt64_FMT, is2id);
438: PetscObjectQuery((PetscObject)is1, composeStr, (PetscObject *)isect);
439: if (*isect == NULL) {
440: ISIntersect(is1, is2, isect);
441: PetscObjectCompose((PetscObject)is1, composeStr, (PetscObject)*isect);
442: } else {
443: PetscObjectReference((PetscObject)*isect);
444: }
445: }
446: return 0;
447: }
449: /*@
450: ISConcatenate - Forms a new `IS` by locally concatenating the indices from an `IS` list without reordering.
452: Collective.
454: Input Parameters:
455: + comm - communicator of the concatenated IS.
456: . len - size of islist array (nonnegative)
457: - islist - array of index sets
459: Output Parameters:
460: . isout - The concatenated index set; empty, if len == 0.
462: Level: intermediate
464: Notes:
465: The semantics of calling this on comm imply that the comms of the members of islist also contain this rank.
467: .seealso: [](sec_scatter), `IS`, `ISDifference()`, `ISSum()`, `ISExpand()`, `ISIntersect()`
468: @*/
469: PetscErrorCode ISConcatenate(MPI_Comm comm, PetscInt len, const IS islist[], IS *isout)
470: {
471: PetscInt i, n, N;
472: const PetscInt *iidx;
473: PetscInt *idx;
476: if (PetscDefined(USE_DEBUG)) {
477: for (i = 0; i < len; ++i)
479: }
481: if (!len) {
482: ISCreateGeneral(comm, 0, NULL, PETSC_OWN_POINTER, isout);
483: return 0;
484: }
486: N = 0;
487: for (i = 0; i < len; ++i) {
488: if (islist[i]) {
489: ISGetLocalSize(islist[i], &n);
490: N += n;
491: }
492: }
493: PetscMalloc1(N, &idx);
494: N = 0;
495: for (i = 0; i < len; ++i) {
496: if (islist[i]) {
497: ISGetLocalSize(islist[i], &n);
498: ISGetIndices(islist[i], &iidx);
499: PetscArraycpy(idx + N, iidx, n);
500: ISRestoreIndices(islist[i], &iidx);
501: N += n;
502: }
503: }
504: ISCreateGeneral(comm, N, idx, PETSC_OWN_POINTER, isout);
505: return 0;
506: }
508: /*@
509: ISListToPair - Convert an `IS` list to a pair of `IS` of equal length defining an equivalent integer multimap.
510: Each `IS` on the input list is assigned an integer j so that all of the indices of that `IS` are
511: mapped to j.
513: Collective.
515: Input Parameters:
516: + comm - `MPI_Comm`
517: . listlen - `IS` list length
518: - islist - `IS` list
520: Output Parameters:
521: + xis - domain `IS`
522: - yis - range `I`S
524: Level: developer
526: Notes:
527: The global integers assigned to the `IS` of the local input list might not correspond to the
528: local numbers of the `IS` on that list, but the two *orderings* are the same: the global
529: integers assigned to the `IS` on the local list form a strictly increasing sequence.
531: The `IS` on the input list can belong to subcommunicators of comm, and the subcommunicators
532: on the input `IS` list are assumed to be in a "deadlock-free" order.
534: Local lists of `PetscObject`s (or their subcomms) on a comm are "deadlock-free" if subcomm1
535: preceeds subcomm2 on any local list, then it preceeds subcomm2 on all ranks.
536: Equivalently, the local numbers of the subcomms on each local list are drawn from some global
537: numbering. This is ensured, for example, by `ISPairToList()`.
539: .seealso: `IS`, `ISPairToList()`
540: @*/
541: PetscErrorCode ISListToPair(MPI_Comm comm, PetscInt listlen, IS islist[], IS *xis, IS *yis)
542: {
543: PetscInt ncolors, *colors, i, leni, len, *xinds, *yinds, k, j;
544: const PetscInt *indsi;
546: PetscMalloc1(listlen, &colors);
547: PetscObjectsListGetGlobalNumbering(comm, listlen, (PetscObject *)islist, &ncolors, colors);
548: len = 0;
549: for (i = 0; i < listlen; ++i) {
550: ISGetLocalSize(islist[i], &leni);
551: len += leni;
552: }
553: PetscMalloc1(len, &xinds);
554: PetscMalloc1(len, &yinds);
555: k = 0;
556: for (i = 0; i < listlen; ++i) {
557: ISGetLocalSize(islist[i], &leni);
558: ISGetIndices(islist[i], &indsi);
559: for (j = 0; j < leni; ++j) {
560: xinds[k] = indsi[j];
561: yinds[k] = colors[i];
562: ++k;
563: }
564: }
565: PetscFree(colors);
566: ISCreateGeneral(comm, len, xinds, PETSC_OWN_POINTER, xis);
567: ISCreateGeneral(comm, len, yinds, PETSC_OWN_POINTER, yis);
568: return 0;
569: }
571: /*@
572: ISPairToList - Convert an `IS` pair encoding an integer map to a list of `IS`.
573: Each `IS` on the output list contains the preimage for each index on the second input `IS`.
574: The `IS` on the output list are constructed on the subcommunicators of the input `IS` pair.
575: Each subcommunicator corresponds to the preimage of some index j -- this subcomm contains
576: exactly the ranks that assign some indices i to j. This is essentially the inverse of
577: `ISListToPair()`.
579: Collective
581: Input Parameters:
582: + xis - domain IS
583: - yis - range IS
585: Output Parameters:
586: + listlen - length of islist
587: - islist - list of ISs breaking up indis by color
589: Level: developer
591: Notes:
592: xis and yis must be of the same length and have congruent communicators.
594: The resulting `IS` have subcommunicators in a "deadlock-free" order (see `ISListToPair()`).
596: .seealso: `IS`, `ISListToPair()`
597: @*/
598: PetscErrorCode ISPairToList(IS xis, IS yis, PetscInt *listlen, IS **islist)
599: {
600: IS indis = xis, coloris = yis;
601: PetscInt *inds, *colors, llen, ilen, lstart, lend, lcount, l;
602: PetscMPIInt rank, size, llow, lhigh, low, high, color, subsize;
603: const PetscInt *ccolors, *cinds;
604: MPI_Comm comm, subcomm;
611: PetscObjectGetComm((PetscObject)xis, &comm);
612: MPI_Comm_rank(comm, &rank);
613: MPI_Comm_rank(comm, &size);
614: /* Extract, copy and sort the local indices and colors on the color. */
615: ISGetLocalSize(coloris, &llen);
616: ISGetLocalSize(indis, &ilen);
618: ISGetIndices(coloris, &ccolors);
619: ISGetIndices(indis, &cinds);
620: PetscMalloc2(ilen, &inds, llen, &colors);
621: PetscArraycpy(inds, cinds, ilen);
622: PetscArraycpy(colors, ccolors, llen);
623: PetscSortIntWithArray(llen, colors, inds);
624: /* Determine the global extent of colors. */
625: llow = 0;
626: lhigh = -1;
627: lstart = 0;
628: lcount = 0;
629: while (lstart < llen) {
630: lend = lstart + 1;
631: while (lend < llen && colors[lend] == colors[lstart]) ++lend;
632: llow = PetscMin(llow, colors[lstart]);
633: lhigh = PetscMax(lhigh, colors[lstart]);
634: ++lcount;
635: }
636: MPIU_Allreduce(&llow, &low, 1, MPI_INT, MPI_MIN, comm);
637: MPIU_Allreduce(&lhigh, &high, 1, MPI_INT, MPI_MAX, comm);
638: *listlen = 0;
639: if (low <= high) {
640: if (lcount > 0) {
641: *listlen = lcount;
642: if (!*islist) PetscMalloc1(lcount, islist);
643: }
644: /*
645: Traverse all possible global colors, and participate in the subcommunicators
646: for the locally-supported colors.
647: */
648: lcount = 0;
649: lstart = 0;
650: lend = 0;
651: for (l = low; l <= high; ++l) {
652: /*
653: Find the range of indices with the same color, which is not smaller than l.
654: Observe that, since colors is sorted, and is a subsequence of [low,high],
655: as soon as we find a new color, it is >= l.
656: */
657: if (lstart < llen) {
658: /* The start of the next locally-owned color is identified. Now look for the end. */
659: if (lstart == lend) {
660: lend = lstart + 1;
661: while (lend < llen && colors[lend] == colors[lstart]) ++lend;
662: }
663: /* Now check whether the identified color segment matches l. */
665: }
666: color = (PetscMPIInt)(colors[lstart] == l);
667: /* Check whether a proper subcommunicator exists. */
668: MPIU_Allreduce(&color, &subsize, 1, MPI_INT, MPI_SUM, comm);
670: if (subsize == 1) subcomm = PETSC_COMM_SELF;
671: else if (subsize == size) subcomm = comm;
672: else {
673: /* a proper communicator is necessary, so we create it. */
674: MPI_Comm_split(comm, color, rank, &subcomm);
675: }
676: if (colors[lstart] == l) {
677: /* If we have l among the local colors, we create an IS to hold the corresponding indices. */
678: ISCreateGeneral(subcomm, lend - lstart, inds + lstart, PETSC_COPY_VALUES, *islist + lcount);
679: /* Position lstart at the beginning of the next local color. */
680: lstart = lend;
681: /* Increment the counter of the local colors split off into an IS. */
682: ++lcount;
683: }
684: if (subsize > 0 && subsize < size) {
685: /*
686: Irrespective of color, destroy the split off subcomm:
687: a subcomm used in the IS creation above is duplicated
688: into a proper PETSc comm.
689: */
690: MPI_Comm_free(&subcomm);
691: }
692: } /* for (l = low; l < high; ++l) */
693: } /* if (low <= high) */
694: PetscFree2(inds, colors);
695: return 0;
696: }
698: /*@
699: ISEmbed - Embed `IS` a into `IS` b by finding the locations in b that have the same indices as in a.
700: If c is the `IS` of these locations, we have a = b*c, regarded as a composition of the
701: corresponding `ISLocalToGlobalMap`s.
703: Not collective.
705: Input Parameters:
706: + a - `IS` to embed
707: . b - `IS` to embed into
708: - drop - flag indicating whether to drop a's indices that are not in b.
710: Output Parameters:
711: . c - local embedding indices
713: Level: developer
715: Notes:
716: If some of a's global indices are not among b's indices the embedding is impossible. The local indices of a
717: corresponding to these global indices are either mapped to -1 (if !drop) or are omitted (if drop). In the former
718: case the size of c is that same as that of a, in the latter case c's size may be smaller.
720: The resulting `IS` is sequential, since the index substitution it encodes is purely local.
722: .seealso: `IS`, `ISLocalToGlobalMapping`
723: @*/
724: PetscErrorCode ISEmbed(IS a, IS b, PetscBool drop, IS *c)
725: {
726: ISLocalToGlobalMapping ltog;
727: ISGlobalToLocalMappingMode gtoltype = IS_GTOLM_DROP;
728: PetscInt alen, clen, *cindices, *cindices2;
729: const PetscInt *aindices;
734: ISLocalToGlobalMappingCreateIS(b, <og);
735: ISGetLocalSize(a, &alen);
736: ISGetIndices(a, &aindices);
737: PetscMalloc1(alen, &cindices);
738: if (!drop) gtoltype = IS_GTOLM_MASK;
739: ISGlobalToLocalMappingApply(ltog, gtoltype, alen, aindices, &clen, cindices);
740: ISRestoreIndices(a, &aindices);
741: ISLocalToGlobalMappingDestroy(<og);
742: if (clen != alen) {
743: cindices2 = cindices;
744: PetscMalloc1(clen, &cindices);
745: PetscArraycpy(cindices, cindices2, clen);
746: PetscFree(cindices2);
747: }
748: ISCreateGeneral(PETSC_COMM_SELF, clen, cindices, PETSC_OWN_POINTER, c);
749: return 0;
750: }
752: /*@
753: ISSortPermutation - calculate the permutation of the indices into a nondecreasing order.
755: Not collective.
757: Input Parameters:
758: + f - `IS` to sort
759: - always - build the permutation even when f's indices are nondecreasing.
761: Output Parameter:
762: . h - permutation or NULL, if f is nondecreasing and always == `PETSC_FALSE`.
764: Level: advanced
766: Notes:
767: Indices in f are unchanged. f[h[i]] is the i-th smallest f index.
769: If always == `PETSC_FALSE`, an extra check is performed to see whether
770: the f indices are nondecreasing. h is built on `PETSC_COMM_SELF`, since
771: the permutation has a local meaning only.
773: .seealso: `IS`, `ISLocalToGlobalMapping`, `ISSort()`
774: @*/
775: PetscErrorCode ISSortPermutation(IS f, PetscBool always, IS *h)
776: {
777: const PetscInt *findices;
778: PetscInt fsize, *hindices, i;
779: PetscBool isincreasing;
783: ISGetLocalSize(f, &fsize);
784: ISGetIndices(f, &findices);
785: *h = NULL;
786: if (!always) {
787: isincreasing = PETSC_TRUE;
788: for (i = 1; i < fsize; ++i) {
789: if (findices[i] <= findices[i - 1]) {
790: isincreasing = PETSC_FALSE;
791: break;
792: }
793: }
794: if (isincreasing) {
795: ISRestoreIndices(f, &findices);
796: return 0;
797: }
798: }
799: PetscMalloc1(fsize, &hindices);
800: for (i = 0; i < fsize; ++i) hindices[i] = i;
801: PetscSortIntWithPermutation(fsize, findices, hindices);
802: ISRestoreIndices(f, &findices);
803: ISCreateGeneral(PETSC_COMM_SELF, fsize, hindices, PETSC_OWN_POINTER, h);
804: ISSetInfo(*h, IS_PERMUTATION, IS_LOCAL, PETSC_FALSE, PETSC_TRUE);
805: return 0;
806: }