[eigen] Re: Eigen solver usage: simple questions

[William Oquendo <woquendo@xxxxxxxxx>]

Dear all,

It seems that .eigenvalues() assumes the matrix is selfadjoint since it is calling the self adjoint mehotd (as shown by the compiler in the assertion and after looking the source code), and of course produces bad results when compiling -DNDEBUG and disabling the assertion

Is there a way to explicitly indicates that a given matrix is not selfadjoint, and in that case to use the more general procedure ?

Best regards / Cordialmente,

--
William Oquendo
Phd Candidate
Simulation Of Physical Systems Group
National University of Colombia
Linux User # 321481
*********************
Este correo puede carecer de tildes o eñes ya que el teclado no contiene estos caracteres. Presento excusas por eso.

*********************

On Tue, Jul 20, 2010 at 4:06 PM, William Oquendo <woquendo@xxxxxxxxx> wrote:

Dear all, 

First of all I apologize in advance if my questions are too simple, I read the docs but the info about eigen solvers is small. Second, I would like to thank you for the great work in Eigen!

I have been using eigen  (currently 2.0.15, compiler g++ 4.4.4, Slackware64 13.1) to solve eigenvalue problems. But now I am concerned with the way I am using the library. I have a 3x3 array (simulating the matrix) and I need to compute, by now, only the eigenvalues. When I execute the code (please see a reduced example copied at the end and also attached), the program dies with the following error assertion (in general my matrix is not self adjoint):

a.out: /usr/local/include/eigen2/Eigen/src/QR/SelfAdjointEigenSolver.h:296: Eigen::Matrix<typename Eigen::NumTraits<typename Eigen::ei_traits<T>::Scalar>::Real, Eigen::ei_traits::ColsAtCompileTime, 1, 2, Eigen::ei_traits::ColsAtCompileTime, 1> Eigen::MatrixBase<Derived>::eigenvalues() const [with Derived = Eigen::Map<Eigen::Matrix<double, 3, 3, 2, 3, 3>, 1>]: Assertion `Flags&SelfAdjointBit' failed.

But, if I compile with the flag -DNDEBUG, the programs ends successfully with the correct data (even with null matrices), as far as I have checked. How can I specify to Eigen that my matrix is not selfAdjoint? (if that is the problem ...

Let's assume now that I need both the eigen values and eigen vectors, and I am using the full eigen solver: 

EigenSolver<Matrix3d> solver(Meigen);

Is it possible to sort the eigen values by norm and, of course, their corresponding eigen vectors by means of some internal eigen function? I can write a function to do that, but I would like to know if Eigen can do it already.

Thanks in advance for your kind help and reply.

Best regards / Cordialmente,

--
William Oquendo
Phd Candidate
Simulation Of Physical Systems Group
National University of Colombia
Linux User # 321481
*********************
Este correo puede carecer de tildes o eñes ya que el teclado no contiene estos caracteres. Presento excusas por eso.

*********************

CODE

#include <iostream>

#include <eigen2/Eigen/Core>

#include <eigen2/Eigen/Array>

#include <eigen2/Eigen/QR>

using namespace Eigen;

using namespace std;

void ComputeEigenValues(double Matrix[][3], double & l0, double & l1, double & l2);

int main()

{

  double M[3][3] = { 1, 0, 0,

    0, 1, 0,

    0, 0, 1 };

  double e0, e1, e2;

  

  ComputeEigenValues(M, e0, e1, e2);

  

  clog << "The eigen values are " << endl

       << e0 << endl

       << e1 << endl

       << e2 << endl

       << endl;

  

  return 0;

}

void ComputeEigenValues(double Matrix[][3], double & l0, double & l1, double & l2)

{

  l0 = l1 = l2 = 0;

  Map<Matrix3d> Meigen(&Matrix[0][0]);

  Vector3d veval =  Meigen.eigenvalues();

  l0 = veval..real()[0];

  l1 = veval.real()[1];

  l2 = veval.real()[2];

}