Actual source code: general.c

  1: /*
  2:      Provides the functions for index sets (IS) defined by a list of integers.
  3: */
  4: #include <../src/vec/is/is/impls/general/general.h>
  5: #include <petsc/private/viewerhdf5impl.h>

  7: static PetscErrorCode ISDuplicate_General(IS is, IS *newIS)
  8: {
  9:   IS_General *sub = (IS_General *)is->data;
 10:   PetscInt    n, bs;

 12:   PetscFunctionBegin;
 13:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
 14:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is), n, sub->idx, PETSC_COPY_VALUES, newIS));
 15:   PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
 16:   PetscCall(PetscLayoutSetBlockSize((*newIS)->map, bs));
 17:   PetscFunctionReturn(PETSC_SUCCESS);
 18: }

 20: static PetscErrorCode ISDestroy_General(IS is)
 21: {
 22:   IS_General *is_general = (IS_General *)is->data;

 24:   PetscFunctionBegin;
 25:   if (is_general->allocated) PetscCall(PetscFree(is_general->idx));
 26:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISGeneralSetIndices_C", NULL));
 27:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISGeneralFilter_C", NULL));
 28:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISGeneralSetIndicesFromMask_C", NULL));
 29:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISShift_C", NULL));
 30:   PetscCall(PetscFree(is->data));
 31:   PetscFunctionReturn(PETSC_SUCCESS);
 32: }

 34: static PetscErrorCode ISCopy_General(IS is, IS isy)
 35: {
 36:   IS_General *is_general = (IS_General *)is->data, *isy_general = (IS_General *)isy->data;
 37:   PetscInt    n;

 39:   PetscFunctionBegin;
 40:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
 41:   PetscCall(PetscArraycpy(isy_general->idx, is_general->idx, n));
 42:   PetscFunctionReturn(PETSC_SUCCESS);
 43: }

 45: static PetscErrorCode ISShift_General(IS is, PetscInt shift, IS isy)
 46: {
 47:   IS_General *is_general = (IS_General *)is->data, *isy_general = (IS_General *)isy->data;
 48:   PetscInt    i, n;

 50:   PetscFunctionBegin;
 51:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
 52:   for (i = 0; i < n; i++) isy_general->idx[i] = is_general->idx[i] + shift;
 53:   PetscFunctionReturn(PETSC_SUCCESS);
 54: }

 56: static PetscErrorCode ISOnComm_General(IS is, MPI_Comm comm, PetscCopyMode mode, IS *newis)
 57: {
 58:   IS_General *sub = (IS_General *)is->data;
 59:   PetscInt    n;

 61:   PetscFunctionBegin;
 62:   PetscCheck(mode != PETSC_OWN_POINTER, comm, PETSC_ERR_ARG_WRONG, "Cannot use PETSC_OWN_POINTER");
 63:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
 64:   PetscCall(ISCreateGeneral(comm, n, sub->idx, mode, newis));
 65:   PetscFunctionReturn(PETSC_SUCCESS);
 66: }

 68: static PetscErrorCode ISSetBlockSize_General(IS is, PetscInt bs)
 69: {
 70:   PetscFunctionBegin;
 71:   PetscCall(PetscLayoutSetBlockSize(is->map, bs));
 72:   PetscFunctionReturn(PETSC_SUCCESS);
 73: }

 75: static PetscErrorCode ISContiguousLocal_General(IS is, PetscInt gstart, PetscInt gend, PetscInt *start, PetscBool *contig)
 76: {
 77:   IS_General *sub = (IS_General *)is->data;
 78:   PetscInt    n, i, p;

 80:   PetscFunctionBegin;
 81:   *start  = 0;
 82:   *contig = PETSC_TRUE;
 83:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
 84:   if (!n) PetscFunctionReturn(PETSC_SUCCESS);
 85:   p = sub->idx[0];
 86:   if (p < gstart) goto nomatch;
 87:   *start = p - gstart;
 88:   if (n > gend - p) goto nomatch;
 89:   for (i = 1; i < n; i++, p++) {
 90:     if (sub->idx[i] != p + 1) goto nomatch;
 91:   }
 92:   PetscFunctionReturn(PETSC_SUCCESS);
 93: nomatch:
 94:   *start  = -1;
 95:   *contig = PETSC_FALSE;
 96:   PetscFunctionReturn(PETSC_SUCCESS);
 97: }

 99: static PetscErrorCode ISLocate_General(IS is, PetscInt key, PetscInt *location)
100: {
101:   IS_General *sub = (IS_General *)is->data;
102:   PetscInt    numIdx, i;
103:   PetscBool   sorted;

105:   PetscFunctionBegin;
106:   PetscCall(PetscLayoutGetLocalSize(is->map, &numIdx));
107:   PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, &sorted));
108:   if (sorted) PetscCall(PetscFindInt(key, numIdx, sub->idx, location));
109:   else {
110:     const PetscInt *idx = sub->idx;

112:     *location = -1;
113:     for (i = 0; i < numIdx; i++) {
114:       if (idx[i] == key) {
115:         *location = i;
116:         PetscFunctionReturn(PETSC_SUCCESS);
117:       }
118:     }
119:   }
120:   PetscFunctionReturn(PETSC_SUCCESS);
121: }

