GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-04-21 00:59:36 +08:00
Code Issues Packages Projects Releases Wiki Activity

Improve robustness of 2x2 eigenvalue with shifting and scaling

Browse Source
This commit is contained in:
Gael Guennebaud 2016-07-26 14:43:54 +02:00
parent 7f7e84aa36
commit 95113cb15c
1 changed files with 13 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
View File

@ -740,14 +740,18 @@ struct direct_selfadjoint_eigenvalues<SolverType,2,false>
EigenvectorsType& eivecs = solver.m_eivec; EigenvectorsType& eivecs = solver.m_eivec;
VectorType& eivals = solver.m_eivalues; VectorType& eivals = solver.m_eivalues;
// map the matrix coefficients to [-1:1] to avoid over- and underflow. // Shift the matrix to the mean eigenvalue and map the matrix coefficients to [-1:1] to avoid over- and underflow.
Scalar scale = mat.cwiseAbs().maxCoeff(); Scalar shift = mat.trace() / Scalar(2);
scale = numext::maxi(scale,Scalar(1)); MatrixType scaledMat = mat;
MatrixType scaledMat = mat / scale; scaledMat.coeffRef(0,1) = mat.coeff(1,0);
scaledMat.diagonal().array() -= shift;
Scalar scale = scaledMat.cwiseAbs().maxCoeff();
if(scale > Scalar(0))
scaledMat /= scale;
// Compute the eigenvalues // Compute the eigenvalues
computeRoots(scaledMat,eivals); computeRoots(scaledMat,eivals);
// compute the eigen vectors // compute the eigen vectors
if(computeEigenvectors) if(computeEigenvectors)
{ {
@ -775,10 +779,11 @@ struct direct_selfadjoint_eigenvalues<SolverType,2,false>
eivecs.col(0) << eivecs.col(1).unitOrthogonal(); eivecs.col(0) << eivecs.col(1).unitOrthogonal();
} }
} }
// Rescale back to the original size. // Rescale back to the original size.
eivals *= scale; eivals *= scale;
eivals.array() += shift;
solver.m_info = Success; solver.m_info = Success;
solver.m_isInitialized = true; solver.m_isInitialized = true;
solver.m_eigenvectorsOk = computeEigenvectors; solver.m_eigenvectorsOk = computeEigenvectors;