eigen/doc/InsideEigenExample.dox

namespace Eigen {

/** \page TopicInsideEigenExample What happens inside Eigen, on a simple example

\eigenAutoToc

<hr>


Consider the following example program:

\code
#include<Eigen/Core>

int main()
{
  int size = 50;
  // VectorXf is a vector of floats, with dynamic size.
  Eigen::VectorXf u(size), v(size), w(size);
  u = v + w;
}
\endcode

The goal of this page is to understand how Eigen compiles it, assuming that SSE2 vectorization is enabled (GCC option -msse2).

\section WhyInteresting Why it's interesting

Maybe you think, that the above example program is so simple, that compiling it shouldn't involve anything interesting. So before starting, let us explain what is nontrivial in compiling it correctly -- that is, producing optimized code -- so that the complexity of Eigen, that we'll explain here, is really useful.

Look at the line of code
\code
  u = v + w;   //   (*)
\endcode

The first important thing about compiling it, is that the arrays should be traversed only once, like
\code
  for(int i = 0; i < size; i++) u[i] = v[i] + w[i];
\endcode
The problem is that if we make a naive C++ library where the VectorXf class has an operator+ returning a VectorXf, then the line of code (*) will amount to:
\code
  VectorXf tmp = v + w;
  VectorXf u = tmp;
\endcode
Obviously, the introduction of the temporary \a tmp here is useless. It has a very bad effect on performance, first because the creation of \a tmp requires a dynamic memory allocation in this context, and second as there are now two for loops:
\code
  for(int i = 0; i < size; i++) tmp[i] = v[i] + w[i];
  for(int i = 0; i < size; i++) u[i] = tmp[i];
\endcode
Traversing the arrays twice instead of once is terrible for performance, as it means that we do many redundant memory accesses.

The second important thing about compiling the above program, is to make correct use of SSE2 instructions. Notice that Eigen also supports AltiVec and that all the discussion that we make here applies also to AltiVec.

SSE2, like AltiVec, is a set of instructions allowing to perform computations on packets of 128 bits at once. Since a float is 32 bits, this means that SSE2 instructions can handle 4 floats at once. This means that, if correctly used, they can make our computation go up to 4x faster.

However, in the above program, we have chosen size=50, so our vectors consist of 50 float's, and 50 is not a multiple of 4. This means that we cannot hope to do all of that computation using SSE2 instructions. The second best thing, to which we should aim, is to handle the 48 first coefficients with SSE2 instructions, since 48 is the biggest multiple of 4 below 50, and then handle separately, without SSE2, the 49th and 50th coefficients. Something like this:

\code
  for(int i = 0; i < 4*(size/4); i+=4) u.packet(i)  = v.packet(i) + w.packet(i);
  for(int i = 4*(size/4); i < size; i++) u[i] = v[i] + w[i];
\endcode

So let us look line by line at our example program, and let's follow Eigen as it compiles it.

\section ConstructingVectors Constructing vectors

Let's analyze the first line:

\code
  Eigen::VectorXf u(size), v(size), w(size);
\endcode

First of all, VectorXf is the following typedef:
\code
  typedef Matrix<float, Dynamic, 1> VectorXf;
\endcode

The class template Matrix is declared in src/Core/util/ForwardDeclarations.h with 6 template parameters, but the last 3 are automatically determined by the first 3. So you don't need to worry about them for now. Here, Matrix\<float, Dynamic, 1\> means a matrix of floats, with a dynamic number of rows and 1 column.

The Matrix class inherits a base class, MatrixBase. Don't worry about it, for now it suffices to say that MatrixBase is what unifies matrices/vectors and all the expressions types -- more on that below.

When we do
\code
  Eigen::VectorXf u(size);
\endcode
the constructor that is called is Matrix::Matrix(Index), in src/Core/Matrix.h. Besides some assertions, all it does is to construct the \a m_storage member, which is of type DenseStorage\<float, Dynamic, Dynamic, 1, 0\> (the last parameter carries the Options).

You may wonder, isn't it overengineering to have the storage in a separate class? The reason is that the Matrix class template covers all kinds of matrices and vector: both fixed-size and dynamic-size. The storage method is not the same in these two cases. For fixed-size, the matrix coefficients are stored as a plain member array. For dynamic-size, the coefficients will be stored as a pointer to a dynamically-allocated array. Because of this, we need to abstract storage away from the Matrix class. That's DenseStorage.

Let's look at this constructor, in src/Core/DenseStorage.h. You can see that there are many partial template specializations of the DenseStorage implementation here, treating separately the cases where dimensions are Dynamic or fixed at compile-time. The partial specialization that we are looking at is the one for a dynamic size with a fixed number of columns:
\code
template<typename T, int Cols, int Options> class DenseStorage_impl<T, Dynamic, Dynamic, Cols, Options>
\endcode

Here, the constructor called is DenseStorage_impl(Index size, Index rows, Index cols)
with size=50, rows=50, cols=1:
\code
constexpr DenseStorage_impl(Index size, Index rows, Index /*cols*/)
    : m_data(conditional_aligned_new_auto<T, Align>(size)), m_rows(rows) {}
