[OpenMP] 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?

Paste this command into your terminal:

The password to connect is

Question 1: Sparse matrix-vector multiplication

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

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

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

typedef struct            sparse_matrix  sparse_matrix_t;

/* ---------------------------------------------------------------------
 * 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

  ( double*                 lhs,
    const sparse_matrix_t*  mat,
    const double*           rhs,
    int chunksize )
{
  const int  nrows = mat->nrows;

  #pragma omp parallel shared(mat, lhs, rhs)
  {
    #pragma omp for schedule(dynamic, chunksize)
    for (int irow = 0; irow < nrows; irow++ )
    {
      double sum  = 0.0;
      #pragma omp simd reduction(+:sum)
      for (int i = mat->offset[irow]; i < mat->offset[irow + 1]; i++ )
      {
        sum += mat->value[i] * rhs[mat->index[i]];
      }
      lhs[irow] = sum;
    }
  }
}

/* ---------------------------------------------------------------------
 * 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";
  double*          lhs   = NULL;
  double*          rhs   = NULL;

  sparse_matrix_t  matrix;
  int              iiter;

  printf ( "----- Matrix-vector multiplication: OpenMP -----\n" );
  // time
  t0_read = wtime ();

  int max_num_threads = omp_get_max_threads();
  //max_num_threads = 3;
  omp_set_num_threads(max_num_threads);
  printf("used number of threads: %d \n", max_num_threads);

  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(*lhs) );
  rhs = (double*) malloc ( matrix.ncols * sizeof(*rhs) );

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

  // set up initial right hand side
  vector_init(rhs,matrix.ncols);

  // time
  t1_read = wtime() - t0_read;

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

  // chunksize for parallel for loop
  int chunksize = (matrix.nrows / (max_num_threads*10)) +1;

  // time
  t0_matvec = wtime ();

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

    matvec ( lhs, &matrix, rhs, chunksize );

    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] );
    }
    */
  }

  // time
  t1_matvec = wtime() - t0_matvec;

  free   ( lhs );
  free   ( rhs );

  sparse_matrix_destroy ( &matrix );

  // time
  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.

Author(s)	Erik-Jan Lingen, Matthias Möller
Deadline	Keine Frist
Submission limit	No limitation

Information

Sign in