Матрицы

#include <vector>
#include <iostream>
#include <chrono>
#include <utility>
#include <stdexcept>
using namespace std;

class Matrix {
   int n,m;
   vector<vector<double>> data;
public:
   Matrix(int n, int m):n(n),m(m),data(n,vector<double>(m,0.0)) {}
   friend Matrix operator*(const Matrix &a, const Matrix &b);
   pair<int,int> size() const { return {n,m}; }
};
int main()
{
   Matrix a(1000,1000),b(1000,1000);
   auto start = std::chrono::high_resolution_clock::now();
   Matrix c=a*b;
   auto end = std::chrono::high_resolution_clock::now();
   std::cout << "Time " << std::chrono::duration_cast<std::chrono::milliseconds>(end-start).count() << " ms\n";
}

Реализация по определению операции:

Matrix operator*(const Matrix &a, const Matrix &b) {
   if (a.m!=b.n) throw runtime_error("неверные размеры");
   Matrix c(a.n,b.m);
   for (int i = 0; i < c.n; i++)
      for (int j = 0; j < c.m; j++)
         for (int k = 0; k < a.m; k++)
            c.data[i][j]+=a.data[i][k]*b.data[k][j];
   return c;
}

Перестановка циклов позволяет компилятору вынести некоторые выражения из цикла. Также на время влияет близость элементов в памяти из-за кэша процессора.

Matrix operator*(const Matrix &a, const Matrix &b) {
   if (a.m!=b.n) throw runtime_error("неверные размеры");
   Matrix c(a.n,b.m);
   for (int i = 0; i < c.n; i++)
     for (int k = 0; k < a.m; k++)
       for (int j = 0; j < c.m; j++)
            c.data[i][j]+=a.data[i][k]*b.data[k][j];
   return c;
}

Компилятор оптимизирует эти циклы так:

Matrix operator*(const Matrix &a, const Matrix &b) {
   if (a.m!=b.n) throw runtime_error("неверные размеры");
   Matrix c(a.n,b.m);
   for (int i = 0; i < c.n; i++) 
   {  auto &ci=c.data[i];
      for (int k = 0; k < a.m; k++)
      { double aik=a.data[i][k];
        auto &bk=b.data[k];
        for (int j = 0; j < c.m; j++)
          ci[j]+=aik*bk[j];
      }
   }
   return c;
}

Для разрежённых матриц может оказаться полезным выполнить предварительное транспонирование перед умножением

// транспонирование
SparseMatrix SparseMatrix::transp() const {
  SparseMatrix r(m,n);
  for(int i=0;i<n;++i)
    for(auto &p: data[i])
      r.data[p.first].emplace_back(i,p.second); // добавить в конец
  return r;
}
// умножение (friend для SparseMatrix)
SparseMatrix operator*(const SparseMatrix & a, const SparseMatrix & b) { 
  if (a.m!=b.n) throw runtime_error("неверные размеры");
  SparseMatrix c(a.n,b.m);
  SparseMatrix bt(b.transp());
  for (int i = 0; i < c.n; i++)
    for (int j = 0; j < c.m; j++) {
      double r=0;
// слияние за O(q), где q - количество ненулевых элементов в строке, аналог for(int k = 0; k < a.m; k++) r+=a[i,k]*bt[j,k];
      auto it1=a.data[i].begin();
      auto it2=bt.data[j].begin();
      while(it1!=a.data[i].end() && it2!=bt.data[j].end()) {
        if(it1->first==it2->first) { r+=it1->second*it2->second; ++it1; ++it2; }
        else if(it1->first<it2->first) ++it1;
        else ++it2;
      }
      if(r)
        c.data[i].emplace_back(j,r);
  }
  return c;
}

Для плотных матриц транспонирование также немного улучшает время за счет последовательного доступа к памяти, но изменение порядка циклов улучшает больше.