[Threads] Parallel sparse matrix-vector multiplication

This material is part of the course Practical Introduction to Parallel Programming and meant for educational purposes only.

Implement a program that reads a sparse matrix \(\mathbf{A}\) from a file, stores it in a Compressed Sparse Row (CSR) format, and performs a series of parallel, sparse matrix-vector products of the form

\begin{equation*} \mathbf{x}_{i} = \mathbf{A} \, \mathbf{x}_{i-1} \end{equation*}

The first vector \(\mathbf{x}_{0}\) is generated randomly.

This code given below implements the sequential version of the program. It defines a struct sparse_matrix_t that represents a sparse matrix.

The program reads a sparse matrix from a file named matrix.dat.

You should (try to) implement three parallel versions of sparse matrix-vector product: one with MPI; one with OpenMP and one with threads.

Hints

Partition the matrix row wise and assign each batch of rows to a different thread or process.

Read the matrix on the master thread/process.

Start with the parallel implementation based on OpenMP as this is the most straightforward.

Compare the parallel performance of each implementation. Is there a significant difference?

Copiez-collez cette commande dans votre terminal :

Le mot de pase pour se connecter est

Question 1: Sparse matrix-vector multiplication

#include <stdio.h>
#include <assert.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>
#include <pthread.h>
#include <math.h>

/* ---------------------------------------------------------------------
 * types
 * ------------------------------------------------------------------ */

struct                    sparse_matrix
{
  int                       nrows;
  int                       ncols;
  int*                      offset;
  int*                      index;
  double*                   value;
};

typedef struct            sparse_matrix  sparse_matrix_t;


struct                    matvec_arg
{
  int row_start;
  int row_end;
  pthread_mutex_t *mutex;
};

typedef struct matvec_arg       matvec_arg_t;


/* ---------------------------------------------------------------------
 * global variables
 * ------------------------------------------------------------------ */

sparse_matrix_t matrix;
double* lhs;
double* rhs;

/* ---------------------------------------------------------------------
 * wtime()
 *
 *   Returns the wall-clock time in seconds since some moment in
 *   the past.
 * ------------------------------------------------------------------ */

double                    wtime

  ( void )

{
  static int             ini = 0;
  static struct timeval  tv0;

  struct timeval         tv;

  gettimeofday ( &tv, NULL );

  if ( ! ini )
  {
    tv0 = tv;
    ini = 1;
  }

  if ( tv0.tv_usec > tv.tv_usec )
  {
    tv.tv_usec  = 1000000 - (tv0.tv_usec - tv.tv_usec);
    tv.tv_sec  -= 1;
  }
  else
  {
    tv.tv_usec -= tv0.tv_usec;
  }

  return ((double) tv.tv_sec + 1.0e-6 * tv.tv_usec);
}


/* ---------------------------------------------------------------------
 * sparse_matrix_init()
 *
 *   Initialises a sparse matrix, given the number of rows, columns
 *   and non-zero values.
 * ------------------------------------------------------------------ */

void                      sparse_matrix_init

  ( sparse_matrix_t*        matrix,
    int                     nrows,
    int                     ncols,
    int                     nnz )

{
  assert ( nrows >= 0 );
  assert ( ncols >= 0 );
  assert ( nnz   >= 0 );

  memset ( matrix, 0x0, sizeof(*matrix) );

  matrix->nrows  = nrows;
  matrix->ncols  = ncols;
  matrix->offset = (int*)malloc ( (nrows + 1) * sizeof(*matrix->offset) );
  matrix->index  = (int*)malloc (  nnz        * sizeof(*matrix->index) );
  matrix->value  = (double*)malloc (  nnz        * sizeof(*matrix->value) );

  if ( (matrix->offset == NULL) ||
       (matrix->index  == NULL) ||
       (matrix->value  == NULL) )
  {
    fprintf ( stderr,
              "ERROR: failed to allocate a sparse matrix.\n" );
    exit    ( 1 );
  }

  memset ( matrix->offset, 0x0,
           (nrows + 1) * sizeof(*matrix->offset) );
}

/* ---------------------------------------------------------------------
 * sparse_matrix_destroy()
 *
 *   Destroys a sparse matrix.
 * ------------------------------------------------------------------ */

void                      sparse_matrix_destroy

  ( sparse_matrix_t*        matrix )

{
  free   ( matrix->offset );
  free   ( matrix->index );
  free   ( matrix->value );
  memset ( matrix, 0x0, sizeof(*matrix) );
}

/* ---------------------------------------------------------------------
 * sparse_matrix_read()
 *
 *   Reads a sparse matrix from a file in the Matrix Market
 *   coordinate format.
 * ------------------------------------------------------------------ */

void                      sparse_matrix_read

  ( sparse_matrix_t*        matrix,
    const char*             fname )

