Designing a matrix header library - using expression templates (ET)

[Quasar Chunawala]

Hi guys,

I would like to design a Matrix class that uses expression templates and avoids temporaries. I would like to start with something very simple, so I decided to just implement matrix addition, first.

I created a class AddExp to represent a compound expression and overloaded the () operator. And the copy-assignment operator Matrix& Matrix::operator=(AddExp& addExp) triggers the evaluation, by the recursively invoking () on each of the elements in this expression.

However, the compiler complains that it can't convert AddExp<Matrix<int,3,3>,Matrix<int,3,3>> to a Matrix<int,3,3>. It puzzles me, why the copy-assignment operator of the Matrix class can't do that.

Note that, if I replace Matrix C by auto C, the code will execute fine, but it's important to be able to write statements of the form Matrix C = <something>.

Compiler Explorer - C++ (x86-64 gcc 12.2)

I would like to ask for some help; if you guys played with expression templates, or find something fundamentally wrong with my code, it would be nice to correct it, and learn from it.

Code:
[CODE lang="cpp" title="Compile-time matrix using ET"]#include <array>
#include <iostream>


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) { 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);
}


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(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);
}

private:
std::array<T, ROWS * COLS> data;
int m_rows;
int m_cols;
};

int main() {
Matrix A{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};

Matrix B{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};

Matrix<int, 3, 3> C = A + B; //Does not compile

return 0;
}[/CODE]

 

[Quasar Chunawala]

There were two problems above:

• Initialization v/s Assignment. We need to add a copy-constructor in the Matrix class that accepts a const RHS& exp, because Matrix C = A + B calls the copy constructor when building the object C, and not the assignment operator operator=.
• const correctness. The overloaded method () should be made const.

The below code now builds and executes as expected.

[CODE lang="cpp" title="Lazy Evaluation of coefficient wise matrix operations"]
#include <array>
#include <cassert>
#include <iostream>
#include <format>

/**
* Compile time fixed-size matrix.
*/

template <typename LHS, typename RHS, typename T, int M, int N>
class AddExp;

template <typename LHS, typename RHS, typename T, int M, int N>
class SubExp;

template <typename T, int ROWS, int COLS>
class Matrix;

template <typename T, typename RHS, int ROWS, int COLS>
class ScalarMult
{
template<typename L, typename R>
ScalarMult(T k, AddExp<L, R, T, ROWS, COLS>& rhsExp) : scalarConstant{ k }, m_rhsExp{ rhsExp } {}

template<typename L, typename R>
ScalarMult(T k, SubExp<L, R, T, ROWS, COLS>& rhsExp) : scalarConstant{ k }, m_rhsExp{ rhsExp } {}

template<typename R>
ScalarMult(T k, ScalarMult<T, R, ROWS, COLS>& rhsExp) : scalarConstant{ k }, m_rhsExp{ rhsExp } {}

ScalarMult(T k, Matrix<T, ROWS, COLS>& rhsExp) : scalarConstant{ k }, m_rhsExp{ rhsExp } {}

T operator()(int i, int j) const { return scalarConstant * m_rhsExp(i, j); }

private:
T scalarConstant;
RHS& m_rhsExp;
};

template <typename LHS, typename RHS, typename T, int M, int N>
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, M, N>, R, T, M, N> operator+(R& exp) {
return AddExp<AddExp<LHS, RHS, T, M, N>, R, T, M, N>(*this, exp);
}

template <typename R>
SubExp<AddExp<LHS, RHS, T, M, N>, R, T, M, N> operator-(R& exp) {
return SubExp<AddExp<LHS, RHS, T, M, N>, R, T, M, N>(*this, exp);
}

