1 /*
2 Provides the functions for index sets (IS) defined by a list of integers.
3 These are for blocks of data, each block is indicated with a single integer.
4 */
5 #include <petsc/private/isimpl.h> /*I "petscis.h" I*/
6 #include <petscviewer.h>
7
8 typedef struct {
9 PetscBool sorted; /* are the blocks sorted? */
10 PetscBool allocated; /* did we allocate the index array ourselves? */
11 PetscInt *idx;
12 } IS_Block;
13
ISDestroy_Block(IS is)14 static PetscErrorCode ISDestroy_Block(IS is)
15 {
16 IS_Block *sub = (IS_Block *)is->data;
17
18 PetscFunctionBegin;
19 if (sub->allocated) PetscCall(PetscFree(sub->idx));
20 PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockSetIndices_C", NULL));
21 PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockGetIndices_C", NULL));
22 PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockRestoreIndices_C", NULL));
23 PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockGetSize_C", NULL));
24 PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockGetLocalSize_C", NULL));
25 PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISShift_C", NULL));
26 PetscCall(PetscFree(is->data));
27 PetscFunctionReturn(PETSC_SUCCESS);
28 }
29
ISLocate_Block(IS is,PetscInt key,PetscInt * location)30 static PetscErrorCode ISLocate_Block(IS is, PetscInt key, PetscInt *location)
31 {
32 IS_Block *sub = (IS_Block *)is->data;
33 PetscInt numIdx, i, bs, bkey, mkey;
34 PetscBool sorted;
35
36 PetscFunctionBegin;
37 PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
38 PetscCall(PetscLayoutGetLocalSize(is->map, &numIdx));
39 numIdx /= bs;
40 bkey = key / bs;
41 mkey = key % bs;
42 if (mkey < 0) {
43 bkey--;
44 mkey += bs;
45 }
46 PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, &sorted));
47 if (sorted) {
48 PetscCall(PetscFindInt(bkey, numIdx, sub->idx, location));
49 } else {
50 const PetscInt *idx = sub->idx;
51
52 *location = -1;
53 for (i = 0; i < numIdx; i++) {
54 if (idx[i] == bkey) {
55 *location = i;
56 break;
57 }
58 }
59 }
60 if (*location >= 0) *location = *location * bs + mkey;
61 PetscFunctionReturn(PETSC_SUCCESS);
62 }
63
ISGetIndices_Block(IS in,const PetscInt * idx[])64 static PetscErrorCode ISGetIndices_Block(IS in, const PetscInt *idx[])
65 {
66 IS_Block *sub = (IS_Block *)in->data;
67 PetscInt i, j, k, bs, n, *ii, *jj;
68
69 PetscFunctionBegin;
70 PetscCall(PetscLayoutGetBlockSize(in->map, &bs));
71 PetscCall(PetscLayoutGetLocalSize(in->map, &n));
72 n /= bs;
73 if (bs == 1) *idx = sub->idx;
74 else {
75 if (n) {
76 PetscCall(PetscMalloc1(bs * n, &jj));
77 *idx = jj;
78 k = 0;
79 ii = sub->idx;
80 for (i = 0; i < n; i++)
81 for (j = 0; j < bs; j++) jj[k++] = bs * ii[i] + j;
82 } else {
83 /* do not malloc for zero size because F90Array1dCreate() inside ISRestoreArray() does not keep array when zero length array */
84 *idx = NULL;
85 }
86 }
87 PetscFunctionReturn(PETSC_SUCCESS);
88 }
89
ISRestoreIndices_Block(IS is,const PetscInt * idx[])90 static PetscErrorCode ISRestoreIndices_Block(IS is, const PetscInt *idx[])
91 {
92 IS_Block *sub = (IS_Block *)is->data;
93 PetscInt bs;
94
95 PetscFunctionBegin;
96 PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
97 if (bs != 1) {
98 PetscCall(PetscFree(*(void **)idx));
99 } else {
100 /* F90Array1dCreate() inside ISRestoreArray() does not keep array when zero length array */
101 PetscCheck(is->map->n <= 0 || *idx == sub->idx, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Must restore with value from ISGetIndices()");
102 }
103 PetscFunctionReturn(PETSC_SUCCESS);
104 }
105
ISInvertPermutation_Block(IS is,PetscInt nlocal,IS * isout)106 static PetscErrorCode ISInvertPermutation_Block(IS is, PetscInt nlocal, IS *isout)
107 {
108 IS_Block *sub = (IS_Block *)is->data;
109 PetscInt i, *ii, bs, n, *idx = sub->idx;
110 PetscMPIInt size;
111
112 PetscFunctionBegin;
113 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is), &size));
114 PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
115 PetscCall(PetscLayoutGetLocalSize(is->map, &n));
116 n /= bs;
117 PetscCheck(size == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "No inversion written yet for block IS");
118 PetscCall(PetscMalloc1(n, &ii));
119 for (i = 0; i < n; i++) ii[idx[i]] = i;
120 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, n, ii, PETSC_OWN_POINTER, isout));
121 PetscCall(ISSetPermutation(*isout));
122 PetscFunctionReturn(PETSC_SUCCESS);
123 }
124
ISView_Block(IS is,PetscViewer viewer)125 static PetscErrorCode ISView_Block(IS is, PetscViewer viewer)
126 {
127 IS_Block *sub = (IS_Block *)is->data;
128 PetscInt i, bs, n, *idx = sub->idx;
129 PetscBool isascii, ibinary;
130
131 PetscFunctionBegin;
132 PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
133 PetscCall(PetscLayoutGetLocalSize(is->map, &n));
134 n /= bs;
135 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
136 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &ibinary));
137 if (isascii) {
138 PetscViewerFormat fmt;
139
140 PetscCall(PetscViewerGetFormat(viewer, &fmt));
141 if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
142 IS ist;
143 const char *name;
144 const PetscInt *idx;
145 PetscInt n;
146
147 PetscCall(PetscObjectGetName((PetscObject)is, &name));
148 PetscCall(ISGetLocalSize(is, &n));
149 PetscCall(ISGetIndices(is, &idx));
150 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is), n, idx, PETSC_USE_POINTER, &ist));
151 PetscCall(PetscObjectSetName((PetscObject)ist, name));
152 PetscCall(ISView(ist, viewer));
153 PetscCall(ISDestroy(&ist));
154 PetscCall(ISRestoreIndices(is, &idx));
155 } else {
156 PetscBool isperm;
157
158 PetscCall(ISGetInfo(is, IS_PERMUTATION, IS_GLOBAL, PETSC_FALSE, &isperm));
159 if (isperm) PetscCall(PetscViewerASCIIPrintf(viewer, "Block Index set is permutation\n"));
160 PetscCall(PetscViewerASCIIPushSynchronized(viewer));
161 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Block size %" PetscInt_FMT "\n", bs));
162 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Number of block indices in set %" PetscInt_FMT "\n", n));
163 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "The first indices of each block are\n"));
164 for (i = 0; i < n; i++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Block %" PetscInt_FMT " Index %" PetscInt_FMT "\n", i, idx[i]));
165 PetscCall(PetscViewerFlush(viewer));
166 PetscCall(PetscViewerASCIIPopSynchronized(viewer));
167 }
168 } else if (ibinary) PetscCall(ISView_Binary(is, viewer));
169 PetscFunctionReturn(PETSC_SUCCESS);
170 }
171
ISSort_Block(IS is)172 static PetscErrorCode ISSort_Block(IS is)
173 {
174 IS_Block *sub = (IS_Block *)is->data;
175 PetscInt bs, n;
176
177 PetscFunctionBegin;
178 PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
179 PetscCall(PetscLayoutGetLocalSize(is->map, &n));
180 PetscCall(PetscIntSortSemiOrdered(n / bs, sub->idx));
181 PetscFunctionReturn(PETSC_SUCCESS);
182 }
183
ISSortRemoveDups_Block(IS is)184 static PetscErrorCode ISSortRemoveDups_Block(IS is)
185 {
186 IS_Block *sub = (IS_Block *)is->data;
187 PetscInt bs, n, nb;
188 PetscBool sorted;
189
190 PetscFunctionBegin;
191 PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
192 PetscCall(PetscLayoutGetLocalSize(is->map, &n));
193 nb = n / bs;
194 PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, &sorted));
195 if (sorted) {
196 PetscCall(PetscSortedRemoveDupsInt(&nb, sub->idx));
197 } else {
198 PetscCall(PetscSortRemoveDupsInt(&nb, sub->idx));
199 }
200 PetscCall(PetscLayoutDestroy(&is->map));
201 PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), nb * bs, PETSC_DECIDE, bs, &is->map));
202 PetscFunctionReturn(PETSC_SUCCESS);
203 }
204
ISSorted_Block(IS is,PetscBool * flg)205 static PetscErrorCode ISSorted_Block(IS is, PetscBool *flg)
206 {
207 PetscFunctionBegin;
208 PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, flg));
209 PetscFunctionReturn(PETSC_SUCCESS);
210 }
211
ISSortedLocal_Block(IS is,PetscBool * flg)212 static PetscErrorCode ISSortedLocal_Block(IS is, PetscBool *flg)
213 {
214 IS_Block *sub = (IS_Block *)is->data;
215 PetscInt n, bs, i, *idx;
216
217 PetscFunctionBegin;
218 PetscCall(PetscLayoutGetLocalSize(is->map, &n));
219 PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
220 n /= bs;
221 idx = sub->idx;
222 for (i = 1; i < n; i++)
223 if (idx[i] < idx[i - 1]) break;
224 if (i < n) *flg = PETSC_FALSE;
225 else *flg = PETSC_TRUE;
226 PetscFunctionReturn(PETSC_SUCCESS);
227 }
228
ISUniqueLocal_Block(IS is,PetscBool * flg)229 static PetscErrorCode ISUniqueLocal_Block(IS is, PetscBool *flg)
230 {
231 IS_Block *sub = (IS_Block *)is->data;
232 PetscInt n, bs, i, *idx, *idxcopy = NULL;
233 PetscBool sortedLocal;
234
235 PetscFunctionBegin;
236 PetscCall(PetscLayoutGetLocalSize(is->map, &n));
237 PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
238 n /= bs;
239 idx = sub->idx;
240 PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, &sortedLocal));
241 if (!sortedLocal) {
242 PetscCall(PetscMalloc1(n, &idxcopy));
243 PetscCall(PetscArraycpy(idxcopy, idx, n));
244 PetscCall(PetscIntSortSemiOrdered(n, idxcopy));
245 idx = idxcopy;
246 }
247 for (i = 1; i < n; i++)
248 if (idx[i] == idx[i - 1]) break;
249 if (i < n) *flg = PETSC_FALSE;
250 else *flg = PETSC_TRUE;
251 PetscCall(PetscFree(idxcopy));
252 PetscFunctionReturn(PETSC_SUCCESS);
253 }
254
ISPermutationLocal_Block(IS is,PetscBool * flg)255 static PetscErrorCode ISPermutationLocal_Block(IS is, PetscBool *flg)
256 {
257 IS_Block *sub = (IS_Block *)is->data;
258 PetscInt n, bs, i, *idx, *idxcopy = NULL;
259 PetscBool sortedLocal;
260
261 PetscFunctionBegin;
262 PetscCall(PetscLayoutGetLocalSize(is->map, &n));
263 PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
264 n /= bs;
265 idx = sub->idx;
266 PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, &sortedLocal));
267 if (!sortedLocal) {
268 PetscCall(PetscMalloc1(n, &idxcopy));
269 PetscCall(PetscArraycpy(idxcopy, idx, n));
270 PetscCall(PetscIntSortSemiOrdered(n, idxcopy));
271 idx = idxcopy;
272 }
273 for (i = 0; i < n; i++)
274 if (idx[i] != i) break;
275 if (i < n) *flg = PETSC_FALSE;
276 else *flg = PETSC_TRUE;
277 PetscCall(PetscFree(idxcopy));
278 PetscFunctionReturn(PETSC_SUCCESS);
279 }
280
ISIntervalLocal_Block(IS is,PetscBool * flg)281 static PetscErrorCode ISIntervalLocal_Block(IS is, PetscBool *flg)
282 {
283 IS_Block *sub = (IS_Block *)is->data;
284 PetscInt n, bs, i, *idx;
285
286 PetscFunctionBegin;
287 PetscCall(PetscLayoutGetLocalSize(is->map, &n));
288 PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
289 n /= bs;
290 idx = sub->idx;
291 for (i = 1; i < n; i++)
292 if (idx[i] != idx[i - 1] + 1) break;
293 if (i < n) *flg = PETSC_FALSE;
294 else *flg = PETSC_TRUE;
295 PetscFunctionReturn(PETSC_SUCCESS);
296 }
297
ISDuplicate_Block(IS is,IS * newIS)298 static PetscErrorCode ISDuplicate_Block(IS is, IS *newIS)
299 {
300 IS_Block *sub = (IS_Block *)is->data;
301 PetscInt bs, n;
302
303 PetscFunctionBegin;
304 PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
305 PetscCall(PetscLayoutGetLocalSize(is->map, &n));
306 n /= bs;
307 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)is), bs, n, sub->idx, PETSC_COPY_VALUES, newIS));
308 PetscFunctionReturn(PETSC_SUCCESS);
309 }
310
ISCopy_Block(IS is,IS isy)311 static PetscErrorCode ISCopy_Block(IS is, IS isy)
312 {
313 IS_Block *is_block = (IS_Block *)is->data, *isy_block = (IS_Block *)isy->data;
314 PetscInt bs, n;
315
316 PetscFunctionBegin;
317 PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
318 PetscCall(PetscLayoutGetLocalSize(is->map, &n));
319 PetscCall(PetscArraycpy(isy_block->idx, is_block->idx, n / bs));
320 PetscFunctionReturn(PETSC_SUCCESS);
321 }
322
ISOnComm_Block(IS is,MPI_Comm comm,PetscCopyMode mode,IS * newis)323 static PetscErrorCode ISOnComm_Block(IS is, MPI_Comm comm, PetscCopyMode mode, IS *newis)
324 {
325 IS_Block *sub = (IS_Block *)is->data;
326 PetscInt bs, n;
327
328 PetscFunctionBegin;
329 PetscCheck(mode != PETSC_OWN_POINTER, comm, PETSC_ERR_ARG_WRONG, "Cannot use PETSC_OWN_POINTER");
330 PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
331 PetscCall(PetscLayoutGetLocalSize(is->map, &n));
332 PetscCall(ISCreateBlock(comm, bs, n / bs, sub->idx, mode, newis));
333 PetscFunctionReturn(PETSC_SUCCESS);
334 }
335
ISShift_Block(IS is,PetscInt shift,IS isy)336 static PetscErrorCode ISShift_Block(IS is, PetscInt shift, IS isy)
337 {
338 IS_Block *isb = (IS_Block *)is->data;
339 IS_Block *isby = (IS_Block *)isy->data;
340 PetscInt i, n, bs;
341
342 PetscFunctionBegin;
343 PetscCall(PetscLayoutGetLocalSize(is->map, &n));
344 PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
345 shift /= bs;
346 for (i = 0; i < n / bs; i++) isby->idx[i] = isb->idx[i] + shift;
347 PetscFunctionReturn(PETSC_SUCCESS);
348 }
349
ISSetBlockSize_Block(IS is,PetscInt bs)350 static PetscErrorCode ISSetBlockSize_Block(IS is, PetscInt bs)
351 {
352 PetscFunctionBegin;
353 PetscCheck(is->map->bs <= 0 || bs == is->map->bs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot change blocksize %" PetscInt_FMT " (to %" PetscInt_FMT ") if ISType is ISBLOCK", is->map->bs, bs);
354 PetscCall(PetscLayoutSetBlockSize(is->map, bs));
355 PetscFunctionReturn(PETSC_SUCCESS);
356 }
357
ISToGeneral_Block(IS inis)358 static PetscErrorCode ISToGeneral_Block(IS inis)
359 {
360 IS_Block *sub = (IS_Block *)inis->data;
361 PetscInt bs, n;
362 const PetscInt *idx;
363
364 PetscFunctionBegin;
365 PetscCall(ISGetBlockSize(inis, &bs));
366 PetscCall(ISGetLocalSize(inis, &n));
367 PetscCall(ISGetIndices(inis, &idx));
368 if (bs == 1) {
369 PetscCopyMode mode = sub->allocated ? PETSC_OWN_POINTER : PETSC_USE_POINTER;
370 sub->allocated = PETSC_FALSE; /* prevent deallocation when changing the subtype*/
371 PetscCall(ISSetType(inis, ISGENERAL));
372 PetscCall(ISGeneralSetIndices(inis, n, idx, mode));
373 } else {
374 PetscCall(ISSetType(inis, ISGENERAL));
375 PetscCall(ISGeneralSetIndices(inis, n, idx, PETSC_OWN_POINTER));
376 }
377 PetscFunctionReturn(PETSC_SUCCESS);
378 }
379
380 // clang-format off
381 static const struct _ISOps myops = {
382 PetscDesignatedInitializer(getindices, ISGetIndices_Block),
383 PetscDesignatedInitializer(restoreindices, ISRestoreIndices_Block),
384 PetscDesignatedInitializer(invertpermutation, ISInvertPermutation_Block),
385 PetscDesignatedInitializer(sort, ISSort_Block),
386 PetscDesignatedInitializer(sortremovedups, ISSortRemoveDups_Block),
387 PetscDesignatedInitializer(sorted, ISSorted_Block),
388 PetscDesignatedInitializer(duplicate, ISDuplicate_Block),
389 PetscDesignatedInitializer(destroy, ISDestroy_Block),
390 PetscDesignatedInitializer(view, ISView_Block),
391 PetscDesignatedInitializer(load, ISLoad_Default),
392 PetscDesignatedInitializer(copy, ISCopy_Block),
393 PetscDesignatedInitializer(togeneral, ISToGeneral_Block),
394 PetscDesignatedInitializer(oncomm, ISOnComm_Block),
395 PetscDesignatedInitializer(setblocksize, ISSetBlockSize_Block),
396 PetscDesignatedInitializer(contiguous, NULL),
397 PetscDesignatedInitializer(locate, ISLocate_Block),
398 /* we can have specialized local routines for determining properties,
399 * but unless the block size is the same on each process (which is not guaranteed at
400 * the moment), then trying to do something specialized for global properties is too
401 * complicated */
402 PetscDesignatedInitializer(sortedlocal, ISSortedLocal_Block),
403 PetscDesignatedInitializer(sortedglobal, NULL),
404 PetscDesignatedInitializer(uniquelocal, ISUniqueLocal_Block),
405 PetscDesignatedInitializer(uniqueglobal, NULL),
406 PetscDesignatedInitializer(permlocal, ISPermutationLocal_Block),
407 PetscDesignatedInitializer(permglobal, NULL),
408 PetscDesignatedInitializer(intervallocal, ISIntervalLocal_Block),
409 PetscDesignatedInitializer(intervalglobal, NULL)
410 };
411 // clang-format on
412
413 /*@
414 ISBlockSetIndices - Set integers representing blocks of indices in an index set of `ISType` `ISBLOCK`
415
416 Collective
417
418 Input Parameters:
419 + is - the index set
420 . bs - number of elements in each block
421 . n - the length of the index set (the number of blocks)
422 . idx - the list of integers, one for each block, the integers contain the index of the first index of each block divided by the block size
423 - mode - see `PetscCopyMode`, only `PETSC_COPY_VALUES` and `PETSC_OWN_POINTER` are supported
424
425 Level: beginner
426
427 Notes:
428 When the communicator is not `MPI_COMM_SELF`, the operations on the
429 index sets, IS, are NOT conceptually the same as `MPI_Group` operations.
430 The index sets are then distributed sets of indices and thus certain operations
431 on them are collective.
432
433 The convenience routine `ISCreateBlock()` allows one to create the `IS` and provide the blocks in a single function call.
434
435 Example:
436 If you wish to index the values {0,1,4,5}, then use
437 a block size of 2 and idx of {0,2}.
438
439 .seealso: [](sec_scatter), `IS`, `ISCreateStride()`, `ISCreateGeneral()`, `ISAllGather()`, `ISCreateBlock()`, `ISBLOCK`, `ISGeneralSetIndices()`
440 @*/
ISBlockSetIndices(IS is,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode)441 PetscErrorCode ISBlockSetIndices(IS is, PetscInt bs, PetscInt n, const PetscInt idx[], PetscCopyMode mode)
442 {
443 PetscFunctionBegin;
444 PetscCall(ISClearInfoCache(is, PETSC_FALSE));
445 PetscUseMethod(is, "ISBlockSetIndices_C", (IS, PetscInt, PetscInt, const PetscInt[], PetscCopyMode), (is, bs, n, idx, mode));
446 PetscFunctionReturn(PETSC_SUCCESS);
447 }
448
ISBlockSetIndices_Block(IS is,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode)449 static PetscErrorCode ISBlockSetIndices_Block(IS is, PetscInt bs, PetscInt n, const PetscInt idx[], PetscCopyMode mode)
450 {
451 PetscInt i, min, max;
452 IS_Block *sub = (IS_Block *)is->data;
453 PetscLayout map;
454
455 PetscFunctionBegin;
456 PetscCheck(bs >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "block size < 1");
457 PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length < 0");
458 if (n) PetscAssertPointer(idx, 4);
459
460 PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), n * bs, is->map->N, bs, &map));
461 PetscCall(PetscLayoutDestroy(&is->map));
462 is->map = map;
463
464 if (sub->allocated) PetscCall(PetscFree(sub->idx));
465 if (mode == PETSC_COPY_VALUES) {
466 PetscCall(PetscMalloc1(n, &sub->idx));
467 PetscCall(PetscArraycpy(sub->idx, idx, n));
468 sub->allocated = PETSC_TRUE;
469 } else if (mode == PETSC_OWN_POINTER) {
470 sub->idx = (PetscInt *)idx;
471 sub->allocated = PETSC_TRUE;
472 } else if (mode == PETSC_USE_POINTER) {
473 sub->idx = (PetscInt *)idx;
474 sub->allocated = PETSC_FALSE;
475 }
476
477 if (n) {
478 min = max = idx[0];
479 for (i = 1; i < n; i++) {
480 if (idx[i] < min) min = idx[i];
481 if (idx[i] > max) max = idx[i];
482 }
483 is->min = bs * min;
484 is->max = bs * max + bs - 1;
485 } else {
486 is->min = PETSC_INT_MAX;
487 is->max = PETSC_INT_MIN;
488 }
489 PetscFunctionReturn(PETSC_SUCCESS);
490 }
491
492 /*@
493 ISCreateBlock - Creates a data structure for an index set containing
494 a list of integers. Each integer represents a fixed block size set of indices.
495
496 Collective
497
498 Input Parameters:
499 + comm - the MPI communicator
500 . bs - number of elements in each block
501 . n - the length of the index set (the number of blocks)
502 . idx - the list of integers, one for each block, the integers contain the index of the first entry of each block divided by the block size
503 - mode - see `PetscCopyMode`, only `PETSC_COPY_VALUES` and `PETSC_OWN_POINTER` are supported in this routine
504
505 Output Parameter:
506 . is - the new index set
507
508 Level: beginner
509
510 Notes:
511 When the communicator is not `MPI_COMM_SELF`, the operations on the
512 index sets, `IS`, are NOT conceptually the same as `MPI_Group` operations.
513 The index sets are then distributed sets of indices and thus certain operations
514 on them are collective.
515
516 The routine `ISBlockSetIndices()` can be used to provide the indices to a preexisting block `IS`
517
518 Example:
519 If you wish to index the values {0,1,6,7}, then use
520 a block size of 2 and idx of {0,3}.
521
522 .seealso: [](sec_scatter), `IS`, `ISCreateStride()`, `ISCreateGeneral()`, `ISAllGather()`, `ISBlockSetIndices()`, `ISBLOCK`, `ISGENERAL`
523 @*/
ISCreateBlock(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode,IS * is)524 PetscErrorCode ISCreateBlock(MPI_Comm comm, PetscInt bs, PetscInt n, const PetscInt idx[], PetscCopyMode mode, IS *is)
525 {
526 PetscFunctionBegin;
527 PetscAssertPointer(is, 6);
528 PetscCheck(bs >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "block size < 1");
529 PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length < 0");
530 if (n) PetscAssertPointer(idx, 4);
531
532 PetscCall(ISCreate(comm, is));
533 PetscCall(ISSetType(*is, ISBLOCK));
534 PetscCall(ISBlockSetIndices(*is, bs, n, idx, mode));
535 PetscFunctionReturn(PETSC_SUCCESS);
536 }
537
ISBlockGetIndices_Block(IS is,const PetscInt * idx[])538 static PetscErrorCode ISBlockGetIndices_Block(IS is, const PetscInt *idx[])
539 {
540 IS_Block *sub = (IS_Block *)is->data;
541
542 PetscFunctionBegin;
543 *idx = sub->idx;
544 PetscFunctionReturn(PETSC_SUCCESS);
545 }
546
ISBlockRestoreIndices_Block(IS is,const PetscInt * idx[])547 static PetscErrorCode ISBlockRestoreIndices_Block(IS is, const PetscInt *idx[])
548 {
549 PetscFunctionBegin;
550 PetscFunctionReturn(PETSC_SUCCESS);
551 }
552
553 /*@C
554 ISBlockGetIndices - Gets the indices associated with each block in an `ISBLOCK`
555
556 Not Collective
557
558 Input Parameter:
559 . is - the index set
560
561 Output Parameter:
562 . idx - the integer indices, one for each block and count of block not indices
563
564 Level: intermediate
565
566 Note:
567 Call `ISBlockRestoreIndices()` when you no longer need access to the indices
568
569 Fortran Note:
570 .vb
571 PetscInt, pointer :: idx(:)
572 .ve
573
574 .seealso: [](sec_scatter), `IS`, `ISBLOCK`, `ISGetIndices()`, `ISBlockRestoreIndices()`, `ISBlockSetIndices()`, `ISCreateBlock()`
575 @*/
ISBlockGetIndices(IS is,const PetscInt * idx[])576 PetscErrorCode ISBlockGetIndices(IS is, const PetscInt *idx[])
577 {
578 PetscFunctionBegin;
579 PetscUseMethod(is, "ISBlockGetIndices_C", (IS, const PetscInt *[]), (is, idx));
580 PetscFunctionReturn(PETSC_SUCCESS);
581 }
582
583 /*@C
584 ISBlockRestoreIndices - Restores the indices associated with each block in an `ISBLOCK` obtained with `ISBlockGetIndices()`
585
586 Not Collective
587
588 Input Parameter:
589 . is - the index set
590
591 Output Parameter:
592 . idx - the integer indices
593
594 Level: intermediate
595
596 Fortran Note:
597 .vb
598 PetscInt, pointer :: idx(:)
599 .ve
600
601 .seealso: [](sec_scatter), `IS`, `ISBLOCK`, `ISRestoreIndices()`, `ISBlockGetIndices()`
602 @*/
ISBlockRestoreIndices(IS is,const PetscInt * idx[])603 PetscErrorCode ISBlockRestoreIndices(IS is, const PetscInt *idx[])
604 {
605 PetscFunctionBegin;
606 PetscUseMethod(is, "ISBlockRestoreIndices_C", (IS, const PetscInt *[]), (is, idx));
607 PetscFunctionReturn(PETSC_SUCCESS);
608 }
609
610 /*@
611 ISBlockGetLocalSize - Returns the local number of blocks in the index set of `ISType` `ISBLOCK`
612
613 Not Collective
614
615 Input Parameter:
616 . is - the index set
617
618 Output Parameter:
619 . size - the local number of blocks
620
621 Level: intermediate
622
623 .seealso: [](sec_scatter), `IS`, `ISGetBlockSize()`, `ISBlockGetSize()`, `ISGetSize()`, `ISCreateBlock()`, `ISBLOCK`
624 @*/
ISBlockGetLocalSize(IS is,PetscInt * size)625 PetscErrorCode ISBlockGetLocalSize(IS is, PetscInt *size)
626 {
627 PetscFunctionBegin;
628 PetscUseMethod(is, "ISBlockGetLocalSize_C", (IS, PetscInt *), (is, size));
629 PetscFunctionReturn(PETSC_SUCCESS);
630 }
631
ISBlockGetLocalSize_Block(IS is,PetscInt * size)632 static PetscErrorCode ISBlockGetLocalSize_Block(IS is, PetscInt *size)
633 {
634 PetscInt bs, n;
635
636 PetscFunctionBegin;
637 PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
638 PetscCall(PetscLayoutGetLocalSize(is->map, &n));
639 *size = n / bs;
640 PetscFunctionReturn(PETSC_SUCCESS);
641 }
642
643 /*@
644 ISBlockGetSize - Returns the global number of blocks in parallel in the index set of `ISType` `ISBLOCK`
645
646 Not Collective
647
648 Input Parameter:
649 . is - the index set
650
651 Output Parameter:
652 . size - the global number of blocks
653
654 Level: intermediate
655
656 .seealso: [](sec_scatter), `IS`, `ISGetBlockSize()`, `ISBlockGetLocalSize()`, `ISGetSize()`, `ISCreateBlock()`, `ISBLOCK`
657 @*/
ISBlockGetSize(IS is,PetscInt * size)658 PetscErrorCode ISBlockGetSize(IS is, PetscInt *size)
659 {
660 PetscFunctionBegin;
661 PetscUseMethod(is, "ISBlockGetSize_C", (IS, PetscInt *), (is, size));
662 PetscFunctionReturn(PETSC_SUCCESS);
663 }
664
ISBlockGetSize_Block(IS is,PetscInt * size)665 static PetscErrorCode ISBlockGetSize_Block(IS is, PetscInt *size)
666 {
667 PetscInt bs, N;
668
669 PetscFunctionBegin;
670 PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
671 PetscCall(PetscLayoutGetSize(is->map, &N));
672 *size = N / bs;
673 PetscFunctionReturn(PETSC_SUCCESS);
674 }
675
ISCreate_Block(IS is)676 PETSC_INTERN PetscErrorCode ISCreate_Block(IS is)
677 {
678 IS_Block *sub;
679
680 PetscFunctionBegin;
681 PetscCall(PetscNew(&sub));
682 is->data = (void *)sub;
683 is->ops[0] = myops;
684 PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockSetIndices_C", ISBlockSetIndices_Block));
685 PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockGetIndices_C", ISBlockGetIndices_Block));
686 PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockRestoreIndices_C", ISBlockRestoreIndices_Block));
687 PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockGetSize_C", ISBlockGetSize_Block));
688 PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockGetLocalSize_C", ISBlockGetLocalSize_Block));
689 PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISShift_C", ISShift_Block));
690 PetscFunctionReturn(PETSC_SUCCESS);
691 }
692