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