Actual source code: isblock.c
petsc-3.11.4 2019-09-28
1: /* Routines to be used by MatIncreaseOverlap() for BAIJ and SBAIJ matrices */
2: #include <petscis.h>
3: #include <petscbt.h>
4: #include <petscctable.h>
6: /*@
7: ISCompressIndicesGeneral - convert the indices into block indices
8: Input Parameters:
9: + n - maximum possible length of the index set
10: . nkeys - expected number of keys when PETSC_USE_CTABLE
11: . bs - the size of block
12: . imax - the number of index sets
13: - is_in - the non-blocked array of index sets
15: Output Parameter:
16: . is_out - the blocked new index set
18: Level: intermediate
20: .seealso: ISExpandIndicesGeneral()
21: @*/
22: PetscErrorCode ISCompressIndicesGeneral(PetscInt n,PetscInt nkeys,PetscInt bs,PetscInt imax,const IS is_in[],IS is_out[])
23: {
24: PetscErrorCode ierr;
25: PetscInt isz,len,i,j,ival,Nbs;
26: const PetscInt *idx;
27: #if defined(PETSC_USE_CTABLE)
28: PetscTable gid1_lid1;
29: PetscInt tt, gid1, *nidx,Nkbs;
30: PetscTablePosition tpos;
31: #else
32: PetscInt *nidx;
33: PetscBT table;
34: #endif
37: Nbs = n/bs;
38: #if defined(PETSC_USE_CTABLE)
39: Nkbs = nkeys/bs;
40: PetscTableCreate(Nkbs,Nbs,&gid1_lid1);
41: #else
42: PetscMalloc1(Nbs,&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,INSERT_VALUES);
60: isz++;
61: }
62: #else
63: if (ival>Nbs) SETERRQ(PETSC_COMM_SELF,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);
69: #if defined(PETSC_USE_CTABLE)
70: PetscMalloc1(isz,&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_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"index greater than array-dim");
76: nidx[tt] = gid1 - 1;
77: j++;
78: }
79: if (j != isz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"table error: jj != isz");
80: ISCreateGeneral(PetscObjectComm((PetscObject)is_in[i]),isz,nidx,PETSC_OWN_POINTER,(is_out+i));
81: #else
82: ISCreateGeneral(PetscObjectComm((PetscObject)is_in[i]),isz,nidx,PETSC_COPY_VALUES,(is_out+i));
83: #endif
84: }
85: #if defined(PETSC_USE_CTABLE)
86: PetscTableDestroy(&gid1_lid1);
87: #else
88: PetscBTDestroy(&table);
89: PetscFree(nidx);
90: #endif
91: return(0);
92: }
94: PetscErrorCode ISCompressIndicesSorted(PetscInt n,PetscInt bs,PetscInt imax,const IS is_in[],IS is_out[])
95: {
97: PetscInt i,j,k,val,len,*nidx,bbs;
98: const PetscInt *idx,*idx_local;
99: PetscBool flg,isblock;
100: #if defined(PETSC_USE_CTABLE)
101: PetscInt maxsz;
102: #else
103: PetscInt Nbs=n/bs;
104: #endif
107: for (i=0; i<imax; i++) {
108: ISSorted(is_in[i],&flg);
109: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Indices are not sorted");
110: }
112: #if defined(PETSC_USE_CTABLE)
113: /* Now check max size */
114: for (i=0,maxsz=0; i<imax; i++) {
115: ISGetLocalSize(is_in[i],&len);
116: if (len%bs !=0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Indices are not block ordered");
117: len = len/bs; /* The reduced index size */
118: if (len > maxsz) maxsz = len;
119: }
120: PetscMalloc1(maxsz,&nidx);
121: #else
122: PetscMalloc1(Nbs,&nidx);
123: #endif
124: /* Now check if the indices are in block order */
125: for (i=0; i<imax; i++) {
126: ISGetLocalSize(is_in[i],&len);
128: /* special case where IS is already block IS of the correct size */
129: PetscObjectTypeCompare((PetscObject)is_in[i],ISBLOCK,&isblock);
130: if (isblock) {
131: ISBlockGetLocalSize(is_in[i],&bbs);
132: if (bs == bbs) {
133: len = len/bs;
134: ISBlockGetIndices(is_in[i],&idx);
135: ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,(is_out+i));
136: ISBlockRestoreIndices(is_in[i],&idx);
137: continue;
138: }
139: }
140: ISGetIndices(is_in[i],&idx);
141: if (len%bs !=0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Indices are not block ordered");
143: len = len/bs; /* The reduced index size */
144: idx_local = idx;
145: for (j=0; j<len; j++) {
146: val = idx_local[0];
147: if (val%bs != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Indices are not block ordered");
148: for (k=0; k<bs; k++) {
149: if (val+k != idx_local[k]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Indices are not block ordered");
150: }
151: nidx[j] = val/bs;
152: idx_local += bs;
153: }
154: ISRestoreIndices(is_in[i],&idx);
155: ISCreateGeneral(PetscObjectComm((PetscObject)is_in[i]),len,nidx,PETSC_COPY_VALUES,(is_out+i));
156: }
157: PetscFree(nidx);
158: return(0);
159: }
161: /*@C
162: ISExpandIndicesGeneral - convert the indices into non-block indices
163: Input Parameters:
164: + n - the length of the index set (not being used)
165: . nkeys - expected number of keys when PETSC_USE_CTABLE (not being used)
166: . bs - the size of block
167: . imax - the number of index sets
168: - is_in - the blocked array of index sets
170: Output Parameter:
171: . is_out - the non-blocked new index set
173: Level: intermediate
175: .seealso: ISCompressIndicesGeneral()
176: @*/
177: PetscErrorCode ISExpandIndicesGeneral(PetscInt n,PetscInt nkeys,PetscInt bs,PetscInt imax,const IS is_in[],IS is_out[])
178: {
180: PetscInt len,i,j,k,*nidx;
181: const PetscInt *idx;
182: PetscInt maxsz;
185: /* Check max size of is_in[] */
186: maxsz = 0;
187: for (i=0; i<imax; i++) {
188: ISGetLocalSize(is_in[i],&len);
189: if (len > maxsz) maxsz = len;
190: }
191: PetscMalloc1(maxsz*bs,&nidx);
193: for (i=0; i<imax; i++) {
194: ISGetLocalSize(is_in[i],&len);
195: ISGetIndices(is_in[i],&idx);
196: for (j=0; j<len ; ++j) {
197: for (k=0; k<bs; k++) nidx[j*bs+k] = idx[j]*bs+k;
198: }
199: ISRestoreIndices(is_in[i],&idx);
200: ISCreateGeneral(PetscObjectComm((PetscObject)is_in[i]),len*bs,nidx,PETSC_COPY_VALUES,is_out+i);
201: }
202: PetscFree(nidx);
203: return(0);
204: }