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

  9:    Input Parameters:
 10: +    n - maximum possible length of the index set
 11: .    nkeys - expected number of keys when PETSC_USE_CTABLE
 12: .    bs - the size of block
 13: .    imax - the number of index sets
 14: -    is_in - the non-blocked array of index sets

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

 19:    Level: intermediate

 21: .seealso: ISExpandIndicesGeneral()
 22: @*/
 23: PetscErrorCode  ISCompressIndicesGeneral(PetscInt n,PetscInt nkeys,PetscInt bs,PetscInt imax,const IS is_in[],IS is_out[])
 24: {
 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

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

 68: #if defined(PETSC_USE_CTABLE)
 69:     PetscMalloc1(isz,&nidx);
 70:     PetscTableGetHeadPosition(gid1_lid1,&tpos);
 71:     j    = 0;
 72:     while (tpos) {
 73:       PetscTableGetNext(gid1_lid1,&tpos,&gid1,&tt);
 75:       nidx[tt] = gid1 - 1;
 76:       j++;
 77:     }
 79:     ISCreateGeneral(PetscObjectComm((PetscObject)is_in[i]),isz,nidx,PETSC_OWN_POINTER,(is_out+i));
 80: #else
 81:     ISCreateGeneral(PetscObjectComm((PetscObject)is_in[i]),isz,nidx,PETSC_COPY_VALUES,(is_out+i));
 82: #endif
 83:   }
 84: #if defined(PETSC_USE_CTABLE)
 85:   PetscTableDestroy(&gid1_lid1);
 86: #else
 87:   PetscBTDestroy(&table);
 88:   PetscFree(nidx);
 89: #endif
 90:   return 0;
 91: }

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

104:   for (i=0; i<imax; i++) {
105:     ISSorted(is_in[i],&flg);
107:   }

109: #if defined(PETSC_USE_CTABLE)
110:   /* Now check max size */
111:   for (i=0,maxsz=0; i<imax; i++) {
112:     ISGetLocalSize(is_in[i],&len);
114:     len = len/bs; /* The reduced index size */
115:     if (len > maxsz) maxsz = len;
116:   }
117:   PetscMalloc1(maxsz,&nidx);
118: #else
119:   PetscMalloc1(Nbs,&nidx);
120: #endif
121:   /* Now check if the indices are in block order */
122:   for (i=0; i<imax; i++) {
123:     ISGetLocalSize(is_in[i],&len);

125:     /* special case where IS is already block IS of the correct size */
126:     PetscObjectTypeCompare((PetscObject)is_in[i],ISBLOCK,&isblock);
127:     if (isblock) {
128:       ISBlockGetLocalSize(is_in[i],&bbs);
129:       if (bs == bbs) {
130:         len  = len/bs;
131:         ISBlockGetIndices(is_in[i],&idx);
132:         ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,(is_out+i));
133:         ISBlockRestoreIndices(is_in[i],&idx);
134:         continue;
135:       }
136:     }
137:     ISGetIndices(is_in[i],&idx);

140:     len       = len/bs; /* The reduced index size */
141:     idx_local = idx;
142:     for (j=0; j<len; j++) {
143:       val = idx_local[0];
145:       for (k=0; k<bs; k++) {
147:       }
148:       nidx[j]    = val/bs;
149:       idx_local += bs;
150:     }
151:     ISRestoreIndices(is_in[i],&idx);
152:     ISCreateGeneral(PetscObjectComm((PetscObject)is_in[i]),len,nidx,PETSC_COPY_VALUES,(is_out+i));
153:   }
154:   PetscFree(nidx);
155:   return 0;
156: }

158: /*@C
159:    ISExpandIndicesGeneral - convert the indices into non-block indices

161:    Input Parameters:
162: +    n - the length of the index set (not being used)
163: .    nkeys - expected number of keys when PETSC_USE_CTABLE (not being used)
164: .    bs - the size of block
165: .    imax - the number of index sets
166: -    is_in - the blocked array of index sets

168:    Output Parameter:
169: .    is_out - the non-blocked new index set

171:    Level: intermediate

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

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

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