C++ Expression Templates

Can expression templates using references to temporaries be re-useable?

The comment above has a very effective way to check the problem with the dangling reference. Note that if you try to print the values from the main function in your example the program will still work because the object that will have the dangling reference bound to it will be created also on the stack space of main. I tried to move the code which is assigned to expr inside a function and the program crashed as expected (the temporary object will be in another stack frame):

auto makeExpr1(Vec const& v0, Vec const& v1, Vec const& v2) {
    return v0 + v1 + v2;
}
// ... in main:
auto expr = makeExpr1(v0, v1, v2);

The problem you highlighted here appears in the cases of creating an expression that can be lazily evaluated in languages like C++. A somehow similar situation can occur in the context of range expressions (C++20 ranges).
Below is my quick attempt to fix that code and make it work with lvalues and rvalues added with the operator + (I apologise for the ugly parts and possible mistakes). This will store copy of their operands only when they are going to be out of scope and will result in dangling references in the old code.

Regarding re-usability: as long as you define a type for every operation and a corresponding operator '?' function ('?' being the simbol of the operation) this approch should give you a starting point for any binary operation on such a vector.

#include <cassert>
#include <vector>
#include <utility>
#include <iostream>

/*
 * Passes lvalues and stores rvalues
 */
template <typename T> class Wrapper;

template <typename T> class Wrapper<T&> {
    private:
        T& ref;

    public:
        Wrapper(T& ref) : ref(ref) {}
        T& get() { return ref; }
        const T& get() const { return ref; }
};

template <typename T> class Wrapper<T&&> {
    private:
        T value;

    public:
        Wrapper(T&& ref) : value(std::move(ref)) {}
        T& get() { return value; }
        const T& get() const { return value; }
};

template <typename E>
class VecExpression {
  public:
    double operator[](size_t i) const 
    {
        // Delegation to the actual expression type. This avoids dynamic polymorphism (a.k.a. virtual functions in C++)
        return static_cast<E const&>(*this)[i];
    }
    size_t size()               const { return static_cast<E const&>(*this).size(); }
};

/*
 * Forwards the reference and const qualifiers
 *  of the expression type to the expression itself
 */
template <typename E> constexpr E& forwardRef(VecExpression<E>& ve) {
    return static_cast<E&>(ve);
}

template <typename E> constexpr const E& forwardRef(const VecExpression<E>& ve) {
    return static_cast<const E&>(ve);
}

template <typename E> constexpr E&& forwardRef(VecExpression<E>&& ve) {
    return static_cast<E&&>(ve);
}

class Vec : public VecExpression<Vec> {
    std::vector<double> elems;

  public:
    double operator[](size_t i) const { return elems[i]; }
    double &operator[](size_t i)      { return elems[i]; }
    size_t size() const               { return elems.size(); }

    Vec(size_t n) : elems(n) {}

    // construct vector using initializer list 
    Vec(std::initializer_list<double> init) : elems(init) {}

    // A Vec can be constructed from any VecExpression, forcing its evaluation.
    template <typename E>
    Vec(VecExpression<E> const& expr) : elems(expr.size()) {
        std::cout << "Expr ctor\n"; // Very quick test
        for (size_t i = 0; i != expr.size(); ++i) {
            elems[i] = expr[i];
        }
    }

    // Move ctor added for checking
    Vec(Vec&& vec) : elems(std::move(vec.elems)) {
        std::cout << "Move ctor\n";  // Very quick test
    }
};

/*
 * Now VecSum is a sum between possibly const - qualified
 *  and referenced expression types
 */
template <typename E1, typename E2>
class VecSum : public VecExpression<VecSum<E1, E2>> {

    Wrapper<E1> _u;
    Wrapper<E2> _v;

public:

    VecSum(E1 u, E2 v) : _u(static_cast<E1>(u)), _v(static_cast<E2>(v)) {
        assert(_u.get().size() == _v.get().size());
    }

    double operator[](size_t i) const { return _u.get()[i] + _v.get()[i]; }
    size_t size()               const { return _v.get().size(); }
};

