Actual source code: iscoloring.c


  2: #include <petsc/private/isimpl.h>
  3: #include <petscviewer.h>
  4: #include <petscsf.h>

  6: const char *const ISColoringTypes[] = {"global", "ghosted", "ISColoringType", "IS_COLORING_", NULL};

  8: PetscErrorCode ISColoringReference(ISColoring coloring)
  9: {
 10:   coloring->refct++;
 11:   return 0;
 12: }

 14: /*@C
 15:    ISColoringSetType - indicates if the coloring is for the local representation (including ghost points) or the global representation of a `Mat`

 17:    Collective

 19:    Input Parameters:
 20: +    coloring - the coloring object
 21: -    type - either `IS_COLORING_LOCAL` or `IS_COLORING_GLOBAL`

 23:    Notes:
 24:    `IS_COLORING_LOCAL` can lead to faster computations since parallel ghost point updates are not needed for each color

 26:    With `IS_COLORING_LOCAL` the coloring is in the numbering of the local vector, for `IS_COLORING_GLOBAL` it is in the numbering of the global vector

 28:    Level: intermediate

 30: .seealso: `MatFDColoringCreate()`, `ISColoring`, `ISColoringType`, `ISColoringCreate()`, `IS_COLORING_LOCAL`, `IS_COLORING_GLOBAL`, `ISColoringGetType()`
 31: @*/
 32: PetscErrorCode ISColoringSetType(ISColoring coloring, ISColoringType type)
 33: {
 34:   coloring->ctype = type;
 35:   return 0;
 36: }

 38: /*@C

 40:     ISColoringGetType - gets if the coloring is for the local representation (including ghost points) or the global representation

 42:    Collective

 44:    Input Parameter:
 45: .   coloring - the coloring object

 47:    Output Parameter:
 48: .    type - either `IS_COLORING_LOCAL` or `IS_COLORING_GLOBAL`

 50:    Level: intermediate

 52: .seealso: `MatFDColoringCreate()`, `ISColoring`, `ISColoringType`, `ISColoringCreate()`, `IS_COLORING_LOCAL`, `IS_COLORING_GLOBAL`, `ISColoringSetType()`
 53: @*/
 54: PetscErrorCode ISColoringGetType(ISColoring coloring, ISColoringType *type)
 55: {
 56:   *type = coloring->ctype;
 57:   return 0;
 58: }

 60: /*@
 61:    ISColoringDestroy - Destroys an `ISColoring` coloring context.

 63:    Collective

 65:    Input Parameter:
 66: .  iscoloring - the coloring context

 68:    Level: advanced

 70: .seealso: `ISColoring`, `ISColoringView()`, `MatColoring`
 71: @*/
 72: PetscErrorCode ISColoringDestroy(ISColoring *iscoloring)
 73: {
 74:   PetscInt i;

 76:   if (!*iscoloring) return 0;
 78:   if (--(*iscoloring)->refct > 0) {
 79:     *iscoloring = NULL;
 80:     return 0;
 81:   }

 83:   if ((*iscoloring)->is) {
 84:     for (i = 0; i < (*iscoloring)->n; i++) ISDestroy(&(*iscoloring)->is[i]);
 85:     PetscFree((*iscoloring)->is);
 86:   }
 87:   if ((*iscoloring)->allocated) PetscFree((*iscoloring)->colors);
 88:   PetscCommDestroy(&(*iscoloring)->comm);
 89:   PetscFree((*iscoloring));
 90:   return 0;
 91: }

 93: /*@C
 94:   ISColoringViewFromOptions - Processes command line options to determine if/how an `ISColoring` object is to be viewed.

 96:   Collective

 98:   Input Parameters:
 99: + obj   - the `ISColoring` object
100: . prefix - prefix to use for viewing, or NULL to use prefix of 'mat'
101: - optionname - option to activate viewing

103:   Level: intermediate

105:   Developer Note:
106:   This cannot use `PetscObjectViewFromOptions()` because `ISColoring` is not a `PetscObject`

108: .seealso: `ISColoring`, `ISColoringView()`
109: @*/
110: PetscErrorCode ISColoringViewFromOptions(ISColoring obj, PetscObject bobj, const char optionname[])
111: {
112:   PetscViewer       viewer;
113:   PetscBool         flg;
114:   PetscViewerFormat format;
115:   char             *prefix;

117:   prefix = bobj ? bobj->prefix : NULL;
118:   PetscOptionsGetViewer(obj->comm, NULL, prefix, optionname, &viewer, &format, &flg);
119:   if (flg) {
120:     PetscViewerPushFormat(viewer, format);
121:     ISColoringView(obj, viewer);
122:     PetscViewerPopFormat(viewer);
123:     PetscViewerDestroy(&viewer);
124:   }
125:   return 0;
126: }

