[eigen] Issue in SelfAdjointEigenSolver for 3-by-3 matrices

[Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx>]

Hi all,

There seems to be an issue with SelfAdjointEigenSolver when applied to 
3-by-3 matrices, because the unit test eigensolver_selfadjoint_1 fails 
occasionally on my computer.

I think I boiled it down to this example:

------------------------------------------------

#include <iostream>
#include <Eigen/Core>
#include <Eigen/Eigenvalues>

using namespace Eigen;
using namespace std;

int main()
{
  Matrix3f A = Matrix3f::Identity(3,3);
  A(0,2) = 0.001;
  SelfAdjointEigenSolver<Matrix3f> eigOfA(A);
  cout << eigOfA.eigenvalues() << endl;
}

------------------------------------------------

The output of this program is

  1
  1
  1

This program computes the eigenvalues of the matrix

  1      0     0.001
  0      1     0
  0.001  0     1

They are easily computed by hand; the exact eigenvalues are 0.999, 1, and 
1.001. I think this is rather too big an error.

The issue seems to be in

  ei_tridiagonalization_inplace_selector<MatrixType,3>::run()

in Tridiagonalization.h, where we have

    RealScalar v1norm2 = ei_abs2(mat(2,0));
    if (ei_isMuchSmallerThan(v1norm2, RealScalar(1)))
    {
	...

If the condition in the if statement is true, then mat(2,0) is simply 
ignored. This does not seem the right condition, and I'm not even sure why 
we need two branches here. Can somebody please help me here? I'd rather 
not touch code I don't understand ...

This is all old code and the same issue is present in Eigen 2, so perhaps 
we want to backport the fix before the next version is released?

Cheers,
Jitse