eigen2 porting tips

Hi,

just a few notes on differences between Eigen1 and 2 which affect Avogadro.

vector.ortho() becomes vector.someOrthogonal()
class MatrixP becomes class Transform
computeFittingHyperplane() becomes fitHyperplane()
and moreover the parameter “points” becomes an array of pointers
to vectors. This part of the API is rather likely to evolve but
since it’s only one call throughout avogadro, it should be OK :slight_smile:
someAntecedent() becomes lu().solve()

for the rotations, the API is now much richer: see AngleAxis and Quaternion.
In particular we have orthogonal-basis-vector expressions with UnitX() etc.
See the documentation of AngleAxis.

In Eigen1, the operators +,-,* were inefficient as they returned objects by
value. So we provided C-style functions like multiply() to avoid paying that
cost. No longer in Eigen2 : they just return abstract expressions, there is
no overhead anymore. (Luckily, I seem to remember that in Avogadro we never
really bothered to use the faster C-style functions. All the better then,
porting should be easier and you’ll see a bigger benefit from moving to
eigen2).

The storage order of matrices is still column-major by default but you can
change that per-matrix if you need to (OpenGL uses column-major so in general
that’s what you want).

If you need to treat a plain array as a vector or matrix, no need anymore to
construct a Eigen vector or matrix, just use Eigen::Map to wrap the existing
array as an Eigen expression.

If you want to perform coefficient access without assertions (so you keep high
performance even with asserts enabled at compile-time) use coeff() and
coeffRef() instead of operator() and [].

Don’t hesitate to use Eigen xprs to shorten the code considerably – this also
gives Eigen more opportunities for optimization. Here is an example. When
porting this code from glpainter.cpp:

matrix(0, 0) = ortho1(0);
matrix(1, 0) = ortho1(1);
matrix(2, 0) = ortho1(2);
matrix(3, 0) = 0.0;

The best thing to do would be to change ortho1 to be a 4D vector, with last
coord manually set to 0, so you could then do:

matrix.col(0) = ortho1;

and automatically benefit from vectorization (since a 4D vector is perfectly
vectorizable). In fact I don’t for sure if Eigen is clever enough to
vectorize that one, but if it’s not then it’s a bug that we shall fix (this
case is definitely vectorizable as you column is column-major).

A few words about vectorization. It is important to understand that Intel CPUs
vectorization works with 128-bit packets. Consequences:

  • with floats, about the maximum speedup you can hope for is x4 (in some
    special cases like large dot products one can get slightly above that)
  • with doubles, it is x2
  • so vectorization is quite a big incentive to favor floats over doubles.
  • if the vector or matrix size is not suitable to work with 128-bit packet,
    vectorization just can’t happen. That’s the case with 3d vectors and 3x3
    matrices, unfortunately. With large enough dynamic-size matrices, Eigen will
    try to vectorize any size by finishing the last few coefficients without
    vectorization. This is why, with dynamic size, vectorization almost always
    brings a nice speedup. But for fixed size 3, this would not work too well
    (there are a few things we could still attempt, but not much). However with
    fixed-size 4x4 matrices, vectorization works really well, which is nice as
    it’s what OpenGL uses.

Cheers,
Benoit