Actual source code: block.c


  2: /*
  3:      Provides the functions for index sets (IS) defined by a list of integers.
  4:    These are for blocks of data, each block is indicated with a single integer.
  5: */
  6: #include <petsc/private/isimpl.h>
  7: #include <petscviewer.h>

  9: typedef struct {
 10:   PetscBool sorted;    /* are the blocks sorted? */
 11:   PetscBool allocated; /* did we allocate the index array ourselves? */
 12:   PetscInt *idx;
 13: } IS_Block;

 15: static PetscErrorCode ISDestroy_Block(IS is)
 16: {
 17:   IS_Block *sub = (IS_Block *)is->data;

 19:   PetscFunctionBegin;
 20:   if (sub->allocated) PetscCall(PetscFree(sub->idx));
 21:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockSetIndices_C", NULL));
 22:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockGetIndices_C", NULL));
 23:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockRestoreIndices_C", NULL));
 24:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockGetSize_C", NULL));
 25:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockGetLocalSize_C", NULL));
 26:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISShift_C", NULL));
 27:   PetscCall(PetscFree(is->data));
 28:   PetscFunctionReturn(PETSC_SUCCESS);
 29: }

 31: static PetscErrorCode ISLocate_Block(IS is, PetscInt key, PetscInt *location)
 32: {
 33:   IS_Block *sub = (IS_Block *)is->data;
 34:   PetscInt  numIdx, i, bs, bkey, mkey;
 35:   PetscBool sorted;

 37:   PetscFunctionBegin;
 38:   PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
 39:   PetscCall(PetscLayoutGetLocalSize(is->map, &numIdx));
 40:   numIdx /= bs;
 41:   bkey = key / bs;
 42:   mkey = key % bs;
 43:   if (mkey < 0) {
 44:     bkey--;
 45:     mkey += bs;
 46:   }
 47:   PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, &sorted));
 48:   if (sorted) {
 49:     PetscCall(PetscFindInt(bkey, numIdx, sub->idx, location));
 50:   } else {
 51:     const PetscInt *idx = sub->idx;

 53:     *location = -1;
 54:     for (i = 0; i < numIdx; i++) {
 55:       if (idx[i] == bkey) {
 56:         *location = i;
 57:         break;
 58:       }
 59:     }
 60:   }
 61:   if (*location >= 0) *location = *location * bs + mkey;
 62:   PetscFunctionReturn(PETSC_SUCCESS);
 63: }

 65: static PetscErrorCode ISGetIndices_Block(IS in, const PetscInt *idx[])
 66: {
 67:   IS_Block *sub = (IS_Block *)in->data;
 68:   PetscInt  i, j, k, bs, n, *ii, *jj;

 70:   PetscFunctionBegin;
 71:   PetscCall(PetscLayoutGetBlockSize(in->map, &bs));
 72:   PetscCall(PetscLayoutGetLocalSize(in->map, &n));
 73:   n /= bs;
 74:   if (bs == 1) *idx = sub->idx;
 75:   else {
 76:     if (n) {
 77:       PetscCall(PetscMalloc1(bs * n, &jj));
 78:       *idx = jj;
 79:       k    = 0;
 80:       ii   = sub->idx;
 81:       for (i = 0; i < n; i++)
 82:         for (j = 0; j < bs; j++) jj[k++] = bs * ii[i] + j;
 83:     } else {
 84:       /* do not malloc for zero size because F90Array1dCreate() inside ISRestoreArrayF90() does not keep array when zero length array */
 85:       *idx = NULL;
 86:     }
 87:   }
 88:   PetscFunctionReturn(PETSC_SUCCESS);
 89: }

 91: static PetscErrorCode ISRestoreIndices_Block(IS is, const PetscInt *idx[])
 92: {
 93:   IS_Block *sub = (IS_Block *)is->data;
 94:   PetscInt  bs;

 96:   PetscFunctionBegin;
 97:   PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
 98:   if (bs != 1) {
 99:     PetscCall(PetscFree(*(void **)idx));
100:   } else {
101:     /* F90Array1dCreate() inside ISRestoreArrayF90() does not keep array when zero length array */
102:     PetscCheck(is->map->n <= 0 || *idx == sub->idx, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Must restore with value from ISGetIndices()");
103:   }
104:   PetscFunctionReturn(PETSC_SUCCESS);
105: }

107: static PetscErrorCode ISInvertPermutation_Block(IS is, PetscInt nlocal, IS *isout)
108: {
109:   IS_Block   *sub = (IS_Block *)is->data;
110:   PetscInt    i, *ii, bs, n, *idx = sub->idx;
111:   PetscMPIInt size;

113:   PetscFunctionBegin;
114:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is), &size));
115:   PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
116:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
117:   n /= bs;
118:   if (size == 1) {
119:     PetscCall(PetscMalloc1(n, &ii));
120:     for (i = 0; i < n; i++) ii[idx[i]] = i;
121:     PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, n, ii, PETSC_OWN_POINTER, isout));
122:     PetscCall(ISSetPermutation(*isout));
123:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No inversion written yet for block IS");
124:   PetscFunctionReturn(PETSC_SUCCESS);
125: }

