Actual source code: isblock.c
petsc-3.3-p7 2013-05-11
2: /* Routines to be used by MatIncreaseOverlap() for BAIJ and SBAIJ matrices */
3: #include <petscis.h>
4: #include <petscbt.h>
9: /*@C
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: {
27: PetscErrorCode ierr;
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: PetscMalloc(Nbs*sizeof(PetscInt),&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);
71:
72: #if defined (PETSC_USE_CTABLE)
73: PetscMalloc(isz*sizeof(PetscInt),&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: PetscMalloc(maxsz*sizeof(PetscInt),&nidx);
126: #else
127: PetscMalloc(Nbs*sizeof(PetscInt),&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);
164: return(0);
165: }
169: /*@C
170: ISExpandIndicesGeneral - convert the indices into non-block indices
171: Input Parameters:
172: + n - the length of the index set (not being used)
173: . nkeys - expected number of keys when PETSC_USE_CTABLE (not being used)
174: . bs - the size of block
175: . imax - the number of index sets
176: - is_in - the blocked array of index sets
178: Output Parameter:
179: . is_out - the non-blocked new index set
181: Level: intermediate
183: .seealso: ISCompressIndicesGeneral()
184: @*/
185: PetscErrorCode ISExpandIndicesGeneral(PetscInt n,PetscInt nkeys,PetscInt bs,PetscInt imax,const IS is_in[],IS is_out[])
186: {
188: PetscInt len,i,j,k,*nidx;
189: const PetscInt *idx;
190: PetscInt maxsz;
193: /* Check max size of is_in[] */
194: maxsz=0;
195: for (i=0; i<imax; i++) {
196: ISGetLocalSize(is_in[i],&len);
197: if (len > maxsz) maxsz = len;
198: }
199: PetscMalloc(maxsz*bs*sizeof(PetscInt),&nidx);
201: for (i=0; i<imax; i++) {
202: ISGetLocalSize(is_in[i],&len);
203: ISGetIndices(is_in[i],&idx);
204: for (j=0; j<len ; ++j){
205: for (k=0; k<bs; k++) nidx[j*bs+k] = idx[j]*bs+k;
206: }
207: ISRestoreIndices(is_in[i],&idx);
208: ISCreateGeneral(PETSC_COMM_SELF,len*bs,nidx,PETSC_COPY_VALUES,is_out+i);
209: }
210: PetscFree(nidx);
211: return(0);
212: }