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:   PetscErrorCode     ierr;
 26:   PetscInt           isz,len,i,j,ival,Nbs;
 27:   const PetscInt     *idx;
 28: #if defined(PETSC_USE_CTABLE)
 29:   PetscTable         gid1_lid1;
 30:   PetscInt           tt, gid1, *nidx,Nkbs;
 31:   PetscTablePosition tpos;
 32: #else
 33:   PetscInt           *nidx;
 34:   PetscBT            table;
 35: #endif

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

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

 95: PetscErrorCode  ISCompressIndicesSorted(PetscInt n,PetscInt bs,PetscInt imax,const IS is_in[],IS is_out[])
 96: {
 98:   PetscInt       i,j,k,val,len,*nidx,bbs;
 99:   const PetscInt *idx,*idx_local;
100:   PetscBool      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_COMM_SELF,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_COMM_SELF,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:   PetscMalloc1(maxsz,&nidx);
122: #else
123:   PetscMalloc1(Nbs,&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:     PetscObjectTypeCompare((PetscObject)is_in[i],ISBLOCK,&isblock);
131:     if (isblock) {
132:       ISBlockGetLocalSize(is_in[i],&bbs);
133:       if (bs == bbs) {
134:         len  = len/bs;
135:         ISBlockGetIndices(is_in[i],&idx);
136:         ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,(is_out+i));
137:         ISBlockRestoreIndices(is_in[i],&idx);
138:         continue;
139:       }
140:     }
141:     ISGetIndices(is_in[i],&idx);
142:     if (len%bs !=0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Indices are not block ordered");

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

162: /*@C
163:    ISExpandIndicesGeneral - convert the indices into non-block indices

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

172:    Output Parameter:
173: .    is_out - the non-blocked new index set

175:    Level: intermediate

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

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

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