Actual source code: general.c


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

184: #if defined(PETSC_HAVE_HDF5)
185: static PetscErrorCode ISView_General_HDF5(IS is, PetscViewer viewer)
186: {
187:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5 *)viewer->data;
188:   hid_t             filespace;  /* file dataspace identifier */
189:   hid_t             chunkspace; /* chunk dataset property identifier */
190:   hid_t             dset_id;    /* dataset identifier */
191:   hid_t             memspace;   /* memory dataspace identifier */
192:   hid_t             inttype;    /* int type (H5T_NATIVE_INT or H5T_NATIVE_LLONG) */
193:   hid_t             file_id, group;
194:   hsize_t           dim, maxDims[3], dims[3], chunkDims[3], count[3], offset[3];
195:   PetscBool         timestepping;
196:   PetscInt          bs, N, n, timestep = PETSC_MIN_INT, low;
197:   hsize_t           chunksize;
198:   const PetscInt   *ind;
199:   const char       *isname;

201:   PetscFunctionBegin;
202:   PetscCall(ISGetBlockSize(is, &bs));
203:   bs = PetscMax(bs, 1); /* If N = 0, bs  = 0 as well */
204:   PetscCall(PetscViewerHDF5OpenGroup(viewer, NULL, &file_id, &group));
205:   PetscCall(PetscViewerHDF5IsTimestepping(viewer, &timestepping));
206:   if (timestepping) PetscCall(PetscViewerHDF5GetTimestep(viewer, &timestep));

208:   /* Create the dataspace for the dataset.
209:    *
210:    * dims - holds the current dimensions of the dataset
211:    *
212:    * maxDims - holds the maximum dimensions of the dataset (unlimited
213:    * for the number of time steps with the current dimensions for the
214:    * other dimensions; so only additional time steps can be added).
215:    *
216:    * chunkDims - holds the size of a single time step (required to
217:    * permit extending dataset).
218:    */
219:   dim       = 0;
220:   chunksize = 1;
221:   if (timestep >= 0) {
222:     dims[dim]      = timestep + 1;
223:     maxDims[dim]   = H5S_UNLIMITED;
224:     chunkDims[dim] = 1;
225:     ++dim;
226:   }
227:   PetscCall(ISGetSize(is, &N));
228:   PetscCall(ISGetLocalSize(is, &n));
229:   PetscCall(PetscHDF5IntCast(N / bs, dims + dim));

231:   maxDims[dim]   = dims[dim];
232:   chunkDims[dim] = PetscMax(1, dims[dim]);
233:   chunksize *= chunkDims[dim];
234:   ++dim;
235:   if (bs >= 1) {
236:     dims[dim]      = bs;
237:     maxDims[dim]   = dims[dim];
238:     chunkDims[dim] = dims[dim];
239:     chunksize *= chunkDims[dim];
240:     ++dim;
241:   }
242:   /* hdf5 chunks must be less than 4GB */
243:   if (chunksize > PETSC_HDF5_MAX_CHUNKSIZE / 64) {
244:     if (bs >= 1) {
245:       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));
246:       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));
247:     } else {
248:       chunkDims[dim - 1] = PETSC_HDF5_MAX_CHUNKSIZE / 64;
249:     }
250:   }
251:   PetscCallHDF5Return(filespace, H5Screate_simple, (dim, dims, maxDims));

253:   #if defined(PETSC_USE_64BIT_INDICES)
254:   inttype = H5T_NATIVE_LLONG;
255:   #else
256:   inttype = H5T_NATIVE_INT;
257:   #endif

259:   /* Create the dataset with default properties and close filespace */
260:   PetscCall(PetscObjectGetName((PetscObject)is, &isname));
261:   if (!H5Lexists(group, isname, H5P_DEFAULT)) {
262:     /* Create chunk */
263:     PetscCallHDF5Return(chunkspace, H5Pcreate, (H5P_DATASET_CREATE));
264:     PetscCallHDF5(H5Pset_chunk, (chunkspace, dim, chunkDims));

266:     PetscCallHDF5Return(dset_id, H5Dcreate2, (group, isname, inttype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT));
267:     PetscCallHDF5(H5Pclose, (chunkspace));
268:   } else {
269:     PetscCallHDF5Return(dset_id, H5Dopen2, (group, isname, H5P_DEFAULT));
270:     PetscCallHDF5(H5Dset_extent, (dset_id, dims));
271:   }
272:   PetscCallHDF5(H5Sclose, (filespace));

