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: #if defined(PETSC_HAVE_HDF5)
184: static PetscErrorCode ISView_General_HDF5(IS is, PetscViewer viewer)
185: {
186: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5 *)viewer->data;
187: hid_t filespace; /* file dataspace identifier */
188: hid_t chunkspace; /* chunk dataset property identifier */
189: hid_t dset_id; /* dataset identifier */
190: hid_t memspace; /* memory dataspace identifier */
191: hid_t inttype; /* int type (H5T_NATIVE_INT or H5T_NATIVE_LLONG) */
192: hid_t file_id, group;
193: hsize_t dim, maxDims[3], dims[3], chunkDims[3], count[3], offset[3];
194: PetscBool timestepping;
195: PetscInt bs, N, n, timestep = PETSC_MIN_INT, low;
196: hsize_t chunksize;
197: const PetscInt *ind;
198: const char *isname;
200: PetscFunctionBegin;
201: PetscCall(ISGetBlockSize(is, &bs));
202: bs = PetscMax(bs, 1); /* If N = 0, bs = 0 as well */
203: PetscCall(PetscViewerHDF5OpenGroup(viewer, NULL, &file_id, &group));
204: PetscCall(PetscViewerHDF5IsTimestepping(viewer, ×tepping));
205: if (timestepping) PetscCall(PetscViewerHDF5GetTimestep(viewer, ×tep));
207: /* Create the dataspace for the dataset.
208: *
209: * dims - holds the current dimensions of the dataset
210: *
211: * maxDims - holds the maximum dimensions of the dataset (unlimited
212: * for the number of time steps with the current dimensions for the
213: * other dimensions; so only additional time steps can be added).
214: *
215: * chunkDims - holds the size of a single time step (required to
216: * permit extending dataset).
217: */
218: dim = 0;
219: chunksize = 1;
220: if (timestep >= 0) {
221: dims[dim] = timestep + 1;
222: maxDims[dim] = H5S_UNLIMITED;
223: chunkDims[dim] = 1;
224: ++dim;
225: }
226: PetscCall(ISGetSize(is, &N));
227: PetscCall(ISGetLocalSize(is, &n));
228: PetscCall(PetscHDF5IntCast(N / bs, dims + dim));
230: maxDims[dim] = dims[dim];
231: chunkDims[dim] = PetscMax(1, dims[dim]);
232: chunksize *= chunkDims[dim];
233: ++dim;
234: if (bs >= 1) {
235: dims[dim] = bs;
236: maxDims[dim] = dims[dim];
237: chunkDims[dim] = dims[dim];
238: chunksize *= chunkDims[dim];
239: ++dim;
240: }
241: /* hdf5 chunks must be less than 4GB */
242: if (chunksize > PETSC_HDF5_MAX_CHUNKSIZE / 64) {
243: if (bs >= 1) {
244: 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));
245: 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));
246: } else {
247: chunkDims[dim - 1] = PETSC_HDF5_MAX_CHUNKSIZE / 64;
248: }
249: }
250: PetscCallHDF5Return(filespace, H5Screate_simple, (dim, dims, maxDims));
252: #if defined(PETSC_USE_64BIT_INDICES)
253: inttype = H5T_NATIVE_LLONG;
254: #else
255: inttype = H5T_NATIVE_INT;
256: #endif
258: /* Create the dataset with default properties and close filespace */
259: PetscCall(PetscObjectGetName((PetscObject)is, &isname));
260: if (!H5Lexists(group, isname, H5P_DEFAULT)) {
261: /* Create chunk */
262: PetscCallHDF5Return(chunkspace, H5Pcreate, (H5P_DATASET_CREATE));
263: PetscCallHDF5(H5Pset_chunk, (chunkspace, dim, chunkDims));
265: PetscCallHDF5Return(dset_id, H5Dcreate2, (group, isname, inttype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT));
266: PetscCallHDF5(H5Pclose, (chunkspace));
267: } else {
268: PetscCallHDF5Return(dset_id, H5Dopen2, (group, isname, H5P_DEFAULT));
269: PetscCallHDF5(H5Dset_extent, (dset_id, dims));
270: }
271: PetscCallHDF5(H5Sclose, (filespace));
273: /* Each process defines a dataset and writes it to the hyperslab in the file */
274: dim = 0;
275: if (timestep >= 0) {
276: count[dim] = 1;
277: ++dim;
278: }
279: PetscCall(PetscHDF5IntCast(n / bs, count + dim));
280: ++dim;
281: if (bs >= 1) {
282: count[dim] = bs;
283: ++dim;
284: }
285: if (n > 0 || H5_VERSION_GE(1, 10, 0)) {
286: PetscCallHDF5Return(memspace, H5Screate_simple, (dim, count, NULL));
287: } else {
288: /* Can't create dataspace with zero for any dimension, so create null dataspace. */
289: PetscCallHDF5Return(memspace, H5Screate, (H5S_NULL));
290: }
292: /* Select hyperslab in the file */
293: PetscCall(PetscLayoutGetRange(is->map, &low, NULL));
294: dim = 0;
295: if (timestep >= 0) {
296: offset[dim] = timestep;
297: ++dim;
298: }
299: PetscCall(PetscHDF5IntCast(low / bs, offset + dim));
300: ++dim;
301: if (bs >= 1) {
302: offset[dim] = 0;
303: ++dim;
304: }
305: if (n > 0 || H5_VERSION_GE(1, 10, 0)) {
306: PetscCallHDF5Return(filespace, H5Dget_space, (dset_id));
307: PetscCallHDF5(H5Sselect_hyperslab, (filespace, H5S_SELECT_SET, offset, NULL, count, NULL));
308: } else {
309: /* Create null filespace to match null memspace. */
310: PetscCallHDF5Return(filespace, H5Screate, (H5S_NULL));
311: }
313: PetscCall(ISGetIndices(is, &ind));
314: PetscCallHDF5(H5Dwrite, (dset_id, inttype, memspace, filespace, hdf5->dxpl_id, ind));
315: PetscCallHDF5(H5Fflush, (file_id, H5F_SCOPE_GLOBAL));
316: PetscCall(ISRestoreIndices(is, &ind));
318: /* Close/release resources */
319: PetscCallHDF5(H5Gclose, (group));
320: PetscCallHDF5(H5Sclose, (filespace));
321: PetscCallHDF5(H5Sclose, (memspace));
322: PetscCallHDF5(H5Dclose, (dset_id));
324: if (timestepping) PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)is, "timestepping", PETSC_BOOL, ×tepping));
325: PetscCall(PetscInfo(is, "Wrote IS object with name %s\n", isname));
326: PetscFunctionReturn(PETSC_SUCCESS);
327: }
328: #endif
330: static PetscErrorCode ISView_General(IS is, PetscViewer viewer)
331: {
332: IS_General *sub = (IS_General *)is->data;
333: PetscInt i, n, *idx = sub->idx;
334: PetscBool iascii, isbinary, ishdf5;
336: PetscFunctionBegin;
337: PetscCall(PetscLayoutGetLocalSize(is->map, &n));
338: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
339: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
340: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
341: if (iascii) {
342: MPI_Comm comm;
343: PetscMPIInt rank, size;
344: PetscViewerFormat fmt;
345: PetscBool isperm;
347: PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
348: PetscCallMPI(MPI_Comm_rank(comm, &rank));
349: PetscCallMPI(MPI_Comm_size(comm, &size));
351: PetscCall(PetscViewerGetFormat(viewer, &fmt));
352: PetscCall(ISGetInfo(is, IS_PERMUTATION, IS_LOCAL, PETSC_FALSE, &isperm));
353: if (isperm && fmt != PETSC_VIEWER_ASCII_MATLAB) PetscCall(PetscViewerASCIIPrintf(viewer, "Index set is permutation\n"));
354: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
355: if (size > 1) {
356: if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
357: const char *name;
359: PetscCall(PetscObjectGetName((PetscObject)is, &name));
360: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%s_%d = [...\n", name, rank));
361: for (i = 0; i < n; i++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%" PetscInt_FMT "\n", idx[i] + 1));
362: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "];\n"));
363: } else {
364: PetscInt st = 0;
366: if (fmt == PETSC_VIEWER_ASCII_INDEX) st = is->map->rstart;
367: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Number of indices in set %" PetscInt_FMT "\n", rank, n));
368: for (i = 0; i < n; i++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] %" PetscInt_FMT " %" PetscInt_FMT "\n", rank, i + st, idx[i]));
369: }
370: } else {
371: if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
372: const char *name;
374: PetscCall(PetscObjectGetName((PetscObject)is, &name));
375: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%s = [...\n", name));
376: for (i = 0; i < n; i++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%" PetscInt_FMT "\n", idx[i] + 1));
377: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "];\n"));
378: } else {
379: PetscInt st = 0;
381: if (fmt == PETSC_VIEWER_ASCII_INDEX) st = is->map->rstart;
382: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Number of indices in set %" PetscInt_FMT "\n", n));
383: for (i = 0; i < n; i++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT "\n", i + st, idx[i]));
384: }
385: }
386: PetscCall(PetscViewerFlush(viewer));
387: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
388: } else if (isbinary) {
389: PetscCall(ISView_Binary(is, viewer));
390: } else if (ishdf5) {
391: #if defined(PETSC_HAVE_HDF5)
392: PetscCall(ISView_General_HDF5(is, viewer));
393: #endif
394: }
395: PetscFunctionReturn(PETSC_SUCCESS);
396: }
398: static PetscErrorCode ISSort_General(IS is)
399: {
400: IS_General *sub = (IS_General *)is->data;
401: PetscInt n;
403: PetscFunctionBegin;
404: PetscCall(PetscLayoutGetLocalSize(is->map, &n));
405: PetscCall(PetscIntSortSemiOrdered(n, sub->idx));
406: PetscFunctionReturn(PETSC_SUCCESS);
407: }
409: static PetscErrorCode ISSortRemoveDups_General(IS is)
410: {
411: IS_General *sub = (IS_General *)is->data;
412: PetscLayout map;
413: PetscInt n;
414: PetscBool sorted;
416: PetscFunctionBegin;
417: PetscCall(PetscLayoutGetLocalSize(is->map, &n));
418: PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, &sorted));
419: if (sorted) {
420: PetscCall(PetscSortedRemoveDupsInt(&n, sub->idx));
421: } else {
422: PetscCall(PetscSortRemoveDupsInt(&n, sub->idx));
423: }
424: PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), n, PETSC_DECIDE, is->map->bs, &map));
425: PetscCall(PetscLayoutDestroy(&is->map));
426: is->map = map;
427: PetscFunctionReturn(PETSC_SUCCESS);
428: }
430: static PetscErrorCode ISSorted_General(IS is, PetscBool *flg)
431: {
432: PetscFunctionBegin;
433: PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, flg));
434: PetscFunctionReturn(PETSC_SUCCESS);
435: }
437: static PetscErrorCode ISToGeneral_General(IS is)
438: {
439: PetscFunctionBegin;
440: PetscFunctionReturn(PETSC_SUCCESS);
441: }
443: // clang-format off
444: static const struct _ISOps myops = {
445: PetscDesignatedInitializer(getindices, ISGetIndices_General),
446: PetscDesignatedInitializer(restoreindices, ISRestoreIndices_General),
447: PetscDesignatedInitializer(invertpermutation, ISInvertPermutation_General),
448: PetscDesignatedInitializer(sort, ISSort_General),
449: PetscDesignatedInitializer(sortremovedups, ISSortRemoveDups_General),
450: PetscDesignatedInitializer(sorted, ISSorted_General),
451: PetscDesignatedInitializer(duplicate, ISDuplicate_General),
452: PetscDesignatedInitializer(destroy, ISDestroy_General),
453: PetscDesignatedInitializer(view, ISView_General),
454: PetscDesignatedInitializer(load, ISLoad_Default),
455: PetscDesignatedInitializer(copy, ISCopy_General),
456: PetscDesignatedInitializer(togeneral, ISToGeneral_General),
457: PetscDesignatedInitializer(oncomm, ISOnComm_General),
458: PetscDesignatedInitializer(setblocksize, ISSetBlockSize_General),
459: PetscDesignatedInitializer(contiguous, ISContiguousLocal_General),
460: PetscDesignatedInitializer(locate, ISLocate_General),
461: /* no specializations of {sorted,unique,perm,interval}{local,global} because the default
462: * checks in ISGetInfo_XXX in index.c are exactly what we would do for ISGeneral */
463: PetscDesignatedInitializer(sortedlocal, NULL),
464: PetscDesignatedInitializer(sortedglobal, NULL),
465: PetscDesignatedInitializer(uniquelocal, NULL),
466: PetscDesignatedInitializer(uniqueglobal, NULL),
467: PetscDesignatedInitializer(permlocal, NULL),
468: PetscDesignatedInitializer(permglobal, NULL),
469: PetscDesignatedInitializer(intervallocal, NULL),
470: PetscDesignatedInitializer(intervalglobal, NULL)
471: };
472: // clang-format on
474: PETSC_INTERN PetscErrorCode ISSetUp_General(IS is)
475: {
476: IS_General *sub = (IS_General *)is->data;
477: const PetscInt *idx = sub->idx;
478: PetscInt n, i, min, max;
480: PetscFunctionBegin;
481: PetscCall(PetscLayoutGetLocalSize(is->map, &n));
483: if (n) {
484: min = max = idx[0];
485: for (i = 1; i < n; i++) {
486: if (idx[i] < min) min = idx[i];
487: if (idx[i] > max) max = idx[i];
488: }
489: is->min = min;
490: is->max = max;
491: } else {
492: is->min = PETSC_MAX_INT;
493: is->max = PETSC_MIN_INT;
494: }
495: PetscFunctionReturn(PETSC_SUCCESS);
496: }
498: /*@
499: ISCreateGeneral - Creates a data structure for an index set containing a list of integers.
501: Collective
503: Input Parameters:
504: + comm - the MPI communicator
505: . n - the length of the index set
506: . idx - the list of integers
507: - mode - `PETSC_COPY_VALUES`, `PETSC_OWN_POINTER`, or `PETSC_USE_POINTER`; see `PetscCopyMode` for meaning of this flag.
509: Output Parameter:
510: . is - the new index set
512: Level: beginner
514: Notes:
515: When the communicator is not `MPI_COMM_SELF`, the operations on IS are NOT
516: conceptually the same as `MPI_Group` operations. The `IS` are then
517: distributed sets of indices and thus certain operations on them are
518: collective.
520: Use `ISGeneralSetIndices()` to provide indices to an already existing `IS` of `ISType` `ISGENERAL`
522: .seealso: [](sec_scatter), `IS`, `ISGENERAL`, `ISCreateStride()`, `ISCreateBlock()`, `ISAllGather()`, `PETSC_COPY_VALUES`, `PETSC_OWN_POINTER`,
523: `PETSC_USE_POINTER`, `PetscCopyMode`, `ISGeneralSetIndicesFromMask()`
524: @*/
525: PetscErrorCode ISCreateGeneral(MPI_Comm comm, PetscInt n, const PetscInt idx[], PetscCopyMode mode, IS *is)
526: {
527: PetscFunctionBegin;
528: PetscCall(ISCreate(comm, is));
529: PetscCall(ISSetType(*is, ISGENERAL));
530: PetscCall(ISGeneralSetIndices(*is, n, idx, mode));
531: PetscFunctionReturn(PETSC_SUCCESS);
532: }
534: /*@
535: ISGeneralSetIndices - Sets the indices for an `ISGENERAL` index set
537: Logically Collective
539: Input Parameters:
540: + is - the index set
541: . n - the length of the index set
542: . idx - the list of integers
543: - mode - see `PetscCopyMode` for meaning of this flag.
545: Level: beginner
547: Note:
548: Use `ISCreateGeneral()` to create the `IS` and set its indices in a single function call
550: .seealso: [](sec_scatter), `IS`, `ISBLOCK`, `ISCreateGeneral()`, `ISGeneralSetIndicesFromMask()`, `ISBlockSetIndices()`, `ISGENERAL`, `PetscCopyMode`
551: @*/
552: PetscErrorCode ISGeneralSetIndices(IS is, PetscInt n, const PetscInt idx[], PetscCopyMode mode)
553: {
554: PetscFunctionBegin;
556: if (n) PetscAssertPointer(idx, 3);
557: PetscCall(ISClearInfoCache(is, PETSC_FALSE));
558: PetscUseMethod(is, "ISGeneralSetIndices_C", (IS, PetscInt, const PetscInt[], PetscCopyMode), (is, n, idx, mode));
559: PetscFunctionReturn(PETSC_SUCCESS);
560: }
562: static PetscErrorCode ISGeneralSetIndices_General(IS is, PetscInt n, const PetscInt idx[], PetscCopyMode mode)
563: {
564: PetscLayout map;
565: IS_General *sub = (IS_General *)is->data;
567: PetscFunctionBegin;
568: PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length < 0");
569: if (n) PetscAssertPointer(idx, 3);
571: PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), n, PETSC_DECIDE, is->map->bs, &map));
572: PetscCall(PetscLayoutDestroy(&is->map));
573: is->map = map;
575: if (sub->allocated) PetscCall(PetscFree(sub->idx));
576: if (mode == PETSC_COPY_VALUES) {
577: PetscCall(PetscMalloc1(n, &sub->idx));
578: PetscCall(PetscArraycpy(sub->idx, idx, n));
579: sub->allocated = PETSC_TRUE;
580: } else if (mode == PETSC_OWN_POINTER) {
581: sub->idx = (PetscInt *)idx;
582: sub->allocated = PETSC_TRUE;
583: } else {
584: sub->idx = (PetscInt *)idx;
585: sub->allocated = PETSC_FALSE;
586: }
588: PetscCall(ISSetUp_General(is));
589: PetscCall(ISViewFromOptions(is, NULL, "-is_view"));
590: PetscFunctionReturn(PETSC_SUCCESS);
591: }
593: /*@
594: ISGeneralSetIndicesFromMask - Sets the indices for an `ISGENERAL` index set using a boolean mask
596: Collective
598: Input Parameters:
599: + is - the index set
600: . rstart - the range start index (inclusive)
601: . rend - the range end index (exclusive)
602: - mask - the boolean mask array of length rend-rstart, indices will be set for each `PETSC_TRUE` value in the array
604: Level: beginner
606: Note:
607: The mask array may be freed by the user after this call.
609: Example:
610: .vb
611: PetscBool mask[] = {PETSC_FALSE, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, PETSC_TRUE};
612: ISGeneralSetIndicesFromMask(is,10,15,mask);
613: .ve
614: will feed the `IS` with indices
615: .vb
616: {11, 14}
617: .ve
618: locally.
620: .seealso: [](sec_scatter), `IS`, `ISCreateGeneral()`, `ISGeneralSetIndices()`, `ISGENERAL`
621: @*/
622: PetscErrorCode ISGeneralSetIndicesFromMask(IS is, PetscInt rstart, PetscInt rend, const PetscBool mask[])
623: {
624: PetscFunctionBegin;
626: if (rend - rstart) PetscAssertPointer(mask, 4);
627: PetscCall(ISClearInfoCache(is, PETSC_FALSE));
628: PetscUseMethod(is, "ISGeneralSetIndicesFromMask_C", (IS, PetscInt, PetscInt, const PetscBool[]), (is, rstart, rend, mask));
629: PetscFunctionReturn(PETSC_SUCCESS);
630: }
632: static PetscErrorCode ISGeneralSetIndicesFromMask_General(IS is, PetscInt rstart, PetscInt rend, const PetscBool mask[])
633: {
634: PetscInt i, nidx;
635: PetscInt *idx;
637: PetscFunctionBegin;
638: for (i = 0, nidx = 0; i < rend - rstart; i++)
639: if (mask[i]) nidx++;
640: PetscCall(PetscMalloc1(nidx, &idx));
641: for (i = 0, nidx = 0; i < rend - rstart; i++) {
642: if (mask[i]) {
643: idx[nidx] = i + rstart;
644: nidx++;
645: }
646: }
647: PetscCall(ISGeneralSetIndices_General(is, nidx, idx, PETSC_OWN_POINTER));
648: PetscFunctionReturn(PETSC_SUCCESS);
649: }
651: static PetscErrorCode ISGeneralFilter_General(IS is, PetscInt start, PetscInt end)
652: {
653: IS_General *sub = (IS_General *)is->data;
654: PetscInt *idx = sub->idx, *idxnew;
655: PetscInt i, n = is->map->n, nnew = 0, o;
657: PetscFunctionBegin;
658: for (i = 0; i < n; ++i)
659: if (idx[i] >= start && idx[i] < end) nnew++;
660: PetscCall(PetscMalloc1(nnew, &idxnew));
661: for (o = 0, i = 0; i < n; i++) {
662: if (idx[i] >= start && idx[i] < end) idxnew[o++] = idx[i];
663: }
664: PetscCall(ISGeneralSetIndices_General(is, nnew, idxnew, PETSC_OWN_POINTER));
665: PetscFunctionReturn(PETSC_SUCCESS);
666: }
668: /*@
669: ISGeneralFilter - Remove all indices outside of [start, end) from an `ISGENERAL`
671: Collective
673: Input Parameters:
674: + is - the index set
675: . start - the lowest index kept
676: - end - one more than the highest index kept
678: Level: beginner
680: .seealso: [](sec_scatter), `IS`, `ISGENERAL`, `ISCreateGeneral()`, `ISGeneralSetIndices()`
681: @*/
682: PetscErrorCode ISGeneralFilter(IS is, PetscInt start, PetscInt end)
683: {
684: PetscFunctionBegin;
686: PetscCall(ISClearInfoCache(is, PETSC_FALSE));
687: PetscUseMethod(is, "ISGeneralFilter_C", (IS, PetscInt, PetscInt), (is, start, end));
688: PetscFunctionReturn(PETSC_SUCCESS);
689: }
691: PETSC_INTERN PetscErrorCode ISCreate_General(IS is)
692: {
693: IS_General *sub;
695: PetscFunctionBegin;
696: PetscCall(PetscNew(&sub));
697: is->data = (void *)sub;
698: is->ops[0] = myops;
699: PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISGeneralSetIndices_C", ISGeneralSetIndices_General));
700: PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISGeneralSetIndicesFromMask_C", ISGeneralSetIndicesFromMask_General));
701: PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISGeneralFilter_C", ISGeneralFilter_General));
702: PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISShift_C", ISShift_General));
703: PetscFunctionReturn(PETSC_SUCCESS);
704: }