\endcode

Here, the \a m_data member is the actual array of coefficients of the matrix. As you see, it is dynamically allocated. Rather than calling new[] or malloc(), we have our own internal::conditional_aligned_new_auto defined in src/Core/util/Memory.h. What it does is that if vectorization is enabled, then it allocates a suitably aligned array (e.g. 16-byte-aligned for SSE2, more for AVX or AVX512), as that is very useful for vectorization. If no alignment is required, it amounts to the standard new[].

As you can see, the constructor also sets the \a m_rows member to \a rows. Notice that there is no \a m_cols member: indeed, in this partial specialization of DenseStorage_impl, we know the number of columns at compile-time, since the Cols template parameter is different from Dynamic. Namely, in our case, Cols is 1, which is to say that our vector is just a matrix with 1 column. Hence, there is no need to store the number of columns as a runtime variable.

When you call VectorXf::data() to get the pointer to the array of coefficients, it returns DenseStorage::data() which returns the \a m_data member.

When you call VectorXf::size() to get the size of the vector, this is actually a method in the base class MatrixBase. It determines that the vector is a column-vector, since ColsAtCompileTime==1 (this comes from the template parameters in the typedef VectorXf). It deduces that the size is the number of rows, so it returns VectorXf::rows(), which returns DenseStorage::rows(), which returns the \a m_rows member, which was set to \a size by the constructor.

\section ConstructionOfSumXpr Construction of the sum expression

Now that our vectors are constructed, let's move on to the next line:

\code
u = v + w;
\endcode

The executive summary is that operator+ returns a "sum of vectors" expression, but doesn't actually perform the computation. It is the operator=, whose call occurs thereafter, that does the computation.

Let us now see what Eigen does when it sees this:

\code
v + w
\endcode

Here, v and w are of type VectorXf, which is a typedef for a specialization of Matrix (as we explained above), which is a subclass of MatrixBase. So what is being called is

\code
MatrixBase::operator+(const MatrixBase&)
\endcode

The return type of this operator is
\code
CwiseBinaryOp<internal::scalar_sum_op<float, float>, const VectorXf, const VectorXf>
\endcode
The CwiseBinaryOp class is our first encounter with an expression template. As we said, the operator+ doesn't by itself perform any computation, it just returns an abstract "sum of vectors" expression. Since there are also "difference of vectors" and "coefficient-wise product of vectors" expressions, we unify them all as "coefficient-wise binary operations", which we abbreviate as "CwiseBinaryOp". "Coefficient-wise" means that the operations is performed coefficient by coefficient. "binary" means that there are two operands -- we are adding two vectors with one another.

Now you might ask, what if we did something like

\code
v + w + u;
\endcode

The first v + w would return a CwiseBinaryOp as above, so in order for this to compile, we'd need to define an operator+ also in the class CwiseBinaryOp... at this point it starts looking like a nightmare: are we going to have to define all operators in each of the expression classes (as you guessed, CwiseBinaryOp is only one of many) ? This looks like a dead end!

