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 <petsc/private/hashmapi.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: Note:
22: When `is_in[i]` already has the requested block size and needs no compression,
23: `is_out[i]` may be `is_in[i]` itself with an incremented reference count.
25: .seealso: [](sec_scatter), `IS`, `ISGENERAL`, `ISExpandIndicesGeneral()`
26: @*/
27: PetscErrorCode ISCompressIndicesGeneral(PetscInt n, PetscInt nkeys, PetscInt bs, PetscInt imax, const IS is_in[], IS is_out[])
28: {
29: PetscInt isz, len, i, j, ival, bbs;
30: const PetscInt *idx;
31: PetscBool flg;
32: #if PetscDefined(USE_CTABLE)
33: PetscHMapI gid1_lid1 = NULL;
34: PetscInt tt, gid1, *nidx;
35: PetscHashIter tpos;
36: #else
37: PetscInt *nidx;
38: PetscInt Nbs;
39: PetscBT table;
40: #endif
42: PetscFunctionBegin;
43: #if PetscDefined(USE_CTABLE)
44: PetscCall(PetscHMapICreateWithSize(nkeys / bs, &gid1_lid1));
45: #else
46: Nbs = n / bs;
47: PetscCall(PetscMalloc1(Nbs, &nidx));
48: PetscCall(PetscBTCreate(Nbs, &table));
49: #endif
50: for (i = 0; i < imax; i++) {
51: PetscCall(ISGetLocalSize(is_in[i], &len));
52: PetscCall(ISGetBlockSize(is_in[i], &bbs));
53: /* special cases where IS already has the correct block size */
54: if (bs == bbs) {
55: PetscCall(PetscObjectTypeCompare((PetscObject)is_in[i], ISBLOCK, &flg));
56: if (flg) {
57: len = len / bs;
58: PetscCall(ISBlockGetIndices(is_in[i], &idx));
59: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is_in[i]), len, idx, PETSC_COPY_VALUES, is_out + i));
60: PetscCall(ISBlockRestoreIndices(is_in[i], &idx));
61: continue;
62: } else {
63: PetscCall(PetscObjectTypeCompare((PetscObject)is_in[i], ISSTRIDE, &flg));
64: if (flg) {
65: PetscInt first, nblocks;
66: PetscInt *idx;
68: PetscCall(ISStrideGetInfo(is_in[i], &first, &j));
69: if (j == 1) {
70: nblocks = len ? (first + len - 1) / bs - first / bs + 1 : 0;
71: PetscCall(PetscMalloc1(nblocks, &idx));
72: for (j = 0; j < nblocks; ++j) idx[j] = first / bs + j;
73: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is_in[i]), nblocks, idx, PETSC_OWN_POINTER, is_out + i));
74: continue;
75: }
76: } else if (bs == 1) { /* ISGENERAL */
77: is_out[i] = is_in[i];
78: PetscCall(PetscObjectReference((PetscObject)is_in[i]));
79: continue;
80: }
81: }
82: }
83: isz = 0;
84: #if PetscDefined(USE_CTABLE)
85: PetscCall(PetscHMapIClear(gid1_lid1));
86: #else
87: PetscCall(PetscBTMemzero(Nbs, table));
88: #endif
89: PetscCall(ISGetIndices(is_in[i], &idx));
90: for (j = 0; j < len; j++) {
91: ival = idx[j] / bs; /* convert the indices into block indices */
92: #if PetscDefined(USE_CTABLE)
93: PetscCall(PetscHMapIGetWithDefault(gid1_lid1, ival + 1, 0, &tt));
94: if (!tt) {
95: PetscCall(PetscHMapISet(gid1_lid1, ival + 1, isz + 1));
96: isz++;
97: }
98: #else
99: PetscCheck(ival <= Nbs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index greater than mat-dim");
100: if (!PetscBTLookupSet(table, ival)) nidx[isz++] = ival;
101: #endif
102: }
103: PetscCall(ISRestoreIndices(is_in[i], &idx));
105: #if PetscDefined(USE_CTABLE)
106: PetscCall(PetscMalloc1(isz, &nidx));
107: PetscHashIterBegin(gid1_lid1, tpos);
108: j = 0;
109: while (!PetscHashIterAtEnd(gid1_lid1, tpos)) {
110: PetscHashIterGetKey(gid1_lid1, tpos, gid1);
111: PetscHashIterGetVal(gid1_lid1, tpos, tt);
112: PetscCheck(tt-- <= isz, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index greater than array-dim");
113: nidx[tt] = gid1 - 1;
114: j++;
115: PetscHashIterNext(gid1_lid1, tpos);
116: }
117: PetscCheck(j == isz, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "table error: jj != isz");
118: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is_in[i]), isz, nidx, PETSC_OWN_POINTER, is_out + i));
119: #else
120: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is_in[i]), isz, nidx, PETSC_COPY_VALUES, is_out + i));
121: #endif
122: }
123: #if PetscDefined(USE_CTABLE)
124: PetscCall(PetscHMapIDestroy(&gid1_lid1));
125: #else
126: PetscCall(PetscBTDestroy(&table));
127: PetscCall(PetscFree(nidx));
128: #endif
129: PetscFunctionReturn(PETSC_SUCCESS);
130: }
132: /*@
133: ISExpandIndicesGeneral - convert the indices of an array `IS` into non-block indices in an array of `ISGENERAL`
135: Input Parameters:
136: + n - the length of the index set (not being used)
137: . nkeys - expected number of keys when `PETSC_USE_CTABLE` is used
138: . bs - the size of block
139: . imax - the number of index sets
140: - is_in - the blocked array of index sets, must be as large as `imax`
142: Output Parameter:
143: . is_out - the non-blocked new index set, as `ISGENERAL`, must be as large as `imax`
145: Level: intermediate
147: .seealso: [](sec_scatter), `IS`, `ISGENERAL`, `ISCompressIndicesGeneral()`
148: @*/
149: PetscErrorCode ISExpandIndicesGeneral(PetscInt n, PetscInt nkeys, PetscInt bs, PetscInt imax, const IS is_in[], IS is_out[])
150: {
151: PetscInt len, i, j, k, *nidx;
152: const PetscInt *idx;
153: PetscInt maxsz;
155: PetscFunctionBegin;
156: /* Check max size of is_in[] */
157: maxsz = 0;
158: for (i = 0; i < imax; i++) {
159: PetscCall(ISGetLocalSize(is_in[i], &len));
160: if (len > maxsz) maxsz = len;
161: }
162: PetscCall(PetscMalloc1(maxsz * bs, &nidx));
164: for (i = 0; i < imax; i++) {
165: PetscCall(ISGetLocalSize(is_in[i], &len));
166: PetscCall(ISGetIndices(is_in[i], &idx));
167: for (j = 0; j < len; ++j) {
168: for (k = 0; k < bs; k++) nidx[j * bs + k] = idx[j] * bs + k;
169: }
170: PetscCall(ISRestoreIndices(is_in[i], &idx));
171: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is_in[i]), len * bs, nidx, PETSC_COPY_VALUES, is_out + i));
172: }
173: PetscCall(PetscFree(nidx));
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }