Expression templates

Expression templates is a C++ template metaprogramming technique in which templates are used to represent part of an expression. Typically, the template itself represents a particular kind of operation, while the parameters represent the operands to which the operation applies. The expression template can then be evaluated at a later time, or passed to a function. The technique was invented independently by Todd Veldhuizen and David Vandevoorde.

For example, consider a library representing vectors with a class Vec. It is natural to want to overload operator- and operator* so you could write Vec x = alpha*(u - v); where alpha is a scalar and u and v are Vecs. A naive implementation would have operator- and operator* return Vecs. However, then the above expression would mean creating a temporary for u-v then another temporary for alpha times that first temporary, then assigning that to x. Even with the return value optimization this will allocate memory at least twice: once for the temporary u-v and once for the result of the overall expression.

Expression templates delay evaluation so the expression Vec x = alpha*(u-v); essentially generates at compile time a new Vec constructor. It is as if this constructor takes a scalar and two Vecs by reference; it allocates the necessary memory and then performs the computation. Thus only one memory allocation is performed.

An example implementation of expression templates is as follows (using the curiously recurring template pattern as is used by Boost.uBLAS):

#include <vector>
#include <cassert>

/* The template parameter is named 'E' for 'Expression' */

template <typename E>
// A CRTP base class for Vecs with a size and indexing:
class VecExpression {
public:
  typedef std::vector<double>         container_type;
  typedef typename container_type::size_type   size_type;
  typedef typename container_type::value_type  value_type;
  typedef typename container_type::reference   reference;

  size_type  size()                  const { return static_cast<E const&>(*this).size(); }
  value_type operator[](size_type i) const { return static_cast<E const&>(*this)[i];     }

  operator E&()             { return static_cast<      E&>(*this); }
  operator E const&() const { return static_cast<const E&>(*this); }
};

// The actual Vec class:
class Vec : public VecExpression<Vec> {
   container_type _data;
public:
   reference  operator[](size_type i)       { return _data[i]; }
   value_type operator[](size_type i) const { return _data[i]; }
   size_type  size()                  const { return _data.size(); }

   Vec(size_type n) : _data(n) {} // Construct a given size:

   // Construct from any VecExpression:
   template <typename E>
   Vec(VecExpression<E> const& vec) {
     E const& v = vec;
     _data.resize(v.size());
     for (size_type i = 0; i != v.size(); ++i) {
       _data[i] = v[i];
     }
   }
};

template <typename E1, typename E2>
class VecDifference : public VecExpression<VecDifference<E1, E2> > {
  E1 const& _u;
  E2 const& _v;
public:
  typedef Vec::size_type size_type;
  typedef Vec::value_type value_type;
  VecDifference(VecExpression<E1> const& u, VecExpression<E2> const& v) : _u(u), _v(v) {
    assert(u.size() == v.size());
  }
  size_type size() const { return _v.size(); }
  value_type operator[](Vec::size_type i) const { return _u[i] - _v[i]; }
};
  
template <typename E>
class VecScaled : public VecExpression<VecScaled<E> > {
  double _alpha; 
  E const& _v;
public:
  VecScaled(double alpha, VecExpression<E> const& v) : _alpha(alpha), _v(v) {}
  Vec::size_type size() const { return _v.size(); }
  Vec::value_type operator[](Vec::size_type i) const { return _alpha * _v[i]; }
};

// Now we can overload operators:

template <typename E1, typename E2>
VecDifference<E1,E2> const
operator-(VecExpression<E1> const& u, VecExpression<E2> const& v) {
   return VecDifference<E1,E2>(u,v);
}

template <typename E>
VecScaled<E> const
operator*(double alpha, VecExpression<E> const& v) {
   return VecScaled<E>(alpha,v);
}

With the above definitions, the expression alpha*(u-v) is of type

VecScaled<VecDifference<Vec,Vec> >

so calling Vec x = alpha * (u-v) calls the constructor that takes a VecExpression<VecScaled<VecDifference<Vec,Vec> > >. Each line of the for loop then expands from

_data[i] = v[i];

to essentially

_data[i] = alpha * (u[i] - v[i]);

with no temporaries needed and only one pass through each memory block.

Libraries using Expression templates

• Boost uBLAS
• Eigen
• Blitz++