Actual source code: iscoloring.c

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

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

  7: PetscErrorCode ISColoringReference(ISColoring coloring)
  8: {
  9:   PetscFunctionBegin;
 10:   coloring->refct++;
 11:   PetscFunctionReturn(PETSC_SUCCESS);
 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:   Level: intermediate

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

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

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

 39: /*@C

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

 43:   Collective

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

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

 51:   Level: intermediate

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

 62: /*@
 63:   ISColoringDestroy - Destroys an `ISColoring` coloring context.

 65:   Collective

 67:   Input Parameter:
 68: . iscoloring - the coloring context

 70:   Level: advanced

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

 78:   PetscFunctionBegin;
 79:   if (!*iscoloring) PetscFunctionReturn(PETSC_SUCCESS);
 80:   PetscAssertPointer(*iscoloring, 1);
 81:   if (--(*iscoloring)->refct > 0) {
 82:     *iscoloring = NULL;
 83:     PetscFunctionReturn(PETSC_SUCCESS);
 84:   }

 86:   if ((*iscoloring)->is) {
 87:     for (i = 0; i < (*iscoloring)->n; i++) PetscCall(ISDestroy(&(*iscoloring)->is[i]));
 88:     PetscCall(PetscFree((*iscoloring)->is));
 89:   }
 90:   if ((*iscoloring)->allocated) PetscCall(PetscFree((*iscoloring)->colors));
 91:   PetscCall(PetscCommDestroy(&(*iscoloring)->comm));
 92:   PetscCall(PetscFree(*iscoloring));
 93:   PetscFunctionReturn(PETSC_SUCCESS);
 94: }

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

 99:   Collective

101:   Input Parameters:
102: + obj        - the `ISColoring` object
103: . bobj       - prefix to use for viewing, or `NULL` to use prefix of `mat`
104: - optionname - option to activate viewing

106:   Level: intermediate

108:   Developer Notes:
109:   This cannot use `PetscObjectViewFromOptions()` because `ISColoring` is not a `PetscObject`

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

120:   PetscFunctionBegin;
121:   prefix = bobj ? bobj->prefix : NULL;
122:   PetscCall(PetscOptionsGetViewer(obj->comm, NULL, prefix, optionname, &viewer, &format, &flg));
123:   if (flg) {
124:     PetscCall(PetscViewerPushFormat(viewer, format));
125:     PetscCall(ISColoringView(obj, viewer));
126:     PetscCall(PetscViewerPopFormat(viewer));
127:     PetscCall(PetscOptionsRestoreViewer(&viewer));
128:   }
129:   PetscFunctionReturn(PETSC_SUCCESS);
130: }

132: /*@C
133:   ISColoringView - Views an `ISColoring` coloring context.

135:   Collective

137:   Input Parameters:
138: + iscoloring - the coloring context
139: - viewer     - the viewer

141:   Level: advanced

143: .seealso: `ISColoring()`, `ISColoringViewFromOptions()`, `ISColoringDestroy()`, `ISColoringGetIS()`, `MatColoring`
144: @*/
145: PetscErrorCode ISColoringView(ISColoring iscoloring, PetscViewer viewer)
146: {
147:   PetscInt  i;
148:   PetscBool iascii;
149:   IS       *is;

151:   PetscFunctionBegin;
152:   PetscAssertPointer(iscoloring, 1);
153:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(iscoloring->comm, &viewer));

156:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
157:   if (iascii) {
158:     MPI_Comm    comm;
159:     PetscMPIInt size, rank;

161:     PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
162:     PetscCallMPI(MPI_Comm_size(comm, &size));
163:     PetscCallMPI(MPI_Comm_rank(comm, &rank));
164:     PetscCall(PetscViewerASCIIPrintf(viewer, "ISColoring Object: %d MPI processes\n", size));
165:     PetscCall(PetscViewerASCIIPrintf(viewer, "ISColoringType: %s\n", ISColoringTypes[iscoloring->ctype]));
166:     PetscCall(PetscViewerASCIIPushSynchronized(viewer));
167:     PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Number of colors %" PetscInt_FMT "\n", rank, iscoloring->n));
168:     PetscCall(PetscViewerFlush(viewer));
169:     PetscCall(PetscViewerASCIIPopSynchronized(viewer));
170:   }

172:   PetscCall(ISColoringGetIS(iscoloring, PETSC_USE_POINTER, PETSC_IGNORE, &is));
173:   for (i = 0; i < iscoloring->n; i++) PetscCall(ISView(iscoloring->is[i], viewer));
174:   PetscCall(ISColoringRestoreIS(iscoloring, PETSC_USE_POINTER, &is));
175:   PetscFunctionReturn(PETSC_SUCCESS);
176: }

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

