Actual source code: block.c

  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>
  6: #include <petscviewer.h>

  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;

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

 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: }

 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;

 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;

 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: }

 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;

 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 ISRestoreArrayF90() does not keep array when zero length array */
 84:       *idx = NULL;
 85:     }
 86:   }
 87:   PetscFunctionReturn(PETSC_SUCCESS);
 88: }

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

 95:   PetscFunctionBegin;
 96:   PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
 97:   if (bs != 1) {
 98:     PetscCall(PetscFree(*(void **)idx));
 99:   } else {
100:     /* F90Array1dCreate() inside ISRestoreArrayF90() 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: }

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;

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:   if (size == 1) {
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:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No inversion written yet for block IS");
123:   PetscFunctionReturn(PETSC_SUCCESS);
124: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

381: // clang-format off
382: static const struct _ISOps myops = {
383:   PetscDesignatedInitializer(getindices, ISGetIndices_Block),
384:   PetscDesignatedInitializer(restoreindices, ISRestoreIndices_Block),
385:   PetscDesignatedInitializer(invertpermutation, ISInvertPermutation_Block),
386:   PetscDesignatedInitializer(sort, ISSort_Block),
387:   PetscDesignatedInitializer(sortremovedups, ISSortRemoveDups_Block),
388:   PetscDesignatedInitializer(sorted, ISSorted_Block),
389:   PetscDesignatedInitializer(duplicate, ISDuplicate_Block),
390:   PetscDesignatedInitializer(destroy, ISDestroy_Block),
391:   PetscDesignatedInitializer(view, ISView_Block),
392:   PetscDesignatedInitializer(load, ISLoad_Default),
393:   PetscDesignatedInitializer(copy, ISCopy_Block),
394:   PetscDesignatedInitializer(togeneral, ISToGeneral_Block),
395:   PetscDesignatedInitializer(oncomm, ISOnComm_Block),
396:   PetscDesignatedInitializer(setblocksize, ISSetBlockSize_Block),
397:   PetscDesignatedInitializer(contiguous, NULL),
398:   PetscDesignatedInitializer(locate, ISLocate_Block),
399:   /* we can have specialized local routines for determining properties,
400:    * but unless the block size is the same on each process (which is not guaranteed at
401:    * the moment), then trying to do something specialized for global properties is too
402:    * complicated */
403:   PetscDesignatedInitializer(sortedlocal, ISSortedLocal_Block),
404:   PetscDesignatedInitializer(sortedglobal, NULL),
405:   PetscDesignatedInitializer(uniquelocal, ISUniqueLocal_Block),
406:   PetscDesignatedInitializer(uniqueglobal, NULL),
407:   PetscDesignatedInitializer(permlocal, ISPermutationLocal_Block),
408:   PetscDesignatedInitializer(permglobal, NULL),
409:   PetscDesignatedInitializer(intervallocal, ISIntervalLocal_Block),
410:   PetscDesignatedInitializer(intervalglobal, NULL)
411: };
412: // clang-format on

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

417:   Collective

419:   Input Parameters:
420: + is   - the index set
421: . bs   - number of elements in each block
422: . n    - the length of the index set (the number of blocks)
423: . 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
424: - mode - see `PetscCopyMode`, only `PETSC_COPY_VALUES` and `PETSC_OWN_POINTER` are supported

426:   Level: beginner

428:   Notes:
429:   When the communicator is not `MPI_COMM_SELF`, the operations on the
430:   index sets, IS, are NOT conceptually the same as `MPI_Group` operations.
431:   The index sets are then distributed sets of indices and thus certain operations
432:   on them are collective.

434:   The convenience routine `ISCreateBlock()` allows one to create the `IS` and provide the blocks in a single function call.

436:   Example:
437:   If you wish to index the values {0,1,4,5}, then use
438:   a block size of 2 and idx of {0,2}.

440: .seealso: [](sec_scatter), `IS`, `ISCreateStride()`, `ISCreateGeneral()`, `ISAllGather()`, `ISCreateBlock()`, `ISBLOCK`, `ISGeneralSetIndices()`
441: @*/
442: PetscErrorCode ISBlockSetIndices(IS is, PetscInt bs, PetscInt n, const PetscInt idx[], PetscCopyMode mode)
443: {
444:   PetscFunctionBegin;
445:   PetscCall(ISClearInfoCache(is, PETSC_FALSE));
446:   PetscUseMethod(is, "ISBlockSetIndices_C", (IS, PetscInt, PetscInt, const PetscInt[], PetscCopyMode), (is, bs, n, idx, mode));
447:   PetscFunctionReturn(PETSC_SUCCESS);
448: }