{
  FILE*   file = fopen ( fname, "r" );
  int     ierr = 0;

  int     nrows;
  int     ncols;
  int     nnz;
  int     irow;
  int     jcol;
  double  val;
  int     i, k, n;

  if ( ! file )
  {
    fprintf ( stderr, "ERROR: failed to open the file `%s\'.\n",
              fname );
    exit    ( 1 );
  }

  if ( fscanf( file, " %d %d %d", &nrows, &ncols, &nnz ) != 3 )
  {
    ierr = 1;
    goto error;
  }

  if ( (nrows < 0) || (ncols < 0) || (nnz < 0) )
  {
    ierr = 2;
    goto error;
  }

  sparse_matrix_init ( matrix, nrows, ncols, nnz );

  /* First pass: count the number of entries per row. */

  for ( i = 0; i < nnz; i++ )
  {
    if ( fscanf( file, " %d %d %le", &irow, &jcol, &val ) != 3 )
    {
      ierr = 3;
      goto error;
    }

    if ( (irow < 1) || (irow > nrows) ||
         (jcol < 1) || (jcol > ncols) )
    {
      ierr = 4;
      goto error;
    }

    matrix->offset[irow - 1]++;
  }

  fclose ( file );

  file = NULL;

  /* Prepare the offset array. */

  for ( irow = n = 0; irow <= nrows; irow++ )
  {
    k                    = matrix->offset[irow];
    matrix->offset[irow] = n;
    n                   += k;
  }

  if ( n != nnz )
  {
    ierr = 5;
    goto error;
  }

  /* Second pass: store the values in the sparse matrix. */

  file = fopen ( fname, "r" );

  if ( ! file )
  {
    fprintf ( stderr, "ERROR: failed to open the file `%s\'.\n",
              fname );
    exit    ( 1 );
  }

  if ( fscanf( file, " %d %d %d", &nrows, &ncols, &nnz ) != 3 )
  {
    ierr = 6;
    goto error;
  }

  if ( (nrows != matrix->nrows) ||
       (ncols != matrix->ncols) ||
       (nnz   != matrix->offset[nrows]) )
  {
    ierr = 7;
    goto error;
  }

  for ( i = 0; i < nnz; i++ )
  {
    if ( fscanf( file, " %d %d %le", &irow, &jcol, &val ) != 3 )
    {
      ierr = 8;
      goto error;
    }

    if ( (irow < 1) || (irow > nrows) ||
         (jcol < 1) || (jcol > ncols) )
    {
      ierr = 9;
      goto error;
    }

    k                = matrix->offset[irow - 1]++;
    matrix->index[k] = jcol - 1;
    matrix->value[k] = val;
  }

  fclose ( file );

  file = NULL;

  /* Correct the offset array. */

  for ( irow = n = 0; irow < nrows; irow++ )
  {
    k                    = matrix->offset[irow];
    matrix->offset[irow] = n;
    n                    = k;
  }

  if ( n != nnz )
  {
    ierr = 10;
    goto error;
  }

  /* No error occurred if this point is reached. */

  return;

error:

  fprintf ( stderr, "ERROR: failed to read the file `%s\' "
            "(ierr = %d).\n",
            fname, ierr );
  exit    ( 1 );
}

/* ---------------------------------------------------------------------
 * vector_init()
 *
 *   Fills a vector with random values.
 * ------------------------------------------------------------------ */

void                      vector_init

  ( double*                 vec,
    int                     len )

{
  int  i;

  for ( i = 0; i < len; i++ )
  {
    //vec[i] = rand() / (double) RAND_MAX;
    vec[i] = (double)i / 1;
  }
}

/* ---------------------------------------------------------------------
 * matvec()
 *
 *   Executes a sparse matrix-vector product.
 * ------------------------------------------------------------------ */

void*                     matvec

  ( void*                   arg )

{
  matvec_arg_t* matvec_arg = (matvec_arg_t*) arg;
  int row_start = matvec_arg->row_start;
  int row_end= matvec_arg->row_end;

  const int  nrows = row_end - row_start;

  for (int irow = row_start; irow < row_end; irow++ )
  {
    double sum  = 0.0;
    for (int i = matrix.offset[irow]; i < matrix.offset[irow + 1]; i++ )
    {
      sum += matrix.value[i] * rhs[matrix.index[i]];
    }
    lhs[irow] = sum;
  }

  // ------------------------ with mutex ------------------

  /*
  matvec_arg_t* matvec_arg = (matvec_arg_t*) arg;
  int row_start = matvec_arg->row_start;
  int row_end= matvec_arg->row_end;

  const int  nrows = row_end - row_start;


  double* lhs_loc =  (double*)malloc ( nrows * sizeof(double) );

  for (int irow = row_start; irow < row_end; irow++ )
  {
    double sum  = 0.0;
    for (int i = matrix.offset[irow]; i < matrix.offset[irow + 1]; i++ )
    {
      sum += matrix.value[i] * rhs[matrix.index[i]];
    }
    lhs_loc[irow - row_start] = sum;
  }

  pthread_mutex_lock ( matvec_arg->mutex );
  for (int irow = row_start; irow < row_end; irow++ )
  {
    lhs[irow] = lhs_loc[irow - row_start];
  }
  pthread_mutex_unlock ( matvec_arg->mutex );
  */

  return NULL;
}