123: static PetscErrorCode ISGetIndices_General(IS in, const PetscInt *idx[])
124: {
125:   IS_General *sub = (IS_General *)in->data;

127:   PetscFunctionBegin;
128:   *idx = sub->idx;
129:   PetscFunctionReturn(PETSC_SUCCESS);
130: }

132: static PetscErrorCode ISRestoreIndices_General(IS in, const PetscInt *idx[])
133: {
134:   IS_General *sub = (IS_General *)in->data;

136:   PetscFunctionBegin;
137:   /* F90Array1dCreate() inside ISRestoreArrayF90() does not keep array when zero length array */
138:   PetscCheck(in->map->n <= 0 || *idx == sub->idx, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Must restore with value from ISGetIndices()");
139:   PetscFunctionReturn(PETSC_SUCCESS);
140: }

142: static PetscErrorCode ISInvertPermutation_General(IS is, PetscInt nlocal, IS *isout)
143: {
144:   IS_General     *sub = (IS_General *)is->data;
145:   PetscInt        i, *ii, n, nstart;
146:   const PetscInt *idx = sub->idx;
147:   PetscMPIInt     size;
148:   IS              istmp, nistmp;

150:   PetscFunctionBegin;
151:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
152:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is), &size));
153:   if (size == 1) {
154:     PetscCall(PetscMalloc1(n, &ii));
155:     for (i = 0; i < n; i++) ii[idx[i]] = i;
156:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, ii, PETSC_OWN_POINTER, isout));
157:     PetscCall(ISSetPermutation(*isout));
158:   } else {
159:     /* crude, nonscalable get entire IS on each processor */
160:     PetscCall(ISAllGather(is, &istmp));
161:     PetscCall(ISSetPermutation(istmp));
162:     PetscCall(ISInvertPermutation(istmp, PETSC_DECIDE, &nistmp));
163:     PetscCall(ISDestroy(&istmp));
164:     /* get the part we need */
165:     if (nlocal == PETSC_DECIDE) nlocal = n;
166:     PetscCallMPI(MPI_Scan(&nlocal, &nstart, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)is)));
167:     if (PetscDefined(USE_DEBUG)) {
168:       PetscInt    N;
169:       PetscMPIInt rank;
170:       PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)is), &rank));
171:       PetscCall(PetscLayoutGetSize(is->map, &N));
172:       PetscCheck((rank != size - 1) || (nstart == N), PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Sum of nlocal lengths %" PetscInt_FMT " != total IS length %" PetscInt_FMT, nstart, N);
173:     }
174:     nstart -= nlocal;
175:     PetscCall(ISGetIndices(nistmp, &idx));
176:     PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is), nlocal, idx + nstart, PETSC_COPY_VALUES, isout));
177:     PetscCall(ISRestoreIndices(nistmp, &idx));
178:     PetscCall(ISDestroy(&nistmp));
179:   }
180:   PetscFunctionReturn(PETSC_SUCCESS);
181: }

183: /*
184:   FindRun_Private - Determine whether a run exists in the sequence

186:   Input Parameters;
187: + indices   - An array of indices
188: . len       - The length of `indices`
189: . off       - The offset into `indices`
190: - minRunLen - The minimum size of a run

192:   Output Parameters:
193: + runLen - The length of the run, with 0 meaning no run was found
194: . step   - The step between consecutive run values
195: - val    - The initial run value

197:   Note:
198:   We define a run as an arithmetic progression of at least `minRunLen` values, where zero is a valid step.
199: */
200: static PetscErrorCode ISFindRun_Private(const PetscInt indices[], PetscInt len, PetscInt off, PetscInt minRunLen, PetscInt *runLen, PetscInt *step, PetscInt *val)
201: {
202:   PetscInt o = off, s, l;

204:   PetscFunctionBegin;
205:   *runLen = 0;
206:   // We need at least 2 values for a run
207:   if (len - off < PetscMax(minRunLen, 2)) PetscFunctionReturn(PETSC_SUCCESS);
208:   s = indices[o + 1] - indices[o];
209:   l = 2;
210:   ++o;
211:   while (o < len - 1 && (s == indices[o + 1] - indices[o])) {
212:     ++l;
213:     ++o;
214:   }
215:   if (runLen) *runLen = l;
216:   if (step) *step = s;
217:   if (val) *val = indices[off];
218:   PetscFunctionReturn(PETSC_SUCCESS);
219: }

