Parallel Matrix-Vector Multiplication with MPI Collectives

1. Decomposing the Problem

Let's first recall the sequential algorithm for matrix-vector multiplication, where A is an $m\times n$ matrix and x is an $n\times 1$ vector:

$y_i={\sum }_{j=0}^{n-1}{A}_{i,j}\cdot x_j\text{for }i=0,\dots ,m-1$

The key insight for parallelization is that the calculation of each element $y_i$ of the result vector is independent of all other elements $y_k$ (where $k\ne i$). This is a perfect example of data parallelism.

Our approach will be to give each process a "stripe" of rows from the matrix A. This is known as row-wise block-striped decomposition.

This leads to a critical question: what data does each process need to perform its work?

• Unique Data: Each process needs its own unique block of rows from matrix A.
• Shared Data: To compute the inner product for any given row, a process needs access to the entire vector x. Therefore, the vector x must be available to all processes.

With this understanding, we can outline a simple and effective parallel algorithm.

As you saw in the slides, the algorithm is:

1. The root process (e.g., rank 0) initializes the full matrix A and vector x.
2. The root scatters the rows of matrix A among all processes.
3. The root broadcasts the vector x to all processes.
4. Each process computes the matrix-vector product for its local rows, producing a partial result vector.
5. The root gathers the partial results from all processes.
6. The root now holds the complete final result vector y.

2. Mapping the Algorithm to MPI Collectives

Now, let's map this algorithm to the MPI functions you know. Your proficiency in C++ and understanding of memory layouts will be helpful here, especially when dealing with the 2D matrix.

Step 1: Handling the Matrix Data

A 2D matrix in C++ is often represented as a pointer-to-pointer (int**) or allocated as a single contiguous block of memory (a 1D array). For MPI, a single contiguous block is strongly preferred because MPI scatter/gather operations work most easily with flat memory buffers.

You can allocate an $m\times n$ matrix as a 1D array and access element $(i,j)$ using A[i * n + j].

Step 2: Distributing the Data (MPI_Scatter & MPI_Bcast)

• Distributing Matrix A: The root process has the full matrix. We want to send a unique block of rows to each process. This is a perfect job for MPI_Scatter. If we have p processes and an m x n matrix, each process will receive m/p rows, which corresponds to (m/p) * n elements from the flattened 1D array.
• Distributing Vector x: Every process needs an identical copy of the entire vector x. The right tool for this is MPI_Bcast.

Step 3: Local Computation

Each process now has its local rows of A (let's call it A_local) and the full vector x. It can now perform a standard sequential matrix-vector multiplication to produce its piece of the result, y_local.

Step 4: Collecting the Results (MPI_Gather)

Each process holds a partial result vector y_local. To form the final vector y, we need to collect these pieces and concatenate them in rank order. This is the exact function of MPI_Gather. The root process will receive all the y_local vectors into its main y buffer.

This complete process is a beautiful demonstration of how different MPI collectives work together to solve a larger problem.

3. Hands-on Practice: Implementing the Algorithm

Now it's time to translate this design into a working C++ program. Your task is to implement parallel matrix-vector multiplication using the scatter-broadcast-gather pattern.

Your Task:

Assume you have a matrix A of size M x N and a vector x of size N. For simplicity, ensure that M is perfectly divisible by the number of MPI processes.

1. Initialization (Root Process):

   • Let the root process (rank 0) define the matrix dimensions (e.g., M=12, N=4).
   • Dynamically allocate a 1D array for matrix A of size M * N and a 1D array for vector x of size N.
   • Initialize A and x with some values (e.g., A[i*N + j] = i and x[j] = 1.0).
   • Allocate a 1D array for the final result vector y of size M.

2. Data Distribution (All Processes):

   • Calculate rows_per_proc = M / world_size.
   • Each process allocates a local buffer A_local of size rows_per_proc * N to receive its rows of the matrix.
   • Each process allocates a buffer for the full vector x of size N.
   • Use MPI_Scatter to distribute A from the root. The sendcount and recvcount will be rows_per_proc * N.
   • Use MPI_Bcast to broadcast x from the root. The count will be N.

3. Local Computation (All Processes):

   • Each process allocates a local result buffer y_local of size rows_per_proc.
   • Perform the matrix-vector multiplication using A_local and x, storing the result in y_local.

   // Pseudocode for local computation
   for (int i = 0; i < rows_per_proc; ++i) {
       y_local[i] = 0.0;
       for (int j = 0; j < N; ++j) {
           y_local[i] += A_local[i * N + j] * x[j];
       }
   }
   

4. Result Collection (All Processes):

   • Use MPI_Gather to collect the y_local results from all processes into the y vector on the root. The sendcount and recvcount will be rows_per_proc.

5. Verification (Root Process):

   • The root process should print the final vector y to verify the result.

Take your time with this exercise. Managing the buffer sizes and counts for the different MPI calls is the main challenge and a crucial skill to develop.

4. A Glimpse at a More Scalable Approach

The row-wise method is simple and effective, but it has a scalability limitation. As you increase the number of processes, every single process still receives a full copy of vector x. The MPI_Bcast operation and the memory required to store x on every node can become a bottleneck for very large-scale problems.

A more advanced technique, called 2D checkerboard decomposition, partitions both the matrix A and the vectors x and y.

This 2D approach is more complex to implement, requiring you to organize processes into a grid using MPI's Cartesian communicator features (MPI_Cart_create, MPI_Comm_split). However, it is far more scalable and is the standard for high-performance linear algebra libraries.

Conclusion

Excellent work today! You have successfully designed and implemented a parallel algorithm for a non-trivial and important computational kernel. This exercise ties together much of what you've learned about MPI collectives.

Key Takeaways:

• Algorithm is Key: Parallel programming starts with decomposing the problem. For matrix-vector multiplication, a row-wise decomposition is a natural starting point.
• Map Algorithm to MPI Tools: We mapped our communication needs directly to MPI functions: MPI_Scatter for decomposing the matrix, MPI_Bcast for replicating shared data, and MPI_Gather for aggregating the results.
• Memory Layout Matters: Using a contiguous 1D array for your matrix simplifies communication with MPI.
• Scalability is a Goal: While our approach works, we noted its limitations and saw a preview of a more scalable 2D decomposition strategy.

Preview of the Next Lesson

In our discussion of the 2D decomposition, we mentioned that it reduces communication overhead. But how do we quantify this? In the next lesson, we will learn how to model communication overhead and analyze the performance of our parallel applications. We'll explore strategies like message aggregation to see how we can further optimize our MPI programs by minimizing the time spent sending and receiving data.