Hello Eigen community!

I am currently working on a piece of code that somebody else has written, and I have recently found a major bug in it. The code can be found here: https://github.com/alkamid/QCLD/blob/ma ... ce/QCL.cpp

I think I found the place where the bug occurs and it's line #619:

Code: Select all

PHI_VECTOR = -(z_step)*(z_step)*POISSON_M.inverse()*DENSITY_VECTOR;

If I print POISSON_M.inverse() and DENSITY_VECTOR separately, their values seem to be okay (I can check the order of magnitude if you want), but POISSON_M.inverse()*DENSITY_VECTOR just goes crazy -- its values are >10e80, when they should be more like 10e15. This completely destroys further calculations.

The strange thing is, it works just fine on Windows, and the program gives reasonable results. The problem occurs on Linux using gcc 4.8.1, as opposed to gcc 4.4.1 via MinGW on Windows. In both cases I used the latest stable Eigen, 3.2.0. Also, the multiplication in the first iteration of the main loop (the one starting at #228) seems to work fine on Linux.

I would be very, very grateful if you could point me into the right direction with this bug, as I have already spent a few days locating it. If there is anything else I can provide you with, please let me know!

I would recommend you to run your program under valgrind to detect the real memory error.

Then, if POISSON_M is not very small (>4), then doing POISSON_M.lu().solve(DENSITY_VECTOR) will be both faster and numerically more accurate than explicitly computing the inverse.

Thanks ggael. Although the program runs much faster (the matrix is >500x500), the error persists. I will try valgrind.

EDIT: valgrind gives me so many warnings that it will take at least a few days before I sort them out.

EDIT2: Okay, could you please help me start with valgrind? The very first warning I get is:

Code: Select all

==21941== Conditional jump or move depends on uninitialised value(s)
==21941==    at 0x418854: void Eigen::MatrixBase<Eigen::Block<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, -1, 1, true>, -1, 1, false> >::makeHouseholder<Eigen::VectorBlock<Eigen::Block<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, -1, 1, true>, -1, 1, false>, -1> >(Eigen::VectorBlock<Eigen::Block<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, -1, 1, true>, -1, 1, false>, -1>&, double&, double&) const (Householder.h:79)
==21941==    by 0x418CC3: Eigen::MatrixBase<Eigen::Block<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, -1, 1, true>, -1, 1, false> >::makeHouseholderInPlace(double&, double&) (Householder.h:45)
==21941==    by 0x421226: Eigen::HessenbergDecomposition<Eigen::Matrix<double, -1, -1, 0, -1, -1> >::_compute(Eigen::Matrix<double, -1, -1, 0, -1, -1>&, Eigen::Matrix<double, -1, 1, 0, -1, 1>&, Eigen::Matrix<double, 1, -1, 1, 1, -1>&) (HessenbergDecomposition.h:303)
==21941==    by 0x4236D1: Eigen::RealSchur<Eigen::Matrix<double, -1, -1, 0, -1, -1> >::compute(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, bool) (HessenbergDecomposition.h:159)
==21941==    by 0x4237E0: Eigen::EigenSolver<Eigen::Matrix<double, -1, -1, 0, -1, -1> >::compute(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, bool) (EigenSolver.h:372)
==21941==    by 0x41188C: solve_eigen(Eigen::EigenSolver<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, std::vector<layer, std::allocator<layer> >&, std::vector<double, std::allocator<double> >&, unsigned int&, unsigned int&, unsigned int&, unsigned int&, unsigned int&, double&, double&, double&, double&, double&, double&, double&, double&, double&, std::string&, std::string&, std::string&, double&, double&, unsigned int&) (EigenSolver.h:155)
==21941==    by 0x4039C0: main (QCL.cpp:247)
==21941==  Uninitialised value was created by a heap allocation
==21941==    at 0x4C2A2DB: malloc (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)
==21941==    by 0x415483: Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >::resize(long, long) (Memory.h:220)
==21941==    by 0x4114DC: solve_eigen(Eigen::EigenSolver<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, std::vector<layer, std::allocator<layer> >&, std::vector<double, std::allocator<double> >&, unsigned int&, unsigned int&, unsigned int&, unsigned int&, unsigned int&, double&, double&, double&, double&, double&, double&, double&, double&, double&, std::string&, std::string&, std::string&, double&, double&, unsigned int&) (PlainObjectBase.h:627)
==21941==    by 0x4039C0: main (QCL.cpp:247)


Line 247 is:

Code: Select all

solve_eigen         ( H , wafer, phi, VARS );

and because valgrind is complaining about a Matrix object, I suspect something is wrong with H initialisation:

Code: Select all

EigenSolver<MatrixXd>  H (Eigen_dim);

Is it not how I am supposed to initialise it? Can you see what is wrong here?

hm, could you show more code... What solve_eigen is? When you do EigenSolver<MatrixXd> H (Eigen_dim);, this creates an eigenvalue solveur and preallocate memory for matrices of size Eigen_dim x Eigen_dim. However, before using it, you must call H.compute(A); with some matrix A.

ggael wrote:hm, could you show more code... What solve_eigen is?

The whole code is here: https://github.com/alkamid/QCLD/blob/ma ... ce/QCL.cpp

instead of:

H = EigenSolver<MatrixXd> ( H_INITIALISE ,true )

you should do:

H.compute( H_INITIALISE ,true );

maybe there are more issues though...

Also make sure solve_eigen is not inlined to help debugging with valgrind.

It seems to work now, thanks. But what was the issue, what I was doing wrong?

Good question actually. Maybe this is hiding another memory issue (is valgrind completely green now?) or maybe your using an old version of Eigen which had and issue regarding that? The following example is running fine with valgrind:

Code: Select all

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;

int n = 100;

EIGEN_DONT_INLINE void foo(EigenSolver<MatrixXd>& H) {
  MatrixXd M(n,n);
  M.setRandom();
  H = EigenSolver<MatrixXd>(M,true);
}

int main() {
  EigenSolver<MatrixXd> eig0(n);
  foo(eig0);
  std::cout << eig0.eigenvalues() << "\n";
  std::cout << eig0.eigenvectors() << "\n";
}


hm wait, the problem is that not all entries in H_INITIALISE are set. You should call:

H_INITIALISE.setZero();

before starting filling it. By default Eigen's matrices are not intialized to zero, they are not initialized at all.

No, valgrind is far from being green, but at least it got a few more iterations right.

The first valgrind error (warning?) looks exactly the same to the one I posted in my second post, but I have no idea what is wrong there. I am compiling with -O2 flag, so all should be inline.

ggael wrote:hm wait, the problem is that not all entries in H_INITIALISE are set. You should call:

H_INITIALISE.setZero();

before starting filling it. By default Eigen's matrices are not intialized to zero, they are not initialized at all.

Then the program slows down very significantly, but I understand what you mean. It is a tridiagonal matrix -- is there any way that Eigen can solve it faster?

adding .setZero() should not be slower since you are using a dense matrix anyway ! Is your matrix symmetric? If so you could EigenSelfadjointSolver and call computeFromTridiagonal(...).

Oh, I just panicked about the execution time -- it is in fact valgrind that takes so long, not the program itself. Actually, you solved two major issues with my program, thanks a lot for that! Now it seems to run just fine.

The matrix is not symmetric, unfortunately, so I guess I have to use dense matrices.

Now when I run valgrind the first error is:

Code: Select all

==21725== Conditional jump or move depends on uninitialised value(s)
==21725==    at 0x41DE7D: Eigen::internal::partial_lu_impl<double, 0, int>::unblocked_lu(Eigen::Block<Eigen::Map<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 0, Eigen::Stride<0, 0> >, -1, -1, false>&, int*, int&) (Visitor.h:149)
==21725==    by 0x41F0CD: Eigen::internal::partial_lu_impl<double, 0, int>::blocked_lu(long, long, double*, long, int*, int&, long) (PartialPivLU.h:308)
==21725==    by 0x41EA13: Eigen::internal::partial_lu_impl<double, 0, int>::blocked_lu(long, long, double*, long, int*, int&, long) (PartialPivLU.h:343)
==21725==    by 0x41EA13: Eigen::internal::partial_lu_impl<double, 0, int>::blocked_lu(long, long, double*, long, int*, int&, long) (PartialPivLU.h:343)
==21725==    by 0x41FAEE: Eigen::PartialPivLU<Eigen::Matrix<double, -1, -1, 0, -1, -1> >::compute(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&) (PartialPivLU.h:379)
==21725==    by 0x41FE02: Eigen::PartialPivLU<Eigen::Matrix<double, -1, -1, 0, -1, -1> >::PartialPivLU(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&) (PartialPivLU.h:209)
==21725==    by 0x41FF25: Eigen::Matrix<double, -1, -1, 0, -1, -1>::Matrix<Eigen::internal::inverse_impl<Eigen::Matrix<double, -1, -1, 0, -1, -1> > >(Eigen::ReturnByValue<Eigen::internal::inverse_impl<Eigen::Matrix<double, -1, -1, 0, -1, -1> > > const&) (PartialPivLU.h:479)
==21725==    by 0x408FD0: solve_poisson(double, double, Eigen::EigenSolver<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, std::vector<layer, std::allocator<layer> >&, std::vector<double, std::allocator<double> >&, std::vector<double, std::allocator<double> >&, int, unsigned int&, unsigned int&, unsigned int&, unsigned int&, unsigned int&, double&, double&, double&, double&, double&, double&, double&, double&, double&, std::string&, std::string&, std::string&, double&, double&, unsigned int&) (CwiseUnaryOp.h:68)
==21725==    by 0x40FFF6: main (QCL.cpp:246)
==21725==  Uninitialised value was created by a heap allocation
==21725==    at 0x4C274A8: malloc (vg_replace_malloc.c:236)
==21725==    by 0x41F867: Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >::resize(long, long) (Memory.h:220)
==21725==    by 0x407BD0: solve_poisson(double, double, Eigen::EigenSolver<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, std::vector<layer, std::allocator<layer> >&, std::vector<double, std::allocator<double> >&, std::vector<double, std::allocator<double> >&, int, unsigned int&, unsigned int&, unsigned int&, unsigned int&, unsigned int&, double&, double&, double&, double&, double&, double&, double&, double&, double&, std::string&, std::string&, std::string&, double&, double&, unsigned int&) (PlainObjectBase.h:627)
==21725==    by 0x40FFF6: main (QCL.cpp:246)


Line #246 is:

Code: Select all

solve_poisson(wafer[0].get_z_min(),wafer[layer_no - 1].get_z_max(),H,wafer,phi,phi_old2,iter, VARS);

It's probably the same. You should check that you're not assuming a matrix is automatically initialized by zero.

We plane to add a computeFromHessenberg shortcut in EigenSolver. That should save you some times.

Maybe the problem is in "H"? This is how I initialise it: (line 177)

Code: Select all

    EigenSolver<MatrixXd>  H (Eigen_dim);

But is this the right way?