128: /*@C
129:    ISColoringView - Views an `ISColoring` coloring context.

131:    Collective

133:    Input Parameters:
134: +  iscoloring - the coloring context
135: -  viewer - the viewer

137:    Level: advanced

139: .seealso: `ISColoring()`, `ISColoringViewFromOptions()`, `ISColoringDestroy()`, `ISColoringGetIS()`, `MatColoring`
140: @*/
141: PetscErrorCode ISColoringView(ISColoring iscoloring, PetscViewer viewer)
142: {
143:   PetscInt  i;
144:   PetscBool iascii;
145:   IS       *is;

148:   if (!viewer) PetscViewerASCIIGetStdout(iscoloring->comm, &viewer);

151:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
152:   if (iascii) {
153:     MPI_Comm    comm;
154:     PetscMPIInt size, rank;

156:     PetscObjectGetComm((PetscObject)viewer, &comm);
157:     MPI_Comm_size(comm, &size);
158:     MPI_Comm_rank(comm, &rank);
159:     PetscViewerASCIIPrintf(viewer, "ISColoring Object: %d MPI processes\n", size);
160:     PetscViewerASCIIPrintf(viewer, "ISColoringType: %s\n", ISColoringTypes[iscoloring->ctype]);
161:     PetscViewerASCIIPushSynchronized(viewer);
162:     PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Number of colors %" PetscInt_FMT "\n", rank, iscoloring->n);
163:     PetscViewerFlush(viewer);
164:     PetscViewerASCIIPopSynchronized(viewer);
165:   }

167:   ISColoringGetIS(iscoloring, PETSC_USE_POINTER, PETSC_IGNORE, &is);
168:   for (i = 0; i < iscoloring->n; i++) ISView(iscoloring->is[i], viewer);
169:   ISColoringRestoreIS(iscoloring, PETSC_USE_POINTER, &is);
170:   return 0;
171: }

173: /*@C
174:    ISColoringGetColors - Returns an array with the color for each local node

176:    Not Collective

178:    Input Parameter:
179: .  iscoloring - the coloring context

181:    Output Parameters:
182: +  n - number of nodes
183: .  nc - number of colors
184: -  colors - color for each node

186:    Level: advanced

188:    Notes:
189:    Do not free the colors[] array.

191:    The colors[] array will only be valid for the lifetime of the `ISColoring`

193: .seealso: `ISColoring`, `ISColoringValue`, `ISColoringRestoreIS()`, `ISColoringView()`, `ISColoringGetIS()`
194: @*/
195: PetscErrorCode ISColoringGetColors(ISColoring iscoloring, PetscInt *n, PetscInt *nc, const ISColoringValue **colors)
196: {

199:   if (n) *n = iscoloring->N;
200:   if (nc) *nc = iscoloring->n;
201:   if (colors) *colors = iscoloring->colors;
202:   return 0;
203: }