127: static PetscErrorCode ISView_Block(IS is, PetscViewer viewer)
128: {
129:   IS_Block *sub = (IS_Block *)is->data;
130:   PetscInt  i, bs, n, *idx = sub->idx;
131:   PetscBool iascii, ibinary;

133:   PetscFunctionBegin;
134:   PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
135:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
136:   n /= bs;
137:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
138:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &ibinary));
139:   if (iascii) {
140:     PetscViewerFormat fmt;

142:     PetscCall(PetscViewerGetFormat(viewer, &fmt));
143:     if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
144:       IS              ist;
145:       const char     *name;
146:       const PetscInt *idx;
147:       PetscInt        n;

149:       PetscCall(PetscObjectGetName((PetscObject)is, &name));
150:       PetscCall(ISGetLocalSize(is, &n));
151:       PetscCall(ISGetIndices(is, &idx));
152:       PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is), n, idx, PETSC_USE_POINTER, &ist));
153:       PetscCall(PetscObjectSetName((PetscObject)ist, name));
154:       PetscCall(ISView(ist, viewer));
155:       PetscCall(ISDestroy(&ist));
156:       PetscCall(ISRestoreIndices(is, &idx));
157:     } else {
158:       PetscBool isperm;

160:       PetscCall(ISGetInfo(is, IS_PERMUTATION, IS_GLOBAL, PETSC_FALSE, &isperm));
161:       if (isperm) PetscCall(PetscViewerASCIIPrintf(viewer, "Block Index set is permutation\n"));
162:       PetscCall(PetscViewerASCIIPushSynchronized(viewer));
163:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Block size %" PetscInt_FMT "\n", bs));
164:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Number of block indices in set %" PetscInt_FMT "\n", n));
165:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "The first indices of each block are\n"));
166:       for (i = 0; i < n; i++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Block %" PetscInt_FMT " Index %" PetscInt_FMT "\n", i, idx[i]));
167:       PetscCall(PetscViewerFlush(viewer));
168:       PetscCall(PetscViewerASCIIPopSynchronized(viewer));
169:     }
170:   } else if (ibinary) PetscCall(ISView_Binary(is, viewer));
171:   PetscFunctionReturn(PETSC_SUCCESS);
172: }

174: static PetscErrorCode ISSort_Block(IS is)
175: {
176:   IS_Block *sub = (IS_Block *)is->data;
177:   PetscInt  bs, n;

179:   PetscFunctionBegin;
180:   PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
181:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
182:   PetscCall(PetscIntSortSemiOrdered(n / bs, sub->idx));
183:   PetscFunctionReturn(PETSC_SUCCESS);
184: }