221: static PetscErrorCode ISGeneralCheckCompress(IS is, PetscBool *compress)
222: {
223:   const PetscInt  minRun    = 8;
224:   PetscBool       lcompress = PETSC_TRUE, test = PETSC_FALSE;
225:   const PetscInt *idx;
226:   PetscInt        n, off = 0;

228:   PetscFunctionBegin;
229:   *compress = PETSC_FALSE;
230:   PetscCall(PetscOptionsGetBool(NULL, is->hdr.prefix, "-is_view_compress", &test, NULL));
231:   if (!test) PetscFunctionReturn(PETSC_SUCCESS);
232:   PetscCall(ISGetIndices(is, &idx));
233:   PetscCall(ISGetLocalSize(is, &n));
234:   while (off < n) {
235:     PetscInt len;

237:     PetscCall(ISFindRun_Private(idx, n, off, minRun, &len, NULL, NULL));
238:     if (!len) {
239:       lcompress = PETSC_FALSE;
240:       break;
241:     }
242:     off += len;
243:   }
244:   PetscCallMPI(MPI_Allreduce(&lcompress, compress, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
245:   PetscFunctionReturn(PETSC_SUCCESS);
246: }

248: #if defined(PETSC_HAVE_HDF5)
249: // Run length encode the IS
250: static PetscErrorCode ISGeneralCompress(IS is, IS *cis)
251: {
252:   const PetscInt *idx;
253:   const PetscInt  minRun = 8;
254:   PetscInt        runs   = 0;
255:   PetscInt        off    = 0;
256:   PetscInt       *cidx, n, r;
257:   const char     *name;
258:   MPI_Comm        comm;

260:   PetscFunctionBegin;
261:   PetscCall(PetscObjectGetComm((PetscObject)is, &comm));
262:   PetscCall(ISGetIndices(is, &idx));
263:   PetscCall(ISGetLocalSize(is, &n));
264:   if (!n) runs = 0;
265:   // Count runs
266:   while (off < n) {
267:     PetscInt len;

269:     PetscCall(ISFindRun_Private(idx, n, off, minRun, &len, NULL, NULL));
270:     PetscCheck(len, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Compression failed to find a run at offset %" PetscInt_FMT, off);
271:     ++runs;
272:     off += len;
273:   }
274:   // Construct compressed set
275:   PetscCall(PetscMalloc1(runs * 3, &cidx));
276:   off = 0;
277:   r   = 0;
278:   while (off < n) {
279:     PetscCall(ISFindRun_Private(idx, n, off, minRun, &cidx[r * 3 + 0], &cidx[r * 3 + 1], &cidx[r * 3 + 2]));
280:     PetscCheck(cidx[r * 3 + 0], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Compression failed to find a run at offset %" PetscInt_FMT, off);
281:     off += cidx[r * 3 + 0];
282:     ++r;
283:   }
284:   PetscCheck(r == runs, comm, PETSC_ERR_PLIB, "The number of runs %" PetscInt_FMT " != %" PetscInt_FMT " computed chunks", runs, r);
285:   PetscCheck(off == n, comm, PETSC_ERR_PLIB, "The local size %" PetscInt_FMT " != %" PetscInt_FMT " total encoded size", n, off);
286:   PetscCall(ISCreateGeneral(comm, runs * 3, cidx, PETSC_OWN_POINTER, cis));
287:   PetscCall(PetscLayoutSetBlockSize((*cis)->map, 3));
288:   PetscCall(PetscObjectGetName((PetscObject)is, &name));
289:   PetscCall(PetscObjectSetName((PetscObject)*cis, name));
290:   PetscFunctionReturn(PETSC_SUCCESS);
291: }

293: static PetscErrorCode ISView_General_HDF5(IS is, PetscViewer viewer)
294: {
295:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5 *)viewer->data;
296:   hid_t             filespace;  /* file dataspace identifier */
297:   hid_t             chunkspace; /* chunk dataset property identifier */
298:   hid_t             dset_id;    /* dataset identifier */
299:   hid_t             memspace;   /* memory dataspace identifier */
300:   hid_t             inttype;    /* int type (H5T_NATIVE_INT or H5T_NATIVE_LLONG) */
301:   hid_t             file_id, group;
302:   hsize_t           dim, maxDims[3], dims[3], chunkDims[3], count[3], offset[3];
303:   PetscBool         timestepping;
304:   PetscInt          bs, N, n, timestep = PETSC_INT_MIN, low;
305:   hsize_t           chunksize;
306:   const PetscInt   *ind;
307:   const char       *isname;

309:   PetscFunctionBegin;
310:   PetscCall(ISGetBlockSize(is, &bs));
311:   bs = PetscMax(bs, 1); /* If N = 0, bs  = 0 as well */
312:   PetscCall(PetscViewerHDF5OpenGroup(viewer, NULL, &file_id, &group));
313:   PetscCall(PetscViewerHDF5IsTimestepping(viewer, &timestepping));
314:   if (timestepping) PetscCall(PetscViewerHDF5GetTimestep(viewer, &timestep));

316:   /* Create the dataspace for the dataset.
317:    *
318:    * dims - holds the current dimensions of the dataset
319:    *
320:    * maxDims - holds the maximum dimensions of the dataset (unlimited
321:    * for the number of time steps with the current dimensions for the
322:    * other dimensions; so only additional time steps can be added).
323:    *
324:    * chunkDims - holds the size of a single time step (required to
325:    * permit extending dataset).
326:    */
327:   dim       = 0;
328:   chunksize = 1;
329:   if (timestep >= 0) {
330:     dims[dim]      = timestep + 1;
331:     maxDims[dim]   = H5S_UNLIMITED;
332:     chunkDims[dim] = 1;
333:     ++dim;
334:   }
335:   PetscCall(ISGetSize(is, &N));
336:   PetscCall(ISGetLocalSize(is, &n));
337:   PetscCall(PetscHDF5IntCast(N / bs, dims + dim));

339:   maxDims[dim]   = dims[dim];
340:   chunkDims[dim] = PetscMax(1, dims[dim]);
341:   chunksize *= chunkDims[dim];
342:   ++dim;
343:   if (bs >= 1) {
344:     dims[dim]      = bs;
345:     maxDims[dim]   = dims[dim];
346:     chunkDims[dim] = dims[dim];
347:     chunksize *= chunkDims[dim];
348:     ++dim;
349:   }
350:   /* hdf5 chunks must be less than 4GB */
351:   if (chunksize > PETSC_HDF5_MAX_CHUNKSIZE / 64) {
352:     if (bs >= 1) {
353:       if (chunkDims[dim - 2] > (hsize_t)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE / 64))) chunkDims[dim - 2] = (hsize_t)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE / 64));
354:       if (chunkDims[dim - 1] > (hsize_t)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE / 64))) chunkDims[dim - 1] = (hsize_t)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE / 64));
355:     } else {
356:       chunkDims[dim - 1] = PETSC_HDF5_MAX_CHUNKSIZE / 64;
357:     }
358:   }
359:   PetscCallHDF5Return(filespace, H5Screate_simple, ((int)dim, dims, maxDims));

