#### Strassen Matrix Multiplication

a
N
T
O
o
A
and 7 more contributors
```/**
* @brief [Strassen's
* algorithm](https://en.wikipedia.org/wiki/Strassen_algorithm) is one of the
* methods for multiplying two matrices. It is one of the faster algorithms for
* larger matrices than naive multiplication method.
*
* It involves dividing each matrices into 4 blocks, given they are evenly
* divisible, and are combined with new defined matrices involving 7 matrix
* multiplications instead of eight, yielding O(n^2.8073) complexity.
*
* @author [AshishYUO](https://github.com/AshishYUO)
*/
#include <cassert>   /// For assert operation
#include <chrono>    /// For std::chrono; time measurement
#include <iostream>  /// For I/O operations
#include <tuple>     /// For std::tuple
#include <vector>    /// For creating dynamic arrays

/**
* @namespace divide_and_conquer
* @brief Divide and Conquer algorithms
*/
namespace divide_and_conquer {

/**
* @namespace strassens_multiplication
* @brief Namespace for performing strassen's multiplication
*/
namespace strassens_multiplication {

/// Complement of 0 is a max integer.
constexpr size_t MAX_SIZE = ~0ULL;
/**
* @brief Matrix class.
*/
template <typename T,
typename = typename std::enable_if<
std::is_integral<T>::value || std::is_floating_point<T>::value,
bool>::type>
class Matrix {
std::vector<std::vector<T>> _mat;

public:
/**
* @brief Constructor
* @tparam Integer ensuring integers are being evaluated and not other
* data types.
* @param size denoting the size of Matrix as size x size
*/
template <typename Integer,
typename = typename std::enable_if<
std::is_integral<Integer>::value, Integer>::type>
explicit Matrix(const Integer size) {
for (size_t i = 0; i < size; ++i) {
_mat.emplace_back(std::vector<T>(size, 0));
}
}

/**
* @brief Constructor
* @tparam Integer ensuring integers are being evaluated and not other
* data types.
* @param rows denoting the total rows of Matrix
* @param cols denoting the total elements in each row of Matrix
*/
template <typename Integer,
typename = typename std::enable_if<
std::is_integral<Integer>::value, Integer>::type>
Matrix(const Integer rows, const Integer cols) {
for (size_t i = 0; i < rows; ++i) {
_mat.emplace_back(std::vector<T>(cols, 0));
}
}

/**
* @brief Get the matrix shape
* @returns pair of integer denoting total rows and columns
*/
inline std::pair<size_t, size_t> size() const {
return {_mat.size(), _mat[0].size()};
}

/**
* @brief returns the address of the element at ith place
* (here ith row of the matrix)
* @tparam Integer any valid integer
* @param index index which is requested
* @returns the address of the element (here ith row or array)
*/
template <typename Integer,
typename = typename std::enable_if<
std::is_integral<Integer>::value, Integer>::type>
inline std::vector<T> &operator[](const Integer index) {
return _mat[index];
}

/**
* @brief Creates a new matrix and returns a part of it.
* @param row_start start of the row
* @param row_end end of the row
* @param col_start start of the col
* @param col_end end of the column
* @returns A slice of (row_end - row_start) x (col_end - col_start) size of
* array starting from row_start row and col_start column
*/
Matrix slice(const size_t row_start, const size_t row_end = MAX_SIZE,
const size_t col_start = MAX_SIZE,
const size_t col_end = MAX_SIZE) const {
const size_t h_size =
(row_end != MAX_SIZE ? row_end : _mat.size()) - row_start;
const size_t v_size = (col_end != MAX_SIZE ? col_end : _mat[0].size()) -
(col_start != MAX_SIZE ? col_start : 0);
Matrix result = Matrix<T>(h_size, v_size);

const size_t v_start = (col_start != MAX_SIZE ? col_start : 0);
for (size_t i = 0; i < h_size; ++i) {
for (size_t j = 0; j < v_size; ++j) {
result._mat[i][j] = _mat[i + row_start][j + v_start];
}
}
return result;
}

/**
* @brief Horizontally stack the matrix (one after the other)
* @tparam Number any type of number
* @param other the other matrix: note that this array is not modified
* @returns void, but modifies the current array
*/
template <typename Number, typename = typename std::enable_if<
std::is_integral<Number>::value ||
std::is_floating_point<Number>::value,
Number>::type>
void h_stack(const Matrix<Number> &other) {
assert(_mat.size() == other._mat.size());
for (size_t i = 0; i < other._mat.size(); ++i) {
for (size_t j = 0; j < other._mat[i].size(); ++j) {
_mat[i].push_back(other._mat[i][j]);
}
}
}

/**
* @brief Horizontally stack the matrix (current matrix above the other)
* @tparam Number any type of number (Integer or floating point)
* @param other the other matrix: note that this array is not modified
* @returns void, but modifies the current array
*/
template <typename Number, typename = typename std::enable_if<
std::is_integral<Number>::value ||
std::is_floating_point<Number>::value,
Number>::type>
void v_stack(const Matrix<Number> &other) {
assert(_mat[0].size() == other._mat[0].size());
for (size_t i = 0; i < other._mat.size(); ++i) {
_mat.emplace_back(std::vector<T>(other._mat[i].size()));
for (size_t j = 0; j < other._mat[i].size(); ++j) {
_mat.back()[j] = other._mat[i][j];
}
}
}

/**
* @brief Add two matrices and returns a new matrix
* @tparam Number any real value to add
* @param other Other matrix to add to this
* @returns new matrix
*/
template <typename Number, typename = typename std::enable_if<
std::is_integral<Number>::value ||
std::is_floating_point<Number>::value,
bool>::type>
Matrix operator+(const Matrix<Number> &other) const {
assert(this->size() == other.size());
Matrix C = Matrix<Number>(_mat.size(), _mat[0].size());
for (size_t i = 0; i < _mat.size(); ++i) {
for (size_t j = 0; j < _mat[i].size(); ++j) {
C._mat[i][j] = _mat[i][j] + other._mat[i][j];
}
}
return C;
}

/**
* @brief Add another matrices to current matrix
* @tparam Number any real value to add
* @param other Other matrix to add to this
* @returns reference of current matrix
*/
template <typename Number, typename = typename std::enable_if<
std::is_integral<Number>::value ||
std::is_floating_point<Number>::value,
bool>::type>
Matrix &operator+=(const Matrix<Number> &other) const {
assert(this->size() == other.size());
for (size_t i = 0; i < _mat.size(); ++i) {
for (size_t j = 0; j < _mat[i].size(); ++j) {
_mat[i][j] += other._mat[i][j];
}
}
return this;
}

/**
* @brief Subtract two matrices and returns a new matrix
* @tparam Number any real value to multiply
* @param other Other matrix to subtract to this
* @returns new matrix
*/
template <typename Number, typename = typename std::enable_if<
std::is_integral<Number>::value ||
std::is_floating_point<Number>::value,
bool>::type>
Matrix operator-(const Matrix<Number> &other) const {
assert(this->size() == other.size());
Matrix C = Matrix<Number>(_mat.size(), _mat[0].size());
for (size_t i = 0; i < _mat.size(); ++i) {
for (size_t j = 0; j < _mat[i].size(); ++j) {
C._mat[i][j] = _mat[i][j] - other._mat[i][j];
}
}
return C;
}

/**
* @brief Subtract another matrices to current matrix
* @tparam Number any real value to Subtract
* @param other Other matrix to Subtract to this
* @returns reference of current matrix
*/
template <typename Number, typename = typename std::enable_if<
std::is_integral<Number>::value ||
std::is_floating_point<Number>::value,
bool>::type>
Matrix &operator-=(const Matrix<Number> &other) const {
assert(this->size() == other.size());
for (size_t i = 0; i < _mat.size(); ++i) {
for (size_t j = 0; j < _mat[i].size(); ++j) {
_mat[i][j] -= other._mat[i][j];
}
}
return this;
}

/**
* @brief Multiply two matrices and returns a new matrix
* @tparam Number any real value to multiply
* @param other Other matrix to multiply to this
* @returns new matrix
*/
template <typename Number, typename = typename std::enable_if<
std::is_integral<Number>::value ||
std::is_floating_point<Number>::value,
bool>::type>
inline Matrix operator*(const Matrix<Number> &other) const {
assert(_mat[0].size() == other._mat.size());
auto size = this->size();
const size_t row = size.first, col = size.second;
// Main condition for applying strassen's method:
// 1: matrix should be a square matrix
// 2: matrix should be of even size (mat.size() % 2 == 0)
return (row == col && (row & 1) == 0)
? this->strassens_multiplication(other)
: this->naive_multiplication(other);
}

/**
* @brief Multiply matrix with a number and returns a new matrix
* @tparam Number any real value to multiply
* @param other Other real number to multiply to current matrix
* @returns new matrix
*/
template <typename Number, typename = typename std::enable_if<
std::is_integral<Number>::value ||
std::is_floating_point<Number>::value,
bool>::type>
inline Matrix operator*(const Number other) const {
Matrix C = Matrix<Number>(_mat.size(), _mat[0].size());
for (size_t i = 0; i < _mat.size(); ++i) {
for (size_t j = 0; j < _mat[i].size(); ++j) {
C._mat[i][j] = _mat[i][j] * other;
}
}
return C;
}

/**
* @brief Multiply a number to current matrix
* @tparam Number any real value to multiply
* @param other Other matrix to multiply to this
* @returns reference of current matrix
*/
template <typename Number, typename = typename std::enable_if<
std::is_integral<Number>::value ||
std::is_floating_point<Number>::value,
bool>::type>
Matrix &operator*=(const Number other) const {
for (size_t i = 0; i < _mat.size(); ++i) {
for (size_t j = 0; j < _mat[i].size(); ++j) {
_mat[i][j] *= other;
}
}
return this;
}

/**
* @brief Naive multiplication performed on this
* @tparam Number any real value to multiply
* @param other Other matrix to multiply to this
* @returns new matrix
*/
template <typename Number, typename = typename std::enable_if<
std::is_integral<Number>::value ||
std::is_floating_point<Number>::value,
bool>::type>
Matrix naive_multiplication(const Matrix<Number> &other) const {
Matrix C = Matrix<Number>(_mat.size(), other._mat[0].size());

for (size_t i = 0; i < _mat.size(); ++i) {
for (size_t k = 0; k < _mat[0].size(); ++k) {
for (size_t j = 0; j < other._mat[0].size(); ++j) {
C._mat[i][j] += _mat[i][k] * other._mat[k][j];
}
}
}
return C;
}

/**
* @brief Strassens method of multiplying two matrices
* References: https://en.wikipedia.org/wiki/Strassen_algorithm
* @tparam Number any real value to multiply
* @param other Other matrix to multiply to this
* @returns new matrix
*/
template <typename Number, typename = typename std::enable_if<
std::is_integral<Number>::value ||
std::is_floating_point<Number>::value,
bool>::type>
Matrix strassens_multiplication(const Matrix<Number> &other) const {
const size_t size = _mat.size();
// Base case: when a matrix is small enough for faster naive
// multiplication, or the matrix is of odd size, then go with the naive
// multiplication route;
// else; go with the strassen's method.
if (size <= 64ULL || (size & 1ULL)) {
return this->naive_multiplication(other);
} else {
const Matrix<Number>
A = this->slice(0ULL, size >> 1, 0ULL, size >> 1),
B = this->slice(0ULL, size >> 1, size >> 1, size),
C = this->slice(size >> 1, size, 0ULL, size >> 1),
D = this->slice(size >> 1, size, size >> 1, size),
E = other.slice(0ULL, size >> 1, 0ULL, size >> 1),
F = other.slice(0ULL, size >> 1, size >> 1, size),
G = other.slice(size >> 1, size, 0ULL, size >> 1),
H = other.slice(size >> 1, size, size >> 1, size);

Matrix P1 = A.strassens_multiplication(F - H);
Matrix P2 = (A + B).strassens_multiplication(H);
Matrix P3 = (C + D).strassens_multiplication(E);
Matrix P4 = D.strassens_multiplication(G - E);
Matrix P5 = (A + D).strassens_multiplication(E + H);
Matrix P6 = (B - D).strassens_multiplication(G + H);
Matrix P7 = (A - C).strassens_multiplication(E + F);

// Building final matrix C11 would be
//     [      |      ]
//     [ C11  |  C12 ]
// C = [ ____ | ____ ]
//     [      |      ]
//     [ C21  |  C22 ]
//     [      |      ]

Matrix C11 = P5 + P4 - P2 + P6;
Matrix C12 = P1 + P2;
Matrix C21 = P3 + P4;
Matrix C22 = P1 + P5 - P3 - P7;

C21.h_stack(C22);
C11.h_stack(C12);
C11.v_stack(C21);

return C11;
}
}

/**
* @brief Compares two matrices if each of them are equal or not
* @param other other matrix to compare
* @returns whether they are equal or not
*/
bool operator==(const Matrix<T> &other) const {
if (_mat.size() != other._mat.size() ||
_mat[0].size() != other._mat[0].size()) {
return false;
}
for (size_t i = 0; i < _mat.size(); ++i) {
for (size_t j = 0; j < _mat[i].size(); ++j) {
if (_mat[i][j] != other._mat[i][j]) {
return false;
}
}
}
return true;
}

friend std::ostream &operator<<(std::ostream &out, const Matrix<T> &mat) {
for (auto &row : mat._mat) {
for (auto &elem : row) {
out << elem << " ";
}
out << "\n";
}
return out << "\n";
}
};

}  // namespace strassens_multiplication

}  // namespace divide_and_conquer

/**
* @brief Self-test implementations
* @returns void
*/
static void test() {
const size_t s = 512;
auto matrix_demo =
divide_and_conquer::strassens_multiplication::Matrix<size_t>(s, s);

for (size_t i = 0; i < s; ++i) {
for (size_t j = 0; j < s; ++j) {
matrix_demo[i][j] = i + j;
}
}

auto matrix_demo2 =
divide_and_conquer::strassens_multiplication::Matrix<size_t>(s, s);
for (size_t i = 0; i < s; ++i) {
for (size_t j = 0; j < s; ++j) {
matrix_demo2[i][j] = 2 + i + j;
}
}

auto start = std::chrono::system_clock::now();
auto Mat3 = matrix_demo2 * matrix_demo;
auto end = std::chrono::system_clock::now();

std::chrono::duration<double> time = (end - start);
std::cout << "Strassen time: " << time.count() << "s" << std::endl;

start = std::chrono::system_clock::now();
auto conf = matrix_demo2.naive_multiplication(matrix_demo);
end = std::chrono::system_clock::now();

time = end - start;
std::cout << "Normal time: " << time.count() << "s" << std::endl;

// std::cout << Mat3 << conf << std::endl;
assert(Mat3 == conf);
}

/**
* @brief main function
* @returns 0 on exit
*/
int main() {
test();  // run self-test implementation
return 0;
}
```