A Simple, Minimal, Example of C++ Expression Templates

I am preparing a lecture on efficient numerical computation in C++ and this obviously touches upon the technique expression template. I came up with the following, very stripped-down, example which still manages to illustrate expression templates are in general far more efficient than normal C++ value semantics.

The idea here is to illustrate addition of vectors (in this case I'm using plain old std::vector). The expression template mechanism below is in fact so simple here that it only supports the addition operation and there are several shortcomings compared to a proper (and much longer) implementation. However it shows three clear advantages that are critical in high-performance code:

  • No creation of temporaries of length order N (where N is the number of elements)
  • Only a single iteration need be made through all of the vectors being added
  • If only a subset of elements of the result is requested, then only these elements are computed

If you are interested in examples of the technique in practice have a look at Boost.UBLAS, Armadillo, Eigen and Boost.Proto.

// Bojan Nikolic <bojan@bnikolic.co.uk>
// A minimal expression template example

#include <vector>
#include <iostream>

/// For comparision, this is an addition operator on stardard
/// std::vectors the old-fashioned way
template< class T>
std::vector<T> operator+(const std::vector<T> &a,
                         const std::vector<T> &b)
  std::vector<T> res(a.size());
  for(size_t i=0; i<a.size(); ++i)
  return res;

// This is trait class for ops which actually just represent the value
// of a single argument
struct valueop {};
// This is for the addition operation
struct addop {};

template<class E1=std::vector<double>,
         class E2=std::vector<double>,
         class op=valueop>
struct binop
  const E1 &left;
  const E2 &right;

  binop(const E1 &left,
        const E2 &right,
        op opval):

  /// Construction with just one value and default operator is taken
  /// to mean that this is a "terminal" value
  binop(const E1 &val);

  ::binop(const std::vector<double> &val):

/// Evaluete the i-th element of an expression
template<class E1, class E2, class op>
double eval(const binop<E1, E2, op> &o, size_t i)
  // Default implementation returns 0 as that is identity for
  // addition/substraction
  return 0;

/// Evaluation of a "valueop" just returns the value
template<class T>
double eval(const binop<T, T,  valueop> &o,
            size_t i)
  return o.left[i];

/// Evaluation of "addop" adds the corresponding elements
template<class E1, class E2>
double eval(const binop<E1, E2, addop> &o, size_t i)
  return eval(o.left,i)+eval(o.right,i);

template<class E1, class E2>
binop<E1, E2, addop>  operator+ (const E1 &left,
                                 const E2 &right)
  return binop<E1, E2, addop>(left, right, addop());

template <class T>
binop<T, T, valueop> mkVal(T &val)
  return binop<T, T, valueop>(val);

int main(void)
  /// Some vectors to use in the example
  std::vector<double> a(10, 1.0), b(10, 2.0);

  // Addition the old vay, entire vector is computed

  // Now do it the expression template way
  binop<> ba(a), bb(b);
  // If you've got C++ 0x features then you can do:
  // auto bc=mkVal(a);

  // Print 3rd element of the sum, only the sum for the third element
  // is computed
  std::cout<<eval(ba+bb+ba+bb, 3)<<std::endl;