Actual source code: isdiff.c
petsc-3.3-p7 2013-05-11
2: #include <petscis.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 {
58: imin = imax = 0;
59: }
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);
75:
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: PetscMalloc(nout*sizeof(PetscInt),&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: }
97: /*@
98: ISSum - Computes the sum (union) of two index sets.
100: Only sequential version (at the moment)
102: Input Parameter:
103: + is1 - index set to be extended
104: - is2 - index values to be added
106: Output Parameter:
107: . is3 - the sum; this can not be is1 or is2
109: Notes:
110: If n1 and n2 are the sizes of the sets, this takes O(n1+n2) time;
112: Both index sets need to be sorted on input.
114: Level: intermediate
116: .seealso: ISDestroy(), ISView(), ISDifference(), ISExpand()
118: Concepts: index sets^union
119: Concepts: IS^union
121: @*/
122: PetscErrorCode ISSum(IS is1,IS is2,IS *is3)
123: {
124: MPI_Comm comm;
125: PetscBool f;
126: PetscMPIInt size;
127: const PetscInt *i1,*i2;
128: PetscInt n1,n2,n3, p1,p2, *iout;
134: PetscObjectGetComm((PetscObject)(is1),&comm);
135: MPI_Comm_size(comm,&size);
136: if (size>1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently only for uni-processor IS");
138: ISSorted(is1,&f);
139: if (!f) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Arg 1 is not sorted");
140: ISSorted(is2,&f);
141: if (!f) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Arg 2 is not sorted");
143: ISGetLocalSize(is1,&n1);
144: ISGetLocalSize(is2,&n2);
145: if (!n2) return(0);
146: ISGetIndices(is1,&i1);
147: ISGetIndices(is2,&i2);
149: p1 = 0; p2 = 0; n3 = 0;
150: do {
151: if (p1==n1) { /* cleanup for is2 */ n3 += n2-p2; break;
152: } else {
153: while (p2<n2 && i2[p2]<i1[p1]) {n3++; p2++;}
154: if (p2==n2) { /* cleanup for is1 */ n3 += n1-p1; break;
155: } else {
156: if (i2[p2]==i1[p1]) {n3++; p1++; p2++;}
157: }
158: }
159: if (p2==n2) { /* cleanup for is1 */ n3 += n1-p1; break;
160: } else {
161: while (p1<n1 && i1[p1]<i2[p2]) {n3++; p1++;}
162: if (p1==n1) { /* clean up for is2 */ n3 += n2-p2; break;
163: } else {
164: if (i1[p1]==i2[p2]) {n3++; p1++; p2++;}
165: }
166: }
167: } while (p1<n1 || p2<n2);
169: if (n3==n1) { /* no new elements to be added */
170: ISRestoreIndices(is1,&i1);
171: ISRestoreIndices(is2,&i2);
172: ISDuplicate(is1,is3);
173: return(0);
174: }
175: PetscMalloc(n3*sizeof(PetscInt),&iout);
177: p1 = 0; p2 = 0; n3 = 0;
178: do {
179: if (p1==n1) { /* cleanup for is2 */
180: while (p2<n2) iout[n3++] = i2[p2++];
181: break;
182: } else {
183: while (p2<n2 && i2[p2]<i1[p1]) iout[n3++] = i2[p2++];
184: if (p2==n2) { /* cleanup for is1 */
185: while (p1<n1) iout[n3++] = i1[p1++];
186: break;
187: } else {
188: if (i2[p2]==i1[p1]) {iout[n3++] = i1[p1++]; p2++;}
189: }
190: }
191: if (p2==n2) { /* cleanup for is1 */
192: while (p1<n1) iout[n3++] = i1[p1++];
193: break;
194: } else {
195: while (p1<n1 && i1[p1]<i2[p2]) iout[n3++] = i1[p1++];
196: if (p1==n1) { /* clean up for is2 */
197: while (p2<n2) iout[n3++] = i2[p2++];
198: break;
199: } else {
200: if (i1[p1]==i2[p2]) {iout[n3++] = i1[p1++]; p2++;}
201: }
202: }
203: } while (p1<n1 || p2<n2);
205: ISRestoreIndices(is1,&i1);
206: ISRestoreIndices(is2,&i2);
207: ISCreateGeneral(comm,n3,iout,PETSC_OWN_POINTER,is3);
208: return(0);
209: }
213: /*@
214: ISExpand - Computes the union of two index sets, by concatenating 2 lists and
215: removing duplicates.
217: Collective on IS
219: Input Parameter:
220: + is1 - first index set
221: - is2 - index values to be added
223: Output Parameters:
224: . isout - is1 + is2 The index set is2 is appended to is1 removing duplicates
226: Notes:
227: Negative values are removed from the lists. This requires O(imax-imin)
228: memory and O(imax-imin) work, where imin and imax are the bounds on the
229: indices in is1 and is2.
231: The IS's do not need to be sorted.
233: Level: intermediate
235: .seealso: ISDestroy(), ISView(), ISDifference(), ISSum()
237: Concepts: index sets^difference
238: Concepts: IS^difference
240: @*/
241: PetscErrorCode ISExpand(IS is1,IS is2,IS *isout)
242: {
244: PetscInt i,n1,n2,imin,imax,nout,*iout;
245: const PetscInt *i1,*i2;
246: PetscBT mask;
247: MPI_Comm comm;
254: ISGetIndices(is1,&i1);
255: ISGetLocalSize(is1,&n1);
256: ISGetIndices(is2,&i2);
257: ISGetLocalSize(is2,&n2);
259: /* Create a bit mask array to contain required values */
260: if (n1 || n2) {
261: imin = PETSC_MAX_INT;
262: imax = 0;
263: for (i=0; i<n1; i++) {
264: if (i1[i] < 0) continue;
265: imin = PetscMin(imin,i1[i]);
266: imax = PetscMax(imax,i1[i]);
267: }
268: for (i=0; i<n2; i++) {
269: if (i2[i] < 0) continue;
270: imin = PetscMin(imin,i2[i]);
271: imax = PetscMax(imax,i2[i]);
272: }
273: } else {
274: imin = imax = 0;
275: }
276: PetscMalloc((n1+n2)*sizeof(PetscInt),&iout);
277: nout = 0;
278: PetscBTCreate(imax-imin,&mask);
279: /* Put the values from is1 */
280: for (i=0; i<n1; i++) {
281: if (i1[i] < 0) continue;
282: if (!PetscBTLookupSet(mask,i1[i] - imin)) {
283: iout[nout++] = i1[i];
284: }
285: }
286: ISRestoreIndices(is1,&i1);
287: /* Put the values from is2 */
288: for (i=0; i<n2; i++) {
289: if (i2[i] < 0) continue;
290: if (!PetscBTLookupSet(mask,i2[i] - imin)) {
291: iout[nout++] = i2[i];
292: }
293: }
294: ISRestoreIndices(is2,&i2);
296: /* create the new IS containing the sum */
297: PetscObjectGetComm((PetscObject)is1,&comm);
298: ISCreateGeneral(comm,nout,iout,PETSC_OWN_POINTER,isout);
300: PetscBTDestroy(&mask);
301: return(0);
302: }
306: /*@
307: ISConcatenate - Forms a new IS by locally concatenating the indices from an IS list without reordering.
308:
310: Collective on comm.
312: Input Parameter:
313: + comm - communicator of the concatenated IS.
314: . len - size of islist array (nonnegative)
315: - islist - array of index sets
319: Output Parameters:
320: . isout - The concatenated index set; empty, if len == 0.
322: Notes: The semantics of calling this on comm imply that the comms of the members if islist also contain this rank.
324: Level: intermediate
326: .seealso: ISDifference(), ISSum(), ISExpand()
328: Concepts: index sets^concatenation
329: Concepts: IS^concatenation
331: @*/
332: PetscErrorCode ISConcatenate(MPI_Comm comm, PetscInt len, const IS islist[], IS *isout)
333: {
335: PetscInt i,n,N;
336: const PetscInt *iidx;
337: PetscInt *idx;
341: #if defined(PETSC_USE_DEBUG)
342: for(i = 0; i < len; ++i) {
344: }
345: #endif
347: if(!len) {
348: ISCreateStride(comm, 0,0,0, isout);
349: return(0);
350: }
351: if(len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Negative array length: %D", len);
352: N = 0;
353: for(i = 0; i < len; ++i) {
354: ISGetLocalSize(islist[i], &n);
355: N += n;
356: }
357: PetscMalloc(sizeof(PetscInt)*N, &idx);
358: N = 0;
359: for(i = 0; i < len; ++i) {
360: ISGetLocalSize(islist[i], &n);
361: ISGetIndices(islist[i], &iidx);
362: PetscMemcpy(idx+N,iidx, sizeof(PetscInt)*n);
363: N += n;
364: }
365: ISCreateGeneral(comm, N, idx, PETSC_OWN_POINTER, isout);
366: return(0);
367: }
369: /*@
370: ISListToMap - convert an IS list to a pair of ISs of equal length defining an equivalent integer multimap.
371: Each IS on the input list is assigned an integer j so that all of the indices of that IS are
372: mapped to j.
375: Collective on comm.
377: Input arguments:
378: + comm - MPI_Comm
379: . listlen - IS list length
380: - islist - IS list
382: Output arguments:
383: + xis - domain IS
384: - yis - range IS
386: Level: advanced
388: Notes:
389: The global integers assigned to the ISs of the local input list might not correspond to the
390: local numbers of the ISs on that list, but the two *orderings* are the same: the global
391: integers assigned to the ISs on the local list form a strictly increasing sequence.
393: The ISs on the input list can belong to subcommunicators of comm, and the subcommunicators
394: on the input IS list are assumed to be in a "deadlock-free" order.
396: Local lists of PetscObjects (or their subcommes) on a comm are "deadlock-free" if subcomm1
397: preceeds subcomm2 on any local list, then it preceeds subcomm2 on all ranks.
398: Equivalently, the local numbers of the subcomms on each local list are drawn from some global
399: numbering. This is ensured, for example, by ISMapToList().
401: .seealso ISMapToList()
402: @*/
403: #undef __FUNCT__
405: PetscErrorCode ISListToMap(MPI_Comm comm, PetscInt listlen, IS islist[], IS *xis, IS *yis)
406: {
408: PetscInt ncolors, *colors,i, leni,len,*xinds, *yinds,k,j;
409: const PetscInt *indsi;
411: PetscMalloc(listlen*sizeof(PetscInt), &colors);
412: PetscObjectsGetGlobalNumbering(comm, listlen, (PetscObject*)islist,&ncolors, colors);
413: len = 0;
414: for(i = 0; i < listlen; ++i) {
415: ISGetLocalSize(islist[i], &leni);
416: len += leni;
417: }
418: PetscMalloc(len*sizeof(PetscInt), &xinds);
419: PetscMalloc(len*sizeof(PetscInt), &yinds);
420: k = 0;
421: for(i = 0; i < listlen; ++i) {
422: ISGetLocalSize(islist[i], &leni);
423: ISGetIndices(islist[i],&indsi);
424: for(j = 0; j < leni; ++j) {
425: xinds[k] = indsi[j];
426: yinds[k] = colors[i];
427: ++k;
428: }
429: }
430: PetscFree(colors);
431: ISCreateGeneral(comm,len,xinds,PETSC_OWN_POINTER,xis);
432: ISCreateGeneral(comm,len,yinds,PETSC_OWN_POINTER,yis);
433: return(0);
434: }
437: /*@
438: ISMapToList - convert an IS pair encoding an integer map to a list of ISs.
439: Each IS on the output list contains the preimage for each index on the second input IS.
440: The ISs on the output list are constructed on the subcommunicators of the input IS pair.
441: Each subcommunicator corresponds to the preimage of some index j -- this subcomm contains
442: exactly the ranks that assign some indices i to j. This is essentially the inverse of
443: ISListToMap().
445: Collective on indis.
447: Input arguments:
448: + xis - domain IS
449: - yis - range IS
451: Output arguments:
452: + listlen - length of islist
453: - islist - list of ISs breaking up indis by color
455: Note:
456: + xis and yis must be of the same length and have congruent communicators.
457: - The resulting ISs have subcommunicators in a "deadlock-free" order (see ISListToMap()).
459: Level: advanced
461: .seealso ISListToMap()
462: @*/
463: #undef __FUNCT__
465: PetscErrorCode ISMapToList(IS xis, IS yis, PetscInt *listlen, IS **islist)
466: {
468: IS indis = xis, coloris = yis;
469: PetscInt *inds, *colors, llen, ilen, lstart, lend, lcount,l;
470: PetscMPIInt rank, size, llow, lhigh, low, high,color,subsize;
471: const PetscInt *ccolors, *cinds;
472: MPI_Comm comm, subcomm;
479: comm = ((PetscObject)xis)->comm;
480: MPI_Comm_rank(comm, &rank);
481: MPI_Comm_rank(comm, &size);
482: /* Extract, copy and sort the local indices and colors on the color. */
483: ISGetLocalSize(coloris, &llen);
484: ISGetLocalSize(indis, &ilen);
485: if(llen != ilen) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Incompatible IS sizes: %D and %D", ilen, llen);
486: ISGetIndices(coloris, &ccolors);
487: ISGetIndices(indis, &cinds);
488: PetscMalloc2(ilen,PetscInt,&inds,llen,PetscInt,&colors);
489: PetscMemcpy(inds,cinds,ilen*sizeof(PetscInt));
490: PetscMemcpy(colors,ccolors,llen*sizeof(PetscInt));
491: PetscSortIntWithArray(llen, colors, inds);
492: /* Determine the global extent of colors. */
493: llow = 0; lhigh = -1;
494: lstart = 0; lcount = 0;
495: while(lstart < llen) {
496: lend = lstart+1;
497: while(lend < llen && colors[lend] == colors[lstart]) ++lend;
498: llow = PetscMin(llow,colors[lstart]);
499: lhigh = PetscMax(lhigh,colors[lstart]);
500: ++lcount;
501: }
502: MPI_Allreduce(&llow,&low,1,MPI_INT,MPI_MIN,comm);
503: MPI_Allreduce(&lhigh,&high,1,MPI_INT,MPI_MAX,comm);
504: *listlen = 0;
505: if(low <= high) {
506: if(lcount > 0) {
507: *listlen = lcount;
508: if(!*islist) {
509: PetscMalloc(sizeof(IS)*lcount, islist);
510: }
511: }
512: /*
513: Traverse all possible global colors, and participate in the subcommunicators
514: for the locally-supported colors.
515: */
516: lcount = 0;
517: lstart = 0; lend = 0;
518: for(l = low; l <= high; ++l) {
519: /*
520: Find the range of indices with the same color, which is not smaller than l.
521: Observe that, since colors is sorted, and is a subsequence of [low,high],
522: as soon as we find a new color, it is >= l.
523: */
524: if(lstart < llen) {
525: /* The start of the next locally-owned color is identified. Now look for the end. */
526: if(lstart == lend) {
527: lend = lstart+1;
528: while(lend < llen && colors[lend] == colors[lstart]) ++lend;
529: }
530: /* Now check whether the identified color segment matches l. */
531: 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);
532: }
533: color = (PetscMPIInt)(colors[lstart] == l);
534: /* Check whether a proper subcommunicator exists. */
535: MPI_Allreduce(&color,&subsize,1,MPI_INT,MPI_SUM,comm);
536:
537: if(subsize == 1) subcomm = PETSC_COMM_SELF;
538: else if(subsize == size) subcomm = comm;
539: else {
540: /* a proper communicator is necessary, so we create it. */
541: MPI_Comm_split(comm, color, rank, &subcomm);
542: }
543: if(colors[lstart] == l) {
544: /* If we have l among the local colors, we create an IS to hold the corresponding indices. */
545: ISCreateGeneral(subcomm, lend-lstart,inds+lstart,PETSC_COPY_VALUES,*islist+lcount);
546: /* Position lstart at the beginning of the next local color. */
547: lstart = lend;
548: /* Increment the counter of the local colors split off into an IS. */
549: ++lcount;
550: }
551: if(subsize > 0 && subsize < size) {
552: /*
553: Irrespective of color, destroy the split off subcomm:
554: a subcomm used in the IS creation above is duplicated
555: into a proper PETSc comm.
556: */
557: MPI_Comm_free(&subcomm);
558: }
559: }/* for(l = low; l < high; ++l) */
560: }/* if(low <= high) */
561: PetscFree2(inds,colors);
562: return(0);
563: }
566: /*@
567: ISMapFactorRight - for a pair of ISs a and b, regarded as local-to-global index maps, compute IS c such that
568: a = b*c as a composition of maps. In other words, find a substitution of local indices c
569: such that a factors through c (and b). Another way to look at this is as finding the right
570: factor for b in a (b is the left factor).
572: Not collective.
574: Input arguments:
575: + a - IS to factor
576: . b - left factor
577: - drop - flag indicating whether to drop a's indices that can't factor through b.
579: Output arguments:
580: . c - right local factor
582: Note:
583: If some of a's global indices are not among b's indices the factorization is impossible. The local indices of a
584: corresponding to these global indices are either mapped to -1 (if !drop) or are omitted (if drop). In former
585: case the size of c is that same as that of a, in the latter case c's size may be smaller.
587: The resulting IS is sequential, since the index substition it encodes is purely local.
589: Level: advanced
591: .seealso ISLocalToGlobalMapping
592: @*/
593: #undef __FUNCT__
595: PetscErrorCode ISMapFactorRight(IS a, IS b, PetscBool drop, IS *c)
596: {
598: ISLocalToGlobalMapping ltog;
599: ISGlobalToLocalMappingType gtoltype = IS_GTOLM_DROP;
600: PetscInt alen, clen, *cindices, *cindices2;
601: const PetscInt *aindices;
607: ISLocalToGlobalMappingCreateIS(b, <og);
608: ISGetLocalSize(a, &alen);
609: ISGetIndices(a, &aindices);
610: PetscMalloc(alen*sizeof(PetscInt), &cindices);
611: if(!drop) gtoltype = IS_GTOLM_MASK;
612: ISGlobalToLocalMappingApply(ltog,gtoltype,alen,aindices,&clen,cindices);
613: if(clen != alen) {
614: cindices2 = cindices;
615: PetscMalloc(clen*sizeof(PetscInt), &cindices);
616: PetscMemcpy(cindices,cindices2,clen*sizeof(PetscInt));
617: PetscFree(cindices2);
618: }
619: ISCreateGeneral(PETSC_COMM_SELF, clen, cindices, PETSC_OWN_POINTER, c);
620: return(0);
621: }