274:   /* Each process defines a dataset and writes it to the hyperslab in the file */
275:   dim = 0;
276:   if (timestep >= 0) {
277:     count[dim] = 1;
278:     ++dim;
279:   }
280:   PetscCall(PetscHDF5IntCast(n / bs, count + dim));
281:   ++dim;
282:   if (bs >= 1) {
283:     count[dim] = bs;
284:     ++dim;
285:   }
286:   if (n > 0 || H5_VERSION_GE(1, 10, 0)) {
287:     PetscCallHDF5Return(memspace, H5Screate_simple, (dim, count, NULL));
288:   } else {
289:     /* Can't create dataspace with zero for any dimension, so create null dataspace. */
290:     PetscCallHDF5Return(memspace, H5Screate, (H5S_NULL));
291:   }

293:   /* Select hyperslab in the file */
294:   PetscCall(PetscLayoutGetRange(is->map, &low, NULL));
295:   dim = 0;
296:   if (timestep >= 0) {
297:     offset[dim] = timestep;
298:     ++dim;
299:   }
300:   PetscCall(PetscHDF5IntCast(low / bs, offset + dim));
301:   ++dim;
302:   if (bs >= 1) {
303:     offset[dim] = 0;
304:     ++dim;
305:   }
306:   if (n > 0 || H5_VERSION_GE(1, 10, 0)) {
307:     PetscCallHDF5Return(filespace, H5Dget_space, (dset_id));
308:     PetscCallHDF5(H5Sselect_hyperslab, (filespace, H5S_SELECT_SET, offset, NULL, count, NULL));
309:   } else {
310:     /* Create null filespace to match null memspace. */
311:     PetscCallHDF5Return(filespace, H5Screate, (H5S_NULL));
312:   }

314:   PetscCall(ISGetIndices(is, &ind));
315:   PetscCallHDF5(H5Dwrite, (dset_id, inttype, memspace, filespace, hdf5->dxpl_id, ind));
316:   PetscCallHDF5(H5Fflush, (file_id, H5F_SCOPE_GLOBAL));
317:   PetscCall(ISRestoreIndices(is, &ind));

319:   /* Close/release resources */
320:   PetscCallHDF5(H5Gclose, (group));
321:   PetscCallHDF5(H5Sclose, (filespace));
322:   PetscCallHDF5(H5Sclose, (memspace));
323:   PetscCallHDF5(H5Dclose, (dset_id));

325:   if (timestepping) PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)is, "timestepping", PETSC_BOOL, &timestepping));
326:   PetscCall(PetscInfo(is, "Wrote IS object with name %s\n", isname));
327:   PetscFunctionReturn(PETSC_SUCCESS);
328: }
329: #endif

331: static PetscErrorCode ISView_General(IS is, PetscViewer viewer)
332: {
333:   IS_General *sub = (IS_General *)is->data;
334:   PetscInt    i, n, *idx = sub->idx;
335:   PetscBool   iascii, isbinary, ishdf5;

337:   PetscFunctionBegin;
338:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
339:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
340:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
341:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
342:   if (iascii) {
343:     MPI_Comm          comm;
344:     PetscMPIInt       rank, size;
345:     PetscViewerFormat fmt;
346:     PetscBool         isperm;

348:     PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
349:     PetscCallMPI(MPI_Comm_rank(comm, &rank));
350:     PetscCallMPI(MPI_Comm_size(comm, &size));

352:     PetscCall(PetscViewerGetFormat(viewer, &fmt));
353:     PetscCall(ISGetInfo(is, IS_PERMUTATION, IS_LOCAL, PETSC_FALSE, &isperm));
354:     if (isperm && fmt != PETSC_VIEWER_ASCII_MATLAB) PetscCall(PetscViewerASCIIPrintf(viewer, "Index set is permutation\n"));
355:     PetscCall(PetscViewerASCIIPushSynchronized(viewer));
356:     if (size > 1) {
357:       if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
358:         const char *name;

360:         PetscCall(PetscObjectGetName((PetscObject)is, &name));
361:         PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%s_%d = [...\n", name, rank));
362:         for (i = 0; i < n; i++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%" PetscInt_FMT "\n", idx[i] + 1));
363:         PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "];\n"));
364:       } else {
365:         PetscInt st = 0;

367:         if (fmt == PETSC_VIEWER_ASCII_INDEX) st = is->map->rstart;
368:         PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Number of indices in set %" PetscInt_FMT "\n", rank, n));
369:         for (i = 0; i < n; i++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] %" PetscInt_FMT " %" PetscInt_FMT "\n", rank, i + st, idx[i]));
370:       }
371:     } else {
372:       if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
373:         const char *name;

375:         PetscCall(PetscObjectGetName((PetscObject)is, &name));
376:         PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%s = [...\n", name));
377:         for (i = 0; i < n; i++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%" PetscInt_FMT "\n", idx[i] + 1));
378:         PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "];\n"));
379:       } else {
380:         PetscInt st = 0;

382:         if (fmt == PETSC_VIEWER_ASCII_INDEX) st = is->map->rstart;
383:         PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Number of indices in set %" PetscInt_FMT "\n", n));
384:         for (i = 0; i < n; i++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT "\n", i + st, idx[i]));
385:       }
386:     }
387:     PetscCall(PetscViewerFlush(viewer));
388:     PetscCall(PetscViewerASCIIPopSynchronized(viewer));
389:   } else if (isbinary) {
390:     PetscCall(ISView_Binary(is, viewer));
391:   } else if (ishdf5) {
392: #if defined(PETSC_HAVE_HDF5)
393:     PetscCall(ISView_General_HDF5(is, viewer));
394: #endif
395:   }
396:   PetscFunctionReturn(PETSC_SUCCESS);
397: }

