mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-12-27 07:29:52 +08:00
Code Issues Packages Projects Releases Wiki Activity

bug #790: fix overflow in real_2x2_jacobi_svd

Browse Source
This commit is contained in:
Gael Guennebaud 2014-04-14 13:52:16 +02:00
parent 0587db8bf5
commit 148acf8e4f
1 changed files with 6 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
Eigen/src/SVD/JacobiSVD.h
View File

@ -415,6 +415,7 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q,
JacobiRotation<RealScalar> *j_right) JacobiRotation<RealScalar> *j_right)
{ {
using std::sqrt; using std::sqrt;
using std::abs;
Matrix<RealScalar,2,2> m; Matrix<RealScalar,2,2> m;
m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)), m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)),
numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q)); numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q));
@ -428,9 +429,11 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q,
} }
else else
{ {
RealScalar u = d / t; RealScalar t2d2 = numext::hypot(t,d);
rot1.c() = RealScalar(1) / sqrt(RealScalar(1) + numext::abs2(u)); rot1.c() = abs(t)/t2d2;
rot1.s() = rot1.c() * u; rot1.s() = d/t2d2;
if(t<RealScalar(0))
rot1.s() = -rot1.s();
} }
m.applyOnTheLeft(0,1,rot1); m.applyOnTheLeft(0,1,rot1);
j_right->makeJacobi(m,0,1); j_right->makeJacobi(m,0,1);