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> /*I "petscis.h" I*/
5 #include <petsc/private/viewerhdf5impl.h>
6
ISDuplicate_General(IS is,IS * newIS)7 static PetscErrorCode ISDuplicate_General(IS is, IS *newIS)
8 {
9 IS_General *sub = (IS_General *)is->data;
10 PetscInt n, bs;
11
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 }
19
ISDestroy_General(IS is)20 static PetscErrorCode ISDestroy_General(IS is)
21 {
22 IS_General *is_general = (IS_General *)is->data;
23
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 }
33
ISCopy_General(IS is,IS isy)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;
38
39 PetscFunctionBegin;
40 PetscCall(PetscLayoutGetLocalSize(is->map, &n));
41 PetscCall(PetscArraycpy(isy_general->idx, is_general->idx, n));
42 PetscFunctionReturn(PETSC_SUCCESS);
43 }
44
ISShift_General(IS is,PetscInt shift,IS isy)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;
49
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 }
55
ISOnComm_General(IS is,MPI_Comm comm,PetscCopyMode mode,IS * newis)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;
60
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 }
67
ISSetBlockSize_General(IS is,PetscInt bs)68 static PetscErrorCode ISSetBlockSize_General(IS is, PetscInt bs)
69 {
70 PetscFunctionBegin;
71 PetscCall(PetscLayoutSetBlockSize(is->map, bs));
72 PetscFunctionReturn(PETSC_SUCCESS);
73 }
74
ISContiguousLocal_General(IS is,PetscInt gstart,PetscInt gend,PetscInt * start,PetscBool * contig)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;
79
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 }
98
ISLocate_General(IS is,PetscInt key,PetscInt * location)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;
104
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;
111
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 }
122
ISGetIndices_General(IS in,const PetscInt * idx[])123 static PetscErrorCode ISGetIndices_General(IS in, const PetscInt *idx[])
124 {
125 IS_General *sub = (IS_General *)in->data;
126
127 PetscFunctionBegin;
128 *idx = sub->idx;
129 PetscFunctionReturn(PETSC_SUCCESS);
130 }
131
ISRestoreIndices_General(IS in,const PetscInt * idx[])132 static PetscErrorCode ISRestoreIndices_General(IS in, const PetscInt *idx[])
133 {
134 IS_General *sub = (IS_General *)in->data;
135
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 }
141
ISInvertPermutation_General(IS is,PetscInt nlocal,IS * isout)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;
149
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 }
182
183 /*
184 FindRun_Private - Determine whether a run exists in the sequence
185
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
191
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
196
197 Note:
198 We define a run as an arithmetic progression of at least `minRunLen` values, where zero is a valid step.
199 */
ISFindRun_Private(const PetscInt indices[],PetscInt len,PetscInt off,PetscInt minRunLen,PetscInt * runLen,PetscInt * step,PetscInt * val)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;
203
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 }
220
ISGeneralCheckCompress(IS is,PetscBool * compress)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;
227
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;
236
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(MPIU_Allreduce(&lcompress, compress, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
245 PetscFunctionReturn(PETSC_SUCCESS);
246 }
247
248 #if defined(PETSC_HAVE_HDF5)
249 // Run length encode the IS
ISGeneralCompress(IS is,IS * cis)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;
259
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;
268
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 }
292
ISView_General_HDF5(IS is,PetscViewer viewer)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;
308
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));
315
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));
338
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));
360
361 #if defined(PETSC_USE_64BIT_INDICES)
362 inttype = H5T_NATIVE_LLONG;
363 #else
364 inttype = H5T_NATIVE_INT;
365 #endif
366
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));
373
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));
381
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 }
400
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 }
421
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));
426
427 /* Close/release resources */
428 PetscCallHDF5(H5Gclose, (group));
429 PetscCallHDF5(H5Sclose, (filespace));
430 PetscCallHDF5(H5Sclose, (memspace));
431 PetscCallHDF5(H5Dclose, (dset_id));
432
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 }
437
ISView_General_HDF5_Compressed(IS is,PetscViewer viewer)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;
444
445 PetscFunctionBegin;
446 PetscCall(PetscViewerHDF5IsTimestepping(viewer, ×tepping));
447 PetscCheck(!timestepping, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_WRONG, "Timestepping not supported for compressed writing");
448
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));
453
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
459
ISView_General(IS is,PetscViewer viewer)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 isascii, isbinary, ishdf5, compress;
465
466 PetscFunctionBegin;
467 PetscCall(PetscLayoutGetLocalSize(is->map, &n));
468 PetscCall(ISGeneralCheckCompress(is, &compress));
469 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
470 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
471 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
472 if (isascii) {
473 MPI_Comm comm;
474 PetscMPIInt rank, size;
475 PetscViewerFormat fmt;
476 PetscBool isperm;
477
478 PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
479 PetscCallMPI(MPI_Comm_rank(comm, &rank));
480 PetscCallMPI(MPI_Comm_size(comm, &size));
481
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;
489
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;
496
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;
504
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;
511
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;
524
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 }
532
ISSort_General(IS is)533 static PetscErrorCode ISSort_General(IS is)
534 {
535 IS_General *sub = (IS_General *)is->data;
536 PetscInt n;
537
538 PetscFunctionBegin;
539 PetscCall(PetscLayoutGetLocalSize(is->map, &n));
540 PetscCall(PetscIntSortSemiOrdered(n, sub->idx));
541 PetscFunctionReturn(PETSC_SUCCESS);
542 }
543
ISSortRemoveDups_General(IS is)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;
550
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 }
564
ISSorted_General(IS is,PetscBool * flg)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 }
571
ISToGeneral_General(IS is)572 static PetscErrorCode ISToGeneral_General(IS is)
573 {
574 PetscFunctionBegin;
575 PetscFunctionReturn(PETSC_SUCCESS);
576 }
577
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
608
ISSetUp_General(IS is)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;
614
615 PetscFunctionBegin;
616 PetscCall(PetscLayoutGetLocalSize(is->map, &n));
617
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 }
632
633 /*@
634 ISCreateGeneral - Creates a data structure for an index set containing a list of integers.
635
636 Collective
637
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.
643
644 Output Parameter:
645 . is - the new index set
646
647 Level: beginner
648
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.
654
655 Use `ISGeneralSetIndices()` to provide indices to an already existing `IS` of `ISType` `ISGENERAL`
656
657 .seealso: [](sec_scatter), `IS`, `ISGENERAL`, `ISCreateStride()`, `ISCreateBlock()`, `ISAllGather()`, `PETSC_COPY_VALUES`, `PETSC_OWN_POINTER`,
658 `PETSC_USE_POINTER`, `PetscCopyMode`, `ISGeneralSetIndicesFromMask()`
659 @*/
ISCreateGeneral(MPI_Comm comm,PetscInt n,const PetscInt idx[],PetscCopyMode mode,IS * is)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 }
668
669 /*@
670 ISGeneralSetIndices - Sets the indices for an `ISGENERAL` index set
671
672 Logically Collective
673
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.
679
680 Level: beginner
681
682 Note:
683 Use `ISCreateGeneral()` to create the `IS` and set its indices in a single function call
684
685 .seealso: [](sec_scatter), `IS`, `ISBLOCK`, `ISCreateGeneral()`, `ISGeneralSetIndicesFromMask()`, `ISBlockSetIndices()`, `ISGENERAL`, `PetscCopyMode`
686 @*/
ISGeneralSetIndices(IS is,PetscInt n,const PetscInt idx[],PetscCopyMode mode)687 PetscErrorCode ISGeneralSetIndices(IS is, PetscInt n, const PetscInt idx[], PetscCopyMode mode)
688 {
689 PetscFunctionBegin;
690 PetscValidHeaderSpecific(is, IS_CLASSID, 1);
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 }
697
ISGeneralSetIndices_General(IS is,PetscInt n,const PetscInt idx[],PetscCopyMode mode)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;
702
703 PetscFunctionBegin;
704 PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), n, PETSC_DECIDE, is->map->bs, &map));
705 PetscCall(PetscLayoutDestroy(&is->map));
706 is->map = map;
707
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 }
720
721 PetscCall(ISSetUp_General(is));
722 PetscCall(ISViewFromOptions(is, NULL, "-is_view"));
723 PetscFunctionReturn(PETSC_SUCCESS);
724 }
725
726 /*@
727 ISGeneralSetIndicesFromMask - Sets the indices for an `ISGENERAL` index set using a boolean mask
728
729 Collective
730
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
736
737 Level: beginner
738
739 Note:
740 The mask array may be freed by the user after this call.
741
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.
752
753 .seealso: [](sec_scatter), `IS`, `ISCreateGeneral()`, `ISGeneralSetIndices()`, `ISGENERAL`
754 @*/
ISGeneralSetIndicesFromMask(IS is,PetscInt rstart,PetscInt rend,const PetscBool mask[])755 PetscErrorCode ISGeneralSetIndicesFromMask(IS is, PetscInt rstart, PetscInt rend, const PetscBool mask[])
756 {
757 PetscFunctionBegin;
758 PetscValidHeaderSpecific(is, IS_CLASSID, 1);
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 }
765
ISGeneralSetIndicesFromMask_General(IS is,PetscInt rstart,PetscInt rend,const PetscBool mask[])766 static PetscErrorCode ISGeneralSetIndicesFromMask_General(IS is, PetscInt rstart, PetscInt rend, const PetscBool mask[])
767 {
768 PetscInt i, nidx;
769 PetscInt *idx;
770
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 }
784
ISGeneralFilter_General(IS is,PetscInt start,PetscInt end)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;
790
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 }
801
802 /*@
803 ISGeneralFilter - Remove all indices outside of [start, end) from an `ISGENERAL`
804
805 Collective
806
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`
811
812 Level: beginner
813
814 .seealso: [](sec_scatter), `IS`, `ISGENERAL`, `ISCreateGeneral()`, `ISGeneralSetIndices()`
815 @*/
ISGeneralFilter(IS is,PetscInt start,PetscInt end)816 PetscErrorCode ISGeneralFilter(IS is, PetscInt start, PetscInt end)
817 {
818 PetscFunctionBegin;
819 PetscValidHeaderSpecific(is, IS_CLASSID, 1);
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 }
825
ISCreate_General(IS is)826 PETSC_INTERN PetscErrorCode ISCreate_General(IS is)
827 {
828 IS_General *sub;
829
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 }
840