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:   /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
 18:   PetscObjectTypeCompare((PetscObject)mat, MATSEQBAIJ, &flg1);
 19:   PetscObjectTypeCompare((PetscObject)mat, MATMPIBAIJ, &flg2);
 20:   if (flg1 || flg2) MatGetBlockSize(mat, &bs);

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

 29:   MatGetSize(mat_seq, &n, NULL);
 30:   MatGetOwnershipRange(mat_seq, &start, &end);
 31:   n = n / bs;

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

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

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

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

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

 63:   Level: beginner

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

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