/*
 * Used to create a VecSum by capturing also the reference kind
 *  of the arguments (will be used by the Wrapper inside VecSum)
 */
template <typename E1, typename E2>
auto makeVecSum(E1&& e1, E2&& e2) {
    return VecSum<E1&&, E2&&>(std::forward<E1>(e1), std::forward<E2>(e2));
}

/*
 * Now the operator+ takes the vector expressions by universal references
 */
template <typename VE1, typename VE2>
auto operator+(VE1&& ve1, VE2&& ve2) {
   return makeVecSum(forwardRef(std::forward<VE1>(ve1)), forwardRef(std::forward<VE2>(ve2)));
}

// Now this will work
auto makeExpr1(Vec const& v0, Vec const& v1, Vec const& v2) {
    return v0 + v1 + v2;
}

// This will also work - the rvalue is stored in the
//  expression itself and both will have the same lifetime
auto makeExpr2(Vec const& v0, Vec const& v1) {
    return v0 + v1 + Vec({1.0, 1.0, 1.0, 1.0});
}

int main() {

    Vec v0 = {23.4,12.5,144.56,90.56};
    Vec v1 = {67.12,34.8,90.34,89.30};
    Vec v2 = {34.90,111.9,45.12,90.5};

    auto expr = makeExpr1(v0, v1, v2);
    Vec v1_ = expr;
    Vec v2_ = expr;
    auto expr_ = makeExpr2(v0, v1);

    for (size_t i = 0; i < v1_.size(); ++i)
        std::cout << v1_[i] << " ";
    std::cout << std::endl;

    for (size_t i = 0; i < v2_.size(); ++i)
        std::cout << v2_[i] << " ";
    std::cout << std::endl;

    for (size_t i = 0; i < expr.size(); ++i)
        std::cout << expr[i] << " ";
    std::cout << std::endl;

    for (size_t i = 0; i < expr_.size(); ++i)
        std::cout << expr_[i] << " ";
    std::cout << std::endl;
}

Expression templates and C++11

I was wondering whether lambdas together with move semantics or any other new feature can do as good as ETs. Any thoughts?

Quick Answer

Move semantics are not a total panacea on their own --techniques such as expression templates (ETs) are still needed in C++11 to eliminate overheads such as moving data around! So, to answer your question quickly before diving into the rest of my answer, move semantics, etc. doesn't completely replace ETs as my answer illustrates below.

Detailed Answer

ETs typically return proxy objects to defer evaluation until later, so there is no immediate apparent benefit of C++11 language features until the code that triggers the computation. That said, one would not want to write ET code, however, that triggers run-time code generation during the building of the expression tree with the proxies. Nicely, C++11's move semantics and perfect forwarding can help avoid such overheads should that otherwise occur. (Such would not have been possible in C++03.)

Essentially, when writing ETs one wants to exploit the language features in a way to generate optimal code once the member function(s) of the involved proxy objects are invoked. In C++11 this will include using perfect forwarding, move semantics over copying, etc. if such is actually still needed over and above what the compiler can already do. The name of the game is to minimize the run-time code generated and/or maximize the run-time speed and/or minimize the run-time overhead.

I wanted to actually try some ETs with C++11 features to see if I could elide ALL intermediate temporary instance types with a a = b + c + d; expression. (As this was just a fun break from my normal activities so I did not compare it to or write ET code purely using C++03. Also I did not worry about all aspects of code polishing that appears below.)

To start with, I did not use lambdas --as I preferred to use explicit types and functions-- so I won't argue for/against lambdas with respect to your question. My guess is that they would be similar to using functors and performing no better than the non-ET code below (i.e., moves would be required) --at least until compilers can automatically optimize lambdas using their own internal ETs for such. The code I wrote, however, exploits move semantics and perfect forwarding. Here's what I did starting with the results and then finally presenting the code.

I created a math_vector<N> class where N==3 and it defines an internal private instance of std::array<long double, N>. The members are a default constructor, copy and move constructors and assignments, an initializer list constructor, a destructor, a swap() member, operator [] to access elements of the vector and operator +=. Used without any expression templates, this code:

