Extending MatrixBase

In this section we will see how to add custom methods to MatrixBase. Since all expressions and matrix types inherit MatrixBase, adding a method to MatrixBase make it immediately available to all expressions ! A typical use case is, for instance, to make Eigen compatible with another API.

You certainly know that in C++ it is not possible to add methods to an extending class. So how that's possible ? Here the trick is to include in the declaration of MatrixBase a file defined by the preprocessor token EIGEN_MATRIXBASE_PLUGIN:

class MatrixBase {
  // ...
  #ifdef EIGEN_MATRIXBASE_PLUGIN
  #include EIGEN_MATRIXBASE_PLUGIN
  #endif
};

Therefore to extend MatrixBase with you own methods you just have to create a file with your method declaration and define EIGEN_MATRIXBASE_PLUGIN before you include any Eigen's header file.

Here is an example of such an extension file:
MatrixBaseAddons.h

inline Scalar at(uint i, uint j) const { return this->operator()(i,j); }
inline Scalar& at(uint i, uint j) { return this->operator()(i,j); }
inline Scalar at(uint i) const { return this->operator[](i); }
inline Scalar& at(uint i) { return this->operator[](i); }
inline RealScalar squaredLength() const { return squaredNorm(); }
inline RealScalar length() const { return norm(); }
inline RealScalar invLength(void) const { return fast_inv_sqrt(squaredNorm()); }
template<typename OtherDerived>
inline Scalar squaredDistanceTo(const MatrixBase<OtherDerived>& other) const
{ return (derived() - other.derived()).squaredNorm(); }
template<typename OtherDerived>
inline RealScalar distanceTo(const MatrixBase<OtherDerived>& other) const
{ return ei_sqrt(derived().squaredDistanceTo(other)); }
inline void scaleTo(RealScalar l) { RealScalar vl = norm(); if (vl>1e-9) derived() *= (l/vl); }
inline Transpose<Derived> transposed() {return transpose();}
inline const Transpose<Derived> transposed() const {return transpose();}
inline uint minComponentId(void) const  { int i; minCoeff(&i); return i; }
inline uint maxComponentId(void) const  { int i; maxCoeff(&i); return i; }
template<typename OtherDerived>
void makeFloor(const MatrixBase<OtherDerived>& other) { derived() = derived().cwise().min(other.derived()); }
template<typename OtherDerived>
void makeCeil(const MatrixBase<OtherDerived>& other) { derived() = derived().cwise().max(other.derived()); }
const const CwiseUnaryOp<ei_scalar_add_op<Scalar>, Derived>
operator+(const Scalar& scalar) const { return cwise() + scalar; }
friend const CwiseUnaryOp<ei_scalar_add_op<Scalar>, Derived>
operator+(const Scalar& scalar, const MatrixBase<Derived>& mat) { return mat + scalar; }

Then one can the following declaration in the config.h or whatever prerequisites header file of his project:

#define EIGEN_MATRIXBASE_PLUGIN "MatrixBaseAddons.h"

Using custom scalar types

By default, Eigen currently supports the following scalar types: int, float, double, std::complex<float>, std::complex<double>, long double, long long int (64 bits integers), and bool. The long double is especially useful on x86-64 systems or when the SSE2 instruction set is enabled because it enforces the use of x87 registers with extended accuracy.

In order to add support for a custom type T you need: 1 - make sure the common operator (+,-,*,/,etc.) are supported by the type T 2 - add a specialization of struct Eigen::NumTraits<T> (see NumTraits) 3 - define a couple of math functions for your type such as: ei_sqrt, ei_abs, etc... (see the file Eigen/src/Core/MathFunctions.h)

Here is a concrete example adding support for the Adolc's adouble type. Adolc is an automatic differentiation library. The type adouble is basically a real value tracking the values of any number of partial derivatives.

#ifndef ADLOCSUPPORT_H
#define ADLOCSUPPORT_H
#define ADOLC_TAPELESS
#include <adolc/adouble.h>
#include <Eigen/Core>
namespace Eigen {
template<> struct NumTraits<adtl::adouble>
{
  typedef adtl::adouble Real;
  typedef adtl::adouble FloatingPoint;
  enum {
    IsComplex = 0,
    HasFloatingPoint = 1,
    ReadCost = 1,
    AddCost = 1,
    MulCost = 1
  };
};
}
// the Adolc's type adouble is defined in the adtl namespace
// therefore, the following ei_* functions *must* be defined
// in the same namespace
namespace adtl {
  inline const adouble& ei_conj(const adouble& x)  { return x; }
  inline const adouble& ei_real(const adouble& x)  { return x; }
  inline adouble ei_imag(const adouble&)    { return 0.; }
  inline adouble ei_abs(const adouble&  x)  { return fabs(x); }
  inline adouble ei_abs2(const adouble& x)  { return x*x; }
  inline adouble ei_sqrt(const adouble& x)  { return sqrt(x); }
  inline adouble ei_exp(const adouble&  x)  { return exp(x); }
  inline adouble ei_log(const adouble&  x)  { return log(x); }
  inline adouble ei_sin(const adouble&  x)  { return sin(x); }
  inline adouble ei_cos(const adouble&  x)  { return cos(x); }
  inline adouble ei_pow(const adouble& x, adouble y)  { return pow(x, y); }
}
#endif // ADLOCSUPPORT_H

Preprocessor directives

You can control some aspects of Eigen by defining the following preprocessor tokens them before including any of Eigen's headers.

EIGEN_NO_DEBUG disables Eigen assertions. Like NDEBUG but only affects Eigen's assertions.
EIGEN_DONT_VECTORIZE disables explicit vectorization when defined.
EIGEN_UNROLLING_LIMIT defines the maximal instruction counts to enable meta unrolling of loops. Set it to zero to disable unrolling. The default is 100.
EIGEN_DEFAULT_TO_ROW_MAJOR the default storage order for matrices becomes row-major instead of column-major.
EIGEN_TUNE_FOR_CPU_CACHE_SIZE represents the maximal size in Bytes of L2 blocks. Since several blocks have to stay concurently in L2 cache, this value should correspond to at most 1/4 of the size of L2 cache.
EIGEN_NO_STATIC_ASSERT replaces compile time static assertions by runtime assertions
EIGEN_MATRIXBASE_PLUGIN see Extending MatrixBase