181:   Not Collective

183:   Input Parameter:
184: . iscoloring - the coloring context

186:   Output Parameters:
187: + n      - number of nodes
188: . nc     - number of colors
189: - colors - color for each node

191:   Level: advanced

193:   Notes:
194:   Do not free the `colors` array.

196:   The `colors` array will only be valid for the lifetime of the `ISColoring`

198: .seealso: `ISColoring`, `ISColoringValue`, `ISColoringRestoreIS()`, `ISColoringView()`, `ISColoringGetIS()`
199: @*/
200: PetscErrorCode ISColoringGetColors(ISColoring iscoloring, PetscInt *n, PetscInt *nc, const ISColoringValue **colors)
201: {
202:   PetscFunctionBegin;
203:   PetscAssertPointer(iscoloring, 1);

205:   if (n) *n = iscoloring->N;
206:   if (nc) *nc = iscoloring->n;
207:   if (colors) *colors = iscoloring->colors;
208:   PetscFunctionReturn(PETSC_SUCCESS);
209: }

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

214:   Collective

216:   Input Parameters:
217: + iscoloring - the coloring context
218: - 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

220:   Output Parameters:
221: + nn   - number of index sets in the coloring context
222: - isis - array of index sets

224:   Level: advanced

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

229: .seealso: `ISColoring`, `IS`, `ISColoringRestoreIS()`, `ISColoringView()`, `ISColoringGetColoring()`, `ISColoringGetColors()`
230: @*/
231: PetscErrorCode ISColoringGetIS(ISColoring iscoloring, PetscCopyMode mode, PetscInt *nn, IS *isis[])
232: {
233:   PetscFunctionBegin;
234:   PetscAssertPointer(iscoloring, 1);

236:   if (nn) *nn = iscoloring->n;
237:   if (isis) {
238:     if (!iscoloring->is) {
239:       PetscInt        *mcolors, **ii, nc = iscoloring->n, i, base, n = iscoloring->N;
240:       ISColoringValue *colors = iscoloring->colors;
241:       IS              *is;

243:       if (PetscDefined(USE_DEBUG)) {
244:         for (i = 0; i < n; i++) PetscCheck(((PetscInt)colors[i]) < nc, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coloring is our of range index %d value %d number colors %d", (int)i, (int)colors[i], (int)nc);
245:       }

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

251:       PetscCall(PetscMalloc1(nc, &ii));
252:       PetscCall(PetscMalloc1(n, &ii[0]));
253:       for (i = 1; i < nc; i++) ii[i] = ii[i - 1] + mcolors[i - 1];
254:       PetscCall(PetscArrayzero(mcolors, nc));

256:       if (iscoloring->ctype == IS_COLORING_GLOBAL) {
257:         PetscCallMPI(MPI_Scan(&iscoloring->N, &base, 1, MPIU_INT, MPI_SUM, iscoloring->comm));
258:         base -= iscoloring->N;
259:         for (i = 0; i < n; i++) ii[colors[i]][mcolors[colors[i]]++] = i + base; /* global idx */
260:       } else if (iscoloring->ctype == IS_COLORING_LOCAL) {
261:         for (i = 0; i < n; i++) ii[colors[i]][mcolors[colors[i]]++] = i; /* local idx */
262:       } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not provided for this ISColoringType type");

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

267:       if (mode != PETSC_OWN_POINTER) iscoloring->is = is;
268:       *isis = is;
269:       PetscCall(PetscFree(ii[0]));
270:       PetscCall(PetscFree(ii));
271:       PetscCall(PetscFree(mcolors));
272:     } else {
273:       *isis = iscoloring->is;
274:     }
275:   }
276:   PetscFunctionReturn(PETSC_SUCCESS);
277: }

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

282:   Collective

284:   Input Parameters:
285: + iscoloring - the coloring context
286: . mode       - who retains ownership of the is
287: - is         - array of index sets

289:   Level: advanced

291: .seealso: `ISColoring()`, `IS`, `ISColoringGetIS()`, `ISColoringView()`, `PetscCopyMode`
292: @*/
293: PetscErrorCode ISColoringRestoreIS(ISColoring iscoloring, PetscCopyMode mode, IS *is[])
294: {
295:   PetscFunctionBegin;
296:   PetscAssertPointer(iscoloring, 1);

298:   /* currently nothing is done here */
299:   PetscFunctionReturn(PETSC_SUCCESS);
300: }

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

305:   Collective

307:   Input Parameters:
308: + comm    - communicator for the processors creating the coloring
309: . ncolors - max color value
310: . n       - number of nodes on this processor
311: . colors  - array containing the colors for this MPI rank, color numbers begin at 0, for each local node
312: - mode    - see `PetscCopyMode` for meaning of this flag.

314:   Output Parameter:
315: . iscoloring - the resulting coloring data structure

317:   Options Database Key:
318: . -is_coloring_view - Activates `ISColoringView()`

320:   Level: advanced

322:   Notes:
323:   By default sets coloring type to  `IS_COLORING_GLOBAL`

325: .seealso: `ISColoring`, `ISColoringValue`, `MatColoringCreate()`, `ISColoringView()`, `ISColoringDestroy()`, `ISColoringSetType()`
326: @*/
327: PetscErrorCode ISColoringCreate(MPI_Comm comm, PetscInt ncolors, PetscInt n, const ISColoringValue colors[], PetscCopyMode mode, ISColoring *iscoloring)
328: {
329:   PetscMPIInt size, rank, tag;
330:   PetscInt    base, top, i;
331:   PetscInt    nc, ncwork;
332:   MPI_Status  status;

334:   PetscFunctionBegin;
335:   if (ncolors != PETSC_DECIDE && ncolors > IS_COLORING_MAX) {
336:     PetscCheck(ncolors <= PETSC_MAX_UINT16, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Max color value exceeds %d limit. This number is unrealistic. Perhaps a bug in code? Current max: %d user requested: %" PetscInt_FMT, PETSC_MAX_UINT16, PETSC_IS_COLORING_MAX, ncolors);
337:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Max color value exceeds limit. Perhaps reconfigure PETSc with --with-is-color-value-type=short? Current max: %d user requested: %" PetscInt_FMT, PETSC_IS_COLORING_MAX, ncolors);
338:   }
339:   PetscCall(PetscNew(iscoloring));
340:   PetscCall(PetscCommDuplicate(comm, &(*iscoloring)->comm, &tag));
341:   comm = (*iscoloring)->comm;

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

346:   /* should use MPI_Scan() */
347:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
348:   if (rank == 0) {
349:     base = 0;
350:     top  = n;
351:   } else {
352:     PetscCallMPI(MPI_Recv(&base, 1, MPIU_INT, rank - 1, tag, comm, &status));
353:     top = base + n;
354:   }
355:   if (rank < size - 1) PetscCallMPI(MPI_Send(&top, 1, MPIU_INT, rank + 1, tag, comm));

357:   /* compute the total number of colors */
358:   ncwork = 0;
359:   for (i = 0; i < n; i++) {
360:     if (ncwork < colors[i]) ncwork = colors[i];
361:   }
362:   ncwork++;
363:   PetscCall(MPIU_Allreduce(&ncwork, &nc, 1, MPIU_INT, MPI_MAX, comm));
364:   PetscCheck(nc <= ncolors, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of colors passed in %" PetscInt_FMT " is less than the actual number of colors in array %" PetscInt_FMT, ncolors, nc);
365:   (*iscoloring)->n     = nc;
366:   (*iscoloring)->is    = NULL;
367:   (*iscoloring)->N     = n;
368:   (*iscoloring)->refct = 1;
369:   (*iscoloring)->ctype = IS_COLORING_GLOBAL;
370:   if (mode == PETSC_COPY_VALUES) {
371:     PetscCall(PetscMalloc1(n, &(*iscoloring)->colors));
372:     PetscCall(PetscArraycpy((*iscoloring)->colors, colors, n));
373:     (*iscoloring)->allocated = PETSC_TRUE;
374:   } else if (mode == PETSC_OWN_POINTER) {
375:     (*iscoloring)->colors    = (ISColoringValue *)colors;
376:     (*iscoloring)->allocated = PETSC_TRUE;
377:   } else {
378:     (*iscoloring)->colors    = (ISColoringValue *)colors;
379:     (*iscoloring)->allocated = PETSC_FALSE;
380:   }
381:   PetscCall(ISColoringViewFromOptions(*iscoloring, NULL, "-is_coloring_view"));
382:   PetscCall(PetscInfo(0, "Number of colors %" PetscInt_FMT "\n", nc));
383:   PetscFunctionReturn(PETSC_SUCCESS);
384: }

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

390:   Collective

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

396:   Output Parameter:
397: . rows - contains new numbers from remote or local

399:   Level: advanced

401:   Developer Notes:
402:   This manual page is incomprehensible and still needs to be fixed

404: .seealso: [](sec_scatter), `IS`, `MatPartitioningCreate()`, `ISPartitioningToNumbering()`, `ISPartitioningCount()`
405: @*/
406: PetscErrorCode ISBuildTwoSided(IS ito, IS toindx, IS *rows)
407: {
408:   const PetscInt *ito_indices, *toindx_indices;
409:   PetscInt       *send_indices, rstart, *recv_indices, nrecvs, nsends;
410:   PetscInt       *tosizes, *fromsizes, i, j, *tosizes_tmp, *tooffsets_tmp, ito_ln;
411:   PetscMPIInt    *toranks, *fromranks, size, target_rank, *fromperm_newtoold, nto, nfrom;
412:   PetscLayout     isrmap;
413:   MPI_Comm        comm;
414:   PetscSF         sf;
415:   PetscSFNode    *iremote;

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

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

498:   Collective

500:   Input Parameter:
501: . part - a partitioning as generated by `MatPartitioningApply()` or `MatPartitioningApplyND()`

503:   Output Parameter:
504: . is - on each processor the index set that defines the global numbers
505:          (in the new numbering) for all the nodes currently (before the partitioning)
506:          on that processor

508:   Level: advanced

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

513: .seealso: [](sec_scatter), `IS`, `MatPartitioningCreate()`, `AOCreateBasic()`, `ISPartitioningCount()`
514: @*/
515: PetscErrorCode ISPartitioningToNumbering(IS part, IS *is)
516: {
517:   MPI_Comm        comm;
518:   IS              ndorder;
519:   PetscInt        i, np, npt, n, *starts = NULL, *sums = NULL, *lsizes = NULL, *newi = NULL;
520:   const PetscInt *indices = NULL;

522:   PetscFunctionBegin;
524:   PetscAssertPointer(is, 2);
525:   /* see if the partitioning comes from nested dissection */
526:   PetscCall(PetscObjectQuery((PetscObject)part, "_petsc_matpartitioning_ndorder", (PetscObject *)&ndorder));
527:   if (ndorder) {
528:     PetscCall(PetscObjectReference((PetscObject)ndorder));
529:     *is = ndorder;
530:     PetscFunctionReturn(PETSC_SUCCESS);
531:   }

533:   PetscCall(PetscObjectGetComm((PetscObject)part, &comm));
534:   /* count the number of partitions, i.e., virtual processors */
535:   PetscCall(ISGetLocalSize(part, &n));
536:   PetscCall(ISGetIndices(part, &indices));
537:   np = 0;
538:   for (i = 0; i < n; i++) np = PetscMax(np, indices[i]);
539:   PetscCall(MPIU_Allreduce(&np, &npt, 1, MPIU_INT, MPI_MAX, comm));
540:   np = npt + 1; /* so that it looks like a MPI_Comm_size output */

542:   /*
543:         lsizes - number of elements of each partition on this particular processor
544:         sums - total number of "previous" nodes for any particular partition
545:         starts - global number of first element in each partition on this processor
546:   */
547:   PetscCall(PetscMalloc3(np, &lsizes, np, &starts, np, &sums));
548:   PetscCall(PetscArrayzero(lsizes, np));
549:   for (i = 0; i < n; i++) lsizes[indices[i]]++;
550:   PetscCall(MPIU_Allreduce(lsizes, sums, np, MPIU_INT, MPI_SUM, comm));
551:   PetscCallMPI(MPI_Scan(lsizes, starts, np, MPIU_INT, MPI_SUM, comm));
552:   for (i = 0; i < np; i++) starts[i] -= lsizes[i];
553:   for (i = 1; i < np; i++) {
554:     sums[i] += sums[i - 1];
555:     starts[i] += sums[i - 1];
556:   }

558:   /*
559:       For each local index give it the new global number
560:   */
561:   PetscCall(PetscMalloc1(n, &newi));
562:   for (i = 0; i < n; i++) newi[i] = starts[indices[i]]++;
563:   PetscCall(PetscFree3(lsizes, starts, sums));

565:   PetscCall(ISRestoreIndices(part, &indices));
566:   PetscCall(ISCreateGeneral(comm, n, newi, PETSC_OWN_POINTER, is));
567:   PetscCall(ISSetPermutation(*is));
568:   PetscFunctionReturn(PETSC_SUCCESS);
569: }

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

575:   Collective

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

581:   Output Parameter:
582: . count - array of length size, to contain the number of elements assigned
583:         to each partition, where size is the number of partitions generated
584:          (see notes below).

586:   Level: advanced

588:   Notes:
589:   By default the number of partitions generated (and thus the length
590:   of count) is the size of the communicator associated with `IS`,
591:   but it can be set by `MatPartitioningSetNParts()`.

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

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

597: .seealso: [](sec_scatter), `IS`, `MatPartitioningCreate()`, `AOCreateBasic()`, `ISPartitioningToNumbering()`,
598:           `MatPartitioningSetNParts()`, `MatPartitioningApply()`, `MatPartitioningApplyND()`
599: @*/
600: PetscErrorCode ISPartitioningCount(IS part, PetscInt len, PetscInt count[])
601: {
602:   MPI_Comm        comm;
603:   PetscInt        i, n, *lsizes;
604:   const PetscInt *indices;
605:   PetscMPIInt     npp;

607:   PetscFunctionBegin;
608:   PetscCall(PetscObjectGetComm((PetscObject)part, &comm));
609:   if (len == PETSC_DEFAULT) {
610:     PetscMPIInt size;
611:     PetscCallMPI(MPI_Comm_size(comm, &size));
612:     len = (PetscInt)size;
613:   }

615:   /* count the number of partitions */
616:   PetscCall(ISGetLocalSize(part, &n));
617:   PetscCall(ISGetIndices(part, &indices));
618:   if (PetscDefined(USE_DEBUG)) {
619:     PetscInt np = 0, npt;
620:     for (i = 0; i < n; i++) np = PetscMax(np, indices[i]);
621:     PetscCall(MPIU_Allreduce(&np, &npt, 1, MPIU_INT, MPI_MAX, comm));
622:     np = npt + 1; /* so that it looks like a MPI_Comm_size output */
623:     PetscCheck(np <= len, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Length of count array %" PetscInt_FMT " is less than number of partitions %" PetscInt_FMT, len, np);
624:   }

626:   /*
627:         lsizes - number of elements of each partition on this particular processor
628:         sums - total number of "previous" nodes for any particular partition
629:         starts - global number of first element in each partition on this processor
630:   */
631:   PetscCall(PetscCalloc1(len, &lsizes));
632:   for (i = 0; i < n; i++) {
633:     if (indices[i] > -1) lsizes[indices[i]]++;
634:   }
635:   PetscCall(ISRestoreIndices(part, &indices));
636:   PetscCall(PetscMPIIntCast(len, &npp));
637:   PetscCall(MPIU_Allreduce(lsizes, count, npp, MPIU_INT, MPI_SUM, comm));
638:   PetscCall(PetscFree(lsizes));
639:   PetscFunctionReturn(PETSC_SUCCESS);
640: }

642: /*@
643:   ISAllGather - Given an index set `IS` on each processor, generates a large
644:   index set (same on each processor) by concatenating together each
645:   processors index set.

647:   Collective

649:   Input Parameter:
650: . is - the distributed index set

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

655:   Level: intermediate

657:   Notes:
658:   `ISAllGather()` is clearly not scalable for large index sets.

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

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

668: .seealso: [](sec_scatter), `IS`, `ISCreateGeneral()`, `ISCreateStride()`, `ISCreateBlock()`
669: @*/
670: PetscErrorCode ISAllGather(IS is, IS *isout)
671: {
672:   PetscInt       *indices, n, i, N, step, first;
673:   const PetscInt *lindices;
674:   MPI_Comm        comm;
675:   PetscMPIInt     size, *sizes = NULL, *offsets = NULL, nn;
676:   PetscBool       stride;

678:   PetscFunctionBegin;
680:   PetscAssertPointer(isout, 2);

682:   PetscCall(PetscObjectGetComm((PetscObject)is, &comm));
683:   PetscCallMPI(MPI_Comm_size(comm, &size));
684:   PetscCall(ISGetLocalSize(is, &n));
685:   PetscCall(PetscObjectTypeCompare((PetscObject)is, ISSTRIDE, &stride));
686:   if (size == 1 && stride) { /* should handle parallel ISStride also */
687:     PetscCall(ISStrideGetInfo(is, &first, &step));
688:     PetscCall(ISCreateStride(PETSC_COMM_SELF, n, first, step, isout));
689:   } else {
690:     PetscCall(PetscMalloc2(size, &sizes, size, &offsets));

692:     PetscCall(PetscMPIIntCast(n, &nn));
693:     PetscCallMPI(MPI_Allgather(&nn, 1, MPI_INT, sizes, 1, MPI_INT, comm));
694:     offsets[0] = 0;
695:     for (i = 1; i < size; i++) {
696:       PetscInt s = offsets[i - 1] + sizes[i - 1];
697:       PetscCall(PetscMPIIntCast(s, &offsets[i]));
698:     }
699:     N = offsets[size - 1] + sizes[size - 1];

701:     PetscCall(PetscMalloc1(N, &indices));
702:     PetscCall(ISGetIndices(is, &lindices));
703:     PetscCallMPI(MPI_Allgatherv((void *)lindices, nn, MPIU_INT, indices, sizes, offsets, MPIU_INT, comm));
704:     PetscCall(ISRestoreIndices(is, &lindices));
705:     PetscCall(PetscFree2(sizes, offsets));

707:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, N, indices, PETSC_OWN_POINTER, isout));
708:   }
709:   PetscFunctionReturn(PETSC_SUCCESS);
710: }

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

716:   Collective

718:   Input Parameters:
719: + comm     - communicator to share the indices
720: . n        - local size of set
721: - lindices - local colors

723:   Output Parameters:
724: + outN       - total number of indices
725: - outindices - all of the colors

727:   Level: intermediate

729:   Note:
730:   `ISAllGatherColors()` is clearly not scalable for large index sets.

732: .seealso: `ISCOloringValue`, `ISColoring()`, `ISCreateGeneral()`, `ISCreateStride()`, `ISCreateBlock()`, `ISAllGather()`
733: @*/
734: PetscErrorCode ISAllGatherColors(MPI_Comm comm, PetscInt n, ISColoringValue *lindices, PetscInt *outN, ISColoringValue *outindices[])
735: {
736:   ISColoringValue *indices;
737:   PetscInt         i, N;
738:   PetscMPIInt      size, *offsets = NULL, *sizes = NULL, nn = n;

740:   PetscFunctionBegin;
741:   PetscCallMPI(MPI_Comm_size(comm, &size));
742:   PetscCall(PetscMalloc2(size, &sizes, size, &offsets));

744:   PetscCallMPI(MPI_Allgather(&nn, 1, MPI_INT, sizes, 1, MPI_INT, comm));
745:   offsets[0] = 0;
746:   for (i = 1; i < size; i++) offsets[i] = offsets[i - 1] + sizes[i - 1];
747:   N = offsets[size - 1] + sizes[size - 1];
748:   PetscCall(PetscFree2(sizes, offsets));

750:   PetscCall(PetscMalloc1(N + 1, &indices));
751:   PetscCallMPI(MPI_Allgatherv(lindices, (PetscMPIInt)n, MPIU_COLORING_VALUE, indices, sizes, offsets, MPIU_COLORING_VALUE, comm));

753:   *outindices = indices;
754:   if (outN) *outN = N;
755:   PetscFunctionReturn(PETSC_SUCCESS);
756: }

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

762:   Collective

764:   Input Parameters:
765: + is   - the index set
766: . nmin - the first index desired in the local part of the complement
767: - 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`)

769:   Output Parameter:
770: . isout - the complement

772:   Level: intermediate

774:   Notes:
775:   The communicator for `isout` is the same as for the input `is`

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

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

782: .seealso: [](sec_scatter), `IS`, `ISCreateGeneral()`, `ISCreateStride()`, `ISCreateBlock()`, `ISAllGather()`
783: @*/
784: PetscErrorCode ISComplement(IS is, PetscInt nmin, PetscInt nmax, IS *isout)
785: {
786:   const PetscInt *indices;
787:   PetscInt        n, i, j, unique, cnt, *nindices;
788:   PetscBool       sorted;

790:   PetscFunctionBegin;
792:   PetscAssertPointer(isout, 4);
793:   PetscCheck(nmin >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nmin %" PetscInt_FMT " cannot be negative", nmin);
794:   PetscCheck(nmin <= nmax, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nmin %" PetscInt_FMT " cannot be greater than nmax %" PetscInt_FMT, nmin, nmax);
795:   PetscCall(ISSorted(is, &sorted));
796:   PetscCheck(sorted, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Index set must be sorted");

798:   PetscCall(ISGetLocalSize(is, &n));
799:   PetscCall(ISGetIndices(is, &indices));
800:   if (PetscDefined(USE_DEBUG)) {
801:     for (i = 0; i < n; i++) {
802:       PetscCheck(indices[i] >= nmin, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Index %" PetscInt_FMT "'s value %" PetscInt_FMT " is smaller than minimum given %" PetscInt_FMT, i, indices[i], nmin);
803:       PetscCheck(indices[i] < nmax, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Index %" PetscInt_FMT "'s value %" PetscInt_FMT " is larger than maximum given %" PetscInt_FMT, i, indices[i], nmax);
804:     }
805:   }
806:   /* Count number of unique entries */
807:   unique = (n > 0);
808:   for (i = 0; i < n - 1; i++) {
809:     if (indices[i + 1] != indices[i]) unique++;
810:   }
811:   PetscCall(PetscMalloc1(nmax - nmin - unique, &nindices));
812:   cnt = 0;
813:   for (i = nmin, j = 0; i < nmax; i++) {
814:     if (j < n && i == indices[j]) do {
815:         j++;
816:       } while (j < n && i == indices[j]);
817:     else nindices[cnt++] = i;
818:   }
819:   PetscCheck(cnt == nmax - nmin - unique, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of entries found in complement %" PetscInt_FMT " does not match expected %" PetscInt_FMT, cnt, nmax - nmin - unique);
820:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is), cnt, nindices, PETSC_OWN_POINTER, isout));
821:   PetscCall(ISRestoreIndices(is, &indices));
822:   PetscFunctionReturn(PETSC_SUCCESS);
823: }