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