Actual source code: isblock.c

petsc-3.5.2 2014-09-08
Report Typos and Errors
  2: /* Routines to be used by MatIncreaseOverlap() for BAIJ and SBAIJ matrices */
  3: #include <petscis.h>                       /*I "petscis.h"  I*/
  4: #include <petscbt.h>
  5: #include <petscctable.h>

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

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

 21:    Level: intermediate

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

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

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

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

112:   for (i=0; i<imax; i++) {
113:     ISSorted(is_in[i],&flg);
114:     if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Indices are not sorted");
115:   }

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

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

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

168: /*@C
169:    ISExpandIndicesGeneral - convert the indices into non-block indices
170:    Input Parameters:
171: +  n - the length of the index set   (not being used)
172: .  nkeys - expected number of keys when PETSC_USE_CTABLE (not being used)
173: .  bs - the size of block
174: .  imax - the number of index sets
175: -  is_in - the blocked array of index sets

177:    Output Parameter:
178: .  is_out - the non-blocked new index set

180:    Level: intermediate

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

192:   /* Check max size of is_in[] */
193:   maxsz=0;
194:   for (i=0; i<imax; i++) {
195:     ISGetLocalSize(is_in[i],&len);
196:     if (len > maxsz) maxsz = len;
197:   }
198:   PetscMalloc1(maxsz*bs,&nidx);

200:   for (i=0; i<imax; i++) {
201:     ISGetLocalSize(is_in[i],&len);
202:     ISGetIndices(is_in[i],&idx);
203:     for (j=0; j<len ; ++j) {
204:       for (k=0; k<bs; k++) nidx[j*bs+k] = idx[j]*bs+k;
205:     }
206:     ISRestoreIndices(is_in[i],&idx);
207:     ISCreateGeneral(PETSC_COMM_SELF,len*bs,nidx,PETSC_COPY_VALUES,is_out+i);
208:   }
209:   PetscFree(nidx);
210:   return(0);
211: }