450: static PetscErrorCode ISBlockSetIndices_Block(IS is, PetscInt bs, PetscInt n, const PetscInt idx[], PetscCopyMode mode)
451: {
452:   PetscInt    i, min, max;
453:   IS_Block   *sub = (IS_Block *)is->data;
454:   PetscLayout map;

456:   PetscFunctionBegin;
457:   PetscCheck(bs >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "block size < 1");
458:   PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length < 0");
459:   if (n) PetscAssertPointer(idx, 4);

461:   PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), n * bs, is->map->N, bs, &map));
462:   PetscCall(PetscLayoutDestroy(&is->map));
463:   is->map = map;

465:   if (sub->allocated) PetscCall(PetscFree(sub->idx));
466:   if (mode == PETSC_COPY_VALUES) {
467:     PetscCall(PetscMalloc1(n, &sub->idx));
468:     PetscCall(PetscArraycpy(sub->idx, idx, n));
469:     sub->allocated = PETSC_TRUE;
470:   } else if (mode == PETSC_OWN_POINTER) {
471:     sub->idx       = (PetscInt *)idx;
472:     sub->allocated = PETSC_TRUE;
473:   } else if (mode == PETSC_USE_POINTER) {
474:     sub->idx       = (PetscInt *)idx;
475:     sub->allocated = PETSC_FALSE;
476:   }

478:   if (n) {
479:     min = max = idx[0];
480:     for (i = 1; i < n; i++) {
481:       if (idx[i] < min) min = idx[i];
482:       if (idx[i] > max) max = idx[i];
483:     }
484:     is->min = bs * min;
485:     is->max = bs * max + bs - 1;
486:   } else {
487:     is->min = PETSC_MAX_INT;
488:     is->max = PETSC_MIN_INT;
489:   }
490:   PetscFunctionReturn(PETSC_SUCCESS);
491: }

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

497:   Collective

499:   Input Parameters:
500: + comm - the MPI communicator
501: . bs   - number of elements in each block
502: . n    - the length of the index set (the number of blocks)
503: . 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
504: - mode - see `PetscCopyMode`, only `PETSC_COPY_VALUES` and `PETSC_OWN_POINTER` are supported in this routine

506:   Output Parameter:
507: . is - the new index set

509:   Level: beginner

511:   Notes:
512:   When the communicator is not `MPI_COMM_SELF`, the operations on the
513:   index sets, `IS`, are NOT conceptually the same as `MPI_Group` operations.
514:   The index sets are then distributed sets of indices and thus certain operations
515:   on them are collective.

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

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

523: .seealso: [](sec_scatter), `IS`, `ISCreateStride()`, `ISCreateGeneral()`, `ISAllGather()`, `ISBlockSetIndices()`, `ISBLOCK`, `ISGENERAL`
524: @*/
525: PetscErrorCode ISCreateBlock(MPI_Comm comm, PetscInt bs, PetscInt n, const PetscInt idx[], PetscCopyMode mode, IS *is)
526: {
527:   PetscFunctionBegin;
528:   PetscAssertPointer(is, 6);
529:   PetscCheck(bs >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "block size < 1");
530:   PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "length < 0");
531:   if (n) PetscAssertPointer(idx, 4);

533:   PetscCall(ISCreate(comm, is));
534:   PetscCall(ISSetType(*is, ISBLOCK));
535:   PetscCall(ISBlockSetIndices(*is, bs, n, idx, mode));
536:   PetscFunctionReturn(PETSC_SUCCESS);
537: }

539: static PetscErrorCode ISBlockGetIndices_Block(IS is, const PetscInt *idx[])
540: {
541:   IS_Block *sub = (IS_Block *)is->data;

543:   PetscFunctionBegin;
544:   *idx = sub->idx;
545:   PetscFunctionReturn(PETSC_SUCCESS);
546: }

548: static PetscErrorCode ISBlockRestoreIndices_Block(IS is, const PetscInt *idx[])
549: {
550:   PetscFunctionBegin;
551:   PetscFunctionReturn(PETSC_SUCCESS);
552: }

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

