Actual source code: isdiff.c
1: #define PETSCVEC_DLL
3: #include petscis.h
4: #include petscbt.h
8: /*@
9: ISDifference - Computes the difference between two index sets.
11: Collective on IS
13: Input Parameter:
14: + is1 - first index, to have items removed from it
15: - is2 - index values to be removed
17: Output Parameters:
18: . isout - is1 - is2
20: Notes:
21: Negative values are removed from the lists. is2 may have values
22: that are not in is1. This requires O(imax-imin) memory and O(imax-imin)
23: work, where imin and imax are the bounds on the indices in is1.
25: Level: intermediate
27: Concepts: index sets^difference
28: Concepts: IS^difference
30: .seealso: ISDestroy(), ISView(), ISSum(), ISExpand()
32: @*/
33: PetscErrorCode ISDifference(IS is1,IS is2,IS *isout)
34: {
36: PetscInt i,*i1,*i2,n1,n2,imin,imax,nout,*iout;
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,isout);
90: PetscFree(iout);
92: PetscBTDestroy(mask);
93: return(0);
94: }
98: /*@
99: ISSum - Computes the sum (union) of two index sets in place. Note that
100: is1 is an existing IS, not merely a pointer.
102: Only sequential version (at the moment)
104: Input Parameter:
105: + is1 - index set to be extended
106: - is2 - index values to be added
108: Notes:
109: If n1 and n2 are the sizes of the sets, this takes O(n1+n2) time;
110: if is2 is a subset of is1, is1 is left unchanged, otherwise is1
111: is reallocated.
112: Both index sets need to be sorted on input.
114: Level: intermediate
116: .seealso: ISDestroy(), ISView(), ISDifference(), ISSum()
118: Concepts: index sets^union
119: Concepts: IS^union
121: @*/
122: PetscErrorCode ISSum(IS *is1,IS is2)
123: {
124: MPI_Comm comm;
125: PetscTruth f;
126: PetscMPIInt size;
127: PetscInt *i1,*i2,n1,n2,n3, p1,p2, *iout;
133: PetscObjectGetComm((PetscObject)(*is1),&comm);
134: MPI_Comm_size(comm,&size);
135: if (size>1) SETERRQ(PETSC_ERR_SUP,"Currently only for uni-processor IS");
137: ISSorted(*is1,&f);
138: if (!f) SETERRQ(PETSC_ERR_ARG_INCOMP,"Arg 1 is not sorted");
139: ISSorted(is2,&f);
140: if (!f) SETERRQ(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]) {n3++; p2++;}
153: if (p2==n2) { /* cleanup for is1 */ n3 += n1-p1; break;
154: } else {
155: if (i2[p2]==i1[p1]) {n3++; p1++; p2++;}
156: }
157: }
158: if (p2==n2) { /* cleanup for is1 */ n3 += n1-p1; break;
159: } else {
160: while (p1<n1 && i1[p1]<i2[p2]) {n3++; p1++;}
161: if (p1==n1) { /* clean up for is2 */ n3 += n2-p2; break;
162: } else {
163: if (i1[p1]==i2[p2]) {n3++; p1++; p2++;}
164: }
165: }
166: } while (p1<n1 || p2<n2);
168: if (n3==n1) { /* no new elements to be added */
169: ISRestoreIndices(*is1,&i1);
170: ISRestoreIndices(is2,&i2);
171: return(0);
172: }
173: PetscMalloc(n3*sizeof(PetscInt),&iout);
175: p1 = 0; p2 = 0; n3 = 0;
176: do {
177: if (p1==n1) { /* cleanup for is2 */
178: while (p2<n2) iout[n3++] = i2[p2++];
179: break;
180: } else {
181: while (p2<n2 && i2[p2]<i1[p1]) iout[n3++] = i2[p2++];
182: if (p2==n2) { /* cleanup for is1 */
183: while (p1<n1) iout[n3++] = i1[p1++];
184: break;
185: } else {
186: if (i2[p2]==i1[p1]) {iout[n3++] = i1[p1++]; p2++;}
187: }
188: }
189: if (p2==n2) { /* cleanup for is1 */
190: while (p1<n1) iout[n3++] = i1[p1++];
191: break;
192: } else {
193: while (p1<n1 && i1[p1]<i2[p2]) iout[n3++] = i1[p1++];
194: if (p1==n1) { /* clean up for is2 */
195: while (p2<n2) iout[n3++] = i2[p2++];
196: break;
197: } else {
198: if (i1[p1]==i2[p2]) {iout[n3++] = i1[p1++]; p2++;}
199: }
200: }
201: } while (p1<n1 || p2<n2);
203: ISRestoreIndices(*is1,&i1);
204: ISRestoreIndices(is2,&i2);
205: ISDestroy(*is1);
206: ISCreateGeneral(PETSC_COMM_SELF,n3,iout,is1);
207: PetscFree(iout);
209: return(0);
210: }
214: /*@
215: ISExpand - Computes the union of two index sets, by concatenating 2 lists and
216: removing duplicates.
218: Collective on IS
220: Input Parameter:
221: + is1 - first index set
222: - is2 - index values to be added
224: Output Parameters:
225: . isout - is1 + is2 The index set is2 is appended to is1 removing duplicates
227: Notes:
228: Negative values are removed from the lists. This requires O(imax-imin)
229: memory and O(imax-imin) work, where imin and imax are the bounds on the
230: indices in is1 and is2.
232: Level: intermediate
234: .seealso: ISDestroy(), ISView(), ISDifference(), ISSum()
236: Concepts: index sets^difference
237: Concepts: IS^difference
239: @*/
240: PetscErrorCode ISExpand(IS is1,IS is2,IS *isout)
241: {
243: PetscInt i,*i1,*i2,n1,n2,imin,imax,nout,*iout;
244: PetscBT mask;
245: MPI_Comm comm;
252: ISGetIndices(is1,&i1);
253: ISGetLocalSize(is1,&n1);
254: ISGetIndices(is2,&i2);
255: ISGetLocalSize(is2,&n2);
257: /* Create a bit mask array to contain required values */
258: if (n1 || n2) {
259: imin = PETSC_MAX_INT;
260: imax = 0;
261: for (i=0; i<n1; i++) {
262: if (i1[i] < 0) continue;
263: imin = PetscMin(imin,i1[i]);
264: imax = PetscMax(imax,i1[i]);
265: }
266: for (i=0; i<n2; i++) {
267: if (i2[i] < 0) continue;
268: imin = PetscMin(imin,i2[i]);
269: imax = PetscMax(imax,i2[i]);
270: }
271: } else {
272: imin = imax = 0;
273: }
274: PetscMalloc((n1+n2)*sizeof(PetscInt),&iout);
275: nout = 0;
276: PetscBTCreate(imax-imin,mask);
277: /* Put the values from is1 */
278: for (i=0; i<n1; i++) {
279: if (i1[i] < 0) continue;
280: if (!PetscBTLookupSet(mask,i1[i] - imin)) {
281: iout[nout++] = i1[i];
282: }
283: }
284: ISRestoreIndices(is1,&i1);
285: /* Put the values from is2 */
286: for (i=0; i<n2; i++) {
287: if (i2[i] < 0) continue;
288: if (!PetscBTLookupSet(mask,i2[i] - imin)) {
289: iout[nout++] = i2[i];
290: }
291: }
292: ISRestoreIndices(is2,&i2);
294: /* create the new IS containing the sum */
295: PetscObjectGetComm((PetscObject)is1,&comm);
296: ISCreateGeneral(comm,nout,iout,isout);
297: PetscFree(iout);
299: PetscBTDestroy(mask);
300: return(0);
301: }