205: /*@C
206:    ISColoringGetIS - Extracts index sets from the coloring context. Each is contains the nodes of one color

208:    Collective

210:    Input Parameters:
211: +  iscoloring - the coloring context
212: -  mode - if this value is `PETSC_OWN_POINTER` then the caller owns the pointer and must free the array of `IS` and each `IS` in the array

214:    Output Parameters:
215: +  nn - number of index sets in the coloring context
216: -  is - array of index sets

218:    Level: advanced

220:    Note:
221:    If mode is `PETSC_USE_POINTER` then `ISColoringRestoreIS()` must be called when the `IS` are no longer needed

223: .seealso: `ISColoring`, `IS`, `ISColoringRestoreIS()`, `ISColoringView()`, `ISColoringGetColoring()`, `ISColoringGetColors()`
224: @*/
225: PetscErrorCode ISColoringGetIS(ISColoring iscoloring, PetscCopyMode mode, PetscInt *nn, IS *isis[])
226: {

229:   if (nn) *nn = iscoloring->n;
230:   if (isis) {
231:     if (!iscoloring->is) {
232:       PetscInt        *mcolors, **ii, nc = iscoloring->n, i, base, n = iscoloring->N;
233:       ISColoringValue *colors = iscoloring->colors;
234:       IS              *is;

236:       if (PetscDefined(USE_DEBUG)) {
238:       }

240:       /* generate the lists of nodes for each color */
241:       PetscCalloc1(nc, &mcolors);
242:       for (i = 0; i < n; i++) mcolors[colors[i]]++;

244:       PetscMalloc1(nc, &ii);
245:       PetscMalloc1(n, &ii[0]);
246:       for (i = 1; i < nc; i++) ii[i] = ii[i - 1] + mcolors[i - 1];
247:       PetscArrayzero(mcolors, nc);

249:       if (iscoloring->ctype == IS_COLORING_GLOBAL) {
250:         MPI_Scan(&iscoloring->N, &base, 1, MPIU_INT, MPI_SUM, iscoloring->comm);
251:         base -= iscoloring->N;
252:         for (i = 0; i < n; i++) ii[colors[i]][mcolors[colors[i]]++] = i + base; /* global idx */
253:       } else if (iscoloring->ctype == IS_COLORING_LOCAL) {
254:         for (i = 0; i < n; i++) ii[colors[i]][mcolors[colors[i]]++] = i; /* local idx */
255:       } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not provided for this ISColoringType type");

257:       PetscMalloc1(nc, &is);
258:       for (i = 0; i < nc; i++) ISCreateGeneral(iscoloring->comm, mcolors[i], ii[i], PETSC_COPY_VALUES, is + i);

260:       if (mode != PETSC_OWN_POINTER) iscoloring->is = is;
261:       *isis = is;
262:       PetscFree(ii[0]);
263:       PetscFree(ii);
264:       PetscFree(mcolors);
265:     } else {
266:       *isis = iscoloring->is;
267:     }
268:   }
269:   return 0;
270: }

272: /*@C
273:    ISColoringRestoreIS - Restores the index sets extracted from the coloring context with `ISColoringGetIS()` using `PETSC_USE_POINTER`

275:    Collective

277:    Input Parameters:
278: +  iscoloring - the coloring context
279: .  mode - who retains ownership of the is
280: -  is - array of index sets

282:    Level: advanced

284: .seealso: `ISColoring()`, `IS`, `ISColoringGetIS()`, `ISColoringView()`, `PetscCopyMode`
285: @*/
286: PetscErrorCode ISColoringRestoreIS(ISColoring iscoloring, PetscCopyMode mode, IS *is[])
287: {

290:   /* currently nothing is done here */
291:   return 0;
292: }

