SALSA Analysis Modules
Functions
jpl.c File Reference

Functions for computing the JPL multicolour structure of a matrix. More...

#include <stdlib.h>
#include "anamod.h"
#include "anamodsalsamodules.h"
#include "anamatrix.h"
#include "petscerror.h"
#include "petscmat.h"
#include "petscis.h"
#include "src/mat/impls/aij/mpi/mpiaij.h"
#include "src/mat/impls/aij/seq/aij.h"
Include dependency graph for jpl.c:

Go to the source code of this file.

Functions

static PetscErrorCode LookAtUncolouredVar (int compare, int this_var, Mat A, PetscScalar *clr_array, PetscScalar *ran_array, PetscReal this_rand, int *neighb, int *nneigh, PetscBool *colour_now)
static PetscErrorCode JonesPlassmannColouring (Mat A, PetscBool *flg)
static PetscErrorCode NColours (AnaModNumericalProblem prob, AnalysisItem *rv, int *lv, PetscBool *flg)
static PetscErrorCode ColourSizes (AnaModNumericalProblem prob, AnalysisItem *rv, int *lv, PetscBool *flg)
static PetscErrorCode ColourOffsets (AnaModNumericalProblem prob, AnalysisItem *rv, int *lv, PetscBool *flg)
static PetscErrorCode Colours (AnaModNumericalProblem prob, AnalysisItem *rv, int *lv, PetscBool *flg)
PetscErrorCode RegisterJPLModules (void)

Detailed Description

Functions for computing the JPL multicolour structure of a matrix.

Jones-Plassmann multi-colouring

Usage

Activate this module with

PetscErrorCode RegisterJPLModules();

Compute these elements with

ComputeQuantity("jpl","element",A,(void*)&res,&flg);

Available elements are:

This routine works in parallel; the stored result is distributed over all participating processors.

References

@article{JoPl:heuristic,
author = {M.T. Jones and P.E. Plassmann},
title = {A parallel graph coloring heuristic},
journal = SJSSC,
volume = {14},
year = {1993},
keywords = {incomplete factorization, multicolouring, matrix graph}
}

