Consider the following:

Code: Select all

   Eigen::Matrix<double,3,3> M;
   M << 1,2,3,4,5,6,7,8,9;
   Eigen::Matrix<double,3,1> v1, v2;
   v1 << 1,2,3;
   v2 << 3,4,5;
   const Eigen::Matrix<double,3,1> & N (M*(-(v1-v2)));
   std::cout << N << "\n";


This particular piece of code works fine, but in a more complex context matrix N appears to be non-initialized. If instead of

Code: Select all

   const Eigen::Matrix<double,3,1> & N (M*(-(v1-v2)));


I put

Code: Select all

   const auto & N (M*(-(v1-v2)));


or

Code: Select all

   const Eigen::Matrix<double,3,1> N (M*(-(v1-v2)));


then it works. What is the problem here? Is Eigen unable to handle constant references to temporaries? Is it guaranteed to work with auto, or do I have to copy the temporary?

Hi vitke,

an exact example which actually causes an issue would be great. At the moment, I can only suggest to remove the reference and refrain from using auto.

You should in general refrain from using auto in context with Eigen. Almost every operation implemented in Eigen results in a template expression and the type behind your auto declared variable will be most surprising if you look at it. The expression templates also do not evaluate anything until they are assigned to a Matrix. The std::cout works with auto declared types because the implementation of the streaming operator calls .eval() on the dense expression you are passing which creates a temporary object which is finally evaluated.

Regards,
Hauke

Hi Hauke, thanks for your reply. I cannot reproduce this in an isolated example and the whole code is 40,000+ lines, so I can't post an exact example.

This is probably related: I have two matrices R and psi of type Eigen::Matrix<double,3,3> and I do

Code: Select all

         const auto & M = R.transpose() * psi * R;


If I print M it is uninitialized, if I run Valgrind it reports that M is uninitialized. If I write instead

Code: Select all

         const Eigen::Matrix<double,3,3> & M = R.transpose() * psi * R;


then M is printed just fine and Valgrind does not complain.

I understand that the exact type of a product is compicated, but uninitialized memory!? This does not sound like a template problem. And in the first example I had an opposite case - it was fine with auto and it wasn't fine if I specified the type. I can't see how could this possibly matter for memory initialization. If I print the same product in the following line

Code: Select all

         const auto & M = R.transpose() * psi * R;
         std::cout << M << " " << R.transpose() * psi * R << "\n";


then everything is printed fine and there are no uninitialized numbers. This is most definitely a memory problem.

Hi again,

the only suggestion I can give you then for now is once again to remove the reference (ampersamp) on your variable declarations. It makes no sense in the context you provided and is the only part of your code which looks a bit dodgy.

Regards,
Hauke

Curiously, this happens only with g++-4.8.1. If I use g++-4.7.3 then it does not happen. Here is what gcc 4.8 release notes say:

Several scalability bottle-necks have been removed from GCC's optimization passes. Compilation of extremely large functions, e.g. due to the use of the flatten attribute in the "Eigen" C++ linear algebra templates library, is significantly faster than previous releases of GCC.

Could these changes have introduced a bug in gcc? Is there a way to turn these optimizations off? My functions are not longer than a few hundred lines, but they are mostly inline, so I suppose that at some point in the compilation process a single huge function is created.

Not realated at all, and the fact it's "working" with gcc 4.7 is just luck. I'm pretty sure valgrind will detect memory errors even with gcc 4.7. Your references makes no sense, recall that the type of "A*B" is a abstract representation of a matrix product which basically only stores references to A and B. The fact your code with references is compiling is due to internal implementation details. It won't compile with A+B.

Here is an isolated example:

Code: Select all

#include <Dense>
#include <iostream>

template<int i>
struct A
{
   inline void f() const
   {
      A<i-1> a0, a1;
      a0.f();
      a1.f();
      Eigen::Matrix<double,3,3> M, R;
      M << 1,2,3,4,5,6,7,8,9;
      R << 9,8,7,6,5,4,3,2,1;
      const Eigen::Matrix<double,3,3> & N = M * R;
      std::cout << i << " " << N << "\n";
   }
};

template<>
struct A<0>
{
   inline void f() const {}
};

int main()
{
   A<2> a;
   a.f();
}

I compile it with

Code: Select all

g++-4.8 -Wall -I/usr/local/src/eigen/Eigen/ -O2 -o testEigen testEigen.cc

The output is

1 6.9532e-310 6.94606e-310 6.94605e-310
6.9532e-310 6.94593e-310 6.94606e-310
0 0 6.94606e-310
1 6.9532e-310 6.94606e-310 6.94605e-310
6.9532e-310 6.94593e-310 6.94606e-310
0 0 6.94606e-310
2 3.11661e-317 6.94605e-310 3.11621e-317
2.07452e-317 4.94066e-324 2.07743e-317
3.11703e-317 9.88131e-324 6.94606e-310

The optimization is essential: with -O0 it's all right. Something similar also happens with g++ 4.7.3, but not with g++ 4.6.4. Valgrind gives many errors with g++ 4.8.1 and with g++ 4.7.3, but no errors with g++ 4.6.4. I tried it on three machines running Ubuntu and on one machine running Arch, with at least two different and not too old versions of Eigen, including a version two days old. Valgrind outputs are too long to be pasted here, and I can't find a way to attach them, but you will probably be able to generate them.

