1: #include <petsc/private/matimpl.h>
  2: #include <petsc/private/isimpl.h>

  4: static PetscErrorCode MatColoringApply_Natural(MatColoring mc, ISColoring *iscoloring)
  5: {
  6:   PetscInt         start, end, i, bs = 1, n;
  7:   ISColoringValue *colors;
  8:   MPI_Comm         comm;
  9:   PetscBool        flg1, flg2;
 10:   Mat              mat     = mc->mat;
 11:   Mat              mat_seq = mc->mat;
 12:   PetscMPIInt      size;
 13:   ISColoring       iscoloring_seq;
 14:   ISColoringValue *colors_loc;
 15:   PetscInt         rstart, rend, N_loc, nc;

 17:   PetscFunctionBegin;
 18:   /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
 19:   PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATSEQBAIJ, &flg1));
 20:   PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATMPIBAIJ, &flg2));
 21:   if (flg1 || flg2) PetscCall(MatGetBlockSize(mat, &bs));

 23:   PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
 24:   PetscCallMPI(MPI_Comm_size(comm, &size));
 25:   if (size > 1) {
 26:     /* create a sequential iscoloring on all processors */
 27:     PetscCall(MatGetSeqNonzeroStructure(mat, &mat_seq));
 28:   }

 30:   PetscCall(MatGetSize(mat_seq, &n, NULL));
 31:   PetscCall(MatGetOwnershipRange(mat_seq, &start, &end));
 32:   n = n / bs;
 33:   PetscCheck(n <= IS_COLORING_MAX - 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Maximum color size exceeded");

 35:   start = start / bs;
 36:   end   = end / bs;
 37:   PetscCall(PetscMalloc1(end - start + 1, &colors));
 38:   for (i = start; i < end; i++) colors[i - start] = (ISColoringValue)i;
 39:   PetscCall(ISColoringCreate(comm, n, end - start, colors, PETSC_OWN_POINTER, iscoloring));

 41:   if (size > 1) {
 42:     PetscCall(MatDestroySeqNonzeroStructure(&mat_seq));

 44:     /* convert iscoloring_seq to a parallel iscoloring */
 45:     iscoloring_seq = *iscoloring;
 46:     rstart         = mat->rmap->rstart / bs;
 47:     rend           = mat->rmap->rend / bs;
 48:     N_loc          = rend - rstart; /* number of local nodes */

 50:     /* get local colors for each local node */
 51:     PetscCall(PetscMalloc1(N_loc + 1, &colors_loc));
 52:     for (i = rstart; i < rend; i++) colors_loc[i - rstart] = iscoloring_seq->colors[i];
 53:     /* create a parallel iscoloring */
 54:     nc = iscoloring_seq->n;
 55:     PetscCall(ISColoringCreate(comm, nc, N_loc, colors_loc, PETSC_OWN_POINTER, iscoloring));
 56:     PetscCall(ISColoringDestroy(&iscoloring_seq));
 57:   }
 58:   PetscFunctionReturn(PETSC_SUCCESS);
 59: }

 61: /*MC
 62:   MATCOLORINGNATURAL - implements a trivial coloring routine with one color per column

 64:   Level: beginner

 66:   Note:
 67:   Using this coloring would be extremely inefficient but it is useful for testing

 69: .seealso: `MatColoring`, `MatColoringType`, `MatColoringCreate()`, `MatColoring`, `MatColoringSetType()`, `MatColoringType`
 70: M*/
 71: PETSC_EXTERN PetscErrorCode MatColoringCreate_Natural(MatColoring mc)
 72: {
 73:   PetscFunctionBegin;
 74:   mc->data                = NULL;
 75:   mc->ops->apply          = MatColoringApply_Natural;
 76:   mc->ops->view           = NULL;
 77:   mc->ops->destroy        = NULL;
 78:   mc->ops->setfromoptions = NULL;
 79:   PetscFunctionReturn(PETSC_SUCCESS);
 80: }