399: static PetscErrorCode ISSort_General(IS is)
400: {
401:   IS_General *sub = (IS_General *)is->data;
402:   PetscInt    n;

404:   PetscFunctionBegin;
405:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
406:   PetscCall(PetscIntSortSemiOrdered(n, sub->idx));
407:   PetscFunctionReturn(PETSC_SUCCESS);
408: }

410: static PetscErrorCode ISSortRemoveDups_General(IS is)
411: {
412:   IS_General *sub = (IS_General *)is->data;
413:   PetscLayout map;
414:   PetscInt    n;
415:   PetscBool   sorted;

417:   PetscFunctionBegin;
418:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
419:   PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, &sorted));
420:   if (sorted) {
421:     PetscCall(PetscSortedRemoveDupsInt(&n, sub->idx));
422:   } else {
423:     PetscCall(PetscSortRemoveDupsInt(&n, sub->idx));
424:   }
425:   PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), n, PETSC_DECIDE, is->map->bs, &map));
426:   PetscCall(PetscLayoutDestroy(&is->map));
427:   is->map = map;
428:   PetscFunctionReturn(PETSC_SUCCESS);
429: }

431: static PetscErrorCode ISSorted_General(IS is, PetscBool *flg)
432: {
433:   PetscFunctionBegin;
434:   PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, flg));
435:   PetscFunctionReturn(PETSC_SUCCESS);
436: }

438: PetscErrorCode ISToGeneral_General(IS is)
439: {
440:   PetscFunctionBegin;
441:   PetscFunctionReturn(PETSC_SUCCESS);
442: }

444: static struct _ISOps myops = {ISGetIndices_General, ISRestoreIndices_General, ISInvertPermutation_General, ISSort_General, ISSortRemoveDups_General, ISSorted_General, ISDuplicate_General, ISDestroy_General, ISView_General, ISLoad_Default, ISCopy_General, ISToGeneral_General, ISOnComm_General, ISSetBlockSize_General, ISContiguousLocal_General, ISLocate_General,
445:                               /* no specializations of {sorted,unique,perm,interval}{local,global}
446:                                 * because the default checks in ISGetInfo_XXX in index.c are exactly
447:                                 * what we would do for ISGeneral */
448:                               NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL};

450: PETSC_INTERN PetscErrorCode ISSetUp_General(IS);

452: PetscErrorCode ISSetUp_General(IS is)
453: {
454:   IS_General     *sub = (IS_General *)is->data;
455:   const PetscInt *idx = sub->idx;
456:   PetscInt        n, i, min, max;

458:   PetscFunctionBegin;
459:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));

461:   if (n) {
462:     min = max = idx[0];
463:     for (i = 1; i < n; i++) {
464:       if (idx[i] < min) min = idx[i];
465:       if (idx[i] > max) max = idx[i];
466:     }
467:     is->min = min;
468:     is->max = max;
469:   } else {
470:     is->min = PETSC_MAX_INT;
471:     is->max = PETSC_MIN_INT;
472:   }
473:   PetscFunctionReturn(PETSC_SUCCESS);
474: }

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

