Actual source code: isblock.c

      1: #define PETSCVEC_DLL

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

     11: /*@C
     12:    ISCompressIndicesGeneral - convert the indices into block indices
     13:    Input Parameters:
     14: +  n - the length of the index set
     15: .  bs - the size of block
     16: .  imax - the number of index sets
     17: -  is_in - the non-blocked array of index sets

     19:    Output Parameter:
     20: .  is_out - the blocked new index set

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

     38:   Nbs =n/bs;
     39: #if defined (PETSC_USE_CTABLE)
     40:   PetscTableCreate(Nbs,&gid1_lid1);
     41: #else
     42:   PetscMalloc(Nbs*sizeof(PetscInt),&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);
     60:         isz++;
     61:       }
     62: #else
     63:       if (ival>Nbs) SETERRQ(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);
     68: #if defined (PETSC_USE_CTABLE)
     69:     PetscMalloc(isz*sizeof(PetscInt),&nidx);
     70:     PetscTableGetHeadPosition(gid1_lid1,&tpos);
     71:     j = 0;
     72:     while (tpos) {
     73:       PetscTableGetNext(gid1_lid1,&tpos,&gid1,&tt);
     74:       if (tt-- > isz) { SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"index greater than array-dim"); }
     75:       nidx[tt] = gid1 - 1;
     76:       j++;
     77:     }
     78:     if (j != isz) { SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"table error: jj != isz"); }
     79:     ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is_out+i));
     80:     PetscFree(nidx);
     81: #else
     82:     ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is_out+i));
     83: #endif
     84:   }
     85: #if defined (PETSC_USE_CTABLE)
     86:   PetscTableDelete(gid1_lid1);
     87: #else
     88:   PetscBTDestroy(table);
     89:   PetscFree(nidx);
     90: #endif
     91:   return(0);
     92: }

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

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

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

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

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

    161:   return(0);
    162: }

    166: PetscErrorCode PETSCVEC_DLLEXPORT ISExpandIndicesGeneral(PetscInt n,PetscInt bs,PetscInt imax,const IS is_in[],IS is_out[])
    167: {
    169:   PetscInt       len,i,j,k,*idx,*nidx;
    170: #if defined (PETSC_USE_CTABLE)
    171:   PetscInt       maxsz;
    172: #else
    173:   PetscInt       Nbs;
    174: #endif

    177: #if defined (PETSC_USE_CTABLE)
    178:   /* Now check max size */
    179:   for (i=0,maxsz=0; i<imax; i++) {
    180:     ISGetIndices(is_in[i],&idx);
    181:     ISGetLocalSize(is_in[i],&len);
    182:     if (len*bs > maxsz) maxsz = len*bs;
    183:   }
    184:   PetscMalloc(maxsz*sizeof(PetscInt),&nidx);
    185: #else
    186:   Nbs = n/bs;
    187:   PetscMalloc(Nbs*bs*sizeof(PetscInt),&nidx);
    188: #endif

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