From 3eb5253ca1352391f4ea31c4a2dd06c34c4a33e7 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 2 Sep 2014 14:41:14 +0200 Subject: Optimization: "matrix * real" did not call the special path and the real was converted to a complex. Add respective unit test to avoid future regression. --- test/linearstructure.cpp | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) (limited to 'test') diff --git a/test/linearstructure.cpp b/test/linearstructure.cpp index 618984d5c..b627915ce 100644 --- a/test/linearstructure.cpp +++ b/test/linearstructure.cpp @@ -2,11 +2,16 @@ // for linear algebra. // // Copyright (C) 2006-2008 Benoit Jacob +// Copyright (C) 2014 Gael Guennebaud // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. +static bool g_called; + +#define EIGEN_SPECIAL_SCALAR_MULTIPLE_PLUGIN { g_called = true; } + #include "main.h" template void linearStructure(const MatrixType& m) @@ -68,6 +73,24 @@ template void linearStructure(const MatrixType& m) VERIFY_IS_APPROX(m1.block(0,0,rows,cols) * s1, m1 * s1); } +// Make sure that complex * real and real * complex are properly optimized +template void real_complex(DenseIndex rows = MatrixType::RowsAtCompileTime, DenseIndex cols = MatrixType::ColsAtCompileTime) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + + RealScalar s = internal::random(); + MatrixType m1 = MatrixType::Random(rows, cols); + + g_called = false; + VERIFY_IS_APPROX(s*m1, Scalar(s)*m1); + VERIFY(g_called && "real * matrix not properly optimized"); + + g_called = false; + VERIFY_IS_APPROX(m1*s, m1*Scalar(s)); + VERIFY(g_called && "matrix * real not properly optimized"); +} + void test_linearstructure() { for(int i = 0; i < g_repeat; i++) { @@ -80,5 +103,8 @@ void test_linearstructure() CALL_SUBTEST_7( linearStructure(MatrixXi (internal::random(1,EIGEN_TEST_MAX_SIZE), internal::random(1,EIGEN_TEST_MAX_SIZE))) ); CALL_SUBTEST_8( linearStructure(MatrixXcd(internal::random(1,EIGEN_TEST_MAX_SIZE/2), internal::random(1,EIGEN_TEST_MAX_SIZE/2))) ); CALL_SUBTEST_9( linearStructure(ArrayXXf (internal::random(1,EIGEN_TEST_MAX_SIZE), internal::random(1,EIGEN_TEST_MAX_SIZE))) ); + + CALL_SUBTEST_10( real_complex() ); + CALL_SUBTEST_10( real_complex(10,10) ); } } -- cgit v1.2.3 From 18fbe7e7d4cdb737ef5775bbb32fe62b6f8ef70e Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 2 Sep 2014 14:49:23 +0200 Subject: Fix stableNorm() with respect to NaN and inf, and add respective unit tests. blueNorm() and hypotNorm() are broken wrt to NaN/inf --- Eigen/src/Core/StableNorm.h | 11 +++++++- test/stable_norm.cpp | 69 ++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 78 insertions(+), 2 deletions(-) (limited to 'test') diff --git a/Eigen/src/Core/StableNorm.h b/Eigen/src/Core/StableNorm.h index 525620bad..07a021505 100644 --- a/Eigen/src/Core/StableNorm.h +++ b/Eigen/src/Core/StableNorm.h @@ -20,7 +20,7 @@ inline void stable_norm_kernel(const ExpressionType& bl, Scalar& ssq, Scalar& sc using std::max; Scalar maxCoeff = bl.cwiseAbs().maxCoeff(); - if (maxCoeff>scale) + if(maxCoeff>scale) { ssq = ssq * numext::abs2(scale/maxCoeff); Scalar tmp = Scalar(1)/maxCoeff; @@ -29,12 +29,21 @@ inline void stable_norm_kernel(const ExpressionType& bl, Scalar& ssq, Scalar& sc invScale = NumTraits::highest(); scale = Scalar(1)/invScale; } + else if(maxCoeff>NumTraits::highest()) // we got a INF + { + invScale = Scalar(1); + scale = maxCoeff; + } else { scale = maxCoeff; invScale = tmp; } } + else if(maxCoeff!=maxCoeff) // we got a NaN + { + scale = maxCoeff; + } // TODO if the maxCoeff is much much smaller than the current scale, // then we can neglect this sub vector diff --git a/test/stable_norm.cpp b/test/stable_norm.cpp index 549f91fbf..88cab7aa3 100644 --- a/test/stable_norm.cpp +++ b/test/stable_norm.cpp @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2009 Gael Guennebaud +// Copyright (C) 2009-2014 Gael Guennebaud // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed @@ -14,6 +14,21 @@ template bool isNotNaN(const T& x) return x==x; } +template bool isNaN(const T& x) +{ + return x!=x; +} + +template bool isInf(const T& x) +{ + return x > NumTraits::highest(); +} + +template bool isMinusInf(const T& x) +{ + return x < NumTraits::lowest(); +} + // workaround aggressive optimization in ICC template EIGEN_DONT_INLINE T sub(T a, T b) { return a - b; } @@ -106,6 +121,58 @@ template void stable_norm(const MatrixType& m) VERIFY_IS_APPROX(vrand.rowwise().stableNorm(), vrand.rowwise().norm()); VERIFY_IS_APPROX(vrand.rowwise().blueNorm(), vrand.rowwise().norm()); VERIFY_IS_APPROX(vrand.rowwise().hypotNorm(), vrand.rowwise().norm()); + + // test NaN, +inf, -inf + MatrixType v; + Index i = internal::random(0,rows-1); + Index j = internal::random(0,cols-1); + + // NaN + { + v = vrand; + v(i,j) = RealScalar(0)/RealScalar(0); + VERIFY(!isFinite(v.squaredNorm())); VERIFY(isNaN(v.squaredNorm())); + VERIFY(!isFinite(v.norm())); VERIFY(isNaN(v.norm())); + VERIFY(!isFinite(v.stableNorm())); VERIFY(isNaN(v.stableNorm())); + VERIFY(!isFinite(v.blueNorm())); VERIFY(isNaN(v.blueNorm())); +// VERIFY(!isFinite(v.hypotNorm())); //VERIFY(isNaN(v.hypotNorm())); + } + + // +inf + { + v = vrand; + v(i,j) = RealScalar(1)/RealScalar(0); + VERIFY(!isFinite(v.squaredNorm())); VERIFY(isInf(v.squaredNorm())); + VERIFY(!isFinite(v.norm())); VERIFY(isInf(v.norm())); + VERIFY(!isFinite(v.stableNorm())); VERIFY(isInf(v.stableNorm())); +// VERIFY(!isFinite(v.blueNorm())); //VERIFY(isInf(v.blueNorm())); +// VERIFY(!isFinite(v.hypotNorm())); //VERIFY(isInf(v.hypotNorm())); + } + + // -inf + { + v = vrand; + v(i,j) = RealScalar(-1)/RealScalar(0); + VERIFY(!isFinite(v.squaredNorm())); VERIFY(isInf(v.squaredNorm())); + VERIFY(!isFinite(v.norm())); VERIFY(isInf(v.norm())); + VERIFY(!isFinite(v.stableNorm())); VERIFY(isInf(v.stableNorm())); +// VERIFY(!isFinite(v.blueNorm())); VERIFY(isInf(v.blueNorm())); +// VERIFY(!isFinite(v.hypotNorm())); VERIFY(isInf(v.hypotNorm())); + } + + // mix + { + Index i2 = internal::random(0,rows-1); + Index j2 = internal::random(0,cols-1); + v = vrand; + v(i,j) = RealScalar(-1)/RealScalar(0); + v(i2,j2) = RealScalar(0)/RealScalar(0); + VERIFY(!isFinite(v.squaredNorm())); VERIFY(isNaN(v.squaredNorm())); + VERIFY(!isFinite(v.norm())); VERIFY(isNaN(v.norm())); + VERIFY(!isFinite(v.stableNorm())); VERIFY(isNaN(v.stableNorm())); +// VERIFY(!isFinite(v.blueNorm())); //VERIFY(isNaN(v.blueNorm())); +// VERIFY(!isFinite(v.hypotNorm())); VERIFY(isNaN(v.hypotNorm())); + } } void test_stable_norm() -- cgit v1.2.3 From a44a343f034eb915f9ad6dcd69e4be534d48bbe5 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 2 Sep 2014 15:06:24 +0200 Subject: Fix blueNorm wrt NaN/INF. --- Eigen/src/Core/StableNorm.h | 11 +++++------ test/stable_norm.cpp | 6 +++--- 2 files changed, 8 insertions(+), 9 deletions(-) (limited to 'test') diff --git a/Eigen/src/Core/StableNorm.h b/Eigen/src/Core/StableNorm.h index 07a021505..64d43e1b1 100644 --- a/Eigen/src/Core/StableNorm.h +++ b/Eigen/src/Core/StableNorm.h @@ -64,7 +64,7 @@ blueNorm_impl(const EigenBase& _vec) using std::abs; const Derived& vec(_vec.derived()); static bool initialized = false; - static RealScalar b1, b2, s1m, s2m, overfl, rbig, relerr; + static RealScalar b1, b2, s1m, s2m, rbig, relerr; if(!initialized) { int ibeta, it, iemin, iemax, iexp; @@ -93,7 +93,6 @@ blueNorm_impl(const EigenBase& _vec) iexp = - ((iemax+it)/2); s2m = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp))); // scaling factor for upper range - overfl = rbig*s2m; // overflow boundary for abig eps = RealScalar(pow(double(ibeta), 1-it)); relerr = sqrt(eps); // tolerance for neglecting asml initialized = true; @@ -110,13 +109,13 @@ blueNorm_impl(const EigenBase& _vec) else if(ax < b1) asml += numext::abs2(ax*s1m); else amed += numext::abs2(ax); } + if(amed!=amed) + return amed; // we got a NaN if(abig > RealScalar(0)) { abig = sqrt(abig); - if(abig > overfl) - { - return rbig; - } + if(abig > rbig) // overflow, or *this contains INF values + return abig; // return INF if(amed > RealScalar(0)) { abig = abig/s2m; diff --git a/test/stable_norm.cpp b/test/stable_norm.cpp index 88cab7aa3..119b5b424 100644 --- a/test/stable_norm.cpp +++ b/test/stable_norm.cpp @@ -145,7 +145,7 @@ template void stable_norm(const MatrixType& m) VERIFY(!isFinite(v.squaredNorm())); VERIFY(isInf(v.squaredNorm())); VERIFY(!isFinite(v.norm())); VERIFY(isInf(v.norm())); VERIFY(!isFinite(v.stableNorm())); VERIFY(isInf(v.stableNorm())); -// VERIFY(!isFinite(v.blueNorm())); //VERIFY(isInf(v.blueNorm())); + VERIFY(!isFinite(v.blueNorm())); VERIFY(isInf(v.blueNorm())); // VERIFY(!isFinite(v.hypotNorm())); //VERIFY(isInf(v.hypotNorm())); } @@ -156,7 +156,7 @@ template void stable_norm(const MatrixType& m) VERIFY(!isFinite(v.squaredNorm())); VERIFY(isInf(v.squaredNorm())); VERIFY(!isFinite(v.norm())); VERIFY(isInf(v.norm())); VERIFY(!isFinite(v.stableNorm())); VERIFY(isInf(v.stableNorm())); -// VERIFY(!isFinite(v.blueNorm())); VERIFY(isInf(v.blueNorm())); + VERIFY(!isFinite(v.blueNorm())); VERIFY(isInf(v.blueNorm())); // VERIFY(!isFinite(v.hypotNorm())); VERIFY(isInf(v.hypotNorm())); } @@ -170,7 +170,7 @@ template void stable_norm(const MatrixType& m) VERIFY(!isFinite(v.squaredNorm())); VERIFY(isNaN(v.squaredNorm())); VERIFY(!isFinite(v.norm())); VERIFY(isNaN(v.norm())); VERIFY(!isFinite(v.stableNorm())); VERIFY(isNaN(v.stableNorm())); -// VERIFY(!isFinite(v.blueNorm())); //VERIFY(isNaN(v.blueNorm())); + VERIFY(!isFinite(v.blueNorm())); VERIFY(isNaN(v.blueNorm())); // VERIFY(!isFinite(v.hypotNorm())); VERIFY(isNaN(v.hypotNorm())); } } -- cgit v1.2.3 From 42e27d41a2014043799b44b68b0d5644629279ba Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 2 Sep 2014 16:09:39 +0200 Subject: Fix hypot() and hypotNorm() wrt NaN and INF values. --- Eigen/src/Core/MathFunctions.h | 15 +++++++++++---- Eigen/src/Core/functors/BinaryFunctors.h | 14 +++++++++++--- test/stable_norm.cpp | 8 ++++---- 3 files changed, 26 insertions(+), 11 deletions(-) (limited to 'test') diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index 20fc2be74..5df1fbfec 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h @@ -308,10 +308,17 @@ struct hypot_impl using std::sqrt; RealScalar _x = abs(x); RealScalar _y = abs(y); - RealScalar p = (max)(_x, _y); - if(p==RealScalar(0)) return 0; - RealScalar q = (min)(_x, _y); - RealScalar qp = q/p; + Scalar p, qp; + if(_x>_y) + { + p = _x; + qp = _y / p; + } + else + { + p = _y; + qp = _x / p; + } return p * sqrt(RealScalar(1) + qp*qp); } }; diff --git a/Eigen/src/Core/functors/BinaryFunctors.h b/Eigen/src/Core/functors/BinaryFunctors.h index ba094f5d1..157d075a7 100644 --- a/Eigen/src/Core/functors/BinaryFunctors.h +++ b/Eigen/src/Core/functors/BinaryFunctors.h @@ -167,9 +167,17 @@ template struct scalar_hypot_op { EIGEN_USING_STD_MATH(max); EIGEN_USING_STD_MATH(min); using std::sqrt; - Scalar p = (max)(_x, _y); - Scalar q = (min)(_x, _y); - Scalar qp = q/p; + Scalar p, qp; + if(_x>_y) + { + p = _x; + qp = _y / p; + } + else + { + p = _y; + qp = _x / p; + } return p * sqrt(Scalar(1) + qp*qp); } }; diff --git a/test/stable_norm.cpp b/test/stable_norm.cpp index 119b5b424..f76919af4 100644 --- a/test/stable_norm.cpp +++ b/test/stable_norm.cpp @@ -135,7 +135,7 @@ template void stable_norm(const MatrixType& m) VERIFY(!isFinite(v.norm())); VERIFY(isNaN(v.norm())); VERIFY(!isFinite(v.stableNorm())); VERIFY(isNaN(v.stableNorm())); VERIFY(!isFinite(v.blueNorm())); VERIFY(isNaN(v.blueNorm())); -// VERIFY(!isFinite(v.hypotNorm())); //VERIFY(isNaN(v.hypotNorm())); + VERIFY(!isFinite(v.hypotNorm())); VERIFY(isNaN(v.hypotNorm())); } // +inf @@ -146,7 +146,7 @@ template void stable_norm(const MatrixType& m) VERIFY(!isFinite(v.norm())); VERIFY(isInf(v.norm())); VERIFY(!isFinite(v.stableNorm())); VERIFY(isInf(v.stableNorm())); VERIFY(!isFinite(v.blueNorm())); VERIFY(isInf(v.blueNorm())); -// VERIFY(!isFinite(v.hypotNorm())); //VERIFY(isInf(v.hypotNorm())); + VERIFY(!isFinite(v.hypotNorm())); VERIFY(isInf(v.hypotNorm())); } // -inf @@ -157,7 +157,7 @@ template void stable_norm(const MatrixType& m) VERIFY(!isFinite(v.norm())); VERIFY(isInf(v.norm())); VERIFY(!isFinite(v.stableNorm())); VERIFY(isInf(v.stableNorm())); VERIFY(!isFinite(v.blueNorm())); VERIFY(isInf(v.blueNorm())); -// VERIFY(!isFinite(v.hypotNorm())); VERIFY(isInf(v.hypotNorm())); + VERIFY(!isFinite(v.hypotNorm())); VERIFY(isInf(v.hypotNorm())); } // mix @@ -171,7 +171,7 @@ template void stable_norm(const MatrixType& m) VERIFY(!isFinite(v.norm())); VERIFY(isNaN(v.norm())); VERIFY(!isFinite(v.stableNorm())); VERIFY(isNaN(v.stableNorm())); VERIFY(!isFinite(v.blueNorm())); VERIFY(isNaN(v.blueNorm())); -// VERIFY(!isFinite(v.hypotNorm())); VERIFY(isNaN(v.hypotNorm())); + VERIFY(!isFinite(v.hypotNorm())); VERIFY(isNaN(v.hypotNorm())); } } -- cgit v1.2.3 From 51b3f558bb76c11149fc64971db786798f1b692c Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 8 Sep 2014 10:21:22 +0200 Subject: Fix bug #822: outer products needed linear access, and add respective unit tests --- Eigen/src/Core/GeneralProduct.h | 4 ++-- test/product.h | 8 ++++++++ 2 files changed, 10 insertions(+), 2 deletions(-) (limited to 'test') diff --git a/Eigen/src/Core/GeneralProduct.h b/Eigen/src/Core/GeneralProduct.h index 624b8b6e8..7179eb124 100644 --- a/Eigen/src/Core/GeneralProduct.h +++ b/Eigen/src/Core/GeneralProduct.h @@ -231,7 +231,7 @@ EIGEN_DONT_INLINE void outer_product_selector_run(const ProductType& prod, Dest& // FIXME not very good if rhs is real and lhs complex while alpha is real too const Index cols = dest.cols(); for (Index j=0; j diff --git a/test/product.h b/test/product.h index 856b234ac..0b3abe402 100644 --- a/test/product.h +++ b/test/product.h @@ -139,4 +139,12 @@ template void product(const MatrixType& m) // inner product Scalar x = square2.row(c) * square2.col(c2); VERIFY_IS_APPROX(x, square2.row(c).transpose().cwiseProduct(square2.col(c2)).sum()); + + // outer product + VERIFY_IS_APPROX(m1.col(c) * m1.row(r), m1.block(0,c,rows,1) * m1.block(r,0,1,cols)); + VERIFY_IS_APPROX(m1.row(r).transpose() * m1.col(c).transpose(), m1.block(r,0,1,cols).transpose() * m1.block(0,c,rows,1).transpose()); + VERIFY_IS_APPROX(m1.block(0,c,rows,1) * m1.row(r), m1.block(0,c,rows,1) * m1.block(r,0,1,cols)); + VERIFY_IS_APPROX(m1.col(c) * m1.block(r,0,1,cols), m1.block(0,c,rows,1) * m1.block(r,0,1,cols)); + VERIFY_IS_APPROX(m1.leftCols(1) * m1.row(r), m1.block(0,0,rows,1) * m1.block(r,0,1,cols)); + VERIFY_IS_APPROX(m1.col(c) * m1.topRows(1), m1.block(0,c,rows,1) * m1.block(0,0,1,cols)); } -- cgit v1.2.3 From d6236d3b26f6b652c452d884c440099892fdcdba Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 10 Sep 2014 11:54:20 +0200 Subject: Fix bug #791: infinite loop in JacobiSVD in the presence of NaN. --- Eigen/src/SVD/JacobiSVD.h | 3 ++- test/jacobisvd.cpp | 16 ++++++++++++---- 2 files changed, 14 insertions(+), 5 deletions(-) (limited to 'test') diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index 6f3907f5d..6ff689de3 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -741,7 +741,8 @@ JacobiSVD::compute(const MatrixType& matrix, unsig EIGEN_USING_STD_MATH(max); RealScalar threshold = (max)(considerAsZero, precision * (max)(abs(m_workMatrix.coeff(p,p)), abs(m_workMatrix.coeff(q,q)))); - if((max)(abs(m_workMatrix.coeff(p,q)),abs(m_workMatrix.coeff(q,p))) > threshold) + // We compare both values to threshold instead of calling max to be robust to NaN (See bug 791) + if(abs(m_workMatrix.coeff(p,q))>threshold || abs(m_workMatrix.coeff(q,p)) > threshold) { finished = false; diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp index 36721b496..422d46695 100644 --- a/test/jacobisvd.cpp +++ b/test/jacobisvd.cpp @@ -315,16 +315,23 @@ void jacobisvd_inf_nan() VERIFY(sub(some_inf, some_inf) != sub(some_inf, some_inf)); svd.compute(MatrixType::Constant(10,10,some_inf), ComputeFullU | ComputeFullV); - Scalar some_nan = zero() / zero(); - VERIFY(some_nan != some_nan); - svd.compute(MatrixType::Constant(10,10,some_nan), ComputeFullU | ComputeFullV); + Scalar nan = std::numeric_limits::quiet_NaN(); + VERIFY(nan != nan); + svd.compute(MatrixType::Constant(10,10,nan), ComputeFullU | ComputeFullV); MatrixType m = MatrixType::Zero(10,10); m(internal::random(0,9), internal::random(0,9)) = some_inf; svd.compute(m, ComputeFullU | ComputeFullV); m = MatrixType::Zero(10,10); - m(internal::random(0,9), internal::random(0,9)) = some_nan; + m(internal::random(0,9), internal::random(0,9)) = nan; + svd.compute(m, ComputeFullU | ComputeFullV); + + // regression test for bug 791 + m.resize(3,3); + m << 0, 2*NumTraits::epsilon(), 0.5, + 0, -0.5, 0, + nan, 0, 0; svd.compute(m, ComputeFullU | ComputeFullV); } @@ -437,6 +444,7 @@ void test_jacobisvd() // Test on inf/nan matrix CALL_SUBTEST_7( jacobisvd_inf_nan() ); + CALL_SUBTEST_10( jacobisvd_inf_nan() ); } CALL_SUBTEST_7(( jacobisvd(MatrixXf(internal::random(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2), internal::random(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2))) )); -- cgit v1.2.3 From 84a7ead059917964cc8719f372d7d153bb2cad53 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 10 Sep 2014 11:59:45 +0200 Subject: Add one more regression test for bug #791. --- test/jacobisvd.cpp | 7 +++++++ 1 file changed, 7 insertions(+) (limited to 'test') diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp index 422d46695..7dbb29c59 100644 --- a/test/jacobisvd.cpp +++ b/test/jacobisvd.cpp @@ -333,6 +333,13 @@ void jacobisvd_inf_nan() 0, -0.5, 0, nan, 0, 0; svd.compute(m, ComputeFullU | ComputeFullV); + + m.resize(4,4); + m << 1, 0, 0, 0, + 0, 3, 1, 2e-308, + 1, 0, 1, nan, + 0, nan, nan, 0; + svd.compute(m, ComputeFullU | ComputeFullV); } // Regression test for bug 286: JacobiSVD loops indefinitely with some -- cgit v1.2.3 From 5e890d3ad78a7e5c491a43202993d617fffb964a Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 10 Sep 2014 23:11:58 +0200 Subject: Improve further the accuracy of JacobiSVD wrt under/overflow while improving speed for small matrices (hypot is very slow). --- Eigen/src/SVD/JacobiSVD.h | 17 +++++++++-------- test/jacobisvd.cpp | 29 ++++++++++++++++++++++++++--- 2 files changed, 35 insertions(+), 11 deletions(-) (limited to 'test') diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index 6ff689de3..3ab8a4c8a 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -425,18 +425,19 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, JacobiRotation rot1; RealScalar t = m.coeff(0,0) + m.coeff(1,1); RealScalar d = m.coeff(1,0) - m.coeff(0,1); - if(t == RealScalar(0)) + + if(d == RealScalar(0)) { - rot1.c() = RealScalar(0); - rot1.s() = d > RealScalar(0) ? RealScalar(1) : RealScalar(-1); + rot1.s() = RealScalar(0); + rot1.c() = RealScalar(1); } else { - RealScalar t2d2 = numext::hypot(t,d); - rot1.c() = abs(t)/t2d2; - rot1.s() = d/t2d2; - if(tmakeJacobi(m,0,1); diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp index 7dbb29c59..cd04db5be 100644 --- a/test/jacobisvd.cpp +++ b/test/jacobisvd.cpp @@ -354,11 +354,33 @@ void jacobisvd_underoverflow() Matrix2d M; M << -7.90884e-313, -4.94e-324, 0, 5.60844e-313; + JacobiSVD svd; + svd.compute(M,ComputeFullU|ComputeFullV); + jacobisvd_check_full(M,svd); + + VectorXd value_set(9); + value_set << 0, 1, -1, 5.60844e-313, -5.60844e-313, 4.94e-324, -4.94e-324, -4.94e-223, 4.94e-223; + Array4i id(0,0,0,0); + int k = 0; + do + { + M << value_set(id(0)), value_set(id(1)), value_set(id(2)), value_set(id(3)); + svd.compute(M,ComputeFullU|ComputeFullV); + jacobisvd_check_full(M,svd); + + id(k)++; + if(id(k)>=value_set.size()) + { + while(k<3 && id(k)>=value_set.size()) id(++k)++; + id.head(k).setZero(); + k=0; + } + + } while((id svd; - svd.compute(M); // just check we don't loop indefinitely // Check for overflow: Matrix3d M3; @@ -367,7 +389,8 @@ void jacobisvd_underoverflow() -8.7190887618028355e+307, -7.3453213709232193e+307, -2.4367363684472105e+307; JacobiSVD svd3; - svd3.compute(M3); // just check we don't loop indefinitely + svd3.compute(M3,ComputeFullU|ComputeFullV); // just check we don't loop indefinitely + jacobisvd_check_full(M3,svd3); } void jacobisvd_preallocate() -- cgit v1.2.3 From 9452eb38f812194a676edc1b9eb9d08b7bc0f297 Mon Sep 17 00:00:00 2001 From: Jitse Niesen Date: Fri, 12 Sep 2014 14:52:35 +0100 Subject: Make UpperBidiagonalization accept row-major matrices (bug #769) * Give temporary workspace the same storage order as original matrix * Take storage order into account when determining inner stride of rows and columns * Change one test to use a row-major matrix. --- Eigen/src/SVD/UpperBidiagonalization.h | 29 ++++++++++++++++++++++------- test/upperbidiagonalization.cpp | 2 +- 2 files changed, 23 insertions(+), 8 deletions(-) (limited to 'test') diff --git a/Eigen/src/SVD/UpperBidiagonalization.h b/Eigen/src/SVD/UpperBidiagonalization.h index 225b19e3c..64906bf0c 100644 --- a/Eigen/src/SVD/UpperBidiagonalization.h +++ b/Eigen/src/SVD/UpperBidiagonalization.h @@ -154,14 +154,19 @@ void upperbidiagonalization_blocked_helper(MatrixType& A, typename MatrixType::RealScalar *diagonal, typename MatrixType::RealScalar *upper_diagonal, typename MatrixType::Index bs, - Ref > X, - Ref > Y) + Ref::Flags & RowMajorBit> > X, + Ref::Flags & RowMajorBit> > Y) { typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; - typedef Ref > SubColumnType; - typedef Ref, 0, InnerStride<> > SubRowType; - typedef Ref > SubMatType; + enum { StorageOrder = traits::Flags & RowMajorBit }; + typedef InnerStride ColInnerStride; + typedef InnerStride RowInnerStride; + typedef Ref, 0, ColInnerStride> SubColumnType; + typedef Ref, 0, RowInnerStride> SubRowType; + typedef Ref > SubMatType; Index brows = A.rows(); Index bcols = A.cols(); @@ -288,8 +293,18 @@ void upperbidiagonalization_inplace_blocked(MatrixType& A, BidiagType& bidiagona Index cols = A.cols(); Index size = (std::min)(rows, cols); - Matrix X(rows,maxBlockSize); - Matrix Y(cols,maxBlockSize); + // X and Y are work space + enum { StorageOrder = traits::Flags & RowMajorBit }; + Matrix X(rows,maxBlockSize); + Matrix Y(cols,maxBlockSize); Index blockSize = (std::min)(maxBlockSize,size); Index k = 0; diff --git a/test/upperbidiagonalization.cpp b/test/upperbidiagonalization.cpp index d15bf588b..847b34b55 100644 --- a/test/upperbidiagonalization.cpp +++ b/test/upperbidiagonalization.cpp @@ -35,7 +35,7 @@ void test_upperbidiagonalization() CALL_SUBTEST_1( upperbidiag(MatrixXf(3,3)) ); CALL_SUBTEST_2( upperbidiag(MatrixXd(17,12)) ); CALL_SUBTEST_3( upperbidiag(MatrixXcf(20,20)) ); - CALL_SUBTEST_4( upperbidiag(MatrixXcd(16,15)) ); + CALL_SUBTEST_4( upperbidiag(Matrix,Dynamic,Dynamic,RowMajor>(16,15)) ); CALL_SUBTEST_5( upperbidiag(Matrix()) ); CALL_SUBTEST_6( upperbidiag(Matrix()) ); CALL_SUBTEST_7( upperbidiag(Matrix()) ); -- cgit v1.2.3