Hi,

we are switching our project to Eigen and therefore I started out comparing it with our old matrix class, which was seldon ( http://seldon.sourceforge.net/).

The first thing I noticed that Eigen takes quite some time to simply construct matrices without initializing them. I simply constructed a 30000 times 30000 matrix, which takes Eigen about 20 seconds, whereas Seldon is finished within milliseconds.
Since we fill the matrices up with parallized routines the construction time matters.
Is there a way to do a faster initialization of Eigen matrices?

Thanks in advance

Eigen::MatrixXd(rows,cols) ctor does only a call to new double[rows*cols]. Nothing else. Remark that a 30k^2 matrix requires about 7GB of memory (for doubles) ! If it takes so long then I guess that's because your system does not have such amount of memory left?

It's definitely not swapping, since the cluster I tested this on has 256 GB of memory.
I tested a little bit more and found the reason. Actually I didn't fill up the matrix with double but with complex<double>. Eigen obviously now calls new complex<double>[rows*cols] which is rather slow.
For double everything works fine and the matrix is constructed instantly.
Our old matrix class seems to do something else for complex<double>.

I see, the pb is that std::complex default ctor initialize the values to 0:

complex(const _Tp& __r = _Tp(), const _Tp& __i = _Tp())
: _M_real(__r), _M_imag(__i) { }

Your previous library probably called malloc instead of the operator new. However, some custom types really requires an initialization so we really have to use operator new. Now the question is shall we do an exception for complex<float> and complex<double>?

In the meantime, if you are using Eigen3, you can try something like (not tested):

Code: Select all

namespace Eigen { namespace internal {

template<> inline std::complex<float>* aligned_new(size_t size)
{ return reinterpret_cast<T*>(aligned_malloc(sizeof(std::complex<float>)*size)); }

template<> inline std::complex<double>* aligned_new(size_t size)
{ return reinterpret_cast<T*>(aligned_malloc(sizeof(std::complex<double>)*size)); }

}}

Thanks for your reply. You are right, my previous library used malloc. In this case Eigen3 uses constraint_aligned_new, which initializes the data with malloc and then applies the construct_elements_of_array function. I have to admit I just can guess what this functions does since I am not familiar with the syntax ::new (ptr+i). Still I was able to speed it up by defining

Code: Select all

namespace Eigen { namespace internal {

template<> inline std::complex<double>* construct_elements_of_array(std::complex<double> *ptr, size_t size)
{
 //for (size_t i=0; i < size; ++i) ::new (ptr + i) T;
  return ptr;
}
}}


This means that all complex<> valued matrices are automatically initialized with zero, whether I want it or not. On the other hand matrices with types double or float are not initialized with zero. Do I understand this correctly?
This also means the EIGEN_INITIALIZE_MATRICES_BY_ZERO just creates a useless overhead for complex<> valued matrices.

In my opinion you should create an exception for complex<float> and complex<double>, since otherwise in my opinion it's not really consistent.

yes this is exactly what I said in my previous reply, and I've also posted a message on the mailing list about this issue, You check it there : http://listengine.tuxfamily.org/lists.t ... llist.html

The archive is not updated in realtime.