294: /*@
295:     ISColoringCreate - Generates an `ISColoring` context from lists (provided by each MPI rank) of colors for each node.

297:     Collective

299:     Input Parameters:
300: +   comm - communicator for the processors creating the coloring
301: .   ncolors - max color value
302: .   n - number of nodes on this processor
303: .   colors - array containing the colors for this MPI rank, color numbers begin at 0, for each local node
304: -   mode - see `PetscCopyMode` for meaning of this flag.

306:     Output Parameter:
307: .   iscoloring - the resulting coloring data structure

309:     Options Database Key:
310: .   -is_coloring_view - Activates `ISColoringView()`

312:    Level: advanced

314:     Notes:
315:     By default sets coloring type to  `IS_COLORING_GLOBAL`

317: .seealso: `ISColoring`, `ISColoringValue`, `MatColoringCreate()`, `ISColoringView()`, `ISColoringDestroy()`, `ISColoringSetType()`
318: @*/
319: PetscErrorCode ISColoringCreate(MPI_Comm comm, PetscInt ncolors, PetscInt n, const ISColoringValue colors[], PetscCopyMode mode, ISColoring *iscoloring)
320: {
321:   PetscMPIInt size, rank, tag;
322:   PetscInt    base, top, i;
323:   PetscInt    nc, ncwork;
324:   MPI_Status  status;

326:   if (ncolors != PETSC_DECIDE && ncolors > IS_COLORING_MAX) {
328:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Max color value exceeds limit. Perhaps reconfigure PETSc with --with-is-color-value-type=short?\n Current max: %d user requested: %" PetscInt_FMT, PETSC_IS_COLORING_MAX, ncolors);
329:   }
330:   PetscNew(iscoloring);
331:   PetscCommDuplicate(comm, &(*iscoloring)->comm, &tag);
332:   comm = (*iscoloring)->comm;

334:   /* compute the number of the first node on my processor */
335:   MPI_Comm_size(comm, &size);

337:   /* should use MPI_Scan() */
338:   MPI_Comm_rank(comm, &rank);
339:   if (rank == 0) {
340:     base = 0;
341:     top  = n;
342:   } else {
343:     MPI_Recv(&base, 1, MPIU_INT, rank - 1, tag, comm, &status);
344:     top = base + n;
345:   }
346:   if (rank < size - 1) MPI_Send(&top, 1, MPIU_INT, rank + 1, tag, comm);

348:   /* compute the total number of colors */
349:   ncwork = 0;
350:   for (i = 0; i < n; i++) {
351:     if (ncwork < colors[i]) ncwork = colors[i];
352:   }
353:   ncwork++;
354:   MPIU_Allreduce(&ncwork, &nc, 1, MPIU_INT, MPI_MAX, comm);
356:   (*iscoloring)->n     = nc;
357:   (*iscoloring)->is    = NULL;
358:   (*iscoloring)->N     = n;
359:   (*iscoloring)->refct = 1;
360:   (*iscoloring)->ctype = IS_COLORING_GLOBAL;
361:   if (mode == PETSC_COPY_VALUES) {
362:     PetscMalloc1(n, &(*iscoloring)->colors);
363:     PetscArraycpy((*iscoloring)->colors, colors, n);
364:     (*iscoloring)->allocated = PETSC_TRUE;
365:   } else if (mode == PETSC_OWN_POINTER) {
366:     (*iscoloring)->colors    = (ISColoringValue *)colors;
367:     (*iscoloring)->allocated = PETSC_TRUE;
368:   } else {
369:     (*iscoloring)->colors    = (ISColoringValue *)colors;
370:     (*iscoloring)->allocated = PETSC_FALSE;
371:   }
372:   ISColoringViewFromOptions(*iscoloring, NULL, "-is_coloring_view");
373:   PetscInfo(0, "Number of colors %" PetscInt_FMT "\n", nc);
374:   return 0;
375: }

