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, ×tepping));
206: if (timestepping) PetscCall(PetscViewerHDF5GetTimestep(viewer, ×tep));
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, ×tepping));
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: }