The solution is that CwiseBinaryOp itself, as well as Matrix and all the other expression types, is a subclass of MatrixBase. So it is enough to define once and for all the operators in class MatrixBase.

Since MatrixBase is the common base class of different subclasses, the aspects that depend on the subclass must be abstracted from MatrixBase. This is called polymorphism.

The classical approach to polymorphism in C++ is by means of virtual functions. This is dynamic polymorphism. Here we don't want dynamic polymorphism because the whole design of Eigen is based around the assumption that all the complexity, all the abstraction, gets resolved at compile-time. This is crucial: if the abstraction can't get resolved at compile-time, Eigen's compile-time optimization mechanisms become useless, not to mention that if that abstraction has to be resolved at runtime it'll incur an overhead by itself.

Here, what we want is to have a single class MatrixBase as the base of many subclasses, in such a way that each MatrixBase object (be it a matrix, or vector, or any kind of expression) knows at compile-time (as opposed to run-time) of which particular subclass it is an object (i.e. whether it is a matrix, or an expression, and what kind of expression).

The solution is the <a href="http://en.wikipedia.org/wiki/Curiously_Recurring_Template_Pattern">Curiously Recurring Template Pattern</a>. Let's do the break now. Hopefully you can read this wikipedia page during the break if needed, but it won't be allowed during the exam.

In short, MatrixBase takes a template parameter \a Derived. Whenever we define a subclass Subclass, we actually make Subclass inherit MatrixBase\<Subclass\>. The point is that different subclasses inherit different MatrixBase types. Thanks to this, whenever we have an object of a subclass, and we call on it some MatrixBase method, we still remember even from inside the MatrixBase method which particular subclass we're talking about.

This means that we can put almost all the methods and operators in the base class MatrixBase, and have only the bare minimum in the subclasses. If you look at the subclasses in Eigen, like for instance the CwiseBinaryOp class, they have very few methods. There are coeff() and sometimes coeffRef() methods for access to the coefficients, there are rows() and cols() methods returning the number of rows and columns, but there isn't much more than that. All the meat is in MatrixBase, so it only needs to be coded once for all kinds of expressions, matrices, and vectors.

So let's end this digression and come back to the piece of code from our example program that we were currently analyzing,

\code
v + w
\endcode

Now that MatrixBase is a good friend, let's write fully the prototype of the operator+ that gets called here (the actual declaration is generated by the EIGEN_MAKE_CWISE_BINARY_OP macro in the plugin file src/plugins/CommonCwiseBinaryOps.inc, included from within MatrixBase; slightly simplified, it amounts to):

\code
template<typename Derived>
class MatrixBase
{
  // ...

  template<typename OtherDerived>
  const CwiseBinaryOp<internal::scalar_sum_op<typename internal::traits<Derived>::Scalar,
                                              typename internal::traits<OtherDerived>::Scalar>,
                      const Derived, const OtherDerived>
  operator+(const MatrixBase<OtherDerived> &other) const;

  // ...
};
\endcode

Here of course, \a Derived and \a OtherDerived are VectorXf.

As we said, CwiseBinaryOp is also used for other operations such as subtraction, so it takes another template parameter determining the operation that will be applied to coefficients. This template parameter is a functor, that is, a class in which we have an operator() so it behaves like a function. Here, the functor used is internal::scalar_sum_op. It is defined in src/Core/functors/BinaryFunctors.h.

