[MPI] 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 <mpi.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;

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

{
  const int  nrows = mat->nrows;

    for (int irow = 0; irow < nrows; irow++ )
    {
      double sum  = 0.0;
      //iend = mat->offset[irow + 1];
      for (int i = mat->offset[irow]; i < mat->offset[irow + 1]; i++ )
      {
        //jcol = mat->index[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;

  double*          lhs_loc = NULL;
  double*          lhs     = NULL;
  double*          rhs     = NULL;


  sparse_matrix_t  matrix_loc;
  int              iiter;

  int              proc_count;
  int              my_rank;

  int             nrows_global, ncols;
  int             row_start=0, row_end=0, nrows=0, nnz, step;


  // MPI initialization
  MPI_Init      ( NULL, NULL );
  MPI_Comm_size ( MPI_COMM_WORLD, &proc_count );
  MPI_Comm_rank ( MPI_COMM_WORLD, &my_rank );

  int             root = proc_count-1;


  // values for all other processes; only use and write values as root
  int* nnz_loc       = (int*) malloc ( proc_count * sizeof(int) );
  int* row_start_loc = (int*) malloc ( proc_count * sizeof(int) );
  int* row_end_loc   = (int*) malloc ( proc_count * sizeof(int) );
  int* nrows_loc     = (int*) malloc ( proc_count * sizeof(int) );

  t0_read = wtime ();

  // read matrix as root and send parts to other processes
  if (my_rank == root)
  {
    printf ( "----- Matrix-vector multiplication: MPI -----\n" );
    printf("used number of processes: %d \n", proc_count);

    const char*      fname = "matrix.dat";
    sparse_matrix_t  matrix;
    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] );

    nrows_global = matrix.nrows;
    ncols = nrows_global;


    step = round(nrows_global/proc_count);

    // send
    for (int i = 0; i < proc_count; i++)
    {

      if (i != root )
      {
        row_start_loc[i] = i * step;
        row_end_loc[i] = (i + 1) * step;
        nrows_loc[i] = row_end_loc[i] - row_start_loc[i];
        nnz_loc[i] = matrix.offset[row_end_loc[i]] - matrix.offset[row_start_loc[i]];

        MPI_Send(&nrows_global, 1, MPI_INT, i, 0, MPI_COMM_WORLD);
        MPI_Send(&nnz_loc[i], 1, MPI_INT, i, 0, MPI_COMM_WORLD);
      }
      else
      {
        row_start = row_start_loc[i] = i * step;
        row_end = row_end_loc[i] = nrows_global;
        nrows = nrows_loc[i] = row_end_loc[i] - row_start_loc[i];
        nnz = nnz_loc[i] = matrix.offset[row_end_loc[i]] - matrix.offset[row_start_loc[i]];
      }
    }

    // send parts to other processes

    for (int i = 0; i < proc_count; i++)
    {
      int*     mat_index   = (int*)    malloc ( nnz_loc[i]  * sizeof(int) );
      double*  mat_value   = (double*) malloc ( nnz_loc[i]  * sizeof(double) );
      int*     mat_offset  = (int*)    malloc ( nrows_loc[i]     * sizeof(int) );

      for (int j = 0; j < nnz_loc[i]; j++)
      {
        mat_index[j] = matrix.index[ matrix.offset[row_start_loc[i]] + j];
        mat_value[j] = matrix.value[ matrix.offset[row_start_loc[i]] + j];
      }

      for (int j = 0; j < nrows_loc[i]; j++)
      {
        mat_offset[j] = matrix.offset[row_start_loc[i] + j] -  matrix.offset[row_start_loc[i]];
      }

      if (i != root)
      {
        MPI_Send(mat_index, nnz_loc[i], MPI_INT, i, 0, MPI_COMM_WORLD);
        MPI_Send(mat_value, nnz_loc[i], MPI_DOUBLE, i, 0, MPI_COMM_WORLD);
        MPI_Send(mat_offset, nrows_loc[i], MPI_INT, i, 0, MPI_COMM_WORLD);
      }
      else
      {
        sparse_matrix_init ( &matrix_loc, nrows, ncols, nnz );
        matrix_loc.index = mat_index;
        matrix_loc.value = mat_value;
        matrix_loc.offset = mat_offset;
      }
    }

    sparse_matrix_destroy ( &matrix );

    //printf("global matrix destroyed at root\n");

  }

  // receive parts of global matrices from root
  else
  {
    MPI_Recv(&nrows_global, 1, MPI_INT, root, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
    MPI_Recv(&nnz, 1, MPI_INT, root, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);

    // compute own row counts
    step = round(nrows_global/proc_count);
    row_start = my_rank * step;
    row_end = (my_rank + 1) * step;
    nrows = row_end - row_start;
    ncols = nrows_global;

    // initialize own matrix
    sparse_matrix_init ( &matrix_loc, nrows, ncols, nnz );

    // receive matrix values from root
    MPI_Recv(matrix_loc.index, nnz, MPI_INT, root, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
    MPI_Recv(matrix_loc.value, nnz, MPI_DOUBLE, root, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
    MPI_Recv(matrix_loc.offset, step, MPI_INT, root, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);

  }

  printf("proc %d - row_start: %d; row_end: %d; nrows: %d; nnz: %d\n", my_rank, row_start, row_end, nrows, nnz);



  // set up rhs vectors

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

  lhs_loc = (double*) malloc ( nrows        * sizeof(double) );
  lhs     = (double*) malloc ( nrows_global * sizeof(double) );
  rhs     = (double*) malloc ( nrows_global * sizeof(double) );

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

  // initialize rhs vector at root
  if (my_rank == root)
  {
    vector_init(rhs,ncols);
  }

  // broadcast rhs vector to all processes
  MPI_Bcast ( rhs, ncols, MPI_DOUBLE, root, MPI_COMM_WORLD );

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

  // time
  t1_read = wtime () - t0_read;
  t0_matvec = wtime ();

  // do matrix-vector multiplication
  for ( iiter = 0; iiter < niter; iiter++ )
  {
    double*  tmp = rhs;

    matvec ( lhs_loc, &matrix_loc, rhs );


    if (my_rank == root)
    {
      // receive all lhs parts of other processes and write to global lhs
      for (int i=0; i<proc_count-1; i++)
      {
        double*  lhs_part   = (double*) malloc ( nrows_loc[i] * sizeof(double) );
        MPI_Recv(lhs_part, nrows_loc[i], MPI_DOUBLE, i, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
        for (int j=0;j<nrows_loc[i]; j++)
        {
          lhs[row_start_loc[i]+j] = lhs_part[j];
        }
        free(lhs_part);
      }
      // write own lhs part to global lhs
      for (int j=0;j<nrows; j++)
      {
        lhs[row_start+j] = lhs_loc[j];
      }

    }
    else
    {
      // send lhs part to root
      MPI_Send(lhs_loc, nrows, MPI_DOUBLE, root, 0, MPI_COMM_WORLD);
    }

    // broadcast full lhs vector to all processes
    MPI_Bcast ( lhs, nrows_global, MPI_DOUBLE, root, MPI_COMM_WORLD );

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

    /*
    if (iiter <10 && my_rank == root)
    {
      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;


  // free memory
  sparse_matrix_destroy ( &matrix_loc );
  free   ( lhs_loc );
  free   ( lhs );
  free   ( rhs );
  free   ( nnz_loc );
  free   ( row_start_loc );
  free   ( row_end_loc );
  free   ( nrows_loc );

  MPI_Finalize ();

  t1_total = wtime () - t0_total;

  if (my_rank == root)
  {
    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	No deadline
Submission limit	No limitation

Information

Sign in