Re: [eigen] How to subtract a diagonal matrix from a matrix in one expression?

[ Thread Index | Date Index | More Archives ]

On Thu, Jan 12, 2017 at 4:50 PM, Christoph Hertzberg <chtz@xxxxxxxxxxxxxxxxxxxxxxxx> wrote:
On 2017-01-11 17:51, Adrien Escande wrote:
Dear Alexander,

if you are using fixed size matrices, I don't see any problem with having
temporary objects as they don't induce any memory allocation.
The following one-liner ther work:
const Matrix2d M2 = M - Matrix2d(v.asDiagonal());
and I guess the compiler can be quite smart to avoid any overhead.

Yes, especially for small fixed sized expressions, this is often the case. Unfortunately, it is slightly sub-optimal here.

Depends a lot on the compiler, with clang 3.9 -O3:

M2 = M - Matrix2d(v.asDiagonal());
movsd (%rdx), %xmm0
movsd 8(%rdx), %xmm1
pslldq $8, %xmm1
movapd (%rsi), %xmm2
subpd %xmm0, %xmm2
movapd %xmm2, (%rdi)
movapd 16(%rsi), %xmm0
subpd %xmm1, %xmm0
movapd %xmm0, 16(%rdi)

M2 = M; M2.diagonal()-=v;
movapd (%rsi), %xmm0
movapd %xmm0, (%rdi)
movapd 16(%rsi), %xmm1
movapd %xmm1, 16(%rdi)
subsd (%rdx), %xmm0
movsd %xmm0, (%rdi)
shufpd $1, %xmm1, %xmm1
subsd 8(%rdx), %xmm1
movsd %xmm1, 24(%rdi)

Yes, clang is pretty good at understanding intrinsics and re-vectorizing them.

I was trying this:
  void foo(Matrix2d& M2, const Matrix2d& M, const Vector2d& v) {
      M2 = M - Matrix2d(v.asDiagonal());

I'm getting something like this when compiling it with -O2:
        pxor    %xmm0, %xmm0
        movaps  %xmm0, (%rsp)
        movaps  %xmm0, 16(%rsp)
        movsd   (%rdx), %xmm0
        movsd   %xmm0, (%rsp)
        movsd   8(%rdx), %xmm0
        movsd   %xmm0, 24(%rsp)
        movapd  (%rsi), %xmm0
        subpd   (%rsp), %xmm0
        movaps  %xmm0, (%rdi)
        movapd  16(%rsi), %xmm0
        subpd   16(%rsp), %xmm0
        movaps  %xmm0, 16(%rdi)

That means there is an actual temporary allocated (on the stack).

In this case, I'm getting slightly better code when compiling with -DEIGEN_DONT_VECTORIZE here:
        movsd   (%rsi), %xmm0
        movsd   8(%rdx), %xmm1
        subsd   (%rdx), %xmm0
        movsd   %xmm0, (%rdi)
        movsd   8(%rsi), %xmm0
        movsd   %xmm0, 8(%rdi)
        movsd   16(%rsi), %xmm0
        movsd   %xmm0, 16(%rdi)
        movsd   24(%rsi), %xmm0
        subsd   %xmm1, %xmm0
        movsd   %xmm0, 24(%rdi)

I'm pretty sure, if M and v were known at compile time this would entirely be evaluated at compile time.

Btw: Notice that x-0.0 can be simplified to x, whereas x+0.0 cannot, because (-0.0+0.0) == +0.0. (For any other x, including NaNs or inf, the result is x.)

If one changes the _expression_ to:
  M2 = M; M2.diagonal()-=v;
I'm getting this:
        movapd  (%rsi), %xmm0
        movaps  %xmm0, (%rdi)
        movapd  16(%rsi), %xmm0
        movaps  %xmm0, 16(%rdi)
        movsd   (%rdi), %xmm0
        subsd   (%rdx), %xmm0
        movsd   %xmm0, (%rdi)
        movsd   24(%rdi), %xmm0
        subsd   8(%rdx), %xmm0
        movsd   %xmm0, 24(%rdi)

A a side remarks, be very mindful of the auto keyword with Eigen
If you do
const auto M2 = M - Matrix2d(v.asDiagonal());
M2 will be a CwiseBinaryOp, not a Matrix2d, meaning that you subtraction
will be recomputed each time you are using M2.

Actually, it is much worse in this case:
Matrix2d(v.asDiagonal()) will be a temporary, which gets destructed after the assignment to M2. But the CwiseBinaryOp stored in M2 only holds a const reference to this temporary, therefore the _expression_ actually becomes invalid/undefined behavior (it still sometimes works, because sometimes compilers delay the actual destruction of temporaries).
Generally, one should not store the result of an Eigen operation inside an `auto` variable, unless you really want the _expression_-tree and you are sure that it only references sub-expressions which are valid as long as the _expression_ tree is used.


Best regards,

On Thu, Jan 12, 2017 at 12:44 AM, Alexander Voigt <> wrote:

Dear Eigen developers,

is there a way to subtract a diagonal matrix (which has been constructed
from a vector v using v.asDiagonal()) from a square matrix in one
_expression_?  (I'm asking for one _expression_, because I'd like to make
the result const and avoid a temporary.)


include <Eigen/Core>
using namespace Eigen;

int main() {
   Matrix<double,2,2> M;
   Matrix<double,2,1> v;

   const auto M2 = M - v.asDiagonal(); // does not compile

The marked line yields a compiler error with g++ 4.9.2.
I know that it is possible to write instead:

   Matrix<double,2,2> M2(M);
   M2 -= v.asDiagonal();     // works

However, this latter approach has the disadvantage that M2 cannot be
made const, which I'd like to achieve.  An approach that works with one
_expression_ would be

   const Matrix<double,2,2> M2 = [&M,&v]{
      Matrix<double,2,2> tmp(M);
      tmp -= v.asDiagonal();
      return tmp;

but this has the disadvantage that it is more complicated and needs a
temporary object.

Many thanks and best regards,

Dr. Alexander Voigt
Institute for Theoretical Particle Physics and Cosmology
RWTH Aachen
Sommerfeldstr. 14
52074 Aachen
Room: 28A408
Phone: +49 (0)241 80 27049
Fax  : +49 (0)241 80 22187

 Dipl. Inf., Dipl. Math. Christoph Hertzberg

 Universität Bremen
 FB 3 - Mathematik und Informatik
 AG Robotik
 Robert-Hooke-Straße 1
 28359 Bremen, Germany

 Zentrale: +49 421 178 45-6611

 Besuchsadresse der Nebengeschäftsstelle:
 Robert-Hooke-Straße 5
 28359 Bremen, Germany

 Tel.:    +49 421 178 45-4021
 Empfang: +49 421 178 45-6600
 Fax:     +49 421 178 45-4150
 E-Mail:  chtz@xxxxxxxxxxxxxxxxxxxxxxxx

 Weitere Informationen:

Mail converted by MHonArc 2.6.19+