186: static PetscErrorCode ISSortRemoveDups_Block(IS is)
187: {
188:   IS_Block *sub = (IS_Block *)is->data;
189:   PetscInt  bs, n, nb;
190:   PetscBool sorted;

192:   PetscFunctionBegin;
193:   PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
194:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
195:   nb = n / bs;
196:   PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, &sorted));
197:   if (sorted) {
198:     PetscCall(PetscSortedRemoveDupsInt(&nb, sub->idx));
199:   } else {
200:     PetscCall(PetscSortRemoveDupsInt(&nb, sub->idx));
201:   }
202:   PetscCall(PetscLayoutDestroy(&is->map));
203:   PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), nb * bs, PETSC_DECIDE, bs, &is->map));
204:   PetscFunctionReturn(PETSC_SUCCESS);
205: }

207: static PetscErrorCode ISSorted_Block(IS is, PetscBool *flg)
208: {
209:   PetscFunctionBegin;
210:   PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, flg));
211:   PetscFunctionReturn(PETSC_SUCCESS);
212: }

214: static PetscErrorCode ISSortedLocal_Block(IS is, PetscBool *flg)
215: {
216:   IS_Block *sub = (IS_Block *)is->data;
217:   PetscInt  n, bs, i, *idx;

219:   PetscFunctionBegin;
220:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
221:   PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
222:   n /= bs;
223:   idx = sub->idx;
224:   for (i = 1; i < n; i++)
225:     if (idx[i] < idx[i - 1]) break;
226:   if (i < n) *flg = PETSC_FALSE;
227:   else *flg = PETSC_TRUE;
228:   PetscFunctionReturn(PETSC_SUCCESS);
229: }

231: static PetscErrorCode ISUniqueLocal_Block(IS is, PetscBool *flg)
232: {
233:   IS_Block *sub = (IS_Block *)is->data;
234:   PetscInt  n, bs, i, *idx, *idxcopy = NULL;
235:   PetscBool sortedLocal;

237:   PetscFunctionBegin;
238:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
239:   PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
240:   n /= bs;
241:   idx = sub->idx;
242:   PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, &sortedLocal));
243:   if (!sortedLocal) {
244:     PetscCall(PetscMalloc1(n, &idxcopy));
245:     PetscCall(PetscArraycpy(idxcopy, idx, n));
246:     PetscCall(PetscIntSortSemiOrdered(n, idxcopy));
247:     idx = idxcopy;
248:   }
249:   for (i = 1; i < n; i++)
250:     if (idx[i] == idx[i - 1]) break;
251:   if (i < n) *flg = PETSC_FALSE;
252:   else *flg = PETSC_TRUE;
253:   PetscCall(PetscFree(idxcopy));
254:   PetscFunctionReturn(PETSC_SUCCESS);
255: }

257: static PetscErrorCode ISPermutationLocal_Block(IS is, PetscBool *flg)
258: {
259:   IS_Block *sub = (IS_Block *)is->data;
260:   PetscInt  n, bs, i, *idx, *idxcopy = NULL;
261:   PetscBool sortedLocal;

263:   PetscFunctionBegin;
264:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
265:   PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
266:   n /= bs;
267:   idx = sub->idx;
268:   PetscCall(ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, &sortedLocal));
269:   if (!sortedLocal) {
270:     PetscCall(PetscMalloc1(n, &idxcopy));
271:     PetscCall(PetscArraycpy(idxcopy, idx, n));
272:     PetscCall(PetscIntSortSemiOrdered(n, idxcopy));
273:     idx = idxcopy;
274:   }
275:   for (i = 0; i < n; i++)
276:     if (idx[i] != i) break;
277:   if (i < n) *flg = PETSC_FALSE;
278:   else *flg = PETSC_TRUE;
279:   PetscCall(PetscFree(idxcopy));
280:   PetscFunctionReturn(PETSC_SUCCESS);
281: }