template <int COLS>
Matrix<T, M, COLS> operator*(Matrix<T, N, COLS>& exp) {
Matrix<T, M, COLS> result;
if (getCols() != exp.getRows())
throw std::runtime_error("dimension mismatch!");


int kMax = exp.getRows();
for (int i{}; i < M; ++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 COLS, typename L, typename R>
Matrix<T, M, COLS> operator*(const AddExp<L, R, T, N, COLS>& exp2) {
Matrix<T, M, COLS> result;
if (getCols() != exp2.getRows())
throw std::runtime_error("dimension mismatch!");

int kMax = exp2.getRows();
for (int i{}; i < M; ++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, int M, int N>
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, M, N>, R, T, M, N> operator-(R& exp) {
return SubExp<SubExp<LHS, RHS, T, M, N>, RHS, T, M, N>(*this, exp);
}

template <typename R>
AddExp<SubExp<LHS, RHS, T, M, N>, R, T, M, N > operator+(R& exp) {
return AddExp<SubExp<LHS, RHS, T, M, N>, R, T, M, N>(*this, exp);
}

template <int COLS>
Matrix<T, M, COLS> operator*(Matrix<T, N, COLS>& exp) {
Matrix<T, M, COLS> result;
if (getCols() != exp.getRows())
throw std::runtime_error("dimension mismatch!");


int kMax = exp.getRows();
for (int i{}; i < M; ++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 COLS, typename L, typename R>
Matrix<T, M, COLS> operator*(const SubExp<L, R, T, N, COLS>& exp2) {
Matrix<T, M, COLS> result;
if (getCols() != exp2.getRows())
throw std::runtime_error("dimension mismatch!");

int kMax = exp2.getRows();
for (int i{}; i < M; ++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 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_cols * i + j] = v;
++j;
}
++i;
}
}

T& operator()(int i, int j) { return data[COLS * i + j]; }
T operator()(int i, int j) const { return data[COLS * i + j]; }

template <typename RHS>
AddExp<Matrix<T, ROWS, COLS>, RHS, T, ROWS, COLS> operator+(RHS& exp) {
return AddExp<Matrix<T, ROWS, COLS>, RHS, T, ROWS, COLS>(*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, int P>
Matrix<T, ROWS, P> operator*(const AddExp<L, R, T, COLS, P>& exp) {
Matrix<T, ROWS, P> result;
if (getCols() != exp.getRows())
throw std::runtime_error("dimensions mismatch");

int kMax = exp.getRows();
for (int i{}; i < ROWS; ++i) {
for (int j{}; j < P; ++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, int P>
Matrix<T, ROWS, COLS> operator*(const SubExp<L, R, T, COLS, P>& exp) {
Matrix<T, ROWS, P> result;
if (getCols() != exp.getRows())
throw std::runtime_error("dimensions mismatch");

int kMax = exp.getRows();
for (int i{}; i < m_rows; ++i) {
for (int j{}; j < P; ++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, P> result;
if (getCols() != exp.getRows())
throw std::runtime_error("dimensions mismatch");

int kMax = exp.getRows();
for (int i{}; i < m_rows; ++i) {
for (int j{}; j < P; ++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, int M, int N>
std::ostream& operator<<(std::ostream& stream, const AddExp<LHS, RHS, T, M, N>& 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;
}

using Vector1d = Matrix<double, 1, 1>;
using Vector2d = Matrix<double, 2, 1>;
using Vector3d = Matrix<double, 3, 1>;
using Vector4d = Matrix<double, 4, 1>;


using Matrix2d = Matrix<double, 2, 2>;
using Matrix3d = Matrix<double, 3, 3>;
using Matrix4d = Matrix<double, 4, 4>;

using Vector1f = Matrix<float, 1, 1>;
using Vector2f = Matrix<float, 2, 1>;
using Vector3f = Matrix<float, 3, 1>;
using Vector4f = Matrix<float, 4, 1>;

using Matrix2f = Matrix<float, 2, 2>;
using Matrix3f = Matrix<float, 3, 3>;
using Matrix4f = Matrix<float, 4, 4>;

using Matrix2i = Matrix<int, 2, 2>;
using Matrix3i = Matrix<int, 3, 3>;
using Matrix4i = Matrix<int, 4, 4>;


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, 1> x = { {1,2,3} };

Matrix<int, 3, 1> E = A * x;

Matrix<int, 3, 3> F = (C + C) * (A + B);
std::cout << F;

return 0;
}


[/CODE]

 

[Quasar Chunawala]

I added support for matrix products. Whenever a matrix product is encountered, this is evaluated eagerly, as it is more efficient.

Right now, there's a lot of code duplication - so I hope to refresh my inheritance and polymorphism skills, and create a BinaryExp base class, and have AddExp, SubExp and ScalarMult as its sub-types.

 

[Daniel Duffy]

Why, oh why?

 

[Quasar Chunawala]

Why, oh why?

Hi Daniel, do you mean inheritance with templates complicates design, and makes the code unreadable? Did you mean to suggest, that the design should be kept simple...

I just picked up a couple of C++ books from the library, and would like to apply those ideas, to some small mathematical-finance tasks. e.g. write my own AAD code.

 

[Daniel Duffy]

No, I mean

1. My feeling is your C++ depth is not sufficient.
2. All this stuff has already been done.
3. It's reinvention of the wheel.

 

[Quasar Chunawala]

No, I mean

1. My feeling is your C++ depth is not sufficient.
2. All this stuff has already been done.
3. It's reinvention of the wheel.

For sure. Out of college, I was only trained on C++ 98. I finished reviewing the first 12 topics from this book - Beginning C++ 20 from novice to professional. I am encountering several new ideas. I am not as proficient, as I'd like to be.

On the bright side, I do refer from time to time, to en.cppreference, or after finishing a small task, ask for code reviews on stackoverflow.

When it comes to the choice of tasks, my time would perhaps be better spent on working on some new problem/taking something forward.

But, working on a new task, and learning language features at the same time sounded to me like a double whammy. So, I'd like to atleast get my programming skills to a green-belt/blue-belt level.

 

[Daniel Duffy]

What about writing some ML stuff in C++20.

No one writes home-grown matrices any more.

get my programming skills to a green-belt/blue-belt level.
In judo, you do stuff at that level. No free lunch.

1st, check how sensei does it.

More C++ Idioms/Expression-template - Wikibooks, open books for an open world

en.wikibooks.org

 

[Daniel Duffy]

I added support for matrix products. Whenever a matrix product is encountered, this is evaluated eagerly, as it is more efficient.

Right now, there's a lot of code duplication - so I hope to refresh my inheritance and polymorphism skills, and create a BinaryExp base class, and have AddExp, SubExp and ScalarMult as its sub-types.

This approach is a bit outdated TBH ... the functional language Haskell can do it in its sleep and has been an influence on C++20.