How to combine lazy evaluation with auto in C++

Question

Some background on what I try to do: I am trying to implement a library doing quantum mechanics. As quantum mechanics is basically just linear algebra, I'm using the armadillo linear algebra library underneath. Armadillo uses lazy evaluation to do some smart tricks with matrices, which gives a pretty good abstraction from what is actually going on and looks close to matlab code.

I want to do something similar, but I also want to be able to use auto, which is not possible with armadillo (or eigen).

I have been looking around a little, and this answer contains what I think is the typical way of implementing this: https://stackoverflow.com/a/414260/6306265

The problem with this approach is that when you write

auto C = A+B;

you get a C that is a matrix_add, not a matrix. Even if matrix_add behaves similarly enough to matrix, the fact that matrix_add contains references to A and B makes it awkward to carry around. E.g.

auto A = matrix(2,2,{0,1,0,1});
auto B = matrix(2,2,{1,0,1,0});
auto C = A+B;
C.printmatrix(); // 1,1 ; 1,1

but

auto A = matrix(2,2,{0,1,0,1});
auto B = matrix(2,2,{1,0,1,0});
auto C = A+B;
A(0,0) = 1;
C.printmatrix(); // 2,1 ; 1,1

which is counter-intuitive. As mathematically intuitive behaviour is what I want to achieve, that is a problem.

Even worse is when I do

auto sumMatrices(const matrix& A, const matrix& B)
{
    return A+B;
}

which returns a matrix_add with references to local memory.

I would really like to be able to have the nice, overloaded behaviour but also be able to use auto. My idea was to make a wrapper that can hold either a reference or an instance:

template<class T>
class maybe_reference
{
public:
    maybe_reference(const T& t):
        ptr_(std::make_unique<T>(t)),
        t_(*ptr_)
    {}
    maybe_reference(std::reference_wrapper<const T> t): 
        t_(t.get())
    {}

    const T& get(){return t_;}
private:
    unique_ptr<T> ptr_;
    const T& t_;
}

It may not be implemented exactly this way, but the general idea is to have two constructors that can be clearly distinguished to ensure that get() returns either the referenced object or the one in the unique_ptr.

Modified matrix_add:

class matrix_add {
public:
    friend matrix_add operator+(const matrix& A, const matrix& B);

    matrix_add(matrix_add&& other): A_(other.A_.get()), B_(other.B_.get()){}
private:
    matrix_add(const matrix& A, const matrix& B): A_(std::ref(A)), B_(std::ref(B)){}

    maybe_reference<matrix> A_;
    maybe_reference<matrix> B_;
};

I have left out all the parts that make matrix_add behave like a matrix. The idea is to have the object refer to the outside objects A&B as long as it was constructed with A+B, but when it is move-constructed, it would own copies.

My question is basically: does this work?

I have been thinking that the move-constructor may be elided in some or all cases, which might be devastating.

Also, is there an alternative to achieve the same thing? I have been looking, but it seems that for linear algebra at least its either lazy or auto.

EDIT: Thanks to being reminded of the term "expression templates", my google search was a lot more fruitful. I found this reddit-post: https://www.reddit.com/r/cpp/comments/4puabu/news_about_operator_auto/
and the referenced papers, which allow specification of "casts" to auto. That would be the feature that really would make all of this work.

Answer 1

I think, your basic problem is, that lazy evaluation does not mix well with changing state. I see two possible routes out of this:

1. Make your matrices immutable. If you "modify" a matrix, you actually create a copy with the incorporated change, the original remains intact. This works well semantically (any math works exactly as you expect it to do), however it may incur an intolerable runtime overhead if you are setting your matrices value by value.

   This allows your implementation of matrix_add to silently replace itself with a matrix object when it is evaluated, ensuring that each evaluation is only performed at most once.

2. Make your functions explicit. Don't create matrix_add objects that act as if they were matrices themselves, but create matrix_function objects that operate on some input matrices to yield some result. This allows you to explicitly perform the evaluation where you see fit, and to reuse the functions that you define. However, this approach will lead to a lot of additional code complexity.

I don't think it's a good idea to try to work around this problem by introducing implicit points of forced evaluation: You'll loose large parts of what can be achieved by lazy evaluation, so why bother in the first place? Just my two cents.

Answer 2

You could write a template function evaluate which by default is a NOP, and then overload as necessary.

#include <utility>
#include <type_traits>

struct matrix {};
struct matrix_add {
  matrix operator()() const;
};

matrix_add operator + (matrix const& a, matrix const& b);


template<class T> decltype(auto) evaluate(T&& val) { return std::forward<T>(val); }
matrix evaluate(matrix_add const& lazy) { return lazy(); }
matrix evaluate(matrix_add & lazy) { return lazy(); }
matrix evaluate(matrix_add && lazy) { return lazy(); }



int main()
{
  auto a = matrix();
  auto b = matrix();

  auto c = evaluate(a + b);
  auto d = evaluate(1 + 2);

  static_assert(std::is_same<decltype(c), matrix>::value, "");
  static_assert(std::is_same<decltype(d), int>::value, "");

}

Answer 3

I will define a new operator: eager_eval, like this:

namespace lazy {
  template<class T>
  void eager_eval(T const volatile&)=delete;

  template<class T>
  struct expression {
    template<class D,
      std::enable_if_t<std::is_base_of<expression, std::decay_t<D>>{}, int> =0
    >
    friend T eager_eval( D&& d ) { return std::forward<D>(d); }
  };
}

Whenever you want something to be evaluatable in an eager manner, define eager_eval in its namespace, or derive it from lazy::lazy_expression<target_type>.

So we modify your matrix_add to (A) derive from it with the lazy-produced type you want, and (B) have an operator matrix:

struct matrix_add:
  lazy::expression<matrix>
{
  matrix_add(matrix const& a, matrix const& b) : a(a), b(b) { }

  operator matrix() && { // rvalue ref qualified as it should be.
    matrix result;
    // Do the addition.
    return result;
  }
private:
  matrix const& a, b;
};

and now, anyone can do:

auto e = eager_eval( a+b );

and ADL finds the right type to eager evaluate the lazy expression to.

live example.

You could, optionally, implement a default eager_eval that returns its argument:

  template<class T, class...Ts>
  T eager_eval(T&& t, Ts&&...)  { return std::forward<T>(t); }

then

using lazy::eager_eval;
auto x = eager_eval( 1+2 );

lets you be agnostic to the type you pass to eager_eval; if it is a type that is aware of being lazy via an eager_eval overload, it converts, and if not it does not convert.

The pack in lazy::eager_eval above is to ensure that it has the lowest priority as an overload.

Answer 4

with c++17 class template argument deduction, you may write

struct matrix_expr_foo {};
struct matrix_expr_bar {};

template< typename L, typename R >
struct matrix_add {
    // ...
};

matrix_add<matrix_expr_foo,matrix_expr_bar> operator + (matrix_expr_foo const& a, matrix_expr_bar const& b);

template< typename T >
struct expr {
    expr( T const& expr ){
        // evaluate expr ( to be stored in an appropriate member )
    }
    // ...
};

int main()
{
  auto a = matrix_expr_foo();
  auto b = matrix_expr_bar();
  expr c = a + b;

  /* different naming ?
  auto_ c = a + b;
  ...
  */
}

where expr is meant to act as an auto for expression templates ...