283: static PetscErrorCode ISIntervalLocal_Block(IS is, PetscBool *flg)
284: {
285:   IS_Block *sub = (IS_Block *)is->data;
286:   PetscInt  n, bs, i, *idx;

288:   PetscFunctionBegin;
289:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
290:   PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
291:   n /= bs;
292:   idx = sub->idx;
293:   for (i = 1; i < n; i++)
294:     if (idx[i] != idx[i - 1] + 1) break;
295:   if (i < n) *flg = PETSC_FALSE;
296:   else *flg = PETSC_TRUE;
297:   PetscFunctionReturn(PETSC_SUCCESS);
298: }

300: static PetscErrorCode ISDuplicate_Block(IS is, IS *newIS)
301: {
302:   IS_Block *sub = (IS_Block *)is->data;
303:   PetscInt  bs, n;

305:   PetscFunctionBegin;
306:   PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
307:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
308:   n /= bs;
309:   PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)is), bs, n, sub->idx, PETSC_COPY_VALUES, newIS));
310:   PetscFunctionReturn(PETSC_SUCCESS);
311: }

313: static PetscErrorCode ISCopy_Block(IS is, IS isy)
314: {
315:   IS_Block *is_block = (IS_Block *)is->data, *isy_block = (IS_Block *)isy->data;
316:   PetscInt  bs, n;

318:   PetscFunctionBegin;
319:   PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
320:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
321:   PetscCall(PetscArraycpy(isy_block->idx, is_block->idx, n / bs));
322:   PetscFunctionReturn(PETSC_SUCCESS);
323: }

325: static PetscErrorCode ISOnComm_Block(IS is, MPI_Comm comm, PetscCopyMode mode, IS *newis)
326: {
327:   IS_Block *sub = (IS_Block *)is->data;
328:   PetscInt  bs, n;

330:   PetscFunctionBegin;
331:   PetscCheck(mode != PETSC_OWN_POINTER, comm, PETSC_ERR_ARG_WRONG, "Cannot use PETSC_OWN_POINTER");
332:   PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
333:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
334:   PetscCall(ISCreateBlock(comm, bs, n / bs, sub->idx, mode, newis));
335:   PetscFunctionReturn(PETSC_SUCCESS);
336: }

338: static PetscErrorCode ISShift_Block(IS is, PetscInt shift, IS isy)
339: {
340:   IS_Block *isb  = (IS_Block *)is->data;
341:   IS_Block *isby = (IS_Block *)isy->data;
342:   PetscInt  i, n, bs;

344:   PetscFunctionBegin;
345:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
346:   PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
347:   shift /= bs;
348:   for (i = 0; i < n / bs; i++) isby->idx[i] = isb->idx[i] + shift;
349:   PetscFunctionReturn(PETSC_SUCCESS);
350: }

352: static PetscErrorCode ISSetBlockSize_Block(IS is, PetscInt bs)
353: {
354:   PetscFunctionBegin;
355:   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);
356:   PetscCall(PetscLayoutSetBlockSize(is->map, bs));
357:   PetscFunctionReturn(PETSC_SUCCESS);
358: }

360: static PetscErrorCode ISToGeneral_Block(IS inis)
361: {
362:   IS_Block       *sub = (IS_Block *)inis->data;
363:   PetscInt        bs, n;
364:   const PetscInt *idx;

366:   PetscFunctionBegin;
367:   PetscCall(ISGetBlockSize(inis, &bs));
368:   PetscCall(ISGetLocalSize(inis, &n));
369:   PetscCall(ISGetIndices(inis, &idx));
370:   if (bs == 1) {
371:     PetscCopyMode mode = sub->allocated ? PETSC_OWN_POINTER : PETSC_USE_POINTER;
372:     sub->allocated     = PETSC_FALSE; /* prevent deallocation when changing the subtype*/
373:     PetscCall(ISSetType(inis, ISGENERAL));
374:     PetscCall(ISGeneralSetIndices(inis, n, idx, mode));
375:   } else {
376:     PetscCall(ISSetType(inis, ISGENERAL));
377:     PetscCall(ISGeneralSetIndices(inis, n, idx, PETSC_OWN_POINTER));
378:   }
379:   PetscFunctionReturn(PETSC_SUCCESS);
380: }

