Actual source code: aoreduced.c

      1: #define PETSCDM_DLL

     3:  #include src/dm/ao/aoimpl.h
     4:  #include petscsys.h
     5:  #include petscbt.h

      9: PetscErrorCode AODataSegmentGetReduced_Basic(AOData ao,const char name[],const char segname[],PetscInt n,PetscInt *keys,IS *is)
     10: {
     11:   AODataSegment  *segment;
     12:   AODataKey      *key;
     14:   PetscInt       dsize,i,bs,*found,count,imin,imax,*out;
     15:   char           *idata,*odata;
     16:   PetscBT        mask;
     17:   PetscTruth     flag;

     20:   /* find the correct segment */
     21:   AODataSegmentFind_Private(ao,name,segname,&flag,&key,&segment);
     22:   if (!flag) SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot locate segment");

     24:   if (segment->datatype != PETSC_INT) SETERRQ(PETSC_ERR_ARG_WRONG,"Only for PETSC_INT data");

     26:   /*
     27:      Copy the found values into a contiguous location, keeping them in the
     28:      order of the requested keys
     29:   */
     30:   PetscDataTypeGetSize(segment->datatype,&dsize);
     31:   bs    = segment->bs;
     32:   PetscMalloc((n+1)*bs*dsize,&odata);
     33:   idata = (char*)segment->data;
     34:   for (i=0; i<n; i++) {
     35:     PetscMemcpy(odata + i*bs*dsize,idata + keys[i]*bs*dsize,bs*dsize);
     36:   }

     38:   found = (PetscInt*)odata;
     39:   n     = n*bs;

     41:   /*  Determine the max and min values */
     42:   if (n) {
     43:     imin = PETSC_MAX_INT;
     44:     imax = 0;
     45:     for (i=0; i<n; i++) {
     46:       if (found[i] < 0) continue;
     47:       imin = PetscMin(imin,found[i]);
     48:       imax = PetscMax(imax,found[i]);
     49:     }
     50:   } else {
     51:     imin = imax = 0;
     52:   }
     53:   PetscBTCreate(imax-imin,mask);
     54:   /* Put the values into the mask and count them */
     55:   count = 0;
     56:   for (i=0; i<n; i++) {
     57:     if (found[i] < 0) continue;
     58:     if (!PetscBTLookupSet(mask,found[i] - imin)) count++;
     59:   }
     60:   PetscBTMemzero(imax-imin,mask);
     61:   PetscMalloc((count+1)*sizeof(PetscInt),&out);
     62:   count = 0;
     63:   for (i=0; i<n; i++) {
     64:     if (found[i] < 0) continue;
     65:     if (!PetscBTLookupSet(mask,found[i] - imin)) {out[count++] = found[i];}
     66:   }
     67:   PetscBTDestroy(mask);
     68:   PetscFree(found);

     70:   ISCreateGeneral(ao->comm,count,out,is);
     71:   PetscFree(out);
     72:   return(0);
     73: }