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) {
 21:     MatGetBlockSize(mat,&bs);
 22:   }

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

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

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

 44:   if (size > 1) {
 45:     MatDestroySeqNonzeroStructure(&mat_seq);

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

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

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

 69:   Level: beginner

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