361:   #if defined(PETSC_USE_64BIT_INDICES)
362:   inttype = H5T_NATIVE_LLONG;
363:   #else
364:   inttype = H5T_NATIVE_INT;
365:   #endif

367:   /* Create the dataset with default properties and close filespace */
368:   PetscCall(PetscObjectGetName((PetscObject)is, &isname));
369:   if (!H5Lexists(group, isname, H5P_DEFAULT)) {
370:     /* Create chunk */
371:     PetscCallHDF5Return(chunkspace, H5Pcreate, (H5P_DATASET_CREATE));
372:     PetscCallHDF5(H5Pset_chunk, (chunkspace, (int)dim, chunkDims));

374:     PetscCallHDF5Return(dset_id, H5Dcreate2, (group, isname, inttype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT));
375:     PetscCallHDF5(H5Pclose, (chunkspace));
376:   } else {
377:     PetscCallHDF5Return(dset_id, H5Dopen2, (group, isname, H5P_DEFAULT));
378:     PetscCallHDF5(H5Dset_extent, (dset_id, dims));
379:   }
380:   PetscCallHDF5(H5Sclose, (filespace));

382:   /* Each process defines a dataset and writes it to the hyperslab in the file */
383:   dim = 0;
384:   if (timestep >= 0) {
385:     count[dim] = 1;
386:     ++dim;
387:   }
388:   PetscCall(PetscHDF5IntCast(n / bs, count + dim));
389:   ++dim;
390:   if (bs >= 1) {
391:     count[dim] = bs;
392:     ++dim;
393:   }
394:   if (n > 0 || H5_VERSION_GE(1, 10, 0)) {
395:     PetscCallHDF5Return(memspace, H5Screate_simple, ((int)dim, count, NULL));
396:   } else {
397:     /* Can't create dataspace with zero for any dimension, so create null dataspace. */
398:     PetscCallHDF5Return(memspace, H5Screate, (H5S_NULL));
399:   }

401:   /* Select hyperslab in the file */
402:   PetscCall(PetscLayoutGetRange(is->map, &low, NULL));
403:   dim = 0;
404:   if (timestep >= 0) {
405:     offset[dim] = timestep;
406:     ++dim;
407:   }
408:   PetscCall(PetscHDF5IntCast(low / bs, offset + dim));
409:   ++dim;
410:   if (bs >= 1) {
411:     offset[dim] = 0;
412:     ++dim;
413:   }
414:   if (n > 0 || H5_VERSION_GE(1, 10, 0)) {
415:     PetscCallHDF5Return(filespace, H5Dget_space, (dset_id));
416:     PetscCallHDF5(H5Sselect_hyperslab, (filespace, H5S_SELECT_SET, offset, NULL, count, NULL));
417:   } else {
418:     /* Create null filespace to match null memspace. */
419:     PetscCallHDF5Return(filespace, H5Screate, (H5S_NULL));
420:   }

