Actual source code: natural.c

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