Another example, with no reference: just replace

Code: Select all

const Eigen::Matrix<double,3,3> & N = M * R;


with

Code: Select all

const auto N = R.transpose() * M * R;


and add the compiler option -std=c++0x. Similar output is obtained.

And let me just add: Valgrind does not detect memory errors, but it it repeats this error over and over:

Code: Select all

==6832== Conditional jump or move depends on uninitialised value(s)
==6832==    at 0x539A683: __printf_fp (printf_fp.c:406)
==6832==    by 0x53975B7: vfprintf (vfprintf.c:1629)
==6832==    by 0x53BF441: vsnprintf (vsnprintf.c:120)
==6832==    by 0x4EB648F: ??? (in /usr/lib/x86_64-linux-gnu/libstdc++.so.6.0.18)
==6832==    by 0x4EBC8C5: std::ostreambuf_iterator<char, std::char_traits<char> > std::num_put<char, std::ostreambuf_iterator<char, std::char_traits<char> > >::_M_insert_float<double>(std::ostreambuf_iterator<char, std::char_traits<char> >, std::ios_base&, char, char, double) const (in /usr/lib/x86_64-linux-gnu/libstdc++.so.6.0.18)
==6832==    by 0x4EBCBAF: std::num_put<char, std::ostreambuf_iterator<char, std::char_traits<char> > >::do_put(std::ostreambuf_iterator<char, std::char_traits<char> >, std::ios_base&, char, double) const (in /usr/lib/x86_64-linux-gnu/libstdc++.so.6.0.18)
==6832==    by 0x4EC7EC4: std::ostream& std::ostream::_M_insert<double>(double) (in /usr/lib/x86_64-linux-gnu/libstdc++.so.6.0.18)
==6832==    by 0x401A9B: std::ostream& Eigen::internal::print_matrix<Eigen::Matrix<double, 3, 3, 0, 3, 3> >(std::ostream&, Eigen::Matrix<double, 3, 3, 0, 3, 3> const&, Eigen::IOFormat const&) (ostream:221)
==6832==    by 0x40213F: std::ostream& Eigen::operator<< <Eigen::Matrix<double, 3, 3, 0, 3, 3> >(std::ostream&, Eigen::DenseBase<Eigen::Matrix<double, 3, 3, 0, 3, 3> > const&) (IO.h:244)
==6832==    by 0x402815: A<1>::f() const (testEigen.cc:16)
==6832==    by 0x40283E: A<2>::f() const (testEigen.cc:10)
==6832==    by 0x40135D: main (testEigen.cc:29)


And finally I understood what you were saying about references I understand that it makes no sense to write auto & A = B*C in Eigen context, but why it doesn't make sense to write Matrix<double,3,3> & A = B*C? If I ommit the reference, then the product is evaluated i.e. converted to Matrix and then copied to A. If I put the reference, then copy operation is not performed. Am I right?

You're wrong. Matrix<double,3,3> & A = B*C is not correct because B*C is not a Matrix<double,3,3>. The fact it's compiling is because, internally, B*C stores a Matrix<double,3,3> in the case it is plugged into a more complex expression and needs to be evaluated into a temporary. So at the end, A is storing a reference to an object that is deleted right after this statement, and A becomes invalid. The correct way is to write:

Matrix<double,3,3> A = B*C ;

and thanks to expression templates, there is no temporary or extra copy at all.

OK, I will omit the reference. But I think such behavior violates the standard (§12.2):

4 There are two contexts in which temporaries are destroyed at a different point than the end of the full-
expression. The first context is when a default constructor is called to initialize an element of an array. If
the constructor has one or more default arguments, the destruction of every temporary created in a default
argument is sequenced before the construction of the next array element, if any.
5 The second context is when a reference is bound to a temporary. The temporary to which the reference is
bound or the temporary that is the complete object of a subobject to which the reference is bound persists
for the lifetime of the reference except:
— A temporary bound to a reference member in a constructor’s ctor-initializer (12.6.2) persists until the
constructor exits.
— A temporary bound to a reference parameter in a function call (5.2.2) persists until the completion of
the full-expression containing the call

Thus the lifetime of a temporary to which a reference is bound should be extended for the lifetime of that reference.

The standard gives the following example:

Code: Select all

struct S {
   S();
   S(int);
   friend S operator+(const S&, const S&);
   ~S();
};
S obj1;
const S& cr = S(16)+S(23);
S obj2

and later it says

The temporary T3 bound to the reference cr is destroyed at the end of cr’s lifetime, that is, at the end of the program

So what is the difference with that case and our case? A constant reference is bound to an expression result in both cases.

I think I understand: the standard says

The lifetime of a temporary bound to the returned value in a function return statement (6.6.3) is not
extended; the temporary is destroyed at the end of the full-expression in the return statement.

So since the Eigen matrix expression result is not a matrix, it has to be converted by implicitly calling the eval function, the result of which is a temporary to which I bind my reference. But this temporary is destroyed. It does not violate the standard, although it violates my intuition. Thanks for opening my eyes!

Thus the lifetime of a temporary to which a reference is bound should be extended for the lifetime of that reference.

hm, ok so it might still persists but it won't be initialized anyway. This is coherent with what valgrind told you. In the future this won't compile at all.