{
  cout << "CASE 1:\n";
  math_vector<3> a{1.0, 1.1, 1.2};
  math_vector<3> b{2.0, 2.1, 2.2};
  math_vector<3> c{3.0, 3.1, 3.2};
  math_vector<3> d{4.0, 4.1, 4.2};
  math_vector<3> result = a + b + c + d;
  cout << '[' << &result << "]: " << result << "\n";
}

outputs (when compiled with clang++ 3.1 or g++ 4.8 with -std=c++11 -O3):

CASE 1:
0x7fff8d6edf50: math_vector(initlist)
0x7fff8d6edef0: math_vector(initlist)
0x7fff8d6ede90: math_vector(initlist)
0x7fff8d6ede30: math_vector(initlist)
0x7fff8d6edd70: math_vector(copy: 0x7fff8d6edf50)
0x7fff8d6edda0: math_vector(move: 0x7fff8d6edd70)
0x7fff8d6eddd0: math_vector(move: 0x7fff8d6edda0)
0x7fff8d6edda0: ~math_vector()
0x7fff8d6edd70: ~math_vector()
[0x7fff8d6eddd0]: (10,10.4,10.8)
0x7fff8d6eddd0: ~math_vector()
0x7fff8d6ede30: ~math_vector()
0x7fff8d6ede90: ~math_vector()
0x7fff8d6edef0: ~math_vector()
0x7fff8d6edf50: ~math_vector()

i.e., the four explicit constructed instances using initializer lists (i.e., the initlist items), the result variable (i.e., 0x7fff8d6eddd0), and, also makes an additional three objects copying and moving.

To only focus on temporaries and moving, I created a second case that only creates result as a named variable --all others are rvalues:

{
  cout << "CASE 2:\n";
  math_vector<3> result =
    math_vector<3>{1.0, 1.1, 1.2} +
    math_vector<3>{2.0, 2.1, 2.2} +
    math_vector<3>{3.0, 3.1, 3.2} +
    math_vector<3>{4.0, 4.1, 4.2}
  ;
  cout << '[' << &result << "]: " << result << "\n";
}

which outputs this (again when ETs are NOT used):

CASE 2:
0x7fff8d6edcb0: math_vector(initlist)
0x7fff8d6edc50: math_vector(initlist)
0x7fff8d6edce0: math_vector(move: 0x7fff8d6edcb0)
0x7fff8d6edbf0: math_vector(initlist)
0x7fff8d6edd10: math_vector(move: 0x7fff8d6edce0)
0x7fff8d6edb90: math_vector(initlist)
0x7fff8d6edd40: math_vector(move: 0x7fff8d6edd10)
0x7fff8d6edb90: ~math_vector()
0x7fff8d6edd10: ~math_vector()
0x7fff8d6edbf0: ~math_vector()
0x7fff8d6edce0: ~math_vector()
0x7fff8d6edc50: ~math_vector()
0x7fff8d6edcb0: ~math_vector()
[0x7fff8d6edd40]: (10,10.4,10.8)
0x7fff8d6edd40: ~math_vector()

which is better: only extra move objects are created.

But I wanted better: I wanted zero extra temporaries and to have the code as if I hard-coded it with the one normal coding caveat: all explicitly instantiated types would still be created (i.e., the four initlist constructors and result). To accomplish this I then added expression template code as follows:

1. a proxy math_vector_expr<LeftExpr,BinaryOp,RightExpr> class was created to hold an expression not computed yet,
2. a proxy plus_op class was created to hold the addition operation,
3. a constructor was added to math_vector to accept a math_vector_expr object, and,
4. "starter" member functions were added to trigger the creation of the expression template.

The results using ETs are wonderful: no extra temporaries in either case! The previous two cases above now output:

CASE 1:
0x7fffe7180c60: math_vector(initlist)
0x7fffe7180c90: math_vector(initlist)
0x7fffe7180cc0: math_vector(initlist)
0x7fffe7180cf0: math_vector(initlist)
0x7fffe7180d20: math_vector(expr: 0x7fffe7180d90)
[0x7fffe7180d20]: (10,10.4,10.8)
0x7fffe7180d20: ~math_vector()
0x7fffe7180cf0: ~math_vector()
0x7fffe7180cc0: ~math_vector()
0x7fffe7180c90: ~math_vector()
0x7fffe7180c60: ~math_vector()

