diff options
author | Gael Guennebaud <g.gael@free.fr> | 2016-07-26 14:43:54 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2016-07-26 14:43:54 +0200 |
commit | 95113cb15cbb3efcde54bcbba59d86caa3890d66 (patch) | |
tree | 41280eee4df6be6d92a52dc14c8e72cc10b1336c /Eigen/src/Eigenvalues | |
parent | 7f7e84aa3619e5fadf9a32fddbfd6dcd17db72e9 (diff) |
Improve robustness of 2x2 eigenvalue with shifting and scaling
Diffstat (limited to 'Eigen/src/Eigenvalues')
-rw-r--r-- | Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h | 21 |
1 files changed, 13 insertions, 8 deletions
diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 9f0583fd2..8258fadbf 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -740,14 +740,18 @@ struct direct_selfadjoint_eigenvalues<SolverType,2,false> EigenvectorsType& eivecs = solver.m_eivec; VectorType& eivals = solver.m_eivalues; - // map the matrix coefficients to [-1:1] to avoid over- and underflow. - Scalar scale = mat.cwiseAbs().maxCoeff(); - scale = numext::maxi(scale,Scalar(1)); - MatrixType scaledMat = mat / scale; - + // Shift the matrix to the mean eigenvalue and map the matrix coefficients to [-1:1] to avoid over- and underflow. + Scalar shift = mat.trace() / Scalar(2); + MatrixType scaledMat = mat; + 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 computeRoots(scaledMat,eivals); - + // compute the eigen vectors if(computeEigenvectors) { @@ -775,10 +779,11 @@ struct direct_selfadjoint_eigenvalues<SolverType,2,false> eivecs.col(0) << eivecs.col(1).unitOrthogonal(); } } - + // Rescale back to the original size. eivals *= scale; - + eivals.array() += shift; + solver.m_info = Success; solver.m_isInitialized = true; solver.m_eigenvectorsOk = computeEigenvectors; |