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: }