422:   PetscCall(ISGetIndices(is, &ind));
423:   PetscCallHDF5(H5Dwrite, (dset_id, inttype, memspace, filespace, hdf5->dxpl_id, ind));
424:   PetscCallHDF5(H5Fflush, (file_id, H5F_SCOPE_GLOBAL));
425:   PetscCall(ISRestoreIndices(is, &ind));

427:   /* Close/release resources */
428:   PetscCallHDF5(H5Gclose, (group));
429:   PetscCallHDF5(H5Sclose, (filespace));
430:   PetscCallHDF5(H5Sclose, (memspace));
431:   PetscCallHDF5(H5Dclose, (dset_id));

433:   if (timestepping) PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)is, "timestepping", PETSC_BOOL, &timestepping));
434:   PetscCall(PetscInfo(is, "Wrote IS object with name %s\n", isname));
435:   PetscFunctionReturn(PETSC_SUCCESS);
436: }

438: static PetscErrorCode ISView_General_HDF5_Compressed(IS is, PetscViewer viewer)
439: {
440:   const PetscBool compressed = PETSC_TRUE;
441:   const char     *isname;
442:   IS              cis;
443:   PetscBool       timestepping;

445:   PetscFunctionBegin;
446:   PetscCall(PetscViewerHDF5IsTimestepping(viewer, &timestepping));
447:   PetscCheck(!timestepping, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_WRONG, "Timestepping not supported for compressed writing");

449:   PetscCall(ISGeneralCompress(is, &cis));
450:   PetscCall(ISView_General_HDF5(cis, viewer));
451:   PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)cis, "compressed", PETSC_BOOL, &compressed));
452:   PetscCall(ISDestroy(&cis));

454:   PetscCall(PetscObjectGetName((PetscObject)is, &isname));
455:   PetscCall(PetscInfo(is, "Wrote compressed IS object with name %s\n", isname));
456:   PetscFunctionReturn(PETSC_SUCCESS);
457: }
458: #endif

460: static PetscErrorCode ISView_General(IS is, PetscViewer viewer)
461: {
462:   IS_General *sub = (IS_General *)is->data;
463:   PetscInt   *idx = sub->idx, n;
464:   PetscBool   iascii, isbinary, ishdf5, compress;

466:   PetscFunctionBegin;
467:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
468:   PetscCall(ISGeneralCheckCompress(is, &compress));
469:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
470:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
471:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
472:   if (iascii) {
473:     MPI_Comm          comm;
474:     PetscMPIInt       rank, size;
475:     PetscViewerFormat fmt;
476:     PetscBool         isperm;

478:     PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
479:     PetscCallMPI(MPI_Comm_rank(comm, &rank));
480:     PetscCallMPI(MPI_Comm_size(comm, &size));

482:     PetscCall(PetscViewerGetFormat(viewer, &fmt));
483:     PetscCall(ISGetInfo(is, IS_PERMUTATION, IS_LOCAL, PETSC_FALSE, &isperm));
484:     if (isperm && fmt != PETSC_VIEWER_ASCII_MATLAB) PetscCall(PetscViewerASCIIPrintf(viewer, "Index set is permutation\n"));
485:     PetscCall(PetscViewerASCIIPushSynchronized(viewer));
486:     if (size > 1) {
487:       if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
488:         const char *name;

490:         PetscCall(PetscObjectGetName((PetscObject)is, &name));
491:         PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%s_%d = [...\n", name, rank));
492:         for (PetscInt i = 0; i < n; i++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%" PetscInt_FMT "\n", idx[i] + 1));
493:         PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "];\n"));
494:       } else {
495:         PetscInt st = 0;

497:         if (fmt == PETSC_VIEWER_ASCII_INDEX) st = is->map->rstart;
498:         PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Number of indices in set %" PetscInt_FMT "\n", rank, n));
499:         for (PetscInt i = 0; i < n; i++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] %" PetscInt_FMT " %" PetscInt_FMT "\n", rank, i + st, idx[i]));
500:       }
501:     } else {
502:       if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
503:         const char *name;

505:         PetscCall(PetscObjectGetName((PetscObject)is, &name));
506:         PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%s = [...\n", name));
507:         for (PetscInt i = 0; i < n; i++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%" PetscInt_FMT "\n", idx[i] + 1));
508:         PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "];\n"));
509:       } else {
510:         PetscInt st = 0;

512:         if (fmt == PETSC_VIEWER_ASCII_INDEX) st = is->map->rstart;
513:         PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Number of indices in set %" PetscInt_FMT "\n", n));
514:         for (PetscInt i = 0; i < n; i++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT "\n", i + st, idx[i]));
515:       }
516:     }
517:     PetscCall(PetscViewerFlush(viewer));
518:     PetscCall(PetscViewerASCIIPopSynchronized(viewer));
519:   } else if (isbinary) {
520:     PetscCall(ISView_Binary(is, viewer));
521:   } else if (ishdf5) {
522: #if defined(PETSC_HAVE_HDF5)
523:     if (compress) PetscCall(ISView_General_HDF5_Compressed(is, viewer));
524:     else PetscCall(ISView_General_HDF5(is, viewer));
525: #endif
526:   }
527:   PetscFunctionReturn(PETSC_SUCCESS);
528: }

