| Re: [eigen] How to subtract a diagonal matrix from a matrix in one expression? |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
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. 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
expressions.
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.
Cheers,
Christoph
Best regards,
Adrien
On Thu, Jan 12, 2017 at 12:44 AM, Alexander Voigt <
alexander.voigt@xxxxxxxxxxxxxxxxxxxxx> 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.)
Example:
=====================================================
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,
Alex
--
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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
E-Mail: alexander.voigt@xxxxxxxxxxxxxxxxxxxxx
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
--
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: http://www.informatik.uni-bremen.de/robotik