Multiplying matrices in C++

I have recently come across a (click-bait) post on X claiming that multiplying two 1000×1000 matrices in Python was several times faster than in C++ or other languages. So how do you multiply matrices in C++? Let’s see.

First of all, how do you multiply matrices? It’s pretty simple:

Here, c11 is computed as follows:

c11 = a11*b11 + a12*b21 + a13*b31

Generally, we have:

The simplest way to represent a matrix in C++ is using std::vector. For a matrix of NxN size, accessing element at row i and column j means accessing the vector element at index i*N+j. So we can start with this:

constexpr std::size_t N = 1000;

std::vector<double> A_data(N * N);
std::vector<double> B_data(N * N);
std::vector<double> C_data(N * N);

std::random_device rd{};
std::mt19937 rng(rd());
std::uniform_real_distribution<double> dist(-1000, 1000);

for (auto& v : A_data) v = dist(rng);
for (auto& v : B_data) v = dist(rng);

The naive approach

The simplest way to multiply two matrices is implementing the mathematical operation:

void naive_multiplication(std::vector<double> const & A, std::vector<double> const & B, std::vector<double> & C, std::size_t const N)
{
   for (std::size_t i = 0; i < N; ++i)
   {
      for (std::size_t j = 0; j < N; ++j)
      {
         C[i * N + j] = 0;
         for (std::size_t k = 0; k < N; ++k)
         {
            C[i * N + j] += A[i * N + k] * B[k * N + j];
         }
      }
   }
}

On my machine, running a x64 Release build with Visual Studio 2022 17.20, this takes about 850ms. (This is significantly higher than the number in the X post, which means I just have a slower machine.)

The optimized naive approach

It turns out there is a little optimization you can do to the naive approach: instead of writing to the C matrix after every multiplication, you can compute the sum in a separate variable and only assign that to c(i,j) at the end.

void naive_optimized_multiplication(std::vector<double> const& A, std::vector<double> const& B, std::vector<double>& C, std::size_t const N)
{
   for (std::size_t i = 0; i < N; ++i)
   {
      for (std::size_t j = 0; j < N; ++j)
      {
         double sum = 0;
         for (std::size_t k = 0; k < N; ++k)
         {
            sum += A[i * N + k] * B[k * N + j];
         }

         C[i * N + j] = sum;
      }
   }
}

For me, this takes about 700ms (on average) to execute. This means something around 15% improvement. Why is this better? Because the sum variable is stored in a registry which is very fast to access. On the other hand, writing the the C vector in every iteration forces the compiler to load the existing value from RAM into a registry, add to it, and then store back to RAM. For every element of C there are N (1000 in our case) load and store operations instead of just one.

The parallelized approach

It is possible to parallelize the computation. Parallel algorithms are always more complex but for large data the improvement can be significant. We can compute in parallel the values of each row in the result matrix without worrying about data races.

void parallelized_multiplication(std::vector<double> const& A, std::vector<double> const& B, std::vector<double>& C, std::size_t const N)
{
   std::vector<std::size_t> rows(N, 0);

   std::for_each(
      std::execution::par, 
      rows.begin(), rows.end(),
      [&](std::size_t i)
      {
         for (std::size_t j = 0; j < N; ++j)
         {
            double sum = 0;
            for (std::size_t k = 0; k < N; ++k)
               sum += A[i * N + k] * B[k * N + j];
            C[i * N + j] = sum;
         }
      });
}

On my machine, this takes about 170ms, which is 5 times faster than the naive implementation.

The BLAS approach

BLAS, which stands for Basic Linear Algebra Subprograms, is a specification for low level routines for performing common linear algebra operations such as vector addition, scalar multiplication, dot products, linear combinations, and matrix multiplication.

OpenBLAS is an optimized BLAS library that can be used from C++. We can use this to multiply matrixes. You can check the GitRepo as well as the documentation page for installation instructions and code examples.

Here is how we can multiple matrixes using OpenBLAS:

#pragma comment(lib,"libopenblas.lib")
#include "cblas.h"

void openblas_multiplication(std::vector<double> const& A, std::vector<double> const& B, std::vector<double>& C, blasint const N)
{
   cblas_dgemm(
      CblasRowMajor,     // Row-major layout (C/C++ style)
      CblasNoTrans,      // A not transposed
      CblasNoTrans,      // B not transposed
      N, N, N,
      1.0,               // alpha
      A.data(), N,
      B.data(), N,
      0.0,               // beta
      C.data(), N
   );
}

On my machine, this takes about 10ms. This is more than 15 times faster than the parallelized version and more than 80 times faster than the naive approach.

The C++26 BLAS approach

The C++26 standard library provides an implementation for the BLAS specification in the <linalg> header and std::linalg namespace. At the time of writing this article, no compiler supports it, but you can read more about it here.

Theoretically, the multiplication of two matrixes should look like this in C++26:

#include <linalg>
#include <mdspan>
#include <vector>

void linalg_multiplication(std::vector<double> const& A, std::vector<double> const& B, std::vector<double>& C, std::size_t const N)
{
   // Create mdspan views for the matrices
   std::mdspan A_view(A.data(), N, N);
   std::mdspan B_view(B.data(), N, N);
   std::mdspan C_view(C.data(), N, N);
   
   // Perform matrix multiplication: C = A * B
   std::linalg::matrix_product(A_view, B_view, C_view);
}

However, this might be wrong and should be taken more as a concept. Even if it is actually correct, there is no possibility to measure its performance for the moment. We’ll be able to compare it with OpenBLAS at some point in the future.

Conclusions

Writing a naive implementation for matrix multiplication is very simple. But when you’re looking for performance, optimized libraries such as OpenBLAS are a must.