Actual source code: isblock.c
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 of an array of `IS` into an array of `ISGENERAL` of block indices
9: Input Parameters:
10: + n - maximum possible length of the index set
11: . nkeys - expected number of keys when using `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, as `ISGENERAL`, not as `ISBLOCK`
19: Level: intermediate
21: .seealso: [](sec_scatter), `IS`, `ISGENERAL`, `ISExpandIndicesGeneral()`
22: @*/
23: PetscErrorCode ISCompressIndicesGeneral(PetscInt n, PetscInt nkeys, PetscInt bs, PetscInt imax, const IS is_in[], IS is_out[])
24: {
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
36: Nbs = n / bs;
37: #if defined(PETSC_USE_CTABLE)
38: Nkbs = nkeys / bs;
39: PetscTableCreate(Nkbs, Nbs, &gid1_lid1);
40: #else
41: PetscMalloc1(Nbs, &nidx);
42: PetscBTCreate(Nbs, &table);
43: #endif
44: for (i = 0; i < imax; i++) {
45: isz = 0;
46: #if defined(PETSC_USE_CTABLE)
47: PetscTableRemoveAll(gid1_lid1);
48: #else
49: PetscBTMemzero(Nbs, table);
50: #endif
51: ISGetIndices(is_in[i], &idx);
52: ISGetLocalSize(is_in[i], &len);
53: for (j = 0; j < len; j++) {
54: ival = idx[j] / bs; /* convert the indices into block indices */
55: #if defined(PETSC_USE_CTABLE)
56: PetscTableFind(gid1_lid1, ival + 1, &tt);
57: if (!tt) {
58: PetscTableAdd(gid1_lid1, ival + 1, isz + 1, INSERT_VALUES);
59: isz++;
60: }
61: #else
63: if (!PetscBTLookupSet(table, ival)) nidx[isz++] = ival;
64: #endif
65: }
66: ISRestoreIndices(is_in[i], &idx);
68: #if defined(PETSC_USE_CTABLE)
69: PetscMalloc1(isz, &nidx);
70: PetscTableGetHeadPosition(gid1_lid1, &tpos);
71: j = 0;
72: while (tpos) {
73: PetscTableGetNext(gid1_lid1, &tpos, &gid1, &tt);
75: nidx[tt] = gid1 - 1;
76: j++;
77: }
79: ISCreateGeneral(PetscObjectComm((PetscObject)is_in[i]), isz, nidx, PETSC_OWN_POINTER, (is_out + i));
80: #else
81: ISCreateGeneral(PetscObjectComm((PetscObject)is_in[i]), isz, nidx, PETSC_COPY_VALUES, (is_out + i));
82: #endif
83: }
84: #if defined(PETSC_USE_CTABLE)
85: PetscTableDestroy(&gid1_lid1);
86: #else
87: PetscBTDestroy(&table);
88: PetscFree(nidx);
89: #endif
90: return 0;
91: }
93: PetscErrorCode ISCompressIndicesSorted(PetscInt n, PetscInt bs, PetscInt imax, const IS is_in[], IS is_out[])
94: {
95: PetscInt i, j, k, val, len, *nidx, bbs;
96: const PetscInt *idx, *idx_local;
97: PetscBool flg, isblock;
98: #if defined(PETSC_USE_CTABLE)
99: PetscInt maxsz;
100: #else
101: PetscInt Nbs = n / bs;
102: #endif
104: for (i = 0; i < imax; i++) {
105: ISSorted(is_in[i], &flg);
107: }
109: #if defined(PETSC_USE_CTABLE)
110: /* Now check max size */
111: for (i = 0, maxsz = 0; i < imax; i++) {
112: ISGetLocalSize(is_in[i], &len);
114: len = len / bs; /* The reduced index size */
115: if (len > maxsz) maxsz = len;
116: }
117: PetscMalloc1(maxsz, &nidx);
118: #else
119: PetscMalloc1(Nbs, &nidx);
120: #endif
121: /* Now check if the indices are in block order */
122: for (i = 0; i < imax; i++) {
123: ISGetLocalSize(is_in[i], &len);
125: /* special case where IS is already block IS of the correct size */
126: PetscObjectTypeCompare((PetscObject)is_in[i], ISBLOCK, &isblock);
127: if (isblock) {
128: ISBlockGetLocalSize(is_in[i], &bbs);
129: if (bs == bbs) {
130: len = len / bs;
131: ISBlockGetIndices(is_in[i], &idx);
132: ISCreateGeneral(PETSC_COMM_SELF, len, idx, PETSC_COPY_VALUES, (is_out + i));
133: ISBlockRestoreIndices(is_in[i], &idx);
134: continue;
135: }
136: }
137: ISGetIndices(is_in[i], &idx);
140: len = len / bs; /* The reduced index size */
141: idx_local = idx;
142: for (j = 0; j < len; j++) {
143: val = idx_local[0];
146: nidx[j] = val / bs;
147: idx_local += bs;
148: }
149: ISRestoreIndices(is_in[i], &idx);
150: ISCreateGeneral(PetscObjectComm((PetscObject)is_in[i]), len, nidx, PETSC_COPY_VALUES, (is_out + i));
151: }
152: PetscFree(nidx);
153: return 0;
154: }
156: /*@C
157: ISExpandIndicesGeneral - convert the indices of an array `IS` into non-block indices in an array of `ISGENERAL`
159: Input Parameters:
160: + n - the length of the index set (not being used)
161: . nkeys - expected number of keys when `PETSC_USE_CTABLE` is used
162: . bs - the size of block
163: . imax - the number of index sets
164: - is_in - the blocked array of index sets
166: Output Parameter:
167: . is_out - the non-blocked new index set, as `ISGENERAL`
169: Level: intermediate
171: .seealso: [](sec_scatter), `IS`, `ISGENERAL`, `ISCompressIndicesGeneral()`
172: @*/
173: PetscErrorCode ISExpandIndicesGeneral(PetscInt n, PetscInt nkeys, PetscInt bs, PetscInt imax, const IS is_in[], IS is_out[])
174: {
175: PetscInt len, i, j, k, *nidx;
176: const PetscInt *idx;
177: PetscInt maxsz;
179: /* Check max size of is_in[] */
180: maxsz = 0;
181: for (i = 0; i < imax; i++) {
182: ISGetLocalSize(is_in[i], &len);
183: if (len > maxsz) maxsz = len;
184: }
185: PetscMalloc1(maxsz * bs, &nidx);
187: for (i = 0; i < imax; i++) {
188: ISGetLocalSize(is_in[i], &len);
189: ISGetIndices(is_in[i], &idx);
190: for (j = 0; j < len; ++j) {
191: for (k = 0; k < bs; k++) nidx[j * bs + k] = idx[j] * bs + k;
192: }
193: ISRestoreIndices(is_in[i], &idx);
194: ISCreateGeneral(PetscObjectComm((PetscObject)is_in[i]), len * bs, nidx, PETSC_COPY_VALUES, is_out + i);
195: }
196: PetscFree(nidx);
197: return 0;
198: }