Вычислить произведение квадратных матриц

#ifndef MATRIX_H
#define MATRIX_H
 
#include <algorithm>
#include <iostream>
#include <random>
#include <vector>
 
class Matrix {
 
public:
    Matrix() = delete;
    Matrix(std::size_t, std::size_t);
    Matrix(std::size_t, std::size_t, int);
    Matrix(const Matrix &);
    ~Matrix();
 
    void randomFill(int, int);
    
    std::vector<int>& operator[](std::size_t);
    const std::vector<int>& operator[](std::size_t) const;
 
    friend std::ostream& operator<<(std::ostream&, const Matrix &);
    friend std::istream& operator>>(std::istream&, Matrix &);
    friend Matrix operator*(const Matrix &, const Matrix &);
 
private:
    std::size_t m_rows;
    std::size_t m_cols;
    std::vector<std::vector<int>> m_storage;
 
};
 
#endif // !MATRIX_H

#include "matrix.h"
 
Matrix::Matrix(std::size_t rows, std::size_t cols)
    : m_rows(rows), m_cols(cols), m_storage(m_rows, std::vector<int>(m_cols, 0)) {
 
    //...
 
}
 
Matrix::Matrix(std::size_t rows, std::size_t cols, int value)
    : m_rows(rows), m_cols(cols), m_storage(m_rows, std::vector<int>(m_cols, value)) {
 
    //...
 
}
 
Matrix::Matrix(const Matrix &matrix) {
 
    m_rows = matrix.m_rows;
    m_cols = matrix.m_cols;
    m_storage.resize(m_rows);
 
    for (std::size_t i = 0; i < m_rows; ++i)
        m_storage[i].resize(m_cols);
 
    for (std::size_t i = 0; i < m_rows; ++i)
        //for (std::size_t j = 0; j < m_cols; ++j)
        //  m_storage[i][j] = matrix[i][j];
        std::copy(matrix[i].begin(), matrix[i].end(), m_storage[i].begin());
 
}
 
Matrix::~Matrix() {
 
    //...
 
}
 
void Matrix::randomFill(int min, int max) {
 
    if (min > max)
        std::swap(min, max);
 
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_int_distribution<int> dist(min, max);
 
    for (auto &row : m_storage)
        for (auto &el : row)
            el = dist(gen);
 
}
 
std::vector<int>& Matrix::operator[](std::size_t index) {
 
    return m_storage[index];
 
}
 
const std::vector<int>& Matrix::operator[](std::size_t index) const {
 
    return m_storage[index];
 
}
 
std::ostream& operator<<(std::ostream& os, const Matrix &matrix) {
 
    for (std::size_t i = 0; i < matrix.m_rows; ++i) {
 
        for (std::size_t j = 0; j < matrix.m_cols; ++j)
            os << matrix[i][j] << ' ';
 
        os << '\n';
 
    }
        
    return os;
 
}
 
std::istream& operator>>(std::istream &is, Matrix &matrix) {
 
    for (auto &row : matrix.m_storage)
        for (auto &el : row)
            is >> el;
 
    return is;
 
}
 
Matrix operator*(const Matrix &m1, const Matrix &m2) {
 
    if (m1.m_rows != m2.m_cols)
        throw std::runtime_error("Matrices are not agreed");
 
    Matrix result(m1.m_rows, m1.m_rows);
 
    for (std::size_t i = 0; i < m1.m_rows; ++i)
        for (std::size_t j = 0; j < m2.m_cols; ++j)
            for (std::size_t k = 0; k < m1.m_cols; ++k)
                result[i][j] += m1[i][k] * m2[k][j];
 
    return result;
 
}

#include "matrix.h"
 
int main() {
 
    Matrix m1(2, 3);
    Matrix m2(3, 2);
    Matrix m3(2, 2);
 
    //std::cin >> m1;
    //std::cin >> m2;
    //std::cin >> m3;
 
    m1.randomFill(-5, 5);
    m2.randomFill(-5, 5);
    m3.randomFill(-5, 5);
 
    std::cout << m1 << '\n'
        << m2 << '\n'
        << m3 << '\n';
 
    Matrix result = m1 * m2 * m3;
 
    std::cout << result;
 
    system("PAUSE");
    return 0;
 
}