Actual source code: isblock.c

      1: /* Routines to be used by MatIncreaseOverlap() for BAIJ and SBAIJ matrices */
      2: #include <petscis.h>
      3: #include <petscbt.h>
      4: #include <petscctable.h>

      6: /*@
      7:    ISCompressIndicesGeneral - convert the indices into block indices
      8:    Input Parameters:
      9: +  n - maximum possible length of the index set
     10: .  nkeys - expected number of keys when PETSC_USE_CTABLE
     11: .  bs - the size of block
     12: .  imax - the number of index sets
     13: -  is_in - the non-blocked array of index sets

     15:    Output Parameter:
     16: .  is_out - the blocked new index set

     18:    Level: intermediate

     20: .seealso: ISExpandIndicesGeneral()
     21: @*/
     22: PetscErrorCode  ISCompressIndicesGeneral(PetscInt n,PetscInt nkeys,PetscInt bs,PetscInt imax,const IS is_in[],IS is_out[])
     23: {
     24:   PetscErrorCode     ierr;
     25:   PetscInt           isz,len,i,j,ival,Nbs;
     26:   const PetscInt     *idx;
     27: #if defined(PETSC_USE_CTABLE)
     28:   PetscTable         gid1_lid1;
     29:   PetscInt           tt, gid1, *nidx,Nkbs;
     30:   PetscTablePosition tpos;
     31: #else
     32:   PetscInt           *nidx;
     33:   PetscBT            table;
     34: #endif

     37:   Nbs = n/bs;
     38: #if defined(PETSC_USE_CTABLE)
     39:   Nkbs = nkeys/bs;
     40:   PetscTableCreate(Nkbs,Nbs,&gid1_lid1);
     41: #else
     42:   PetscMalloc1(Nbs,&nidx);
     43:   PetscBTCreate(Nbs,&table);
     44: #endif
     45:   for (i=0; i<imax; i++) {
     46:     isz = 0;
     47: #if defined(PETSC_USE_CTABLE)
     48:     PetscTableRemoveAll(gid1_lid1);
     49: #else
     50:     PetscBTMemzero(Nbs,table);
     51: #endif
     52:     ISGetIndices(is_in[i],&idx);
     53:     ISGetLocalSize(is_in[i],&len);
     54:     for (j=0; j<len; j++) {
     55:       ival = idx[j]/bs; /* convert the indices into block indices */
     56: #if defined(PETSC_USE_CTABLE)
     57:       PetscTableFind(gid1_lid1,ival+1,&tt);
     58:       if (!tt) {
     59:         PetscTableAdd(gid1_lid1,ival+1,isz+1,INSERT_VALUES);
     60:         isz++;
     61:       }
     62: #else
     63:       if (ival>Nbs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"index greater than mat-dim");
     64:       if (!PetscBTLookupSet(table,ival)) nidx[isz++] = ival;
     65: #endif
     66:     }
     67:     ISRestoreIndices(is_in[i],&idx);

     69: #if defined(PETSC_USE_CTABLE)
     70:     PetscMalloc1(isz,&nidx);
     71:     PetscTableGetHeadPosition(gid1_lid1,&tpos);
     72:     j    = 0;
     73:     while (tpos) {
     74:       PetscTableGetNext(gid1_lid1,&tpos,&gid1,&tt);
     75:       if (tt-- > isz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"index greater than array-dim");
     76:       nidx[tt] = gid1 - 1;
     77:       j++;
     78:     }
     79:     if (j != isz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"table error: jj != isz");
     80:     ISCreateGeneral(PetscObjectComm((PetscObject)is_in[i]),isz,nidx,PETSC_OWN_POINTER,(is_out+i));
     81: #else
     82:     ISCreateGeneral(PetscObjectComm((PetscObject)is_in[i]),isz,nidx,PETSC_COPY_VALUES,(is_out+i));
     83: #endif
     84:   }
     85: #if defined(PETSC_USE_CTABLE)
     86:   PetscTableDestroy(&gid1_lid1);
     87: #else
     88:   PetscBTDestroy(&table);
     89:   PetscFree(nidx);
     90: #endif
     91:   return(0);
     92: }

     94: PetscErrorCode  ISCompressIndicesSorted(PetscInt n,PetscInt bs,PetscInt imax,const IS is_in[],IS is_out[])
     95: {
     97:   PetscInt       i,j,k,val,len,*nidx,bbs;
     98:   const PetscInt *idx,*idx_local;
     99:   PetscBool      flg,isblock;
    100: #if defined(PETSC_USE_CTABLE)
    101:   PetscInt       maxsz;
    102: #else
    103:   PetscInt       Nbs=n/bs;
    104: #endif

    107:   for (i=0; i<imax; i++) {
    108:     ISSorted(is_in[i],&flg);
    109:     if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Indices are not sorted");
    110:   }

    112: #if defined(PETSC_USE_CTABLE)
    113:   /* Now check max size */
    114:   for (i=0,maxsz=0; i<imax; i++) {
    115:     ISGetLocalSize(is_in[i],&len);
    116:     if (len%bs !=0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Indices are not block ordered");
    117:     len = len/bs; /* The reduced index size */
    118:     if (len > maxsz) maxsz = len;
    119:   }
    120:   PetscMalloc1(maxsz,&nidx);
    121: #else
    122:   PetscMalloc1(Nbs,&nidx);
    123: #endif
    124:   /* Now check if the indices are in block order */
    125:   for (i=0; i<imax; i++) {
    126:     ISGetLocalSize(is_in[i],&len);

    128:     /* special case where IS is already block IS of the correct size */
    129:     PetscObjectTypeCompare((PetscObject)is_in[i],ISBLOCK,&isblock);
    130:     if (isblock) {
    131:       ISBlockGetLocalSize(is_in[i],&bbs);
    132:       if (bs == bbs) {
    133:         len  = len/bs;
    134:         ISBlockGetIndices(is_in[i],&idx);
    135:         ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,(is_out+i));
    136:         ISBlockRestoreIndices(is_in[i],&idx);
    137:         continue;
    138:       }
    139:     }
    140:     ISGetIndices(is_in[i],&idx);
    141:     if (len%bs !=0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Indices are not block ordered");

    143:     len       = len/bs; /* The reduced index size */
    144:     idx_local = idx;
    145:     for (j=0; j<len; j++) {
    146:       val = idx_local[0];
    147:       if (val%bs != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Indices are not block ordered");
    148:       for (k=0; k<bs; k++) {
    149:         if (val+k != idx_local[k]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Indices are not block ordered");
    150:       }
    151:       nidx[j]    = val/bs;
    152:       idx_local += bs;
    153:     }
    154:     ISRestoreIndices(is_in[i],&idx);
    155:     ISCreateGeneral(PetscObjectComm((PetscObject)is_in[i]),len,nidx,PETSC_COPY_VALUES,(is_out+i));
    156:   }
    157:   PetscFree(nidx);
    158:   return(0);
    159: }

    161: /*@C
    162:    ISExpandIndicesGeneral - convert the indices into non-block indices
    163:    Input Parameters:
    164: +  n - the length of the index set   (not being used)
    165: .  nkeys - expected number of keys when PETSC_USE_CTABLE (not being used)
    166: .  bs - the size of block
    167: .  imax - the number of index sets
    168: -  is_in - the blocked array of index sets

    170:    Output Parameter:
    171: .  is_out - the non-blocked new index set

    173:    Level: intermediate

    175: .seealso: ISCompressIndicesGeneral()
    176: @*/
    177: PetscErrorCode  ISExpandIndicesGeneral(PetscInt n,PetscInt nkeys,PetscInt bs,PetscInt imax,const IS is_in[],IS is_out[])
    178: {
    180:   PetscInt       len,i,j,k,*nidx;
    181:   const PetscInt *idx;
    182:   PetscInt       maxsz;

    185:   /* Check max size of is_in[] */
    186:   maxsz = 0;
    187:   for (i=0; i<imax; i++) {
    188:     ISGetLocalSize(is_in[i],&len);
    189:     if (len > maxsz) maxsz = len;
    190:   }
    191:   PetscMalloc1(maxsz*bs,&nidx);

    193:   for (i=0; i<imax; i++) {
    194:     ISGetLocalSize(is_in[i],&len);
    195:     ISGetIndices(is_in[i],&idx);
    196:     for (j=0; j<len ; ++j) {
    197:       for (k=0; k<bs; k++) nidx[j*bs+k] = idx[j]*bs+k;
    198:     }
    199:     ISRestoreIndices(is_in[i],&idx);
    200:     ISCreateGeneral(PetscObjectComm((PetscObject)is_in[i]),len*bs,nidx,PETSC_COPY_VALUES,is_out+i);
    201:   }
    202:   PetscFree(nidx);
    203:   return(0);
    204: }