Actual source code: isblock.c
petsc-3.8.4 2018-03-24
2: /* Routines to be used by MatIncreaseOverlap() for BAIJ and SBAIJ matrices */
3: #include <petscis.h>
4: #include <petscbt.h>
5: #include <petscctable.h>
7: /*@
8: 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: {
26: PetscInt isz,len,i,j,ival,Nbs;
27: const PetscInt *idx;
28: #if defined(PETSC_USE_CTABLE)
29: PetscTable gid1_lid1;
30: PetscInt tt, gid1, *nidx,Nkbs;
31: PetscTablePosition tpos;
32: #else
33: PetscInt *nidx;
34: PetscBT table;
35: #endif
38: Nbs =n/bs;
39: #if defined(PETSC_USE_CTABLE)
40: Nkbs = nkeys/bs;
41: PetscTableCreate(Nkbs,Nbs,&gid1_lid1);
42: #else
43: PetscMalloc1(Nbs,&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,INSERT_VALUES);
61: isz++;
62: }
63: #else
64: if (ival>Nbs) SETERRQ(PETSC_COMM_SELF,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);
70: #if defined(PETSC_USE_CTABLE)
71: PetscMalloc1(isz,&nidx);
72: PetscTableGetHeadPosition(gid1_lid1,&tpos);
73: j = 0;
74: while (tpos) {
75: PetscTableGetNext(gid1_lid1,&tpos,&gid1,&tt);
76: if (tt-- > isz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"index greater than array-dim");
77: nidx[tt] = gid1 - 1;
78: j++;
79: }
80: if (j != isz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"table error: jj != isz");
81: ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_OWN_POINTER,(is_out+i));
82: #else
83: ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(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: }
95: PetscErrorCode ISCompressIndicesSorted(PetscInt n,PetscInt bs,PetscInt imax,const IS is_in[],IS is_out[])
96: {
98: PetscInt i,j,k,val,len,*nidx,bbs;
99: const PetscInt *idx,*idx_local;
100: PetscBool 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_COMM_SELF,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_COMM_SELF,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: PetscMalloc1(maxsz,&nidx);
122: #else
123: PetscMalloc1(Nbs,&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: PetscObjectTypeCompare((PetscObject)is_in[i],ISBLOCK,&isblock);
131: if (isblock) {
132: ISBlockGetLocalSize(is_in[i],&bbs);
133: if (bs == bbs) {
134: len = len/bs;
135: ISBlockGetIndices(is_in[i],&idx);
136: ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,(is_out+i));
137: ISBlockRestoreIndices(is_in[i],&idx);
138: continue;
139: }
140: }
141: ISGetIndices(is_in[i],&idx);
142: if (len%bs !=0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Indices are not block ordered");
144: len = len/bs; /* The reduced index size */
145: idx_local = idx;
146: for (j=0; j<len; j++) {
147: val = idx_local[0];
148: if (val%bs != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Indices are not block ordered");
149: for (k=0; k<bs; k++) {
150: if (val+k != idx_local[k]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Indices are not block ordered");
151: }
152: nidx[j] = val/bs;
153: idx_local += bs;
154: }
155: ISRestoreIndices(is_in[i],&idx);
156: ISCreateGeneral(PETSC_COMM_SELF,len,nidx,PETSC_COPY_VALUES,(is_out+i));
157: }
158: PetscFree(nidx);
159: return(0);
160: }
162: /*@C
163: ISExpandIndicesGeneral - convert the indices into non-block indices
164: Input Parameters:
165: + n - the length of the index set (not being used)
166: . nkeys - expected number of keys when PETSC_USE_CTABLE (not being used)
167: . bs - the size of block
168: . imax - the number of index sets
169: - is_in - the blocked array of index sets
171: Output Parameter:
172: . is_out - the non-blocked new index set
174: Level: intermediate
176: .seealso: ISCompressIndicesGeneral()
177: @*/
178: PetscErrorCode ISExpandIndicesGeneral(PetscInt n,PetscInt nkeys,PetscInt bs,PetscInt imax,const IS is_in[],IS is_out[])
179: {
181: PetscInt len,i,j,k,*nidx;
182: const PetscInt *idx;
183: PetscInt maxsz;
186: /* Check max size of is_in[] */
187: maxsz=0;
188: for (i=0; i<imax; i++) {
189: ISGetLocalSize(is_in[i],&len);
190: if (len > maxsz) maxsz = len;
191: }
192: PetscMalloc1(maxsz*bs,&nidx);
194: for (i=0; i<imax; i++) {
195: ISGetLocalSize(is_in[i],&len);
196: ISGetIndices(is_in[i],&idx);
197: for (j=0; j<len ; ++j) {
198: for (k=0; k<bs; k++) nidx[j*bs+k] = idx[j]*bs+k;
199: }
200: ISRestoreIndices(is_in[i],&idx);
201: ISCreateGeneral(PETSC_COMM_SELF,len*bs,nidx,PETSC_COPY_VALUES,is_out+i);
202: }
203: PetscFree(nidx);
204: return(0);
205: }