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: }