530: static PetscErrorCode ISSort_General(IS is)
531: {
532:   IS_General *sub = (IS_General *)is->data;
533:   PetscInt    n;

535:   PetscFunctionBegin;
536:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
537:   PetscCall(PetscIntSortSemiOrdered(n, sub->idx));
538:   PetscFunctionReturn(PETSC_SUCCESS);
539: }

541: static PetscErrorCode ISSortRemoveDups_General(IS is)
542: {
543:   IS_General *sub = (IS_General *)is->data;
544:   PetscLayout map;
545:   PetscInt    n;
546:   PetscBool   sorted;

548:   PetscFunctionBegin;
549:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
550:   PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, &sorted));
551:   if (sorted) {
552:     PetscCall(PetscSortedRemoveDupsInt(&n, sub->idx));
553:   } else {
554:     PetscCall(PetscSortRemoveDupsInt(&n, sub->idx));
555:   }
556:   PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), n, PETSC_DECIDE, is->map->bs, &map));
557:   PetscCall(PetscLayoutDestroy(&is->map));
558:   is->map = map;
559:   PetscFunctionReturn(PETSC_SUCCESS);
560: }

562: static PetscErrorCode ISSorted_General(IS is, PetscBool *flg)
563: {
564:   PetscFunctionBegin;
565:   PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, flg));
566:   PetscFunctionReturn(PETSC_SUCCESS);
567: }

569: static PetscErrorCode ISToGeneral_General(IS is)
570: {
571:   PetscFunctionBegin;
572:   PetscFunctionReturn(PETSC_SUCCESS);
573: }

575: // clang-format off
576: static const struct _ISOps myops = {
577:   PetscDesignatedInitializer(getindices, ISGetIndices_General),
578:   PetscDesignatedInitializer(restoreindices, ISRestoreIndices_General),
579:   PetscDesignatedInitializer(invertpermutation, ISInvertPermutation_General),
580:   PetscDesignatedInitializer(sort, ISSort_General),
581:   PetscDesignatedInitializer(sortremovedups, ISSortRemoveDups_General),
582:   PetscDesignatedInitializer(sorted, ISSorted_General),
583:   PetscDesignatedInitializer(duplicate, ISDuplicate_General),
584:   PetscDesignatedInitializer(destroy, ISDestroy_General),
585:   PetscDesignatedInitializer(view, ISView_General),
586:   PetscDesignatedInitializer(load, ISLoad_Default),
587:   PetscDesignatedInitializer(copy, ISCopy_General),
588:   PetscDesignatedInitializer(togeneral, ISToGeneral_General),
589:   PetscDesignatedInitializer(oncomm, ISOnComm_General),
590:   PetscDesignatedInitializer(setblocksize, ISSetBlockSize_General),
591:   PetscDesignatedInitializer(contiguous, ISContiguousLocal_General),
592:   PetscDesignatedInitializer(locate, ISLocate_General),
593:   /* no specializations of {sorted,unique,perm,interval}{local,global} because the default
594:    * checks in ISGetInfo_XXX in index.c are exactly what we would do for ISGeneral */
595:   PetscDesignatedInitializer(sortedlocal, NULL),
596:   PetscDesignatedInitializer(sortedglobal, NULL),
597:   PetscDesignatedInitializer(uniquelocal, NULL),
598:   PetscDesignatedInitializer(uniqueglobal, NULL),
599:   PetscDesignatedInitializer(permlocal, NULL),
600:   PetscDesignatedInitializer(permglobal, NULL),
601:   PetscDesignatedInitializer(intervallocal, NULL),
602:   PetscDesignatedInitializer(intervalglobal, NULL)
603: };
604: // clang-format on

606: PETSC_INTERN PetscErrorCode ISSetUp_General(IS is)
607: {
608:   IS_General     *sub = (IS_General *)is->data;
609:   const PetscInt *idx = sub->idx;
610:   PetscInt        n, i, min, max;

612:   PetscFunctionBegin;
613:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));

615:   if (n) {
616:     min = max = idx[0];
617:     for (i = 1; i < n; i++) {
618:       if (idx[i] < min) min = idx[i];
619:       if (idx[i] > max) max = idx[i];
620:     }
621:     is->min = min;
622:     is->max = max;
623:   } else {
624:     is->min = PETSC_INT_MAX;
625:     is->max = PETSC_INT_MIN;
626:   }
627:   PetscFunctionReturn(PETSC_SUCCESS);
628: }

