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,ival,Nbs;
28: const PetscInt *idx;
29: #if defined (PETSC_USE_CTABLE)
30: PetscTable gid1_lid1;
31: PetscInt tt, gid1, *nidx;
32: PetscTablePosition tpos;
33: #else
34: PetscInt *nidx;
35: PetscBT table;
36: #endif
39: Nbs =n/bs;
40: #if defined (PETSC_USE_CTABLE)
41: PetscTableCreate(Nbs,&gid1_lid1);
42: #else
43: PetscMalloc(Nbs*sizeof(PetscInt),&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);
61: isz++;
62: }
63: #else
64: if (ival>Nbs) SETERRQ(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);
69: #if defined (PETSC_USE_CTABLE)
70: PetscMalloc(isz*sizeof(PetscInt),&nidx);
71: PetscTableGetHeadPosition(gid1_lid1,&tpos);
72: j = 0;
73: while (tpos) {
74: PetscTableGetNext(gid1_lid1,&tpos,&gid1,&tt);
75: if (tt-- > isz) { SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"index greater than array-dim"); }
76: nidx[tt] = gid1 - 1;
77: j++;
78: }
79: if (j != isz) { SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"table error: jj != isz"); }
80: ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is_out+i));
81: PetscFree(nidx);
82: #else
83: ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(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: }
97: PetscErrorCode ISCompressIndicesSorted(PetscInt n,PetscInt bs,PetscInt imax,const IS is_in[],IS is_out[])
98: {
100: PetscInt i,j,k,val,len,*nidx,bbs;
101: const PetscInt *idx,*idx_local;
102: PetscTruth flg,isblock;
103: #if defined (PETSC_USE_CTABLE)
104: PetscInt maxsz;
105: #else
106: PetscInt Nbs=n/bs;
107: #endif
110: for (i=0; i<imax; i++) {
111: ISSorted(is_in[i],&flg);
112: if (!flg) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Indices are not sorted");
113: }
115: #if defined (PETSC_USE_CTABLE)
116: /* Now check max size */
117: for (i=0,maxsz=0; i<imax; i++) {
118: ISGetLocalSize(is_in[i],&len);
119: if (len%bs !=0) SETERRQ(PETSC_ERR_ARG_INCOMP,"Indices are not block ordered");
120: len = len/bs; /* The reduced index size */
121: if (len > maxsz) maxsz = len;
122: }
123: PetscMalloc(maxsz*sizeof(PetscInt),&nidx);
124: #else
125: PetscMalloc(Nbs*sizeof(PetscInt),&nidx);
126: #endif
127: /* Now check if the indices are in block order */
128: for (i=0; i<imax; i++) {
129: ISGetLocalSize(is_in[i],&len);
131: /* special case where IS is already block IS of the correct size */
132: ISBlock(is_in[i],&isblock);
133: if (isblock) {
134: ISBlockGetLocalSize(is_in[i],&bbs);
135: if (bs == bbs) {
136: len = len/bs;
137: ISBlockGetIndices(is_in[i],&idx);
138: for (j=0; j<len; j++) nidx[j] = idx[j]/bs;
139: ISBlockRestoreIndices(is_in[i],&idx);
140: ISCreateGeneral(PETSC_COMM_SELF,len,nidx,(is_out+i));
141: continue;
142: }
143: }
144: ISGetIndices(is_in[i],&idx);
145: if (len%bs !=0) SETERRQ(PETSC_ERR_ARG_INCOMP,"Indices are not block ordered");
147: len = len/bs; /* The reduced index size */
148: idx_local = idx;
149: for (j=0; j<len ; j++) {
150: val = idx_local[0];
151: if (val%bs != 0) SETERRQ(PETSC_ERR_ARG_INCOMP,"Indices are not block ordered");
152: for (k=0; k<bs; k++) {
153: if (val+k != idx_local[k]) SETERRQ(PETSC_ERR_ARG_INCOMP,"Indices are not block ordered");
154: }
155: nidx[j] = val/bs;
156: idx_local +=bs;
157: }
158: ISRestoreIndices(is_in[i],&idx);
159: ISCreateGeneral(PETSC_COMM_SELF,len,nidx,(is_out+i));
160: }
161: PetscFree(nidx);
163: return(0);
164: }
168: PetscErrorCode ISExpandIndicesGeneral(PetscInt n,PetscInt bs,PetscInt imax,const IS is_in[],IS is_out[])
169: {
171: PetscInt len,i,j,k,*nidx;
172: const PetscInt *idx;
173: #if defined (PETSC_USE_CTABLE)
174: PetscInt maxsz;
175: #else
176: PetscInt Nbs;
177: #endif
180: #if defined (PETSC_USE_CTABLE)
181: /* Now check max size */
182: for (i=0,maxsz=0; i<imax; i++) {
183: ISGetIndices(is_in[i],&idx);
184: ISGetLocalSize(is_in[i],&len);
185: if (len*bs > maxsz) maxsz = len*bs;
186: }
187: PetscMalloc(maxsz*sizeof(PetscInt),&nidx);
188: #else
189: Nbs = n/bs;
190: PetscMalloc(Nbs*bs*sizeof(PetscInt),&nidx);
191: #endif
193: for (i=0; i<imax; i++) {
194: ISGetIndices(is_in[i],&idx);
195: ISGetLocalSize(is_in[i],&len);
196: for (j=0; j<len ; ++j){
197: for (k=0; k<bs; k++)
198: nidx[j*bs+k] = idx[j]*bs+k;
199: }
200: ISRestoreIndices(is_in[i],&idx);
201: ISCreateGeneral(PETSC_COMM_SELF,len*bs,nidx,is_out+i);
202: }
203: PetscFree(nidx);
204: return(0);
205: }