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 ISRestoreArray() 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;
225: const PetscInt *idx;
226: PetscInt n, off = 0;
228: PetscFunctionBegin;
229: *compress = PETSC_FALSE;
230: PetscCall(ISGetCompressOutput(is, &lcompress));
231: if (!lcompress) 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, ×tepping));
314: if (timestepping) PetscCall(PetscViewerHDF5GetTimestep(viewer, ×tep));
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, ×tepping));
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, ×tepping));
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: PetscBool vcompress;
525: PetscCall(PetscViewerHDF5GetCompress(viewer, &vcompress));
526: if (vcompress && compress) PetscCall(ISView_General_HDF5_Compressed(is, viewer));
527: else PetscCall(ISView_General_HDF5(is, viewer));
528: #endif
529: }
530: PetscFunctionReturn(PETSC_SUCCESS);
531: }
533: static PetscErrorCode ISSort_General(IS is)
534: {
535: IS_General *sub = (IS_General *)is->data;
536: PetscInt n;
538: PetscFunctionBegin;
539: PetscCall(PetscLayoutGetLocalSize(is->map, &n));
540: PetscCall(PetscIntSortSemiOrdered(n, sub->idx));
541: PetscFunctionReturn(PETSC_SUCCESS);
542: }
544: static PetscErrorCode ISSortRemoveDups_General(IS is)
545: {
546: IS_General *sub = (IS_General *)is->data;
547: PetscLayout map;
548: PetscInt n;
549: PetscBool sorted;
551: PetscFunctionBegin;
552: PetscCall(PetscLayoutGetLocalSize(is->map, &n));
553: PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, &sorted));
554: if (sorted) {
555: PetscCall(PetscSortedRemoveDupsInt(&n, sub->idx));
556: } else {
557: PetscCall(PetscSortRemoveDupsInt(&n, sub->idx));
558: }
559: PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), n, PETSC_DECIDE, is->map->bs, &map));
560: PetscCall(PetscLayoutDestroy(&is->map));
561: is->map = map;
562: PetscFunctionReturn(PETSC_SUCCESS);
563: }
565: static PetscErrorCode ISSorted_General(IS is, PetscBool *flg)
566: {
567: PetscFunctionBegin;
568: PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, flg));
569: PetscFunctionReturn(PETSC_SUCCESS);
570: }
572: static PetscErrorCode ISToGeneral_General(IS is)
573: {
574: PetscFunctionBegin;
575: PetscFunctionReturn(PETSC_SUCCESS);
576: }
578: // clang-format off
579: static const struct _ISOps myops = {
580: PetscDesignatedInitializer(getindices, ISGetIndices_General),
581: PetscDesignatedInitializer(restoreindices, ISRestoreIndices_General),
582: PetscDesignatedInitializer(invertpermutation, ISInvertPermutation_General),
583: PetscDesignatedInitializer(sort, ISSort_General),
584: PetscDesignatedInitializer(sortremovedups, ISSortRemoveDups_General),
585: PetscDesignatedInitializer(sorted, ISSorted_General),
586: PetscDesignatedInitializer(duplicate, ISDuplicate_General),
587: PetscDesignatedInitializer(destroy, ISDestroy_General),
588: PetscDesignatedInitializer(view, ISView_General),
589: PetscDesignatedInitializer(load, ISLoad_Default),
590: PetscDesignatedInitializer(copy, ISCopy_General),
591: PetscDesignatedInitializer(togeneral, ISToGeneral_General),
592: PetscDesignatedInitializer(oncomm, ISOnComm_General),
593: PetscDesignatedInitializer(setblocksize, ISSetBlockSize_General),
594: PetscDesignatedInitializer(contiguous, ISContiguousLocal_General),
595: PetscDesignatedInitializer(locate, ISLocate_General),
596: /* no specializations of {sorted,unique,perm,interval}{local,global} because the default
597: * checks in ISGetInfo_XXX in index.c are exactly what we would do for ISGeneral */
598: PetscDesignatedInitializer(sortedlocal, NULL),
599: PetscDesignatedInitializer(sortedglobal, NULL),
600: PetscDesignatedInitializer(uniquelocal, NULL),
601: PetscDesignatedInitializer(uniqueglobal, NULL),
602: PetscDesignatedInitializer(permlocal, NULL),
603: PetscDesignatedInitializer(permglobal, NULL),
604: PetscDesignatedInitializer(intervallocal, NULL),
605: PetscDesignatedInitializer(intervalglobal, NULL)
606: };
607: // clang-format on
609: PETSC_INTERN PetscErrorCode ISSetUp_General(IS is)
610: {
611: IS_General *sub = (IS_General *)is->data;
612: const PetscInt *idx = sub->idx;
613: PetscInt n, i, min, max;
615: PetscFunctionBegin;
616: PetscCall(PetscLayoutGetLocalSize(is->map, &n));
618: if (n) {
619: min = max = idx[0];
620: for (i = 1; i < n; i++) {
621: if (idx[i] < min) min = idx[i];
622: if (idx[i] > max) max = idx[i];
623: }
624: is->min = min;
625: is->max = max;
626: } else {
627: is->min = PETSC_INT_MAX;
628: is->max = PETSC_INT_MIN;
629: }
630: PetscFunctionReturn(PETSC_SUCCESS);
631: }
633: /*@
634: ISCreateGeneral - Creates a data structure for an index set containing a list of integers.
636: Collective
638: Input Parameters:
639: + comm - the MPI communicator
640: . n - the length of the index set
641: . idx - the list of integers
642: - mode - `PETSC_COPY_VALUES`, `PETSC_OWN_POINTER`, or `PETSC_USE_POINTER`; see `PetscCopyMode` for meaning of this flag.
644: Output Parameter:
645: . is - the new index set
647: Level: beginner
649: Notes:
650: When the communicator is not `MPI_COMM_SELF`, the operations on IS are NOT
651: conceptually the same as `MPI_Group` operations. The `IS` are then
652: distributed sets of indices and thus certain operations on them are
653: collective.
655: Use `ISGeneralSetIndices()` to provide indices to an already existing `IS` of `ISType` `ISGENERAL`
657: .seealso: [](sec_scatter), `IS`, `ISGENERAL`, `ISCreateStride()`, `ISCreateBlock()`, `ISAllGather()`, `PETSC_COPY_VALUES`, `PETSC_OWN_POINTER`,
658: `PETSC_USE_POINTER`, `PetscCopyMode`, `ISGeneralSetIndicesFromMask()`
659: @*/
660: PetscErrorCode ISCreateGeneral(MPI_Comm comm, PetscInt n, const PetscInt idx[], PetscCopyMode mode, IS *is)
661: {
662: PetscFunctionBegin;
663: PetscCall(ISCreate(comm, is));
664: PetscCall(ISSetType(*is, ISGENERAL));
665: PetscCall(ISGeneralSetIndices(*is, n, idx, mode));
666: PetscFunctionReturn(PETSC_SUCCESS);
667: }
669: /*@
670: ISGeneralSetIndices - Sets the indices for an `ISGENERAL` index set
672: Logically Collective
674: Input Parameters:
675: + is - the index set
676: . n - the length of the index set
677: . idx - the list of integers
678: - mode - see `PetscCopyMode` for meaning of this flag.
680: Level: beginner
682: Note:
683: Use `ISCreateGeneral()` to create the `IS` and set its indices in a single function call
685: .seealso: [](sec_scatter), `IS`, `ISBLOCK`, `ISCreateGeneral()`, `ISGeneralSetIndicesFromMask()`, `ISBlockSetIndices()`, `ISGENERAL`, `PetscCopyMode`
686: @*/
687: PetscErrorCode ISGeneralSetIndices(IS is, PetscInt n, const PetscInt idx[], PetscCopyMode mode)
688: {
689: PetscFunctionBegin;
691: PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length < 0");
692: if (n) PetscAssertPointer(idx, 3);
693: PetscCall(ISClearInfoCache(is, PETSC_FALSE));
694: PetscUseMethod(is, "ISGeneralSetIndices_C", (IS, PetscInt, const PetscInt[], PetscCopyMode), (is, n, idx, mode));
695: PetscFunctionReturn(PETSC_SUCCESS);
696: }
698: static PetscErrorCode ISGeneralSetIndices_General(IS is, PetscInt n, const PetscInt idx[], PetscCopyMode mode)
699: {
700: PetscLayout map;
701: IS_General *sub = (IS_General *)is->data;
703: PetscFunctionBegin;
704: PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), n, PETSC_DECIDE, is->map->bs, &map));
705: PetscCall(PetscLayoutDestroy(&is->map));
706: is->map = map;
708: if (sub->allocated) PetscCall(PetscFree(sub->idx));
709: if (mode == PETSC_COPY_VALUES) {
710: PetscCall(PetscMalloc1(n, &sub->idx));
711: PetscCall(PetscArraycpy(sub->idx, idx, n));
712: sub->allocated = PETSC_TRUE;
713: } else if (mode == PETSC_OWN_POINTER) {
714: sub->idx = (PetscInt *)idx;
715: sub->allocated = PETSC_TRUE;
716: } else {
717: sub->idx = (PetscInt *)idx;
718: sub->allocated = PETSC_FALSE;
719: }
721: PetscCall(ISSetUp_General(is));
722: PetscCall(ISViewFromOptions(is, NULL, "-is_view"));
723: PetscFunctionReturn(PETSC_SUCCESS);
724: }
726: /*@
727: ISGeneralSetIndicesFromMask - Sets the indices for an `ISGENERAL` index set using a boolean mask
729: Collective
731: Input Parameters:
732: + is - the index set
733: . rstart - the range start index (inclusive)
734: . rend - the range end index (exclusive)
735: - mask - the boolean mask array of length rend-rstart, indices will be set for each `PETSC_TRUE` value in the array
737: Level: beginner
739: Note:
740: The mask array may be freed by the user after this call.
742: Example:
743: .vb
744: PetscBool mask[] = {PETSC_FALSE, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, PETSC_TRUE};
745: ISGeneralSetIndicesFromMask(is,10,15,mask);
746: .ve
747: will feed the `IS` with indices
748: .vb
749: {11, 14}
750: .ve
751: locally.
753: .seealso: [](sec_scatter), `IS`, `ISCreateGeneral()`, `ISGeneralSetIndices()`, `ISGENERAL`
754: @*/
755: PetscErrorCode ISGeneralSetIndicesFromMask(IS is, PetscInt rstart, PetscInt rend, const PetscBool mask[])
756: {
757: PetscFunctionBegin;
759: 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);
760: if (rend - rstart) PetscAssertPointer(mask, 4);
761: PetscCall(ISClearInfoCache(is, PETSC_FALSE));
762: PetscUseMethod(is, "ISGeneralSetIndicesFromMask_C", (IS, PetscInt, PetscInt, const PetscBool[]), (is, rstart, rend, mask));
763: PetscFunctionReturn(PETSC_SUCCESS);
764: }
766: static PetscErrorCode ISGeneralSetIndicesFromMask_General(IS is, PetscInt rstart, PetscInt rend, const PetscBool mask[])
767: {
768: PetscInt i, nidx;
769: PetscInt *idx;
771: PetscFunctionBegin;
772: for (i = 0, nidx = 0; i < rend - rstart; i++)
773: if (mask[i]) nidx++;
774: PetscCall(PetscMalloc1(nidx, &idx));
775: for (i = 0, nidx = 0; i < rend - rstart; i++) {
776: if (mask[i]) {
777: idx[nidx] = i + rstart;
778: nidx++;
779: }
780: }
781: PetscCall(ISGeneralSetIndices_General(is, nidx, idx, PETSC_OWN_POINTER));
782: PetscFunctionReturn(PETSC_SUCCESS);
783: }
785: static PetscErrorCode ISGeneralFilter_General(IS is, PetscInt start, PetscInt end)
786: {
787: IS_General *sub = (IS_General *)is->data;
788: PetscInt *idx = sub->idx, *idxnew;
789: PetscInt i, n = is->map->n, nnew = 0, o;
791: PetscFunctionBegin;
792: for (i = 0; i < n; ++i)
793: if (idx[i] >= start && idx[i] < end) nnew++;
794: PetscCall(PetscMalloc1(nnew, &idxnew));
795: for (o = 0, i = 0; i < n; i++) {
796: if (idx[i] >= start && idx[i] < end) idxnew[o++] = idx[i];
797: }
798: PetscCall(ISGeneralSetIndices_General(is, nnew, idxnew, PETSC_OWN_POINTER));
799: PetscFunctionReturn(PETSC_SUCCESS);
800: }
802: /*@
803: ISGeneralFilter - Remove all indices outside of [start, end) from an `ISGENERAL`
805: Collective
807: Input Parameters:
808: + is - the index set
809: . start - the lowest index kept
810: - end - one more than the highest index kept, `start` $\le$ `end`
812: Level: beginner
814: .seealso: [](sec_scatter), `IS`, `ISGENERAL`, `ISCreateGeneral()`, `ISGeneralSetIndices()`
815: @*/
816: PetscErrorCode ISGeneralFilter(IS is, PetscInt start, PetscInt end)
817: {
818: PetscFunctionBegin;
820: 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);
821: PetscCall(ISClearInfoCache(is, PETSC_FALSE));
822: PetscUseMethod(is, "ISGeneralFilter_C", (IS, PetscInt, PetscInt), (is, start, end));
823: PetscFunctionReturn(PETSC_SUCCESS);
824: }
826: PETSC_INTERN PetscErrorCode ISCreate_General(IS is)
827: {
828: IS_General *sub;
830: PetscFunctionBegin;
831: PetscCall(PetscNew(&sub));
832: is->data = (void *)sub;
833: is->ops[0] = myops;
834: PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISGeneralSetIndices_C", ISGeneralSetIndices_General));
835: PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISGeneralSetIndicesFromMask_C", ISGeneralSetIndicesFromMask_General));
836: PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISGeneralFilter_C", ISGeneralFilter_General));
837: PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISShift_C", ISShift_General));
838: PetscFunctionReturn(PETSC_SUCCESS);
839: }