Actual source code: isblock.c

petsc-3.8.2 2017-11-09
Report Typos and Errors

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

  7: /*@
  8:    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: {
 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(PETSC_COMM_SELF,isz,nidx,PETSC_OWN_POINTER,(is_out+i));
 82: #else
 83:     ISCreateGeneral(PETSC_COMM_SELF,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(PETSC_COMM_SELF,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
164:    Input Parameters:
165: +  n - the length of the index set   (not being used)
166: .  nkeys - expected number of keys when PETSC_USE_CTABLE (not being used)
167: .  bs - the size of block
168: .  imax - the number of index sets
169: -  is_in - the blocked array of index sets

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

174:    Level: intermediate

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

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

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