377: /*@
378:     ISBuildTwoSided - Takes an `IS` that describes where each element will be mapped globally over all ranks.
379:     Generates an `IS` that contains new numbers from remote or local on the `IS`.

381:     Collective

383:     Input Parameters:
384: +   ito - an `IS` describes where each entry will be mapped. Negative target rank will be ignored
385: -   toindx - an `IS` describes what indices should send. NULL means sending natural numbering

387:     Output Parameter:
388: .   rows - contains new numbers from remote or local

390:    Level: advanced

392:    Developer Note:
393:    This manual page is incomprehensible and needs to be fixed

395: .seealso: [](sec_scatter), `IS`, `MatPartitioningCreate()`, `ISPartitioningToNumbering()`, `ISPartitioningCount()`
396: @*/
397: PetscErrorCode ISBuildTwoSided(IS ito, IS toindx, IS *rows)
398: {
399:   const PetscInt *ito_indices, *toindx_indices;
400:   PetscInt       *send_indices, rstart, *recv_indices, nrecvs, nsends;
401:   PetscInt       *tosizes, *fromsizes, i, j, *tosizes_tmp, *tooffsets_tmp, ito_ln;
402:   PetscMPIInt    *toranks, *fromranks, size, target_rank, *fromperm_newtoold, nto, nfrom;
403:   PetscLayout     isrmap;
404:   MPI_Comm        comm;
405:   PetscSF         sf;
406:   PetscSFNode    *iremote;

408:   PetscObjectGetComm((PetscObject)ito, &comm);
409:   MPI_Comm_size(comm, &size);
410:   ISGetLocalSize(ito, &ito_ln);
411:   ISGetLayout(ito, &isrmap);
412:   PetscLayoutGetRange(isrmap, &rstart, NULL);
413:   ISGetIndices(ito, &ito_indices);
414:   PetscCalloc2(size, &tosizes_tmp, size + 1, &tooffsets_tmp);
415:   for (i = 0; i < ito_ln; i++) {
416:     if (ito_indices[i] < 0) continue;
418:     tosizes_tmp[ito_indices[i]]++;
419:   }
420:   nto = 0;
421:   for (i = 0; i < size; i++) {
422:     tooffsets_tmp[i + 1] = tooffsets_tmp[i] + tosizes_tmp[i];
423:     if (tosizes_tmp[i] > 0) nto++;
424:   }
425:   PetscCalloc2(nto, &toranks, 2 * nto, &tosizes);
426:   nto = 0;
427:   for (i = 0; i < size; i++) {
428:     if (tosizes_tmp[i] > 0) {
429:       toranks[nto]         = i;
430:       tosizes[2 * nto]     = tosizes_tmp[i];   /* size */
431:       tosizes[2 * nto + 1] = tooffsets_tmp[i]; /* offset */
432:       nto++;
433:     }
434:   }
435:   nsends = tooffsets_tmp[size];
436:   PetscCalloc1(nsends, &send_indices);
437:   if (toindx) ISGetIndices(toindx, &toindx_indices);
438:   for (i = 0; i < ito_ln; i++) {
439:     if (ito_indices[i] < 0) continue;
440:     target_rank                              = ito_indices[i];
441:     send_indices[tooffsets_tmp[target_rank]] = toindx ? toindx_indices[i] : (i + rstart);
442:     tooffsets_tmp[target_rank]++;
443:   }
444:   if (toindx) ISRestoreIndices(toindx, &toindx_indices);
445:   ISRestoreIndices(ito, &ito_indices);
446:   PetscFree2(tosizes_tmp, tooffsets_tmp);
447:   PetscCommBuildTwoSided(comm, 2, MPIU_INT, nto, toranks, tosizes, &nfrom, &fromranks, &fromsizes);
448:   PetscFree2(toranks, tosizes);
449:   PetscMalloc1(nfrom, &fromperm_newtoold);
450:   for (i = 0; i < nfrom; i++) fromperm_newtoold[i] = i;
451:   PetscSortMPIIntWithArray(nfrom, fromranks, fromperm_newtoold);
452:   nrecvs = 0;
453:   for (i = 0; i < nfrom; i++) nrecvs += fromsizes[i * 2];
454:   PetscCalloc1(nrecvs, &recv_indices);
455:   PetscMalloc1(nrecvs, &iremote);
456:   nrecvs = 0;
457:   for (i = 0; i < nfrom; i++) {
458:     for (j = 0; j < fromsizes[2 * fromperm_newtoold[i]]; j++) {
459:       iremote[nrecvs].rank    = fromranks[i];
460:       iremote[nrecvs++].index = fromsizes[2 * fromperm_newtoold[i] + 1] + j;
461:     }
462:   }
463:   PetscSFCreate(comm, &sf);
464:   PetscSFSetGraph(sf, nsends, nrecvs, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
465:   PetscSFSetType(sf, PETSCSFBASIC);
466:   /* how to put a prefix ? */
467:   PetscSFSetFromOptions(sf);
468:   PetscSFBcastBegin(sf, MPIU_INT, send_indices, recv_indices, MPI_REPLACE);
469:   PetscSFBcastEnd(sf, MPIU_INT, send_indices, recv_indices, MPI_REPLACE);
470:   PetscSFDestroy(&sf);
471:   PetscFree(fromranks);
472:   PetscFree(fromsizes);
473:   PetscFree(fromperm_newtoold);
474:   PetscFree(send_indices);
475:   if (rows) {
476:     PetscSortInt(nrecvs, recv_indices);
477:     ISCreateGeneral(comm, nrecvs, recv_indices, PETSC_OWN_POINTER, rows);
478:   } else {
479:     PetscFree(recv_indices);
480:   }
481:   return 0;
482: }

484: /*@
485:     ISPartitioningToNumbering - Takes an `IS' that represents a partitioning (the MPI rank that each local entry belongs to) and on each MPI rank
486:     generates an `IS` that contains a new global node number in the new ordering for each entry

488:     Collective

490:     Input Parameters:
491: .   partitioning - a partitioning as generated by `MatPartitioningApply()` or `MatPartitioningApplyND()`

493:     Output Parameter:
494: .   is - on each processor the index set that defines the global numbers
495:          (in the new numbering) for all the nodes currently (before the partitioning)
496:          on that processor

498:    Level: advanced

500:    Note:
501:    The resulting `IS` tells where each local entry is mapped to in a new global ordering

503: .seealso: [](sec_scatter), `IS`, `MatPartitioningCreate()`, `AOCreateBasic()`, `ISPartitioningCount()`
504: @*/
505: PetscErrorCode ISPartitioningToNumbering(IS part, IS *is)
506: {
507:   MPI_Comm        comm;
508:   IS              ndorder;
509:   PetscInt        i, np, npt, n, *starts = NULL, *sums = NULL, *lsizes = NULL, *newi = NULL;
510:   const PetscInt *indices = NULL;

514:   /* see if the partitioning comes from nested dissection */
515:   PetscObjectQuery((PetscObject)part, "_petsc_matpartitioning_ndorder", (PetscObject *)&ndorder);
516:   if (ndorder) {
517:     PetscObjectReference((PetscObject)ndorder);
518:     *is = ndorder;
519:     return 0;
520:   }

522:   PetscObjectGetComm((PetscObject)part, &comm);
523:   /* count the number of partitions, i.e., virtual processors */
524:   ISGetLocalSize(part, &n);
525:   ISGetIndices(part, &indices);
526:   np = 0;
527:   for (i = 0; i < n; i++) np = PetscMax(np, indices[i]);
528:   MPIU_Allreduce(&np, &npt, 1, MPIU_INT, MPI_MAX, comm);
529:   np = npt + 1; /* so that it looks like a MPI_Comm_size output */

531:   /*
532:         lsizes - number of elements of each partition on this particular processor
533:         sums - total number of "previous" nodes for any particular partition
534:         starts - global number of first element in each partition on this processor
535:   */
536:   PetscMalloc3(np, &lsizes, np, &starts, np, &sums);
537:   PetscArrayzero(lsizes, np);
538:   for (i = 0; i < n; i++) lsizes[indices[i]]++;
539:   MPIU_Allreduce(lsizes, sums, np, MPIU_INT, MPI_SUM, comm);
540:   MPI_Scan(lsizes, starts, np, MPIU_INT, MPI_SUM, comm);
541:   for (i = 0; i < np; i++) starts[i] -= lsizes[i];
542:   for (i = 1; i < np; i++) {
543:     sums[i] += sums[i - 1];
544:     starts[i] += sums[i - 1];
545:   }

547:   /*
548:       For each local index give it the new global number
549:   */
550:   PetscMalloc1(n, &newi);
551:   for (i = 0; i < n; i++) newi[i] = starts[indices[i]]++;
552:   PetscFree3(lsizes, starts, sums);

554:   ISRestoreIndices(part, &indices);
555:   ISCreateGeneral(comm, n, newi, PETSC_OWN_POINTER, is);
556:   ISSetPermutation(*is);
557:   return 0;
558: }

560: /*@
561:     ISPartitioningCount - Takes a `IS` that represents a partitioning (the MPI rank that each local entry belongs to) and determines the number of
562:     resulting elements on each (partition) rank

564:     Collective

566:     Input Parameters:
567: +   partitioning - a partitioning as generated by `MatPartitioningApply()` or `MatPartitioningApplyND()`
568: -   len - length of the array count, this is the total number of partitions

570:     Output Parameter:
571: .   count - array of length size, to contain the number of elements assigned
572:         to each partition, where size is the number of partitions generated
573:          (see notes below).

575:    Level: advanced

577:     Notes:
578:     By default the number of partitions generated (and thus the length
579:     of count) is the size of the communicator associated with `IS`,
580:     but it can be set by `MatPartitioningSetNParts()`.

582:     The resulting array of lengths can for instance serve as input of `PCBJacobiSetTotalBlocks()`.

584:     If the partitioning has been obtained by `MatPartitioningApplyND()`, the returned count does not include the separators.

586: .seealso: [](sec_scatter), `IS`, `MatPartitioningCreate()`, `AOCreateBasic()`, `ISPartitioningToNumbering()`,
587:           `MatPartitioningSetNParts()`, `MatPartitioningApply()`, `MatPartitioningApplyND()`
588: @*/
589: PetscErrorCode ISPartitioningCount(IS part, PetscInt len, PetscInt count[])
590: {
591:   MPI_Comm        comm;
592:   PetscInt        i, n, *lsizes;
593:   const PetscInt *indices;
594:   PetscMPIInt     npp;

596:   PetscObjectGetComm((PetscObject)part, &comm);
597:   if (len == PETSC_DEFAULT) {
598:     PetscMPIInt size;
599:     MPI_Comm_size(comm, &size);
600:     len = (PetscInt)size;
601:   }

603:   /* count the number of partitions */
604:   ISGetLocalSize(part, &n);
605:   ISGetIndices(part, &indices);
606:   if (PetscDefined(USE_DEBUG)) {
607:     PetscInt np = 0, npt;
608:     for (i = 0; i < n; i++) np = PetscMax(np, indices[i]);
609:     MPIU_Allreduce(&np, &npt, 1, MPIU_INT, MPI_MAX, comm);
610:     np = npt + 1; /* so that it looks like a MPI_Comm_size output */
612:   }

614:   /*
615:         lsizes - number of elements of each partition on this particular processor
616:         sums - total number of "previous" nodes for any particular partition
617:         starts - global number of first element in each partition on this processor
618:   */
619:   PetscCalloc1(len, &lsizes);
620:   for (i = 0; i < n; i++) {
621:     if (indices[i] > -1) lsizes[indices[i]]++;
622:   }
623:   ISRestoreIndices(part, &indices);
624:   PetscMPIIntCast(len, &npp);
625:   MPIU_Allreduce(lsizes, count, npp, MPIU_INT, MPI_SUM, comm);
626:   PetscFree(lsizes);
627:   return 0;
628: }

630: /*@
631:     ISAllGather - Given an index set `IS` on each processor, generates a large
632:     index set (same on each processor) by concatenating together each
633:     processors index set.

635:     Collective

637:     Input Parameter:
638: .   is - the distributed index set

640:     Output Parameter:
641: .   isout - the concatenated index set (same on all processors)

643:     Level: intermediate

645:     Notes:
646:     `ISAllGather()` is clearly not scalable for large index sets.

648:     The `IS` created on each processor must be created with a common
649:     communicator (e.g., `PETSC_COMM_WORLD`). If the index sets were created
650:     with `PETSC_COMM_SELF`, this routine will not work as expected, since
651:     each process will generate its own new `IS` that consists only of
652:     itself.

654:     The communicator for this new `IS` is `PETSC_COMM_SELF`

656: .seealso: [](sec_scatter), `IS`, `ISCreateGeneral()`, `ISCreateStride()`, `ISCreateBlock()`
657: @*/
658: PetscErrorCode ISAllGather(IS is, IS *isout)
659: {
660:   PetscInt       *indices, n, i, N, step, first;
661:   const PetscInt *lindices;
662:   MPI_Comm        comm;
663:   PetscMPIInt     size, *sizes = NULL, *offsets = NULL, nn;
664:   PetscBool       stride;


669:   PetscObjectGetComm((PetscObject)is, &comm);
670:   MPI_Comm_size(comm, &size);
671:   ISGetLocalSize(is, &n);
672:   PetscObjectTypeCompare((PetscObject)is, ISSTRIDE, &stride);
673:   if (size == 1 && stride) { /* should handle parallel ISStride also */
674:     ISStrideGetInfo(is, &first, &step);
675:     ISCreateStride(PETSC_COMM_SELF, n, first, step, isout);
676:   } else {
677:     PetscMalloc2(size, &sizes, size, &offsets);

679:     PetscMPIIntCast(n, &nn);
680:     MPI_Allgather(&nn, 1, MPI_INT, sizes, 1, MPI_INT, comm);
681:     offsets[0] = 0;
682:     for (i = 1; i < size; i++) {
683:       PetscInt s = offsets[i - 1] + sizes[i - 1];
684:       PetscMPIIntCast(s, &offsets[i]);
685:     }
686:     N = offsets[size - 1] + sizes[size - 1];

688:     PetscMalloc1(N, &indices);
689:     ISGetIndices(is, &lindices);
690:     MPI_Allgatherv((void *)lindices, nn, MPIU_INT, indices, sizes, offsets, MPIU_INT, comm);
691:     ISRestoreIndices(is, &lindices);
692:     PetscFree2(sizes, offsets);

694:     ISCreateGeneral(PETSC_COMM_SELF, N, indices, PETSC_OWN_POINTER, isout);
695:   }
696:   return 0;
697: }

699: /*@C
700:     ISAllGatherColors - Given a a set of colors on each processor, generates a large
701:     set (same on each processor) by concatenating together each processors colors

703:     Collective

705:     Input Parameters:
706: +   comm - communicator to share the indices
707: .   n - local size of set
708: -   lindices - local colors

710:     Output Parameters:
711: +   outN - total number of indices
712: -   outindices - all of the colors

714:     Level: intermediate

716:     Note:
717:     `ISAllGatherColors()` is clearly not scalable for large index sets.

719: .seealso: `ISCOloringValue`, `ISColoring()`, `ISCreateGeneral()`, `ISCreateStride()`, `ISCreateBlock()`, `ISAllGather()`
720: @*/
721: PetscErrorCode ISAllGatherColors(MPI_Comm comm, PetscInt n, ISColoringValue *lindices, PetscInt *outN, ISColoringValue *outindices[])
722: {
723:   ISColoringValue *indices;
724:   PetscInt         i, N;
725:   PetscMPIInt      size, *offsets = NULL, *sizes = NULL, nn = n;

727:   MPI_Comm_size(comm, &size);
728:   PetscMalloc2(size, &sizes, size, &offsets);

730:   MPI_Allgather(&nn, 1, MPI_INT, sizes, 1, MPI_INT, comm);
731:   offsets[0] = 0;
732:   for (i = 1; i < size; i++) offsets[i] = offsets[i - 1] + sizes[i - 1];
733:   N = offsets[size - 1] + sizes[size - 1];
734:   PetscFree2(sizes, offsets);

736:   PetscMalloc1(N + 1, &indices);
737:   MPI_Allgatherv(lindices, (PetscMPIInt)n, MPIU_COLORING_VALUE, indices, sizes, offsets, MPIU_COLORING_VALUE, comm);

739:   *outindices = indices;
740:   if (outN) *outN = N;
741:   return 0;
742: }

744: /*@
745:     ISComplement - Given an index set `IS` generates the complement index set. That is all
746:        all indices that are NOT in the given set.

748:     Collective

750:     Input Parameters:
751: +   is - the index set
752: .   nmin - the first index desired in the local part of the complement
753: -   nmax - the largest index desired in the local part of the complement (note that all indices in is must be greater or equal to nmin and less than nmax)

755:     Output Parameter:
756: .   isout - the complement

758:     Level: intermediate

760:     Notes:
761:     The communicator for this new `IS` is the same as for the input IS

763:     For a parallel `IS`, this will generate the local part of the complement on each process

765:     To generate the entire complement (on each process) of a parallel `IS`, first call `ISAllGather()` and then
766:     call this routine.

768: .seealso: [](sec_scatter), `IS`, `ISCreateGeneral()`, `ISCreateStride()`, `ISCreateBlock()`, `ISAllGather()`
769: @*/
770: PetscErrorCode ISComplement(IS is, PetscInt nmin, PetscInt nmax, IS *isout)
771: {
772:   const PetscInt *indices;
773:   PetscInt        n, i, j, unique, cnt, *nindices;
774:   PetscBool       sorted;

780:   ISSorted(is, &sorted);

783:   ISGetLocalSize(is, &n);
784:   ISGetIndices(is, &indices);
785:   if (PetscDefined(USE_DEBUG)) {
786:     for (i = 0; i < n; i++) {
789:     }
790:   }
791:   /* Count number of unique entries */
792:   unique = (n > 0);
793:   for (i = 0; i < n - 1; i++) {
794:     if (indices[i + 1] != indices[i]) unique++;
795:   }
796:   PetscMalloc1(nmax - nmin - unique, &nindices);
797:   cnt = 0;
798:   for (i = nmin, j = 0; i < nmax; i++) {
799:     if (j < n && i == indices[j]) do {
800:         j++;
801:       } while (j < n && i == indices[j]);
802:     else nindices[cnt++] = i;
803:   }
805:   ISCreateGeneral(PetscObjectComm((PetscObject)is), cnt, nindices, PETSC_OWN_POINTER, isout);
806:   ISRestoreIndices(is, &indices);
807:   return 0;
808: }