Let us now explain the internal::traits here. The internal::scalar_sum_op class takes two template parameters: the types of the numbers to add (which may differ, e.g. when adding a real and a complex matrix). Here of course we want to pass the scalar type (a.k.a. numeric type) of VectorXf, which is \c float, twice. How do we determine which is the scalar type of \a Derived ? Throughout Eigen, all matrix and expression types define a typedef \a Scalar which gives its scalar type. For example, VectorXf::Scalar is a typedef for \c float. So here, if life was easy, we could find the numeric type of \a Derived as just
\code
typename Derived::Scalar
\endcode
Unfortunately, we can't do that here, as the compiler would complain that the type Derived hasn't yet been defined. So we use a workaround: in src/Core/util/ForwardDeclarations.h, we declared (not defined!) all our subclasses, like Matrix, and we also declared the following class template:
\code
template<typename T> struct internal::traits;
\endcode
In src/Core/Matrix.h, right \em before the definition of class Matrix, we define a partial specialization of internal::traits for T=Matrix\<any template parameters\>. In this specialization of internal::traits, we define the Scalar typedef. So when we actually define Matrix, it is legal to refer to "typename internal::traits\<Matrix\>::Scalar".

Anyway, we have declared our operator+. In our case, where \a Derived and \a OtherDerived are VectorXf, the above declaration amounts to:
\code
class MatrixBase<VectorXf>
{
  // ...

  const CwiseBinaryOp<internal::scalar_sum_op<float, float>, const VectorXf, const VectorXf>
  operator+(const MatrixBase<VectorXf> &other) const;

  // ...
};
\endcode

Let's now jump to src/Core/CwiseBinaryOp.h to see how it is defined. As you can see there, all it does is to return a CwiseBinaryOp object, and this object is just storing references to the left-hand-side and right-hand-side expressions -- here, these are the vectors \a v and \a w. Well, the CwiseBinaryOp object is also storing an instance of the (empty) functor class, but you shouldn't worry about it as that is a minor implementation detail.

Thus, the operator+ hasn't performed any actual computation. To summarize, the operation \a v + \a w just returned an object of type CwiseBinaryOp which did nothing else than just storing references to \a v and \a w.

\section Assignment The assignment

At this point, the expression \a v + \a w has finished evaluating, so, in the process of compiling the line of code
\code
u = v + w;
\endcode
we now enter the operator=.

What operator= is being called here? The vector u is an object of class VectorXf, i.e. Matrix. In src/Core/Matrix.h, inside the definition of class Matrix, we see this:
\code
    template<typename OtherDerived>
    EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Matrix& operator=(const DenseBase<OtherDerived>& other)
    {
      return Base::_set(other);
    }
\endcode
Here, Base is a typedef for PlainObjectBase\<Matrix\>. The _set() helper, in src/Core/PlainObjectBase.h, forwards to the central assignment entry point:
\code
    internal::call_assignment(this->derived(), other.derived());
\endcode
This function literally reads "assign a sum of two VectorXf's to another VectorXf", and everything from here on lives in src/Core/AssignEvaluator.h. This file implements %Eigen's evaluator-based assignment mechanism (introduced in %Eigen 3.3): expressions describe \em what to compute, while separate \em evaluator classes, one specialization of internal::evaluator\<\> per expression type, describe \em how to compute coefficients and packets of it.

The first stage is a compile-time dispatch on whether the right-hand side must be assumed to alias the destination:
\code
template <typename Dst, typename Src, typename Func, std::enable_if_t<evaluator_assume_aliasing<Src>::value, int> = 0>
constexpr void call_assignment(Dst& dst, const Src& src, const Func& func) {
  typename plain_matrix_type<Src>::type tmp(src);
  call_assignment_no_alias(dst, tmp, func);
}

template <typename Dst, typename Src, typename Func, std::enable_if_t<!evaluator_assume_aliasing<Src>::value, int> = 0>
constexpr void call_assignment(Dst& dst, const Src& src, const Func& func) {
  call_assignment_no_alias(dst, src, func);
}
\endcode
As explained in \ref TopicLazyEvaluation, certain expressions, most notably the matrix Product expression, are evaluated into a temporary before being assigned, in order to avoid aliasing issues when doing "m = m * m;". The internal::evaluator_assume_aliasing trait identifies them. Our CwiseBinaryOp sum is not one of them -- we said since the beginning that we didn't want a temporary to be introduced here -- so we take the second overload and go straight to call_assignment_no_alias. This is also, incidentally, exactly what writing u.noalias() = v + w would have done directly.