479:    Collective

481:    Input Parameters:
482: +  comm - the MPI communicator
483: .  n - the length of the index set
484: .  idx - the list of integers
485: -  mode - `PETSC_COPY_VALUES`, `PETSC_OWN_POINTER`, or `PETSC_USE_POINTER`; see `PetscCopyMode` for meaning of this flag.

487:    Output Parameter:
488: .  is - the new index set

490:    Level: beginner

492:    Notes:
493:    When the communicator is not `MPI_COMM_SELF`, the operations on IS are NOT
494:    conceptually the same as `MPI_Group` operations. The `IS` are then
495:    distributed sets of indices and thus certain operations on them are
496:    collective.

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

500: .seealso: [](sec_scatter), `IS`, `ISGENERAL`, `ISCreateStride()`, `ISCreateBlock()`, `ISAllGather()`, `PETSC_COPY_VALUES`, `PETSC_OWN_POINTER`,
501:           `PETSC_USE_POINTER`, `PetscCopyMode`, `ISGeneralSetIndicesFromMask()`
502: @*/
503: PetscErrorCode ISCreateGeneral(MPI_Comm comm, PetscInt n, const PetscInt idx[], PetscCopyMode mode, IS *is)
504: {
505:   PetscFunctionBegin;
506:   PetscCall(ISCreate(comm, is));
507:   PetscCall(ISSetType(*is, ISGENERAL));
508:   PetscCall(ISGeneralSetIndices(*is, n, idx, mode));
509:   PetscFunctionReturn(PETSC_SUCCESS);
510: }

512: /*@
513:    ISGeneralSetIndices - Sets the indices for an `ISGENERAL` index set

515:    Logically Collective

517:    Input Parameters:
518: +  is - the index set
519: .  n - the length of the index set
520: .  idx - the list of integers
521: -  mode - see `PetscCopyMode` for meaning of this flag.

523:    Level: beginner

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

528: .seealso: [](sec_scatter), `IS`, `ISBLOCK`, `ISCreateGeneral()`, `ISGeneralSetIndicesFromMask()`, `ISBlockSetIndices()`, `ISGENERAL`, `PetscCopyMode`
529: @*/
530: PetscErrorCode ISGeneralSetIndices(IS is, PetscInt n, const PetscInt idx[], PetscCopyMode mode)
531: {
532:   PetscFunctionBegin;
535:   PetscCall(ISClearInfoCache(is, PETSC_FALSE));
536:   PetscUseMethod(is, "ISGeneralSetIndices_C", (IS, PetscInt, const PetscInt[], PetscCopyMode), (is, n, idx, mode));
537:   PetscFunctionReturn(PETSC_SUCCESS);
538: }

540: PetscErrorCode ISGeneralSetIndices_General(IS is, PetscInt n, const PetscInt idx[], PetscCopyMode mode)
541: {
542:   PetscLayout map;
543:   IS_General *sub = (IS_General *)is->data;

545:   PetscFunctionBegin;
546:   PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length < 0");

549:   PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), n, PETSC_DECIDE, is->map->bs, &map));
550:   PetscCall(PetscLayoutDestroy(&is->map));
551:   is->map = map;

553:   if (sub->allocated) PetscCall(PetscFree(sub->idx));
554:   if (mode == PETSC_COPY_VALUES) {
555:     PetscCall(PetscMalloc1(n, &sub->idx));
556:     PetscCall(PetscArraycpy(sub->idx, idx, n));
557:     sub->allocated = PETSC_TRUE;
558:   } else if (mode == PETSC_OWN_POINTER) {
559:     sub->idx       = (PetscInt *)idx;
560:     sub->allocated = PETSC_TRUE;
561:   } else {
562:     sub->idx       = (PetscInt *)idx;
563:     sub->allocated = PETSC_FALSE;
564:   }

566:   PetscCall(ISSetUp_General(is));
567:   PetscCall(ISViewFromOptions(is, NULL, "-is_view"));
568:   PetscFunctionReturn(PETSC_SUCCESS);
569: }

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

574:    Collective

576:    Input Parameters:
577: +  is - the index set
578: .  rstart - the range start index (inclusive)
579: .  rend - the range end index (exclusive)
580: -  mask - the boolean mask array of length rend-rstart, indices will be set for each `PETSC_TRUE` value in the array