630: /*@
631:   ISCreateGeneral - Creates a data structure for an index set containing a list of integers.

633:   Collective

635:   Input Parameters:
636: + comm - the MPI communicator
637: . n    - the length of the index set
638: . idx  - the list of integers
639: - mode - `PETSC_COPY_VALUES`, `PETSC_OWN_POINTER`, or `PETSC_USE_POINTER`; see `PetscCopyMode` for meaning of this flag.

641:   Output Parameter:
642: . is - the new index set

644:   Level: beginner

646:   Notes:
647:   When the communicator is not `MPI_COMM_SELF`, the operations on IS are NOT
648:   conceptually the same as `MPI_Group` operations. The `IS` are then
649:   distributed sets of indices and thus certain operations on them are
650:   collective.

652:   Use `ISGeneralSetIndices()` to provide indices to an already existing `IS` of `ISType` `ISGENERAL`

654: .seealso: [](sec_scatter), `IS`, `ISGENERAL`, `ISCreateStride()`, `ISCreateBlock()`, `ISAllGather()`, `PETSC_COPY_VALUES`, `PETSC_OWN_POINTER`,
655:           `PETSC_USE_POINTER`, `PetscCopyMode`, `ISGeneralSetIndicesFromMask()`
656: @*/
657: PetscErrorCode ISCreateGeneral(MPI_Comm comm, PetscInt n, const PetscInt idx[], PetscCopyMode mode, IS *is)
658: {
659:   PetscFunctionBegin;
660:   PetscCall(ISCreate(comm, is));
661:   PetscCall(ISSetType(*is, ISGENERAL));
662:   PetscCall(ISGeneralSetIndices(*is, n, idx, mode));
663:   PetscFunctionReturn(PETSC_SUCCESS);
664: }

666: /*@
667:   ISGeneralSetIndices - Sets the indices for an `ISGENERAL` index set

669:   Logically Collective

671:   Input Parameters:
672: + is   - the index set
673: . n    - the length of the index set
674: . idx  - the list of integers
675: - mode - see `PetscCopyMode` for meaning of this flag.

677:   Level: beginner

679:   Note:
680:   Use `ISCreateGeneral()` to create the `IS` and set its indices in a single function call

682: .seealso: [](sec_scatter), `IS`, `ISBLOCK`, `ISCreateGeneral()`, `ISGeneralSetIndicesFromMask()`, `ISBlockSetIndices()`, `ISGENERAL`, `PetscCopyMode`
683: @*/
684: PetscErrorCode ISGeneralSetIndices(IS is, PetscInt n, const PetscInt idx[], PetscCopyMode mode)
685: {
686:   PetscFunctionBegin;
688:   PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length < 0");
689:   if (n) PetscAssertPointer(idx, 3);
690:   PetscCall(ISClearInfoCache(is, PETSC_FALSE));
691:   PetscUseMethod(is, "ISGeneralSetIndices_C", (IS, PetscInt, const PetscInt[], PetscCopyMode), (is, n, idx, mode));
692:   PetscFunctionReturn(PETSC_SUCCESS);
693: }

695: static PetscErrorCode ISGeneralSetIndices_General(IS is, PetscInt n, const PetscInt idx[], PetscCopyMode mode)
696: {
697:   PetscLayout map;
698:   IS_General *sub = (IS_General *)is->data;

700:   PetscFunctionBegin;
701:   PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), n, PETSC_DECIDE, is->map->bs, &map));
702:   PetscCall(PetscLayoutDestroy(&is->map));
703:   is->map = map;

705:   if (sub->allocated) PetscCall(PetscFree(sub->idx));
706:   if (mode == PETSC_COPY_VALUES) {
707:     PetscCall(PetscMalloc1(n, &sub->idx));
708:     PetscCall(PetscArraycpy(sub->idx, idx, n));
709:     sub->allocated = PETSC_TRUE;
710:   } else if (mode == PETSC_OWN_POINTER) {
711:     sub->idx       = (PetscInt *)idx;
712:     sub->allocated = PETSC_TRUE;
713:   } else {
714:     sub->idx       = (PetscInt *)idx;
715:     sub->allocated = PETSC_FALSE;
716:   }

718:   PetscCall(ISSetUp_General(is));
719:   PetscCall(ISViewFromOptions(is, NULL, "-is_view"));
720:   PetscFunctionReturn(PETSC_SUCCESS);
721: }

