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