382: static struct _ISOps myops = {ISGetIndices_Block, ISRestoreIndices_Block, ISInvertPermutation_Block, ISSort_Block, ISSortRemoveDups_Block, ISSorted_Block, ISDuplicate_Block, ISDestroy_Block, ISView_Block, ISLoad_Default, ISCopy_Block, ISToGeneral_Block, ISOnComm_Block, ISSetBlockSize_Block, NULL, ISLocate_Block,
383:                               /* we can have specialized local routines for determining properties,
384:                                 * but unless the block size is the same on each process (which is not guaranteed at
385:                                 * the moment), then trying to do something specialized for global properties is too
386:                                 * complicated */
387:                               ISSortedLocal_Block, NULL, ISUniqueLocal_Block, NULL, ISPermutationLocal_Block, NULL, ISIntervalLocal_Block, NULL};

389: /*@
390:    ISBlockSetIndices - Set integers representing blocks of indices in an index set of `ISType` `ISBLOCK`

392:    Collective

394:    Input Parameters:
395: +  is - the index set
396: .  bs - number of elements in each block
397: .   n - the length of the index set (the number of blocks)
398: .  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
399: -  mode - see `PetscCopyMode`, only `PETSC_COPY_VALUES` and `PETSC_OWN_POINTER` are supported

401:    Level: beginner

403:    Notes:
404:    When the communicator is not `MPI_COMM_SELF`, the operations on the
405:    index sets, IS, are NOT conceptually the same as `MPI_Group` operations.
406:    The index sets are then distributed sets of indices and thus certain operations
407:    on them are collective.

409:    The convenience routine `ISCreateBlock()` allows one to create the `IS` and provide the blocks in a single function call.
410:    Example:
411:    If you wish to index the values {0,1,4,5}, then use
412:    a block size of 2 and idx of {0,2}.

414: .seealso: [](sec_scatter), `IS`, `ISCreateStride()`, `ISCreateGeneral()`, `ISAllGather()`, `ISCreateBlock()`, `ISBLOCK`, `ISGeneralSetIndices()`
415: @*/
416: PetscErrorCode ISBlockSetIndices(IS is, PetscInt bs, PetscInt n, const PetscInt idx[], PetscCopyMode mode)
417: {
418:   PetscFunctionBegin;
419:   PetscCall(ISClearInfoCache(is, PETSC_FALSE));
420:   PetscUseMethod(is, "ISBlockSetIndices_C", (IS, PetscInt, PetscInt, const PetscInt[], PetscCopyMode), (is, bs, n, idx, mode));
421:   PetscFunctionReturn(PETSC_SUCCESS);
422: }

424: static PetscErrorCode ISBlockSetIndices_Block(IS is, PetscInt bs, PetscInt n, const PetscInt idx[], PetscCopyMode mode)
425: {
426:   PetscInt    i, min, max;
427:   IS_Block   *sub = (IS_Block *)is->data;
428:   PetscLayout map;

430:   PetscFunctionBegin;
431:   PetscCheck(bs >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "block size < 1");
432:   PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length < 0");

435:   PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), n * bs, is->map->N, bs, &map));
436:   PetscCall(PetscLayoutDestroy(&is->map));
437:   is->map = map;

439:   if (sub->allocated) PetscCall(PetscFree(sub->idx));
440:   if (mode == PETSC_COPY_VALUES) {
441:     PetscCall(PetscMalloc1(n, &sub->idx));
442:     PetscCall(PetscArraycpy(sub->idx, idx, n));
443:     sub->allocated = PETSC_TRUE;
444:   } else if (mode == PETSC_OWN_POINTER) {
445:     sub->idx       = (PetscInt *)idx;
446:     sub->allocated = PETSC_TRUE;
447:   } else if (mode == PETSC_USE_POINTER) {
448:     sub->idx       = (PetscInt *)idx;
449:     sub->allocated = PETSC_FALSE;
450:   }