@inproceedings{JoPl:parallelunstructured,
author = {M.T. Jones and P.E. Plassmann},
title = {Parallel solution of unstructed, sparse systems of linear equations},
booktitle = {Proceedings of the Sixth SIAM conference on
	Parallel Processing for Scientific Computing},
editor = {R.F. Sincovec and D.E. Keyes and M.R. Leuze and L.R. Petzold
	 and D.A. Reed},
publisher = {SIAM},
address = {Philadelphia},
pages = {471--475},
year = {1993}

  \section Disclaimer

  Unabashed use of British spelling of `colour' needs no disclaimer.
}

Definition in file jpl.c.


Function Documentation

static PetscErrorCode ColourOffsets ( AnaModNumericalProblem  prob,
AnalysisItem rv,
int *  lv,
PetscBool *  flg 
) [static]

Compute an array that has the offsets of the locations of the colours, on this processor. If computing the colouring is done in parallel, then the multicolouring is stored in parallel: it is not gathered onto any processor.

Definition at line 451 of file jpl.c.

References GetDataID(), HASTOEXIST, id, AnalysisItem::ii, and JonesPlassmannColouring().

Referenced by RegisterJPLModules().

{
  Mat A = (Mat)prob;
  int llv=0,*v = NULL; PetscBool has; int id,id2; PetscErrorCode ierr;
  PetscFunctionBegin;

  ierr = GetDataID("jpl","n-colours",&id2,&has); CHKERRQ(ierr);
  HASTOEXIST(has);
  ierr = PetscObjectComposedDataGetInt
    ((PetscObject)A,id2,llv,*flg); CHKERRQ(ierr);
  if (*flg) {
    if (lv) *lv = llv;
    ierr = GetDataID("jpl","colour-offsets",&id,&has); CHKERRQ(ierr);
    HASTOEXIST(has);
    ierr = PetscObjectComposedDataGetIntstar
      ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
  }
  if (*flg) rv->ii = v;
  else {
    ierr = JonesPlassmannColouring(A,flg); CHKERRQ(ierr);
    if (*flg) {
      ierr = PetscObjectComposedDataGetIntstar
        ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
      if (*flg) rv->ii = v;
      ierr = PetscObjectComposedDataGetInt
        ((PetscObject)A,id2,llv,*flg); CHKERRQ(ierr);
      if (*flg)
        if (lv) *lv = llv;
    }
  }
  
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

static PetscErrorCode Colours ( AnaModNumericalProblem  prob,
AnalysisItem rv,
int *  lv,
PetscBool *  flg 
) [static]

Compute an array that has the colour sets on this processor. In order to know where a set starts or ends, the result of ColourOffsets() is needed.

If computing the colouring is done in parallel, then the multicolouring is stored in parallel: it is not gathered onto any processor.

Definition at line 497 of file jpl.c.

References GetDataID(), HASTOEXIST, id, AnalysisItem::ii, and JonesPlassmannColouring().

Referenced by RegisterJPLModules().

{
  Mat A = (Mat)prob; int N;
  int *v = NULL; PetscBool has; int id; PetscErrorCode ierr;
  PetscFunctionBegin;
  ierr = MatGetLocalSize(A,&N,PETSC_NULL); CHKERRQ(ierr);
  if (lv) *lv = N;
  ierr = GetDataID("jpl","colours",&id,&has); CHKERRQ(ierr);
  HASTOEXIST(has);
  ierr = PetscObjectComposedDataGetIntstar
    ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
  if (*flg) rv->ii = v;
  else {
    ierr = JonesPlassmannColouring(A,flg); CHKERRQ(ierr);
    if (*flg) {
      ierr = PetscObjectComposedDataGetIntstar
        ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
      if (*flg) rv->ii = v;
    }
  }
  
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

static PetscErrorCode ColourSizes ( AnaModNumericalProblem  prob,
AnalysisItem rv,
int *  lv,
PetscBool *  flg 
) [static]

Compute an array that has the number of points of each colour, on this processor. If computing the colouring is done in parallel, then the multicolouring is stored in parallel: it is not gathered onto any processor.

Definition at line 408 of file jpl.c.

References GetDataID(), HASTOEXIST, id, AnalysisItem::ii, and JonesPlassmannColouring().

Referenced by RegisterJPLModules().

{
  Mat A = (Mat)prob;
  int llv=0,*v = NULL; PetscBool has; int id,id2; PetscErrorCode ierr;
  PetscFunctionBegin;

  ierr = GetDataID("jpl","n-colours",&id2,&has); CHKERRQ(ierr);
  HASTOEXIST(has);
  ierr = PetscObjectComposedDataGetInt
    ((PetscObject)A,id2,llv,*flg); CHKERRQ(ierr);
  if (*flg) {
    if (lv) *lv = llv;
    ierr = GetDataID("jpl","colour-set-sizes",&id,&has); CHKERRQ(ierr);
    HASTOEXIST(has);
    ierr = PetscObjectComposedDataGetIntstar
      ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
  }
  if (*flg) rv->ii = v;
  else {
    ierr = JonesPlassmannColouring(A,flg); CHKERRQ(ierr);
    if (*flg) {
      ierr = PetscObjectComposedDataGetIntstar
        ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
      if (*flg) rv->ii = v;
      ierr = PetscObjectComposedDataGetInt
        ((PetscObject)A,id2,llv,*flg); CHKERRQ(ierr);
      if (*flg)
        if (lv) *lv = llv;
    }
  }
  
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

static PetscErrorCode JonesPlassmannColouring ( Mat  A,
PetscBool *  flg 
) [static]

Compute a multicolour ordering of a matrix, in parallel.

This is the main computational routine.

Definition at line 129 of file jpl.c.

References GetDataID(), id, LookAtUncolouredVar(), and MatrixComputeQuantity().

Referenced by ColourOffsets(), Colours(), ColourSizes(), and NColours().

{
  MPI_Comm comm; const MatType type;
  Mat dia,off; Vec rand,ran_bord, clrs,clrs_bkp,clr_bord; VecScatter sctr;
  PetscScalar *ran_array,*bran_array,*clr_array,*clr_bkp_array,*bclr_array;
  int *neighb,ncoloured,total_coloured,pass,max_clr,
    global_maxclr,local_size,total_size, id;
  PetscBool parallel,sequential; PetscErrorCode ierr;
  
  PetscFunctionBegin;

  *flg = PETSC_FALSE;

  /*
   * What is the matrix type? We only handle AIJ (nonblocked) types
   */
  ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
  ierr = MatGetType(A,&type); CHKERRQ(ierr);
  ierr = PetscStrcmp(type,MATSEQAIJ,&sequential); CHKERRQ(ierr);
  ierr = PetscStrcmp(type,MATMPIAIJ,&parallel); CHKERRQ(ierr);
  if (!(sequential || parallel)) {
    *flg = PETSC_FALSE;
    PetscFunctionReturn(0);
  }

  /*
   * Set up the matrices and vectors
   */
  if (parallel) {
    Mat_MPIAIJ *aij = (Mat_MPIAIJ *) A->data;
    dia=aij->A; off=aij->B; clr_bord = aij->lvec; sctr = aij->Mvctx;
    ierr = VecDuplicate(aij->lvec,&ran_bord); CHKERRQ(ierr);
  } else {
    dia = A; off = NULL; clr_bord = NULL; sctr = NULL;
  }
  ierr = MatGetSize(A,&total_size,PETSC_NULL); CHKERRQ(ierr);
  ierr = MatGetLocalSize(A,&local_size,PETSC_NULL); CHKERRQ(ierr);

  /* A temporary array for storing neighbour information */
  {
    int l1,l2,len; PetscBool f;
    ierr = MatrixComputeQuantity
      (dia,"structure","max-nnzeros-per-row",(AnalysisItem*)&l1,PETSC_NULL,
       &f); CHKERRQ(ierr);
    if (parallel) {
      ierr = MatrixComputeQuantity
        (off,"structure","max-nnzeros-per-row",(AnalysisItem*)&l2,PETSC_NULL,
         &f); CHKERRQ(ierr);
      if (f) len=l1+l2; else len = local_size;
    } else len = l1;
    ierr = PetscMalloc((len+1)*sizeof(int),&neighb); CHKERRQ(ierr);
  }

  /*
   * Create the vectors for random and colours
   */
  {
    if (parallel) {
      ierr = VecCreateMPI(comm,local_size,PETSC_DECIDE,&rand); CHKERRQ(ierr);
    } else {
      ierr = VecCreateSeq(MPI_COMM_SELF,local_size,&rand); CHKERRQ(ierr);
    }
    {
      PetscRandom rctx;
      ierr = PetscRandomCreate(MPI_COMM_SELF,&rctx); CHKERRQ(ierr);
      ierr = PetscRandomSetType(rctx,PETSCRAND48); CHKERRQ(ierr);
      ierr = VecSetRandom(rand,rctx); CHKERRQ(ierr);
      ierr = PetscRandomDestroy(&rctx); CHKERRQ(ierr);
    }
    ierr = VecDuplicate(rand,&clrs); CHKERRQ(ierr);
    ierr = VecDuplicate(rand,&clrs_bkp); CHKERRQ(ierr);
    ierr = VecSet(clrs,0.); CHKERRQ(ierr);
    if (parallel) {
      ierr = VecScatterBegin
        (sctr,rand,ran_bord,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
      ierr = VecScatterEnd
        (sctr,rand,ran_bord,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
      ierr = VecGetArray(ran_bord,&bran_array); CHKERRQ(ierr);
      ierr = VecGetArray(clr_bord,&bclr_array); CHKERRQ(ierr);
    }
    ierr = VecGetArray(rand,&ran_array); CHKERRQ(ierr);
    ierr = VecGetArray(clrs,&clr_array); CHKERRQ(ierr);
    ierr = VecGetArray(clrs_bkp,&clr_bkp_array); CHKERRQ(ierr);
  }

  /*
   * Here is the main loop.
   * Iterate until everything is coloured.
   */
  ncoloured = 0; max_clr = 0;
  for (pass=0; ; pass++) {
    PetscBool busy = PETSC_FALSE,gbusy; int var;
    /*printf(".. pass %d :",pass);*/
    if (max_clr>pass)
      SETERRQ2(MPI_COMM_WORLD,1,"This can not happen: maxclr=%d in pass %d",max_clr,pass);
    ierr = PetscMemcpy
      (clr_bkp_array,clr_array,local_size*sizeof(PetscScalar)); CHKERRQ(ierr);

    /*
     * Loop over all variables in every pass
     */
    for (var=0; var<local_size; var++) {
      int nn = 0; PetscBool now;

      /* Skip already coloured variables */
      if (clr_array[var]>0) continue;

      /* All others, look at all connections */
      ierr = LookAtUncolouredVar
        (var,var,dia,clr_bkp_array,ran_array,ran_array[var],
         neighb,&nn,&now); CHKERRQ(ierr);
      if (now && parallel) {
        if (parallel) {
          ierr = LookAtUncolouredVar
            (-1,var,off,bclr_array,bran_array,ran_array[var],
             neighb,&nn,&now); CHKERRQ(ierr);
        }
      }

      /* If this variable has the highest value among uncoloured neighbours,
       * colour it now. */
      if (now) {
        int clr=1,i;
      attempt:
        for (i=0; i<nn; i++)
          if (neighb[i]==clr) {clr++; goto attempt;}
        clr_array[var] = (PetscScalar) clr;
        if (clr>max_clr) max_clr = clr;
        ncoloured++; busy = PETSC_TRUE; 
        /*printf(" %d->%d",var,clr);*/
      } else {
        /*printf(" %d.",var);*/
      }
    } /*printf("\n");*/

    /* Exit if all points are coloured */
    ierr = MPI_Allreduce
      ((void*)&ncoloured,(void*)&total_coloured,1,MPI_INT,MPI_SUM,comm);
    if (total_coloured==total_size) break;

    /* If nothing was coloured in this pass, we're stuck */
    ierr = MPI_Allreduce
      ((void*)&busy,(void*)&gbusy,1,MPI_INT,MPI_MAX,comm);
    if (!gbusy) SETERRQ(MPI_COMM_WORLD,1,"We are stuck");

    if (parallel) {
      ierr = VecScatterBegin
        (sctr,clrs,clr_bord,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
      ierr = VecScatterEnd
        (sctr,clrs,clr_bord,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
    }
  }

  /*
   * Store the number of colours
   */
  ierr = MPI_Allreduce
    ((void*)&max_clr,(void*)&global_maxclr,1,MPI_INT,MPI_MAX,comm);
  ierr = GetDataID("jpl","n-colours",&id,PETSC_NULL); CHKERRQ(ierr);
  ierr = PetscObjectComposedDataSetInt
    ((PetscObject)A,id,global_maxclr); CHKERRQ(ierr);

  /*
   * Count how many points of each colour
   */
  {
    int *colour_count,*colour_offsets,*ptr,*colours,clr,i;
    ierr = PetscMalloc
      ((global_maxclr+1)*sizeof(int),&colour_count); CHKERRQ(ierr);
    colour_count[0] = global_maxclr;
    for (clr=0; clr<global_maxclr; clr++)
      colour_count[clr+1] = 0;
    for (i=0; i<local_size; i++)
      colour_count[(int)(clr_array[i])]++; /* colours are from 1 */
    /*
    printf("coloured:");
    for (i=0; i<local_size; i++) printf(" %e",clr_array[i]);
    printf("\n");
    printf("counts:");
    for (i=0; i<colour_count[0]; i++) printf(" %d",colour_count[i+1]);
    printf("\n");
    */
    ierr = GetDataID("jpl","colour-set-sizes",&id,PETSC_NULL); CHKERRQ(ierr);
    ierr = PetscObjectComposedDataSetIntstar
      ((PetscObject)A,id,colour_count); CHKERRQ(ierr);

    ierr = PetscMalloc
      ((global_maxclr+1)*sizeof(int),&colour_offsets); CHKERRQ(ierr);
    colour_offsets[0] = global_maxclr;
    colour_offsets[1] = 0;
    for (clr=1; clr<global_maxclr; clr++)
      colour_offsets[1+clr] = colour_offsets[clr]+colour_count[clr];
    /*
    printf("offsets:");
    for (i=0; i<colour_offsets[0]; i++) printf(" %d",colour_offsets[i+1]);
    printf("\n");
    */
    ierr = GetDataID("jpl","colour-offsets",&id,PETSC_NULL); CHKERRQ(ierr);
    ierr = PetscObjectComposedDataSetIntstar
      ((PetscObject)A,id,colour_offsets); CHKERRQ(ierr);

    ierr = PetscMalloc((local_size+1)*sizeof(int),&colours); CHKERRQ(ierr);
    colours[0] = local_size;
    ierr = PetscMalloc(global_maxclr*sizeof(int),&ptr); CHKERRQ(ierr);
    PetscMemcpy(ptr,colour_offsets+1,global_maxclr*sizeof(int));
    for (i=0; i<local_size; i++) {
      int c=(int)(clr_array[i]);
      colours[1+ptr[c-1]++] = i;
    }
    ierr = PetscFree(ptr); CHKERRQ(ierr);
    /*
    printf("colours:");
    for (i=0; i<colours[0]; i++) printf(" %d",colours[i+1]);
    printf("\n");
    */
    ierr = GetDataID("jpl","colours",&id,PETSC_NULL); CHKERRQ(ierr);
    ierr = PetscObjectComposedDataSetIntstar
      ((PetscObject)A,id,colours); CHKERRQ(ierr);
  }

  /* Clean up */
  ierr = PetscFree(neighb); CHKERRQ(ierr);
  ierr = VecRestoreArray(rand,&ran_array); CHKERRQ(ierr);
  ierr = VecDestroy(&rand); CHKERRQ(ierr);
  ierr = VecRestoreArray(clrs,&clr_array); CHKERRQ(ierr);
  ierr = VecDestroy(&clrs); CHKERRQ(ierr);
  ierr = VecRestoreArray(clrs_bkp,&clr_bkp_array); CHKERRQ(ierr);
  ierr = VecDestroy(&clrs_bkp); CHKERRQ(ierr);
  if (parallel) {
    ierr = VecRestoreArray(ran_bord,&bran_array); CHKERRQ(ierr);
    ierr = VecDestroy(&ran_bord); CHKERRQ(ierr);
    ierr = VecRestoreArray(clr_bord,&bclr_array); CHKERRQ(ierr);
    /* clr_bord is the aij->lvec array; not to be destroyed */
  }

  *flg = PETSC_TRUE;


  PetscFunctionReturn(0);
}

Here is the call graph for this function:

static PetscErrorCode LookAtUncolouredVar ( int  compare,
int  this_var,
Mat  A,
PetscScalar *  clr_array,
PetscScalar *  ran_array,
PetscReal  this_rand,
int *  neighb,
int *  nneigh,
PetscBool *  colour_now 
) [static]

Investigate the local and remote connections of a variable to see whether it needs to be coloured, and if so, with what colour.

Definition at line 77 of file jpl.c.

Referenced by JonesPlassmannColouring().

{

  int nnz,j; const int *idx; PetscErrorCode ierr;

  PetscFunctionBegin;

  *colour_now = PETSC_TRUE;

  ierr = MatGetRow(A,this_var,&nnz,&idx,PETSC_NULL); CHKERRQ(ierr);
  for (j=0; j<nnz; j++) {
    /* loop over all points connected to 'this_var' */
    int other,clr;
    other=idx[j]; clr = (int)(clr_array[other]);

    /* no need to look at ourselves */
    if (other==compare) continue;
    /* this really shouldn't happen */
    /*
    if (ran_array[other]==this_rand)
      SETERRQ3(MPI_COMM_WORLD,1,"Use a better random: %e @ %d,%d",this_rand,compare,other);
    */
    /*
     * If the other variable is uncoloured and has a higher random value,
     * we can not colour the current variable.
     */
    if ( (clr==0 && ran_array[other]>this_rand) ||
         (ran_array[other]==this_rand && other>this_var) ) {
      *colour_now = PETSC_FALSE;
      goto exit;
    }
    /*
     * Otherwise, we take note of the colour of the other variable
     * for the computation of the colour of this variable later
     */
    neighb[(*nneigh)++] = clr;
  }
 exit:
  ierr = MatRestoreRow(A,this_var,&nnz,&idx,PETSC_NULL); CHKERRQ(ierr);
  
  PetscFunctionReturn(0);
}
static PetscErrorCode NColours ( AnaModNumericalProblem  prob,
AnalysisItem rv,
int *  lv,
PetscBool *  flg 
) [static]

Compute the global number of colours. Note that any given processor does not need to have all the colours.

Definition at line 377 of file jpl.c.

References GetDataID(), HASTOEXIST, AnalysisItem::i, id, and JonesPlassmannColouring().

Referenced by RegisterJPLModules().

{
  Mat A = (Mat)prob;
  int v = 0; PetscBool has; int id; PetscErrorCode ierr;
  PetscFunctionBegin;
  ierr = GetDataID("jpl","n-colours",&id,&has); CHKERRQ(ierr);
  HASTOEXIST(has);
  ierr = PetscObjectComposedDataGetInt
    ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
  if (*flg) rv->i = v;
  else {
    ierr = JonesPlassmannColouring(A,flg); CHKERRQ(ierr);
    if (*flg) {
      ierr = PetscObjectComposedDataGetInt
        ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
      if (*flg) rv->i = v;
    }
  }
  
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

PetscErrorCode RegisterJPLModules ( void  )

Definition at line 523 of file jpl.c.

References ANALYSISINTARRAY, ANALYSISINTEGER, ColourOffsets(), Colours(), ColourSizes(), NColours(), and RegisterModule().

Referenced by AnaModRegisterSalsaModules(), and AnaModRegisterStandardModules().

{
  PetscErrorCode ierr;
  PetscFunctionBegin;

  ierr = RegisterModule
    ("jpl","n-colours",ANALYSISINTEGER,&NColours); CHKERRQ(ierr);
  ierr = RegisterModule
    ("jpl","colour-set-sizes",ANALYSISINTARRAY,&ColourSizes); CHKERRQ(ierr);
  ierr = RegisterModule
    ("jpl","colour-offsets",ANALYSISINTARRAY,&ColourOffsets); CHKERRQ(ierr);
  ierr = RegisterModule
    ("jpl","colours",ANALYSISINTARRAY,&Colours); CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

Here is the call graph for this function: