Actual source code: isdiff.c
petsc-3.7.5 2017-01-01
2: #include <petsc/private/isimpl.h> /*I "petscis.h" I*/
3: #include <petscbt.h>
7: /*@
8: ISDifference - Computes the difference between two index sets.
10: Collective on IS
12: Input Parameter:
13: + is1 - first index, to have items removed from it
14: - is2 - index values to be removed
16: Output Parameters:
17: . isout - is1 - is2
19: Notes:
20: Negative values are removed from the lists. is2 may have values
21: that are not in is1. This requires O(imax-imin) memory and O(imax-imin)
22: work, where imin and imax are the bounds on the indices in is1.
24: Level: intermediate
26: Concepts: index sets^difference
27: Concepts: IS^difference
29: .seealso: ISDestroy(), ISView(), ISSum(), ISExpand()
31: @*/
32: PetscErrorCode ISDifference(IS is1,IS is2,IS *isout)
33: {
35: PetscInt i,n1,n2,imin,imax,nout,*iout;
36: const PetscInt *i1,*i2;
37: PetscBT mask;
38: MPI_Comm comm;
45: ISGetIndices(is1,&i1);
46: ISGetLocalSize(is1,&n1);
48: /* Create a bit mask array to contain required values */
49: if (n1) {
50: imin = PETSC_MAX_INT;
51: imax = 0;
52: for (i=0; i<n1; i++) {
53: if (i1[i] < 0) continue;
54: imin = PetscMin(imin,i1[i]);
55: imax = PetscMax(imax,i1[i]);
56: }
57: } else imin = imax = 0;
59: PetscBTCreate(imax-imin,&mask);
60: /* Put the values from is1 */
61: for (i=0; i<n1; i++) {
62: if (i1[i] < 0) continue;
63: PetscBTSet(mask,i1[i] - imin);
64: }
65: ISRestoreIndices(is1,&i1);
66: /* Remove the values from is2 */
67: ISGetIndices(is2,&i2);
68: ISGetLocalSize(is2,&n2);
69: for (i=0; i<n2; i++) {
70: if (i2[i] < imin || i2[i] > imax) continue;
71: PetscBTClear(mask,i2[i] - imin);
72: }
73: ISRestoreIndices(is2,&i2);
75: /* Count the number in the difference */
76: nout = 0;
77: for (i=0; i<imax-imin+1; i++) {
78: if (PetscBTLookup(mask,i)) nout++;
79: }
81: /* create the new IS containing the difference */
82: PetscMalloc1(nout,&iout);
83: nout = 0;
84: for (i=0; i<imax-imin+1; i++) {
85: if (PetscBTLookup(mask,i)) iout[nout++] = i + imin;
86: }
87: PetscObjectGetComm((PetscObject)is1,&comm);
88: ISCreateGeneral(comm,nout,iout,PETSC_OWN_POINTER,isout);
90: PetscBTDestroy(&mask);
91: return(0);
92: }
96: /*@
97: ISSum - Computes the sum (union) of two index sets.
99: Only sequential version (at the moment)
101: Input Parameter:
102: + is1 - index set to be extended
103: - is2 - index values to be added
105: Output Parameter:
106: . is3 - the sum; this can not be is1 or is2
108: Notes:
109: If n1 and n2 are the sizes of the sets, this takes O(n1+n2) time;
111: Both index sets need to be sorted on input.
113: Level: intermediate
115: .seealso: ISDestroy(), ISView(), ISDifference(), ISExpand()
117: Concepts: index sets^union
118: Concepts: IS^union
120: @*/
121: PetscErrorCode ISSum(IS is1,IS is2,IS *is3)
122: {
123: MPI_Comm comm;
124: PetscBool f;
125: PetscMPIInt size;
126: const PetscInt *i1,*i2;
127: PetscInt n1,n2,n3, p1,p2, *iout;
133: PetscObjectGetComm((PetscObject)(is1),&comm);
134: MPI_Comm_size(comm,&size);
135: if (size>1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently only for uni-processor IS");
137: ISSorted(is1,&f);
138: if (!f) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Arg 1 is not sorted");
139: ISSorted(is2,&f);
140: if (!f) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Arg 2 is not sorted");
142: ISGetLocalSize(is1,&n1);
143: ISGetLocalSize(is2,&n2);
144: if (!n2) return(0);
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: }
222: /*@
223: ISExpand - Computes the union of two index sets, by concatenating 2 lists and
224: removing duplicates.
226: Collective on IS
228: Input Parameter:
229: + is1 - first index set
230: - is2 - index values to be added
232: Output Parameters:
233: . isout - is1 + is2 The index set is2 is appended to is1 removing duplicates
235: Notes:
236: Negative values are removed from the lists. This requires O(imax-imin)
237: memory and O(imax-imin) work, where imin and imax are the bounds on the
238: indices in is1 and is2.
240: The IS's do not need to be sorted.
242: Level: intermediate
244: .seealso: ISDestroy(), ISView(), ISDifference(), ISSum()
246: Concepts: index sets^difference
247: Concepts: IS^difference
249: @*/
250: PetscErrorCode ISExpand(IS is1,IS is2,IS *isout)
251: {
253: PetscInt i,n1,n2,imin,imax,nout,*iout;
254: const PetscInt *i1,*i2;
255: PetscBT mask;
256: MPI_Comm comm;
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: }
310: /*@
311: ISConcatenate - Forms a new IS by locally concatenating the indices from an IS list without reordering.
314: Collective on comm.
316: Input Parameter:
317: + comm - communicator of the concatenated IS.
318: . len - size of islist array (nonnegative)
319: - islist - array of index sets
323: Output Parameters:
324: . isout - The concatenated index set; empty, if len == 0.
326: Notes: The semantics of calling this on comm imply that the comms of the members if islist also contain this rank.
328: Level: intermediate
330: .seealso: ISDifference(), ISSum(), ISExpand()
332: Concepts: index sets^concatenation
333: Concepts: IS^concatenation
335: @*/
336: PetscErrorCode ISConcatenate(MPI_Comm comm, PetscInt len, const IS islist[], IS *isout)
337: {
339: PetscInt i,n,N;
340: const PetscInt *iidx;
341: PetscInt *idx;
345: #if defined(PETSC_USE_DEBUG)
347: #endif
349: if (!len) {
350: ISCreateStride(comm, 0,0,0, isout);
351: return(0);
352: }
353: if (len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Negative array length: %D", len);
354: N = 0;
355: for (i = 0; i < len; ++i) {
356: ISGetLocalSize(islist[i], &n);
357: N += n;
358: }
359: PetscMalloc(sizeof(PetscInt)*N, &idx);
360: N = 0;
361: for (i = 0; i < len; ++i) {
362: ISGetLocalSize(islist[i], &n);
363: ISGetIndices(islist[i], &iidx);
364: PetscMemcpy(idx+N,iidx, sizeof(PetscInt)*n);
365: ISRestoreIndices(islist[i], &iidx);
366: N += n;
367: }
368: ISCreateGeneral(comm, N, idx, PETSC_OWN_POINTER, isout);
369: return(0);
370: }
372: /*@
373: ISListToPair - convert an IS list to a pair of ISs of equal length defining an equivalent integer multimap.
374: Each IS on the input list is assigned an integer j so that all of the indices of that IS are
375: mapped to j.
378: Collective on comm.
380: Input arguments:
381: + comm - MPI_Comm
382: . listlen - IS list length
383: - islist - IS list
385: Output arguments:
386: + xis - domain IS
387: - yis - range IS
389: Level: advanced
391: Notes:
392: The global integers assigned to the ISs of the local input list might not correspond to the
393: local numbers of the ISs on that list, but the two *orderings* are the same: the global
394: integers assigned to the ISs on the local list form a strictly increasing sequence.
396: The ISs on the input list can belong to subcommunicators of comm, and the subcommunicators
397: on the input IS list are assumed to be in a "deadlock-free" order.
399: Local lists of PetscObjects (or their subcommes) on a comm are "deadlock-free" if subcomm1
400: preceeds subcomm2 on any local list, then it preceeds subcomm2 on all ranks.
401: Equivalently, the local numbers of the subcomms on each local list are drawn from some global
402: numbering. This is ensured, for example, by ISPairToList().
404: .seealso ISPairToList()
405: @*/
406: #undef __FUNCT__
408: PetscErrorCode ISListToPair(MPI_Comm comm, PetscInt listlen, IS islist[], IS *xis, IS *yis)
409: {
411: PetscInt ncolors, *colors,i, leni,len,*xinds, *yinds,k,j;
412: const PetscInt *indsi;
415: PetscMalloc1(listlen, &colors);
416: PetscObjectsListGetGlobalNumbering(comm, listlen, (PetscObject*)islist,&ncolors, colors);
417: len = 0;
418: for (i = 0; i < listlen; ++i) {
419: ISGetLocalSize(islist[i], &leni);
420: len += leni;
421: }
422: PetscMalloc1(len, &xinds);
423: PetscMalloc1(len, &yinds);
424: k = 0;
425: for (i = 0; i < listlen; ++i) {
426: ISGetLocalSize(islist[i], &leni);
427: ISGetIndices(islist[i],&indsi);
428: for (j = 0; j < leni; ++j) {
429: xinds[k] = indsi[j];
430: yinds[k] = colors[i];
431: ++k;
432: }
433: }
434: PetscFree(colors);
435: ISCreateGeneral(comm,len,xinds,PETSC_OWN_POINTER,xis);
436: ISCreateGeneral(comm,len,yinds,PETSC_OWN_POINTER,yis);
437: return(0);
438: }
441: /*@
442: ISPairToList - convert an IS pair encoding an integer map to a list of ISs.
443: Each IS on the output list contains the preimage for each index on the second input IS.
444: The ISs on the output list are constructed on the subcommunicators of the input IS pair.
445: Each subcommunicator corresponds to the preimage of some index j -- this subcomm contains
446: exactly the ranks that assign some indices i to j. This is essentially the inverse of
447: ISListToPair().
449: Collective on indis.
451: Input arguments:
452: + xis - domain IS
453: - yis - range IS
455: Output arguments:
456: + listlen - length of islist
457: - islist - list of ISs breaking up indis by color
459: Note:
460: + xis and yis must be of the same length and have congruent communicators.
461: - The resulting ISs have subcommunicators in a "deadlock-free" order (see ISListToPair()).
463: Level: advanced
465: .seealso ISListToPair()
466: @*/
467: #undef __FUNCT__
469: PetscErrorCode ISPairToList(IS xis, IS yis, PetscInt *listlen, IS **islist)
470: {
472: IS indis = xis, coloris = yis;
473: PetscInt *inds, *colors, llen, ilen, lstart, lend, lcount,l;
474: PetscMPIInt rank, size, llow, lhigh, low, high,color,subsize;
475: const PetscInt *ccolors, *cinds;
476: MPI_Comm comm, subcomm;
484: PetscObjectGetComm((PetscObject)xis,&comm);
485: MPI_Comm_rank(comm, &rank);
486: MPI_Comm_rank(comm, &size);
487: /* Extract, copy and sort the local indices and colors on the color. */
488: ISGetLocalSize(coloris, &llen);
489: ISGetLocalSize(indis, &ilen);
490: if (llen != ilen) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Incompatible IS sizes: %D and %D", ilen, llen);
491: ISGetIndices(coloris, &ccolors);
492: ISGetIndices(indis, &cinds);
493: PetscMalloc2(ilen,&inds,llen,&colors);
494: PetscMemcpy(inds,cinds,ilen*sizeof(PetscInt));
495: PetscMemcpy(colors,ccolors,llen*sizeof(PetscInt));
496: PetscSortIntWithArray(llen, colors, inds);
497: /* Determine the global extent of colors. */
498: llow = 0; lhigh = -1;
499: lstart = 0; lcount = 0;
500: while (lstart < llen) {
501: lend = lstart+1;
502: while (lend < llen && colors[lend] == colors[lstart]) ++lend;
503: llow = PetscMin(llow,colors[lstart]);
504: lhigh = PetscMax(lhigh,colors[lstart]);
505: ++lcount;
506: }
507: MPIU_Allreduce(&llow,&low,1,MPI_INT,MPI_MIN,comm);
508: MPIU_Allreduce(&lhigh,&high,1,MPI_INT,MPI_MAX,comm);
509: *listlen = 0;
510: if (low <= high) {
511: if (lcount > 0) {
512: *listlen = lcount;
513: if (!*islist) {
514: PetscMalloc(sizeof(IS)*lcount, islist);
515: }
516: }
517: /*
518: Traverse all possible global colors, and participate in the subcommunicators
519: for the locally-supported colors.
520: */
521: lcount = 0;
522: lstart = 0; lend = 0;
523: for (l = low; l <= high; ++l) {
524: /*
525: Find the range of indices with the same color, which is not smaller than l.
526: Observe that, since colors is sorted, and is a subsequence of [low,high],
527: as soon as we find a new color, it is >= l.
528: */
529: if (lstart < llen) {
530: /* The start of the next locally-owned color is identified. Now look for the end. */
531: if (lstart == lend) {
532: lend = lstart+1;
533: while (lend < llen && colors[lend] == colors[lstart]) ++lend;
534: }
535: /* Now check whether the identified color segment matches l. */
536: 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);
537: }
538: color = (PetscMPIInt)(colors[lstart] == l);
539: /* Check whether a proper subcommunicator exists. */
540: MPIU_Allreduce(&color,&subsize,1,MPI_INT,MPI_SUM,comm);
542: if (subsize == 1) subcomm = PETSC_COMM_SELF;
543: else if (subsize == size) subcomm = comm;
544: else {
545: /* a proper communicator is necessary, so we create it. */
546: MPI_Comm_split(comm, color, rank, &subcomm);
547: }
548: if (colors[lstart] == l) {
549: /* If we have l among the local colors, we create an IS to hold the corresponding indices. */
550: ISCreateGeneral(subcomm, lend-lstart,inds+lstart,PETSC_COPY_VALUES,*islist+lcount);
551: /* Position lstart at the beginning of the next local color. */
552: lstart = lend;
553: /* Increment the counter of the local colors split off into an IS. */
554: ++lcount;
555: }
556: if (subsize > 0 && subsize < size) {
557: /*
558: Irrespective of color, destroy the split off subcomm:
559: a subcomm used in the IS creation above is duplicated
560: into a proper PETSc comm.
561: */
562: MPI_Comm_free(&subcomm);
563: }
564: } /* for (l = low; l < high; ++l) */
565: } /* if (low <= high) */
566: PetscFree2(inds,colors);
567: return(0);
568: }
571: /*@
572: ISEmbed - embed IS a into IS b by finding the locations in b that have the same indices as in a.
573: If c is the IS of these locations, we have a = b*c, regarded as a composition of the
574: corresponding ISLocalToGlobalMaps.
576: Not collective.
578: Input arguments:
579: + a - IS to embed
580: . b - IS to embed into
581: - drop - flag indicating whether to drop a's indices that are not in b.
583: Output arguments:
584: . c - local embedding indices
586: Note:
587: If some of a's global indices are not among b's indices the embedding is impossible. The local indices of a
588: corresponding to these global indices are either mapped to -1 (if !drop) or are omitted (if drop). In the former
589: case the size of c is that same as that of a, in the latter case c's size may be smaller.
591: The resulting IS is sequential, since the index substition it encodes is purely local.
593: Level: advanced
595: .seealso ISLocalToGlobalMapping
596: @*/
597: #undef __FUNCT__
599: PetscErrorCode ISEmbed(IS a, IS b, PetscBool drop, IS *c)
600: {
601: PetscErrorCode ierr;
602: ISLocalToGlobalMapping ltog;
603: ISGlobalToLocalMappingType gtoltype = IS_GTOLM_DROP;
604: PetscInt alen, clen, *cindices, *cindices2;
605: const PetscInt *aindices;
611: ISLocalToGlobalMappingCreateIS(b, <og);
612: ISGetLocalSize(a, &alen);
613: ISGetIndices(a, &aindices);
614: PetscMalloc1(alen, &cindices);
615: if (!drop) gtoltype = IS_GTOLM_MASK;
616: ISGlobalToLocalMappingApply(ltog,gtoltype,alen,aindices,&clen,cindices);
617: ISLocalToGlobalMappingDestroy(<og);
618: if (clen != alen) {
619: cindices2 = cindices;
620: PetscMalloc1(clen, &cindices);
621: PetscMemcpy(cindices,cindices2,clen*sizeof(PetscInt));
622: PetscFree(cindices2);
623: }
624: ISCreateGeneral(PETSC_COMM_SELF,clen,cindices,PETSC_OWN_POINTER,c);
625: return(0);
626: }
629: /*@
630: ISSortPermutation - calculate the permutation of the indices into a nondecreasing order.
632: Not collective.
634: Input arguments:
635: + f - IS to sort
636: - always - build the permutation even when f's indices are nondecreasin.
638: Output argument:
639: . h - permutation or NULL, if f is nondecreasing and always == PETSC_TRUE.
642: Note: Indices in f are unchanged. f[h[i]] is the i-th smallest f index.
643: If always == PETSC_FALSE, an extra check is peformed to see whether
644: the f indices are nondecreasing. h is built on PETSC_COMM_SELF, since
645: the permutation has a local meaning only.
647: Level: advanced
649: .seealso ISLocalToGlobalMapping, ISSort(), PetscIntSortWithPermutation()
650: @*/
651: #undef __FUNCT__
653: PetscErrorCode ISSortPermutation(IS f,PetscBool always,IS *h)
654: {
655: PetscErrorCode ierr;
656: const PetscInt *findices;
657: PetscInt fsize,*hindices,i;
658: PetscBool isincreasing;
663: ISGetLocalSize(f,&fsize);
664: ISGetIndices(f,&findices);
665: *h = NULL;
666: if (!always) {
667: isincreasing = PETSC_TRUE;
668: for (i = 1; i < fsize; ++i) {
669: if (findices[i] <= findices[i-1]) {
670: isincreasing = PETSC_FALSE;
671: break;
672: }
673: }
674: if (isincreasing) {
675: ISRestoreIndices(f,&findices);
676: return(0);
677: }
678: }
679: PetscMalloc1(fsize,&hindices);
680: for (i = 0; i < fsize; ++i) hindices[i] = i;
681: PetscSortIntWithPermutation(fsize,findices,hindices);
682: ISRestoreIndices(f,&findices);
683: ISCreateGeneral(PETSC_COMM_SELF,fsize,hindices,PETSC_OWN_POINTER,h);
684: return(0);
685: }