452:   if (n) {
453:     min = max = idx[0];
454:     for (i = 1; i < n; i++) {
455:       if (idx[i] < min) min = idx[i];
456:       if (idx[i] > max) max = idx[i];
457:     }
458:     is->min = bs * min;
459:     is->max = bs * max + bs - 1;
460:   } else {
461:     is->min = PETSC_MAX_INT;
462:     is->max = PETSC_MIN_INT;
463:   }
464:   PetscFunctionReturn(PETSC_SUCCESS);
465: }

467: /*@
468:    ISCreateBlock - Creates a data structure for an index set containing
469:    a list of integers. Each integer represents a fixed block size set of indices.

471:    Collective

473:    Input Parameters:
474: +  comm - the MPI communicator
475: .  bs - number of elements in each block
476: .  n - the length of the index set (the number of blocks)
477: .  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
478: -  mode - see `PetscCopyMode`, only `PETSC_COPY_VALUES` and `PETSC_OWN_POINTER` are supported in this routine

480:    Output Parameter:
481: .  is - the new index set

483:    Level: beginner

485:    Notes:
486:    When the communicator is not `MPI_COMM_SELF`, the operations on the
487:    index sets, `IS`, are NOT conceptually the same as `MPI_Group` operations.
488:    The index sets are then distributed sets of indices and thus certain operations
489:    on them are collective.

491:    The routine `ISBlockSetIndices()` can be used to provide the indices to a preexisting block `IS`

493:    Example:
494:    If you wish to index the values {0,1,6,7}, then use
495:    a block size of 2 and idx of {0,3}.

497: .seealso: [](sec_scatter), `IS`, `ISCreateStride()`, `ISCreateGeneral()`, `ISAllGather()`, `ISBlockSetIndices()`, `ISBLOCK`, `ISGENERAL`
498: @*/
499: PetscErrorCode ISCreateBlock(MPI_Comm comm, PetscInt bs, PetscInt n, const PetscInt idx[], PetscCopyMode mode, IS *is)
500: {
501:   PetscFunctionBegin;
503:   PetscCheck(bs >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "block size < 1");
504:   PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length < 0");

507:   PetscCall(ISCreate(comm, is));
508:   PetscCall(ISSetType(*is, ISBLOCK));
509:   PetscCall(ISBlockSetIndices(*is, bs, n, idx, mode));
510:   PetscFunctionReturn(PETSC_SUCCESS);
511: }

513: static PetscErrorCode ISBlockGetIndices_Block(IS is, const PetscInt *idx[])
514: {
515:   IS_Block *sub = (IS_Block *)is->data;

517:   PetscFunctionBegin;
518:   *idx = sub->idx;
519:   PetscFunctionReturn(PETSC_SUCCESS);
520: }

522: static PetscErrorCode ISBlockRestoreIndices_Block(IS is, const PetscInt *idx[])
523: {
524:   PetscFunctionBegin;
525:   PetscFunctionReturn(PETSC_SUCCESS);
526: }

528: /*@C
529:    ISBlockGetIndices - Gets the indices associated with each block in an `ISBLOCK`

531:    Not Collective

533:    Input Parameter:
534: .  is - the index set

536:    Output Parameter:
537: .  idx - the integer indices, one for each block and count of block not indices

539:    Level: intermediate

541:    Note:
542:    Call `ISBlockRestoreIndices()` when you no longer need access to the indices

544: .seealso: [](sec_scatter), `IS`, `ISBLOCK`, `ISGetIndices()`, `ISBlockRestoreIndices()`, `ISBlockSetIndices()`, `ISCreateBlock()`
545: @*/
546: PetscErrorCode ISBlockGetIndices(IS is, const PetscInt *idx[])
547: {
548:   PetscFunctionBegin;
549:   PetscUseMethod(is, "ISBlockGetIndices_C", (IS, const PetscInt *[]), (is, idx));
550:   PetscFunctionReturn(PETSC_SUCCESS);
551: }

