Actual source code: isdiff.c
petsc-master 2020-02-23
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 Parameter:
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: Level: intermediate
25: .seealso: ISDestroy(), ISView(), ISSum(), ISExpand()
26: @*/
27: PetscErrorCode ISDifference(IS is1,IS is2,IS *isout)
28: {
30: PetscInt i,n1,n2,imin,imax,nout,*iout;
31: const PetscInt *i1,*i2;
32: PetscBT mask;
33: MPI_Comm comm;
40: ISGetIndices(is1,&i1);
41: ISGetLocalSize(is1,&n1);
43: /* Create a bit mask array to contain required values */
44: if (n1) {
45: imin = PETSC_MAX_INT;
46: imax = 0;
47: for (i=0; i<n1; i++) {
48: if (i1[i] < 0) continue;
49: imin = PetscMin(imin,i1[i]);
50: imax = PetscMax(imax,i1[i]);
51: }
52: } else imin = imax = 0;
54: PetscBTCreate(imax-imin,&mask);
55: /* Put the values from is1 */
56: for (i=0; i<n1; i++) {
57: if (i1[i] < 0) continue;
58: PetscBTSet(mask,i1[i] - imin);
59: }
60: ISRestoreIndices(is1,&i1);
61: /* Remove the values from is2 */
62: ISGetIndices(is2,&i2);
63: ISGetLocalSize(is2,&n2);
64: for (i=0; i<n2; i++) {
65: if (i2[i] < imin || i2[i] > imax) continue;
66: PetscBTClear(mask,i2[i] - imin);
67: }
68: ISRestoreIndices(is2,&i2);
70: /* Count the number in the difference */
71: nout = 0;
72: for (i=0; i<imax-imin+1; i++) {
73: if (PetscBTLookup(mask,i)) nout++;
74: }
76: /* create the new IS containing the difference */
77: PetscMalloc1(nout,&iout);
78: nout = 0;
79: for (i=0; i<imax-imin+1; i++) {
80: if (PetscBTLookup(mask,i)) iout[nout++] = i + imin;
81: }
82: PetscObjectGetComm((PetscObject)is1,&comm);
83: ISCreateGeneral(comm,nout,iout,PETSC_OWN_POINTER,isout);
85: PetscBTDestroy(&mask);
86: return(0);
87: }
89: /*@
90: ISSum - Computes the sum (union) of two index sets.
92: Only sequential version (at the moment)
94: Input Parameter:
95: + is1 - index set to be extended
96: - is2 - index values to be added
98: Output Parameter:
99: . is3 - the sum; this can not be is1 or is2
101: Notes:
102: If n1 and n2 are the sizes of the sets, this takes O(n1+n2) time;
104: Both index sets need to be sorted on input.
106: Level: intermediate
108: .seealso: ISDestroy(), ISView(), ISDifference(), ISExpand()
111: @*/
112: PetscErrorCode ISSum(IS is1,IS is2,IS *is3)
113: {
114: MPI_Comm comm;
115: PetscBool f;
116: PetscMPIInt size;
117: const PetscInt *i1,*i2;
118: PetscInt n1,n2,n3, p1,p2, *iout;
124: PetscObjectGetComm((PetscObject)(is1),&comm);
125: MPI_Comm_size(comm,&size);
126: if (size>1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently only for uni-processor IS");
128: ISSorted(is1,&f);
129: if (!f) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Arg 1 is not sorted");
130: ISSorted(is2,&f);
131: if (!f) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Arg 2 is not sorted");
133: ISGetLocalSize(is1,&n1);
134: ISGetLocalSize(is2,&n2);
135: if (!n2) {
136: ISDuplicate(is1,is3);
137: return(0);
138: }
139: ISGetIndices(is1,&i1);
140: ISGetIndices(is2,&i2);
142: p1 = 0; p2 = 0; n3 = 0;
143: do {
144: if (p1==n1) { /* cleanup for is2 */ n3 += n2-p2; break;
145: } else {
146: while (p2<n2 && i2[p2]<i1[p1]) {
147: n3++; p2++;
148: }
149: if (p2==n2) {
150: /* cleanup for is1 */
151: n3 += n1-p1; break;
152: } else {
153: if (i2[p2]==i1[p1]) { n3++; p1++; p2++; }
154: }
155: }
156: if (p2==n2) {
157: /* cleanup for is1 */
158: n3 += n1-p1; break;
159: } else {
160: while (p1<n1 && i1[p1]<i2[p2]) {
161: n3++; p1++;
162: }
163: if (p1==n1) {
164: /* clean up for is2 */
165: n3 += n2-p2; break;
166: } else {
167: if (i1[p1]==i2[p2]) { n3++; p1++; p2++; }
168: }
169: }
170: } while (p1<n1 || p2<n2);
172: if (n3==n1) { /* no new elements to be added */
173: ISRestoreIndices(is1,&i1);
174: ISRestoreIndices(is2,&i2);
175: ISDuplicate(is1,is3);
176: return(0);
177: }
178: PetscMalloc1(n3,&iout);
180: p1 = 0; p2 = 0; n3 = 0;
181: do {
182: if (p1==n1) { /* cleanup for is2 */
183: while (p2<n2) iout[n3++] = i2[p2++];
184: break;
185: } else {
186: while (p2<n2 && i2[p2]<i1[p1]) iout[n3++] = i2[p2++];
187: if (p2==n2) { /* cleanup for is1 */
188: while (p1<n1) iout[n3++] = i1[p1++];
189: break;
190: } else {
191: if (i2[p2]==i1[p1]) { iout[n3++] = i1[p1++]; p2++; }
192: }
193: }
194: if (p2==n2) { /* cleanup for is1 */
195: while (p1<n1) iout[n3++] = i1[p1++];
196: break;
197: } else {
198: while (p1<n1 && i1[p1]<i2[p2]) iout[n3++] = i1[p1++];
199: if (p1==n1) { /* clean up for is2 */
200: while (p2<n2) iout[n3++] = i2[p2++];
201: break;
202: } else {
203: if (i1[p1]==i2[p2]) { iout[n3++] = i1[p1++]; p2++; }
204: }
205: }
206: } while (p1<n1 || p2<n2);
208: ISRestoreIndices(is1,&i1);
209: ISRestoreIndices(is2,&i2);
210: ISCreateGeneral(comm,n3,iout,PETSC_OWN_POINTER,is3);
211: return(0);
212: }
214: /*@
215: ISExpand - Computes the union of two index sets, by concatenating 2 lists and
216: removing duplicates.
218: Collective on IS
220: Input Parameter:
221: + is1 - first index set
222: - is2 - index values to be added
224: Output Parameters:
225: . isout - is1 + is2 The index set is2 is appended to is1 removing duplicates
227: Notes:
228: Negative values are removed from the lists. This requires O(imax-imin)
229: memory and O(imax-imin) work, where imin and imax are the bounds on the
230: indices in is1 and is2.
232: The IS's do not need to be sorted.
234: Level: intermediate
236: .seealso: ISDestroy(), ISView(), ISDifference(), ISSum()
239: @*/
240: PetscErrorCode ISExpand(IS is1,IS is2,IS *isout)
241: {
243: PetscInt i,n1,n2,imin,imax,nout,*iout;
244: const PetscInt *i1,*i2;
245: PetscBT mask;
246: MPI_Comm comm;
253: if (!is1 && !is2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Both arguments cannot be NULL");
254: if (!is1) {ISDuplicate(is2, isout);return(0);}
255: if (!is2) {ISDuplicate(is1, isout);return(0);}
256: ISGetIndices(is1,&i1);
257: ISGetLocalSize(is1,&n1);
258: ISGetIndices(is2,&i2);
259: ISGetLocalSize(is2,&n2);
261: /* Create a bit mask array to contain required values */
262: if (n1 || n2) {
263: imin = PETSC_MAX_INT;
264: imax = 0;
265: for (i=0; i<n1; i++) {
266: if (i1[i] < 0) continue;
267: imin = PetscMin(imin,i1[i]);
268: imax = PetscMax(imax,i1[i]);
269: }
270: for (i=0; i<n2; i++) {
271: if (i2[i] < 0) continue;
272: imin = PetscMin(imin,i2[i]);
273: imax = PetscMax(imax,i2[i]);
274: }
275: } else imin = imax = 0;
277: PetscMalloc1(n1+n2,&iout);
278: nout = 0;
279: PetscBTCreate(imax-imin,&mask);
280: /* Put the values from is1 */
281: for (i=0; i<n1; i++) {
282: if (i1[i] < 0) continue;
283: if (!PetscBTLookupSet(mask,i1[i] - imin)) iout[nout++] = i1[i];
284: }
285: ISRestoreIndices(is1,&i1);
286: /* Put the values from is2 */
287: for (i=0; i<n2; i++) {
288: if (i2[i] < 0) continue;
289: if (!PetscBTLookupSet(mask,i2[i] - imin)) iout[nout++] = i2[i];
290: }
291: ISRestoreIndices(is2,&i2);
293: /* create the new IS containing the sum */
294: PetscObjectGetComm((PetscObject)is1,&comm);
295: ISCreateGeneral(comm,nout,iout,PETSC_OWN_POINTER,isout);
297: PetscBTDestroy(&mask);
298: return(0);
299: }
301: /*@
302: ISIntersect - Computes the intersection of two index sets, by sorting and comparing.
304: Collective on IS
306: Input Parameter:
307: + is1 - first index set
308: - is2 - second index set
310: Output Parameters:
311: . isout - the sorted intersection of is1 and is2
313: Notes:
314: Negative values are removed from the lists. This requires O(min(is1,is2))
315: memory and O(max(is1,is2)log(max(is1,is2))) work
317: The IS's do not need to be sorted.
319: Level: intermediate
321: .seealso: ISDestroy(), ISView(), ISDifference(), ISSum(), ISExpand()
322: @*/
323: PetscErrorCode ISIntersect(IS is1,IS is2,IS *isout)
324: {
326: PetscInt i,n1,n2,nout,*iout;
327: const PetscInt *i1,*i2;
328: IS is1sorted = NULL, is2sorted = NULL;
329: PetscBool sorted, lsorted;
330: MPI_Comm comm;
337: PetscObjectGetComm((PetscObject)is1,&comm);
339: ISGetLocalSize(is1,&n1);
340: ISGetLocalSize(is2,&n2);
341: if (n1 < n2) {
342: IS tempis = is1;
343: PetscInt ntemp = n1;
345: is1 = is2;
346: is2 = tempis;
347: n1 = n2;
348: n2 = ntemp;
349: }
350: ISSorted(is1,&lsorted);
351: MPIU_Allreduce(&lsorted,&sorted,1,MPIU_BOOL,MPI_LAND,comm);
352: if (!sorted) {
353: ISDuplicate(is1,&is1sorted);
354: ISSort(is1sorted);
355: ISGetIndices(is1sorted,&i1);
356: } else {
357: is1sorted = is1;
358: PetscObjectReference((PetscObject)is1);
359: ISGetIndices(is1,&i1);
360: }
361: ISSorted(is2,&lsorted);
362: MPIU_Allreduce(&lsorted,&sorted,1,MPIU_BOOL,MPI_LAND,comm);
363: if (!sorted) {
364: ISDuplicate(is2,&is2sorted);
365: ISSort(is2sorted);
366: ISGetIndices(is2sorted,&i2);
367: } else {
368: is2sorted = is2;
369: PetscObjectReference((PetscObject)is2);
370: ISGetIndices(is2,&i2);
371: }
373: PetscMalloc1(n2,&iout);
375: for (i = 0, nout = 0; i < n2; i++) {
376: PetscInt key = i2[i];
377: PetscInt loc;
379: ISLocate(is1sorted,key,&loc);
380: if (loc >= 0) {
381: if (!nout || iout[nout-1] < key) {
382: iout[nout++] = key;
383: }
384: }
385: }
386: PetscRealloc(nout*sizeof(PetscInt),&iout);
388: /* create the new IS containing the sum */
389: ISCreateGeneral(comm,nout,iout,PETSC_OWN_POINTER,isout);
391: ISRestoreIndices(is2sorted,&i2);
392: ISDestroy(&is2sorted);
393: ISRestoreIndices(is1sorted,&i1);
394: ISDestroy(&is1sorted);
395: return(0);
396: }
398: PetscErrorCode ISIntersect_Caching_Internal(IS is1, IS is2, IS *isect)
399: {
403: *isect = NULL;
404: if (is2 && is1) {
405: char composeStr[33] = {0};
406: PetscObjectId is2id;
408: PetscObjectGetId((PetscObject)is2,&is2id);
409: PetscSNPrintf(composeStr,32,"ISIntersect_Caching_%x",is2id);
410: PetscObjectQuery((PetscObject) is1, composeStr, (PetscObject *) isect);
411: if (*isect == NULL) {
412: ISIntersect(is1, is2, isect);
413: PetscObjectCompose((PetscObject) is1, composeStr, (PetscObject) *isect);
414: } else {
415: PetscObjectReference((PetscObject) *isect);
416: }
417: }
418: return(0);
419: }
421: /*@
422: ISConcatenate - Forms a new IS by locally concatenating the indices from an IS list without reordering.
425: Collective.
427: Input Parameter:
428: + comm - communicator of the concatenated IS.
429: . len - size of islist array (nonnegative)
430: - islist - array of index sets
432: Output Parameters:
433: . isout - The concatenated index set; empty, if len == 0.
435: Notes:
436: The semantics of calling this on comm imply that the comms of the members if islist also contain this rank.
438: Level: intermediate
440: .seealso: ISDifference(), ISSum(), ISExpand()
443: @*/
444: PetscErrorCode ISConcatenate(MPI_Comm comm, PetscInt len, const IS islist[], IS *isout)
445: {
447: PetscInt i,n,N;
448: const PetscInt *iidx;
449: PetscInt *idx;
453: #if defined(PETSC_USE_DEBUG)
455: #endif
457: if (!len) {
458: ISCreateStride(comm, 0,0,0, isout);
459: return(0);
460: }
461: if (len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Negative array length: %D", len);
462: N = 0;
463: for (i = 0; i < len; ++i) {
464: if (islist[i]) {
465: ISGetLocalSize(islist[i], &n);
466: N += n;
467: }
468: }
469: PetscMalloc1(N, &idx);
470: N = 0;
471: for (i = 0; i < len; ++i) {
472: if (islist[i]) {
473: ISGetLocalSize(islist[i], &n);
474: ISGetIndices(islist[i], &iidx);
475: PetscArraycpy(idx+N,iidx, n);
476: ISRestoreIndices(islist[i], &iidx);
477: N += n;
478: }
479: }
480: ISCreateGeneral(comm, N, idx, PETSC_OWN_POINTER, isout);
481: return(0);
482: }
484: /*@
485: ISListToPair - convert an IS list to a pair of ISs of equal length defining an equivalent integer multimap.
486: Each IS on the input list is assigned an integer j so that all of the indices of that IS are
487: mapped to j.
490: Collective.
492: Input arguments:
493: + comm - MPI_Comm
494: . listlen - IS list length
495: - islist - IS list
497: Output arguments:
498: + xis - domain IS
499: - yis - range IS
501: Level: advanced
503: Notes:
504: The global integers assigned to the ISs of the local input list might not correspond to the
505: local numbers of the ISs on that list, but the two *orderings* are the same: the global
506: integers assigned to the ISs on the local list form a strictly increasing sequence.
508: The ISs on the input list can belong to subcommunicators of comm, and the subcommunicators
509: on the input IS list are assumed to be in a "deadlock-free" order.
511: Local lists of PetscObjects (or their subcommes) on a comm are "deadlock-free" if subcomm1
512: preceeds subcomm2 on any local list, then it preceeds subcomm2 on all ranks.
513: Equivalently, the local numbers of the subcomms on each local list are drawn from some global
514: numbering. This is ensured, for example, by ISPairToList().
516: .seealso ISPairToList()
517: @*/
518: PetscErrorCode ISListToPair(MPI_Comm comm, PetscInt listlen, IS islist[], IS *xis, IS *yis)
519: {
521: PetscInt ncolors, *colors,i, leni,len,*xinds, *yinds,k,j;
522: const PetscInt *indsi;
525: PetscMalloc1(listlen, &colors);
526: PetscObjectsListGetGlobalNumbering(comm, listlen, (PetscObject*)islist,&ncolors, colors);
527: len = 0;
528: for (i = 0; i < listlen; ++i) {
529: ISGetLocalSize(islist[i], &leni);
530: len += leni;
531: }
532: PetscMalloc1(len, &xinds);
533: PetscMalloc1(len, &yinds);
534: k = 0;
535: for (i = 0; i < listlen; ++i) {
536: ISGetLocalSize(islist[i], &leni);
537: ISGetIndices(islist[i],&indsi);
538: for (j = 0; j < leni; ++j) {
539: xinds[k] = indsi[j];
540: yinds[k] = colors[i];
541: ++k;
542: }
543: }
544: PetscFree(colors);
545: ISCreateGeneral(comm,len,xinds,PETSC_OWN_POINTER,xis);
546: ISCreateGeneral(comm,len,yinds,PETSC_OWN_POINTER,yis);
547: return(0);
548: }
551: /*@
552: ISPairToList - convert an IS pair encoding an integer map to a list of ISs.
553: Each IS on the output list contains the preimage for each index on the second input IS.
554: The ISs on the output list are constructed on the subcommunicators of the input IS pair.
555: Each subcommunicator corresponds to the preimage of some index j -- this subcomm contains
556: exactly the ranks that assign some indices i to j. This is essentially the inverse of
557: ISListToPair().
559: Collective on indis.
561: Input arguments:
562: + xis - domain IS
563: - yis - range IS
565: Output arguments:
566: + listlen - length of islist
567: - islist - list of ISs breaking up indis by color
569: Note:
570: + xis and yis must be of the same length and have congruent communicators.
571: - The resulting ISs have subcommunicators in a "deadlock-free" order (see ISListToPair()).
573: Level: advanced
575: .seealso ISListToPair()
576: @*/
577: PetscErrorCode ISPairToList(IS xis, IS yis, PetscInt *listlen, IS **islist)
578: {
580: IS indis = xis, coloris = yis;
581: PetscInt *inds, *colors, llen, ilen, lstart, lend, lcount,l;
582: PetscMPIInt rank, size, llow, lhigh, low, high,color,subsize;
583: const PetscInt *ccolors, *cinds;
584: MPI_Comm comm, subcomm;
592: PetscObjectGetComm((PetscObject)xis,&comm);
593: MPI_Comm_rank(comm, &rank);
594: MPI_Comm_rank(comm, &size);
595: /* Extract, copy and sort the local indices and colors on the color. */
596: ISGetLocalSize(coloris, &llen);
597: ISGetLocalSize(indis, &ilen);
598: if (llen != ilen) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Incompatible IS sizes: %D and %D", ilen, llen);
599: ISGetIndices(coloris, &ccolors);
600: ISGetIndices(indis, &cinds);
601: PetscMalloc2(ilen,&inds,llen,&colors);
602: PetscArraycpy(inds,cinds,ilen);
603: PetscArraycpy(colors,ccolors,llen);
604: PetscSortIntWithArray(llen, colors, inds);
605: /* Determine the global extent of colors. */
606: llow = 0; lhigh = -1;
607: lstart = 0; lcount = 0;
608: while (lstart < llen) {
609: lend = lstart+1;
610: while (lend < llen && colors[lend] == colors[lstart]) ++lend;
611: llow = PetscMin(llow,colors[lstart]);
612: lhigh = PetscMax(lhigh,colors[lstart]);
613: ++lcount;
614: }
615: MPIU_Allreduce(&llow,&low,1,MPI_INT,MPI_MIN,comm);
616: MPIU_Allreduce(&lhigh,&high,1,MPI_INT,MPI_MAX,comm);
617: *listlen = 0;
618: if (low <= high) {
619: if (lcount > 0) {
620: *listlen = lcount;
621: if (!*islist) {
622: PetscMalloc1(lcount, islist);
623: }
624: }
625: /*
626: Traverse all possible global colors, and participate in the subcommunicators
627: for the locally-supported colors.
628: */
629: lcount = 0;
630: lstart = 0; lend = 0;
631: for (l = low; l <= high; ++l) {
632: /*
633: Find the range of indices with the same color, which is not smaller than l.
634: Observe that, since colors is sorted, and is a subsequence of [low,high],
635: as soon as we find a new color, it is >= l.
636: */
637: if (lstart < llen) {
638: /* The start of the next locally-owned color is identified. Now look for the end. */
639: if (lstart == lend) {
640: lend = lstart+1;
641: while (lend < llen && colors[lend] == colors[lstart]) ++lend;
642: }
643: /* Now check whether the identified color segment matches l. */
644: if (colors[lstart] < l) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Locally owned color %D at location %D is < than the next global color %D", colors[lstart], lcount, l);
645: }
646: color = (PetscMPIInt)(colors[lstart] == l);
647: /* Check whether a proper subcommunicator exists. */
648: MPIU_Allreduce(&color,&subsize,1,MPI_INT,MPI_SUM,comm);
650: if (subsize == 1) subcomm = PETSC_COMM_SELF;
651: else if (subsize == size) subcomm = comm;
652: else {
653: /* a proper communicator is necessary, so we create it. */
654: MPI_Comm_split(comm, color, rank, &subcomm);
655: }
656: if (colors[lstart] == l) {
657: /* If we have l among the local colors, we create an IS to hold the corresponding indices. */
658: ISCreateGeneral(subcomm, lend-lstart,inds+lstart,PETSC_COPY_VALUES,*islist+lcount);
659: /* Position lstart at the beginning of the next local color. */
660: lstart = lend;
661: /* Increment the counter of the local colors split off into an IS. */
662: ++lcount;
663: }
664: if (subsize > 0 && subsize < size) {
665: /*
666: Irrespective of color, destroy the split off subcomm:
667: a subcomm used in the IS creation above is duplicated
668: into a proper PETSc comm.
669: */
670: MPI_Comm_free(&subcomm);
671: }
672: } /* for (l = low; l < high; ++l) */
673: } /* if (low <= high) */
674: PetscFree2(inds,colors);
675: return(0);
676: }
679: /*@
680: ISEmbed - embed IS a into IS b by finding the locations in b that have the same indices as in a.
681: If c is the IS of these locations, we have a = b*c, regarded as a composition of the
682: corresponding ISLocalToGlobalMaps.
684: Not collective.
686: Input arguments:
687: + a - IS to embed
688: . b - IS to embed into
689: - drop - flag indicating whether to drop a's indices that are not in b.
691: Output arguments:
692: . c - local embedding indices
694: Note:
695: If some of a's global indices are not among b's indices the embedding is impossible. The local indices of a
696: corresponding to these global indices are either mapped to -1 (if !drop) or are omitted (if drop). In the former
697: case the size of c is that same as that of a, in the latter case c's size may be smaller.
699: The resulting IS is sequential, since the index substition it encodes is purely local.
701: Level: advanced
703: .seealso ISLocalToGlobalMapping
704: @*/
705: PetscErrorCode ISEmbed(IS a, IS b, PetscBool drop, IS *c)
706: {
707: PetscErrorCode ierr;
708: ISLocalToGlobalMapping ltog;
709: ISGlobalToLocalMappingMode gtoltype = IS_GTOLM_DROP;
710: PetscInt alen, clen, *cindices, *cindices2;
711: const PetscInt *aindices;
717: ISLocalToGlobalMappingCreateIS(b, <og);
718: ISGetLocalSize(a, &alen);
719: ISGetIndices(a, &aindices);
720: PetscMalloc1(alen, &cindices);
721: if (!drop) gtoltype = IS_GTOLM_MASK;
722: ISGlobalToLocalMappingApply(ltog,gtoltype,alen,aindices,&clen,cindices);
723: ISLocalToGlobalMappingDestroy(<og);
724: if (clen != alen) {
725: cindices2 = cindices;
726: PetscMalloc1(clen, &cindices);
727: PetscArraycpy(cindices,cindices2,clen);
728: PetscFree(cindices2);
729: }
730: ISCreateGeneral(PETSC_COMM_SELF,clen,cindices,PETSC_OWN_POINTER,c);
731: return(0);
732: }
735: /*@
736: ISSortPermutation - calculate the permutation of the indices into a nondecreasing order.
738: Not collective.
740: Input arguments:
741: + f - IS to sort
742: - always - build the permutation even when f's indices are nondecreasing.
744: Output argument:
745: . h - permutation or NULL, if f is nondecreasing and always == PETSC_FALSE.
748: Note: Indices in f are unchanged. f[h[i]] is the i-th smallest f index.
749: If always == PETSC_FALSE, an extra check is peformed to see whether
750: the f indices are nondecreasing. h is built on PETSC_COMM_SELF, since
751: the permutation has a local meaning only.
753: Level: advanced
755: .seealso ISLocalToGlobalMapping, ISSort()
756: @*/
757: PetscErrorCode ISSortPermutation(IS f,PetscBool always,IS *h)
758: {
759: PetscErrorCode ierr;
760: const PetscInt *findices;
761: PetscInt fsize,*hindices,i;
762: PetscBool isincreasing;
767: ISGetLocalSize(f,&fsize);
768: ISGetIndices(f,&findices);
769: *h = NULL;
770: if (!always) {
771: isincreasing = PETSC_TRUE;
772: for (i = 1; i < fsize; ++i) {
773: if (findices[i] <= findices[i-1]) {
774: isincreasing = PETSC_FALSE;
775: break;
776: }
777: }
778: if (isincreasing) {
779: ISRestoreIndices(f,&findices);
780: return(0);
781: }
782: }
783: PetscMalloc1(fsize,&hindices);
784: for (i = 0; i < fsize; ++i) hindices[i] = i;
785: PetscSortIntWithPermutation(fsize,findices,hindices);
786: ISRestoreIndices(f,&findices);
787: ISCreateGeneral(PETSC_COMM_SELF,fsize,hindices,PETSC_OWN_POINTER,h);
788: ISSetInfo(*h,IS_PERMUTATION,IS_LOCAL,PETSC_FALSE,PETSC_TRUE);
789: return(0);
790: }