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  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  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  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: }