Inside call_assignment_no_alias, a compile-time constant NeedToTranspose handles the special case of assigning a row-vector to a column-vector (the one exception to the rule that assignment requires matching dimensions); here both sides are column vectors, so it is false. After a few static assertions, the call lands in the Assignment class template, which dispatches on the \em shapes of the two sides. Both sides here have dense shape, so the Dense2Dense specialization runs, which (ignoring a runtime aliasing check in debug builds) simply calls call_dense_assignment_loop:
\code
template <typename DstXprType, typename SrcXprType, typename Functor>
constexpr void call_dense_assignment_loop(DstXprType& dst, const SrcXprType& src, const Functor& func) {
  using DstEvaluatorType = evaluator<DstXprType>;
  using SrcEvaluatorType = evaluator<SrcXprType>;

  SrcEvaluatorType srcEvaluator(src);
  resize_if_allowed(dst, src, func);
  DstEvaluatorType dstEvaluator(dst);

  using Kernel = generic_dense_assignment_kernel<DstEvaluatorType, SrcEvaluatorType, Functor>;
  Kernel kernel(dstEvaluator, srcEvaluator, func, dst.const_cast_derived());

  dense_assignment_loop<Kernel>::run(kernel);
}
\endcode

This is where the evaluators appear. The evaluator of our CwiseBinaryOp expression stores the functor and, recursively, the evaluators of its two operands; the evaluator of a Matrix simply stores the data pointer. The \a kernel object bundles the destination evaluator, the source evaluator, and the assignment functor (here internal::assign_op, meaning plain assignment; compound assignments like += reuse this exact machinery with a different functor). The kernel exposes two essential operations:
\code
  void assignCoeff(Index index) {
    m_functor.assignCoeff(m_dst.coeffRef(index), m_src.coeff(index));
  }
  template <int StoreMode, int LoadMode, typename Packet>
  void assignPacket(Index index) {
    m_functor.template assignPacket<StoreMode>(&m_dst.coeffRef(index), m_src.template packet<LoadMode, Packet>(index));
  }
\endcode

The last line hands the kernel to dense_assignment_loop, whose run() function selects the evaluation strategy. A helper class, internal::copy_using_evaluator_traits, computes two compile-time constants from the evaluators' capabilities: \a Traversal (how to iterate: coefficient by coefficient, linearly, by aligned packets, ...) and \a Unrolling (whether to fully unroll the loop at compile time). In our example \a Traversal is \a LinearVectorizedTraversal -- both sides can be addressed linearly by a single integer index and support packet access -- and \a Unrolling is \a NoUnrolling, which is obvious since our vectors have dynamic size so there's no way to unroll the loop at compile-time.

So the partial specialization that runs is, in slightly simplified form:
\code
template <typename Kernel>
struct dense_assignment_loop_impl<Kernel, LinearVectorizedTraversal, NoUnrolling> {
  static void run(Kernel& kernel) {
    const Index size = kernel.size();
    const Index alignedStart = DstIsAligned ? 0 : first_aligned<Alignment>(kernel.dstDataPtr(), size);
    const Index alignedEnd = alignedStart + numext::round_down(size - alignedStart, PacketSize);

    head_loop::run(kernel, 0, alignedStart);                      // scalar (or partial-packet) prologue

    for (Index index = alignedStart; index < alignedEnd; index += PacketSize)
      kernel.template assignPacket<Alignment, SrcAlignment, PacketType>(index);

    tail_loop::run(kernel, alignedEnd, size);                     // scalar (or partial-packet) epilogue
  }
};
\endcode

Here's how it works. As we said at the beginning, vectorization works with blocks of 4 floats with SSE, so \a PacketSize is 4. There are two potential problems that the loop structure deals with:
\li first, vectorized stores are most efficient when they are aligned to the packet size. So when writing to the coefficients of the destination, we want to group them into packets of 4 such that each packet is 128-bit-aligned. In general, this requires skipping a few coefficients at the beginning, which is the purpose of \a alignedStart and the head loop. In our case the destination is a VectorXf, and remember that in the construction of the vectors we allocated aligned arrays. The evaluator advertises this alignment at compile time, so \a alignedStart is zero and the head loop is avoided altogether without any runtime check.
\li second, the number of coefficients is not in general a multiple of \a PacketSize. Here, there are 50 coefficients and \a PacketSize is 4, so \a alignedEnd is 48 and the last 2 coefficients are handled by the tail loop.