CASE 2:
0x7fffe7180dd0: math_vector(initlist)
0x7fffe7180e20: math_vector(initlist)
0x7fffe7180e70: math_vector(initlist)
0x7fffe7180eb0: math_vector(initlist)
0x7fffe7180d20: math_vector(expr: 0x7fffe7180dc0)
0x7fffe7180eb0: ~math_vector()
0x7fffe7180e70: ~math_vector()
0x7fffe7180e20: ~math_vector()
0x7fffe7180dd0: ~math_vector()
[0x7fffe7180d20]: (10,10.4,10.8)
0x7fffe7180d20: ~math_vector()

i.e., exactly 5 constructor calls and 5 destructor calls in each case. In fact, if you ask the compiler to generate the assembler code between the 4 initlist constructor calls and the outputting of result one gets this beautiful string of assembler code:

fldt    128(%rsp)
leaq    128(%rsp), %rdi
leaq    80(%rsp), %rbp
fldt    176(%rsp)
faddp   %st, %st(1)
fldt    224(%rsp)
faddp   %st, %st(1)
fldt    272(%rsp)
faddp   %st, %st(1)
fstpt   80(%rsp)
fldt    144(%rsp)
fldt    192(%rsp)
faddp   %st, %st(1)
fldt    240(%rsp)
faddp   %st, %st(1)
fldt    288(%rsp)
faddp   %st, %st(1)
fstpt   96(%rsp)
fldt    160(%rsp)
fldt    208(%rsp)
faddp   %st, %st(1)
fldt    256(%rsp)
faddp   %st, %st(1)
fldt    304(%rsp)
faddp   %st, %st(1)
fstpt   112(%rsp)

with g++ and clang++ outputs similar (even smaller) code. No function calls, etc. --just a bunch of adds which is EXACTLY what one wants!

The C++11 code to achieve this follows. Simply #define DONT_USE_EXPR_TEMPL to not use ETs or don't define it at all to use ETs.

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

//#define DONT_USE_EXPR_TEMPL

//===========================================================================

template <std::size_t N> class math_vector;

template <
  typename LeftExpr,
  typename BinaryOp,
  typename RightExpr
>
class math_vector_expr
{
  public:
    math_vector_expr() = delete;

    math_vector_expr(LeftExpr l, RightExpr r) : 
      l_(std::forward<LeftExpr>(l)), 
      r_(std::forward<RightExpr>(r))
    {
    }

    // Prohibit copying...
    math_vector_expr(math_vector_expr const&) = delete;
    math_vector_expr& operator =(math_vector_expr const&) = delete;

    // Allow moves...
    math_vector_expr(math_vector_expr&&) = default;
    math_vector_expr& operator =(math_vector_expr&&) = default;

    template <typename RE>
    auto operator +(RE&& re) const ->
      math_vector_expr<
        math_vector_expr<LeftExpr,BinaryOp,RightExpr> const&,
        BinaryOp,
        decltype(std::forward<RE>(re))
      >
    {
      return 
        math_vector_expr<
          math_vector_expr<LeftExpr,BinaryOp,RightExpr> const&,
          BinaryOp,
          decltype(std::forward<RE>(re))
        >(*this, std::forward<RE>(re))
      ;
    }

    auto le() -> 
      typename std::add_lvalue_reference<LeftExpr>::type
      { return l_; }

    auto le() const ->
      typename std::add_lvalue_reference<
        typename std::add_const<LeftExpr>::type
      >::type
      { return l_; }

    auto re() -> 
      typename std::add_lvalue_reference<RightExpr>::type
      { return r_; }

    auto re() const -> 
      typename std::add_lvalue_reference<
        typename std::add_const<RightExpr>::type
      >::type
      { return r_; }

    auto operator [](std::size_t index) const ->
      decltype(
        BinaryOp::apply(this->le()[index], this->re()[index])
      )
    {
      return BinaryOp::apply(le()[index], re()[index]);
    }

