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

[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]

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