Actual source code: isdiff.c
petsc-3.10.5 2019-03-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) return(0);
141: ISGetIndices(is1,&i1);
142: ISGetIndices(is2,&i2);
144: p1 = 0; p2 = 0; n3 = 0;
145: do {
146: if (p1==n1) { /* cleanup for is2 */ n3 += n2-p2; break;
147: } else {
148: while (p2<n2 && i2[p2]<i1[p1]) {
149: n3++; p2++;
150: }
151: if (p2==n2) {
152: /* cleanup for is1 */
153: n3 += n1-p1; break;
154: } else {
155: if (i2[p2]==i1[p1]) { n3++; p1++; p2++; }
156: }
157: }
158: if (p2==n2) {
159: /* cleanup for is1 */
160: n3 += n1-p1; break;
161: } else {
162: while (p1<n1 && i1[p1]<i2[p2]) {
163: n3++; p1++;
164: }
165: if (p1==n1) {
166: /* clean up for is2 */
167: n3 += n2-p2; break;
168: } else {
169: if (i1[p1]==i2[p2]) { n3++; p1++; p2++; }
170: }
171: }
172: } while (p1<n1 || p2<n2);
174: if (n3==n1) { /* no new elements to be added */
175: ISRestoreIndices(is1,&i1);
176: ISRestoreIndices(is2,&i2);
177: ISDuplicate(is1,is3);
178: return(0);
179: }
180: PetscMalloc1(n3,&iout);
182: p1 = 0; p2 = 0; n3 = 0;
183: do {
184: if (p1==n1) { /* cleanup for is2 */
185: while (p2<n2) iout[n3++] = i2[p2++];
186: break;
187: } else {
188: while (p2<n2 && i2[p2]<i1[p1]) iout[n3++] = i2[p2++];
189: if (p2==n2) { /* cleanup for is1 */
190: while (p1<n1) iout[n3++] = i1[p1++];
191: break;
192: } else {
193: if (i2[p2]==i1[p1]) { iout[n3++] = i1[p1++]; p2++; }
194: }
195: }
196: if (p2==n2) { /* cleanup for is1 */
197: while (p1<n1) iout[n3++] = i1[p1++];
198: break;
199: } else {
200: while (p1<n1 && i1[p1]<i2[p2]) iout[n3++] = i1[p1++];
201: if (p1==n1) { /* clean up for is2 */
202: while (p2<n2) iout[n3++] = i2[p2++];
203: break;
204: } else {
205: if (i1[p1]==i2[p2]) { iout[n3++] = i1[p1++]; p2++; }
206: }
207: }
208: } while (p1<n1 || p2<n2);
210: ISRestoreIndices(is1,&i1);
211: ISRestoreIndices(is2,&i2);
212: ISCreateGeneral(comm,n3,iout,PETSC_OWN_POINTER,is3);
213: return(0);
214: }
216: /*@
217: ISExpand - Computes the union of two index sets, by concatenating 2 lists and
218: removing duplicates.
220: Collective on IS
222: Input Parameter:
223: + is1 - first index set
224: - is2 - index values to be added
226: Output Parameters:
227: . isout - is1 + is2 The index set is2 is appended to is1 removing duplicates
229: Notes:
230: Negative values are removed from the lists. This requires O(imax-imin)
231: memory and O(imax-imin) work, where imin and imax are the bounds on the
232: indices in is1 and is2.
234: The IS's do not need to be sorted.
236: Level: intermediate
238: .seealso: ISDestroy(), ISView(), ISDifference(), ISSum()
240: Concepts: index sets^difference
241: Concepts: IS^difference
243: @*/
244: PetscErrorCode ISExpand(IS is1,IS is2,IS *isout)
245: {
247: PetscInt i,n1,n2,imin,imax,nout,*iout;
248: const PetscInt *i1,*i2;
249: PetscBT mask;
250: MPI_Comm comm;
257: if (!is1 && !is2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Both arguments cannot be NULL");
258: if (!is1) {ISDuplicate(is2, isout);return(0);}
259: if (!is2) {ISDuplicate(is1, isout);return(0);}
260: ISGetIndices(is1,&i1);
261: ISGetLocalSize(is1,&n1);
262: ISGetIndices(is2,&i2);
263: ISGetLocalSize(is2,&n2);
265: /* Create a bit mask array to contain required values */
266: if (n1 || n2) {
267: imin = PETSC_MAX_INT;
268: imax = 0;
269: for (i=0; i<n1; i++) {
270: if (i1[i] < 0) continue;
271: imin = PetscMin(imin,i1[i]);
272: imax = PetscMax(imax,i1[i]);
273: }
274: for (i=0; i<n2; i++) {
275: if (i2[i] < 0) continue;
276: imin = PetscMin(imin,i2[i]);
277: imax = PetscMax(imax,i2[i]);
278: }
279: } else imin = imax = 0;
281: PetscMalloc1(n1+n2,&iout);
282: nout = 0;
283: PetscBTCreate(imax-imin,&mask);
284: /* Put the values from is1 */
285: for (i=0; i<n1; i++) {
286: if (i1[i] < 0) continue;
287: if (!PetscBTLookupSet(mask,i1[i] - imin)) iout[nout++] = i1[i];
288: }
289: ISRestoreIndices(is1,&i1);
290: /* Put the values from is2 */
291: for (i=0; i<n2; i++) {
292: if (i2[i] < 0) continue;
293: if (!PetscBTLookupSet(mask,i2[i] - imin)) iout[nout++] = i2[i];
294: }
295: ISRestoreIndices(is2,&i2);
297: /* create the new IS containing the sum */
298: PetscObjectGetComm((PetscObject)is1,&comm);
299: ISCreateGeneral(comm,nout,iout,PETSC_OWN_POINTER,isout);
301: PetscBTDestroy(&mask);
302: return(0);
303: }
305: /*@
306: ISIntersect - Computes the union of two index sets, by concatenating 2 lists, sorting,
307: and finding duplicates.
309: Collective on IS
311: Input Parameter:
312: + is1 - first index set
313: - is2 - second index set
315: Output Parameters:
316: . isout - the sorted intersection of is1 and is2
318: Notes:
319: Negative values are removed from the lists. This requires O(min(is1,is2))
320: memory and O(max(is1,is2)log(max(is1,is2))) work
322: The IS's do not need to be sorted.
324: Level: intermediate
326: .seealso: ISDestroy(), ISView(), ISDifference(), ISSum(), ISExpand()
328: Concepts: index sets^intersection
329: Concepts: IS^intersection
331: @*/
332: PetscErrorCode ISIntersect(IS is1,IS is2,IS *isout)
333: {
335: PetscInt i,n1,n2,nout,*iout;
336: const PetscInt *i1,*i2;
337: IS is1sorted = NULL, is2sorted = NULL;
338: PetscBool sorted;
339: MPI_Comm comm;
346: ISGetLocalSize(is1,&n1);
347: ISGetLocalSize(is2,&n2);
348: if (n1 < n2) {
349: IS tempis = is1;
350: int ntemp = n1;
352: is1 = is2;
353: is2 = tempis;
354: n1 = n2;
355: n2 = ntemp;
356: }
357: ISSorted(is1,&sorted);
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,&sorted);
368: if (!sorted) {
369: ISDuplicate(is2,&is2sorted);
370: ISSort(is2sorted);
371: ISGetIndices(is2sorted,&i2);
372: } else {
373: is2sorted = is2;
374: PetscObjectReference((PetscObject)is2);
375: ISGetIndices(is2,&i2);
376: }
378: PetscMalloc1(n2,&iout);
380: for (i = 0, nout = 0; i < n2; i++) {
381: PetscInt key = i2[i];
382: PetscInt loc;
384: ISLocate(is1sorted,key,&loc);
385: if (loc >= 0) {
386: if (!nout || iout[nout-1] < key) {
387: iout[nout++] = key;
388: }
389: }
390: }
391: PetscRealloc(sizeof (PetscInt) * (size_t) nout,&iout);
393: /* create the new IS containing the sum */
394: PetscObjectGetComm((PetscObject)is1,&comm);
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 on comm.
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()
448: Concepts: index sets^concatenation
449: Concepts: IS^concatenation
451: @*/
452: PetscErrorCode ISConcatenate(MPI_Comm comm, PetscInt len, const IS islist[], IS *isout)
453: {
455: PetscInt i,n,N;
456: const PetscInt *iidx;
457: PetscInt *idx;
461: #if defined(PETSC_USE_DEBUG)
463: #endif
465: if (!len) {
466: ISCreateStride(comm, 0,0,0, isout);
467: return(0);
468: }
469: if (len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Negative array length: %D", len);
470: N = 0;
471: for (i = 0; i < len; ++i) {
472: if (islist[i]) {
473: ISGetLocalSize(islist[i], &n);
474: N += n;
475: }
476: }
477: PetscMalloc1(N, &idx);
478: N = 0;
479: for (i = 0; i < len; ++i) {
480: if (islist[i]) {
481: ISGetLocalSize(islist[i], &n);
482: ISGetIndices(islist[i], &iidx);
483: PetscMemcpy(idx+N,iidx, sizeof(PetscInt)*n);
484: ISRestoreIndices(islist[i], &iidx);
485: N += n;
486: }
487: }
488: ISCreateGeneral(comm, N, idx, PETSC_OWN_POINTER, isout);
489: return(0);
490: }
492: /*@
493: ISListToPair - convert an IS list to a pair of ISs of equal length defining an equivalent integer multimap.
494: Each IS on the input list is assigned an integer j so that all of the indices of that IS are
495: mapped to j.
498: Collective on comm.
500: Input arguments:
501: + comm - MPI_Comm
502: . listlen - IS list length
503: - islist - IS list
505: Output arguments:
506: + xis - domain IS
507: - yis - range IS
509: Level: advanced
511: Notes:
512: The global integers assigned to the ISs of the local input list might not correspond to the
513: local numbers of the ISs on that list, but the two *orderings* are the same: the global
514: integers assigned to the ISs on the local list form a strictly increasing sequence.
516: The ISs on the input list can belong to subcommunicators of comm, and the subcommunicators
517: on the input IS list are assumed to be in a "deadlock-free" order.
519: Local lists of PetscObjects (or their subcommes) on a comm are "deadlock-free" if subcomm1
520: preceeds subcomm2 on any local list, then it preceeds subcomm2 on all ranks.
521: Equivalently, the local numbers of the subcomms on each local list are drawn from some global
522: numbering. This is ensured, for example, by ISPairToList().
524: .seealso ISPairToList()
525: @*/
526: PetscErrorCode ISListToPair(MPI_Comm comm, PetscInt listlen, IS islist[], IS *xis, IS *yis)
527: {
529: PetscInt ncolors, *colors,i, leni,len,*xinds, *yinds,k,j;
530: const PetscInt *indsi;
533: PetscMalloc1(listlen, &colors);
534: PetscObjectsListGetGlobalNumbering(comm, listlen, (PetscObject*)islist,&ncolors, colors);
535: len = 0;
536: for (i = 0; i < listlen; ++i) {
537: ISGetLocalSize(islist[i], &leni);
538: len += leni;
539: }
540: PetscMalloc1(len, &xinds);
541: PetscMalloc1(len, &yinds);
542: k = 0;
543: for (i = 0; i < listlen; ++i) {
544: ISGetLocalSize(islist[i], &leni);
545: ISGetIndices(islist[i],&indsi);
546: for (j = 0; j < leni; ++j) {
547: xinds[k] = indsi[j];
548: yinds[k] = colors[i];
549: ++k;
550: }
551: }
552: PetscFree(colors);
553: ISCreateGeneral(comm,len,xinds,PETSC_OWN_POINTER,xis);
554: ISCreateGeneral(comm,len,yinds,PETSC_OWN_POINTER,yis);
555: return(0);
556: }
559: /*@
560: ISPairToList - convert an IS pair encoding an integer map to a list of ISs.
561: Each IS on the output list contains the preimage for each index on the second input IS.
562: The ISs on the output list are constructed on the subcommunicators of the input IS pair.
563: Each subcommunicator corresponds to the preimage of some index j -- this subcomm contains
564: exactly the ranks that assign some indices i to j. This is essentially the inverse of
565: ISListToPair().
567: Collective on indis.
569: Input arguments:
570: + xis - domain IS
571: - yis - range IS
573: Output arguments:
574: + listlen - length of islist
575: - islist - list of ISs breaking up indis by color
577: Note:
578: + xis and yis must be of the same length and have congruent communicators.
579: - The resulting ISs have subcommunicators in a "deadlock-free" order (see ISListToPair()).
581: Level: advanced
583: .seealso ISListToPair()
584: @*/
585: PetscErrorCode ISPairToList(IS xis, IS yis, PetscInt *listlen, IS **islist)
586: {
588: IS indis = xis, coloris = yis;
589: PetscInt *inds, *colors, llen, ilen, lstart, lend, lcount,l;
590: PetscMPIInt rank, size, llow, lhigh, low, high,color,subsize;
591: const PetscInt *ccolors, *cinds;
592: MPI_Comm comm, subcomm;
600: PetscObjectGetComm((PetscObject)xis,&comm);
601: MPI_Comm_rank(comm, &rank);
602: MPI_Comm_rank(comm, &size);
603: /* Extract, copy and sort the local indices and colors on the color. */
604: ISGetLocalSize(coloris, &llen);
605: ISGetLocalSize(indis, &ilen);
606: if (llen != ilen) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Incompatible IS sizes: %D and %D", ilen, llen);
607: ISGetIndices(coloris, &ccolors);
608: ISGetIndices(indis, &cinds);
609: PetscMalloc2(ilen,&inds,llen,&colors);
610: PetscMemcpy(inds,cinds,ilen*sizeof(PetscInt));
611: PetscMemcpy(colors,ccolors,llen*sizeof(PetscInt));
612: PetscSortIntWithArray(llen, colors, inds);
613: /* Determine the global extent of colors. */
614: llow = 0; lhigh = -1;
615: lstart = 0; lcount = 0;
616: while (lstart < llen) {
617: lend = lstart+1;
618: while (lend < llen && colors[lend] == colors[lstart]) ++lend;
619: llow = PetscMin(llow,colors[lstart]);
620: lhigh = PetscMax(lhigh,colors[lstart]);
621: ++lcount;
622: }
623: MPIU_Allreduce(&llow,&low,1,MPI_INT,MPI_MIN,comm);
624: MPIU_Allreduce(&lhigh,&high,1,MPI_INT,MPI_MAX,comm);
625: *listlen = 0;
626: if (low <= high) {
627: if (lcount > 0) {
628: *listlen = lcount;
629: if (!*islist) {
630: PetscMalloc1(lcount, islist);
631: }
632: }
633: /*
634: Traverse all possible global colors, and participate in the subcommunicators
635: for the locally-supported colors.
636: */
637: lcount = 0;
638: lstart = 0; lend = 0;
639: for (l = low; l <= high; ++l) {
640: /*
641: Find the range of indices with the same color, which is not smaller than l.
642: Observe that, since colors is sorted, and is a subsequence of [low,high],
643: as soon as we find a new color, it is >= l.
644: */
645: if (lstart < llen) {
646: /* The start of the next locally-owned color is identified. Now look for the end. */
647: if (lstart == lend) {
648: lend = lstart+1;
649: while (lend < llen && colors[lend] == colors[lstart]) ++lend;
650: }
651: /* Now check whether the identified color segment matches l. */
652: 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);
653: }
654: color = (PetscMPIInt)(colors[lstart] == l);
655: /* Check whether a proper subcommunicator exists. */
656: MPIU_Allreduce(&color,&subsize,1,MPI_INT,MPI_SUM,comm);
658: if (subsize == 1) subcomm = PETSC_COMM_SELF;
659: else if (subsize == size) subcomm = comm;
660: else {
661: /* a proper communicator is necessary, so we create it. */
662: MPI_Comm_split(comm, color, rank, &subcomm);
663: }
664: if (colors[lstart] == l) {
665: /* If we have l among the local colors, we create an IS to hold the corresponding indices. */
666: ISCreateGeneral(subcomm, lend-lstart,inds+lstart,PETSC_COPY_VALUES,*islist+lcount);
667: /* Position lstart at the beginning of the next local color. */
668: lstart = lend;
669: /* Increment the counter of the local colors split off into an IS. */
670: ++lcount;
671: }
672: if (subsize > 0 && subsize < size) {
673: /*
674: Irrespective of color, destroy the split off subcomm:
675: a subcomm used in the IS creation above is duplicated
676: into a proper PETSc comm.
677: */
678: MPI_Comm_free(&subcomm);
679: }
680: } /* for (l = low; l < high; ++l) */
681: } /* if (low <= high) */
682: PetscFree2(inds,colors);
683: return(0);
684: }
687: /*@
688: ISEmbed - embed IS a into IS b by finding the locations in b that have the same indices as in a.
689: If c is the IS of these locations, we have a = b*c, regarded as a composition of the
690: corresponding ISLocalToGlobalMaps.
692: Not collective.
694: Input arguments:
695: + a - IS to embed
696: . b - IS to embed into
697: - drop - flag indicating whether to drop a's indices that are not in b.
699: Output arguments:
700: . c - local embedding indices
702: Note:
703: If some of a's global indices are not among b's indices the embedding is impossible. The local indices of a
704: corresponding to these global indices are either mapped to -1 (if !drop) or are omitted (if drop). In the former
705: case the size of c is that same as that of a, in the latter case c's size may be smaller.
707: The resulting IS is sequential, since the index substition it encodes is purely local.
709: Level: advanced
711: .seealso ISLocalToGlobalMapping
712: @*/
713: PetscErrorCode ISEmbed(IS a, IS b, PetscBool drop, IS *c)
714: {
715: PetscErrorCode ierr;
716: ISLocalToGlobalMapping ltog;
717: ISGlobalToLocalMappingMode gtoltype = IS_GTOLM_DROP;
718: PetscInt alen, clen, *cindices, *cindices2;
719: const PetscInt *aindices;
725: ISLocalToGlobalMappingCreateIS(b, <og);
726: ISGetLocalSize(a, &alen);
727: ISGetIndices(a, &aindices);
728: PetscMalloc1(alen, &cindices);
729: if (!drop) gtoltype = IS_GTOLM_MASK;
730: ISGlobalToLocalMappingApply(ltog,gtoltype,alen,aindices,&clen,cindices);
731: ISLocalToGlobalMappingDestroy(<og);
732: if (clen != alen) {
733: cindices2 = cindices;
734: PetscMalloc1(clen, &cindices);
735: PetscMemcpy(cindices,cindices2,clen*sizeof(PetscInt));
736: PetscFree(cindices2);
737: }
738: ISCreateGeneral(PETSC_COMM_SELF,clen,cindices,PETSC_OWN_POINTER,c);
739: return(0);
740: }
743: /*@
744: ISSortPermutation - calculate the permutation of the indices into a nondecreasing order.
746: Not collective.
748: Input arguments:
749: + f - IS to sort
750: - always - build the permutation even when f's indices are nondecreasin.
752: Output argument:
753: . h - permutation or NULL, if f is nondecreasing and always == PETSC_TRUE.
756: Note: Indices in f are unchanged. f[h[i]] is the i-th smallest f index.
757: If always == PETSC_FALSE, an extra check is peformed to see whether
758: the f indices are nondecreasing. h is built on PETSC_COMM_SELF, since
759: the permutation has a local meaning only.
761: Level: advanced
763: .seealso ISLocalToGlobalMapping, ISSort(), PetscIntSortWithPermutation()
764: @*/
765: PetscErrorCode ISSortPermutation(IS f,PetscBool always,IS *h)
766: {
767: PetscErrorCode ierr;
768: const PetscInt *findices;
769: PetscInt fsize,*hindices,i;
770: PetscBool isincreasing;
775: ISGetLocalSize(f,&fsize);
776: ISGetIndices(f,&findices);
777: *h = NULL;
778: if (!always) {
779: isincreasing = PETSC_TRUE;
780: for (i = 1; i < fsize; ++i) {
781: if (findices[i] <= findices[i-1]) {
782: isincreasing = PETSC_FALSE;
783: break;
784: }
785: }
786: if (isincreasing) {
787: ISRestoreIndices(f,&findices);
788: return(0);
789: }
790: }
791: PetscMalloc1(fsize,&hindices);
792: for (i = 0; i < fsize; ++i) hindices[i] = i;
793: PetscSortIntWithPermutation(fsize,findices,hindices);
794: ISRestoreIndices(f,&findices);
795: ISCreateGeneral(PETSC_COMM_SELF,fsize,hindices,PETSC_OWN_POINTER,h);
796: return(0);
797: }