/* ---------------------------------------------------------------------
 * main()
 * ------------------------------------------------------------------ */

int                       main()

{

  double           t0_matvec, t0_read, t0_total, t1_matvec, t1_read, t1_total;
  t0_total = wtime ();

  const int        niter = 1000;
  const char*      fname = "matrix.dat";
  int              iiter;

  const int num_threads = 8;

  printf ( "----- Matrix-vector multiplication: Threads -----\n" );
  printf("used number of threads: %d \n", num_threads);

  t0_read = wtime ();

  printf ( "Reading matrix from `%s\' ...\n\n", fname );

  sparse_matrix_read ( &matrix, fname );

  printf ( "  nrows : %d\n",   matrix.nrows );
  printf ( "  ncols : %d\n",   matrix.ncols );
  printf ( "  nnz   : %d\n\n", matrix.offset[matrix.nrows] );

  printf ( "Initialising the first vector ...\n" );


  lhs = (double*)malloc ( matrix.nrows * sizeof(double) );
  rhs = (double*)malloc ( matrix.ncols * sizeof(double) );

  if ( (lhs == NULL) || (rhs == NULL) )
  {
    fprintf ( stderr, "ERROR: failed to allocate the vectors.\n" );
    return 1;
  }

  vector_init(rhs,matrix.ncols);

  t1_read = wtime() - t0_read;

  t0_matvec = wtime ();

  // set matvec arguments: single mutex, seperate row_start, row_end
  pthread_mutex_t mutex;
  pthread_mutex_init ( &mutex, NULL );

  int step = round(matrix.nrows/num_threads);

  matvec_arg_t matvec_args[num_threads];

  for (int i = 0; i < num_threads; i++)
  {
    int row_start = i * step;
    int row_end;
    if (i!= num_threads -1)
    {
      row_end = (i + 1) * step;
    }
    else
    {
      row_end = matrix.nrows;
    }
    matvec_args[i].row_start = row_start;
    matvec_args[i].row_end = row_end;
    matvec_args[i].mutex = &mutex;
  }

  printf ( "Executing %d matrix-vector products ...\n", niter );

  pthread_t tid[num_threads];


  for ( iiter = 0; iiter < niter; iiter++ )
  {
    double*  tmp = rhs;

    for (int j = 0; j < num_threads; j++)
    {
      pthread_create ( &tid[j], NULL, matvec, &matvec_args[j] );
      //printf ( "thread %d spawned \n", j );
    }

    for (int j = 0; j < num_threads; j++)
    {
      pthread_join ( tid[j], NULL );
      //printf ( "thread %d joined \n", j );
    }

    /* Swap the left-hand side and right-hand side vectors. */
    rhs = lhs;
    lhs = tmp;

    /*
    if (iiter <10)
    {
      printf ( "entry 0     of lhs: \t %g \n", lhs[0] );
      printf ( "entry 7000  of lhs: \t %g \n", lhs[7000] );
      printf ( "entry 17000 of lhs: \t %g \n", lhs[17000] );
    }
    */
  }

  t1_matvec = wtime() - t0_matvec;

  pthread_mutex_destroy ( &mutex );

  free   ( lhs );
  free   ( rhs );

  sparse_matrix_destroy ( &matrix );

  t1_total = wtime() - t0_total;

  printf ( "matrix read: \t\t\t done in %g seconds.\n", t1_read );
  printf ( "matrix vector multiplication: \t done in %g seconds.\n", t1_matvec );
  printf ( "total: \t\t\t\t done in %g seconds.\n", t1_total );
  //printf ( "%g \t %g \t %g \n", t1_read, t1_matvec, t1_total);

  return 0;
}

Question 2: Number of MPI processes

Specify the number of MPI processes (1-64) that should be used for running your program (mpirun -np <nproc>)

Question 3: OpenMP environment variables

Specify the OpenMP environment variables that should be used for running your program (mpirun -x OMP_NUM_THREADS=1 -x ...). Use semicolons to separate multiple variables

Question 4: Command Line Argument

Specify the command line arguments that should be passed to your program when it is run. Arguments must be given in quotes and separated by commas. Leave this filed empty if you do not want to specify command line arguments.

Example: "1","int","arg=2" is interpreted as three arguments 1, int, and arg=2.

Auteur(s)	Erik-Jan Lingen, Matthias Möller
Date limite	Pas de date limite
Limite de soumission	Pas de limite

Informations

Se connecter