582:    Level: beginner

584:    Note:
585:    The mask array may be freed by the user after this call.

587:    Example:
588: .vb
589:    PetscBool mask[] = {PETSC_FALSE, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, PETSC_TRUE};
590:    ISGeneralSetIndicesFromMask(is,10,15,mask);
591: .ve
592:    will feed the `IS` with indices
593: .vb
594:   {11, 14}
595: .ve
596:   locally.

598: .seealso: [](sec_scatter), `IS`, `ISCreateGeneral()`, `ISGeneralSetIndices()`, `ISGENERAL`
599: @*/
600: PetscErrorCode ISGeneralSetIndicesFromMask(IS is, PetscInt rstart, PetscInt rend, const PetscBool mask[])
601: {
602:   PetscFunctionBegin;
605:   PetscCall(ISClearInfoCache(is, PETSC_FALSE));
606:   PetscUseMethod(is, "ISGeneralSetIndicesFromMask_C", (IS, PetscInt, PetscInt, const PetscBool[]), (is, rstart, rend, mask));
607:   PetscFunctionReturn(PETSC_SUCCESS);
608: }

610: PetscErrorCode ISGeneralSetIndicesFromMask_General(IS is, PetscInt rstart, PetscInt rend, const PetscBool mask[])
611: {
612:   PetscInt  i, nidx;
613:   PetscInt *idx;

615:   PetscFunctionBegin;
616:   for (i = 0, nidx = 0; i < rend - rstart; i++)
617:     if (mask[i]) nidx++;
618:   PetscCall(PetscMalloc1(nidx, &idx));
619:   for (i = 0, nidx = 0; i < rend - rstart; i++) {
620:     if (mask[i]) {
621:       idx[nidx] = i + rstart;
622:       nidx++;
623:     }
624:   }
625:   PetscCall(ISGeneralSetIndices_General(is, nidx, idx, PETSC_OWN_POINTER));
626:   PetscFunctionReturn(PETSC_SUCCESS);
627: }

629: static PetscErrorCode ISGeneralFilter_General(IS is, PetscInt start, PetscInt end)
630: {
631:   IS_General *sub = (IS_General *)is->data;
632:   PetscInt   *idx = sub->idx, *idxnew;
633:   PetscInt    i, n = is->map->n, nnew = 0, o;

635:   PetscFunctionBegin;
636:   for (i = 0; i < n; ++i)
637:     if (idx[i] >= start && idx[i] < end) nnew++;
638:   PetscCall(PetscMalloc1(nnew, &idxnew));
639:   for (o = 0, i = 0; i < n; i++) {
640:     if (idx[i] >= start && idx[i] < end) idxnew[o++] = idx[i];
641:   }
642:   PetscCall(ISGeneralSetIndices_General(is, nnew, idxnew, PETSC_OWN_POINTER));
643:   PetscFunctionReturn(PETSC_SUCCESS);
644: }

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

649:    Collective

651:    Input Parameters:
652: +  is - the index set
653: .  start - the lowest index kept
654: -  end - one more than the highest index kept

656:    Level: beginner

658: .seealso: [](sec_scatter), `IS`, `ISGENERAL`, `ISCreateGeneral()`, `ISGeneralSetIndices()`
659: @*/
660: PetscErrorCode ISGeneralFilter(IS is, PetscInt start, PetscInt end)
661: {
662:   PetscFunctionBegin;
664:   PetscCall(ISClearInfoCache(is, PETSC_FALSE));
665:   PetscUseMethod(is, "ISGeneralFilter_C", (IS, PetscInt, PetscInt), (is, start, end));
666:   PetscFunctionReturn(PETSC_SUCCESS);
667: }

669: PETSC_EXTERN PetscErrorCode ISCreate_General(IS is)
670: {
671:   IS_General *sub;

673:   PetscFunctionBegin;
674:   PetscCall(PetscNew(&sub));
675:   is->data = (void *)sub;
676:   PetscCall(PetscMemcpy(is->ops, &myops, sizeof(myops)));
677:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISGeneralSetIndices_C", ISGeneralSetIndices_General));
678:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISGeneralSetIndicesFromMask_C", ISGeneralSetIndicesFromMask_General));
679:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISGeneralFilter_C", ISGeneralFilter_General));
680:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISShift_C", ISShift_General));
681:   PetscFunctionReturn(PETSC_SUCCESS);
682: }