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