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:   PetscFunctionBegin;
 11:   coloring->refct++;
 12:   PetscFunctionReturn(PETSC_SUCCESS);
 13: }

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

 18:    Collective

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

 24:    Level: intermediate

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

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

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

 40: /*@C

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

 44:    Collective

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

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

 52:    Level: intermediate

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

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

 66:    Collective

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

 71:    Level: advanced

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

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

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

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

100:   Collective

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

107:   Level: intermediate

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

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

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

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

136:    Collective

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

142:    Level: advanced

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

152:   PetscFunctionBegin;
154:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(iscoloring->comm, &viewer));

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

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

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

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

182:    Not Collective

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

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

192:    Level: advanced

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

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

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

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

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

215:    Collective

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

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

225:    Level: advanced

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

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

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

244:       if (PetscDefined(USE_DEBUG)) {
245:         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);
246:       }

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

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

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

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

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

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

283:    Collective

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

290:    Level: advanced

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

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

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

306:     Collective

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

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

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

321:    Level: advanced

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

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

335:   PetscFunctionBegin;
336:   if (ncolors != PETSC_DECIDE && ncolors > IS_COLORING_MAX) {
337:     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?\nCurrent max: %d user requested: %" PetscInt_FMT, PETSC_MAX_UINT16, PETSC_IS_COLORING_MAX, ncolors);
338:     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);
339:   }
340:   PetscCall(PetscNew(iscoloring));
341:   PetscCall(PetscCommDuplicate(comm, &(*iscoloring)->comm, &tag));
342:   comm = (*iscoloring)->comm;

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

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

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

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

391:     Collective

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

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

400:    Level: advanced

402:    Developer Note:
403:    This manual page is incomprehensible and still needs to be fixed

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

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

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

499:     Collective

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

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

509:    Level: advanced

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

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

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

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

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

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

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

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

576:     Collective

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

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

587:    Level: advanced

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

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

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

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

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

616:   /* count the number of partitions */
617:   PetscCall(ISGetLocalSize(part, &n));
618:   PetscCall(ISGetIndices(part, &indices));
619:   if (PetscDefined(USE_DEBUG)) {
620:     PetscInt np = 0, npt;
621:     for (i = 0; i < n; i++) np = PetscMax(np, indices[i]);
622:     PetscCall(MPIU_Allreduce(&np, &npt, 1, MPIU_INT, MPI_MAX, comm));
623:     np = npt + 1; /* so that it looks like a MPI_Comm_size output */
624:     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);
625:   }

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

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

648:     Collective

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

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

656:     Level: intermediate

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

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

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

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

679:   PetscFunctionBegin;

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

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

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

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

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

717:     Collective

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

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

728:     Level: intermediate

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

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

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

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

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

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

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

763:     Collective

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

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

773:     Level: intermediate

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

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

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

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

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

799:   PetscCall(ISGetLocalSize(is, &n));
800:   PetscCall(ISGetIndices(is, &indices));
801:   if (PetscDefined(USE_DEBUG)) {
802:     for (i = 0; i < n; i++) {
803:       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);
804:       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);
805:     }
806:   }
807:   /* Count number of unique entries */
808:   unique = (n > 0);
809:   for (i = 0; i < n - 1; i++) {
810:     if (indices[i + 1] != indices[i]) unique++;
811:   }
812:   PetscCall(PetscMalloc1(nmax - nmin - unique, &nindices));
813:   cnt = 0;
814:   for (i = nmin, j = 0; i < nmax; i++) {
815:     if (j < n && i == indices[j]) do {
816:         j++;
817:       } while (j < n && i == indices[j]);
818:     else nindices[cnt++] = i;
819:   }
820:   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);
821:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is), cnt, nindices, PETSC_OWN_POINTER, isout));
822:   PetscCall(ISRestoreIndices(is, &indices));
823:   PetscFunctionReturn(PETSC_SUCCESS);
824: }