Now let's follow one call to kernel.assignPacket(index) all the way down. The assignment functor internal::assign_op, from src/Core/functors/AssignmentFunctors.h, stores the source packet into the destination memory:
\code
template <int Alignment, typename Packet>
void assign_op::assignPacket(DstScalar* a, const Packet& b) const {
  pstoret<DstScalar, Packet, Alignment>(a, b);
}
\endcode
internal::pstoret, in src/Core/GenericPacketMath.h, is a dispatcher selecting between the aligned internal::pstore and the unaligned internal::pstoreu at compile time. Since \a Alignment says the store is aligned, it calls internal::pstore, whose SSE specialization in src/Core/arch/SSE/PacketMath.h is essentially:
\code
template <>
inline void pstore<float>(float* to, const Packet4f& from) { _mm_store_ps(to, from); }
\endcode
Here, Packet4f is %Eigen's name for the SSE type __m128 ("packet of 4 floats"), and _mm_store_ps is an SSE intrinsic representing a single SSE instruction. The packet_traits\<float\> specialization in the same file is what declared, at compile time, that packets of floats have size 4, which is where \a PacketSize came from.

The source packet was produced by m_src.packet(index), where m_src is the evaluator of our sum expression. Its packet() function, in src/Core/CoreEvaluators.h, reads one packet from each operand's evaluator and combines them with the functor:
\code
template <int LoadMode, typename PacketType>
PacketType binary_evaluator<CwiseBinaryOp>::packet(Index index) const {
  return m_d.func().packetOp(m_d.lhsImpl.template packet<LoadMode, PacketType>(index),
                             m_d.rhsImpl.template packet<LoadMode, PacketType>(index));
}
\endcode
Here, \a lhsImpl and \a rhsImpl are the evaluators of the vectors \a v and \a w. Their packet() function loads 4 consecutive floats with internal::ploadt, the load-side analogue of internal::pstoret, which resolves to the _mm_load_ps intrinsic. Finally, the functor's packetOp(), in src/Core/functors/BinaryFunctors.h, is:
\code
template <typename Packet>
Packet scalar_sum_op::packetOp(const Packet& a, const Packet& b) const {
  return internal::padd(a, b);
}
\endcode
and internal::padd's SSE specialization is a single _mm_add_ps instruction.

To summarize, the vectorized loop has been compiled to the following code: for \a index taking the values 0, 4, 8, ..., 44 -- twelve iterations covering the first 48 coefficients -- read the packet of 4 floats starting at v[index] and the packet starting at w[index] using two _mm_load_ps SSE instructions, add them together using an _mm_add_ps instruction, then store the result using an _mm_store_ps instruction.

There remains the tail loop handling the last few (here, the last 2) coefficients. Depending on the instruction set it either uses partial-packet loads and stores, or falls back to kernel.assignCoeff(index), which does the same walk through the evaluators coefficient by coefficient: assignPacket() becomes assignCoeff(), packet() becomes coeff(), and pstore() becomes a plain scalar store. If you followed us this far, you can probably understand this part by yourself.

We see that all the C++ abstraction of Eigen goes away during compilation and that we indeed are precisely controlling which assembly instructions we emit. Such is the beauty of C++! Since we have such precise control over the emitted assembly instructions, but such complex logic to choose the right instructions, we can say that Eigen really behaves like an optimizing compiler. If you prefer, you could say that Eigen behaves like a script for the compiler. In a sense, C++ template metaprogramming is scripting the compiler -- and it's been shown that this scripting language is Turing-complete. See <a href="http://en.wikipedia.org/wiki/Template_metaprogramming"> Wikipedia</a>.

*/

}