Actual source code: isdiff.c
petsc-3.11.4 2019-09-28
2: #include <petsc/private/isimpl.h>
3: #include <petscbt.h>
5: /*@
6: ISDifference - Computes the difference between two index sets.
8: Collective on IS
10: Input Parameter:
11: + is1 - first index, to have items removed from it
12: - is2 - index values to be removed
14: Output Parameters:
15: . isout - is1 - is2
17: Notes:
18: Negative values are removed from the lists. is2 may have values
19: that are not in is1. This requires O(imax-imin) memory and O(imax-imin)
20: work, where imin and imax are the bounds on the indices in is1.
22: Level: intermediate
24: Concepts: index sets^difference
25: Concepts: IS^difference
27: .seealso: ISDestroy(), ISView(), ISSum(), ISExpand()
29: @*/
30: PetscErrorCode ISDifference(IS is1,IS is2,IS *isout)
31: {
33: PetscInt i,n1,n2,imin,imax,nout,*iout;
34: const PetscInt *i1,*i2;
35: PetscBT mask;
36: MPI_Comm comm;
43: ISGetIndices(is1,&i1);
44: ISGetLocalSize(is1,&n1);
46: /* Create a bit mask array to contain required values */
47: if (n1) {
48: imin = PETSC_MAX_INT;
49: imax = 0;
50: for (i=0; i<n1; i++) {
51: if (i1[i] < 0) continue;
52: imin = PetscMin(imin,i1[i]);
53: imax = PetscMax(imax,i1[i]);
54: }
55: } else imin = imax = 0;
57: PetscBTCreate(imax-imin,&mask);
58: /* Put the values from is1 */
59: for (i=0; i<n1; i++) {
60: if (i1[i] < 0) continue;
61: PetscBTSet(mask,i1[i] - imin);
62: }
63: ISRestoreIndices(is1,&i1);
64: /* Remove the values from is2 */
65: ISGetIndices(is2,&i2);
66: ISGetLocalSize(is2,&n2);
67: for (i=0; i<n2; i++) {
68: if (i2[i] < imin || i2[i] > imax) continue;
69: PetscBTClear(mask,i2[i] - imin);
70: }
71: ISRestoreIndices(is2,&i2);
73: /* Count the number in the difference */
74: nout = 0;
75: for (i=0; i<imax-imin+1; i++) {
76: if (PetscBTLookup(mask,i)) nout++;
77: }
79: /* create the new IS containing the difference */
80: PetscMalloc1(nout,&iout);
81: nout = 0;
82: for (i=0; i<imax-imin+1; i++) {
83: if (PetscBTLookup(mask,i)) iout[nout++] = i + imin;
84: }
85: PetscObjectGetComm((PetscObject)is1,&comm);
86: ISCreateGeneral(comm,nout,iout,PETSC_OWN_POINTER,isout);
88: PetscBTDestroy(&mask);
89: return(0);
90: }
92: /*@
93: ISSum - Computes the sum (union) of two index sets.
95: Only sequential version (at the moment)
97: Input Parameter:
98: + is1 - index set to be extended
99: - is2 - index values to be added
101: Output Parameter:
102: . is3 - the sum; this can not be is1 or is2
104: Notes:
105: If n1 and n2 are the sizes of the sets, this takes O(n1+n2) time;
107: Both index sets need to be sorted on input.
109: Level: intermediate
111: .seealso: ISDestroy(), ISView(), ISDifference(), ISExpand()
113: Concepts: index sets^union
114: Concepts: IS^union
116: @*/
117: PetscErrorCode ISSum(IS is1,IS is2,IS *is3)
118: {
119: MPI_Comm comm;
120: PetscBool f;
121: PetscMPIInt size;
122: const PetscInt *i1,*i2;
123: PetscInt n1,n2,n3, p1,p2, *iout;
129: PetscObjectGetComm((PetscObject)(is1),&comm);
130: MPI_Comm_size(comm,&size);
131: if (size>1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently only for uni-processor IS");
133: ISSorted(is1,&f);
134: if (!f) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Arg 1 is not sorted");
135: ISSorted(is2,&f);
136: if (!f) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Arg 2 is not sorted");
138: ISGetLocalSize(is1,&n1);
139: ISGetLocalSize(is2,&n2);
140: if (!n2) {
141: ISDuplicate(is1,is3);
142: return(0);
143: }
144: ISGetIndices(is1,&i1);
145: ISGetIndices(is2,&i2);
147: p1 = 0; p2 = 0; n3 = 0;
148: do {
149: if (p1==n1) { /* cleanup for is2 */ n3 += n2-p2; break;
150: } else {
151: while (p2<n2 && i2[p2]<i1[p1]) {
152: n3++; p2++;
153: }
154: if (p2==n2) {
155: /* cleanup for is1 */
156: n3 += n1-p1; break;
157: } else {
158: if (i2[p2]==i1[p1]) { n3++; p1++; p2++; }
159: }
160: }
161: if (p2==n2) {
162: /* cleanup for is1 */
163: n3 += n1-p1; break;
164: } else {
165: while (p1<n1 && i1[p1]<i2[p2]) {
166: n3++; p1++;
167: }
168: if (p1==n1) {
169: /* clean up for is2 */
170: n3 += n2-p2; break;
171: } else {
172: if (i1[p1]==i2[p2]) { n3++; p1++; p2++; }
173: }
174: }
175: } while (p1<n1 || p2<n2);
177: if (n3==n1) { /* no new elements to be added */
178: ISRestoreIndices(is1,&i1);
179: ISRestoreIndices(is2,&i2);
180: ISDuplicate(is1,is3);
181: return(0);
182: }
183: PetscMalloc1(n3,&iout);
185: p1 = 0; p2 = 0; n3 = 0;
186: do {
187: if (p1==n1) { /* cleanup for is2 */
188: while (p2<n2) iout[n3++] = i2[p2++];
189: break;
190: } else {
191: while (p2<n2 && i2[p2]<i1[p1]) iout[n3++] = i2[p2++];
192: if (p2==n2) { /* cleanup for is1 */
193: while (p1<n1) iout[n3++] = i1[p1++];
194: break;
195: } else {
196: if (i2[p2]==i1[p1]) { iout[n3++] = i1[p1++]; p2++; }
197: }
198: }
199: if (p2==n2) { /* cleanup for is1 */
200: while (p1<n1) iout[n3++] = i1[p1++];
201: break;
202: } else {
203: while (p1<n1 && i1[p1]<i2[p2]) iout[n3++] = i1[p1++];
204: if (p1==n1) { /* clean up for is2 */
205: while (p2<n2) iout[n3++] = i2[p2++];
206: break;
207: } else {
208: if (i1[p1]==i2[p2]) { iout[n3++] = i1[p1++]; p2++; }
209: }
210: }
211: } while (p1<n1 || p2<n2);
213: ISRestoreIndices(is1,&i1);
214: ISRestoreIndices(is2,&i2);
215: ISCreateGeneral(comm,n3,iout,PETSC_OWN_POINTER,is3);
216: return(0);
217: }
219: /*@
220: ISExpand - Computes the union of two index sets, by concatenating 2 lists and
221: removing duplicates.
223: Collective on IS
225: Input Parameter:
226: + is1 - first index set
227: - is2 - index values to be added
229: Output Parameters:
230: . isout - is1 + is2 The index set is2 is appended to is1 removing duplicates
232: Notes:
233: Negative values are removed from the lists. This requires O(imax-imin)
234: memory and O(imax-imin) work, where imin and imax are the bounds on the
235: indices in is1 and is2.
237: The IS's do not need to be sorted.
239: Level: intermediate
241: .seealso: ISDestroy(), ISView(), ISDifference(), ISSum()
243: Concepts: index sets^difference
244: Concepts: IS^difference
246: @*/
247: PetscErrorCode ISExpand(IS is1,IS is2,IS *isout)
248: {
250: PetscInt i,n1,n2,imin,imax,nout,*iout;
251: const PetscInt *i1,*i2;
252: PetscBT mask;
253: MPI_Comm comm;
260: if (!is1 && !is2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Both arguments cannot be NULL");
261: if (!is1) {ISDuplicate(is2, isout);return(0);}
262: if (!is2) {ISDuplicate(is1, isout);return(0);}
263: ISGetIndices(is1,&i1);
264: ISGetLocalSize(is1,&n1);
265: ISGetIndices(is2,&i2);
266: ISGetLocalSize(is2,&n2);
268: /* Create a bit mask array to contain required values */
269: if (n1 || n2) {
270: imin = PETSC_MAX_INT;
271: imax = 0;
272: for (i=0; i<n1; i++) {
273: if (i1[i] < 0) continue;
274: imin = PetscMin(imin,i1[i]);
275: imax = PetscMax(imax,i1[i]);
276: }
277: for (i=0; i<n2; i++) {
278: if (i2[i] < 0) continue;
279: imin = PetscMin(imin,i2[i]);
280: imax = PetscMax(imax,i2[i]);
281: }
282: } else imin = imax = 0;
284: PetscMalloc1(n1+n2,&iout);
285: nout = 0;
286: PetscBTCreate(imax-imin,&mask);
287: /* Put the values from is1 */
288: for (i=0; i<n1; i++) {
289: if (i1[i] < 0) continue;
290: if (!PetscBTLookupSet(mask,i1[i] - imin)) iout[nout++] = i1[i];
291: }
292: ISRestoreIndices(is1,&i1);
293: /* Put the values from is2 */
294: for (i=0; i<n2; i++) {
295: if (i2[i] < 0) continue;
296: if (!PetscBTLookupSet(mask,i2[i] - imin)) iout[nout++] = i2[i];
297: }
298: ISRestoreIndices(is2,&i2);
300: /* create the new IS containing the sum */
301: PetscObjectGetComm((PetscObject)is1,&comm);
302: ISCreateGeneral(comm,nout,iout,PETSC_OWN_POINTER,isout);
304: PetscBTDestroy(&mask);
305: return(0);
306: }
308: /*@
309: ISIntersect - Computes the union of two index sets, by concatenating 2 lists, sorting,
310: and finding duplicates.
312: Collective on IS
314: Input Parameter:
315: + is1 - first index set
316: - is2 - second index set
318: Output Parameters:
319: . isout - the sorted intersection of is1 and is2
321: Notes:
322: Negative values are removed from the lists. This requires O(min(is1,is2))
323: memory and O(max(is1,is2)log(max(is1,is2))) work
325: The IS's do not need to be sorted.
327: Level: intermediate
329: .seealso: ISDestroy(), ISView(), ISDifference(), ISSum(), ISExpand()
331: Concepts: index sets^intersection
332: Concepts: IS^intersection
334: @*/
335: PetscErrorCode ISIntersect(IS is1,IS is2,IS *isout)
336: {
338: PetscInt i,n1,n2,nout,*iout;
339: const PetscInt *i1,*i2;
340: IS is1sorted = NULL, is2sorted = NULL;
341: PetscBool sorted;
342: MPI_Comm comm;
349: ISGetLocalSize(is1,&n1);
350: ISGetLocalSize(is2,&n2);
351: if (n1 < n2) {
352: IS tempis = is1;
353: int ntemp = n1;
355: is1 = is2;
356: is2 = tempis;
357: n1 = n2;
358: n2 = ntemp;
359: }
360: ISSorted(is1,&sorted);
361: if (!sorted) {
362: ISDuplicate(is1,&is1sorted);
363: ISSort(is1sorted);
364: ISGetIndices(is1sorted,&i1);
365: } else {
366: is1sorted = is1;
367: PetscObjectReference((PetscObject)is1);
368: ISGetIndices(is1,&i1);
369: }
370: ISSorted(is2,&sorted);
371: if (!sorted) {
372: ISDuplicate(is2,&is2sorted);
373: ISSort(is2sorted);
374: ISGetIndices(is2sorted,&i2);
375: } else {
376: is2sorted = is2;
377: PetscObjectReference((PetscObject)is2);
378: ISGetIndices(is2,&i2);
379: }
381: PetscMalloc1(n2,&iout);
383: for (i = 0, nout = 0; i < n2; i++) {
384: PetscInt key = i2[i];
385: PetscInt loc;
387: ISLocate(is1sorted,key,&loc);
388: if (loc >= 0) {
389: if (!nout || iout[nout-1] < key) {
390: iout[nout++] = key;
391: }
392: }
393: }
394: PetscRealloc(sizeof (PetscInt) * (size_t) nout,&iout);
396: /* create the new IS containing the sum */
397: PetscObjectGetComm((PetscObject)is1,&comm);
398: ISCreateGeneral(comm,nout,iout,PETSC_OWN_POINTER,isout);
400: ISRestoreIndices(is2sorted,&i2);
401: ISDestroy(&is2sorted);
402: ISRestoreIndices(is1sorted,&i1);
403: ISDestroy(&is1sorted);
404: return(0);
405: }
407: PetscErrorCode ISIntersect_Caching_Internal(IS is1, IS is2, IS *isect)
408: {
412: *isect = NULL;
413: if (is2 && is1) {
414: char composeStr[33] = {0};
415: PetscObjectId is2id;
417: PetscObjectGetId((PetscObject)is2,&is2id);
418: PetscSNPrintf(composeStr,32,"ISIntersect_Caching_%x",is2id);
419: PetscObjectQuery((PetscObject) is1, composeStr, (PetscObject *) isect);
420: if (*isect == NULL) {
421: ISIntersect(is1, is2, isect);
422: PetscObjectCompose((PetscObject) is1, composeStr, (PetscObject) *isect);
423: } else {
424: PetscObjectReference((PetscObject) *isect);
425: }
426: }
427: return(0);
428: }
430: /*@
431: ISConcatenate - Forms a new IS by locally concatenating the indices from an IS list without reordering.
434: Collective on comm.
436: Input Parameter:
437: + comm - communicator of the concatenated IS.
438: . len - size of islist array (nonnegative)
439: - islist - array of index sets
441: Output Parameters:
442: . isout - The concatenated index set; empty, if len == 0.
444: Notes:
445: The semantics of calling this on comm imply that the comms of the members if islist also contain this rank.
447: Level: intermediate
449: .seealso: ISDifference(), ISSum(), ISExpand()
451: Concepts: index sets^concatenation
452: Concepts: IS^concatenation
454: @*/
455: PetscErrorCode ISConcatenate(MPI_Comm comm, PetscInt len, const IS islist[], IS *isout)
456: {
458: PetscInt i,n,N;
459: const PetscInt *iidx;
460: PetscInt *idx;
464: #if defined(PETSC_USE_DEBUG)
466: #endif
468: if (!len) {
469: ISCreateStride(comm, 0,0,0, isout);
470: return(0);
471: }
472: if (len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Negative array length: %D", len);
473: N = 0;
474: for (i = 0; i < len; ++i) {
475: if (islist[i]) {
476: ISGetLocalSize(islist[i], &n);
477: N += n;
478: }
479: }
480: PetscMalloc1(N, &idx);
481: N = 0;
482: for (i = 0; i < len; ++i) {
483: if (islist[i]) {
484: ISGetLocalSize(islist[i], &n);
485: ISGetIndices(islist[i], &iidx);
486: PetscMemcpy(idx+N,iidx, sizeof(PetscInt)*n);
487: ISRestoreIndices(islist[i], &iidx);
488: N += n;
489: }
490: }
491: ISCreateGeneral(comm, N, idx, PETSC_OWN_POINTER, isout);
492: return(0);
493: }
495: /*@
496: ISListToPair - convert an IS list to a pair of ISs of equal length defining an equivalent integer multimap.
497: Each IS on the input list is assigned an integer j so that all of the indices of that IS are
498: mapped to j.
501: Collective on comm.
503: Input arguments:
504: + comm - MPI_Comm
505: . listlen - IS list length
506: - islist - IS list
508: Output arguments:
509: + xis - domain IS
510: - yis - range IS
512: Level: advanced
514: Notes:
515: The global integers assigned to the ISs of the local input list might not correspond to the
516: local numbers of the ISs on that list, but the two *orderings* are the same: the global
517: integers assigned to the ISs on the local list form a strictly increasing sequence.
519: The ISs on the input list can belong to subcommunicators of comm, and the subcommunicators
520: on the input IS list are assumed to be in a "deadlock-free" order.
522: Local lists of PetscObjects (or their subcommes) on a comm are "deadlock-free" if subcomm1
523: preceeds subcomm2 on any local list, then it preceeds subcomm2 on all ranks.
524: Equivalently, the local numbers of the subcomms on each local list are drawn from some global
525: numbering. This is ensured, for example, by ISPairToList().
527: .seealso ISPairToList()
528: @*/
529: PetscErrorCode ISListToPair(MPI_Comm comm, PetscInt listlen, IS islist[], IS *xis, IS *yis)
530: {
532: PetscInt ncolors, *colors,i, leni,len,*xinds, *yinds,k,j;
533: const PetscInt *indsi;
536: PetscMalloc1(listlen, &colors);
537: PetscObjectsListGetGlobalNumbering(comm, listlen, (PetscObject*)islist,&ncolors, colors);
538: len = 0;
539: for (i = 0; i < listlen; ++i) {
540: ISGetLocalSize(islist[i], &leni);
541: len += leni;
542: }
543: PetscMalloc1(len, &xinds);
544: PetscMalloc1(len, &yinds);
545: k = 0;
546: for (i = 0; i < listlen; ++i) {
547: ISGetLocalSize(islist[i], &leni);
548: ISGetIndices(islist[i],&indsi);
549: for (j = 0; j < leni; ++j) {
550: xinds[k] = indsi[j];
551: yinds[k] = colors[i];
552: ++k;
553: }
554: }
555: PetscFree(colors);
556: ISCreateGeneral(comm,len,xinds,PETSC_OWN_POINTER,xis);
557: ISCreateGeneral(comm,len,yinds,PETSC_OWN_POINTER,yis);
558: return(0);
559: }
562: /*@
563: ISPairToList - convert an IS pair encoding an integer map to a list of ISs.
564: Each IS on the output list contains the preimage for each index on the second input IS.
565: The ISs on the output list are constructed on the subcommunicators of the input IS pair.
566: Each subcommunicator corresponds to the preimage of some index j -- this subcomm contains
567: exactly the ranks that assign some indices i to j. This is essentially the inverse of
568: ISListToPair().
570: Collective on indis.
572: Input arguments:
573: + xis - domain IS
574: - yis - range IS
576: Output arguments:
577: + listlen - length of islist
578: - islist - list of ISs breaking up indis by color
580: Note:
581: + xis and yis must be of the same length and have congruent communicators.
582: - The resulting ISs have subcommunicators in a "deadlock-free" order (see ISListToPair()).
584: Level: advanced
586: .seealso ISListToPair()
587: @*/
588: PetscErrorCode ISPairToList(IS xis, IS yis, PetscInt *listlen, IS **islist)
589: {
591: IS indis = xis, coloris = yis;
592: PetscInt *inds, *colors, llen, ilen, lstart, lend, lcount,l;
593: PetscMPIInt rank, size, llow, lhigh, low, high,color,subsize;
594: const PetscInt *ccolors, *cinds;
595: MPI_Comm comm, subcomm;
603: PetscObjectGetComm((PetscObject)xis,&comm);
604: MPI_Comm_rank(comm, &rank);
605: MPI_Comm_rank(comm, &size);
606: /* Extract, copy and sort the local indices and colors on the color. */
607: ISGetLocalSize(coloris, &llen);
608: ISGetLocalSize(indis, &ilen);
609: if (llen != ilen) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Incompatible IS sizes: %D and %D", ilen, llen);
610: ISGetIndices(coloris, &ccolors);
611: ISGetIndices(indis, &cinds);
612: PetscMalloc2(ilen,&inds,llen,&colors);
613: PetscMemcpy(inds,cinds,ilen*sizeof(PetscInt));
614: PetscMemcpy(colors,ccolors,llen*sizeof(PetscInt));
615: PetscSortIntWithArray(llen, colors, inds);
616: /* Determine the global extent of colors. */
617: llow = 0; lhigh = -1;
618: lstart = 0; lcount = 0;
619: while (lstart < llen) {
620: lend = lstart+1;
621: while (lend < llen && colors[lend] == colors[lstart]) ++lend;
622: llow = PetscMin(llow,colors[lstart]);
623: lhigh = PetscMax(lhigh,colors[lstart]);
624: ++lcount;
625: }
626: MPIU_Allreduce(&llow,&low,1,MPI_INT,MPI_MIN,comm);
627: MPIU_Allreduce(&lhigh,&high,1,MPI_INT,MPI_MAX,comm);
628: *listlen = 0;
629: if (low <= high) {
630: if (lcount > 0) {
631: *listlen = lcount;
632: if (!*islist) {
633: PetscMalloc1(lcount, islist);
634: }
635: }
636: /*
637: Traverse all possible global colors, and participate in the subcommunicators
638: for the locally-supported colors.
639: */
640: lcount = 0;
641: lstart = 0; lend = 0;
642: for (l = low; l <= high; ++l) {
643: /*
644: Find the range of indices with the same color, which is not smaller than l.
645: Observe that, since colors is sorted, and is a subsequence of [low,high],
646: as soon as we find a new color, it is >= l.
647: */
648: if (lstart < llen) {
649: /* The start of the next locally-owned color is identified. Now look for the end. */
650: if (lstart == lend) {
651: lend = lstart+1;
652: while (lend < llen && colors[lend] == colors[lstart]) ++lend;
653: }
654: /* Now check whether the identified color segment matches l. */
655: 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);
656: }
657: color = (PetscMPIInt)(colors[lstart] == l);
658: /* Check whether a proper subcommunicator exists. */
659: MPIU_Allreduce(&color,&subsize,1,MPI_INT,MPI_SUM,comm);
661: if (subsize == 1) subcomm = PETSC_COMM_SELF;
662: else if (subsize == size) subcomm = comm;
663: else {
664: /* a proper communicator is necessary, so we create it. */
665: MPI_Comm_split(comm, color, rank, &subcomm);
666: }
667: if (colors[lstart] == l) {
668: /* If we have l among the local colors, we create an IS to hold the corresponding indices. */
669: ISCreateGeneral(subcomm, lend-lstart,inds+lstart,PETSC_COPY_VALUES,*islist+lcount);
670: /* Position lstart at the beginning of the next local color. */
671: lstart = lend;
672: /* Increment the counter of the local colors split off into an IS. */
673: ++lcount;
674: }
675: if (subsize > 0 && subsize < size) {
676: /*
677: Irrespective of color, destroy the split off subcomm:
678: a subcomm used in the IS creation above is duplicated
679: into a proper PETSc comm.
680: */
681: MPI_Comm_free(&subcomm);
682: }
683: } /* for (l = low; l < high; ++l) */
684: } /* if (low <= high) */
685: PetscFree2(inds,colors);
686: return(0);
687: }
690: /*@
691: ISEmbed - embed IS a into IS b by finding the locations in b that have the same indices as in a.
692: If c is the IS of these locations, we have a = b*c, regarded as a composition of the
693: corresponding ISLocalToGlobalMaps.
695: Not collective.
697: Input arguments:
698: + a - IS to embed
699: . b - IS to embed into
700: - drop - flag indicating whether to drop a's indices that are not in b.
702: Output arguments:
703: . c - local embedding indices
705: Note:
706: If some of a's global indices are not among b's indices the embedding is impossible. The local indices of a
707: corresponding to these global indices are either mapped to -1 (if !drop) or are omitted (if drop). In the former
708: case the size of c is that same as that of a, in the latter case c's size may be smaller.
710: The resulting IS is sequential, since the index substition it encodes is purely local.
712: Level: advanced
714: .seealso ISLocalToGlobalMapping
715: @*/
716: PetscErrorCode ISEmbed(IS a, IS b, PetscBool drop, IS *c)
717: {
718: PetscErrorCode ierr;
719: ISLocalToGlobalMapping ltog;
720: ISGlobalToLocalMappingMode gtoltype = IS_GTOLM_DROP;
721: PetscInt alen, clen, *cindices, *cindices2;
722: const PetscInt *aindices;
728: ISLocalToGlobalMappingCreateIS(b, <og);
729: ISGetLocalSize(a, &alen);
730: ISGetIndices(a, &aindices);
731: PetscMalloc1(alen, &cindices);
732: if (!drop) gtoltype = IS_GTOLM_MASK;
733: ISGlobalToLocalMappingApply(ltog,gtoltype,alen,aindices,&clen,cindices);
734: ISLocalToGlobalMappingDestroy(<og);
735: if (clen != alen) {
736: cindices2 = cindices;
737: PetscMalloc1(clen, &cindices);
738: PetscMemcpy(cindices,cindices2,clen*sizeof(PetscInt));
739: PetscFree(cindices2);
740: }
741: ISCreateGeneral(PETSC_COMM_SELF,clen,cindices,PETSC_OWN_POINTER,c);
742: return(0);
743: }
746: /*@
747: ISSortPermutation - calculate the permutation of the indices into a nondecreasing order.
749: Not collective.
751: Input arguments:
752: + f - IS to sort
753: - always - build the permutation even when f's indices are nondecreasin.
755: Output argument:
756: . h - permutation or NULL, if f is nondecreasing and always == PETSC_TRUE.
759: Note: Indices in f are unchanged. f[h[i]] is the i-th smallest f index.
760: If always == PETSC_FALSE, an extra check is peformed to see whether
761: the f indices are nondecreasing. h is built on PETSC_COMM_SELF, since
762: the permutation has a local meaning only.
764: Level: advanced
766: .seealso ISLocalToGlobalMapping, ISSort(), PetscIntSortWithPermutation()
767: @*/
768: PetscErrorCode ISSortPermutation(IS f,PetscBool always,IS *h)
769: {
770: PetscErrorCode ierr;
771: const PetscInt *findices;
772: PetscInt fsize,*hindices,i;
773: PetscBool isincreasing;
778: ISGetLocalSize(f,&fsize);
779: ISGetIndices(f,&findices);
780: *h = NULL;
781: if (!always) {
782: isincreasing = PETSC_TRUE;
783: for (i = 1; i < fsize; ++i) {
784: if (findices[i] <= findices[i-1]) {
785: isincreasing = PETSC_FALSE;
786: break;
787: }
788: }
789: if (isincreasing) {
790: ISRestoreIndices(f,&findices);
791: return(0);
792: }
793: }
794: PetscMalloc1(fsize,&hindices);
795: for (i = 0; i < fsize; ++i) hindices[i] = i;
796: PetscSortIntWithPermutation(fsize,findices,hindices);
797: ISRestoreIndices(f,&findices);
798: ISCreateGeneral(PETSC_COMM_SELF,fsize,hindices,PETSC_OWN_POINTER,h);
799: return(0);
800: }