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: const PetscInt *idx;
225: PetscInt n, off = 0;
227: PetscFunctionBegin;
228: *compress = PETSC_FALSE;
229: PetscCall(ISGetCompressOutput(is, compress));
230: if (!*compress) PetscFunctionReturn(PETSC_SUCCESS);
231: PetscCall(ISGetIndices(is, &idx));
232: PetscCall(ISGetLocalSize(is, &n));
233: while (off < n) {
234: PetscInt len;
236: PetscCall(ISFindRun_Private(idx, n, off, minRun, &len, NULL, NULL));
237: if (!len) {
238: *compress = PETSC_FALSE;
239: break;
240: }
241: off += len;
242: }
243: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, compress, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }
247: #if PetscDefined(HAVE_HDF5)
248: // Run length encode the IS
249: static PetscErrorCode ISGeneralCompress(IS is, IS *cis)
250: {
251: const PetscInt *idx;
252: const PetscInt minRun = 8;
253: PetscInt runs = 0;
254: PetscInt off = 0;
255: PetscInt *cidx, n, r;
256: const char *name;
257: MPI_Comm comm;
259: PetscFunctionBegin;
260: PetscCall(PetscObjectGetComm((PetscObject)is, &comm));
261: PetscCall(ISGetIndices(is, &idx));
262: PetscCall(ISGetLocalSize(is, &n));
263: if (!n) runs = 0;
264: // Count runs
265: while (off < n) {
266: PetscInt len;
268: PetscCall(ISFindRun_Private(idx, n, off, minRun, &len, NULL, NULL));
269: PetscCheck(len, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Compression failed to find a run at offset %" PetscInt_FMT, off);
270: ++runs;
271: off += len;
272: }
273: // Construct compressed set
274: PetscCall(PetscMalloc1(runs * 3, &cidx));
275: off = 0;
276: r = 0;
277: while (off < n) {
278: PetscCall(ISFindRun_Private(idx, n, off, minRun, &cidx[r * 3 + 0], &cidx[r * 3 + 1], &cidx[r * 3 + 2]));
279: PetscCheck(cidx[r * 3 + 0], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Compression failed to find a run at offset %" PetscInt_FMT, off);
280: off += cidx[r * 3 + 0];
281: ++r;
282: }
283: PetscCheck(r == runs, comm, PETSC_ERR_PLIB, "The number of runs %" PetscInt_FMT " != %" PetscInt_FMT " computed chunks", runs, r);
284: PetscCheck(off == n, comm, PETSC_ERR_PLIB, "The local size %" PetscInt_FMT " != %" PetscInt_FMT " total encoded size", n, off);
285: PetscCall(ISCreateGeneral(comm, runs * 3, cidx, PETSC_OWN_POINTER, cis));
286: PetscCall(PetscLayoutSetBlockSize((*cis)->map, 3));
287: PetscCall(PetscObjectGetName((PetscObject)is, &name));
288: PetscCall(PetscObjectSetName((PetscObject)*cis, name));
289: PetscFunctionReturn(PETSC_SUCCESS);
290: }
292: static PetscErrorCode ISView_General_HDF5(IS is, PetscViewer viewer)
293: {
294: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5 *)viewer->data;
295: hid_t filespace; /* file dataspace identifier */
296: hid_t chunkspace; /* chunk dataset property identifier */
297: hid_t dset_id; /* dataset identifier */
298: hid_t memspace; /* memory dataspace identifier */
299: hid_t inttype; /* int type (H5T_NATIVE_INT or H5T_NATIVE_LLONG) */
300: hid_t file_id, group;
301: hsize_t dim, maxDims[3], dims[3], chunkDims[3], count[3], offset[3];
302: PetscBool timestepping;
303: PetscInt bs, N, n, timestep = PETSC_INT_MIN, low;
304: hsize_t chunksize;
305: const PetscInt *ind;
306: const char *isname;
308: PetscFunctionBegin;
309: PetscCall(ISGetBlockSize(is, &bs));
310: bs = PetscMax(bs, 1); /* If N = 0, bs = 0 as well */
311: PetscCall(PetscViewerHDF5OpenGroup(viewer, NULL, &file_id, &group));
312: PetscCall(PetscViewerHDF5IsTimestepping(viewer, ×tepping));
313: if (timestepping) PetscCall(PetscViewerHDF5GetTimestep(viewer, ×tep));
315: /* Create the dataspace for the dataset.
316: *
317: * dims - holds the current dimensions of the dataset
318: *
319: * maxDims - holds the maximum dimensions of the dataset (unlimited
320: * for the number of time steps with the current dimensions for the
321: * other dimensions; so only additional time steps can be added).
322: *
323: * chunkDims - holds the size of a single time step (required to
324: * permit extending dataset).
325: */
326: dim = 0;
327: chunksize = 1;
328: if (timestep >= 0) {
329: dims[dim] = timestep + 1;
330: maxDims[dim] = H5S_UNLIMITED;
331: chunkDims[dim] = 1;
332: ++dim;
333: }
334: PetscCall(ISGetSize(is, &N));
335: PetscCall(ISGetLocalSize(is, &n));
336: PetscCall(PetscHDF5IntCast(N / bs, dims + dim));
338: maxDims[dim] = dims[dim];
339: chunkDims[dim] = PetscMax(1, dims[dim]);
340: chunksize *= chunkDims[dim];
341: ++dim;
342: if (bs >= 1) {
343: dims[dim] = bs;
344: maxDims[dim] = dims[dim];
345: chunkDims[dim] = dims[dim];
346: chunksize *= chunkDims[dim];
347: ++dim;
348: }
349: /* hdf5 chunks must be less than 4GB */
350: if (chunksize > PETSC_HDF5_MAX_CHUNKSIZE / 64) {
351: if (bs >= 1) {
352: 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));
353: 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));
354: } else {
355: chunkDims[dim - 1] = PETSC_HDF5_MAX_CHUNKSIZE / 64;
356: }
357: }
358: PetscCallHDF5Return(filespace, H5Screate_simple, ((int)dim, dims, maxDims));
360: inttype = PetscDefined(USE_64BIT_INDICES) ? H5T_NATIVE_LLONG : H5T_NATIVE_INT;
362: /* Create the dataset with default properties and close filespace */
363: PetscCall(PetscObjectGetName((PetscObject)is, &isname));
364: if (!H5Lexists(group, isname, H5P_DEFAULT)) {
365: /* Create chunk */
366: PetscCallHDF5Return(chunkspace, H5Pcreate, (H5P_DATASET_CREATE));
367: PetscCallHDF5(H5Pset_chunk, (chunkspace, (int)dim, chunkDims));
369: PetscCallHDF5Return(dset_id, H5Dcreate2, (group, isname, inttype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT));
370: PetscCallHDF5(H5Pclose, (chunkspace));
371: } else {
372: PetscCallHDF5Return(dset_id, H5Dopen2, (group, isname, H5P_DEFAULT));
373: PetscCallHDF5(H5Dset_extent, (dset_id, dims));
374: }
375: PetscCallHDF5(H5Sclose, (filespace));
377: /* Each process defines a dataset and writes it to the hyperslab in the file */
378: dim = 0;
379: if (timestep >= 0) {
380: count[dim] = 1;
381: ++dim;
382: }
383: PetscCall(PetscHDF5IntCast(n / bs, count + dim));
384: ++dim;
385: if (bs >= 1) {
386: count[dim] = bs;
387: ++dim;
388: }
389: if (n > 0 || H5_VERSION_GE(1, 10, 0)) {
390: PetscCallHDF5Return(memspace, H5Screate_simple, ((int)dim, count, NULL));
391: } else {
392: /* Can't create dataspace with zero for any dimension, so create null dataspace. */
393: PetscCallHDF5Return(memspace, H5Screate, (H5S_NULL));
394: }
396: /* Select hyperslab in the file */
397: PetscCall(PetscLayoutGetRange(is->map, &low, NULL));
398: dim = 0;
399: if (timestep >= 0) {
400: offset[dim] = timestep;
401: ++dim;
402: }
403: PetscCall(PetscHDF5IntCast(low / bs, offset + dim));
404: ++dim;
405: if (bs >= 1) {
406: offset[dim] = 0;
407: ++dim;
408: }
409: if (n > 0 || H5_VERSION_GE(1, 10, 0)) {
410: PetscCallHDF5Return(filespace, H5Dget_space, (dset_id));
411: PetscCallHDF5(H5Sselect_hyperslab, (filespace, H5S_SELECT_SET, offset, NULL, count, NULL));
412: } else {
413: /* Create null filespace to match null memspace. */
414: PetscCallHDF5Return(filespace, H5Screate, (H5S_NULL));
415: }
417: PetscCall(ISGetIndices(is, &ind));
418: PetscCallHDF5(H5Dwrite, (dset_id, inttype, memspace, filespace, hdf5->dxpl_id, ind));
419: PetscCallHDF5(H5Fflush, (file_id, H5F_SCOPE_GLOBAL));
420: PetscCall(ISRestoreIndices(is, &ind));
422: /* Close/release resources */
423: PetscCallHDF5(H5Gclose, (group));
424: PetscCallHDF5(H5Sclose, (filespace));
425: PetscCallHDF5(H5Sclose, (memspace));
426: PetscCallHDF5(H5Dclose, (dset_id));
428: if (timestepping) PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)is, "timestepping", PETSC_BOOL, ×tepping));
429: PetscCall(PetscInfo(is, "Wrote IS object with name %s\n", isname));
430: PetscFunctionReturn(PETSC_SUCCESS);
431: }
433: static PetscErrorCode ISView_General_HDF5_Compressed(IS is, PetscViewer viewer)
434: {
435: const PetscBool compressed = PETSC_TRUE;
436: const char *isname;
437: IS cis;
438: PetscBool timestepping;
440: PetscFunctionBegin;
441: PetscCall(PetscViewerHDF5IsTimestepping(viewer, ×tepping));
442: PetscCheck(!timestepping, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_WRONG, "Timestepping not supported for compressed writing");
444: PetscCall(ISGeneralCompress(is, &cis));
445: PetscCall(ISView_General_HDF5(cis, viewer));
446: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)cis, "compressed", PETSC_BOOL, &compressed));
447: PetscCall(ISDestroy(&cis));
449: PetscCall(PetscObjectGetName((PetscObject)is, &isname));
450: PetscCall(PetscInfo(is, "Wrote compressed IS object with name %s\n", isname));
451: PetscFunctionReturn(PETSC_SUCCESS);
452: }
453: #endif
455: static PetscErrorCode ISView_General(IS is, PetscViewer viewer)
456: {
457: IS_General *sub = (IS_General *)is->data;
458: PetscInt *idx = sub->idx, n;
459: PetscBool isascii, isbinary, ishdf5, compress;
461: PetscFunctionBegin;
462: PetscCall(PetscLayoutGetLocalSize(is->map, &n));
463: PetscCall(ISGeneralCheckCompress(is, &compress));
464: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
465: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
466: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
467: if (isascii) {
468: MPI_Comm comm;
469: PetscMPIInt rank, size;
470: PetscViewerFormat fmt;
471: PetscBool isperm;
473: PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
474: PetscCallMPI(MPI_Comm_rank(comm, &rank));
475: PetscCallMPI(MPI_Comm_size(comm, &size));
477: PetscCall(PetscViewerGetFormat(viewer, &fmt));
478: PetscCall(ISGetInfo(is, IS_PERMUTATION, IS_LOCAL, PETSC_FALSE, &isperm));
479: if (isperm && fmt != PETSC_VIEWER_ASCII_MATLAB) PetscCall(PetscViewerASCIIPrintf(viewer, "Index set is permutation\n"));
480: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
481: if (size > 1) {
482: if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
483: const char *name;
485: PetscCall(PetscObjectGetName((PetscObject)is, &name));
486: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%s_%d = [...\n", name, rank));
487: for (PetscInt i = 0; i < n; i++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%" PetscInt_FMT "\n", idx[i] + 1));
488: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "];\n"));
489: } else {
490: PetscInt st = 0;
492: if (fmt == PETSC_VIEWER_ASCII_INDEX) st = is->map->rstart;
493: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Number of indices in set %" PetscInt_FMT "\n", rank, n));
494: for (PetscInt i = 0; i < n; i++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] %" PetscInt_FMT " %" PetscInt_FMT "\n", rank, i + st, idx[i]));
495: }
496: } else {
497: if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
498: const char *name;
500: PetscCall(PetscObjectGetName((PetscObject)is, &name));
501: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%s = [...\n", name));
502: for (PetscInt i = 0; i < n; i++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%" PetscInt_FMT "\n", idx[i] + 1));
503: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "];\n"));
504: } else {
505: PetscInt st = 0;
507: if (fmt == PETSC_VIEWER_ASCII_INDEX) st = is->map->rstart;
508: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Number of indices in set %" PetscInt_FMT "\n", n));
509: for (PetscInt i = 0; i < n; i++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT "\n", i + st, idx[i]));
510: }
511: }
512: PetscCall(PetscViewerFlush(viewer));
513: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
514: } else if (isbinary) {
515: PetscCall(ISView_Binary(is, viewer));
516: } else if (ishdf5) {
517: #if PetscDefined(HAVE_HDF5)
518: PetscBool vcompress;
520: PetscCall(PetscViewerHDF5GetCompress(viewer, &vcompress));
521: if (vcompress && compress) PetscCall(ISView_General_HDF5_Compressed(is, viewer));
522: else PetscCall(ISView_General_HDF5(is, viewer));
523: #endif
524: }
525: PetscFunctionReturn(PETSC_SUCCESS);
526: }
528: static PetscErrorCode ISSort_General(IS is)
529: {
530: IS_General *sub = (IS_General *)is->data;
531: PetscInt n;
533: PetscFunctionBegin;
534: PetscCall(PetscLayoutGetLocalSize(is->map, &n));
535: PetscCall(PetscIntSortSemiOrdered(n, sub->idx));
536: PetscFunctionReturn(PETSC_SUCCESS);
537: }
539: static PetscErrorCode ISSortRemoveDups_General(IS is)
540: {
541: IS_General *sub = (IS_General *)is->data;
542: PetscLayout map;
543: PetscInt n;
544: PetscBool sorted;
546: PetscFunctionBegin;
547: PetscCall(PetscLayoutGetLocalSize(is->map, &n));
548: PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, &sorted));
549: if (sorted) {
550: PetscCall(PetscSortedRemoveDupsInt(&n, sub->idx));
551: } else {
552: PetscCall(PetscSortRemoveDupsInt(&n, sub->idx));
553: }
554: PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), n, PETSC_DECIDE, is->map->bs, &map));
555: PetscCall(PetscLayoutDestroy(&is->map));
556: is->map = map;
557: PetscFunctionReturn(PETSC_SUCCESS);
558: }
560: static PetscErrorCode ISSorted_General(IS is, PetscBool *flg)
561: {
562: PetscFunctionBegin;
563: PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, flg));
564: PetscFunctionReturn(PETSC_SUCCESS);
565: }
567: static PetscErrorCode ISToGeneral_General(IS is)
568: {
569: PetscFunctionBegin;
570: PetscFunctionReturn(PETSC_SUCCESS);
571: }
573: // clang-format off
574: static const struct _ISOps myops = {
575: PetscDesignatedInitializer(getindices, ISGetIndices_General),
576: PetscDesignatedInitializer(restoreindices, ISRestoreIndices_General),
577: PetscDesignatedInitializer(invertpermutation, ISInvertPermutation_General),
578: PetscDesignatedInitializer(sort, ISSort_General),
579: PetscDesignatedInitializer(sortremovedups, ISSortRemoveDups_General),
580: PetscDesignatedInitializer(sorted, ISSorted_General),
581: PetscDesignatedInitializer(duplicate, ISDuplicate_General),
582: PetscDesignatedInitializer(destroy, ISDestroy_General),
583: PetscDesignatedInitializer(view, ISView_General),
584: PetscDesignatedInitializer(load, ISLoad_Default),
585: PetscDesignatedInitializer(copy, ISCopy_General),
586: PetscDesignatedInitializer(togeneral, ISToGeneral_General),
587: PetscDesignatedInitializer(oncomm, ISOnComm_General),
588: PetscDesignatedInitializer(setblocksize, ISSetBlockSize_General),
589: PetscDesignatedInitializer(contiguous, ISContiguousLocal_General),
590: PetscDesignatedInitializer(locate, ISLocate_General),
591: /* no specializations of {sorted,unique,perm,interval}{local,global} because the default
592: * checks in ISGetInfo_XXX in index.c are exactly what we would do for ISGeneral */
593: PetscDesignatedInitializer(sortedlocal, NULL),
594: PetscDesignatedInitializer(sortedglobal, NULL),
595: PetscDesignatedInitializer(uniquelocal, NULL),
596: PetscDesignatedInitializer(uniqueglobal, NULL),
597: PetscDesignatedInitializer(permlocal, NULL),
598: PetscDesignatedInitializer(permglobal, NULL),
599: PetscDesignatedInitializer(intervallocal, NULL),
600: PetscDesignatedInitializer(intervalglobal, NULL)
601: };
602: // clang-format on
604: PETSC_INTERN PetscErrorCode ISSetUp_General(IS is)
605: {
606: IS_General *sub = (IS_General *)is->data;
607: const PetscInt *idx = sub->idx;
608: PetscInt n, i, min, max;
610: PetscFunctionBegin;
611: PetscCall(PetscLayoutGetLocalSize(is->map, &n));
613: if (n) {
614: min = max = idx[0];
615: for (i = 1; i < n; i++) {
616: if (idx[i] < min) min = idx[i];
617: if (idx[i] > max) max = idx[i];
618: }
619: is->min = min;
620: is->max = max;
621: } else {
622: is->min = PETSC_INT_MAX;
623: is->max = PETSC_INT_MIN;
624: }
625: PetscFunctionReturn(PETSC_SUCCESS);
626: }
628: /*@
629: ISCreateGeneral - Creates a data structure for an index set containing a list of integers.
631: Collective
633: Input Parameters:
634: + comm - the MPI communicator
635: . n - the length of the index set
636: . idx - the list of integers
637: - mode - `PETSC_COPY_VALUES`, `PETSC_OWN_POINTER`, or `PETSC_USE_POINTER`; see `PetscCopyMode` for meaning of this flag.
639: Output Parameter:
640: . is - the new index set
642: Level: beginner
644: Notes:
645: When the communicator is not `MPI_COMM_SELF`, the operations on IS are NOT
646: conceptually the same as `MPI_Group` operations. The `IS` are then
647: distributed sets of indices and thus certain operations on them are
648: collective.
650: Use `ISGeneralSetIndices()` to provide indices to an already existing `IS` of `ISType` `ISGENERAL`
652: .seealso: [](sec_scatter), `IS`, `ISGENERAL`, `ISCreateStride()`, `ISCreateBlock()`, `ISAllGather()`, `PETSC_COPY_VALUES`, `PETSC_OWN_POINTER`,
653: `PETSC_USE_POINTER`, `PetscCopyMode`, `ISGeneralSetIndicesFromMask()`
654: @*/
655: PetscErrorCode ISCreateGeneral(MPI_Comm comm, PetscInt n, const PetscInt idx[], PetscCopyMode mode, IS *is)
656: {
657: PetscFunctionBegin;
658: PetscCall(ISCreate(comm, is));
659: PetscCall(ISSetType(*is, ISGENERAL));
660: PetscCall(ISGeneralSetIndices(*is, n, idx, mode));
661: PetscFunctionReturn(PETSC_SUCCESS);
662: }
664: /*@
665: ISGeneralSetIndices - Sets the indices for an `ISGENERAL` index set
667: Logically Collective
669: Input Parameters:
670: + is - the index set
671: . n - the length of the index set
672: . idx - the list of integers
673: - mode - see `PetscCopyMode` for meaning of this flag.
675: Level: beginner
677: Note:
678: Use `ISCreateGeneral()` to create the `IS` and set its indices in a single function call
680: .seealso: [](sec_scatter), `IS`, `ISBLOCK`, `ISCreateGeneral()`, `ISGeneralSetIndicesFromMask()`, `ISBlockSetIndices()`, `ISGENERAL`, `PetscCopyMode`
681: @*/
682: PetscErrorCode ISGeneralSetIndices(IS is, PetscInt n, const PetscInt idx[], PetscCopyMode mode)
683: {
684: PetscFunctionBegin;
686: PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length < 0");
687: if (n) PetscAssertPointer(idx, 3);
688: PetscCall(ISClearInfoCache(is, PETSC_FALSE));
689: PetscUseMethod(is, "ISGeneralSetIndices_C", (IS, PetscInt, const PetscInt[], PetscCopyMode), (is, n, idx, mode));
690: PetscFunctionReturn(PETSC_SUCCESS);
691: }
693: static PetscErrorCode ISGeneralSetIndices_General(IS is, PetscInt n, const PetscInt idx[], PetscCopyMode mode)
694: {
695: PetscLayout map;
696: IS_General *sub = (IS_General *)is->data;
698: PetscFunctionBegin;
699: PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), n, PETSC_DECIDE, is->map->bs, &map));
700: PetscCall(PetscLayoutDestroy(&is->map));
701: is->map = map;
703: if (sub->allocated) PetscCall(PetscFree(sub->idx));
704: if (mode == PETSC_COPY_VALUES) {
705: PetscCall(PetscMalloc1(n, &sub->idx));
706: PetscCall(PetscArraycpy(sub->idx, idx, n));
707: sub->allocated = PETSC_TRUE;
708: } else if (mode == PETSC_OWN_POINTER) {
709: sub->idx = (PetscInt *)idx;
710: sub->allocated = PETSC_TRUE;
711: } else {
712: sub->idx = (PetscInt *)idx;
713: sub->allocated = PETSC_FALSE;
714: }
716: PetscCall(ISSetUp_General(is));
717: PetscCall(ISViewFromOptions(is, NULL, "-is_view"));
718: PetscFunctionReturn(PETSC_SUCCESS);
719: }
721: /*@
722: ISGeneralSetIndicesFromMask - Sets the indices for an `ISGENERAL` index set using a boolean mask
724: Collective
726: Input Parameters:
727: + is - the index set
728: . rstart - the range start index (inclusive)
729: . rend - the range end index (exclusive)
730: - mask - the boolean mask array of length rend-rstart, indices will be set for each `PETSC_TRUE` value in the array
732: Level: beginner
734: Note:
735: The mask array may be freed by the user after this call.
737: Example:
738: .vb
739: PetscBool mask[] = {PETSC_FALSE, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, PETSC_TRUE};
740: ISGeneralSetIndicesFromMask(is,10,15,mask);
741: .ve
742: will feed the `IS` with indices
743: .vb
744: {11, 14}
745: .ve
746: locally.
748: .seealso: [](sec_scatter), `IS`, `ISCreateGeneral()`, `ISGeneralSetIndices()`, `ISGENERAL`
749: @*/
750: PetscErrorCode ISGeneralSetIndicesFromMask(IS is, PetscInt rstart, PetscInt rend, const PetscBool mask[])
751: {
752: PetscFunctionBegin;
754: 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);
755: if (rend - rstart) PetscAssertPointer(mask, 4);
756: PetscCall(ISClearInfoCache(is, PETSC_FALSE));
757: PetscUseMethod(is, "ISGeneralSetIndicesFromMask_C", (IS, PetscInt, PetscInt, const PetscBool[]), (is, rstart, rend, mask));
758: PetscFunctionReturn(PETSC_SUCCESS);
759: }
761: static PetscErrorCode ISGeneralSetIndicesFromMask_General(IS is, PetscInt rstart, PetscInt rend, const PetscBool mask[])
762: {
763: PetscInt i, nidx;
764: PetscInt *idx;
766: PetscFunctionBegin;
767: for (i = 0, nidx = 0; i < rend - rstart; i++)
768: if (mask[i]) nidx++;
769: PetscCall(PetscMalloc1(nidx, &idx));
770: for (i = 0, nidx = 0; i < rend - rstart; i++) {
771: if (mask[i]) {
772: idx[nidx] = i + rstart;
773: nidx++;
774: }
775: }
776: PetscCall(ISGeneralSetIndices_General(is, nidx, idx, PETSC_OWN_POINTER));
777: PetscFunctionReturn(PETSC_SUCCESS);
778: }
780: static PetscErrorCode ISGeneralFilter_General(IS is, PetscInt start, PetscInt end)
781: {
782: IS_General *sub = (IS_General *)is->data;
783: PetscInt *idx = sub->idx, *idxnew;
784: PetscInt i, n = is->map->n, nnew = 0, o;
786: PetscFunctionBegin;
787: for (i = 0; i < n; ++i)
788: if (idx[i] >= start && idx[i] < end) nnew++;
789: PetscCall(PetscMalloc1(nnew, &idxnew));
790: for (o = 0, i = 0; i < n; i++) {
791: if (idx[i] >= start && idx[i] < end) idxnew[o++] = idx[i];
792: }
793: PetscCall(ISGeneralSetIndices_General(is, nnew, idxnew, PETSC_OWN_POINTER));
794: PetscFunctionReturn(PETSC_SUCCESS);
795: }
797: /*@
798: ISGeneralFilter - Remove all indices outside of [start, end) from an `ISGENERAL`
800: Collective
802: Input Parameters:
803: + is - the index set
804: . start - the lowest index kept
805: - end - one more than the highest index kept, `start` $\le$ `end`
807: Level: beginner
809: .seealso: [](sec_scatter), `IS`, `ISGENERAL`, `ISCreateGeneral()`, `ISGeneralSetIndices()`
810: @*/
811: PetscErrorCode ISGeneralFilter(IS is, PetscInt start, PetscInt end)
812: {
813: PetscFunctionBegin;
815: 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);
816: PetscCall(ISClearInfoCache(is, PETSC_FALSE));
817: PetscUseMethod(is, "ISGeneralFilter_C", (IS, PetscInt, PetscInt), (is, start, end));
818: PetscFunctionReturn(PETSC_SUCCESS);
819: }
821: PETSC_INTERN PetscErrorCode ISCreate_General(IS is)
822: {
823: IS_General *sub;
825: PetscFunctionBegin;
826: PetscCall(PetscNew(&sub));
827: is->data = (void *)sub;
828: is->ops[0] = myops;
829: PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISGeneralSetIndices_C", ISGeneralSetIndices_General));
830: PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISGeneralSetIndicesFromMask_C", ISGeneralSetIndicesFromMask_General));
831: PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISGeneralFilter_C", ISGeneralFilter_General));
832: PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISShift_C", ISShift_General));
833: PetscFunctionReturn(PETSC_SUCCESS);
834: }