Hi,

I wrap a class around Eigen's API of SelfAdjointEigenSolver.

I.e.:

Code: Select all

   class MatrixEigenSym
   {
      public:
         typedef double TScalar;
         typedef int TIndex;
         typedef Eigen::Matrix <TScalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> Base2;
//          EIGEN_MAKE_ALIGNED_OPERATOR_NEW
      protected:
         /// the matrix object this class is operating on
         Matrix * m_object;
         /// the solver object
         Eigen::SelfAdjointEigenSolver< Base2 > m_eigen;
         /// is true if the eigenvectors were computed
         bool m_eigenvectorsAvailable;
         /// is true if the eigen problem is a general eigen problem
         bool m_isGeneralEigenproblem;
      public:
         //! destructor
         virtual ~MatrixEigenSym();
         //! constructor
         MatrixEigenSym(Matrix & _object, const bool &computevectors=true);
         MatrixEigenSym(Matrix & _object, Matrix & mass, const bool &computevectors=true);

...


Matrix is a user defined matrix type wich wraps Eigen's MatrixXd class.
Before the actual constructor of the eigen solver can be called, we must check the dimensions of the given input data.

This is done by

Code: Select all

   MatrixEigenSym::MatrixEigenSym(Matrix & _object, Matrix & mass,  const bool &computevectors)
   :
   m_object((&_object)),
   m_eigen( std::max(1,_object.rows()) ),
   m_eigenvectorsAvailable(computevectors),
   m_isGeneralEigenproblem(true)
   {
                // do all the checks here

                // factorize:
//       m_eigen = Eigen::SelfAdjointEigenSolver< Base2 >(_object,mass,computevectors);
//       m_eigen.compute(_object,mass,computevectors);
   }




What is critical, that is the line

Code: Select all

   
m_eigen( std::max(1,_object.rows()) ),


It allocates required space for the eigenvectors and eigenvalues, at least of dimension 1 (because we do not know what dimension the input marix will have).

The default constructor does not work in the moment (Assert on the dimension).

Using VALGRIND, I obtain

• "Invalid writes" when setting the values of the boolean members
• "Invalid reads" when accessing the boolean members (for example std::cout<<...)

Using EIGEN_MAKE_ALIGNED_OPERATOR_NEW does not change anything.

But removing the member

Code: Select all

Eigen::SelfAdjointEigenSolver< Base2 > m_eigen

from the class, will let the code work fine. I guess, Eigen is messing around in the memory where it shouldn't...

Thanks for any advice on that...

First of all, it would greatly help here if you could provide a self-contained test program that we can compile and run valgrind on.

From your description, nothing indicates that the invalid read/writes actually come from Eigen itself...

About the fact that the default constructor doesn't work (assert) in the dynamically sized case: this is definitely something we'll fix before the 3.0 release.

Well, I just prepared a small test program - and figured out that the bug does not appear in the test case. Must be somewhere inside of the SWIG bindings...

Anyhow, when applying the generalized eigenvalue problem, I get thousands "uninitialized value" warnings, see code below

Code: Select all

   // the basic math library used here:
   #include <Eigen/Core>
   #include <Eigen/Array> // for cwise access
   #include <Eigen/QR>

#define GLOBAL_THROW(x,m) \
   {\
   std::cerr << m <<"\n"; \
   throw std::exception(); \
   }

   class MatrixEigenSym
   {
      public:
         typedef double TScalar;
         typedef int TIndex;
         //typedef Eigen::Matrix <TScalar, Eigen::Dynamic, Eigen::Dynamic, matrix_flags> Base;
         typedef Eigen::MatrixXd Base;
         typedef Eigen::SelfAdjointEigenSolver< Base > TSolver;
//          EIGEN_MAKE_ALIGNED_OPERATOR_NEW
      protected:
         Base * m_object;
         TSolver m_eigen;
         bool m_eigenvectorsAvailable;
         bool m_isGeneralEigenproblem;
      public:
         virtual ~MatrixEigenSym();
         MatrixEigenSym(Eigen::MatrixXd & _object, const bool &computevectors=true);
         MatrixEigenSym(Eigen::MatrixXd & _object, Eigen::MatrixXd & mass, const bool &computevectors=true);
      protected:
         const Base & object() const{return *m_object;}
         Base & object(){return *m_object;}
         const Eigen::SelfAdjointEigenSolver<Base> & eigen() const {return m_eigen;}
         Eigen::SelfAdjointEigenSolver<Base> & eigen(){return m_eigen;}
         
         
   };

   
   
   

   MatrixEigenSym::~MatrixEigenSym()
   {
   }
   
   MatrixEigenSym::MatrixEigenSym(Base & _object,  const bool &computevectors)
   :
   m_object((&_object)),
   m_eigen( std::max(1,_object.rows()) ),
//    m_eigen( Eigen::MatrixXd::Zero(std::max(1,_object.rows()),std::max(1,_object.rows()) ) ),
   m_eigenvectorsAvailable(computevectors),
   m_isGeneralEigenproblem(false)
   {
      if (object().Base::rows() != object().Base::cols())
         GLOBAL_THROW      (logger_tng_tmath, "Can not perform modal analyses, square matrix expected.")
         
//       m_eigen = Eigen::SelfAdjointEigenSolver< Base2 >(_object,computevectors);
      m_eigen.compute(_object,computevectors);
//       m_eigen = new TSolver(_object,computevectors);
   }
   
   MatrixEigenSym::MatrixEigenSym(Base & _object, Base & mass,  const bool &computevectors)
   :
   m_object((&_object)),
   m_eigen( std::max(1,_object.rows()) ),
//    m_eigen( Eigen::MatrixXd::Zero(std::max(1,_object.rows()),std::max(1,_object.rows()) ) ),
   m_eigenvectorsAvailable(computevectors),
   m_isGeneralEigenproblem(true)
   {
      if (object().Base::rows() != object().Base::cols())
         GLOBAL_THROW      (logger_tng_tmath, "Can not perform modal analyses, square matrix expected.")
      if (object().Base::rows() != mass.Base::rows() || object().Base::cols() != mass.Base::cols())
         GLOBAL_THROW      (logger_tng_tmath, "Can not perform modal analyses, mass matrix has incompatible dimensions.")
         
// std::cout<<"MatrixEigenSym(Matrix & _object, Matrix & mass,  const bool &computevectors)\n";
// std::cout<<"rows: "<<_object.rows();
// std::cout<<", cols: "<<_object.cols();
// std::cout<<", vectors: "<<m_eigenvectorsAvailable;
// std::cout<<", isgeneral: "<<m_isGeneralEigenproblem<<"\n";

      m_eigenvectorsAvailable = computevectors;
      m_isGeneralEigenproblem = true;

//       m_eigen = Eigen::SelfAdjointEigenSolver< Base2 >(_object,mass,computevectors);
      m_eigen.compute(_object,mass,computevectors);
   }




int main()
{
   Eigen::MatrixXd mat(4,4);
   mat.setRandom();
   Eigen::MatrixXd mat2(4,4);
   mat2.setRandom();
   
   
   MatrixEigenSym eigen(mat, mat2, true);
//    MatrixEigenSym eigen(mat);
}


You got me thinking for a long time!

In the critical line:

Code: Select all

      m_eigen.compute(_object,mass,computevectors);


You are calling the variant to solve the generalized eigenproblem. As is explained in the API docs for that method,
http://eigen.tuxfamily.org/dox/classEig ... c4f5c2d11c
the 2nd matrix that you pass to this method is required to be positive definite. Here you are passing a non-positive matrix. Because of that, NaN values are generated.

I modified your test program to let mat2 be positive and suddenly there's no problem anymore.

I'm all for improving eigen to be more tolerant to bad input, but here, the problem happens inside of the LLT decomposition that is internally used by the generalized selfadjoint eigensolver. In that context, checking for positiveness can only be done at the cost of a runtime overhead and even so, it can't be very reliable. So I'm not sure that we should try to do something about it...