553: /*@C
554:    ISBlockRestoreIndices - Restores the indices associated with each block  in an `ISBLOCK` obtained with `ISBlockGetIndices()`

556:    Not Collective

558:    Input Parameter:
559: .  is - the index set

561:    Output Parameter:
562: .  idx - the integer indices

564:    Level: intermediate

566: .seealso: [](sec_scatter), `IS`, `ISBLOCK`, `ISRestoreIndices()`, `ISBlockGetIndices()`
567: @*/
568: PetscErrorCode ISBlockRestoreIndices(IS is, const PetscInt *idx[])
569: {
570:   PetscFunctionBegin;
571:   PetscUseMethod(is, "ISBlockRestoreIndices_C", (IS, const PetscInt *[]), (is, idx));
572:   PetscFunctionReturn(PETSC_SUCCESS);
573: }

575: /*@
576:    ISBlockGetLocalSize - Returns the local number of blocks in the index set of `ISType` `ISBLOCK`

578:    Not Collective

580:    Input Parameter:
581: .  is - the index set

583:    Output Parameter:
584: .  size - the local number of blocks

586:    Level: intermediate

588: .seealso: [](sec_scatter), `IS`, `ISGetBlockSize()`, `ISBlockGetSize()`, `ISGetSize()`, `ISCreateBlock()`, `ISBLOCK`
589: @*/
590: PetscErrorCode ISBlockGetLocalSize(IS is, PetscInt *size)
591: {
592:   PetscFunctionBegin;
593:   PetscUseMethod(is, "ISBlockGetLocalSize_C", (IS, PetscInt *), (is, size));
594:   PetscFunctionReturn(PETSC_SUCCESS);
595: }

597: static PetscErrorCode ISBlockGetLocalSize_Block(IS is, PetscInt *size)
598: {
599:   PetscInt bs, n;

601:   PetscFunctionBegin;
602:   PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
603:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
604:   *size = n / bs;
605:   PetscFunctionReturn(PETSC_SUCCESS);
606: }

608: /*@
609:    ISBlockGetSize - Returns the global number of blocks in parallel in the index set of `ISType` `ISBLOCK`

611:    Not Collective

613:    Input Parameter:
614: .  is - the index set

616:    Output Parameter:
617: .  size - the global number of blocks

619:    Level: intermediate

621: .seealso: [](sec_scatter), `IS`, `ISGetBlockSize()`, `ISBlockGetLocalSize()`, `ISGetSize()`, `ISCreateBlock()`, `ISBLOCK`
622: @*/
623: PetscErrorCode ISBlockGetSize(IS is, PetscInt *size)
624: {
625:   PetscFunctionBegin;
626:   PetscUseMethod(is, "ISBlockGetSize_C", (IS, PetscInt *), (is, size));
627:   PetscFunctionReturn(PETSC_SUCCESS);
628: }

630: static PetscErrorCode ISBlockGetSize_Block(IS is, PetscInt *size)
631: {
632:   PetscInt bs, N;

634:   PetscFunctionBegin;
635:   PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
636:   PetscCall(PetscLayoutGetSize(is->map, &N));
637:   *size = N / bs;
638:   PetscFunctionReturn(PETSC_SUCCESS);
639: }

641: PETSC_EXTERN PetscErrorCode ISCreate_Block(IS is)
642: {
643:   IS_Block *sub;

645:   PetscFunctionBegin;
646:   PetscCall(PetscNew(&sub));
647:   is->data = (void *)sub;
648:   PetscCall(PetscMemcpy(is->ops, &myops, sizeof(myops)));
649:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockSetIndices_C", ISBlockSetIndices_Block));
650:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockGetIndices_C", ISBlockGetIndices_Block));
651:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockRestoreIndices_C", ISBlockRestoreIndices_Block));
652:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockGetSize_C", ISBlockGetSize_Block));
653:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockGetLocalSize_C", ISBlockGetLocalSize_Block));
654:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISShift_C", ISShift_Block));
655:   PetscFunctionReturn(PETSC_SUCCESS);
656: }