Начало
Учебные материалы
Матрицы
Представление матриц Операции над матрицами LU-разложение и вычисление обратной матрицы Умножение матриц Алгоритм умножения матриц Штрассена
Умножение матриц
#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";
}
```
--
Реализация по определению операции:
```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 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;
}
```
--
Перестановка циклов позволяет компилятору вынести некоторые выражения из цикла. Также на время влияет близость элементов в памяти из-за кэша процессора.
```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;
}
```
Компилятор оптимизирует эти циклы так:
```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;
}
```
Для разрежённых матриц может оказаться полезным выполнить предварительное транспонирование перед умножением
```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;
}
```
Для плотных матриц транспонирование также немного улучшает время за счет последовательного доступа к памяти,
но изменение порядка циклов улучшает больше.
