This is a follow-on question to my earlier post on writing a matrix-algebra header. I am using expression templates to avoid the creation of temporaries.
Lazy-evaluation is great for matrix operations performed coefficient-wise such as +, -, +=, -=, ==, scalar-multiplication etc. But, I learnt from the Eigen website, that it's not efficient for matrix products.
My code eagerly evaluates any matrix products. The trouble is, while I can write expressions such as A * B, A * (B+C) and (A+B)*C, I am unable to write (A+B)*(C+D), that is where both the left- and right- operands are AddExp objects.
The compiler fails to find the right signature, even though I have add the relevant method definition.
Could someone take a look at my code, and help me - what definition should I add, to write expressions of the type (A+B)*(C+D)?
Here is the code:
Lazy evaluation using expression templates
#include <array>
#include <iostream>
#include <cassert>
/**
* Compile time fixed-size matrix.
*/
template <typename LHS, typename RHS, typename T>
class SubExp;
template <typename T, int ROWS, int COLS>
class Matrix;
template <typename LHS, typename RHS, typename T>
class AddExp {
public:
AddExp(LHS& lhsExp, RHS& rhsExp) : m_lhsExp{ lhsExp }, m_rhsExp{ rhsExp } {}
T operator()(int i, int j) const { return m_lhsExp(i, j) + m_rhsExp(i, j); }
template <typename R>
AddExp<AddExp<LHS, RHS, T>, R, T> operator+(R& exp) {
return AddExp<AddExp<LHS, RHS, T>, R, T>(*this, exp);
}
template <typename R>
SubExp<AddExp<LHS, RHS, T>, R, T> operator-(R& exp) {
return SubExp<AddExp<LHS, RHS, T>, R, T>(*this, exp);
}
template <int ROWS, int COLS>
Matrix<T, ROWS, COLS> operator*(Matrix<T, ROWS, COLS>& exp) {
Matrix<T, ROWS, COLS> result;
assert(result.getRows() == exp.getCols());
assert(result.getCols() == exp.getRows());
int kMax = exp.getRows();
for (int i{}; i < ROWS; ++i)
{
for (int j{}; j < COLS; ++j)
{
for (int k{}; k < kMax; ++k)
{
result(i, j) += (*this)(i, k) * exp(k, j);
}
//std::cout << "result(" << i << "," << j << ")" << result(i, j) << "\n";
}
}
return result;
}
template <int ROWS, int COLS, typename L, typename R>
Matrix<T, ROWS, COLS> operator*(const AddExp<L, R, T>& exp2) {
Matrix<T, ROWS, COLS> result;
assert(getRows() == exp2.getCols());
assert(getCols() == exp2.getRows());
int kMax = exp2.getRows();
for (int i{}; i < ROWS; ++i)
{
for (int j{}; j < COLS; ++j)
{
for (int k{}; k < kMax; ++k)
{
result(i, j) += (*this)(i, k) * exp2(k, j);
}
//std::cout << "result(" << i << "," << j << ")" << result(i, j) << "\n";
}
}
return result;
}
LHS& getLHS() { return m_lhsExp; }
RHS& getRHS() { return m_rhsExp; }
int getRows() const { return m_lhsExp.getRows(); }
int getCols() const { return m_lhsExp.getCols(); }
private:
LHS& m_lhsExp;
RHS& m_rhsExp;
};
template <typename LHS, typename RHS, typename T>
class SubExp {
public:
SubExp(LHS& lhsExp, RHS& rhsExp) : m_lhsExp{ lhsExp }, m_rhsExp{ rhsExp } {}
T operator()(int i, int j) const { return m_lhsExp(i, j) - m_rhsExp(i, j); }
template <typename R>
SubExp<SubExp<LHS, RHS, T>, R, T> operator-(R& exp) {
return SubExp<SubExp<LHS, RHS, T>, R, T>(*this, exp);
}
template <typename R>
AddExp<SubExp<LHS, RHS, T>, R, T> operator+(R& exp) {
return AddExp<SubExp<LHS, RHS, T>, R, T>(*this, exp);
}
LHS& getLHS() { return m_lhsExp; }
RHS& getRHS() { return m_rhsExp; }
int getRows() const { return m_lhsExp.getRows(); }
int getCols() const { return m_lhsExp.getCols(); }
private:
LHS& m_lhsExp;
RHS& m_rhsExp;
};
template <typename T = double, int ROWS = 3, int COLS = 3>
class Matrix {
public:
Matrix() : m_rows{ ROWS }, m_cols{ COLS }, data{} {}
Matrix(const Matrix& m) : m_rows{ m.m_rows }, m_cols{ m.m_cols }, data{ m.data } {}
template <typename RHS>
Matrix(const RHS& exp) : m_rows{ ROWS }, m_cols{ COLS } {
for (int i{}; i < ROWS; ++i) {
for (int j{}; j < COLS; ++j) {
(*this)(i, j) = exp(i, j);
}
}
}
Matrix(std::initializer_list<std::initializer_list<T>> lst)
: m_rows{ ROWS }, m_cols{ COLS }, data{} {
int i{ 0 };
for (auto& l : lst) {
int j{ 0 };
for (auto& v : l) {
data[m_rows * i + j] = v;
++j;
}
++i;
}
}
T& operator()(int i, int j) { return data[ROWS * i + j]; }
template <typename RHS>
AddExp<Matrix<T, ROWS, COLS>, RHS, T> operator+(RHS& exp) {
return AddExp<Matrix<T, ROWS, COLS>, RHS, T>(*this, exp);
}
template <typename RHS>
Matrix<T, ROWS, COLS>& operator=(const RHS& exp) {
for (int i{}; i < ROWS; ++i) {
for (int j{}; j < COLS; ++j) {
(*this)(i, j) = exp(i, j);
}
}
return (*this);
}
template <typename L, typename R>
Matrix<T, ROWS, COLS> operator*(const AddExp<L, R, T>& exp) {
Matrix<T, ROWS, COLS> result;
assert(result.m_rows == exp.getCols());
assert(result.m_cols == exp.getRows());
int kMax = exp.getRows();
for (int i{}; i < m_rows; ++i)
{
for (int j{}; j < m_cols; ++j)
{
for (int k{}; k < kMax; ++k)
{
result(i, j) += (*this)(i, k) * exp(k, j);
}
//std::cout << "result(" << i << "," << j << ")" << result(i, j) << "\n";
}
}
return result;
}
template <typename L, typename R>
Matrix<T, ROWS, COLS> operator*(const SubExp<L, R, T>& exp) {
Matrix<T, ROWS, COLS> result;
assert(result.m_rows == exp.getCols());
assert(result.m_cols == exp.getRows());
int kMax = exp.getRows();
for (int i{}; i < m_rows; ++i)
{
for (int j{}; j < m_cols; ++j)
{
for (int k{}; k < kMax; ++k)
{
result(i, j) += (*this)(i, k) * exp(k, j);
}
}
}
return result;
}
template <int P>
Matrix<T, ROWS, P> operator*(const Matrix<T, COLS, P>& exp) {
Matrix<T, ROWS, COLS> result;
assert(result.m_rows == exp.getCols());
assert(result.m_cols == exp.getRows());
int kMax = exp.getRows();
for (int i{}; i < m_rows; ++i)
{
for (int j{}; j < m_cols; ++j)
{
for (int k{}; k < kMax; ++k)
{
result(i, j) += (*this)(i, k) * exp(k, j);
}
//std::cout << "result(" << i << "," << j << ")" << result(i, j) << "\n";
}
}
return result;
}
auto getData() { return data; }
int getRows() const { return m_rows; }
int getCols() const { return m_cols; }
private:
std::array<T, ROWS* COLS> data;
int m_rows;
int m_cols;
};
template <typename T, int ROWS, int COLS>
std::ostream& operator<<(std::ostream& stream, Matrix<T, ROWS, COLS>& m)
{
stream << "Matrix<" << m.getRows() << "," << m.getCols() << ">" << "[ \n";
for (int i{ 0 }; i < ROWS; ++i)
{
for (int j{ 0 }; j < COLS; ++j)
{
stream << " " << m(i, j) << " ";
}
stream << "\n";
}
stream << " ]";
return stream;
}
template <typename LHS, typename RHS, typename T>
std::ostream& operator<<(std::ostream& stream, const AddExp<LHS, RHS, T>& m)
{
stream << "Matrix<" << m.getRows() << "," << m.getCols() << ">" << "[ \n";
for (int i{ 0 }; i < m.getRows(); ++i)
{
for (int j{ 0 }; j < m.getCols(); ++j)
{
stream << " " << m(i, j) << " ";
}
stream << "\n";
}
stream << " ]";
return stream;
}
int main() {
Matrix<int, 3, 3> A{ {1, 0, 0}, {0, 1, 0}, {0, 0, 1} };
Matrix<int, 3, 3> B{ {1, 0, 0}, {0, 1, 0}, {0, 0, 1} };
Matrix<int, 3, 3> C = A + B;
Matrix<int, 3, 3> D = C * (A + B);
// Matrix<int, 3, 3> D = (C+C)*(A+B); //does not compile
std::cout << D;
return 0;
}
Note: the following solution uses the implementation from my previous answer.
Determining the dimensions of a AddExpr is hard to do without some modifications resulting in problems with determining the dimensions of the multiplication expression. (The dimensions need to be a compile time constant.) The only way of determining the dimensions would be to recursively inspect the rhs and lhs template params of the expressions, but simply keeping the info as static constexpr member variables as suggested in my previous answer makes it much easier to do this.
As stated in the answer it would be best to implement the operators at namespace scope which allows you to implement them for all types of matrix-like operands only once:
#include <algorithm>
#include <array>
#include <iostream>
#include <type_traits>
#include <stdexcept>
// ##########################################################################
// begin content from answer https://stackoverflow.com/a/76022576/2991525 ###
// ##########################################################################
using IndexType = unsigned;
template<class T, IndexType rows, IndexType columns>
struct MatrixType
{
using value_type = T;
static constexpr IndexType Rows = rows;
static constexpr IndexType Columns = columns;
};
template<class T>
concept MatrixLike = requires(IndexType row, IndexType column, T const mat)
{
mat(row, column);
std::is_base_of_v<MatrixType<typename T::value_type, T::Rows, T::Columns>, std::remove_cvref_t<T>>;
{T::Rows} -> std::convertible_to<IndexType>;
{T::Columns} -> std::convertible_to<IndexType>;
};
template<class T, class TargetValueType, IndexType rows, IndexType columns>
concept SizedMatrixLike = MatrixLike<T> && std::is_convertible_v<typename T::value_type, TargetValueType> && (T::Rows == rows) && (T::Columns == columns);
template<class M1, class M2>
using ValueAddResult_t = decltype(std::declval<typename M1::value_type>() + std::declval<typename M2::value_type>());
/**
* Compile time fixed-size matrix.
*/
template <MatrixLike LHS, MatrixLike RHS>
class AddExp : public MatrixType<ValueAddResult_t<LHS, RHS>, LHS::Rows, LHS::Columns> {
public:
static_assert(LHS::Rows == RHS::Rows, "row dimension mismatch");
static_assert(LHS::Columns == RHS::Columns, "column dimension mismatch");
AddExp(LHS const& lhsExp, RHS const& rhsExp)
: m_lhsExp{ lhsExp }, m_rhsExp{ rhsExp }
{
}
decltype(auto) operator()(IndexType i, IndexType j) const
{
return m_lhsExp(i, j) + m_rhsExp(i, j);
}
// friend class Exp<AddExp,T>;
private:
LHS const& m_lhsExp;
RHS const& m_rhsExp;
};
template <typename T = double, IndexType ROWS = 3, IndexType COLS = 3>
class Matrix : public MatrixType<T, ROWS, COLS> {
public:
Matrix() : m_data{} {}
Matrix(std::initializer_list<std::initializer_list<T>> lst)
: m_data{}
{
if (lst.size() > ROWS)
{
throw std::runtime_error("too many rows provided");
}
IndexType i{ 0 };
for (auto& l : lst) {
if (l.size() > COLS)
{
throw std::runtime_error("one of the rows has to many values");
}
std::copy(l.begin(), l.end(), m_data.begin() + (ROWS * i));
++i;
}
}
// copying the array directly is probably cheaper than using the () operator for each element
Matrix(Matrix const&) = default;
Matrix& operator=(Matrix const&) = default;
template<SizedMatrixLike<T, ROWS, COLS> U>
Matrix(U const& other)
{
*this = other;
}
T& operator()(IndexType column, IndexType row)
{
return m_data[ROWS * column + row];
}
T const& operator()(IndexType column, IndexType row) const
{
return m_data[ROWS * column + row];
}
template <SizedMatrixLike<T, ROWS, COLS> RHS>
Matrix& operator=(const RHS& other) {
auto out = m_data.begin();
for (IndexType c = 0; c != COLS; ++c)
{
for (IndexType r = 0; r != ROWS; ++r)
{
*out = other(c, r);
++out;
}
}
return *this;
}
private:
std::array<T, ROWS * COLS> m_data;
};
template<MatrixLike T, MatrixLike U>
requires (T::Rows == U::Rows) && (T::Columns == U::Columns)
AddExp<T, U> operator+(T const& s1, U const& s2)
{
return AddExp(s1, s2);
}
// ##########################################################################
// end content from old answer ##############################################
// ##########################################################################
template<class M1, class M2>
using ValueMultiplyResult_t = decltype(std::declval<typename M1::value_type>() * std::declval<typename M2::value_type>());
template<MatrixLike F1, MatrixLike F2>
requires (F1::Columns == F2::Rows)
auto operator*(F1 const& f1, F2 const& f2)
{
using ResultType = Matrix<ValueMultiplyResult_t<F1, F2>, F1::Rows, F2::Columns>;
ResultType result;
for (IndexType column{}; column < ResultType::Columns; ++column)
{
for (IndexType row{}; row < ResultType::Rows; ++row)
{
auto& resultElement = result(column, row);
for (IndexType f1Column{}; f1Column < F1::Columns; ++f1Column)
{
resultElement += f1(f1Column, row) * f2(column, f1Column);
}
}
}
return result;
}
template <MatrixLike T>
std::ostream& operator<<(std::ostream& stream, T const& m)
{
stream << "Matrix<" << T::Rows << "," << T::Columns << ">" << "[ \n";
for (IndexType column{ 0 }; column < T::Columns; ++column)
{
for (IndexType row{ 0 }; row < T::Rows; ++row)
{
stream << " " << m(column, row) << " ";
}
stream << "\n";
}
stream << " ]";
return stream;
}
int main() {
Matrix<int, 3, 3> A{ {1, 0, 0}, {0, 1, 0}, {0, 0, 1} };
Matrix<int, 3, 3> B{ {1, 0, 0}, {0, 1, 0}, {0, 0, 1} };
Matrix<int, 3, 3> C = A + B;
Matrix<int, 3, 3> D = C * (A + B);
std::cout << D << '\n';
Matrix<int, 3, 3> E = (C+C)*(A+B);
std::cout << E << '\n';
}