  private:
    LeftExpr l_;
    RightExpr r_;
};

//===========================================================================

template <typename T>
struct plus_op
{
  static T apply(T const& a, T const& b)
  {
    return a + b;
  }

  static T apply(T&& a, T const& b)
  {
    a += b;
    return std::move(a);
  }

  static T apply(T const& a, T&& b)
  {
    b += a;
    return std::move(b);
  }

  static T apply(T&& a, T&& b)
  {
    a += b;
    return std::move(a);
  }
};

//===========================================================================

template <std::size_t N>
class math_vector
{
  using impl_type = std::array<long double, N>;

  public:
    math_vector()
    {
      using namespace std;
      fill(begin(v_), end(v_), impl_type{});
      std::cout << this << ": math_vector()" << endl;
    }

    math_vector(math_vector const& mv) noexcept
    {
      using namespace std;
      copy(begin(mv.v_), end(mv.v_), begin(v_));
      std::cout << this << ": math_vector(copy: " << &mv << ")" << endl;
    }

    math_vector(math_vector&& mv) noexcept
    {
      using namespace std;
      move(begin(mv.v_), end(mv.v_), begin(v_));
      std::cout << this << ": math_vector(move: " << &mv << ")" << endl;
    }

    math_vector(std::initializer_list<typename impl_type::value_type> l)
    {
      using namespace std;
      copy(begin(l), end(l), begin(v_));
      std::cout << this << ": math_vector(initlist)" << endl;
    }

    math_vector& operator =(math_vector const& mv) noexcept
    {
      using namespace std;
      copy(begin(mv.v_), end(mv.v_), begin(v_));
      std::cout << this << ": math_vector op =(copy: " << &mv << ")" << endl;
      return *this;
    }

    math_vector& operator =(math_vector&& mv) noexcept
    {
      using namespace std;
      move(begin(mv.v_), end(mv.v_), begin(v_));
      std::cout << this << ": math_vector op =(move: " << &mv << ")" << endl;
      return *this;
    }

    ~math_vector()
    {
      using namespace std;
      std::cout << this << ": ~math_vector()" << endl;
    }

    void swap(math_vector& mv)
    {
      using namespace std;
      for (std::size_t i = 0; i<N; ++i)
        swap(v_[i], mv[i]);
    }

    auto operator [](std::size_t index) const
      -> typename impl_type::value_type const&
    {
      return v_[index];
    }

    auto operator [](std::size_t index)
      -> typename impl_type::value_type&
    {
      return v_[index];
    }

    math_vector& operator +=(math_vector const& b)
    {
      for (std::size_t i = 0; i<N; ++i)
        v_[i] += b[i];
      return *this;
    }

  #ifndef DONT_USE_EXPR_TEMPL

    template <typename LE, typename Op, typename RE>
    math_vector(math_vector_expr<LE,Op,RE>&& mve)
    {
      for (std::size_t i = 0; i < N; ++i)
        v_[i] = mve[i];
      std::cout << this << ": math_vector(expr: " << &mve << ")" << std::endl;
    }

    template <typename RightExpr>
    math_vector& operator =(RightExpr&& re)
    {
      for (std::size_t i = 0; i<N; ++i)
        v_[i] = re[i];
      return *this;
    }

    template <typename RightExpr>
    math_vector& operator +=(RightExpr&& re)
    {
      for (std::size_t i = 0; i<N; ++i)
        v_[i] += re[i];
      return *this;
    }

    template <typename RightExpr>
    auto operator +(RightExpr&& re) const ->
      math_vector_expr<
        math_vector const&, 
        plus_op<typename impl_type::value_type>,
        decltype(std::forward<RightExpr>(re))
      >
    {
      return 
        math_vector_expr<
          math_vector const&, 
          plus_op<typename impl_type::value_type>, 
          decltype(std::forward<RightExpr>(re))
        >(
          *this, 
          std::forward<RightExpr>(re)
        )
      ;
    }

  #endif // #ifndef DONT_USE_EXPR_TEMPL

  private:
    impl_type v_;
};

//===========================================================================

template <std::size_t N>
inline void swap(math_vector<N>& a, math_vector<N>& b)
{
  a.swap(b);
}

