Actual source code: isdiff.c
petsc-3.14.6 2021-03-30
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: 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: {
32: PetscInt i,n1,n2,imin,imax,nout,*iout;
33: const PetscInt *i1,*i2;
34: PetscBT mask;
35: 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 Parameter:
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: Notes:
108: If n1 and n2 are the sizes of the sets, this takes O(n1+n2) time;
110: Both index sets need to be sorted on input.
112: Level: intermediate
114: .seealso: 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;
130: PetscObjectGetComm((PetscObject)(is1),&comm);
131: MPI_Comm_size(comm,&size);
132: if (size>1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently only for uni-processor IS");
134: ISSorted(is1,&f);
135: if (!f) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Arg 1 is not sorted");
136: ISSorted(is2,&f);
137: if (!f) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Arg 2 is not sorted");
139: ISGetLocalSize(is1,&n1);
140: ISGetLocalSize(is2,&n2);
141: if (!n2) {
142: ISDuplicate(is1,is3);
143: return(0);
144: }
145: ISGetIndices(is1,&i1);
146: ISGetIndices(is2,&i2);
148: p1 = 0; p2 = 0; n3 = 0;
149: do {
150: if (p1==n1) { /* cleanup for is2 */ n3 += n2-p2; break;
151: } else {
152: while (p2<n2 && i2[p2]<i1[p1]) {
153: n3++; p2++;
154: }
155: if (p2==n2) {
156: /* cleanup for is1 */
157: n3 += n1-p1; break;
158: } else {
159: if (i2[p2]==i1[p1]) { n3++; p1++; p2++; }
160: }
161: }
162: if (p2==n2) {
163: /* cleanup for is1 */
164: n3 += n1-p1; break;
165: } else {
166: while (p1<n1 && i1[p1]<i2[p2]) {
167: n3++; p1++;
168: }
169: if (p1==n1) {
170: /* clean up for is2 */
171: n3 += n2-p2; break;
172: } else {
173: if (i1[p1]==i2[p2]) { n3++; p1++; p2++; }
174: }
175: }
176: } while (p1<n1 || p2<n2);
178: if (n3==n1) { /* no new elements to be added */
179: ISRestoreIndices(is1,&i1);
180: ISRestoreIndices(is2,&i2);
181: ISDuplicate(is1,is3);
182: return(0);
183: }
184: PetscMalloc1(n3,&iout);
186: p1 = 0; p2 = 0; n3 = 0;
187: do {
188: if (p1==n1) { /* cleanup for is2 */
189: while (p2<n2) iout[n3++] = i2[p2++];
190: break;
191: } else {
192: while (p2<n2 && i2[p2]<i1[p1]) iout[n3++] = i2[p2++];
193: if (p2==n2) { /* cleanup for is1 */
194: while (p1<n1) iout[n3++] = i1[p1++];
195: break;
196: } else {
197: if (i2[p2]==i1[p1]) { iout[n3++] = i1[p1++]; p2++; }
198: }
199: }
200: if (p2==n2) { /* cleanup for is1 */
201: while (p1<n1) iout[n3++] = i1[p1++];
202: break;
203: } else {
204: while (p1<n1 && i1[p1]<i2[p2]) iout[n3++] = i1[p1++];
205: if (p1==n1) { /* clean up for is2 */
206: while (p2<n2) iout[n3++] = i2[p2++];
207: break;
208: } else {
209: if (i1[p1]==i2[p2]) { iout[n3++] = i1[p1++]; p2++; }
210: }
211: }
212: } while (p1<n1 || p2<n2);
214: ISRestoreIndices(is1,&i1);
215: ISRestoreIndices(is2,&i2);
216: ISCreateGeneral(comm,n3,iout,PETSC_OWN_POINTER,is3);
217: return(0);
218: }
220: /*@
221: ISExpand - Computes the union of two index sets, by concatenating 2 lists and
222: removing duplicates.
224: Collective on IS
226: Input Parameter:
227: + is1 - first index set
228: - is2 - index values to be added
230: Output Parameters:
231: . isout - is1 + is2 The index set is2 is appended to is1 removing duplicates
233: Notes:
234: Negative values are removed from the lists. This requires O(imax-imin)
235: memory and O(imax-imin) work, where imin and imax are the bounds on the
236: indices in is1 and is2.
238: The IS's do not need to be sorted.
240: Level: intermediate
242: .seealso: ISDestroy(), ISView(), ISDifference(), ISSum()
245: @*/
246: PetscErrorCode ISExpand(IS is1,IS is2,IS *isout)
247: {
249: PetscInt i,n1,n2,imin,imax,nout,*iout;
250: const PetscInt *i1,*i2;
251: PetscBT mask;
252: MPI_Comm comm;
259: if (!is1 && !is2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Both arguments cannot be NULL");
260: if (!is1) {ISDuplicate(is2, isout);return(0);}
261: if (!is2) {ISDuplicate(is1, isout);return(0);}
262: ISGetIndices(is1,&i1);
263: ISGetLocalSize(is1,&n1);
264: ISGetIndices(is2,&i2);
265: ISGetLocalSize(is2,&n2);
267: /* Create a bit mask array to contain required values */
268: if (n1 || n2) {
269: imin = PETSC_MAX_INT;
270: imax = 0;
271: for (i=0; i<n1; i++) {
272: if (i1[i] < 0) continue;
273: imin = PetscMin(imin,i1[i]);
274: imax = PetscMax(imax,i1[i]);
275: }
276: for (i=0; i<n2; i++) {
277: if (i2[i] < 0) continue;
278: imin = PetscMin(imin,i2[i]);
279: imax = PetscMax(imax,i2[i]);
280: }
281: } else imin = imax = 0;
283: PetscMalloc1(n1+n2,&iout);
284: nout = 0;
285: PetscBTCreate(imax-imin,&mask);
286: /* Put the values from is1 */
287: for (i=0; i<n1; i++) {
288: if (i1[i] < 0) continue;
289: if (!PetscBTLookupSet(mask,i1[i] - imin)) iout[nout++] = i1[i];
290: }
291: ISRestoreIndices(is1,&i1);
292: /* Put the values from is2 */
293: for (i=0; i<n2; i++) {
294: if (i2[i] < 0) continue;
295: if (!PetscBTLookupSet(mask,i2[i] - imin)) iout[nout++] = i2[i];
296: }
297: ISRestoreIndices(is2,&i2);
299: /* create the new IS containing the sum */
300: PetscObjectGetComm((PetscObject)is1,&comm);
301: ISCreateGeneral(comm,nout,iout,PETSC_OWN_POINTER,isout);
303: PetscBTDestroy(&mask);
304: return(0);
305: }
307: /*@
308: ISIntersect - Computes the intersection of two index sets, by sorting and comparing.
310: Collective on IS
312: Input Parameter:
313: + is1 - first index set
314: - is2 - second index set
316: Output Parameters:
317: . isout - the sorted intersection of is1 and is2
319: Notes:
320: Negative values are removed from the lists. This requires O(min(is1,is2))
321: memory and O(max(is1,is2)log(max(is1,is2))) work
323: The IS's do not need to be sorted.
325: Level: intermediate
327: .seealso: ISDestroy(), ISView(), ISDifference(), ISSum(), ISExpand()
328: @*/
329: PetscErrorCode ISIntersect(IS is1,IS is2,IS *isout)
330: {
332: PetscInt i,n1,n2,nout,*iout;
333: const PetscInt *i1,*i2;
334: IS is1sorted = NULL, is2sorted = NULL;
335: PetscBool sorted, lsorted;
336: MPI_Comm comm;
343: PetscObjectGetComm((PetscObject)is1,&comm);
345: ISGetLocalSize(is1,&n1);
346: ISGetLocalSize(is2,&n2);
347: if (n1 < n2) {
348: IS tempis = is1;
349: PetscInt ntemp = n1;
351: is1 = is2;
352: is2 = tempis;
353: n1 = n2;
354: n2 = ntemp;
355: }
356: ISSorted(is1,&lsorted);
357: MPIU_Allreduce(&lsorted,&sorted,1,MPIU_BOOL,MPI_LAND,comm);
358: if (!sorted) {
359: ISDuplicate(is1,&is1sorted);
360: ISSort(is1sorted);
361: ISGetIndices(is1sorted,&i1);
362: } else {
363: is1sorted = is1;
364: PetscObjectReference((PetscObject)is1);
365: ISGetIndices(is1,&i1);
366: }
367: ISSorted(is2,&lsorted);
368: MPIU_Allreduce(&lsorted,&sorted,1,MPIU_BOOL,MPI_LAND,comm);
369: if (!sorted) {
370: ISDuplicate(is2,&is2sorted);
371: ISSort(is2sorted);
372: ISGetIndices(is2sorted,&i2);
373: } else {
374: is2sorted = is2;
375: PetscObjectReference((PetscObject)is2);
376: ISGetIndices(is2,&i2);
377: }
379: PetscMalloc1(n2,&iout);
381: for (i = 0, nout = 0; i < n2; i++) {
382: PetscInt key = i2[i];
383: PetscInt loc;
385: ISLocate(is1sorted,key,&loc);
386: if (loc >= 0) {
387: if (!nout || iout[nout-1] < key) {
388: iout[nout++] = key;
389: }
390: }
391: }
392: PetscRealloc(nout*sizeof(PetscInt),&iout);
394: /* create the new IS containing the sum */
395: ISCreateGeneral(comm,nout,iout,PETSC_OWN_POINTER,isout);
397: ISRestoreIndices(is2sorted,&i2);
398: ISDestroy(&is2sorted);
399: ISRestoreIndices(is1sorted,&i1);
400: ISDestroy(&is1sorted);
401: return(0);
402: }
404: PetscErrorCode ISIntersect_Caching_Internal(IS is1, IS is2, IS *isect)
405: {
409: *isect = NULL;
410: if (is2 && is1) {
411: char composeStr[33] = {0};
412: PetscObjectId is2id;
414: PetscObjectGetId((PetscObject)is2,&is2id);
415: PetscSNPrintf(composeStr,32,"ISIntersect_Caching_%x",is2id);
416: PetscObjectQuery((PetscObject) is1, composeStr, (PetscObject *) isect);
417: if (*isect == NULL) {
418: ISIntersect(is1, is2, isect);
419: PetscObjectCompose((PetscObject) is1, composeStr, (PetscObject) *isect);
420: } else {
421: PetscObjectReference((PetscObject) *isect);
422: }
423: }
424: return(0);
425: }
427: /*@
428: ISConcatenate - Forms a new IS by locally concatenating the indices from an IS list without reordering.
431: Collective.
433: Input Parameter:
434: + comm - communicator of the concatenated IS.
435: . len - size of islist array (nonnegative)
436: - islist - array of index sets
438: Output Parameters:
439: . isout - The concatenated index set; empty, if len == 0.
441: Notes:
442: The semantics of calling this on comm imply that the comms of the members if islist also contain this rank.
444: Level: intermediate
446: .seealso: ISDifference(), ISSum(), ISExpand()
449: @*/
450: PetscErrorCode ISConcatenate(MPI_Comm comm, PetscInt len, const IS islist[], IS *isout)
451: {
453: PetscInt i,n,N;
454: const PetscInt *iidx;
455: PetscInt *idx;
459: if (PetscDefined(USE_DEBUG)) {
461: }
463: if (!len) {
464: ISCreateStride(comm, 0,0,0, isout);
465: return(0);
466: }
467: if (len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Negative array length: %D", len);
468: N = 0;
469: for (i = 0; i < len; ++i) {
470: if (islist[i]) {
471: ISGetLocalSize(islist[i], &n);
472: N += n;
473: }
474: }
475: PetscMalloc1(N, &idx);
476: N = 0;
477: for (i = 0; i < len; ++i) {
478: if (islist[i]) {
479: ISGetLocalSize(islist[i], &n);
480: ISGetIndices(islist[i], &iidx);
481: PetscArraycpy(idx+N,iidx, n);
482: ISRestoreIndices(islist[i], &iidx);
483: N += n;
484: }
485: }
486: ISCreateGeneral(comm, N, idx, PETSC_OWN_POINTER, isout);
487: return(0);
488: }
490: /*@
491: ISListToPair - convert an IS list to a pair of ISs of equal length defining an equivalent integer multimap.
492: Each IS on the input list is assigned an integer j so that all of the indices of that IS are
493: mapped to j.
496: Collective.
498: Input arguments:
499: + comm - MPI_Comm
500: . listlen - IS list length
501: - islist - IS list
503: Output arguments:
504: + xis - domain IS
505: - yis - range IS
507: Level: advanced
509: Notes:
510: The global integers assigned to the ISs of the local input list might not correspond to the
511: local numbers of the ISs on that list, but the two *orderings* are the same: the global
512: integers assigned to the ISs on the local list form a strictly increasing sequence.
514: The ISs on the input list can belong to subcommunicators of comm, and the subcommunicators
515: on the input IS list are assumed to be in a "deadlock-free" order.
517: Local lists of PetscObjects (or their subcommes) on a comm are "deadlock-free" if subcomm1
518: preceeds subcomm2 on any local list, then it preceeds subcomm2 on all ranks.
519: Equivalently, the local numbers of the subcomms on each local list are drawn from some global
520: numbering. This is ensured, for example, by ISPairToList().
522: .seealso ISPairToList()
523: @*/
524: PetscErrorCode ISListToPair(MPI_Comm comm, PetscInt listlen, IS islist[], IS *xis, IS *yis)
525: {
527: PetscInt ncolors, *colors,i, leni,len,*xinds, *yinds,k,j;
528: const PetscInt *indsi;
531: PetscMalloc1(listlen, &colors);
532: PetscObjectsListGetGlobalNumbering(comm, listlen, (PetscObject*)islist,&ncolors, colors);
533: len = 0;
534: for (i = 0; i < listlen; ++i) {
535: ISGetLocalSize(islist[i], &leni);
536: len += leni;
537: }
538: PetscMalloc1(len, &xinds);
539: PetscMalloc1(len, &yinds);
540: k = 0;
541: for (i = 0; i < listlen; ++i) {
542: ISGetLocalSize(islist[i], &leni);
543: ISGetIndices(islist[i],&indsi);
544: for (j = 0; j < leni; ++j) {
545: xinds[k] = indsi[j];
546: yinds[k] = colors[i];
547: ++k;
548: }
549: }
550: PetscFree(colors);
551: ISCreateGeneral(comm,len,xinds,PETSC_OWN_POINTER,xis);
552: ISCreateGeneral(comm,len,yinds,PETSC_OWN_POINTER,yis);
553: return(0);
554: }
557: /*@
558: ISPairToList - convert an IS pair encoding an integer map to a list of ISs.
559: Each IS on the output list contains the preimage for each index on the second input IS.
560: The ISs on the output list are constructed on the subcommunicators of the input IS pair.
561: Each subcommunicator corresponds to the preimage of some index j -- this subcomm contains
562: exactly the ranks that assign some indices i to j. This is essentially the inverse of
563: ISListToPair().
565: Collective on indis.
567: Input arguments:
568: + xis - domain IS
569: - yis - range IS
571: Output arguments:
572: + listlen - length of islist
573: - islist - list of ISs breaking up indis by color
575: Note:
576: + xis and yis must be of the same length and have congruent communicators.
577: - The resulting ISs have subcommunicators in a "deadlock-free" order (see ISListToPair()).
579: Level: advanced
581: .seealso ISListToPair()
582: @*/
583: PetscErrorCode ISPairToList(IS xis, IS yis, PetscInt *listlen, IS **islist)
584: {
586: IS indis = xis, coloris = yis;
587: PetscInt *inds, *colors, llen, ilen, lstart, lend, lcount,l;
588: PetscMPIInt rank, size, llow, lhigh, low, high,color,subsize;
589: const PetscInt *ccolors, *cinds;
590: MPI_Comm comm, subcomm;
598: PetscObjectGetComm((PetscObject)xis,&comm);
599: MPI_Comm_rank(comm, &rank);
600: MPI_Comm_rank(comm, &size);
601: /* Extract, copy and sort the local indices and colors on the color. */
602: ISGetLocalSize(coloris, &llen);
603: ISGetLocalSize(indis, &ilen);
604: if (llen != ilen) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Incompatible IS sizes: %D and %D", ilen, llen);
605: ISGetIndices(coloris, &ccolors);
606: ISGetIndices(indis, &cinds);
607: PetscMalloc2(ilen,&inds,llen,&colors);
608: PetscArraycpy(inds,cinds,ilen);
609: PetscArraycpy(colors,ccolors,llen);
610: PetscSortIntWithArray(llen, colors, inds);
611: /* Determine the global extent of colors. */
612: llow = 0; lhigh = -1;
613: lstart = 0; lcount = 0;
614: while (lstart < llen) {
615: lend = lstart+1;
616: while (lend < llen && colors[lend] == colors[lstart]) ++lend;
617: llow = PetscMin(llow,colors[lstart]);
618: lhigh = PetscMax(lhigh,colors[lstart]);
619: ++lcount;
620: }
621: MPIU_Allreduce(&llow,&low,1,MPI_INT,MPI_MIN,comm);
622: MPIU_Allreduce(&lhigh,&high,1,MPI_INT,MPI_MAX,comm);
623: *listlen = 0;
624: if (low <= high) {
625: if (lcount > 0) {
626: *listlen = lcount;
627: if (!*islist) {
628: PetscMalloc1(lcount, islist);
629: }
630: }
631: /*
632: Traverse all possible global colors, and participate in the subcommunicators
633: for the locally-supported colors.
634: */
635: lcount = 0;
636: lstart = 0; lend = 0;
637: for (l = low; l <= high; ++l) {
638: /*
639: Find the range of indices with the same color, which is not smaller than l.
640: Observe that, since colors is sorted, and is a subsequence of [low,high],
641: as soon as we find a new color, it is >= l.
642: */
643: if (lstart < llen) {
644: /* The start of the next locally-owned color is identified. Now look for the end. */
645: if (lstart == lend) {
646: lend = lstart+1;
647: while (lend < llen && colors[lend] == colors[lstart]) ++lend;
648: }
649: /* Now check whether the identified color segment matches l. */
650: 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);
651: }
652: color = (PetscMPIInt)(colors[lstart] == l);
653: /* Check whether a proper subcommunicator exists. */
654: MPIU_Allreduce(&color,&subsize,1,MPI_INT,MPI_SUM,comm);
656: if (subsize == 1) subcomm = PETSC_COMM_SELF;
657: else if (subsize == size) subcomm = comm;
658: else {
659: /* a proper communicator is necessary, so we create it. */
660: MPI_Comm_split(comm, color, rank, &subcomm);
661: }
662: if (colors[lstart] == l) {
663: /* If we have l among the local colors, we create an IS to hold the corresponding indices. */
664: ISCreateGeneral(subcomm, lend-lstart,inds+lstart,PETSC_COPY_VALUES,*islist+lcount);
665: /* Position lstart at the beginning of the next local color. */
666: lstart = lend;
667: /* Increment the counter of the local colors split off into an IS. */
668: ++lcount;
669: }
670: if (subsize > 0 && subsize < size) {
671: /*
672: Irrespective of color, destroy the split off subcomm:
673: a subcomm used in the IS creation above is duplicated
674: into a proper PETSc comm.
675: */
676: MPI_Comm_free(&subcomm);
677: }
678: } /* for (l = low; l < high; ++l) */
679: } /* if (low <= high) */
680: PetscFree2(inds,colors);
681: return(0);
682: }
685: /*@
686: ISEmbed - embed IS a into IS b by finding the locations in b that have the same indices as in a.
687: If c is the IS of these locations, we have a = b*c, regarded as a composition of the
688: corresponding ISLocalToGlobalMaps.
690: Not collective.
692: Input arguments:
693: + a - IS to embed
694: . b - IS to embed into
695: - drop - flag indicating whether to drop a's indices that are not in b.
697: Output arguments:
698: . c - local embedding indices
700: Note:
701: If some of a's global indices are not among b's indices the embedding is impossible. The local indices of a
702: corresponding to these global indices are either mapped to -1 (if !drop) or are omitted (if drop). In the former
703: case the size of c is that same as that of a, in the latter case c's size may be smaller.
705: The resulting IS is sequential, since the index substition it encodes is purely local.
707: Level: advanced
709: .seealso ISLocalToGlobalMapping
710: @*/
711: PetscErrorCode ISEmbed(IS a, IS b, PetscBool drop, IS *c)
712: {
713: PetscErrorCode ierr;
714: ISLocalToGlobalMapping ltog;
715: ISGlobalToLocalMappingMode gtoltype = IS_GTOLM_DROP;
716: PetscInt alen, clen, *cindices, *cindices2;
717: const PetscInt *aindices;
723: ISLocalToGlobalMappingCreateIS(b, <og);
724: ISGetLocalSize(a, &alen);
725: ISGetIndices(a, &aindices);
726: PetscMalloc1(alen, &cindices);
727: if (!drop) gtoltype = IS_GTOLM_MASK;
728: ISGlobalToLocalMappingApply(ltog,gtoltype,alen,aindices,&clen,cindices);
729: ISLocalToGlobalMappingDestroy(<og);
730: if (clen != alen) {
731: cindices2 = cindices;
732: PetscMalloc1(clen, &cindices);
733: PetscArraycpy(cindices,cindices2,clen);
734: PetscFree(cindices2);
735: }
736: ISCreateGeneral(PETSC_COMM_SELF,clen,cindices,PETSC_OWN_POINTER,c);
737: return(0);
738: }
741: /*@
742: ISSortPermutation - calculate the permutation of the indices into a nondecreasing order.
744: Not collective.
746: Input arguments:
747: + f - IS to sort
748: - always - build the permutation even when f's indices are nondecreasing.
750: Output argument:
751: . h - permutation or NULL, if f is nondecreasing and always == PETSC_FALSE.
754: Note: Indices in f are unchanged. f[h[i]] is the i-th smallest f index.
755: If always == PETSC_FALSE, an extra check is peformed to see whether
756: the f indices are nondecreasing. h is built on PETSC_COMM_SELF, since
757: the permutation has a local meaning only.
759: Level: advanced
761: .seealso ISLocalToGlobalMapping, ISSort()
762: @*/
763: PetscErrorCode ISSortPermutation(IS f,PetscBool always,IS *h)
764: {
765: PetscErrorCode ierr;
766: const PetscInt *findices;
767: PetscInt fsize,*hindices,i;
768: PetscBool isincreasing;
773: ISGetLocalSize(f,&fsize);
774: ISGetIndices(f,&findices);
775: *h = NULL;
776: if (!always) {
777: isincreasing = PETSC_TRUE;
778: for (i = 1; i < fsize; ++i) {
779: if (findices[i] <= findices[i-1]) {
780: isincreasing = PETSC_FALSE;
781: break;
782: }
783: }
784: if (isincreasing) {
785: ISRestoreIndices(f,&findices);
786: return(0);
787: }
788: }
789: PetscMalloc1(fsize,&hindices);
790: for (i = 0; i < fsize; ++i) hindices[i] = i;
791: PetscSortIntWithPermutation(fsize,findices,hindices);
792: ISRestoreIndices(f,&findices);
793: ISCreateGeneral(PETSC_COMM_SELF,fsize,hindices,PETSC_OWN_POINTER,h);
794: ISSetInfo(*h,IS_PERMUTATION,IS_LOCAL,PETSC_FALSE,PETSC_TRUE);
795: return(0);
796: }