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 IS
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: Notes:
19: Negative values are removed from the lists. is2 may have values
20: that are not in is1. This requires O(imax-imin) memory and O(imax-imin)
21: work, where imin and imax are the bounds on the indices in is1.
23: If is2 is NULL, the result is the same as for an empty IS, i.e., a duplicate of is1.
25: Level: intermediate
27: .seealso: ISDestroy(), ISView(), ISSum(), ISExpand()
28: @*/
29: PetscErrorCode ISDifference(IS is1,IS is2,IS *isout)
30: {
31: PetscInt i,n1,n2,imin,imax,nout,*iout;
32: const PetscInt *i1,*i2;
33: PetscBT mask;
34: MPI_Comm comm;
38: if (!is2) {
39: ISDuplicate(is1, isout);
40: return 0;
41: }
44: ISGetIndices(is1,&i1);
45: ISGetLocalSize(is1,&n1);
47: /* Create a bit mask array to contain required values */
48: if (n1) {
49: imin = PETSC_MAX_INT;
50: imax = 0;
51: for (i=0; i<n1; i++) {
52: if (i1[i] < 0) continue;
53: imin = PetscMin(imin,i1[i]);
54: imax = PetscMax(imax,i1[i]);
55: }
56: } else imin = imax = 0;
58: PetscBTCreate(imax-imin,&mask);
59: /* Put the values from is1 */
60: for (i=0; i<n1; i++) {
61: if (i1[i] < 0) continue;
62: PetscBTSet(mask,i1[i] - imin);
63: }
64: ISRestoreIndices(is1,&i1);
65: /* Remove the values from is2 */
66: ISGetIndices(is2,&i2);
67: ISGetLocalSize(is2,&n2);
68: for (i=0; i<n2; i++) {
69: if (i2[i] < imin || i2[i] > imax) continue;
70: PetscBTClear(mask,i2[i] - imin);
71: }
72: ISRestoreIndices(is2,&i2);
74: /* Count the number in the difference */
75: nout = 0;
76: for (i=0; i<imax-imin+1; i++) {
77: if (PetscBTLookup(mask,i)) nout++;
78: }
80: /* create the new IS containing the difference */
81: PetscMalloc1(nout,&iout);
82: nout = 0;
83: for (i=0; i<imax-imin+1; i++) {
84: if (PetscBTLookup(mask,i)) iout[nout++] = i + imin;
85: }
86: PetscObjectGetComm((PetscObject)is1,&comm);
87: ISCreateGeneral(comm,nout,iout,PETSC_OWN_POINTER,isout);
89: PetscBTDestroy(&mask);
90: return 0;
91: }
93: /*@
94: ISSum - Computes the sum (union) of two index sets.
96: Only sequential version (at the moment)
98: Input Parameters:
99: + is1 - index set to be extended
100: - is2 - index values to be added
102: Output Parameter:
103: . is3 - the sum; this can not be is1 or is2
105: Notes:
106: If n1 and n2 are the sizes of the sets, this takes O(n1+n2) time;
108: Both index sets need to be sorted on input.
110: Level: intermediate
112: .seealso: ISDestroy(), ISView(), ISDifference(), ISExpand()
114: @*/
115: PetscErrorCode ISSum(IS is1,IS is2,IS *is3)
116: {
117: MPI_Comm comm;
118: PetscBool f;
119: PetscMPIInt size;
120: const PetscInt *i1,*i2;
121: PetscInt n1,n2,n3, p1,p2, *iout;
125: PetscObjectGetComm((PetscObject)(is1),&comm);
126: MPI_Comm_size(comm,&size);
129: ISSorted(is1,&f);
131: ISSorted(is2,&f);
134: ISGetLocalSize(is1,&n1);
135: ISGetLocalSize(is2,&n2);
136: if (!n2) {
137: ISDuplicate(is1,is3);
138: return 0;
139: }
140: ISGetIndices(is1,&i1);
141: ISGetIndices(is2,&i2);
143: p1 = 0; p2 = 0; n3 = 0;
144: do {
145: if (p1==n1) { /* cleanup for is2 */ n3 += n2-p2; break;
146: } else {
147: while (p2<n2 && i2[p2]<i1[p1]) {
148: n3++; p2++;
149: }
150: if (p2==n2) {
151: /* cleanup for is1 */
152: n3 += n1-p1; break;
153: } else {
154: if (i2[p2]==i1[p1]) { n3++; p1++; p2++; }
155: }
156: }
157: if (p2==n2) {
158: /* cleanup for is1 */
159: n3 += n1-p1; break;
160: } else {
161: while (p1<n1 && i1[p1]<i2[p2]) {
162: n3++; p1++;
163: }
164: if (p1==n1) {
165: /* clean up for is2 */
166: n3 += n2-p2; break;
167: } else {
168: if (i1[p1]==i2[p2]) { n3++; p1++; p2++; }
169: }
170: }
171: } while (p1<n1 || p2<n2);
173: if (n3==n1) { /* no new elements to be added */
174: ISRestoreIndices(is1,&i1);
175: ISRestoreIndices(is2,&i2);
176: ISDuplicate(is1,is3);
177: return 0;
178: }
179: PetscMalloc1(n3,&iout);
181: p1 = 0; p2 = 0; n3 = 0;
182: do {
183: if (p1==n1) { /* cleanup for is2 */
184: while (p2<n2) iout[n3++] = i2[p2++];
185: break;
186: } else {
187: while (p2<n2 && i2[p2]<i1[p1]) iout[n3++] = i2[p2++];
188: if (p2==n2) { /* cleanup for is1 */
189: while (p1<n1) iout[n3++] = i1[p1++];
190: break;
191: } else {
192: if (i2[p2]==i1[p1]) { iout[n3++] = i1[p1++]; p2++; }
193: }
194: }
195: if (p2==n2) { /* cleanup for is1 */
196: while (p1<n1) iout[n3++] = i1[p1++];
197: break;
198: } else {
199: while (p1<n1 && i1[p1]<i2[p2]) iout[n3++] = i1[p1++];
200: if (p1==n1) { /* clean up for is2 */
201: while (p2<n2) iout[n3++] = i2[p2++];
202: break;
203: } else {
204: if (i1[p1]==i2[p2]) { iout[n3++] = i1[p1++]; p2++; }
205: }
206: }
207: } while (p1<n1 || p2<n2);
209: ISRestoreIndices(is1,&i1);
210: ISRestoreIndices(is2,&i2);
211: ISCreateGeneral(comm,n3,iout,PETSC_OWN_POINTER,is3);
212: return 0;
213: }
215: /*@
216: ISExpand - Computes the union of two index sets, by concatenating 2 lists and
217: removing duplicates.
219: Collective on IS
221: Input Parameters:
222: + is1 - first index set
223: - is2 - index values to be added
225: Output Parameters:
226: . isout - is1 + is2 The index set is2 is appended to is1 removing duplicates
228: Notes:
229: Negative values are removed from the lists. This requires O(imax-imin)
230: memory and O(imax-imin) work, where imin and imax are the bounds on the
231: indices in is1 and is2.
233: The IS's do not need to be sorted.
235: Level: intermediate
237: .seealso: ISDestroy(), ISView(), ISDifference(), ISSum()
239: @*/
240: PetscErrorCode ISExpand(IS is1,IS is2,IS *isout)
241: {
242: PetscInt i,n1,n2,imin,imax,nout,*iout;
243: const PetscInt *i1,*i2;
244: PetscBT mask;
245: MPI_Comm comm;
252: if (!is1) {ISDuplicate(is2, isout));PetscFunctionReturn(0;}
253: if (!is2) {ISDuplicate(is1, isout));PetscFunctionReturn(0;}
254: ISGetIndices(is1,&i1);
255: ISGetLocalSize(is1,&n1);
256: ISGetIndices(is2,&i2);
257: ISGetLocalSize(is2,&n2);
259: /* Create a bit mask array to contain required values */
260: if (n1 || n2) {
261: imin = PETSC_MAX_INT;
262: imax = 0;
263: for (i=0; i<n1; i++) {
264: if (i1[i] < 0) continue;
265: imin = PetscMin(imin,i1[i]);
266: imax = PetscMax(imax,i1[i]);
267: }
268: for (i=0; i<n2; i++) {
269: if (i2[i] < 0) continue;
270: imin = PetscMin(imin,i2[i]);
271: imax = PetscMax(imax,i2[i]);
272: }
273: } else imin = imax = 0;
275: PetscMalloc1(n1+n2,&iout);
276: nout = 0;
277: PetscBTCreate(imax-imin,&mask);
278: /* Put the values from is1 */
279: for (i=0; i<n1; i++) {
280: if (i1[i] < 0) continue;
281: if (!PetscBTLookupSet(mask,i1[i] - imin)) iout[nout++] = i1[i];
282: }
283: ISRestoreIndices(is1,&i1);
284: /* Put the values from is2 */
285: for (i=0; i<n2; i++) {
286: if (i2[i] < 0) continue;
287: if (!PetscBTLookupSet(mask,i2[i] - imin)) iout[nout++] = i2[i];
288: }
289: ISRestoreIndices(is2,&i2);
291: /* create the new IS containing the sum */
292: PetscObjectGetComm((PetscObject)is1,&comm);
293: ISCreateGeneral(comm,nout,iout,PETSC_OWN_POINTER,isout);
295: PetscBTDestroy(&mask);
296: return 0;
297: }
299: /*@
300: ISIntersect - Computes the intersection of two index sets, by sorting and comparing.
302: Collective on IS
304: Input Parameters:
305: + is1 - first index set
306: - is2 - second index set
308: Output Parameters:
309: . isout - the sorted intersection of is1 and is2
311: Notes:
312: Negative values are removed from the lists. This requires O(min(is1,is2))
313: memory and O(max(is1,is2)log(max(is1,is2))) work
315: The IS's do not need to be sorted.
317: Level: intermediate
319: .seealso: ISDestroy(), ISView(), ISDifference(), ISSum(), ISExpand()
320: @*/
321: PetscErrorCode ISIntersect(IS is1,IS is2,IS *isout)
322: {
323: PetscInt i,n1,n2,nout,*iout;
324: const PetscInt *i1,*i2;
325: IS is1sorted = NULL, is2sorted = NULL;
326: PetscBool sorted, lsorted;
327: MPI_Comm comm;
333: PetscObjectGetComm((PetscObject)is1,&comm);
335: ISGetLocalSize(is1,&n1);
336: ISGetLocalSize(is2,&n2);
337: if (n1 < n2) {
338: IS tempis = is1;
339: PetscInt ntemp = n1;
341: is1 = is2;
342: is2 = tempis;
343: n1 = n2;
344: n2 = ntemp;
345: }
346: ISSorted(is1,&lsorted);
347: MPIU_Allreduce(&lsorted,&sorted,1,MPIU_BOOL,MPI_LAND,comm);
348: if (!sorted) {
349: ISDuplicate(is1,&is1sorted);
350: ISSort(is1sorted);
351: ISGetIndices(is1sorted,&i1);
352: } else {
353: is1sorted = is1;
354: PetscObjectReference((PetscObject)is1);
355: ISGetIndices(is1,&i1);
356: }
357: ISSorted(is2,&lsorted);
358: MPIU_Allreduce(&lsorted,&sorted,1,MPIU_BOOL,MPI_LAND,comm);
359: if (!sorted) {
360: ISDuplicate(is2,&is2sorted);
361: ISSort(is2sorted);
362: ISGetIndices(is2sorted,&i2);
363: } else {
364: is2sorted = is2;
365: PetscObjectReference((PetscObject)is2);
366: ISGetIndices(is2,&i2);
367: }
369: PetscMalloc1(n2,&iout);
371: for (i = 0, nout = 0; i < n2; i++) {
372: PetscInt key = i2[i];
373: PetscInt loc;
375: ISLocate(is1sorted,key,&loc);
376: if (loc >= 0) {
377: if (!nout || iout[nout-1] < key) {
378: iout[nout++] = key;
379: }
380: }
381: }
382: PetscRealloc(nout*sizeof(PetscInt),&iout);
384: /* create the new IS containing the sum */
385: ISCreateGeneral(comm,nout,iout,PETSC_OWN_POINTER,isout);
387: ISRestoreIndices(is2sorted,&i2);
388: ISDestroy(&is2sorted);
389: ISRestoreIndices(is1sorted,&i1);
390: ISDestroy(&is1sorted);
391: return 0;
392: }
394: PetscErrorCode ISIntersect_Caching_Internal(IS is1, IS is2, IS *isect)
395: {
396: *isect = NULL;
397: if (is2 && is1) {
398: char composeStr[33] = {0};
399: PetscObjectId is2id;
401: PetscObjectGetId((PetscObject)is2,&is2id);
402: PetscSNPrintf(composeStr,32,"ISIntersect_Caching_%" PetscInt64_FMT,is2id);
403: PetscObjectQuery((PetscObject) is1, composeStr, (PetscObject *) isect);
404: if (*isect == NULL) {
405: ISIntersect(is1, is2, isect);
406: PetscObjectCompose((PetscObject) is1, composeStr, (PetscObject) *isect);
407: } else {
408: PetscObjectReference((PetscObject) *isect);
409: }
410: }
411: return 0;
412: }
414: /*@
415: ISConcatenate - Forms a new IS by locally concatenating the indices from an IS list without reordering.
417: Collective.
419: Input Parameters:
420: + comm - communicator of the concatenated IS.
421: . len - size of islist array (nonnegative)
422: - islist - array of index sets
424: Output Parameters:
425: . isout - The concatenated index set; empty, if len == 0.
427: Notes:
428: The semantics of calling this on comm imply that the comms of the members if islist also contain this rank.
430: Level: intermediate
432: .seealso: ISDifference(), ISSum(), ISExpand()
434: @*/
435: PetscErrorCode ISConcatenate(MPI_Comm comm, PetscInt len, const IS islist[], IS *isout)
436: {
437: PetscInt i,n,N;
438: const PetscInt *iidx;
439: PetscInt *idx;
442: if (PetscDefined(USE_DEBUG)) {
444: }
446: if (!len) {
447: ISCreateStride(comm, 0,0,0, isout);
448: return 0;
449: }
451: N = 0;
452: for (i = 0; i < len; ++i) {
453: if (islist[i]) {
454: ISGetLocalSize(islist[i], &n);
455: N += n;
456: }
457: }
458: PetscMalloc1(N, &idx);
459: N = 0;
460: for (i = 0; i < len; ++i) {
461: if (islist[i]) {
462: ISGetLocalSize(islist[i], &n);
463: ISGetIndices(islist[i], &iidx);
464: PetscArraycpy(idx+N,iidx, n);
465: ISRestoreIndices(islist[i], &iidx);
466: N += n;
467: }
468: }
469: ISCreateGeneral(comm, N, idx, PETSC_OWN_POINTER, isout);
470: return 0;
471: }
473: /*@
474: ISListToPair - convert an IS list to a pair of ISs of equal length defining an equivalent integer multimap.
475: Each IS on the input list is assigned an integer j so that all of the indices of that IS are
476: mapped to j.
478: Collective.
480: Input arguments:
481: + comm - MPI_Comm
482: . listlen - IS list length
483: - islist - IS list
485: Output arguments:
486: + xis - domain IS
487: - yis - range IS
489: Level: advanced
491: Notes:
492: The global integers assigned to the ISs of the local input list might not correspond to the
493: local numbers of the ISs on that list, but the two *orderings* are the same: the global
494: integers assigned to the ISs on the local list form a strictly increasing sequence.
496: The ISs on the input list can belong to subcommunicators of comm, and the subcommunicators
497: on the input IS list are assumed to be in a "deadlock-free" order.
499: Local lists of PetscObjects (or their subcommes) on a comm are "deadlock-free" if subcomm1
500: preceeds subcomm2 on any local list, then it preceeds subcomm2 on all ranks.
501: Equivalently, the local numbers of the subcomms on each local list are drawn from some global
502: numbering. This is ensured, for example, by ISPairToList().
504: .seealso ISPairToList()
505: @*/
506: PetscErrorCode ISListToPair(MPI_Comm comm, PetscInt listlen, IS islist[], IS *xis, IS *yis)
507: {
508: PetscInt ncolors, *colors,i, leni,len,*xinds, *yinds,k,j;
509: const PetscInt *indsi;
511: PetscMalloc1(listlen, &colors);
512: PetscObjectsListGetGlobalNumbering(comm, listlen, (PetscObject*)islist,&ncolors, colors);
513: len = 0;
514: for (i = 0; i < listlen; ++i) {
515: ISGetLocalSize(islist[i], &leni);
516: len += leni;
517: }
518: PetscMalloc1(len, &xinds);
519: PetscMalloc1(len, &yinds);
520: k = 0;
521: for (i = 0; i < listlen; ++i) {
522: ISGetLocalSize(islist[i], &leni);
523: ISGetIndices(islist[i],&indsi);
524: for (j = 0; j < leni; ++j) {
525: xinds[k] = indsi[j];
526: yinds[k] = colors[i];
527: ++k;
528: }
529: }
530: PetscFree(colors);
531: ISCreateGeneral(comm,len,xinds,PETSC_OWN_POINTER,xis);
532: ISCreateGeneral(comm,len,yinds,PETSC_OWN_POINTER,yis);
533: return 0;
534: }
536: /*@
537: ISPairToList - convert an IS pair encoding an integer map to a list of ISs.
538: Each IS on the output list contains the preimage for each index on the second input IS.
539: The ISs on the output list are constructed on the subcommunicators of the input IS pair.
540: Each subcommunicator corresponds to the preimage of some index j -- this subcomm contains
541: exactly the ranks that assign some indices i to j. This is essentially the inverse of
542: ISListToPair().
544: Collective on indis.
546: Input arguments:
547: + xis - domain IS
548: - yis - range IS
550: Output arguments:
551: + listlen - length of islist
552: - islist - list of ISs breaking up indis by color
554: Note:
555: xis and yis must be of the same length and have congruent communicators.
557: The resulting ISs have subcommunicators in a "deadlock-free" order (see ISListToPair()).
559: Level: advanced
561: .seealso ISListToPair()
562: @*/
563: PetscErrorCode ISPairToList(IS xis, IS yis, PetscInt *listlen, IS **islist)
564: {
565: IS indis = xis, coloris = yis;
566: PetscInt *inds, *colors, llen, ilen, lstart, lend, lcount,l;
567: PetscMPIInt rank, size, llow, lhigh, low, high,color,subsize;
568: const PetscInt *ccolors, *cinds;
569: MPI_Comm comm, subcomm;
576: PetscObjectGetComm((PetscObject)xis,&comm);
577: MPI_Comm_rank(comm, &rank);
578: MPI_Comm_rank(comm, &size);
579: /* Extract, copy and sort the local indices and colors on the color. */
580: ISGetLocalSize(coloris, &llen);
581: ISGetLocalSize(indis, &ilen);
583: ISGetIndices(coloris, &ccolors);
584: ISGetIndices(indis, &cinds);
585: PetscMalloc2(ilen,&inds,llen,&colors);
586: PetscArraycpy(inds,cinds,ilen);
587: PetscArraycpy(colors,ccolors,llen);
588: PetscSortIntWithArray(llen, colors, inds);
589: /* Determine the global extent of colors. */
590: llow = 0; lhigh = -1;
591: lstart = 0; lcount = 0;
592: while (lstart < llen) {
593: lend = lstart+1;
594: while (lend < llen && colors[lend] == colors[lstart]) ++lend;
595: llow = PetscMin(llow,colors[lstart]);
596: lhigh = PetscMax(lhigh,colors[lstart]);
597: ++lcount;
598: }
599: MPIU_Allreduce(&llow,&low,1,MPI_INT,MPI_MIN,comm);
600: MPIU_Allreduce(&lhigh,&high,1,MPI_INT,MPI_MAX,comm);
601: *listlen = 0;
602: if (low <= high) {
603: if (lcount > 0) {
604: *listlen = lcount;
605: if (!*islist) {
606: PetscMalloc1(lcount, islist);
607: }
608: }
609: /*
610: Traverse all possible global colors, and participate in the subcommunicators
611: for the locally-supported colors.
612: */
613: lcount = 0;
614: lstart = 0; lend = 0;
615: for (l = low; l <= high; ++l) {
616: /*
617: Find the range of indices with the same color, which is not smaller than l.
618: Observe that, since colors is sorted, and is a subsequence of [low,high],
619: as soon as we find a new color, it is >= l.
620: */
621: if (lstart < llen) {
622: /* The start of the next locally-owned color is identified. Now look for the end. */
623: if (lstart == lend) {
624: lend = lstart+1;
625: while (lend < llen && colors[lend] == colors[lstart]) ++lend;
626: }
627: /* Now check whether the identified color segment matches l. */
629: }
630: color = (PetscMPIInt)(colors[lstart] == l);
631: /* Check whether a proper subcommunicator exists. */
632: MPIU_Allreduce(&color,&subsize,1,MPI_INT,MPI_SUM,comm);
634: if (subsize == 1) subcomm = PETSC_COMM_SELF;
635: else if (subsize == size) subcomm = comm;
636: else {
637: /* a proper communicator is necessary, so we create it. */
638: MPI_Comm_split(comm, color, rank, &subcomm);
639: }
640: if (colors[lstart] == l) {
641: /* If we have l among the local colors, we create an IS to hold the corresponding indices. */
642: ISCreateGeneral(subcomm, lend-lstart,inds+lstart,PETSC_COPY_VALUES,*islist+lcount);
643: /* Position lstart at the beginning of the next local color. */
644: lstart = lend;
645: /* Increment the counter of the local colors split off into an IS. */
646: ++lcount;
647: }
648: if (subsize > 0 && subsize < size) {
649: /*
650: Irrespective of color, destroy the split off subcomm:
651: a subcomm used in the IS creation above is duplicated
652: into a proper PETSc comm.
653: */
654: MPI_Comm_free(&subcomm);
655: }
656: } /* for (l = low; l < high; ++l) */
657: } /* if (low <= high) */
658: PetscFree2(inds,colors);
659: return 0;
660: }
662: /*@
663: ISEmbed - embed IS a into IS b by finding the locations in b that have the same indices as in a.
664: If c is the IS of these locations, we have a = b*c, regarded as a composition of the
665: corresponding ISLocalToGlobalMaps.
667: Not collective.
669: Input arguments:
670: + a - IS to embed
671: . b - IS to embed into
672: - drop - flag indicating whether to drop a's indices that are not in b.
674: Output arguments:
675: . c - local embedding indices
677: Note:
678: If some of a's global indices are not among b's indices the embedding is impossible. The local indices of a
679: corresponding to these global indices are either mapped to -1 (if !drop) or are omitted (if drop). In the former
680: case the size of c is that same as that of a, in the latter case c's size may be smaller.
682: The resulting IS is sequential, since the index substition it encodes is purely local.
684: Level: advanced
686: .seealso ISLocalToGlobalMapping
687: @*/
688: PetscErrorCode ISEmbed(IS a, IS b, PetscBool drop, IS *c)
689: {
690: ISLocalToGlobalMapping ltog;
691: ISGlobalToLocalMappingMode gtoltype = IS_GTOLM_DROP;
692: PetscInt alen, clen, *cindices, *cindices2;
693: const PetscInt *aindices;
698: ISLocalToGlobalMappingCreateIS(b, <og);
699: ISGetLocalSize(a, &alen);
700: ISGetIndices(a, &aindices);
701: PetscMalloc1(alen, &cindices);
702: if (!drop) gtoltype = IS_GTOLM_MASK;
703: ISGlobalToLocalMappingApply(ltog,gtoltype,alen,aindices,&clen,cindices);
704: ISLocalToGlobalMappingDestroy(<og);
705: if (clen != alen) {
706: cindices2 = cindices;
707: PetscMalloc1(clen, &cindices);
708: PetscArraycpy(cindices,cindices2,clen);
709: PetscFree(cindices2);
710: }
711: ISCreateGeneral(PETSC_COMM_SELF,clen,cindices,PETSC_OWN_POINTER,c);
712: return 0;
713: }
715: /*@
716: ISSortPermutation - calculate the permutation of the indices into a nondecreasing order.
718: Not collective.
720: Input arguments:
721: + f - IS to sort
722: - always - build the permutation even when f's indices are nondecreasing.
724: Output argument:
725: . h - permutation or NULL, if f is nondecreasing and always == PETSC_FALSE.
727: Note: Indices in f are unchanged. f[h[i]] is the i-th smallest f index.
728: If always == PETSC_FALSE, an extra check is peformed to see whether
729: the f indices are nondecreasing. h is built on PETSC_COMM_SELF, since
730: the permutation has a local meaning only.
732: Level: advanced
734: .seealso ISLocalToGlobalMapping, ISSort()
735: @*/
736: PetscErrorCode ISSortPermutation(IS f,PetscBool always,IS *h)
737: {
738: const PetscInt *findices;
739: PetscInt fsize,*hindices,i;
740: PetscBool isincreasing;
744: ISGetLocalSize(f,&fsize);
745: ISGetIndices(f,&findices);
746: *h = NULL;
747: if (!always) {
748: isincreasing = PETSC_TRUE;
749: for (i = 1; i < fsize; ++i) {
750: if (findices[i] <= findices[i-1]) {
751: isincreasing = PETSC_FALSE;
752: break;
753: }
754: }
755: if (isincreasing) {
756: ISRestoreIndices(f,&findices);
757: return 0;
758: }
759: }
760: PetscMalloc1(fsize,&hindices);
761: for (i = 0; i < fsize; ++i) hindices[i] = i;
762: PetscSortIntWithPermutation(fsize,findices,hindices);
763: ISRestoreIndices(f,&findices);
764: ISCreateGeneral(PETSC_COMM_SELF,fsize,hindices,PETSC_OWN_POINTER,h);
765: ISSetInfo(*h,IS_PERMUTATION,IS_LOCAL,PETSC_FALSE,PETSC_TRUE);
766: return 0;
767: }