//===========================================================================

#ifdef DONT_USE_EXPR_TEMPL

template <std::size_t N>
inline math_vector<N> operator +(
  math_vector<N> const& a, 
  math_vector<N> const& b
)
{
  math_vector<N> retval(a);
  retval += b;
  return retval;
}

template <std::size_t N>
inline math_vector<N> operator +(
  math_vector<N>&& a, 
  math_vector<N> const& b
)
{
  a += b;
  return std::move(a);
}

template <std::size_t N>
inline math_vector<N> operator +(
  math_vector<N> const& a, 
  math_vector<N>&& b
)
{
  b += a;
  return std::move(b);
}

template <std::size_t N>
inline math_vector<N> operator +(
  math_vector<N>&& a, 
  math_vector<N>&& b
)
{
  a += std::move(b);
  return std::move(a);
}

#endif // #ifdef DONT_USE_EXPR_TEMPL

//===========================================================================

template <std::size_t N>
std::ostream& operator <<(std::ostream& os, math_vector<N> const& mv)
{
  os << '(';
  for (std::size_t i = 0; i < N; ++i)
    os << mv[i] << ((i+1 != N) ? ',' : ')');
  return os;
}

//===========================================================================

int main()
{
  using namespace std;

  try
  {
    {
      cout << "CASE 1:\n";
      math_vector<3> a{1.0, 1.1, 1.2};
      math_vector<3> b{2.0, 2.1, 2.2};
      math_vector<3> c{3.0, 3.1, 3.2};
      math_vector<3> d{4.0, 4.1, 4.2};
      math_vector<3> result = a + b + c + d;
      cout << '[' << &result << "]: " << result << "\n";
    }
    cout << endl;
    {
      cout << "CASE 2:\n";
      math_vector<3> result =
        math_vector<3>{1.0, 1.1, 1.2} +
        math_vector<3>{2.0, 2.1, 2.2} +
        math_vector<3>{3.0, 3.1, 3.2} +
        math_vector<3>{4.0, 4.1, 4.2}
      ;
      cout << '[' << &result << "]: " << result << "\n";
    }
  }
  catch (...)
  {
    return 1;
  }
}

//===========================================================================

C++ Expression Templates

What are C++ Expression Templates in simple terms?

Expression templates are a category of C++ template meta programming which delays evaluation of subexpressions until the full expression is known, so that optimizations (especially the elimination of temporaries) can be applied.

Are there books around that discuss numerical methods/computations using C++ Expression Templates?

I believe ET's were invented by Todd Veldhuizen who published a paper on it 15 years ago. (It seems that many older links to it are dead by now, but currently here is a version of it.) Some material about it is in David Vandevoorde's and Nicolai Josuttis' C++ Templates: The Complete Guide.

In what way, C++ Expression Templates are better than using pure C?

They allow you to write your code in an expressive high level way without losing performance. For example,

void f(const my_array<double> a1, const my_array<double> a2) 
{ 
  my_array<double> a3 = 1.2 * a1 + a1 * a2; 
  // ..
}

can be optimized all the way down to

for( my_array<double>::size_type idx=0; idx<a1.size(); ++idx ) 
  a3[idx] = 1.2*a1[idx] + a1[idx]*a2[idx]; 

which is faster, but harder to understand.

C++ expression templates lifetime

Instead of storing reference here:

template<typename E1, typename E2>
class ExpressionSum : public Expression<ExpressionSum<E1,E2>> {
    E1 const& u_; // <------| These are references
    E2 const& v_; // <------|
public:
    ExpressionSum(E1 const& u, E2 const& v) : u_(u), v_(v)
    { }

    // ...
};

These does not cause lifetime extension. The wikipedia article assume the expression template is never stored and only live in the same statement as the expression.

Store them as value:

template<typename E1, typename E2>
class ExpressionSum : public Expression<ExpressionSum<E1,E2>> {
    E1 u_; // <------| Fixed!
    E2 v_; // <------|
public:
    ExpressionSum(E1 const& u, E2 const& v) : u_(u), v_(v)
    { }

    // ...
};