557:   Not Collective

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

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

565:   Level: intermediate

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

570: .seealso: [](sec_scatter), `IS`, `ISBLOCK`, `ISGetIndices()`, `ISBlockRestoreIndices()`, `ISBlockSetIndices()`, `ISCreateBlock()`
571: @*/
572: PetscErrorCode ISBlockGetIndices(IS is, const PetscInt *idx[])
573: {
574:   PetscFunctionBegin;
575:   PetscUseMethod(is, "ISBlockGetIndices_C", (IS, const PetscInt *[]), (is, idx));
576:   PetscFunctionReturn(PETSC_SUCCESS);
577: }

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

582:   Not Collective

584:   Input Parameter:
585: . is - the index set

587:   Output Parameter:
588: . idx - the integer indices

590:   Level: intermediate

592: .seealso: [](sec_scatter), `IS`, `ISBLOCK`, `ISRestoreIndices()`, `ISBlockGetIndices()`
593: @*/
594: PetscErrorCode ISBlockRestoreIndices(IS is, const PetscInt *idx[])
595: {
596:   PetscFunctionBegin;
597:   PetscUseMethod(is, "ISBlockRestoreIndices_C", (IS, const PetscInt *[]), (is, idx));
598:   PetscFunctionReturn(PETSC_SUCCESS);
599: }

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

604:   Not Collective

606:   Input Parameter:
607: . is - the index set

609:   Output Parameter:
610: . size - the local number of blocks

612:   Level: intermediate

614: .seealso: [](sec_scatter), `IS`, `ISGetBlockSize()`, `ISBlockGetSize()`, `ISGetSize()`, `ISCreateBlock()`, `ISBLOCK`
615: @*/
616: PetscErrorCode ISBlockGetLocalSize(IS is, PetscInt *size)
617: {
618:   PetscFunctionBegin;
619:   PetscUseMethod(is, "ISBlockGetLocalSize_C", (IS, PetscInt *), (is, size));
620:   PetscFunctionReturn(PETSC_SUCCESS);
621: }

623: static PetscErrorCode ISBlockGetLocalSize_Block(IS is, PetscInt *size)
624: {
625:   PetscInt bs, n;

627:   PetscFunctionBegin;
628:   PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
629:   PetscCall(PetscLayoutGetLocalSize(is->map, &n));
630:   *size = n / bs;
631:   PetscFunctionReturn(PETSC_SUCCESS);
632: }

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

637:   Not Collective

639:   Input Parameter:
640: . is - the index set

642:   Output Parameter:
643: . size - the global number of blocks

645:   Level: intermediate

647: .seealso: [](sec_scatter), `IS`, `ISGetBlockSize()`, `ISBlockGetLocalSize()`, `ISGetSize()`, `ISCreateBlock()`, `ISBLOCK`
648: @*/
649: PetscErrorCode ISBlockGetSize(IS is, PetscInt *size)
650: {
651:   PetscFunctionBegin;
652:   PetscUseMethod(is, "ISBlockGetSize_C", (IS, PetscInt *), (is, size));
653:   PetscFunctionReturn(PETSC_SUCCESS);
654: }

656: static PetscErrorCode ISBlockGetSize_Block(IS is, PetscInt *size)
657: {
658:   PetscInt bs, N;

660:   PetscFunctionBegin;
661:   PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
662:   PetscCall(PetscLayoutGetSize(is->map, &N));
663:   *size = N / bs;
664:   PetscFunctionReturn(PETSC_SUCCESS);
665: }

667: PETSC_INTERN PetscErrorCode ISCreate_Block(IS is)
668: {
669:   IS_Block *sub;

671:   PetscFunctionBegin;
672:   PetscCall(PetscNew(&sub));
673:   is->data   = (void *)sub;
674:   is->ops[0] = myops;
675:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockSetIndices_C", ISBlockSetIndices_Block));
676:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockGetIndices_C", ISBlockGetIndices_Block));
677:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockRestoreIndices_C", ISBlockRestoreIndices_Block));
678:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockGetSize_C", ISBlockGetSize_Block));
679:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISBlockGetLocalSize_C", ISBlockGetLocalSize_Block));
680:   PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISShift_C", ISShift_Block));
681:   PetscFunctionReturn(PETSC_SUCCESS);
682: }