Designing a Matrix class - expression template issues

Question

I would like to design a Matrixclass 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 that will help . 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 spend quite some time over this. I would like to ask for some help; if you guys have played with expression templates, or find something fundamentally wrong with my code, it would be nice to correct it, and learn from it.

#include <array>
#include <iostream>

/**
 * Compile time fixed-size 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) { 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);
    }

    // friend class Exp<AddExp,T>;
   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);
    }

    // friend class Exp<Matrix,T>;
   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;

    return 0;
}

Answer 1

I see two problems: assignment vs initialization and const-correctness.

assignment vs initialization

In your main you have:

Matrix<int, 3, 3> C = A + B;

This requires a Matrix constructor which takes a const RHS& exp. Chenge it to:

Matrix<int, 3, 3> C;
C = A + B;

const-correctness

Your T operator()(int i, int j) should be const, since you don't want to change the expression or data when performing a sum.

Check it on Godbolt.

Answer 2

If you declare a variable and assign it in the same statement a constructor is actually used, i.e.

Matrix<int, 3, 3> C = A + B;

does not use the operator= of Matrix<int, 3, 3>, but looks for a constructor the expression A+B could be passed to.

You've only implemented the assignment operator though. You need to implement the a constructor taking a reference to a "matrix-like" type too for the compiler to accept this expression. You could simply do this by using the assignment operator in the constructor body.

Some additional notes:

• You use a fixed-size data structure to store the element data; furthermore you don't have any logic for resizing the matrix, which would only be possible in a very limited way, so I recommend getting rid of the member variables and simply defining static constexpr variables. Those variables could help you with avoiding matrices of incompatible dimensions together.

• operator= and operator() are incompatible; the indices in the implementation of operator= need to be reversed.

  Consider the case where a matrix has 100 rows and 1 column. The innermost loop would use (*this)(99, 0) which would return data[ROWS * 99 + 0], i.e. data[9801] which is clearly out of bounds.

• You don't check, if the matrix dimensions are compatible; furthermore you could determine the element type of AddExp automatically based on LHS and RHS.

• You could implement operator+ at namespace scope to avoid having to repeat the implementation for each expression types; you simply need to restrict the possible template parameters to something that is "matrix-like" using concepts (or SFINAE).

• If you take the parameter of operator= as reference to const, operator() for the parameter must have an implementation of operator() that works for const objects.

The following is an updated version taking these notes into account (C++ 20 required because of concept use):

#include <algorithm>
#include <array>
#include <iostream>
#include <type_traits>
#include <stdexcept>

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

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

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