You can also extend std::tuple to piggyback on it's EBO:

template<typename E1, typename E2>
class ExpressionSum : public Expression<ExpressionSum<E1,E2>>, private std::tuple<E1, E2> {
    auto u_() const -> E1 const& { return std::get<0>(*this); }
    auto v_() const -> E2 const& { return std::get<1>(*this); }
public:
    ExpressionSum(E1 const& u, E2 const& v) : std::tuple<E1, E2>(u, v)
    { }

    // ...
};

How to use expression-templates for specific types?

Simply use SFINAE to check if the type is an arithmethic one and specialize as needed.

Example:

template <typename E1, typename E2, typename Enable = void > class VecSum;

template <typename E1, typename E2>
class VecSum< E1, E2,
      typename std::enable_if_t<!std::is_arithmetic<E1>::value && !std::is_arithmetic<E2>::value>
      > : public VecExpression<VecSum<E1, E2> >
{
    E1 const& _u;
    E2 const& _v;

    public:

    VecSum(E1 const& u, E2 const& v) : _u(u), _v(v)
    {
        assert(u.size() == v.size());
    }

    double operator[](size_t i) const { return _u[i] + _v[i]; }
    size_t size()               const { return _v.size(); }
};

template <typename E1, typename E2>
class VecSum < E1, E2,
      typename std::enable_if_t< std::is_arithmetic<E1>::value && !std::is_arithmetic<E2>::value>
      > : public VecExpression<VecSum<E1, E2> >
{
    E1 const& _u;
    E2 const& _v;

    public:

    VecSum(E1 const& u, E2 const& v) : _u(u), _v(v)
    {
    }

    double operator[](size_t i) const { return _u + _v[i]; }
    size_t size()               const { return _v.size(); }
};

template <typename E1, typename E2>
class VecSum < E1, E2,
      typename std::enable_if_t< !std::is_arithmetic<E1>::value && std::is_arithmetic<E2>::value>
      > : public VecExpression<VecSum<E1, E2> >
{
    E1 const& _u;
    E2 const& _v;

    public:

    VecSum(E1 const& u, E2 const& v) : _u(u), _v(v)
    {
    }

    double operator[](size_t i) const { return _u[i] + _v; }
    size_t size()               const { return _u.size(); }
};

int main(){
    Vec v0 = { 1, 2, 3 ,4 };
    Vec v1 = {10, 20,30,40 };
    Vec v2 = {100,200,300,400 };

    {

        Vec sum = v0+v1+v2;
        Vec v3(4);

        for(int i=0;i<4;++i)
            v3[i]=sum[i];

        for(unsigned int i=0;i<v3.size();++i)
            std::cout << v3[i] << std::endl;
    }

    std::cout << "with lhs skalar" << std::endl;

    {
        Vec sum = 5 + 50 + v1;
        Vec v3(4);

        for(int i=0;i<4;++i)
            v3[i]=sum[i];

        for(unsigned int i=0;i<v3.size();++i)
            std::cout << v3[i] << std::endl;

    }

    std::cout << "with rhs skalar" << std::endl;

    {
        Vec sum = v1 + 5 + 50 ;
        Vec v3(4);

        for(int i=0;i<4;++i)
            v3[i]=sum[i];

        for(unsigned int i=0;i<v3.size();++i)
            std::cout << v3[i] << std::endl;

    }

}

Expression Templates - C++ Templates: The Complete Guide

Consider, for example, the right hand side of

x = 1.2*x + x*y;

What the quote says is that this is composed of two different categories.

The heavy array x and y objects are not defined within this expression, but rather before it:

Array<double> x(1000), y(1000);

So, as you build expressions using them, you don't have to worry whether they're still alive - they were defined beforehand. Since they're heavy, you want to capture them by reference, and, fortunately, their lifetime makes that possible.

Conversely, the lightweight A_Scale objects are generated within the expression (e.g., implicitly by the 1.2 above). Since they're temporaries, you have to worry about their lifetime. Since they're lightweight, it's not a problem.

That's the rationale for the traits class differentiating between them: the former are by reference, and the latter are by value (they are copied).

---

Related Topics