723: /*@
724:   ISGeneralSetIndicesFromMask - Sets the indices for an `ISGENERAL` index set using a boolean mask

726:   Collective

728:   Input Parameters:
729: + is     - the index set
730: . rstart - the range start index (inclusive)
731: . rend   - the range end index (exclusive)
732: - mask   - the boolean mask array of length rend-rstart, indices will be set for each `PETSC_TRUE` value in the array

734:   Level: beginner

736:   Note:
737:   The mask array may be freed by the user after this call.

739:   Example:
740: .vb
741:    PetscBool mask[] = {PETSC_FALSE, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, PETSC_TRUE};
742:    ISGeneralSetIndicesFromMask(is,10,15,mask);
743: .ve
744:   will feed the `IS` with indices
745: .vb
746:   {11, 14}
747: .ve
748:   locally.

750: .seealso: [](sec_scatter), `IS`, `ISCreateGeneral()`, `ISGeneralSetIndices()`, `ISGENERAL`
751: @*/
752: PetscErrorCode ISGeneralSetIndicesFromMask(IS is, PetscInt rstart, PetscInt rend, const PetscBool mask[])
753: {
754:   PetscFunctionBegin;
756:   PetscCheck(rstart <= rend, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "rstart %" PetscInt_FMT " must be less than or equal to rend %" PetscInt_FMT, rstart, rend);
757:   if (rend - rstart) PetscAssertPointer(mask, 4);
758:   PetscCall(ISClearInfoCache(is, PETSC_FALSE));
759:   PetscUseMethod(is, "ISGeneralSetIndicesFromMask_C", (IS, PetscInt, PetscInt, const PetscBool[]), (is, rstart, rend, mask));
760:   PetscFunctionReturn(PETSC_SUCCESS);
761: }

763: static PetscErrorCode ISGeneralSetIndicesFromMask_General(IS is, PetscInt rstart, PetscInt rend, const PetscBool mask[])
764: {
765:   PetscInt  i, nidx;
766:   PetscInt *idx;

768:   PetscFunctionBegin;
769:   for (i = 0, nidx = 0; i < rend - rstart; i++)
770:     if (mask[i]) nidx++;
771:   PetscCall(PetscMalloc1(nidx, &idx));
772:   for (i = 0, nidx = 0; i < rend - rstart; i++) {
773:     if (mask[i]) {
774:       idx[nidx] = i + rstart;
775:       nidx++;
776:     }
777:   }
778:   PetscCall(ISGeneralSetIndices_General(is, nidx, idx, PETSC_OWN_POINTER));
779:   PetscFunctionReturn(PETSC_SUCCESS);
780: }

782: static PetscErrorCode ISGeneralFilter_General(IS is, PetscInt start, PetscInt end)
783: {
784:   IS_General *sub = (IS_General *)is->data;
785:   PetscInt   *idx = sub->idx, *idxnew;
786:   PetscInt    i, n = is->map->n, nnew = 0, o;

788:   PetscFunctionBegin;
789:   for (i = 0; i < n; ++i)
790:     if (idx[i] >= start && idx[i] < end) nnew++;
791:   PetscCall(PetscMalloc1(nnew, &idxnew));
792:   for (o = 0, i = 0; i < n; i++) {
793:     if (idx[i] >= start && idx[i] < end) idxnew[o++] = idx[i];
794:   }
795:   PetscCall(ISGeneralSetIndices_General(is, nnew, idxnew, PETSC_OWN_POINTER));
796:   PetscFunctionReturn(PETSC_SUCCESS);
797: }

799: /*@
800:   ISGeneralFilter - Remove all indices outside of [start, end) from an `ISGENERAL`

802:   Collective

804:   Input Parameters:
805: + is    - the index set
806: . start - the lowest index kept
807: - end   - one more than the highest index kept, `start` $\le$ `end`

809:   Level: beginner

811: .seealso: [](sec_scatter), `IS`, `ISGENERAL`, `ISCreateGeneral()`, `ISGeneralSetIndices()`
812: @*/
813: PetscErrorCode ISGeneralFilter(IS is, PetscInt start, PetscInt end)
814: {
815:   PetscFunctionBegin;
817:   PetscCheck(start <= end, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "start %" PetscInt_FMT " must be less than or equal to end %" PetscInt_FMT, start, end);
818:   PetscCall(ISClearInfoCache(is, PETSC_FALSE));
819:   PetscUseMethod(is, "ISGeneralFilter_C", (IS, PetscInt, PetscInt), (is, start, end));
820:   PetscFunctionReturn(PETSC_SUCCESS);
821: }

823: PETSC_INTERN PetscErrorCode ISCreate_General(IS is)
824: {
825:   IS_General *sub;

827:   PetscFunctionBegin;
828:   PetscCall(PetscNew(&sub));
829:   is->data   = (void *)sub;
830:   is->ops[0] = myops;
831:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISGeneralSetIndices_C", ISGeneralSetIndices_General));
832:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISGeneralSetIndicesFromMask_C", ISGeneralSetIndicesFromMask_General));
833:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISGeneralFilter_C", ISGeneralFilter_General));
834:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISShift_C", ISShift_General));
835:   PetscFunctionReturn(PETSC_SUCCESS);
836: }