diff options
author | Carlos Hernandez <chernand@google.com> | 2014-08-05 17:53:32 -0700 |
---|---|---|
committer | Carlos Hernandez <chernand@google.com> | 2014-08-05 17:54:05 -0700 |
commit | 7faaa9f3f0df9d23790277834d426c3d992ac3ba (patch) | |
tree | b788ae3b96daf9f5a79d8ec434e1e9edd56b3a72 /test | |
parent | 810535bb0c575a003b32076e5352ab8fd3f23a1c (diff) | |
download | eigen-7faaa9f3f0df9d23790277834d426c3d992ac3ba.tar.gz |
Update Eigen to the latest stable release, 3.2.2android-wear-5.1.1_r1android-wear-5.1.0_r1android-wear-5.0.0_r1android-l-preview_r2android-cts-5.1_r9android-cts-5.1_r8android-cts-5.1_r7android-cts-5.1_r6android-cts-5.1_r5android-cts-5.1_r4android-cts-5.1_r3android-cts-5.1_r28android-cts-5.1_r27android-cts-5.1_r26android-cts-5.1_r25android-cts-5.1_r24android-cts-5.1_r23android-cts-5.1_r22android-cts-5.1_r21android-cts-5.1_r20android-cts-5.1_r2android-cts-5.1_r19android-cts-5.1_r18android-cts-5.1_r17android-cts-5.1_r16android-cts-5.1_r15android-cts-5.1_r14android-cts-5.1_r13android-cts-5.1_r10android-cts-5.1_r1android-cts-5.0_r9android-cts-5.0_r8android-cts-5.0_r7android-cts-5.0_r6android-cts-5.0_r5android-cts-5.0_r4android-cts-5.0_r3android-5.1.1_r9android-5.1.1_r8android-5.1.1_r7android-5.1.1_r6android-5.1.1_r5android-5.1.1_r4android-5.1.1_r38android-5.1.1_r37android-5.1.1_r36android-5.1.1_r35android-5.1.1_r34android-5.1.1_r33android-5.1.1_r30android-5.1.1_r3android-5.1.1_r29android-5.1.1_r28android-5.1.1_r26android-5.1.1_r25android-5.1.1_r24android-5.1.1_r23android-5.1.1_r22android-5.1.1_r20android-5.1.1_r2android-5.1.1_r19android-5.1.1_r18android-5.1.1_r17android-5.1.1_r16android-5.1.1_r15android-5.1.1_r14android-5.1.1_r13android-5.1.1_r12android-5.1.1_r10android-5.1.1_r1android-5.1.0_r5android-5.1.0_r4android-5.1.0_r3android-5.1.0_r1android-5.0.2_r3android-5.0.2_r1android-5.0.1_r1android-5.0.0_r7android-5.0.0_r6android-5.0.0_r5.1android-5.0.0_r5android-5.0.0_r4android-5.0.0_r3android-5.0.0_r2android-5.0.0_r1lollipop-wear-releaselollipop-releaselollipop-mr1-wfc-releaselollipop-mr1-releaselollipop-mr1-fi-releaselollipop-mr1-devlollipop-mr1-cts-releaselollipop-devlollipop-cts-releasel-preview
./Eigen/src/Core/util/NonMPL2.h is left untouched, so that
usage of non MPL2 code is disabled.
Change-Id: I86fc9257b3c30d0ca15b268d4ef07bf038bba7ca
Diffstat (limited to 'test')
85 files changed, 2745 insertions, 708 deletions
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 6f8fc4ae3..ccb0fc798 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -86,6 +86,27 @@ else() ei_add_property(EIGEN_MISSING_BACKENDS "PaStiX, ") endif() +if(METIS_FOUND) + add_definitions("-DEIGEN_METIS_SUPPORT") + include_directories(${METIS_INCLUDES}) + ei_add_property(EIGEN_TESTED_BACKENDS "METIS, ") +else() + ei_add_property(EIGEN_MISSING_BACKENDS "METIS, ") +endif() + +find_package(SPQR) +if(SPQR_FOUND AND BLAS_FOUND AND LAPACK_FOUND) + if(CHOLMOD_FOUND) + add_definitions("-DEIGEN_SPQR_SUPPORT") + include_directories(${SPQR_INCLUDES}) + set(SPQR_ALL_LIBS ${SPQR_LIBRARIES} ${CHOLMOD_LIBRARIES} ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES}) + set(SPARSE_LIBS ${SPARSE_LIBS} ${SPQR_ALL_LIBS}) + ei_add_property(EIGEN_TESTED_BACKENDS "SPQR, ") + else(CHOLMOD_FOUND) + ei_add_property(EIGEN_MISSING_BACKENDS "SPQR, ") + endif(CHOLMOD_FOUND) +endif() + option(EIGEN_TEST_NOQT "Disable Qt support in unit tests" OFF) if(NOT EIGEN_TEST_NOQT) find_package(Qt4) @@ -101,6 +122,9 @@ if(TEST_LIB) add_definitions("-DEIGEN_EXTERN_INSTANTIATIONS=1") endif(TEST_LIB) +set_property(GLOBAL PROPERTY EIGEN_CURRENT_SUBPROJECT "Official") +add_custom_target(BuildOfficial) + ei_add_test(meta) ei_add_test(sizeof) ei_add_test(dynalloc) @@ -129,13 +153,14 @@ ei_add_test(diagonal) ei_add_test(miscmatrices) ei_add_test(commainitializer) ei_add_test(smallvectors) -ei_add_test(map) +ei_add_test(mapped_matrix) ei_add_test(mapstride) ei_add_test(mapstaticmethods) ei_add_test(array) ei_add_test(array_for_matrix) ei_add_test(array_replicate) ei_add_test(array_reverse) +ei_add_test(ref) ei_add_test(triangular) ei_add_test(selfadjoint) ei_add_test(product_selfadjoint) @@ -162,6 +187,8 @@ ei_add_test(schur_complex) ei_add_test(eigensolver_selfadjoint) ei_add_test(eigensolver_generic) ei_add_test(eigensolver_complex) +ei_add_test(real_qz) +ei_add_test(eigensolver_generalized_real) ei_add_test(jacobi) ei_add_test(jacobisvd) ei_add_test(geo_orthomethods) @@ -177,9 +204,6 @@ ei_add_test(stdvector_overload) ei_add_test(stdlist) ei_add_test(stddeque) ei_add_test(resize) -if(QT4_FOUND) - ei_add_test(qtvector "" "${QT_QTCORE_LIBRARY}") -endif(QT4_FOUND) ei_add_test(sparse_vector) ei_add_test(sparse_basic) ei_add_test(sparse_product) @@ -190,7 +214,6 @@ ei_add_test(swap) ei_add_test(conservative_resize) ei_add_test(permutationmatrices) ei_add_test(sparse_permutations) -ei_add_test(eigen2support) ei_add_test(nullary) ei_add_test(nesting_ops "${CMAKE_CXX_FLAGS_DEBUG}") ei_add_test(zerosized) @@ -198,11 +221,21 @@ ei_add_test(dontalign) ei_add_test(sizeoverflow) ei_add_test(prec_inverse_4x4) ei_add_test(vectorwiseop) +ei_add_test(special_numbers) ei_add_test(simplicial_cholesky) ei_add_test(conjugate_gradient) ei_add_test(bicgstab) +ei_add_test(sparselu) +ei_add_test(sparseqr) +# ei_add_test(denseLM) + +if(QT4_FOUND) + ei_add_test(qtvector "" "${QT_QTCORE_LIBRARY}") +endif(QT4_FOUND) + +ei_add_test(eigen2support) if(UMFPACK_FOUND) ei_add_test(umfpack_support "" "${UMFPACK_ALL_LIBS}") @@ -224,6 +257,14 @@ if(PASTIX_FOUND AND (SCOTCH_FOUND OR METIS_FOUND)) ei_add_test(pastix_support "" "${PASTIX_ALL_LIBS}") endif() +if(SPQR_FOUND AND CHOLMOD_FOUND) + ei_add_test(spqr_support "" "${SPQR_ALL_LIBS}") +endif() + +if(METIS_FOUND) +ei_add_test(metis_support "" "${METIS_LIBRARIES}") +endif() + string(TOLOWER "${CMAKE_CXX_COMPILER}" cmake_cxx_compiler_tolower) if(cmake_cxx_compiler_tolower MATCHES "qcc") set(CXX_IS_QCC "ON") diff --git a/test/adjoint.cpp b/test/adjoint.cpp index b6cf0a68b..ea36f7841 100644 --- a/test/adjoint.cpp +++ b/test/adjoint.cpp @@ -11,11 +11,54 @@ #include "main.h" +template<bool IsInteger> struct adjoint_specific; + +template<> struct adjoint_specific<true> { + template<typename Vec, typename Mat, typename Scalar> + static void run(const Vec& v1, const Vec& v2, Vec& v3, const Mat& square, Scalar s1, Scalar s2) { + VERIFY(test_isApproxWithRef((s1 * v1 + s2 * v2).dot(v3), numext::conj(s1) * v1.dot(v3) + numext::conj(s2) * v2.dot(v3), 0)); + VERIFY(test_isApproxWithRef(v3.dot(s1 * v1 + s2 * v2), s1*v3.dot(v1)+s2*v3.dot(v2), 0)); + + // check compatibility of dot and adjoint + VERIFY(test_isApproxWithRef(v1.dot(square * v2), (square.adjoint() * v1).dot(v2), 0)); + } +}; + +template<> struct adjoint_specific<false> { + template<typename Vec, typename Mat, typename Scalar> + static void run(const Vec& v1, const Vec& v2, Vec& v3, const Mat& square, Scalar s1, Scalar s2) { + typedef typename NumTraits<Scalar>::Real RealScalar; + using std::abs; + + RealScalar ref = NumTraits<Scalar>::IsInteger ? RealScalar(0) : (std::max)((s1 * v1 + s2 * v2).norm(),v3.norm()); + VERIFY(test_isApproxWithRef((s1 * v1 + s2 * v2).dot(v3), numext::conj(s1) * v1.dot(v3) + numext::conj(s2) * v2.dot(v3), ref)); + VERIFY(test_isApproxWithRef(v3.dot(s1 * v1 + s2 * v2), s1*v3.dot(v1)+s2*v3.dot(v2), ref)); + + VERIFY_IS_APPROX(v1.squaredNorm(), v1.norm() * v1.norm()); + // check normalized() and normalize() + VERIFY_IS_APPROX(v1, v1.norm() * v1.normalized()); + v3 = v1; + v3.normalize(); + VERIFY_IS_APPROX(v1, v1.norm() * v3); + VERIFY_IS_APPROX(v3, v1.normalized()); + VERIFY_IS_APPROX(v3.norm(), RealScalar(1)); + + // check compatibility of dot and adjoint + ref = NumTraits<Scalar>::IsInteger ? 0 : (std::max)((std::max)(v1.norm(),v2.norm()),(std::max)((square * v2).norm(),(square.adjoint() * v1).norm())); + VERIFY(internal::isMuchSmallerThan(abs(v1.dot(square * v2) - (square.adjoint() * v1).dot(v2)), ref, test_precision<Scalar>())); + + // check that Random().normalized() works: tricky as the random xpr must be evaluated by + // normalized() in order to produce a consistent result. + VERIFY_IS_APPROX(Vec::Random(v1.size()).normalized().norm(), RealScalar(1)); + } +}; + template<typename MatrixType> void adjoint(const MatrixType& m) { /* this test covers the following files: Transpose.h Conjugate.h Dot.h */ + using std::abs; typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; @@ -43,45 +86,21 @@ template<typename MatrixType> void adjoint(const MatrixType& m) // check multiplicative behavior VERIFY_IS_APPROX((m1.adjoint() * m2).adjoint(), m2.adjoint() * m1); - VERIFY_IS_APPROX((s1 * m1).adjoint(), internal::conj(s1) * m1.adjoint()); + VERIFY_IS_APPROX((s1 * m1).adjoint(), numext::conj(s1) * m1.adjoint()); - // check basic properties of dot, norm, norm2 - typedef typename NumTraits<Scalar>::Real RealScalar; + // check basic properties of dot, squaredNorm + VERIFY_IS_APPROX(numext::conj(v1.dot(v2)), v2.dot(v1)); + VERIFY_IS_APPROX(numext::real(v1.dot(v1)), v1.squaredNorm()); - RealScalar ref = NumTraits<Scalar>::IsInteger ? RealScalar(0) : (std::max)((s1 * v1 + s2 * v2).norm(),v3.norm()); - VERIFY(test_isApproxWithRef((s1 * v1 + s2 * v2).dot(v3), internal::conj(s1) * v1.dot(v3) + internal::conj(s2) * v2.dot(v3), ref)); - VERIFY(test_isApproxWithRef(v3.dot(s1 * v1 + s2 * v2), s1*v3.dot(v1)+s2*v3.dot(v2), ref)); - VERIFY_IS_APPROX(internal::conj(v1.dot(v2)), v2.dot(v1)); - VERIFY_IS_APPROX(internal::real(v1.dot(v1)), v1.squaredNorm()); - if(!NumTraits<Scalar>::IsInteger) { - VERIFY_IS_APPROX(v1.squaredNorm(), v1.norm() * v1.norm()); - // check normalized() and normalize() - VERIFY_IS_APPROX(v1, v1.norm() * v1.normalized()); - v3 = v1; - v3.normalize(); - VERIFY_IS_APPROX(v1, v1.norm() * v3); - VERIFY_IS_APPROX(v3, v1.normalized()); - VERIFY_IS_APPROX(v3.norm(), RealScalar(1)); - } - VERIFY_IS_MUCH_SMALLER_THAN(internal::abs(vzero.dot(v1)), static_cast<RealScalar>(1)); + adjoint_specific<NumTraits<Scalar>::IsInteger>::run(v1, v2, v3, square, s1, s2); - // check compatibility of dot and adjoint + VERIFY_IS_MUCH_SMALLER_THAN(abs(vzero.dot(v1)), static_cast<RealScalar>(1)); - ref = NumTraits<Scalar>::IsInteger ? 0 : (std::max)((std::max)(v1.norm(),v2.norm()),(std::max)((square * v2).norm(),(square.adjoint() * v1).norm())); - VERIFY(test_isApproxWithRef(v1.dot(square * v2), (square.adjoint() * v1).dot(v2), ref)); - // like in testBasicStuff, test operator() to check const-qualification Index r = internal::random<Index>(0, rows-1), c = internal::random<Index>(0, cols-1); - VERIFY_IS_APPROX(m1.conjugate()(r,c), internal::conj(m1(r,c))); - VERIFY_IS_APPROX(m1.adjoint()(c,r), internal::conj(m1(r,c))); - - if(!NumTraits<Scalar>::IsInteger) - { - // check that Random().normalized() works: tricky as the random xpr must be evaluated by - // normalized() in order to produce a consistent result. - VERIFY_IS_APPROX(VectorType::Random(rows).normalized().norm(), RealScalar(1)); - } + VERIFY_IS_APPROX(m1.conjugate()(r,c), numext::conj(m1(r,c))); + VERIFY_IS_APPROX(m1.adjoint()(c,r), numext::conj(m1(r,c))); // check inplace transpose m3 = m1; diff --git a/test/array.cpp b/test/array.cpp index 3548fa641..c607da631 100644 --- a/test/array.cpp +++ b/test/array.cpp @@ -13,7 +13,6 @@ template<typename ArrayType> void array(const ArrayType& m) { typedef typename ArrayType::Index Index; typedef typename ArrayType::Scalar Scalar; - typedef typename NumTraits<Scalar>::Real RealScalar; typedef Array<Scalar, ArrayType::RowsAtCompileTime, 1> ColVectorType; typedef Array<Scalar, 1, ArrayType::ColsAtCompileTime> RowVectorType; @@ -64,9 +63,12 @@ template<typename ArrayType> void array(const ArrayType& m) VERIFY_IS_APPROX(m1, m3 / m2); // reductions - VERIFY_IS_APPROX(m1.colwise().sum().sum(), m1.sum()); - VERIFY_IS_APPROX(m1.rowwise().sum().sum(), m1.sum()); - if (!internal::isApprox(m1.sum(), (m1+m2).sum(), test_precision<Scalar>())) + VERIFY_IS_APPROX(m1.abs().colwise().sum().sum(), m1.abs().sum()); + VERIFY_IS_APPROX(m1.abs().rowwise().sum().sum(), m1.abs().sum()); + using std::abs; + VERIFY_IS_MUCH_SMALLER_THAN(abs(m1.colwise().sum().sum() - m1.sum()), m1.abs().sum()); + VERIFY_IS_MUCH_SMALLER_THAN(abs(m1.rowwise().sum().sum() - m1.sum()), m1.abs().sum()); + if (!internal::isMuchSmallerThan(abs(m1.sum() - (m1+m2).sum()), m1.abs().sum(), test_precision<Scalar>())) VERIFY_IS_NOT_APPROX(((m1+m2).rowwise().sum()).sum(), m1.sum()); VERIFY_IS_APPROX(m1.colwise().sum(), m1.colwise().redux(internal::scalar_sum_op<Scalar>())); @@ -83,10 +85,10 @@ template<typename ArrayType> void array(const ArrayType& m) template<typename ArrayType> void comparisons(const ArrayType& m) { + using std::abs; typedef typename ArrayType::Index Index; typedef typename ArrayType::Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; - typedef Array<Scalar, ArrayType::RowsAtCompileTime, 1> VectorType; Index rows = m.rows(); Index cols = m.cols(); @@ -120,7 +122,7 @@ template<typename ArrayType> void comparisons(const ArrayType& m) Scalar mid = (m1.cwiseAbs().minCoeff() + m1.cwiseAbs().maxCoeff())/Scalar(2); for (int j=0; j<cols; ++j) for (int i=0; i<rows; ++i) - m3(i,j) = internal::abs(m1(i,j))<mid ? 0 : m1(i,j); + m3(i,j) = abs(m1(i,j))<mid ? 0 : m1(i,j); VERIFY_IS_APPROX( (m1.abs()<ArrayType::Constant(rows,cols,mid)) .select(ArrayType::Zero(rows,cols),m1), m3); // shorter versions: @@ -149,6 +151,8 @@ template<typename ArrayType> void comparisons(const ArrayType& m) template<typename ArrayType> void array_real(const ArrayType& m) { + using std::abs; + using std::sqrt; typedef typename ArrayType::Index Index; typedef typename ArrayType::Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; @@ -157,58 +161,58 @@ template<typename ArrayType> void array_real(const ArrayType& m) Index cols = m.cols(); ArrayType m1 = ArrayType::Random(rows, cols), - m2 = ArrayType::Random(rows, cols), - m3(rows, cols); + m2 = ArrayType::Random(rows, cols), + m3(rows, cols); Scalar s1 = internal::random<Scalar>(); // these tests are mostly to check possible compilation issues. - VERIFY_IS_APPROX(m1.sin(), std::sin(m1)); - VERIFY_IS_APPROX(m1.sin(), internal::sin(m1)); - VERIFY_IS_APPROX(m1.cos(), std::cos(m1)); - VERIFY_IS_APPROX(m1.cos(), internal::cos(m1)); - VERIFY_IS_APPROX(m1.asin(), std::asin(m1)); - VERIFY_IS_APPROX(m1.asin(), internal::asin(m1)); - VERIFY_IS_APPROX(m1.acos(), std::acos(m1)); - VERIFY_IS_APPROX(m1.acos(), internal::acos(m1)); - VERIFY_IS_APPROX(m1.tan(), std::tan(m1)); - VERIFY_IS_APPROX(m1.tan(), internal::tan(m1)); + VERIFY_IS_APPROX(m1.sin(), sin(m1)); + VERIFY_IS_APPROX(m1.cos(), cos(m1)); + VERIFY_IS_APPROX(m1.asin(), asin(m1)); + VERIFY_IS_APPROX(m1.acos(), acos(m1)); + VERIFY_IS_APPROX(m1.tan(), tan(m1)); - VERIFY_IS_APPROX(internal::cos(m1+RealScalar(3)*m2), internal::cos((m1+RealScalar(3)*m2).eval())); - VERIFY_IS_APPROX(std::cos(m1+RealScalar(3)*m2), std::cos((m1+RealScalar(3)*m2).eval())); + VERIFY_IS_APPROX(cos(m1+RealScalar(3)*m2), cos((m1+RealScalar(3)*m2).eval())); - VERIFY_IS_APPROX(m1.abs().sqrt(), std::sqrt(std::abs(m1))); - VERIFY_IS_APPROX(m1.abs().sqrt(), internal::sqrt(internal::abs(m1))); - VERIFY_IS_APPROX(m1.abs(), internal::sqrt(internal::abs2(m1))); + VERIFY_IS_APPROX(m1.abs().sqrt(), sqrt(abs(m1))); + VERIFY_IS_APPROX(m1.abs(), sqrt(numext::abs2(m1))); - VERIFY_IS_APPROX(internal::abs2(internal::real(m1)) + internal::abs2(internal::imag(m1)), internal::abs2(m1)); - VERIFY_IS_APPROX(internal::abs2(std::real(m1)) + internal::abs2(std::imag(m1)), internal::abs2(m1)); + VERIFY_IS_APPROX(numext::abs2(numext::real(m1)) + numext::abs2(numext::imag(m1)), numext::abs2(m1)); + VERIFY_IS_APPROX(numext::abs2(real(m1)) + numext::abs2(imag(m1)), numext::abs2(m1)); if(!NumTraits<Scalar>::IsComplex) - VERIFY_IS_APPROX(internal::real(m1), m1); + VERIFY_IS_APPROX(numext::real(m1), m1); - VERIFY_IS_APPROX(m1.abs().log(), std::log(std::abs(m1))); - VERIFY_IS_APPROX(m1.abs().log(), internal::log(internal::abs(m1))); + // shift argument of logarithm so that it is not zero + Scalar smallNumber = NumTraits<Scalar>::dummy_precision(); + VERIFY_IS_APPROX((m1.abs() + smallNumber).log() , log(abs(m1) + smallNumber)); - VERIFY_IS_APPROX(m1.exp(), std::exp(m1)); - VERIFY_IS_APPROX(m1.exp() * m2.exp(), std::exp(m1+m2)); - VERIFY_IS_APPROX(m1.exp(), internal::exp(m1)); - VERIFY_IS_APPROX(m1.exp() / m2.exp(), std::exp(m1-m2)); + VERIFY_IS_APPROX(m1.exp() * m2.exp(), exp(m1+m2)); + VERIFY_IS_APPROX(m1.exp(), exp(m1)); + VERIFY_IS_APPROX(m1.exp() / m2.exp(),(m1-m2).exp()); VERIFY_IS_APPROX(m1.pow(2), m1.square()); - VERIFY_IS_APPROX(std::pow(m1,2), m1.square()); + VERIFY_IS_APPROX(pow(m1,2), m1.square()); ArrayType exponents = ArrayType::Constant(rows, cols, RealScalar(2)); - VERIFY_IS_APPROX(std::pow(m1,exponents), m1.square()); + VERIFY_IS_APPROX(Eigen::pow(m1,exponents), m1.square()); m3 = m1.abs(); VERIFY_IS_APPROX(m3.pow(RealScalar(0.5)), m3.sqrt()); - VERIFY_IS_APPROX(std::pow(m3,RealScalar(0.5)), m3.sqrt()); + VERIFY_IS_APPROX(pow(m3,RealScalar(0.5)), m3.sqrt()); // scalar by array division - const RealScalar tiny = std::sqrt(std::numeric_limits<RealScalar>::epsilon()); + const RealScalar tiny = sqrt(std::numeric_limits<RealScalar>::epsilon()); s1 += Scalar(tiny); m1 += ArrayType::Constant(rows,cols,Scalar(tiny)); VERIFY_IS_APPROX(s1/m1, s1 * m1.inverse()); + + // check inplace transpose + m3 = m1; + m3.transposeInPlace(); + VERIFY_IS_APPROX(m3,m1.transpose()); + m3.transposeInPlace(); + VERIFY_IS_APPROX(m3,m1); } template<typename ArrayType> void array_complex(const ArrayType& m) @@ -223,11 +227,10 @@ template<typename ArrayType> void array_complex(const ArrayType& m) for (Index i = 0; i < m.rows(); ++i) for (Index j = 0; j < m.cols(); ++j) - m2(i,j) = std::sqrt(m1(i,j)); + m2(i,j) = sqrt(m1(i,j)); VERIFY_IS_APPROX(m1.sqrt(), m2); - VERIFY_IS_APPROX(m1.sqrt(), std::sqrt(m1)); - VERIFY_IS_APPROX(m1.sqrt(), internal::sqrt(m1)); + VERIFY_IS_APPROX(m1.sqrt(), Eigen::sqrt(m1)); } template<typename ArrayType> void min_max(const ArrayType& m) diff --git a/test/array_for_matrix.cpp b/test/array_for_matrix.cpp index 4b637c3a6..9a50f99ab 100644 --- a/test/array_for_matrix.cpp +++ b/test/array_for_matrix.cpp @@ -13,7 +13,6 @@ template<typename MatrixType> void array_for_matrix(const MatrixType& m) { typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits<Scalar>::Real RealScalar; typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVectorType; typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType; @@ -26,10 +25,10 @@ template<typename MatrixType> void array_for_matrix(const MatrixType& m) ColVectorType cv1 = ColVectorType::Random(rows); RowVectorType rv1 = RowVectorType::Random(cols); - + Scalar s1 = internal::random<Scalar>(), s2 = internal::random<Scalar>(); - + // scalar addition VERIFY_IS_APPROX(m1.array() + s1, s1 + m1.array()); VERIFY_IS_APPROX((m1.array() + s1).matrix(), MatrixType::Constant(rows,cols,s1) + m1); @@ -42,10 +41,10 @@ template<typename MatrixType> void array_for_matrix(const MatrixType& m) VERIFY_IS_APPROX(m3, (m1.array() - s1).matrix()); // reductions - VERIFY_IS_MUCH_SMALLER_THAN(m1.colwise().sum().sum() - m1.sum(), m1.cwiseAbs().maxCoeff()); - VERIFY_IS_MUCH_SMALLER_THAN(m1.rowwise().sum().sum() - m1.sum(), m1.cwiseAbs().maxCoeff()); - VERIFY_IS_MUCH_SMALLER_THAN(m1.colwise().sum() + m2.colwise().sum() - (m1+m2).colwise().sum(), (m1+m2).cwiseAbs().maxCoeff()); - VERIFY_IS_MUCH_SMALLER_THAN(m1.rowwise().sum() - m2.rowwise().sum() - (m1-m2).rowwise().sum(), (m1-m2).cwiseAbs().maxCoeff()); + VERIFY_IS_MUCH_SMALLER_THAN(m1.colwise().sum().sum() - m1.sum(), m1.squaredNorm()); + VERIFY_IS_MUCH_SMALLER_THAN(m1.rowwise().sum().sum() - m1.sum(), m1.squaredNorm()); + VERIFY_IS_MUCH_SMALLER_THAN(m1.colwise().sum() + m2.colwise().sum() - (m1+m2).colwise().sum(), (m1+m2).squaredNorm()); + VERIFY_IS_MUCH_SMALLER_THAN(m1.rowwise().sum() - m2.rowwise().sum() - (m1-m2).rowwise().sum(), (m1-m2).squaredNorm()); VERIFY_IS_APPROX(m1.colwise().sum(), m1.colwise().redux(internal::scalar_sum_op<Scalar>())); // vector-wise ops @@ -73,10 +72,10 @@ template<typename MatrixType> void array_for_matrix(const MatrixType& m) template<typename MatrixType> void comparisons(const MatrixType& m) { + using std::abs; typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; - typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; Index rows = m.rows(); Index cols = m.cols(); @@ -110,7 +109,7 @@ template<typename MatrixType> void comparisons(const MatrixType& m) Scalar mid = (m1.cwiseAbs().minCoeff() + m1.cwiseAbs().maxCoeff())/Scalar(2); for (int j=0; j<cols; ++j) for (int i=0; i<rows; ++i) - m3(i,j) = internal::abs(m1(i,j))<mid ? 0 : m1(i,j); + m3(i,j) = abs(m1(i,j))<mid ? 0 : m1(i,j); VERIFY_IS_APPROX( (m1.array().abs()<MatrixType::Constant(rows,cols,mid).array()) .select(MatrixType::Zero(rows,cols),m1), m3); // shorter versions: @@ -133,12 +132,13 @@ template<typename MatrixType> void comparisons(const MatrixType& m) template<typename VectorType> void lpNorm(const VectorType& v) { + using std::sqrt; VectorType u = VectorType::Random(v.size()); VERIFY_IS_APPROX(u.template lpNorm<Infinity>(), u.cwiseAbs().maxCoeff()); VERIFY_IS_APPROX(u.template lpNorm<1>(), u.cwiseAbs().sum()); - VERIFY_IS_APPROX(u.template lpNorm<2>(), internal::sqrt(u.array().abs().square().sum())); - VERIFY_IS_APPROX(internal::pow(u.template lpNorm<5>(), typename VectorType::RealScalar(5)), u.array().abs().pow(5).sum()); + VERIFY_IS_APPROX(u.template lpNorm<2>(), sqrt(u.array().abs().square().sum())); + VERIFY_IS_APPROX(numext::pow(u.template lpNorm<5>(), typename VectorType::RealScalar(5)), u.array().abs().pow(5).sum()); } template<typename MatrixType> void cwise_min_max(const MatrixType& m) @@ -163,11 +163,53 @@ template<typename MatrixType> void cwise_min_max(const MatrixType& m) // min/max with scalar input VERIFY_IS_APPROX(MatrixType::Constant(rows,cols, minM1), m1.cwiseMin( minM1)); - VERIFY_IS_APPROX(m1, m1.cwiseMin( maxM1)); + VERIFY_IS_APPROX(m1, m1.cwiseMin(maxM1)); + VERIFY_IS_APPROX(-m1, (-m1).cwiseMin(-minM1)); + VERIFY_IS_APPROX(-m1.array(), ((-m1).array().min)( -minM1)); VERIFY_IS_APPROX(MatrixType::Constant(rows,cols, maxM1), m1.cwiseMax( maxM1)); - VERIFY_IS_APPROX(m1, m1.cwiseMax( minM1)); + VERIFY_IS_APPROX(m1, m1.cwiseMax(minM1)); + VERIFY_IS_APPROX(-m1, (-m1).cwiseMax(-maxM1)); + VERIFY_IS_APPROX(-m1.array(), ((-m1).array().max)(-maxM1)); + + VERIFY_IS_APPROX(MatrixType::Constant(rows,cols, minM1).array(), (m1.array().min)( minM1)); + VERIFY_IS_APPROX(m1.array(), (m1.array().min)( maxM1)); + + VERIFY_IS_APPROX(MatrixType::Constant(rows,cols, maxM1).array(), (m1.array().max)( maxM1)); + VERIFY_IS_APPROX(m1.array(), (m1.array().max)( minM1)); + +} +template<typename MatrixTraits> void resize(const MatrixTraits& t) +{ + typedef typename MatrixTraits::Index Index; + typedef typename MatrixTraits::Scalar Scalar; + typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType; + typedef Array<Scalar,Dynamic,Dynamic> Array2DType; + typedef Matrix<Scalar,Dynamic,1> VectorType; + typedef Array<Scalar,Dynamic,1> Array1DType; + + Index rows = t.rows(), cols = t.cols(); + + MatrixType m(rows,cols); + VectorType v(rows); + Array2DType a2(rows,cols); + Array1DType a1(rows); + + m.array().resize(rows+1,cols+1); + VERIFY(m.rows()==rows+1 && m.cols()==cols+1); + a2.matrix().resize(rows+1,cols+1); + VERIFY(a2.rows()==rows+1 && a2.cols()==cols+1); + v.array().resize(cols); + VERIFY(v.size()==cols); + a1.matrix().resize(cols); + VERIFY(a1.size()==cols); +} + +void regression_bug_654() +{ + ArrayXf a = RowVectorXf(3); + VectorXf v = Array<float,1,Dynamic>(3); } void test_array_for_matrix() @@ -202,4 +244,10 @@ void test_array_for_matrix() CALL_SUBTEST_5( lpNorm(VectorXf(internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) ); CALL_SUBTEST_4( lpNorm(VectorXcf(internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) ); } + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_4( resize(MatrixXcf(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) ); + CALL_SUBTEST_5( resize(MatrixXf(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) ); + CALL_SUBTEST_6( resize(MatrixXi(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) ); + } + CALL_SUBTEST_6( regression_bug_654() ); } diff --git a/test/array_replicate.cpp b/test/array_replicate.cpp index 94da7425b..f412d1aed 100644 --- a/test/array_replicate.cpp +++ b/test/array_replicate.cpp @@ -16,7 +16,6 @@ template<typename MatrixType> void replicate(const MatrixType& m) */ typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits<Scalar>::Real RealScalar; typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; typedef Matrix<Scalar, Dynamic, Dynamic> MatrixX; typedef Matrix<Scalar, Dynamic, 1> VectorX; diff --git a/test/basicstuff.cpp b/test/basicstuff.cpp index 48db531c1..8c0621ecd 100644 --- a/test/basicstuff.cpp +++ b/test/basicstuff.cpp @@ -52,8 +52,7 @@ template<typename MatrixType> void basicStuff(const MatrixType& m) VERIFY_IS_APPROX( v1, v1); VERIFY_IS_NOT_APPROX( v1, 2*v1); VERIFY_IS_MUCH_SMALLER_THAN( vzero, v1); - if(!NumTraits<Scalar>::IsInteger) - VERIFY_IS_MUCH_SMALLER_THAN( vzero, v1.norm()); + VERIFY_IS_MUCH_SMALLER_THAN( vzero, v1.squaredNorm()); VERIFY_IS_NOT_MUCH_SMALLER_THAN(v1, v1); VERIFY_IS_APPROX( vzero, v1-v1); VERIFY_IS_APPROX( m1, m1); @@ -142,10 +141,10 @@ template<typename MatrixType> void basicStuffComplex(const MatrixType& m) Scalar s1 = internal::random<Scalar>(), s2 = internal::random<Scalar>(); - VERIFY(internal::real(s1)==internal::real_ref(s1)); - VERIFY(internal::imag(s1)==internal::imag_ref(s1)); - internal::real_ref(s1) = internal::real(s2); - internal::imag_ref(s1) = internal::imag(s2); + VERIFY(numext::real(s1)==numext::real_ref(s1)); + VERIFY(numext::imag(s1)==numext::imag_ref(s1)); + numext::real_ref(s1) = numext::real(s2); + numext::imag_ref(s1) = numext::imag(s2); VERIFY(internal::isApprox(s1, s2, NumTraits<RealScalar>::epsilon())); // extended precision in Intel FPUs means that s1 == s2 in the line above is not guaranteed. diff --git a/test/block.cpp b/test/block.cpp index 0969262ca..9ed5d7bc5 100644 --- a/test/block.cpp +++ b/test/block.cpp @@ -10,6 +10,26 @@ #define EIGEN_NO_STATIC_ASSERT // otherwise we fail at compile time on unused paths #include "main.h" +template<typename MatrixType, typename Index, typename Scalar> +typename Eigen::internal::enable_if<!NumTraits<typename MatrixType::Scalar>::IsComplex,typename MatrixType::Scalar>::type +block_real_only(const MatrixType &m1, Index r1, Index r2, Index c1, Index c2, const Scalar& s1) { + // check cwise-Functions: + VERIFY_IS_APPROX(m1.row(r1).cwiseMax(s1), m1.cwiseMax(s1).row(r1)); + VERIFY_IS_APPROX(m1.col(c1).cwiseMin(s1), m1.cwiseMin(s1).col(c1)); + + VERIFY_IS_APPROX(m1.block(r1,c1,r2-r1+1,c2-c1+1).cwiseMin(s1), m1.cwiseMin(s1).block(r1,c1,r2-r1+1,c2-c1+1)); + VERIFY_IS_APPROX(m1.block(r1,c1,r2-r1+1,c2-c1+1).cwiseMax(s1), m1.cwiseMax(s1).block(r1,c1,r2-r1+1,c2-c1+1)); + + return Scalar(0); +} + +template<typename MatrixType, typename Index, typename Scalar> +typename Eigen::internal::enable_if<NumTraits<typename MatrixType::Scalar>::IsComplex,typename MatrixType::Scalar>::type +block_real_only(const MatrixType &, Index, Index, Index, Index, const Scalar&) { + return Scalar(0); +} + + template<typename MatrixType> void block(const MatrixType& m) { typedef typename MatrixType::Index Index; @@ -37,6 +57,8 @@ template<typename MatrixType> void block(const MatrixType& m) Index c1 = internal::random<Index>(0,cols-1); Index c2 = internal::random<Index>(c1,cols-1); + block_real_only(m1, r1, r2, c1, c1, s1); + //check row() and col() VERIFY_IS_EQUAL(m1.col(c1).transpose(), m1.transpose().row(c1)); //check operator(), both constant and non-constant, on row() and col() @@ -51,7 +73,8 @@ template<typename MatrixType> void block(const MatrixType& m) VERIFY_IS_APPROX(m1.col(c1), m1_copy.col(c1) + s1 * m1_copy.col(c2)); m1.col(c1).col(0) += s1 * m1_copy.col(c2); VERIFY_IS_APPROX(m1.col(c1), m1_copy.col(c1) + Scalar(2) * s1 * m1_copy.col(c2)); - + + //check block() Matrix<Scalar,Dynamic,Dynamic> b1(1,1); b1(0,0) = m1(r1,c1); @@ -77,6 +100,12 @@ template<typename MatrixType> void block(const MatrixType& m) // check that fixed block() and block() agree Matrix<Scalar,Dynamic,Dynamic> b = m1.template block<BlockRows,BlockCols>(3,3); VERIFY_IS_EQUAL(b, m1.block(3,3,BlockRows,BlockCols)); + + // same tests with mixed fixed/dynamic size + m1.template block<BlockRows,Dynamic>(1,1,BlockRows,BlockCols) *= s1; + m1.template block<BlockRows,Dynamic>(1,1,BlockRows,BlockCols)(0,3) = m1.template block<2,5>(1,1)(1,2); + Matrix<Scalar,Dynamic,Dynamic> b2 = m1.template block<Dynamic,BlockCols>(3,3,2,5); + VERIFY_IS_EQUAL(b2, m1.block(3,3,BlockRows,BlockCols)); } if (rows>2) @@ -96,11 +125,11 @@ template<typename MatrixType> void block(const MatrixType& m) } // stress some basic stuffs with block matrices - VERIFY(internal::real(ones.col(c1).sum()) == RealScalar(rows)); - VERIFY(internal::real(ones.row(r1).sum()) == RealScalar(cols)); + VERIFY(numext::real(ones.col(c1).sum()) == RealScalar(rows)); + VERIFY(numext::real(ones.row(r1).sum()) == RealScalar(cols)); - VERIFY(internal::real(ones.col(c1).dot(ones.col(c2))) == RealScalar(rows)); - VERIFY(internal::real(ones.row(r1).dot(ones.row(r2))) == RealScalar(cols)); + VERIFY(numext::real(ones.col(c1).dot(ones.col(c2))) == RealScalar(rows)); + VERIFY(numext::real(ones.row(r1).dot(ones.row(r2))) == RealScalar(cols)); // now test some block-inside-of-block. diff --git a/test/cholesky.cpp b/test/cholesky.cpp index 14e01c006..64bcbccc4 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -83,14 +83,14 @@ template<typename MatrixType> void cholesky(const MatrixType& m) symm += a1 * a1.adjoint(); } - SquareMatrixType symmUp = symm.template triangularView<Upper>(); - SquareMatrixType symmLo = symm.template triangularView<Lower>(); - // to test if really Cholesky only uses the upper triangular part, uncomment the following // FIXME: currently that fails !! //symm.template part<StrictlyLower>().setZero(); { + SquareMatrixType symmUp = symm.template triangularView<Upper>(); + SquareMatrixType symmLo = symm.template triangularView<Lower>(); + LLT<SquareMatrixType,Lower> chollo(symmLo); VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix()); vecX = chollo.solve(vecB); @@ -114,6 +114,21 @@ template<typename MatrixType> void cholesky(const MatrixType& m) VERIFY_IS_APPROX(MatrixType(chollo.matrixU().transpose().conjugate()), MatrixType(chollo.matrixL())); VERIFY_IS_APPROX(MatrixType(cholup.matrixL().transpose().conjugate()), MatrixType(cholup.matrixU())); VERIFY_IS_APPROX(MatrixType(cholup.matrixU().transpose().conjugate()), MatrixType(cholup.matrixL())); + + // test some special use cases of SelfCwiseBinaryOp: + MatrixType m1 = MatrixType::Random(rows,cols), m2(rows,cols); + m2 = m1; + m2 += symmLo.template selfadjointView<Lower>().llt().solve(matB); + VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView<Lower>().llt().solve(matB)); + m2 = m1; + m2 -= symmLo.template selfadjointView<Lower>().llt().solve(matB); + VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB)); + m2 = m1; + m2.noalias() += symmLo.template selfadjointView<Lower>().llt().solve(matB); + VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView<Lower>().llt().solve(matB)); + m2 = m1; + m2.noalias() -= symmLo.template selfadjointView<Lower>().llt().solve(matB); + VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB)); } // LDLT @@ -165,22 +180,58 @@ template<typename MatrixType> void cholesky(const MatrixType& m) // restore if(sign == -1) symm = -symm; - } - // test some special use cases of SelfCwiseBinaryOp: - MatrixType m1 = MatrixType::Random(rows,cols), m2(rows,cols); - m2 = m1; - m2 += symmLo.template selfadjointView<Lower>().llt().solve(matB); - VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView<Lower>().llt().solve(matB)); - m2 = m1; - m2 -= symmLo.template selfadjointView<Lower>().llt().solve(matB); - VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB)); - m2 = m1; - m2.noalias() += symmLo.template selfadjointView<Lower>().llt().solve(matB); - VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView<Lower>().llt().solve(matB)); - m2 = m1; - m2.noalias() -= symmLo.template selfadjointView<Lower>().llt().solve(matB); - VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB)); + // check matrices coming from linear constraints with Lagrange multipliers + if(rows>=3) + { + SquareMatrixType A = symm; + int c = internal::random<int>(0,rows-2); + A.bottomRightCorner(c,c).setZero(); + // Make sure a solution exists: + vecX.setRandom(); + vecB = A * vecX; + vecX.setZero(); + ldltlo.compute(A); + VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix()); + vecX = ldltlo.solve(vecB); + VERIFY_IS_APPROX(A * vecX, vecB); + } + + // check non-full rank matrices + if(rows>=3) + { + int r = internal::random<int>(1,rows-1); + Matrix<Scalar,Dynamic,Dynamic> a = Matrix<Scalar,Dynamic,Dynamic>::Random(rows,r); + SquareMatrixType A = a * a.adjoint(); + // Make sure a solution exists: + vecX.setRandom(); + vecB = A * vecX; + vecX.setZero(); + ldltlo.compute(A); + VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix()); + vecX = ldltlo.solve(vecB); + VERIFY_IS_APPROX(A * vecX, vecB); + } + + // check matrices with a wide spectrum + if(rows>=3) + { + RealScalar s = (std::min)(16,std::numeric_limits<RealScalar>::max_exponent10/8); + Matrix<Scalar,Dynamic,Dynamic> a = Matrix<Scalar,Dynamic,Dynamic>::Random(rows,rows); + Matrix<RealScalar,Dynamic,1> d = Matrix<RealScalar,Dynamic,1>::Random(rows); + for(int k=0; k<rows; ++k) + d(k) = d(k)*std::pow(RealScalar(10),internal::random<RealScalar>(-s,s)); + SquareMatrixType A = a * d.asDiagonal() * a.adjoint(); + // Make sure a solution exists: + vecX.setRandom(); + vecB = A * vecX; + vecX.setZero(); + ldltlo.compute(A); + VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix()); + vecX = ldltlo.solve(vecB); + VERIFY_IS_APPROX(A * vecX, vecB); + } + } // update/downdate CALL_SUBTEST(( test_chol_update<SquareMatrixType,LLT>(symm) )); @@ -262,6 +313,45 @@ template<typename MatrixType> void cholesky_bug241(const MatrixType& m) VERIFY_IS_APPROX(matA * vecX, vecB); } +// LDLT is not guaranteed to work for indefinite matrices, but happens to work fine if matrix is diagonal. +// This test checks that LDLT reports correctly that matrix is indefinite. +// See http://forum.kde.org/viewtopic.php?f=74&t=106942 and bug 736 +template<typename MatrixType> void cholesky_definiteness(const MatrixType& m) +{ + eigen_assert(m.rows() == 2 && m.cols() == 2); + MatrixType mat; + { + mat << 1, 0, 0, -1; + LDLT<MatrixType> ldlt(mat); + VERIFY(!ldlt.isNegative()); + VERIFY(!ldlt.isPositive()); + } + { + mat << 1, 2, 2, 1; + LDLT<MatrixType> ldlt(mat); + VERIFY(!ldlt.isNegative()); + VERIFY(!ldlt.isPositive()); + } + { + mat << 0, 0, 0, 0; + LDLT<MatrixType> ldlt(mat); + VERIFY(ldlt.isNegative()); + VERIFY(ldlt.isPositive()); + } + { + mat << 0, 0, 0, 1; + LDLT<MatrixType> ldlt(mat); + VERIFY(!ldlt.isNegative()); + VERIFY(ldlt.isPositive()); + } + { + mat << -1, 0, 0, 0; + LDLT<MatrixType> ldlt(mat); + VERIFY(ldlt.isNegative()); + VERIFY(!ldlt.isPositive()); + } +} + template<typename MatrixType> void cholesky_verify_assert() { MatrixType tmp; @@ -284,11 +374,12 @@ template<typename MatrixType> void cholesky_verify_assert() void test_cholesky() { - int s; + int s = 0; for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1( cholesky(Matrix<double,1,1>()) ); CALL_SUBTEST_3( cholesky(Matrix2d()) ); CALL_SUBTEST_3( cholesky_bug241(Matrix2d()) ); + CALL_SUBTEST_3( cholesky_definiteness(Matrix2d()) ); CALL_SUBTEST_4( cholesky(Matrix3f()) ); CALL_SUBTEST_5( cholesky(Matrix4d()) ); s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE); @@ -306,5 +397,6 @@ void test_cholesky() CALL_SUBTEST_9( LLT<MatrixXf>(10) ); CALL_SUBTEST_9( LDLT<MatrixXf>(10) ); - EIGEN_UNUSED_VARIABLE(s) + TEST_SET_BUT_UNUSED_VARIABLE(s) + TEST_SET_BUT_UNUSED_VARIABLE(nb_temporaries) } diff --git a/test/conservative_resize.cpp b/test/conservative_resize.cpp index 4d11e4075..498421b4c 100644 --- a/test/conservative_resize.cpp +++ b/test/conservative_resize.cpp @@ -60,34 +60,51 @@ void run_matrix_tests() template <typename Scalar> void run_vector_tests() { - typedef Matrix<Scalar, 1, Eigen::Dynamic> MatrixType; + typedef Matrix<Scalar, 1, Eigen::Dynamic> VectorType; - MatrixType m, n; + VectorType m, n; // boundary cases ... - m = n = MatrixType::Random(50); + m = n = VectorType::Random(50); m.conservativeResize(1); VERIFY_IS_APPROX(m, n.segment(0,1)); - m = n = MatrixType::Random(50); + m = n = VectorType::Random(50); m.conservativeResize(50); VERIFY_IS_APPROX(m, n.segment(0,50)); + + m = n = VectorType::Random(50); + m.conservativeResize(m.rows(),1); + VERIFY_IS_APPROX(m, n.segment(0,1)); + + m = n = VectorType::Random(50); + m.conservativeResize(m.rows(),50); + VERIFY_IS_APPROX(m, n.segment(0,50)); // random shrinking ... for (int i=0; i<50; ++i) { const int size = internal::random<int>(1,50); - m = n = MatrixType::Random(50); + m = n = VectorType::Random(50); m.conservativeResize(size); VERIFY_IS_APPROX(m, n.segment(0,size)); + + m = n = VectorType::Random(50); + m.conservativeResize(m.rows(), size); + VERIFY_IS_APPROX(m, n.segment(0,size)); } // random growing with zeroing ... for (int i=0; i<50; ++i) { const int size = internal::random<int>(50,100); - m = n = MatrixType::Random(50); - m.conservativeResizeLike(MatrixType::Zero(size)); + m = n = VectorType::Random(50); + m.conservativeResizeLike(VectorType::Zero(size)); + VERIFY_IS_APPROX(m.segment(0,50), n); + VERIFY( size<=50 || m.segment(50,size-50).sum() == Scalar(0) ); + + m = n = VectorType::Random(50); + m.conservativeResizeLike(Matrix<Scalar,Dynamic,Dynamic>::Zero(1,size)); VERIFY_IS_APPROX(m.segment(0,50), n); VERIFY( size<=50 || m.segment(50,size-50).sum() == Scalar(0) ); } @@ -95,20 +112,23 @@ void run_vector_tests() void test_conservative_resize() { - CALL_SUBTEST_1((run_matrix_tests<int, Eigen::RowMajor>())); - CALL_SUBTEST_1((run_matrix_tests<int, Eigen::ColMajor>())); - CALL_SUBTEST_2((run_matrix_tests<float, Eigen::RowMajor>())); - CALL_SUBTEST_2((run_matrix_tests<float, Eigen::ColMajor>())); - CALL_SUBTEST_3((run_matrix_tests<double, Eigen::RowMajor>())); - CALL_SUBTEST_3((run_matrix_tests<double, Eigen::ColMajor>())); - CALL_SUBTEST_4((run_matrix_tests<std::complex<float>, Eigen::RowMajor>())); - CALL_SUBTEST_4((run_matrix_tests<std::complex<float>, Eigen::ColMajor>())); - CALL_SUBTEST_5((run_matrix_tests<std::complex<double>, Eigen::RowMajor>())); - CALL_SUBTEST_6((run_matrix_tests<std::complex<double>, Eigen::ColMajor>())); - - CALL_SUBTEST_1((run_vector_tests<int>())); - CALL_SUBTEST_2((run_vector_tests<float>())); - CALL_SUBTEST_3((run_vector_tests<double>())); - CALL_SUBTEST_4((run_vector_tests<std::complex<float> >())); - CALL_SUBTEST_5((run_vector_tests<std::complex<double> >())); + for(int i=0; i<g_repeat; ++i) + { + CALL_SUBTEST_1((run_matrix_tests<int, Eigen::RowMajor>())); + CALL_SUBTEST_1((run_matrix_tests<int, Eigen::ColMajor>())); + CALL_SUBTEST_2((run_matrix_tests<float, Eigen::RowMajor>())); + CALL_SUBTEST_2((run_matrix_tests<float, Eigen::ColMajor>())); + CALL_SUBTEST_3((run_matrix_tests<double, Eigen::RowMajor>())); + CALL_SUBTEST_3((run_matrix_tests<double, Eigen::ColMajor>())); + CALL_SUBTEST_4((run_matrix_tests<std::complex<float>, Eigen::RowMajor>())); + CALL_SUBTEST_4((run_matrix_tests<std::complex<float>, Eigen::ColMajor>())); + CALL_SUBTEST_5((run_matrix_tests<std::complex<double>, Eigen::RowMajor>())); + CALL_SUBTEST_6((run_matrix_tests<std::complex<double>, Eigen::ColMajor>())); + + CALL_SUBTEST_1((run_vector_tests<int>())); + CALL_SUBTEST_2((run_vector_tests<float>())); + CALL_SUBTEST_3((run_vector_tests<double>())); + CALL_SUBTEST_4((run_vector_tests<std::complex<float> >())); + CALL_SUBTEST_5((run_vector_tests<std::complex<double> >())); + } } diff --git a/test/corners.cpp b/test/corners.cpp index 4705c5f05..3c64c32a1 100644 --- a/test/corners.cpp +++ b/test/corners.cpp @@ -62,6 +62,16 @@ template<typename MatrixType, int CRows, int CCols, int SRows, int SCols> void c VERIFY_IS_EQUAL((matrix.template bottomLeftCorner<r,c>()), (matrix.template block<r,c>(rows-r,0))); VERIFY_IS_EQUAL((matrix.template bottomRightCorner<r,c>()), (matrix.template block<r,c>(rows-r,cols-c))); + VERIFY_IS_EQUAL((matrix.template topLeftCorner<r,c>()), (matrix.template topLeftCorner<r,Dynamic>(r,c))); + VERIFY_IS_EQUAL((matrix.template topRightCorner<r,c>()), (matrix.template topRightCorner<r,Dynamic>(r,c))); + VERIFY_IS_EQUAL((matrix.template bottomLeftCorner<r,c>()), (matrix.template bottomLeftCorner<r,Dynamic>(r,c))); + VERIFY_IS_EQUAL((matrix.template bottomRightCorner<r,c>()), (matrix.template bottomRightCorner<r,Dynamic>(r,c))); + + VERIFY_IS_EQUAL((matrix.template topLeftCorner<r,c>()), (matrix.template topLeftCorner<Dynamic,c>(r,c))); + VERIFY_IS_EQUAL((matrix.template topRightCorner<r,c>()), (matrix.template topRightCorner<Dynamic,c>(r,c))); + VERIFY_IS_EQUAL((matrix.template bottomLeftCorner<r,c>()), (matrix.template bottomLeftCorner<Dynamic,c>(r,c))); + VERIFY_IS_EQUAL((matrix.template bottomRightCorner<r,c>()), (matrix.template bottomRightCorner<Dynamic,c>(r,c))); + VERIFY_IS_EQUAL((matrix.template topRows<r>()), (matrix.template block<r,cols>(0,0))); VERIFY_IS_EQUAL((matrix.template middleRows<r>(sr)), (matrix.template block<r,cols>(sr,0))); VERIFY_IS_EQUAL((matrix.template bottomRows<r>()), (matrix.template block<r,cols>(rows-r,0))); @@ -74,6 +84,16 @@ template<typename MatrixType, int CRows, int CCols, int SRows, int SCols> void c VERIFY_IS_EQUAL((const_matrix.template bottomLeftCorner<r,c>()), (const_matrix.template block<r,c>(rows-r,0))); VERIFY_IS_EQUAL((const_matrix.template bottomRightCorner<r,c>()), (const_matrix.template block<r,c>(rows-r,cols-c))); + VERIFY_IS_EQUAL((const_matrix.template topLeftCorner<r,c>()), (const_matrix.template topLeftCorner<r,Dynamic>(r,c))); + VERIFY_IS_EQUAL((const_matrix.template topRightCorner<r,c>()), (const_matrix.template topRightCorner<r,Dynamic>(r,c))); + VERIFY_IS_EQUAL((const_matrix.template bottomLeftCorner<r,c>()), (const_matrix.template bottomLeftCorner<r,Dynamic>(r,c))); + VERIFY_IS_EQUAL((const_matrix.template bottomRightCorner<r,c>()), (const_matrix.template bottomRightCorner<r,Dynamic>(r,c))); + + VERIFY_IS_EQUAL((const_matrix.template topLeftCorner<r,c>()), (const_matrix.template topLeftCorner<Dynamic,c>(r,c))); + VERIFY_IS_EQUAL((const_matrix.template topRightCorner<r,c>()), (const_matrix.template topRightCorner<Dynamic,c>(r,c))); + VERIFY_IS_EQUAL((const_matrix.template bottomLeftCorner<r,c>()), (const_matrix.template bottomLeftCorner<Dynamic,c>(r,c))); + VERIFY_IS_EQUAL((const_matrix.template bottomRightCorner<r,c>()), (const_matrix.template bottomRightCorner<Dynamic,c>(r,c))); + VERIFY_IS_EQUAL((const_matrix.template topRows<r>()), (const_matrix.template block<r,cols>(0,0))); VERIFY_IS_EQUAL((const_matrix.template middleRows<r>(sr)), (const_matrix.template block<r,cols>(sr,0))); VERIFY_IS_EQUAL((const_matrix.template bottomRows<r>()), (const_matrix.template block<r,cols>(rows-r,0))); diff --git a/test/cwiseop.cpp b/test/cwiseop.cpp index ca6e4211c..247fa2a09 100644 --- a/test/cwiseop.cpp +++ b/test/cwiseop.cpp @@ -28,17 +28,52 @@ template<typename Scalar> struct AddIfNull { enum { Cost = NumTraits<Scalar>::AddCost }; }; +template<typename MatrixType> +typename Eigen::internal::enable_if<!NumTraits<typename MatrixType::Scalar>::IsInteger,typename MatrixType::Scalar>::type +cwiseops_real_only(MatrixType& m1, MatrixType& m2, MatrixType& m3, MatrixType& mones) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + + VERIFY_IS_APPROX(m1.cwise() / m2, m1.cwise() * (m2.cwise().inverse())); + m3 = m1.cwise().abs().cwise().sqrt(); + VERIFY_IS_APPROX(m3.cwise().square(), m1.cwise().abs()); + VERIFY_IS_APPROX(m1.cwise().square().cwise().sqrt(), m1.cwise().abs()); + VERIFY_IS_APPROX(m1.cwise().abs().cwise().log().cwise().exp() , m1.cwise().abs()); + + VERIFY_IS_APPROX(m1.cwise().pow(2), m1.cwise().square()); + m3 = (m1.cwise().abs().cwise()<=RealScalar(0.01)).select(mones,m1); + VERIFY_IS_APPROX(m3.cwise().pow(-1), m3.cwise().inverse()); + m3 = m1.cwise().abs(); + VERIFY_IS_APPROX(m3.cwise().pow(RealScalar(0.5)), m3.cwise().sqrt()); + +// VERIFY_IS_APPROX(m1.cwise().tan(), m1.cwise().sin().cwise() / m1.cwise().cos()); + VERIFY_IS_APPROX(mones, m1.cwise().sin().cwise().square() + m1.cwise().cos().cwise().square()); + m3 = m1; + m3.cwise() /= m2; + VERIFY_IS_APPROX(m3, m1.cwise() / m2); + + return Scalar(0); +} + +template<typename MatrixType> +typename Eigen::internal::enable_if<NumTraits<typename MatrixType::Scalar>::IsInteger,typename MatrixType::Scalar>::type +cwiseops_real_only(MatrixType& , MatrixType& , MatrixType& , MatrixType& ) +{ + return 0; +} + template<typename MatrixType> void cwiseops(const MatrixType& m) { typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits<Scalar>::Real RealScalar; typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; Index rows = m.rows(); Index cols = m.cols(); MatrixType m1 = MatrixType::Random(rows, cols), + m1bis = m1, m2 = MatrixType::Random(rows, cols), m3(rows, cols), m4(rows, cols), @@ -101,26 +136,6 @@ template<typename MatrixType> void cwiseops(const MatrixType& m) VERIFY_IS_APPROX(m3, m1.cwise() * m2); VERIFY_IS_APPROX(mones, m2.cwise()/m2); - if(!NumTraits<Scalar>::IsInteger) - { - VERIFY_IS_APPROX(m1.cwise() / m2, m1.cwise() * (m2.cwise().inverse())); - m3 = m1.cwise().abs().cwise().sqrt(); - VERIFY_IS_APPROX(m3.cwise().square(), m1.cwise().abs()); - VERIFY_IS_APPROX(m1.cwise().square().cwise().sqrt(), m1.cwise().abs()); - VERIFY_IS_APPROX(m1.cwise().abs().cwise().log().cwise().exp() , m1.cwise().abs()); - - VERIFY_IS_APPROX(m1.cwise().pow(2), m1.cwise().square()); - m3 = (m1.cwise().abs().cwise()<=RealScalar(0.01)).select(mones,m1); - VERIFY_IS_APPROX(m3.cwise().pow(-1), m3.cwise().inverse()); - m3 = m1.cwise().abs(); - VERIFY_IS_APPROX(m3.cwise().pow(RealScalar(0.5)), m3.cwise().sqrt()); - -// VERIFY_IS_APPROX(m1.cwise().tan(), m1.cwise().sin().cwise() / m1.cwise().cos()); - VERIFY_IS_APPROX(mones, m1.cwise().sin().cwise().square() + m1.cwise().cos().cwise().square()); - m3 = m1; - m3.cwise() /= m2; - VERIFY_IS_APPROX(m3, m1.cwise() / m2); - } // check min VERIFY_IS_APPROX( m1.cwise().min(m2), m2.cwise().min(m1) ); @@ -148,8 +163,10 @@ template<typename MatrixType> void cwiseops(const MatrixType& m) VERIFY( (m1.cwise().max(m2).cwise() > (m1-mones)).all() ); VERIFY( (m1.cwise()<m1.unaryExpr(bind2nd(plus<Scalar>(), Scalar(1)))).all() ); - VERIFY( !(m1.cwise()<m1.unaryExpr(bind2nd(minus<Scalar>(), Scalar(1)))).all() ); - VERIFY( !(m1.cwise()>m1.unaryExpr(bind2nd(plus<Scalar>(), Scalar(1)))).any() ); + VERIFY( !(m1.cwise()<m1bis.unaryExpr(bind2nd(minus<Scalar>(), Scalar(1)))).all() ); + VERIFY( !(m1.cwise()>m1bis.unaryExpr(bind2nd(plus<Scalar>(), Scalar(1)))).any() ); + + cwiseops_real_only(m1, m2, m3, mones); } void test_cwiseop() diff --git a/test/denseLM.cpp b/test/denseLM.cpp new file mode 100644 index 000000000..0aa736ea3 --- /dev/null +++ b/test/denseLM.cpp @@ -0,0 +1,190 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Desire Nuentsa <desire.nuentsa_wakam@inria.fr> +// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr> +// +// 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/. + +#include <iostream> +#include <fstream> +#include <iomanip> + +#include "main.h" +#include <Eigen/LevenbergMarquardt> +using namespace std; +using namespace Eigen; + +template<typename Scalar> +struct DenseLM : DenseFunctor<Scalar> +{ + typedef DenseFunctor<Scalar> Base; + typedef typename Base::JacobianType JacobianType; + typedef Matrix<Scalar,Dynamic,1> VectorType; + + DenseLM(int n, int m) : DenseFunctor<Scalar>(n,m) + { } + + VectorType model(const VectorType& uv, VectorType& x) + { + VectorType y; // Should change to use expression template + int m = Base::values(); + int n = Base::inputs(); + eigen_assert(uv.size()%2 == 0); + eigen_assert(uv.size() == n); + eigen_assert(x.size() == m); + y.setZero(m); + int half = n/2; + VectorBlock<const VectorType> u(uv, 0, half); + VectorBlock<const VectorType> v(uv, half, half); + for (int j = 0; j < m; j++) + { + for (int i = 0; i < half; i++) + y(j) += u(i)*std::exp(-(x(j)-i)*(x(j)-i)/(v(i)*v(i))); + } + return y; + + } + void initPoints(VectorType& uv_ref, VectorType& x) + { + m_x = x; + m_y = this->model(uv_ref, x); + } + + int operator()(const VectorType& uv, VectorType& fvec) + { + + int m = Base::values(); + int n = Base::inputs(); + eigen_assert(uv.size()%2 == 0); + eigen_assert(uv.size() == n); + eigen_assert(fvec.size() == m); + int half = n/2; + VectorBlock<const VectorType> u(uv, 0, half); + VectorBlock<const VectorType> v(uv, half, half); + for (int j = 0; j < m; j++) + { + fvec(j) = m_y(j); + for (int i = 0; i < half; i++) + { + fvec(j) -= u(i) *std::exp(-(m_x(j)-i)*(m_x(j)-i)/(v(i)*v(i))); + } + } + + return 0; + } + int df(const VectorType& uv, JacobianType& fjac) + { + int m = Base::values(); + int n = Base::inputs(); + eigen_assert(n == uv.size()); + eigen_assert(fjac.rows() == m); + eigen_assert(fjac.cols() == n); + int half = n/2; + VectorBlock<const VectorType> u(uv, 0, half); + VectorBlock<const VectorType> v(uv, half, half); + for (int j = 0; j < m; j++) + { + for (int i = 0; i < half; i++) + { + fjac.coeffRef(j,i) = -std::exp(-(m_x(j)-i)*(m_x(j)-i)/(v(i)*v(i))); + fjac.coeffRef(j,i+half) = -2.*u(i)*(m_x(j)-i)*(m_x(j)-i)/(std::pow(v(i),3)) * std::exp(-(m_x(j)-i)*(m_x(j)-i)/(v(i)*v(i))); + } + } + return 0; + } + VectorType m_x, m_y; //Data Points +}; + +template<typename FunctorType, typename VectorType> +int test_minimizeLM(FunctorType& functor, VectorType& uv) +{ + LevenbergMarquardt<FunctorType> lm(functor); + LevenbergMarquardtSpace::Status info; + + info = lm.minimize(uv); + + VERIFY_IS_EQUAL(info, 1); + //FIXME Check other parameters + return info; +} + +template<typename FunctorType, typename VectorType> +int test_lmder(FunctorType& functor, VectorType& uv) +{ + typedef typename VectorType::Scalar Scalar; + LevenbergMarquardtSpace::Status info; + LevenbergMarquardt<FunctorType> lm(functor); + info = lm.lmder1(uv); + + VERIFY_IS_EQUAL(info, 1); + //FIXME Check other parameters + return info; +} + +template<typename FunctorType, typename VectorType> +int test_minimizeSteps(FunctorType& functor, VectorType& uv) +{ + LevenbergMarquardtSpace::Status info; + LevenbergMarquardt<FunctorType> lm(functor); + info = lm.minimizeInit(uv); + if (info==LevenbergMarquardtSpace::ImproperInputParameters) + return info; + do + { + info = lm.minimizeOneStep(uv); + } while (info==LevenbergMarquardtSpace::Running); + + VERIFY_IS_EQUAL(info, 1); + //FIXME Check other parameters + return info; +} + +template<typename T> +void test_denseLM_T() +{ + typedef Matrix<T,Dynamic,1> VectorType; + + int inputs = 10; + int values = 1000; + DenseLM<T> dense_gaussian(inputs, values); + VectorType uv(inputs),uv_ref(inputs); + VectorType x(values); + + // Generate the reference solution + uv_ref << -2, 1, 4 ,8, 6, 1.8, 1.2, 1.1, 1.9 , 3; + + //Generate the reference data points + x.setRandom(); + x = 10*x; + x.array() += 10; + dense_gaussian.initPoints(uv_ref, x); + + // Generate the initial parameters + VectorBlock<VectorType> u(uv, 0, inputs/2); + VectorBlock<VectorType> v(uv, inputs/2, inputs/2); + + // Solve the optimization problem + + //Solve in one go + u.setOnes(); v.setOnes(); + test_minimizeLM(dense_gaussian, uv); + + //Solve until the machine precision + u.setOnes(); v.setOnes(); + test_lmder(dense_gaussian, uv); + + // Solve step by step + v.setOnes(); u.setOnes(); + test_minimizeSteps(dense_gaussian, uv); + +} + +void test_denseLM() +{ + CALL_SUBTEST_2(test_denseLM_T<double>()); + + // CALL_SUBTEST_2(test_sparseLM_T<std::complex<double>()); +} diff --git a/test/determinant.cpp b/test/determinant.cpp index e93f2f297..758f3afbb 100644 --- a/test/determinant.cpp +++ b/test/determinant.cpp @@ -39,7 +39,7 @@ template<typename MatrixType> void determinant(const MatrixType& m) m2.col(i).swap(m2.col(j)); VERIFY_IS_APPROX(m2.determinant(), -m1.determinant()); VERIFY_IS_APPROX(m2.determinant(), m2.transpose().determinant()); - VERIFY_IS_APPROX(internal::conj(m2.determinant()), m2.adjoint().determinant()); + VERIFY_IS_APPROX(numext::conj(m2.determinant()), m2.adjoint().determinant()); m2 = m1; m2.row(i) += x*m2.row(j); VERIFY_IS_APPROX(m2.determinant(), m1.determinant()); @@ -53,8 +53,8 @@ template<typename MatrixType> void determinant(const MatrixType& m) void test_determinant() { - int s; for(int i = 0; i < g_repeat; i++) { + int s = 0; CALL_SUBTEST_1( determinant(Matrix<float, 1, 1>()) ); CALL_SUBTEST_2( determinant(Matrix<double, 2, 2>()) ); CALL_SUBTEST_3( determinant(Matrix<double, 3, 3>()) ); @@ -62,6 +62,6 @@ void test_determinant() CALL_SUBTEST_5( determinant(Matrix<std::complex<double>, 10, 10>()) ); s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4); CALL_SUBTEST_6( determinant(MatrixXd(s, s)) ); + TEST_SET_BUT_UNUSED_VARIABLE(s) } - EIGEN_UNUSED_VARIABLE(s) } diff --git a/test/diagonal.cpp b/test/diagonal.cpp index 95cd10372..53814a588 100644 --- a/test/diagonal.cpp +++ b/test/diagonal.cpp @@ -13,9 +13,6 @@ template<typename MatrixType> void diagonal(const MatrixType& m) { typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; - typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; - typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType; Index rows = m.rows(); Index cols = m.cols(); @@ -31,17 +28,14 @@ template<typename MatrixType> void diagonal(const MatrixType& m) if (rows>2) { enum { - N1 = MatrixType::RowsAtCompileTime>1 ? 1 : 0, - N2 = MatrixType::RowsAtCompileTime>2 ? -2 : 0 + N1 = MatrixType::RowsAtCompileTime>2 ? 2 : 0, + N2 = MatrixType::RowsAtCompileTime>1 ? -1 : 0 }; // check sub/super diagonal - if(m1.template diagonal<N1>().RowsAtCompileTime!=Dynamic) + if(MatrixType::SizeAtCompileTime!=Dynamic) { VERIFY(m1.template diagonal<N1>().RowsAtCompileTime == m1.diagonal(N1).size()); - } - if(m1.template diagonal<N2>().RowsAtCompileTime!=Dynamic) - { VERIFY(m1.template diagonal<N2>().RowsAtCompileTime == m1.diagonal(N2).size()); } @@ -56,12 +50,12 @@ template<typename MatrixType> void diagonal(const MatrixType& m) VERIFY_IS_APPROX(m2.template diagonal<N2>()[0], static_cast<Scalar>(6) * m1.template diagonal<N2>()[0]); m2.diagonal(N1) = 2 * m1.diagonal(N1); - VERIFY_IS_APPROX(m2.diagonal<N1>(), static_cast<Scalar>(2) * m1.diagonal(N1)); + VERIFY_IS_APPROX(m2.template diagonal<N1>(), static_cast<Scalar>(2) * m1.diagonal(N1)); m2.diagonal(N1)[0] *= 3; VERIFY_IS_APPROX(m2.diagonal(N1)[0], static_cast<Scalar>(6) * m1.diagonal(N1)[0]); m2.diagonal(N2) = 2 * m1.diagonal(N2); - VERIFY_IS_APPROX(m2.diagonal<N2>(), static_cast<Scalar>(2) * m1.diagonal(N2)); + VERIFY_IS_APPROX(m2.template diagonal<N2>(), static_cast<Scalar>(2) * m1.diagonal(N2)); m2.diagonal(N2)[0] *= 3; VERIFY_IS_APPROX(m2.diagonal(N2)[0], static_cast<Scalar>(6) * m1.diagonal(N2)[0]); } diff --git a/test/diagonalmatrices.cpp b/test/diagonalmatrices.cpp index 3f5776dfc..149f1db2f 100644 --- a/test/diagonalmatrices.cpp +++ b/test/diagonalmatrices.cpp @@ -13,7 +13,6 @@ template<typename MatrixType> void diagonalmatrices(const MatrixType& m) { typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime }; typedef Matrix<Scalar, Rows, 1> VectorType; typedef Matrix<Scalar, 1, Cols> RowVectorType; @@ -32,6 +31,8 @@ template<typename MatrixType> void diagonalmatrices(const MatrixType& m) rv2 = RowVectorType::Random(cols); LeftDiagonalMatrix ldm1(v1), ldm2(v2); RightDiagonalMatrix rdm1(rv1), rdm2(rv2); + + Scalar s1 = internal::random<Scalar>(); SquareMatrixType sq_m1 (v1.asDiagonal()); VERIFY_IS_APPROX(sq_m1, v1.asDiagonal().toDenseMatrix()); @@ -76,6 +77,13 @@ template<typename MatrixType> void diagonalmatrices(const MatrixType& m) big.block(i,j,rows,cols) = big.block(i,j,rows,cols) * rv1.asDiagonal(); VERIFY_IS_APPROX((big.block(i,j,rows,cols)) , m1 * rv1.asDiagonal() ); + + // scalar multiple + VERIFY_IS_APPROX(LeftDiagonalMatrix(ldm1*s1).diagonal(), ldm1.diagonal() * s1); + VERIFY_IS_APPROX(LeftDiagonalMatrix(s1*ldm1).diagonal(), s1 * ldm1.diagonal()); + + VERIFY_IS_APPROX(m1 * (rdm1 * s1), (m1 * rdm1) * s1); + VERIFY_IS_APPROX(m1 * (s1 * rdm1), (m1 * rdm1) * s1); } void test_diagonalmatrices() diff --git a/test/dynalloc.cpp b/test/dynalloc.cpp index c7ccbffef..7e41bfa97 100644 --- a/test/dynalloc.cpp +++ b/test/dynalloc.cpp @@ -53,7 +53,7 @@ void check_aligned_new() void check_aligned_stack_alloc() { - for(int i = 1; i < 1000; i++) + for(int i = 1; i < 400; i++) { ei_declare_aligned_stack_constructed_variable(float,p,i,0); VERIFY(size_t(p)%ALIGNMENT==0); @@ -82,6 +82,7 @@ class MyClassA template<typename T> void check_dynaligned() { T* obj = new T; + VERIFY(T::NeedsToAlign==1); VERIFY(size_t(obj)%ALIGNMENT==0); delete obj; } diff --git a/test/eigen2support.cpp b/test/eigen2support.cpp index 7e02bdf5b..ad1d98091 100644 --- a/test/eigen2support.cpp +++ b/test/eigen2support.cpp @@ -43,7 +43,9 @@ template<typename MatrixType> void eigen2support(const MatrixType& m) VERIFY_IS_EQUAL((m1.col(0).end(1)), (m1.col(0).segment(rows-1,1))); VERIFY_IS_EQUAL((m1.col(0).template end<1>()), (m1.col(0).segment(rows-1,1))); - using namespace internal; + using std::cos; + using numext::real; + using numext::abs2; VERIFY_IS_EQUAL(ei_cos(s1), cos(s1)); VERIFY_IS_EQUAL(ei_real(s1), real(s1)); VERIFY_IS_EQUAL(ei_abs2(s1), abs2(s1)); diff --git a/test/eigensolver_complex.cpp b/test/eigensolver_complex.cpp index 0c2059512..c9d8c0877 100644 --- a/test/eigensolver_complex.cpp +++ b/test/eigensolver_complex.cpp @@ -41,9 +41,6 @@ template<typename MatrixType> void eigensolver(const MatrixType& m) typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; - typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; - typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealVectorType; - typedef typename std::complex<typename NumTraits<typename MatrixType::Scalar>::Real> Complex; MatrixType a = MatrixType::Random(rows,cols); MatrixType symmA = a.adjoint() * a; @@ -59,6 +56,17 @@ template<typename MatrixType> void eigensolver(const MatrixType& m) // another algorithm so results may differ slightly verify_is_approx_upto_permutation(a.eigenvalues(), ei1.eigenvalues()); + ComplexEigenSolver<MatrixType> ei2; + ei2.setMaxIterations(ComplexSchur<MatrixType>::m_maxIterationsPerRow * rows).compute(a); + VERIFY_IS_EQUAL(ei2.info(), Success); + VERIFY_IS_EQUAL(ei2.eigenvectors(), ei1.eigenvectors()); + VERIFY_IS_EQUAL(ei2.eigenvalues(), ei1.eigenvalues()); + if (rows > 2) { + ei2.setMaxIterations(1).compute(a); + VERIFY_IS_EQUAL(ei2.info(), NoConvergence); + VERIFY_IS_EQUAL(ei2.getMaxIterations(), 1); + } + ComplexEigenSolver<MatrixType> eiNoEivecs(a, false); VERIFY_IS_EQUAL(eiNoEivecs.info(), Success); VERIFY_IS_APPROX(ei1.eigenvalues(), eiNoEivecs.eigenvalues()); @@ -93,7 +101,7 @@ template<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m void test_eigensolver_complex() { - int s; + int s = 0; for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1( eigensolver(Matrix4cf()) ); s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4); @@ -101,7 +109,6 @@ void test_eigensolver_complex() CALL_SUBTEST_3( eigensolver(Matrix<std::complex<float>, 1, 1>()) ); CALL_SUBTEST_4( eigensolver(Matrix3f()) ); } - CALL_SUBTEST_1( eigensolver_verify_assert(Matrix4cf()) ); s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4); CALL_SUBTEST_2( eigensolver_verify_assert(MatrixXcd(s,s)) ); @@ -109,7 +116,7 @@ void test_eigensolver_complex() CALL_SUBTEST_4( eigensolver_verify_assert(Matrix3f()) ); // Test problem size constructors - CALL_SUBTEST_5(ComplexEigenSolver<MatrixXf>(s)); + CALL_SUBTEST_5(ComplexEigenSolver<MatrixXf> tmp(s)); - EIGEN_UNUSED_VARIABLE(s) + TEST_SET_BUT_UNUSED_VARIABLE(s) } diff --git a/test/eigensolver_generalized_real.cpp b/test/eigensolver_generalized_real.cpp new file mode 100644 index 000000000..566a4bdc6 --- /dev/null +++ b/test/eigensolver_generalized_real.cpp @@ -0,0 +1,59 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr> +// +// 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/. + +#include "main.h" +#include <limits> +#include <Eigen/Eigenvalues> + +template<typename MatrixType> void generalized_eigensolver_real(const MatrixType& m) +{ + typedef typename MatrixType::Index Index; + /* this test covers the following files: + GeneralizedEigenSolver.h + */ + Index rows = m.rows(); + Index cols = m.cols(); + + typedef typename MatrixType::Scalar Scalar; + typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; + + MatrixType a = MatrixType::Random(rows,cols); + MatrixType b = MatrixType::Random(rows,cols); + MatrixType a1 = MatrixType::Random(rows,cols); + MatrixType b1 = MatrixType::Random(rows,cols); + MatrixType spdA = a.adjoint() * a + a1.adjoint() * a1; + MatrixType spdB = b.adjoint() * b + b1.adjoint() * b1; + + // lets compare to GeneralizedSelfAdjointEigenSolver + GeneralizedSelfAdjointEigenSolver<MatrixType> symmEig(spdA, spdB); + GeneralizedEigenSolver<MatrixType> eig(spdA, spdB); + + VERIFY_IS_EQUAL(eig.eigenvalues().imag().cwiseAbs().maxCoeff(), 0); + + VectorType realEigenvalues = eig.eigenvalues().real(); + std::sort(realEigenvalues.data(), realEigenvalues.data()+realEigenvalues.size()); + VERIFY_IS_APPROX(realEigenvalues, symmEig.eigenvalues()); +} + +void test_eigensolver_generalized_real() +{ + for(int i = 0; i < g_repeat; i++) { + int s = 0; + CALL_SUBTEST_1( generalized_eigensolver_real(Matrix4f()) ); + s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4); + CALL_SUBTEST_2( generalized_eigensolver_real(MatrixXd(s,s)) ); + + // some trivial but implementation-wise tricky cases + CALL_SUBTEST_2( generalized_eigensolver_real(MatrixXd(1,1)) ); + CALL_SUBTEST_2( generalized_eigensolver_real(MatrixXd(2,2)) ); + CALL_SUBTEST_3( generalized_eigensolver_real(Matrix<double,1,1>()) ); + CALL_SUBTEST_4( generalized_eigensolver_real(Matrix2d()) ); + TEST_SET_BUT_UNUSED_VARIABLE(s) + } +} diff --git a/test/eigensolver_generic.cpp b/test/eigensolver_generic.cpp index 0b55ccd93..005af81eb 100644 --- a/test/eigensolver_generic.cpp +++ b/test/eigensolver_generic.cpp @@ -2,7 +2,7 @@ // for linear algebra. // // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> -// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk> +// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk> // // 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 @@ -23,7 +23,6 @@ template<typename MatrixType> void eigensolver(const MatrixType& m) typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; - typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealVectorType; typedef typename std::complex<typename NumTraits<typename MatrixType::Scalar>::Real> Complex; @@ -45,6 +44,17 @@ template<typename MatrixType> void eigensolver(const MatrixType& m) VERIFY_IS_APPROX(ei1.eigenvectors().colwise().norm(), RealVectorType::Ones(rows).transpose()); VERIFY_IS_APPROX(a.eigenvalues(), ei1.eigenvalues()); + EigenSolver<MatrixType> ei2; + ei2.setMaxIterations(RealSchur<MatrixType>::m_maxIterationsPerRow * rows).compute(a); + VERIFY_IS_EQUAL(ei2.info(), Success); + VERIFY_IS_EQUAL(ei2.eigenvectors(), ei1.eigenvectors()); + VERIFY_IS_EQUAL(ei2.eigenvalues(), ei1.eigenvalues()); + if (rows > 2) { + ei2.setMaxIterations(1).compute(a); + VERIFY_IS_EQUAL(ei2.info(), NoConvergence); + VERIFY_IS_EQUAL(ei2.getMaxIterations(), 1); + } + EigenSolver<MatrixType> eiNoEivecs(a, false); VERIFY_IS_EQUAL(eiNoEivecs.info(), Success); VERIFY_IS_APPROX(ei1.eigenvalues(), eiNoEivecs.eigenvalues()); @@ -78,7 +88,7 @@ template<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m void test_eigensolver_generic() { - int s; + int s = 0; for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1( eigensolver(Matrix4f()) ); s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4); @@ -98,7 +108,7 @@ void test_eigensolver_generic() CALL_SUBTEST_4( eigensolver_verify_assert(Matrix2d()) ); // Test problem size constructors - CALL_SUBTEST_5(EigenSolver<MatrixXf>(s)); + CALL_SUBTEST_5(EigenSolver<MatrixXf> tmp(s)); // regression test for bug 410 CALL_SUBTEST_2( @@ -111,5 +121,5 @@ void test_eigensolver_generic() } ); - EIGEN_UNUSED_VARIABLE(s) + TEST_SET_BUT_UNUSED_VARIABLE(s) } diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp index 02dbdb429..5c6ecd875 100644 --- a/test/eigensolver_selfadjoint.cpp +++ b/test/eigensolver_selfadjoint.cpp @@ -23,9 +23,6 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m) typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; - typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; - typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealVectorType; - typedef typename std::complex<typename NumTraits<typename MatrixType::Scalar>::Real> Complex; RealScalar largerEps = 10*test_precision<RealScalar>(); @@ -113,7 +110,7 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m) void test_eigensolver_selfadjoint() { - int s; + int s = 0; for(int i = 0; i < g_repeat; i++) { // very important to test 3x3 and 2x2 matrices since we provide special paths for them CALL_SUBTEST_1( selfadjointeigensolver(Matrix2d()) ); @@ -138,9 +135,9 @@ void test_eigensolver_selfadjoint() // Test problem size constructors s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4); - CALL_SUBTEST_8(SelfAdjointEigenSolver<MatrixXf>(s)); - CALL_SUBTEST_8(Tridiagonalization<MatrixXf>(s)); + CALL_SUBTEST_8(SelfAdjointEigenSolver<MatrixXf> tmp1(s)); + CALL_SUBTEST_8(Tridiagonalization<MatrixXf> tmp2(s)); - EIGEN_UNUSED_VARIABLE(s) + TEST_SET_BUT_UNUSED_VARIABLE(s) } diff --git a/test/exceptions.cpp b/test/exceptions.cpp index 8c48b2f7b..b83fb82ba 100644 --- a/test/exceptions.cpp +++ b/test/exceptions.cpp @@ -69,6 +69,10 @@ class ScalarWithExceptions static int countdown; }; +ScalarWithExceptions real(const ScalarWithExceptions &x) { return x; } +ScalarWithExceptions imag(const ScalarWithExceptions & ) { return 0; } +ScalarWithExceptions conj(const ScalarWithExceptions &x) { return x; } + int ScalarWithExceptions::instances = 0; int ScalarWithExceptions::countdown = 0; diff --git a/test/geo_alignedbox.cpp b/test/geo_alignedbox.cpp index 5886f9181..8e36adbe3 100644 --- a/test/geo_alignedbox.cpp +++ b/test/geo_alignedbox.cpp @@ -15,6 +15,10 @@ #include<iostream> using namespace std; +template<typename T> EIGEN_DONT_INLINE +void kill_extra_precision(T& x) { eigen_assert(&x != 0); } + + template<typename BoxType> void alignedbox(const BoxType& _box) { /* this test covers the following files: @@ -36,6 +40,10 @@ template<typename BoxType> void alignedbox(const BoxType& _box) BoxType b0(dim); BoxType b1(VectorType::Random(dim),VectorType::Random(dim)); BoxType b2; + + kill_extra_precision(b1); + kill_extra_precision(p0); + kill_extra_precision(p1); b0.extend(p0); b0.extend(p1); @@ -71,7 +79,6 @@ void alignedboxCastTests(const BoxType& _box) // casting typedef typename BoxType::Index Index; typedef typename BoxType::Scalar Scalar; - typedef typename NumTraits<Scalar>::Real RealScalar; typedef Matrix<Scalar, BoxType::AmbientDimAtCompileTime, 1> VectorType; const Index dim = _box.dim(); @@ -109,7 +116,7 @@ void specificTest1() VERIFY_IS_APPROX( 14.0f, box.volume() ); VERIFY_IS_APPROX( 53.0f, box.diagonal().squaredNorm() ); - VERIFY_IS_APPROX( internal::sqrt( 53.0f ), box.diagonal().norm() ); + VERIFY_IS_APPROX( std::sqrt( 53.0f ), box.diagonal().norm() ); VERIFY_IS_APPROX( m, box.corner( BoxType::BottomLeft ) ); VERIFY_IS_APPROX( M, box.corner( BoxType::TopRight ) ); diff --git a/test/geo_eulerangles.cpp b/test/geo_eulerangles.cpp index 9bf149d2a..b4830bd41 100644 --- a/test/geo_eulerangles.cpp +++ b/test/geo_eulerangles.cpp @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2008-2012 Gael Guennebaud <gael.guennebaud@inria.fr> // // 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 @@ -12,10 +12,55 @@ #include <Eigen/LU> #include <Eigen/SVD> -template<typename Scalar> void eulerangles(void) + +template<typename Scalar> +void verify_euler(const Matrix<Scalar,3,1>& ea, int i, int j, int k) { typedef Matrix<Scalar,3,3> Matrix3; typedef Matrix<Scalar,3,1> Vector3; + typedef AngleAxis<Scalar> AngleAxisx; + using std::abs; + Matrix3 m(AngleAxisx(ea[0], Vector3::Unit(i)) * AngleAxisx(ea[1], Vector3::Unit(j)) * AngleAxisx(ea[2], Vector3::Unit(k))); + Vector3 eabis = m.eulerAngles(i, j, k); + Matrix3 mbis(AngleAxisx(eabis[0], Vector3::Unit(i)) * AngleAxisx(eabis[1], Vector3::Unit(j)) * AngleAxisx(eabis[2], Vector3::Unit(k))); + VERIFY_IS_APPROX(m, mbis); + /* If I==K, and ea[1]==0, then there no unique solution. */ + /* The remark apply in the case where I!=K, and |ea[1]| is close to pi/2. */ + if( (i!=k || ea[1]!=0) && (i==k || !internal::isApprox(abs(ea[1]),Scalar(M_PI/2),test_precision<Scalar>())) ) + VERIFY((ea-eabis).norm() <= test_precision<Scalar>()); + + // approx_or_less_than does not work for 0 + VERIFY(0 < eabis[0] || test_isMuchSmallerThan(eabis[0], Scalar(1))); + VERIFY_IS_APPROX_OR_LESS_THAN(eabis[0], Scalar(M_PI)); + VERIFY_IS_APPROX_OR_LESS_THAN(-Scalar(M_PI), eabis[1]); + VERIFY_IS_APPROX_OR_LESS_THAN(eabis[1], Scalar(M_PI)); + VERIFY_IS_APPROX_OR_LESS_THAN(-Scalar(M_PI), eabis[2]); + VERIFY_IS_APPROX_OR_LESS_THAN(eabis[2], Scalar(M_PI)); +} + +template<typename Scalar> void check_all_var(const Matrix<Scalar,3,1>& ea) +{ + verify_euler(ea, 0,1,2); + verify_euler(ea, 0,1,0); + verify_euler(ea, 0,2,1); + verify_euler(ea, 0,2,0); + + verify_euler(ea, 1,2,0); + verify_euler(ea, 1,2,1); + verify_euler(ea, 1,0,2); + verify_euler(ea, 1,0,1); + + verify_euler(ea, 2,0,1); + verify_euler(ea, 2,0,2); + verify_euler(ea, 2,1,0); + verify_euler(ea, 2,1,2); +} + +template<typename Scalar> void eulerangles() +{ + typedef Matrix<Scalar,3,3> Matrix3; + typedef Matrix<Scalar,3,1> Vector3; + typedef Array<Scalar,3,1> Array3; typedef Quaternion<Scalar> Quaternionx; typedef AngleAxis<Scalar> AngleAxisx; @@ -24,25 +69,38 @@ template<typename Scalar> void eulerangles(void) q1 = AngleAxisx(a, Vector3::Random().normalized()); Matrix3 m; m = q1; - - #define VERIFY_EULER(I,J,K, X,Y,Z) { \ - Vector3 ea = m.eulerAngles(I,J,K); \ - VERIFY_IS_APPROX(m, Matrix3(AngleAxisx(ea[0], Vector3::Unit##X()) * AngleAxisx(ea[1], Vector3::Unit##Y()) * AngleAxisx(ea[2], Vector3::Unit##Z()))); \ - } - VERIFY_EULER(0,1,2, X,Y,Z); - VERIFY_EULER(0,1,0, X,Y,X); - VERIFY_EULER(0,2,1, X,Z,Y); - VERIFY_EULER(0,2,0, X,Z,X); - - VERIFY_EULER(1,2,0, Y,Z,X); - VERIFY_EULER(1,2,1, Y,Z,Y); - VERIFY_EULER(1,0,2, Y,X,Z); - VERIFY_EULER(1,0,1, Y,X,Y); - - VERIFY_EULER(2,0,1, Z,X,Y); - VERIFY_EULER(2,0,2, Z,X,Z); - VERIFY_EULER(2,1,0, Z,Y,X); - VERIFY_EULER(2,1,2, Z,Y,Z); + + Vector3 ea = m.eulerAngles(0,1,2); + check_all_var(ea); + ea = m.eulerAngles(0,1,0); + check_all_var(ea); + + // Check with purely random Quaternion: + q1.coeffs() = Quaternionx::Coefficients::Random().normalized(); + m = q1; + ea = m.eulerAngles(0,1,2); + check_all_var(ea); + ea = m.eulerAngles(0,1,0); + check_all_var(ea); + + // Check with random angles in range [0:pi]x[-pi:pi]x[-pi:pi]. + ea = (Array3::Random() + Array3(1,0,0))*Scalar(M_PI)*Array3(0.5,1,1); + check_all_var(ea); + + ea[2] = ea[0] = internal::random<Scalar>(0,Scalar(M_PI)); + check_all_var(ea); + + ea[0] = ea[1] = internal::random<Scalar>(0,Scalar(M_PI)); + check_all_var(ea); + + ea[1] = 0; + check_all_var(ea); + + ea.head(2).setZero(); + check_all_var(ea); + + ea.setZero(); + check_all_var(ea); } void test_geo_eulerangles() diff --git a/test/geo_hyperplane.cpp b/test/geo_hyperplane.cpp index 3fc80c4c7..f26fc1329 100644 --- a/test/geo_hyperplane.cpp +++ b/test/geo_hyperplane.cpp @@ -22,7 +22,6 @@ template<typename HyperplaneType> void hyperplane(const HyperplaneType& _plane) const Index dim = _plane.dim(); enum { Options = HyperplaneType::Options }; typedef typename HyperplaneType::Scalar Scalar; - typedef typename NumTraits<Scalar>::Real RealScalar; typedef Matrix<Scalar, HyperplaneType::AmbientDimAtCompileTime, 1> VectorType; typedef Matrix<Scalar, HyperplaneType::AmbientDimAtCompileTime, HyperplaneType::AmbientDimAtCompileTime> MatrixType; @@ -79,6 +78,7 @@ template<typename HyperplaneType> void hyperplane(const HyperplaneType& _plane) template<typename Scalar> void lines() { + using std::abs; typedef Hyperplane<Scalar, 2> HLine; typedef ParametrizedLine<Scalar, 2> PLine; typedef Matrix<Scalar,2,1> Vector; @@ -90,7 +90,7 @@ template<typename Scalar> void lines() Vector u = Vector::Random(); Vector v = Vector::Random(); Scalar a = internal::random<Scalar>(); - while (internal::abs(a-1) < 1e-4) a = internal::random<Scalar>(); + while (abs(a-1) < 1e-4) a = internal::random<Scalar>(); while (u.norm() < 1e-4) u = Vector::Random(); while (v.norm() < 1e-4) v = Vector::Random(); diff --git a/test/geo_parametrizedline.cpp b/test/geo_parametrizedline.cpp index 4e1f845dd..f0462d40a 100644 --- a/test/geo_parametrizedline.cpp +++ b/test/geo_parametrizedline.cpp @@ -18,13 +18,12 @@ template<typename LineType> void parametrizedline(const LineType& _line) /* this test covers the following files: ParametrizedLine.h */ + using std::abs; typedef typename LineType::Index Index; const Index dim = _line.dim(); typedef typename LineType::Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; typedef Matrix<Scalar, LineType::AmbientDimAtCompileTime, 1> VectorType; - typedef Matrix<Scalar, LineType::AmbientDimAtCompileTime, - LineType::AmbientDimAtCompileTime> MatrixType; typedef Hyperplane<Scalar,LineType::AmbientDimAtCompileTime> HyperplaneType; VectorType p0 = VectorType::Random(dim); @@ -35,7 +34,7 @@ template<typename LineType> void parametrizedline(const LineType& _line) LineType l0(p0, d0); Scalar s0 = internal::random<Scalar>(); - Scalar s1 = internal::abs(internal::random<Scalar>()); + Scalar s1 = abs(internal::random<Scalar>()); VERIFY_IS_MUCH_SMALLER_THAN( l0.distance(p0), RealScalar(1) ); VERIFY_IS_MUCH_SMALLER_THAN( l0.distance(p0+s0*d0), RealScalar(1) ); diff --git a/test/geo_quaternion.cpp b/test/geo_quaternion.cpp index 6e6922864..1694b32c7 100644 --- a/test/geo_quaternion.cpp +++ b/test/geo_quaternion.cpp @@ -23,22 +23,22 @@ template<typename T> T bounded_acos(T v) template<typename QuatType> void check_slerp(const QuatType& q0, const QuatType& q1) { + using std::abs; typedef typename QuatType::Scalar Scalar; - typedef Matrix<Scalar,3,1> VectorType; typedef AngleAxis<Scalar> AA; Scalar largeEps = test_precision<Scalar>(); Scalar theta_tot = AA(q1*q0.inverse()).angle(); if(theta_tot>M_PI) - theta_tot = 2.*M_PI-theta_tot; - for(Scalar t=0; t<=1.001; t+=0.1) + theta_tot = Scalar(2.*M_PI)-theta_tot; + for(Scalar t=0; t<=Scalar(1.001); t+=Scalar(0.1)) { QuatType q = q0.slerp(t,q1); Scalar theta = AA(q*q0.inverse()).angle(); - VERIFY(internal::abs(q.norm() - 1) < largeEps); + VERIFY(abs(q.norm() - 1) < largeEps); if(theta_tot==0) VERIFY(theta_tot==0); - else VERIFY(internal::abs(theta/theta_tot - t) < largeEps); + else VERIFY(abs(theta - t * theta_tot) < largeEps); } } @@ -47,8 +47,7 @@ template<typename Scalar, int Options> void quaternion(void) /* this test covers the following files: Quaternion.h */ - - typedef Matrix<Scalar,3,3> Matrix3; + using std::abs; typedef Matrix<Scalar,3,1> Vector3; typedef Matrix<Scalar,4,1> Vector4; typedef Quaternion<Scalar,Options> Quaternionx; @@ -82,13 +81,13 @@ template<typename Scalar, int Options> void quaternion(void) q2 = AngleAxisx(a, v1.normalized()); // angular distance - Scalar refangle = internal::abs(AngleAxisx(q1.inverse()*q2).angle()); + Scalar refangle = abs(AngleAxisx(q1.inverse()*q2).angle()); if (refangle>Scalar(M_PI)) refangle = Scalar(2)*Scalar(M_PI) - refangle; if((q1.coeffs()-q2.coeffs()).norm() > 10*largeEps) { - VERIFY_IS_MUCH_SMALLER_THAN(internal::abs(q1.angularDistance(q2) - refangle), Scalar(1)); + VERIFY_IS_MUCH_SMALLER_THAN(abs(q1.angularDistance(q2) - refangle), Scalar(1)); } // rotation matrix conversion @@ -109,7 +108,7 @@ template<typename Scalar, int Options> void quaternion(void) // Do not execute the test if the rotation angle is almost zero, or // the rotation axis and v1 are almost parallel. - if (internal::abs(aa.angle()) > 5*test_precision<Scalar>() + if (abs(aa.angle()) > 5*test_precision<Scalar>() && (aa.axis() - v1.normalized()).norm() < 1.99 && (aa.axis() + v1.normalized()).norm() < 1.99) { @@ -157,7 +156,7 @@ template<typename Scalar, int Options> void quaternion(void) check_slerp(q1,q2); q1 = AngleAxisx(b, v1.normalized()); - q2 = AngleAxisx(b+M_PI, v1.normalized()); + q2 = AngleAxisx(b+Scalar(M_PI), v1.normalized()); check_slerp(q1,q2); q1 = AngleAxisx(b, v1.normalized()); @@ -171,23 +170,36 @@ template<typename Scalar, int Options> void quaternion(void) template<typename Scalar> void mapQuaternion(void){ typedef Map<Quaternion<Scalar>, Aligned> MQuaternionA; + typedef Map<const Quaternion<Scalar>, Aligned> MCQuaternionA; typedef Map<Quaternion<Scalar> > MQuaternionUA; typedef Map<const Quaternion<Scalar> > MCQuaternionUA; typedef Quaternion<Scalar> Quaternionx; + typedef Matrix<Scalar,3,1> Vector3; + typedef AngleAxis<Scalar> AngleAxisx; + + Vector3 v0 = Vector3::Random(), + v1 = Vector3::Random(); + Scalar a = internal::random<Scalar>(-Scalar(M_PI), Scalar(M_PI)); EIGEN_ALIGN16 Scalar array1[4]; EIGEN_ALIGN16 Scalar array2[4]; EIGEN_ALIGN16 Scalar array3[4+1]; Scalar* array3unaligned = array3+1; + + MQuaternionA mq1(array1); + MCQuaternionA mcq1(array1); + MQuaternionA mq2(array2); + MQuaternionUA mq3(array3unaligned); + MCQuaternionUA mcq3(array3unaligned); // std::cerr << array1 << " " << array2 << " " << array3 << "\n"; - MQuaternionA(array1).coeffs().setRandom(); - (MQuaternionA(array2)) = MQuaternionA(array1); - (MQuaternionUA(array3unaligned)) = MQuaternionA(array1); + mq1 = AngleAxisx(a, v0.normalized()); + mq2 = mq1; + mq3 = mq1; - Quaternionx q1 = MQuaternionA(array1); - Quaternionx q2 = MQuaternionA(array2); - Quaternionx q3 = MQuaternionUA(array3unaligned); + Quaternionx q1 = mq1; + Quaternionx q2 = mq2; + Quaternionx q3 = mq3; Quaternionx q4 = MCQuaternionUA(array3unaligned); VERIFY_IS_APPROX(q1.coeffs(), q2.coeffs()); @@ -197,6 +209,23 @@ template<typename Scalar> void mapQuaternion(void){ if(internal::packet_traits<Scalar>::Vectorizable) VERIFY_RAISES_ASSERT((MQuaternionA(array3unaligned))); #endif + + VERIFY_IS_APPROX(mq1 * (mq1.inverse() * v1), v1); + VERIFY_IS_APPROX(mq1 * (mq1.conjugate() * v1), v1); + + VERIFY_IS_APPROX(mcq1 * (mcq1.inverse() * v1), v1); + VERIFY_IS_APPROX(mcq1 * (mcq1.conjugate() * v1), v1); + + VERIFY_IS_APPROX(mq3 * (mq3.inverse() * v1), v1); + VERIFY_IS_APPROX(mq3 * (mq3.conjugate() * v1), v1); + + VERIFY_IS_APPROX(mcq3 * (mcq3.inverse() * v1), v1); + VERIFY_IS_APPROX(mcq3 * (mcq3.conjugate() * v1), v1); + + VERIFY_IS_APPROX(mq1*mq2, q1*q2); + VERIFY_IS_APPROX(mq3*mq2, q3*q2); + VERIFY_IS_APPROX(mcq1*mq2, q1*q2); + VERIFY_IS_APPROX(mcq3*mq2, q3*q2); } template<typename Scalar> void quaternionAlignment(void){ diff --git a/test/geo_transformations.cpp b/test/geo_transformations.cpp index f4d65aabc..ee3030b5d 100644 --- a/test/geo_transformations.cpp +++ b/test/geo_transformations.cpp @@ -17,22 +17,11 @@ template<typename Scalar, int Mode, int Options> void non_projective_only() /* this test covers the following files: Cross.h Quaternion.h, Transform.cpp */ - typedef Matrix<Scalar,2,2> Matrix2; - typedef Matrix<Scalar,3,3> Matrix3; - typedef Matrix<Scalar,4,4> Matrix4; - typedef Matrix<Scalar,2,1> Vector2; typedef Matrix<Scalar,3,1> Vector3; - typedef Matrix<Scalar,4,1> Vector4; typedef Quaternion<Scalar> Quaternionx; typedef AngleAxis<Scalar> AngleAxisx; - typedef Transform<Scalar,2,Mode,Options> Transform2; typedef Transform<Scalar,3,Mode,Options> Transform3; - typedef Transform<Scalar,2,Isometry,Options> Isometry2; - typedef Transform<Scalar,3,Isometry,Options> Isometry3; - typedef typename Transform3::MatrixType MatrixType; - typedef DiagonalMatrix<Scalar,2> AlignedScaling2; typedef DiagonalMatrix<Scalar,3> AlignedScaling3; - typedef Translation<Scalar,2> Translation2; typedef Translation<Scalar,3> Translation3; Vector3 v0 = Vector3::Random(), @@ -88,7 +77,8 @@ template<typename Scalar, int Mode, int Options> void transformations() /* this test covers the following files: Cross.h Quaternion.h, Transform.cpp */ - typedef Matrix<Scalar,2,2> Matrix2; + using std::cos; + using std::abs; typedef Matrix<Scalar,3,3> Matrix3; typedef Matrix<Scalar,4,4> Matrix4; typedef Matrix<Scalar,2,1> Vector2; @@ -98,10 +88,7 @@ template<typename Scalar, int Mode, int Options> void transformations() typedef AngleAxis<Scalar> AngleAxisx; typedef Transform<Scalar,2,Mode,Options> Transform2; typedef Transform<Scalar,3,Mode,Options> Transform3; - typedef Transform<Scalar,2,Isometry,Options> Isometry2; - typedef Transform<Scalar,3,Isometry,Options> Isometry3; typedef typename Transform3::MatrixType MatrixType; - typedef DiagonalMatrix<Scalar,2> AlignedScaling2; typedef DiagonalMatrix<Scalar,3> AlignedScaling3; typedef Translation<Scalar,2> Translation2; typedef Translation<Scalar,3> Translation3; @@ -115,7 +102,7 @@ template<typename Scalar, int Mode, int Options> void transformations() VERIFY_IS_APPROX(v0, AngleAxisx(a, v0.normalized()) * v0); VERIFY_IS_APPROX(-v0, AngleAxisx(Scalar(M_PI), v0.unitOrthogonal()) * v0); - VERIFY_IS_APPROX(internal::cos(a)*v0.squaredNorm(), v0.dot(AngleAxisx(a, v0.unitOrthogonal()) * v0)); + VERIFY_IS_APPROX(cos(a)*v0.squaredNorm(), v0.dot(AngleAxisx(a, v0.unitOrthogonal()) * v0)); m = AngleAxisx(a, v0.normalized()).toRotationMatrix().adjoint(); VERIFY_IS_APPROX(Matrix3::Identity(), m * AngleAxisx(a, v0.normalized())); VERIFY_IS_APPROX(Matrix3::Identity(), AngleAxisx(a, v0.normalized()) * m); @@ -155,7 +142,7 @@ template<typename Scalar, int Mode, int Options> void transformations() // Transform // TODO complete the tests ! a = 0; - while (internal::abs(a)<Scalar(0.1)) + while (abs(a)<Scalar(0.1)) a = internal::random<Scalar>(-Scalar(0.4)*Scalar(M_PI), Scalar(0.4)*Scalar(M_PI)); q1 = AngleAxisx(a, v0.normalized()); Transform3 t0, t1, t2; @@ -249,7 +236,7 @@ template<typename Scalar, int Mode, int Options> void transformations() Vector2 v20 = Vector2::Random(); Vector2 v21 = Vector2::Random(); for (int k=0; k<2; ++k) - if (internal::abs(v21[k])<Scalar(1e-3)) v21[k] = Scalar(1e-3); + if (abs(v21[k])<Scalar(1e-3)) v21[k] = Scalar(1e-3); t21.setIdentity(); t21.linear() = Rotation2D<Scalar>(a).toRotationMatrix(); VERIFY_IS_APPROX(t20.fromPositionOrientationScale(v20,a,v21).matrix(), @@ -292,6 +279,13 @@ template<typename Scalar, int Mode, int Options> void transformations() t1 = Eigen::Scaling(s0,s0,s0) * t1; VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + t0 = t3; + t0.scale(s0); + t1 = t3 * Eigen::Scaling(s0); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + t0.prescale(s0); + t1 = Eigen::Scaling(s0) * t1; + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); t0.setIdentity(); t0.prerotate(q1).prescale(v0).pretranslate(v0); @@ -400,8 +394,8 @@ template<typename Scalar, int Mode, int Options> void transformations() Rotation2D<double> r2d1d = r2d1.template cast<double>(); VERIFY_IS_APPROX(r2d1d.template cast<Scalar>(),r2d1); - t20 = Translation2(v20) * (Rotation2D<Scalar>(s0) * Scaling(s0)); - t21 = Translation2(v20) * Rotation2D<Scalar>(s0) * Scaling(s0); + t20 = Translation2(v20) * (Rotation2D<Scalar>(s0) * Eigen::Scaling(s0)); + t21 = Translation2(v20) * Rotation2D<Scalar>(s0) * Eigen::Scaling(s0); VERIFY_IS_APPROX(t20,t21); } diff --git a/test/householder.cpp b/test/householder.cpp index 203dce46c..c5f6b5e4f 100644 --- a/test/householder.cpp +++ b/test/householder.cpp @@ -29,8 +29,6 @@ template<typename MatrixType> void householder(const MatrixType& m) typedef Matrix<Scalar, Dynamic, MatrixType::ColsAtCompileTime> HBlockMatrixType; typedef Matrix<Scalar, Dynamic, 1> HCoeffsVectorType; - typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> RightSquareMatrixType; - typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, Dynamic> VBlockMatrixType; typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime> TMatrixType; Matrix<Scalar, EIGEN_SIZE_MAX(MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime), 1> _tmp((std::max)(rows,cols)); @@ -62,8 +60,8 @@ template<typename MatrixType> void householder(const MatrixType& m) m1.applyHouseholderOnTheLeft(essential,beta,tmp); VERIFY_IS_APPROX(m1.norm(), m2.norm()); if(rows>=2) VERIFY_IS_MUCH_SMALLER_THAN(m1.block(1,0,rows-1,cols).norm(), m1.norm()); - VERIFY_IS_MUCH_SMALLER_THAN(internal::imag(m1(0,0)), internal::real(m1(0,0))); - VERIFY_IS_APPROX(internal::real(m1(0,0)), alpha); + VERIFY_IS_MUCH_SMALLER_THAN(numext::imag(m1(0,0)), numext::real(m1(0,0))); + VERIFY_IS_APPROX(numext::real(m1(0,0)), alpha); v1 = VectorType::Random(rows); if(even) v1.tail(rows-1).setZero(); @@ -74,8 +72,8 @@ template<typename MatrixType> void householder(const MatrixType& m) m3.applyHouseholderOnTheRight(essential,beta,tmp); VERIFY_IS_APPROX(m3.norm(), m4.norm()); if(rows>=2) VERIFY_IS_MUCH_SMALLER_THAN(m3.block(0,1,rows,rows-1).norm(), m3.norm()); - VERIFY_IS_MUCH_SMALLER_THAN(internal::imag(m3(0,0)), internal::real(m3(0,0))); - VERIFY_IS_APPROX(internal::real(m3(0,0)), alpha); + VERIFY_IS_MUCH_SMALLER_THAN(numext::imag(m3(0,0)), numext::real(m3(0,0))); + VERIFY_IS_APPROX(numext::real(m3(0,0)), alpha); // test householder sequence on the left with a shift @@ -91,12 +89,29 @@ template<typename MatrixType> void householder(const MatrixType& m) hseq.setLength(hc.size()).setShift(shift); VERIFY(hseq.length() == hc.size()); VERIFY(hseq.shift() == shift); - + MatrixType m5 = m2; m5.block(shift,0,brows,cols).template triangularView<StrictlyLower>().setZero(); VERIFY_IS_APPROX(hseq * m5, m1); // test applying hseq directly m3 = hseq; VERIFY_IS_APPROX(m3 * m5, m1); // test evaluating hseq to a dense matrix, then applying + + SquareMatrixType hseq_mat = hseq; + SquareMatrixType hseq_mat_conj = hseq.conjugate(); + SquareMatrixType hseq_mat_adj = hseq.adjoint(); + SquareMatrixType hseq_mat_trans = hseq.transpose(); + SquareMatrixType m6 = SquareMatrixType::Random(rows, rows); + VERIFY_IS_APPROX(hseq_mat.adjoint(), hseq_mat_adj); + VERIFY_IS_APPROX(hseq_mat.conjugate(), hseq_mat_conj); + VERIFY_IS_APPROX(hseq_mat.transpose(), hseq_mat_trans); + VERIFY_IS_APPROX(hseq_mat * m6, hseq_mat * m6); + VERIFY_IS_APPROX(hseq_mat.adjoint() * m6, hseq_mat_adj * m6); + VERIFY_IS_APPROX(hseq_mat.conjugate() * m6, hseq_mat_conj * m6); + VERIFY_IS_APPROX(hseq_mat.transpose() * m6, hseq_mat_trans * m6); + VERIFY_IS_APPROX(m6 * hseq_mat, m6 * hseq_mat); + VERIFY_IS_APPROX(m6 * hseq_mat.adjoint(), m6 * hseq_mat_adj); + VERIFY_IS_APPROX(m6 * hseq_mat.conjugate(), m6 * hseq_mat_conj); + VERIFY_IS_APPROX(m6 * hseq_mat.transpose(), m6 * hseq_mat_trans); // test householder sequence on the right with a shift diff --git a/test/inverse.cpp b/test/inverse.cpp index cff42dd8d..8187b088d 100644 --- a/test/inverse.cpp +++ b/test/inverse.cpp @@ -13,6 +13,7 @@ template<typename MatrixType> void inverse(const MatrixType& m) { + using std::abs; typedef typename MatrixType::Index Index; /* this test covers the following files: Inverse.h @@ -21,8 +22,6 @@ template<typename MatrixType> void inverse(const MatrixType& m) Index cols = m.cols(); typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits<Scalar>::Real RealScalar; - typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType; MatrixType m1(rows, cols), m2(rows, cols), @@ -42,6 +41,9 @@ template<typename MatrixType> void inverse(const MatrixType& m) VERIFY_IS_APPROX(MatrixType(m1.transpose().inverse()), MatrixType(m1.inverse().transpose())); #if !defined(EIGEN_TEST_PART_5) && !defined(EIGEN_TEST_PART_6) + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType; + //computeInverseAndDetWithCheck tests //First: an invertible matrix bool invertible; @@ -63,7 +65,7 @@ template<typename MatrixType> void inverse(const MatrixType& m) MatrixType m3 = v3*v3.transpose(), m4(rows,cols); m3.computeInverseAndDetWithCheck(m4, det, invertible); VERIFY( rows==1 ? invertible : !invertible ); - VERIFY_IS_MUCH_SMALLER_THAN(internal::abs(det-m3.determinant()), RealScalar(1)); + VERIFY_IS_MUCH_SMALLER_THAN(abs(det-m3.determinant()), RealScalar(1)); m3.computeInverseWithCheck(m4, invertible); VERIFY( rows==1 ? invertible : !invertible ); #endif @@ -84,19 +86,19 @@ template<typename MatrixType> void inverse(const MatrixType& m) void test_inverse() { - int s; + int s = 0; for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1( inverse(Matrix<double,1,1>()) ); CALL_SUBTEST_2( inverse(Matrix2d()) ); CALL_SUBTEST_3( inverse(Matrix3f()) ); CALL_SUBTEST_4( inverse(Matrix4f()) ); CALL_SUBTEST_4( inverse(Matrix<float,4,4,DontAlign>()) ); - s = internal::random<int>(50,320); + s = internal::random<int>(50,320); CALL_SUBTEST_5( inverse(MatrixXf(s,s)) ); s = internal::random<int>(25,100); CALL_SUBTEST_6( inverse(MatrixXcd(s,s)) ); CALL_SUBTEST_7( inverse(Matrix4d()) ); CALL_SUBTEST_7( inverse(Matrix<double,4,4,DontAlign>()) ); } - EIGEN_UNUSED_VARIABLE(s) + TEST_SET_BUT_UNUSED_VARIABLE(s) } diff --git a/test/jacobi.cpp b/test/jacobi.cpp index f64f5d08f..7ccd4124b 100644 --- a/test/jacobi.cpp +++ b/test/jacobi.cpp @@ -14,7 +14,6 @@ template<typename MatrixType, typename JacobiScalar> void jacobi(const MatrixType& m = MatrixType()) { - typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Index Index; Index rows = m.rows(); Index cols = m.cols(); @@ -41,8 +40,8 @@ void jacobi(const MatrixType& m = MatrixType()) MatrixType b = a; b.applyOnTheLeft(p, q, rot); - VERIFY_IS_APPROX(b.row(p), c * a.row(p) + internal::conj(s) * a.row(q)); - VERIFY_IS_APPROX(b.row(q), -s * a.row(p) + internal::conj(c) * a.row(q)); + VERIFY_IS_APPROX(b.row(p), c * a.row(p) + numext::conj(s) * a.row(q)); + VERIFY_IS_APPROX(b.row(q), -s * a.row(p) + numext::conj(c) * a.row(q)); } { @@ -55,7 +54,7 @@ void jacobi(const MatrixType& m = MatrixType()) MatrixType b = a; b.applyOnTheRight(p, q, rot); VERIFY_IS_APPROX(b.col(p), c * a.col(p) - s * a.col(q)); - VERIFY_IS_APPROX(b.col(q), internal::conj(s) * a.col(p) + internal::conj(c) * a.col(q)); + VERIFY_IS_APPROX(b.col(q), numext::conj(s) * a.col(p) + numext::conj(c) * a.col(q)); } } @@ -75,7 +74,8 @@ void test_jacobi() // complex<float> is really important to test as it is the only way to cover conjugation issues in certain unaligned paths CALL_SUBTEST_6(( jacobi<MatrixXcf, float>(MatrixXcf(r,c)) )); CALL_SUBTEST_6(( jacobi<MatrixXcf, std::complex<float> >(MatrixXcf(r,c)) )); - (void) r; - (void) c; + + TEST_SET_BUT_UNUSED_VARIABLE(r); + TEST_SET_BUT_UNUSED_VARIABLE(c); } } diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp index f6c567829..7c21f0ab3 100644 --- a/test/jacobisvd.cpp +++ b/test/jacobisvd.cpp @@ -27,11 +27,8 @@ void jacobisvd_check_full(const MatrixType& m, const JacobiSVD<MatrixType, QRPre }; typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits<Scalar>::Real RealScalar; typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixUType; typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime> MatrixVType; - typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType; - typedef Matrix<Scalar, ColsAtCompileTime, 1> InputVectorType; MatrixType sigma = MatrixType::Zero(rows,cols); sigma.diagonal() = svd.singularValues().template cast<Scalar>(); @@ -70,6 +67,7 @@ template<typename MatrixType, int QRPreconditioner> void jacobisvd_solve(const MatrixType& m, unsigned int computationOptions) { typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::Index Index; Index rows = m.rows(); Index cols = m.cols(); @@ -84,9 +82,90 @@ void jacobisvd_solve(const MatrixType& m, unsigned int computationOptions) RhsType rhs = RhsType::Random(rows, internal::random<Index>(1, cols)); JacobiSVD<MatrixType, QRPreconditioner> svd(m, computationOptions); + + if(internal::is_same<RealScalar,double>::value) svd.setThreshold(1e-8); + else if(internal::is_same<RealScalar,float>::value) svd.setThreshold(1e-4); + SolutionType x = svd.solve(rhs); + + RealScalar residual = (m*x-rhs).norm(); + // Check that there is no significantly better solution in the neighborhood of x + if(!test_isMuchSmallerThan(residual,rhs.norm())) + { + // If the residual is very small, then we have an exact solution, so we are already good. + for(int k=0;k<x.rows();++k) + { + SolutionType y(x); + y.row(k).array() += 2*NumTraits<RealScalar>::epsilon(); + RealScalar residual_y = (m*y-rhs).norm(); + VERIFY( test_isApprox(residual_y,residual) || residual < residual_y ); + + y.row(k) = x.row(k).array() - 2*NumTraits<RealScalar>::epsilon(); + residual_y = (m*y-rhs).norm(); + VERIFY( test_isApprox(residual_y,residual) || residual < residual_y ); + } + } + // evaluate normal equation which works also for least-squares solutions - VERIFY_IS_APPROX(m.adjoint()*m*x,m.adjoint()*rhs); + if(internal::is_same<RealScalar,double>::value) + { + // This test is not stable with single precision. + // This is probably because squaring m signicantly affects the precision. + VERIFY_IS_APPROX(m.adjoint()*m*x,m.adjoint()*rhs); + } + + // check minimal norm solutions + { + // generate a full-rank m x n problem with m<n + enum { + RankAtCompileTime2 = ColsAtCompileTime==Dynamic ? Dynamic : (ColsAtCompileTime)/2+1, + RowsAtCompileTime3 = ColsAtCompileTime==Dynamic ? Dynamic : ColsAtCompileTime+1 + }; + typedef Matrix<Scalar, RankAtCompileTime2, ColsAtCompileTime> MatrixType2; + typedef Matrix<Scalar, RankAtCompileTime2, 1> RhsType2; + typedef Matrix<Scalar, ColsAtCompileTime, RankAtCompileTime2> MatrixType2T; + Index rank = RankAtCompileTime2==Dynamic ? internal::random<Index>(1,cols) : Index(RankAtCompileTime2); + MatrixType2 m2(rank,cols); + int guard = 0; + do { + m2.setRandom(); + } while(m2.jacobiSvd().setThreshold(test_precision<Scalar>()).rank()!=rank && (++guard)<10); + VERIFY(guard<10); + RhsType2 rhs2 = RhsType2::Random(rank); + // use QR to find a reference minimal norm solution + HouseholderQR<MatrixType2T> qr(m2.adjoint()); + Matrix<Scalar,Dynamic,1> tmp = qr.matrixQR().topLeftCorner(rank,rank).template triangularView<Upper>().adjoint().solve(rhs2); + tmp.conservativeResize(cols); + tmp.tail(cols-rank).setZero(); + SolutionType x21 = qr.householderQ() * tmp; + // now check with SVD + JacobiSVD<MatrixType2, ColPivHouseholderQRPreconditioner> svd2(m2, computationOptions); + SolutionType x22 = svd2.solve(rhs2); + VERIFY_IS_APPROX(m2*x21, rhs2); + VERIFY_IS_APPROX(m2*x22, rhs2); + VERIFY_IS_APPROX(x21, x22); + + // Now check with a rank deficient matrix + typedef Matrix<Scalar, RowsAtCompileTime3, ColsAtCompileTime> MatrixType3; + typedef Matrix<Scalar, RowsAtCompileTime3, 1> RhsType3; + Index rows3 = RowsAtCompileTime3==Dynamic ? internal::random<Index>(rank+1,2*cols) : Index(RowsAtCompileTime3); + Matrix<Scalar,RowsAtCompileTime3,Dynamic> C = Matrix<Scalar,RowsAtCompileTime3,Dynamic>::Random(rows3,rank); + MatrixType3 m3 = C * m2; + RhsType3 rhs3 = C * rhs2; + JacobiSVD<MatrixType3, ColPivHouseholderQRPreconditioner> svd3(m3, computationOptions); + SolutionType x3 = svd3.solve(rhs3); + if(svd3.rank()!=rank) { + std::cout << m3 << "\n\n"; + std::cout << svd3.singularValues().transpose() << "\n"; + std::cout << svd3.rank() << " == " << rank << "\n"; + std::cout << x21.norm() << " == " << x3.norm() << "\n"; + } +// VERIFY_IS_APPROX(m3*x3, rhs3); + VERIFY_IS_APPROX(m3*x21, rhs3); + VERIFY_IS_APPROX(m2*x3, rhs2); + + VERIFY_IS_APPROX(x21, x3); + } } template<typename MatrixType, int QRPreconditioner> @@ -95,27 +174,30 @@ void jacobisvd_test_all_computation_options(const MatrixType& m) if (QRPreconditioner == NoQRPreconditioner && m.rows() != m.cols()) return; JacobiSVD<MatrixType, QRPreconditioner> fullSvd(m, ComputeFullU|ComputeFullV); - - jacobisvd_check_full(m, fullSvd); - jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeFullV); - + CALL_SUBTEST(( jacobisvd_check_full(m, fullSvd) )); + CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeFullV) )); + + #if defined __INTEL_COMPILER + // remark #111: statement is unreachable + #pragma warning disable 111 + #endif if(QRPreconditioner == FullPivHouseholderQRPreconditioner) return; - jacobisvd_compare_to_full(m, ComputeFullU, fullSvd); - jacobisvd_compare_to_full(m, ComputeFullV, fullSvd); - jacobisvd_compare_to_full(m, 0, fullSvd); + CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeFullU, fullSvd) )); + CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeFullV, fullSvd) )); + CALL_SUBTEST(( jacobisvd_compare_to_full(m, 0, fullSvd) )); if (MatrixType::ColsAtCompileTime == Dynamic) { // thin U/V are only available with dynamic number of columns - jacobisvd_compare_to_full(m, ComputeFullU|ComputeThinV, fullSvd); - jacobisvd_compare_to_full(m, ComputeThinV, fullSvd); - jacobisvd_compare_to_full(m, ComputeThinU|ComputeFullV, fullSvd); - jacobisvd_compare_to_full(m, ComputeThinU , fullSvd); - jacobisvd_compare_to_full(m, ComputeThinU|ComputeThinV, fullSvd); - jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeThinV); - jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeFullV); - jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeThinV); + CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeFullU|ComputeThinV, fullSvd) )); + CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinV, fullSvd) )); + CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinU|ComputeFullV, fullSvd) )); + CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinU , fullSvd) )); + CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinU|ComputeThinV, fullSvd) )); + CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeThinV) )); + CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeFullV) )); + CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeThinV) )); // test reconstruction typedef typename MatrixType::Index Index; @@ -128,12 +210,29 @@ void jacobisvd_test_all_computation_options(const MatrixType& m) template<typename MatrixType> void jacobisvd(const MatrixType& a = MatrixType(), bool pickrandom = true) { - MatrixType m = pickrandom ? MatrixType::Random(a.rows(), a.cols()) : a; + MatrixType m = a; + if(pickrandom) + { + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef typename MatrixType::Index Index; + Index diagSize = (std::min)(a.rows(), a.cols()); + RealScalar s = std::numeric_limits<RealScalar>::max_exponent10/4; + s = internal::random<RealScalar>(1,s); + Matrix<RealScalar,Dynamic,1> d = Matrix<RealScalar,Dynamic,1>::Random(diagSize); + for(Index k=0; k<diagSize; ++k) + d(k) = d(k)*std::pow(RealScalar(10),internal::random<RealScalar>(-s,s)); + m = Matrix<Scalar,Dynamic,Dynamic>::Random(a.rows(),diagSize) * d.asDiagonal() * Matrix<Scalar,Dynamic,Dynamic>::Random(diagSize,a.cols()); + // cancel some coeffs + Index n = internal::random<Index>(0,m.size()-1); + for(Index i=0; i<n; ++i) + m(internal::random<Index>(0,m.rows()-1), internal::random<Index>(0,m.cols()-1)) = Scalar(0); + } - jacobisvd_test_all_computation_options<MatrixType, FullPivHouseholderQRPreconditioner>(m); - jacobisvd_test_all_computation_options<MatrixType, ColPivHouseholderQRPreconditioner>(m); - jacobisvd_test_all_computation_options<MatrixType, HouseholderQRPreconditioner>(m); - jacobisvd_test_all_computation_options<MatrixType, NoQRPreconditioner>(m); + CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, FullPivHouseholderQRPreconditioner>(m) )); + CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, ColPivHouseholderQRPreconditioner>(m) )); + CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, HouseholderQRPreconditioner>(m) )); + CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, NoQRPreconditioner>(m) )); } template<typename MatrixType> void jacobisvd_verify_assert(const MatrixType& m) @@ -260,7 +359,7 @@ void jacobisvd_preallocate() MatrixXf m = v.asDiagonal(); internal::set_is_malloc_allowed(false); - VERIFY_RAISES_ASSERT(VectorXf v(10);) + VERIFY_RAISES_ASSERT(VectorXf tmp(10);) JacobiSVD<MatrixXf> svd; internal::set_is_malloc_allowed(true); svd.compute(m); @@ -323,6 +422,11 @@ void test_jacobisvd() int r = internal::random<int>(1, 30), c = internal::random<int>(1, 30); + + TEST_SET_BUT_UNUSED_VARIABLE(r) + TEST_SET_BUT_UNUSED_VARIABLE(c) + + CALL_SUBTEST_10(( jacobisvd<MatrixXd>(MatrixXd(r,c)) )); CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(r,c)) )); CALL_SUBTEST_8(( jacobisvd<MatrixXcd>(MatrixXcd(r,c)) )); (void) r; diff --git a/test/linearstructure.cpp b/test/linearstructure.cpp index fd071c995..618984d5c 100644 --- a/test/linearstructure.cpp +++ b/test/linearstructure.cpp @@ -11,6 +11,7 @@ template<typename MatrixType> void linearStructure(const MatrixType& m) { + using std::abs; /* this test covers the following files: CwiseUnaryOp.h, CwiseBinaryOp.h, SelfCwiseBinaryOp.h */ @@ -27,7 +28,7 @@ template<typename MatrixType> void linearStructure(const MatrixType& m) m3(rows, cols); Scalar s1 = internal::random<Scalar>(); - while (internal::abs(s1)<1e-3) s1 = internal::random<Scalar>(); + while (abs(s1)<1e-3) s1 = internal::random<Scalar>(); Index r = internal::random<Index>(0, rows-1), c = internal::random<Index>(0, cols-1); diff --git a/test/lu.cpp b/test/lu.cpp index 6cbcb0a95..25f86755a 100644 --- a/test/lu.cpp +++ b/test/lu.cpp @@ -14,7 +14,6 @@ using namespace std; template<typename MatrixType> void lu_non_invertible() { typedef typename MatrixType::Index Index; - typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; /* this test covers the following files: LU.h @@ -100,7 +99,6 @@ template<typename MatrixType> void lu_invertible() /* this test covers the following files: LU.h */ - typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; int size = internal::random<int>(1,EIGEN_TEST_MAX_SIZE); @@ -132,8 +130,6 @@ template<typename MatrixType> void lu_partial_piv() PartialPivLU.h */ typedef typename MatrixType::Index Index; - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; Index rows = internal::random<Index>(1,4); Index cols = rows; diff --git a/test/main.h b/test/main.h index 2da327c17..14f0d2f78 100644 --- a/test/main.h +++ b/test/main.h @@ -14,6 +14,7 @@ #include <iostream> #include <fstream> #include <string> +#include <sstream> #include <vector> #include <typeinfo> #include <limits> @@ -30,6 +31,10 @@ // B0 is defined in POSIX header termios.h #define B0 FORBIDDEN_IDENTIFIER + +// shuts down ICC's remark #593: variable "XXX" was set but never used +#define TEST_SET_BUT_UNUSED_VARIABLE(X) X = X + 0; + // the following file is automatically generated by cmake #include "split_test_helper.h" @@ -37,6 +42,11 @@ #undef NDEBUG #endif +// On windows CE, NDEBUG is automatically defined <assert.h> if NDEBUG is not defined. +#ifndef DEBUG +#define DEBUG +#endif + // bounds integer values for AltiVec #ifdef __ALTIVEC__ #define EIGEN_MAKING_DOCS @@ -48,11 +58,6 @@ #define DEFAULT_REPEAT 10 -#ifdef __ICC -// disable warning #279: controlling expression is constant -#pragma warning disable 279 -#endif - namespace Eigen { static std::vector<std::string> g_test_stack; @@ -169,12 +174,17 @@ namespace Eigen #define EIGEN_INTERNAL_DEBUGGING #include <Eigen/QR> // required for createRandomPIMatrixOfRank -static void verify_impl(bool condition, const char *testname, const char *file, int line, const char *condition_as_string) +inline void verify_impl(bool condition, const char *testname, const char *file, int line, const char *condition_as_string) { if (!condition) { - std::cerr << "Test " << testname << " failed in " << file << " (" << line << ")" \ - << std::endl << " " << condition_as_string << std::endl << std::endl; \ + std::cerr << "Test " << testname << " failed in " << file << " (" << line << ")" + << std::endl << " " << condition_as_string << std::endl; + std::cerr << "Stack:\n"; + const int test_stack_size = static_cast<int>(Eigen::g_test_stack.size()); + for(int i=test_stack_size-1; i>=0; --i) + std::cerr << " - " << Eigen::g_test_stack[i] << "\n"; + std::cerr << "\n"; abort(); } } @@ -260,6 +270,7 @@ inline bool test_isApprox(const Type1& a, const Type2& b) // The idea behind this function is to compare the two scalars a and b where // the scalar ref is a hint about the expected order of magnitude of a and b. +// WARNING: the scalar a and b must be positive // Therefore, if for some reason a and b are very small compared to ref, // we won't issue a false negative. // This test could be: abs(a-b) <= eps * ref @@ -290,6 +301,10 @@ inline bool test_isUnitary(const MatrixBase<Derived>& m) return m.isUnitary(test_precision<typename internal::traits<Derived>::Scalar>()); } +// Forward declaration to avoid ICC warning +template<typename T, typename U> +bool test_is_equal(const T& actual, const U& expected); + template<typename T, typename U> bool test_is_equal(const T& actual, const U& expected) { @@ -307,6 +322,9 @@ bool test_is_equal(const T& actual, const U& expected) * A partial isometry is a matrix all of whose singular values are either 0 or 1. * This is very useful to test rank-revealing algorithms. */ +// Forward declaration to avoid ICC warning +template<typename MatrixType> +void createRandomPIMatrixOfRank(typename MatrixType::Index desired_rank, typename MatrixType::Index rows, typename MatrixType::Index cols, MatrixType& m); template<typename MatrixType> void createRandomPIMatrixOfRank(typename MatrixType::Index desired_rank, typename MatrixType::Index rows, typename MatrixType::Index cols, MatrixType& m) { @@ -345,6 +363,9 @@ void createRandomPIMatrixOfRank(typename MatrixType::Index desired_rank, typenam m = qra.householderQ() * d * qrb.householderQ(); } +// Forward declaration to avoid ICC warning +template<typename PermutationVectorType> +void randomPermutationVector(PermutationVectorType& v, typename PermutationVectorType::Index size); template<typename PermutationVectorType> void randomPermutationVector(PermutationVectorType& v, typename PermutationVectorType::Index size) { @@ -371,20 +392,22 @@ template<> struct GetDifferentType<double> { typedef float type; }; template<typename T> struct GetDifferentType<std::complex<T> > { typedef std::complex<typename GetDifferentType<T>::type> type; }; -template<typename T> std::string type_name() { return "other"; } -template<> std::string type_name<float>() { return "float"; } -template<> std::string type_name<double>() { return "double"; } -template<> std::string type_name<int>() { return "int"; } -template<> std::string type_name<std::complex<float> >() { return "complex<float>"; } +// Forward declaration to avoid ICC warning +template<typename T> std::string type_name(); +template<typename T> std::string type_name() { return "other"; } +template<> std::string type_name<float>() { return "float"; } +template<> std::string type_name<double>() { return "double"; } +template<> std::string type_name<int>() { return "int"; } +template<> std::string type_name<std::complex<float> >() { return "complex<float>"; } template<> std::string type_name<std::complex<double> >() { return "complex<double>"; } -template<> std::string type_name<std::complex<int> >() { return "complex<int>"; } +template<> std::string type_name<std::complex<int> >() { return "complex<int>"; } // forward declaration of the main test function void EIGEN_CAT(test_,EIGEN_TEST_FUNC)(); using namespace Eigen; -void set_repeat_from_string(const char *str) +inline void set_repeat_from_string(const char *str) { errno = 0; g_repeat = int(strtoul(str, 0, 10)); @@ -396,10 +419,10 @@ void set_repeat_from_string(const char *str) g_has_set_repeat = true; } -void set_seed_from_string(const char *str) +inline void set_seed_from_string(const char *str) { errno = 0; - g_seed = strtoul(str, 0, 10); + g_seed = int(strtoul(str, 0, 10)); if(errno || g_seed == 0) { std::cout << "Invalid seed value " << str << std::endl; @@ -462,11 +485,25 @@ int main(int argc, char *argv[]) if(!g_has_set_repeat) g_repeat = DEFAULT_REPEAT; std::cout << "Initializing random number generator with seed " << g_seed << std::endl; + std::stringstream ss; + ss << "Seed: " << g_seed; + g_test_stack.push_back(ss.str()); srand(g_seed); std::cout << "Repeating each test " << g_repeat << " times" << std::endl; - Eigen::g_test_stack.push_back(EI_PP_MAKE_STRING(EIGEN_TEST_FUNC)); + Eigen::g_test_stack.push_back(std::string(EI_PP_MAKE_STRING(EIGEN_TEST_FUNC))); EIGEN_CAT(test_,EIGEN_TEST_FUNC)(); return 0; } + +// These warning are disabled here such that they are still ON when parsing Eigen's header files. +#if defined __INTEL_COMPILER + // remark #383: value copied to temporary, reference to temporary used + // -> this warning is raised even for legal usage as: g_test_stack.push_back("foo"); where g_test_stack is a std::vector<std::string> + // remark #1418: external function definition with no prior declaration + // -> this warning is raised for all our test functions. Declaring them static would fix the issue. + // warning #279: controlling expression is constant + // remark #1572: floating-point equality and inequality comparisons are unreliable + #pragma warning disable 279 383 1418 1572 +#endif diff --git a/test/map.cpp b/test/mapped_matrix.cpp index fe983e802..de9dbbde3 100644 --- a/test/map.cpp +++ b/test/mapped_matrix.cpp @@ -102,9 +102,6 @@ template<typename VectorType> void map_static_methods(const VectorType& m) template<typename PlainObjectType> void check_const_correctness(const PlainObjectType&) { - typedef typename PlainObjectType::Index Index; - typedef typename PlainObjectType::Scalar Scalar; - // there's a lot that we can't test here while still having this test compile! // the only possible approach would be to run a script trying to compile stuff and checking that it fails. // CMake can help with that. @@ -117,7 +114,7 @@ template<typename PlainObjectType> void check_const_correctness(const PlainObjec VERIFY( !(Map<ConstPlainObjectType, Aligned>::Flags & LvalueBit) ); } -void test_map() +void test_mapped_matrix() { for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1( map_class_vector(Matrix<float, 1, 1>()) ); diff --git a/test/mapstride.cpp b/test/mapstride.cpp index fe35b9d23..b1dc9de2a 100644 --- a/test/mapstride.cpp +++ b/test/mapstride.cpp @@ -116,7 +116,7 @@ template<int Alignment,typename MatrixType> void map_class_matrix(const MatrixTy void test_mapstride() { for(int i = 0; i < g_repeat; i++) { - EIGEN_UNUSED int maxn = 30; + int maxn = 30; CALL_SUBTEST_1( map_class_vector<Aligned>(Matrix<float, 1, 1>()) ); CALL_SUBTEST_1( map_class_vector<Unaligned>(Matrix<float, 1, 1>()) ); CALL_SUBTEST_2( map_class_vector<Aligned>(Vector4d()) ); @@ -142,5 +142,7 @@ void test_mapstride() CALL_SUBTEST_5( map_class_matrix<Unaligned>(MatrixXi(internal::random<int>(1,maxn),internal::random<int>(1,maxn))) ); CALL_SUBTEST_6( map_class_matrix<Aligned>(MatrixXcd(internal::random<int>(1,maxn),internal::random<int>(1,maxn))) ); CALL_SUBTEST_6( map_class_matrix<Unaligned>(MatrixXcd(internal::random<int>(1,maxn),internal::random<int>(1,maxn))) ); + + TEST_SET_BUT_UNUSED_VARIABLE(maxn); } } diff --git a/test/meta.cpp b/test/meta.cpp index dc1d128d5..3302c5887 100644 --- a/test/meta.cpp +++ b/test/meta.cpp @@ -11,9 +11,6 @@ void test_meta() { - typedef float & FloatRef; - typedef const float & ConstFloatRef; - VERIFY((internal::conditional<(3<4),internal::true_type, internal::false_type>::type::value)); VERIFY(( internal::is_same<float,float>::value)); VERIFY((!internal::is_same<float,double>::value)); @@ -56,7 +53,7 @@ void test_meta() VERIFY(( internal::is_same<float,internal::remove_pointer<float* const >::type >::value)); VERIFY(internal::meta_sqrt<1>::ret == 1); - #define VERIFY_META_SQRT(X) VERIFY(internal::meta_sqrt<X>::ret == int(internal::sqrt(double(X)))) + #define VERIFY_META_SQRT(X) VERIFY(internal::meta_sqrt<X>::ret == int(std::sqrt(double(X)))) VERIFY_META_SQRT(2); VERIFY_META_SQRT(3); VERIFY_META_SQRT(4); diff --git a/test/metis_support.cpp b/test/metis_support.cpp new file mode 100644 index 000000000..932b04074 --- /dev/null +++ b/test/metis_support.cpp @@ -0,0 +1,39 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. +#include "sparse_solver.h" +#include <Eigen/SparseLU> +#include <Eigen/MetisSupport> +#include <unsupported/Eigen/SparseExtra> + +template<typename T> void test_metis_T() +{ + SparseLU<SparseMatrix<T, ColMajor>, MetisOrdering<int> > sparselu_metis; + + check_sparse_square_solving(sparselu_metis); +} + +void test_metis_support() +{ + CALL_SUBTEST_1(test_metis_T<double>()); +} diff --git a/test/miscmatrices.cpp b/test/miscmatrices.cpp index af0481cfe..ef20dc749 100644 --- a/test/miscmatrices.cpp +++ b/test/miscmatrices.cpp @@ -17,7 +17,6 @@ template<typename MatrixType> void miscMatrices(const MatrixType& m) typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; - typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType; Index rows = m.rows(); Index cols = m.cols(); diff --git a/test/nesting_ops.cpp b/test/nesting_ops.cpp index 938ebcb7a..1e8523283 100644 --- a/test/nesting_ops.cpp +++ b/test/nesting_ops.cpp @@ -12,17 +12,9 @@ template <typename MatrixType> void run_nesting_ops(const MatrixType& _m) { typename MatrixType::Nested m(_m); - typedef typename MatrixType::Scalar Scalar; -#ifdef NDEBUG - const bool is_debug = false; -#else - const bool is_debug = true; -#endif - - // Make really sure that we are in debug mode! We don't want any type of - // inlining for these tests to pass. - VERIFY(is_debug); + // Make really sure that we are in debug mode! + VERIFY_RAISES_ASSERT(eigen_assert(false)); // The only intention of these tests is to ensure that this code does // not trigger any asserts or segmentation faults... more to come. diff --git a/test/nomalloc.cpp b/test/nomalloc.cpp index d4ffcefcb..cbd02dd21 100644 --- a/test/nomalloc.cpp +++ b/test/nomalloc.cpp @@ -12,6 +12,12 @@ #ifdef __GNUC__ #define throw(X) #endif + +#ifdef __INTEL_COMPILER + // disable "warning #76: argument to macro is empty" produced by the above hack + #pragma warning disable 76 +#endif + // discard stack allocation as that too bypasses malloc #define EIGEN_STACK_ALLOCATION_LIMIT 0 // any heap allocation will raise an assert @@ -30,7 +36,6 @@ template<typename MatrixType> void nomalloc(const MatrixType& m) */ typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; - typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; Index rows = m.rows(); Index cols = m.cols(); diff --git a/test/nullary.cpp b/test/nullary.cpp index 1220e3f97..5408d88b2 100644 --- a/test/nullary.cpp +++ b/test/nullary.cpp @@ -91,6 +91,12 @@ void testVectorType(const VectorType& base) scalar.setLinSpaced(1,low,high); VERIFY_IS_APPROX( scalar, ScalarMatrix::Constant(high) ); VERIFY_IS_APPROX( ScalarMatrix::LinSpaced(1,low,high), ScalarMatrix::Constant(high) ); + + // regression test for bug 526 (linear vectorized transversal) + if (size > 1) { + m.tail(size-1).setLinSpaced(low, high); + VERIFY_IS_APPROX(m(size-1), high); + } } template<typename MatrixType> diff --git a/test/packetmath.cpp b/test/packetmath.cpp index c1464e994..2c0519c41 100644 --- a/test/packetmath.cpp +++ b/test/packetmath.cpp @@ -40,7 +40,7 @@ template<typename Scalar> bool areApprox(const Scalar* a, const Scalar* b, int s { for (int i=0; i<size; ++i) { - if (!internal::isApprox(a[i],b[i])) + if (a[i]!=b[i] && !internal::isApprox(a[i],b[i])) { std::cout << "[" << Map<const Matrix<Scalar,1,Dynamic> >(a,size) << "]" << " != " << Map<const Matrix<Scalar,1,Dynamic> >(b,size) << "\n"; return false; @@ -99,6 +99,7 @@ struct packet_helper<false,Packet> template<typename Scalar> void packetmath() { + using std::abs; typedef typename internal::packet_traits<Scalar>::type Packet; const int PacketSize = internal::packet_traits<Scalar>::size; typedef typename NumTraits<Scalar>::Real RealScalar; @@ -113,7 +114,7 @@ template<typename Scalar> void packetmath() { data1[i] = internal::random<Scalar>()/RealScalar(PacketSize); data2[i] = internal::random<Scalar>()/RealScalar(PacketSize); - refvalue = (std::max)(refvalue,internal::abs(data1[i])); + refvalue = (std::max)(refvalue,abs(data1[i])); } internal::pstore(data2, internal::pload<Packet>(data1)); @@ -144,7 +145,6 @@ template<typename Scalar> void packetmath() for (int i=0; i<PacketSize; ++i) ref[i] = data1[i+offset]; - typedef Matrix<Scalar, PacketSize, 1> Vector; VERIFY(areApprox(ref, data2, PacketSize) && "internal::palign"); } @@ -156,7 +156,7 @@ template<typename Scalar> void packetmath() CHECK_CWISE2(REF_DIV, internal::pdiv); #endif CHECK_CWISE1(internal::negate, internal::pnegate); - CHECK_CWISE1(internal::conj, internal::pconj); + CHECK_CWISE1(numext::conj, internal::pconj); for(int offset=0;offset<3;++offset) { @@ -207,6 +207,7 @@ template<typename Scalar> void packetmath() template<typename Scalar> void packetmath_real() { + using std::abs; typedef typename internal::packet_traits<Scalar>::type Packet; const int PacketSize = internal::packet_traits<Scalar>::size; @@ -217,35 +218,50 @@ template<typename Scalar> void packetmath_real() for (int i=0; i<size; ++i) { - data1[i] = internal::random<Scalar>(-1e3,1e3); - data2[i] = internal::random<Scalar>(-1e3,1e3); + data1[i] = internal::random<Scalar>(-1,1) * std::pow(Scalar(10), internal::random<Scalar>(-3,3)); + data2[i] = internal::random<Scalar>(-1,1) * std::pow(Scalar(10), internal::random<Scalar>(-3,3)); } - CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasSin, internal::sin, internal::psin); - CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasCos, internal::cos, internal::pcos); - CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasTan, internal::tan, internal::ptan); + CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasSin, std::sin, internal::psin); + CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasCos, std::cos, internal::pcos); + CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasTan, std::tan, internal::ptan); for (int i=0; i<size; ++i) { data1[i] = internal::random<Scalar>(-1,1); data2[i] = internal::random<Scalar>(-1,1); } - CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasASin, internal::asin, internal::pasin); - CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasACos, internal::acos, internal::pacos); + CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasASin, std::asin, internal::pasin); + CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasACos, std::acos, internal::pacos); for (int i=0; i<size; ++i) { data1[i] = internal::random<Scalar>(-87,88); data2[i] = internal::random<Scalar>(-87,88); } - CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasExp, internal::exp, internal::pexp); + CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasExp, std::exp, internal::pexp); for (int i=0; i<size; ++i) { - data1[i] = internal::random<Scalar>(0,1e6); - data2[i] = internal::random<Scalar>(0,1e6); + data1[i] = internal::random<Scalar>(0,1) * std::pow(Scalar(10), internal::random<Scalar>(-6,6)); + data2[i] = internal::random<Scalar>(0,1) * std::pow(Scalar(10), internal::random<Scalar>(-6,6)); } - CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasLog, internal::log, internal::plog); - CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasSqrt, internal::sqrt, internal::psqrt); + if(internal::random<float>(0,1)<0.1) + data1[internal::random<int>(0, PacketSize)] = 0; + CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasLog, std::log, internal::plog); + CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasSqrt, std::sqrt, internal::psqrt); +} + +template<typename Scalar> void packetmath_notcomplex() +{ + using std::abs; + typedef typename internal::packet_traits<Scalar>::type Packet; + const int PacketSize = internal::packet_traits<Scalar>::size; + + EIGEN_ALIGN16 Scalar data1[internal::packet_traits<Scalar>::size*4]; + EIGEN_ALIGN16 Scalar data2[internal::packet_traits<Scalar>::size*4]; + EIGEN_ALIGN16 Scalar ref[internal::packet_traits<Scalar>::size*4]; + + Array<Scalar,Dynamic,1>::Map(data1, internal::packet_traits<Scalar>::size*4).setRandom(); ref[0] = data1[0]; for (int i=0; i<PacketSize; ++i) @@ -254,7 +270,7 @@ template<typename Scalar> void packetmath_real() CHECK_CWISE2((std::min), internal::pmin); CHECK_CWISE2((std::max), internal::pmax); - CHECK_CWISE1(internal::abs, internal::pabs); + CHECK_CWISE1(abs, internal::pabs); ref[0] = data1[0]; for (int i=0; i<PacketSize; ++i) @@ -336,6 +352,10 @@ void test_packetmath() CALL_SUBTEST_1( packetmath<std::complex<float> >() ); CALL_SUBTEST_2( packetmath<std::complex<double> >() ); + CALL_SUBTEST_1( packetmath_notcomplex<float>() ); + CALL_SUBTEST_2( packetmath_notcomplex<double>() ); + CALL_SUBTEST_3( packetmath_notcomplex<int>() ); + CALL_SUBTEST_1( packetmath_real<float>() ); CALL_SUBTEST_2( packetmath_real<double>() ); diff --git a/test/pastix_support.cpp b/test/pastix_support.cpp index 0e57227f9..14da0944b 100644 --- a/test/pastix_support.cpp +++ b/test/pastix_support.cpp @@ -41,4 +41,4 @@ void test_pastix_support() CALL_SUBTEST_2(test_pastix_T<double>()); CALL_SUBTEST_3( (test_pastix_T_LU<std::complex<float> >()) ); CALL_SUBTEST_4(test_pastix_T_LU<std::complex<double> >()); -}
\ No newline at end of file +} diff --git a/test/permutationmatrices.cpp b/test/permutationmatrices.cpp index 00f666ccd..7b0dbc763 100644 --- a/test/permutationmatrices.cpp +++ b/test/permutationmatrices.cpp @@ -14,7 +14,6 @@ template<typename MatrixType> void permutationmatrices(const MatrixType& m) { typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime, Options = MatrixType::Options }; typedef PermutationMatrix<Rows> LeftPermutationType; diff --git a/test/prec_inverse_4x4.cpp b/test/prec_inverse_4x4.cpp index f7d0aff70..c4ef2d4bd 100644 --- a/test/prec_inverse_4x4.cpp +++ b/test/prec_inverse_4x4.cpp @@ -14,7 +14,6 @@ template<typename MatrixType> void inverse_permutation_4x4() { typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; Vector4i indices(0,1,2,3); for(int i = 0; i < 24; ++i) { @@ -29,6 +28,7 @@ template<typename MatrixType> void inverse_permutation_4x4() template<typename MatrixType> void inverse_general_4x4(int repeat) { + using std::abs; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; double error_sum = 0., error_max = 0.; @@ -38,7 +38,7 @@ template<typename MatrixType> void inverse_general_4x4(int repeat) RealScalar absdet; do { m = MatrixType::Random(); - absdet = internal::abs(m.determinant()); + absdet = abs(m.determinant()); } while(absdet < NumTraits<Scalar>::epsilon()); MatrixType inv = m.inverse(); double error = double( (m*inv-MatrixType::Identity()).norm() * absdet / NumTraits<Scalar>::epsilon() ); diff --git a/test/product.h b/test/product.h index 4aa9fd56d..856b234ac 100644 --- a/test/product.h +++ b/test/product.h @@ -24,7 +24,6 @@ template<typename MatrixType> void product(const MatrixType& m) */ typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits<Scalar>::NonInteger NonInteger; typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> RowVectorType; typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> ColVectorType; typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> RowSquareMatrixType; diff --git a/test/product_extra.cpp b/test/product_extra.cpp index 9a6bf0792..744a1ef7f 100644 --- a/test/product_extra.cpp +++ b/test/product_extra.cpp @@ -13,7 +13,6 @@ template<typename MatrixType> void product_extra(const MatrixType& m) { typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits<Scalar>::NonInteger NonInteger; typedef Matrix<Scalar, 1, Dynamic> RowVectorType; typedef Matrix<Scalar, Dynamic, 1> ColVectorType; typedef Matrix<Scalar, Dynamic, Dynamic, @@ -43,7 +42,7 @@ template<typename MatrixType> void product_extra(const MatrixType& m) VERIFY_IS_APPROX(m3.noalias() = m1.adjoint() * square.adjoint(), m1.adjoint().eval() * square.adjoint().eval()); VERIFY_IS_APPROX(m3.noalias() = m1.adjoint() * m2, m1.adjoint().eval() * m2); VERIFY_IS_APPROX(m3.noalias() = (s1 * m1.adjoint()) * m2, (s1 * m1.adjoint()).eval() * m2); - VERIFY_IS_APPROX(m3.noalias() = ((s1 * m1).adjoint()) * m2, (internal::conj(s1) * m1.adjoint()).eval() * m2); + VERIFY_IS_APPROX(m3.noalias() = ((s1 * m1).adjoint()) * m2, (numext::conj(s1) * m1.adjoint()).eval() * m2); VERIFY_IS_APPROX(m3.noalias() = (- m1.adjoint() * s1) * (s3 * m2), (- m1.adjoint() * s1).eval() * (s3 * m2).eval()); VERIFY_IS_APPROX(m3.noalias() = (s2 * m1.adjoint() * s1) * m2, (s2 * m1.adjoint() * s1).eval() * m2); VERIFY_IS_APPROX(m3.noalias() = (-m1*s2) * s1*m2.adjoint(), (-m1*s2).eval() * (s1*m2.adjoint()).eval()); @@ -135,6 +134,35 @@ void zero_sized_objects() a*b; } +void unaligned_objects() +{ + // Regression test for the bug reported here: + // http://forum.kde.org/viewtopic.php?f=74&t=107541 + // Recall the matrix*vector kernel avoid unaligned loads by loading two packets and then reassemble then. + // There was a mistake in the computation of the valid range for fully unaligned objects: in some rare cases, + // memory was read outside the allocated matrix memory. Though the values were not used, this might raise segfault. + for(int m=450;m<460;++m) + { + for(int n=8;n<12;++n) + { + MatrixXf M(m, n); + VectorXf v1(n), r1(500); + RowVectorXf v2(m), r2(16); + + M.setRandom(); + v1.setRandom(); + v2.setRandom(); + for(int o=0; o<4; ++o) + { + r1.segment(o,m).noalias() = M * v1; + VERIFY_IS_APPROX(r1.segment(o,m), M * MatrixXf(v1)); + r2.segment(o,n).noalias() = v2 * M; + VERIFY_IS_APPROX(r2.segment(o,n), MatrixXf(v2) * M); + } + } + } +} + void test_product_extra() { for(int i = 0; i < g_repeat; i++) { @@ -143,6 +171,7 @@ void test_product_extra() CALL_SUBTEST_2( mat_mat_scalar_scalar_product() ); CALL_SUBTEST_3( product_extra(MatrixXcf(internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2), internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2))) ); CALL_SUBTEST_4( product_extra(MatrixXcd(internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2), internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2))) ); - CALL_SUBTEST_5( zero_sized_objects() ); } + CALL_SUBTEST_5( zero_sized_objects() ); + CALL_SUBTEST_6( unaligned_objects() ); } diff --git a/test/product_mmtr.cpp b/test/product_mmtr.cpp index 879cfe16a..7d6746800 100644 --- a/test/product_mmtr.cpp +++ b/test/product_mmtr.cpp @@ -19,8 +19,6 @@ template<typename Scalar> void mmtr(int size) { - typedef typename NumTraits<Scalar>::Real RealScalar; - typedef Matrix<Scalar,Dynamic,Dynamic,ColMajor> MatrixColMaj; typedef Matrix<Scalar,Dynamic,Dynamic,RowMajor> MatrixRowMaj; diff --git a/test/product_notemporary.cpp b/test/product_notemporary.cpp index cf9dbdd03..258d238e2 100644 --- a/test/product_notemporary.cpp +++ b/test/product_notemporary.cpp @@ -9,7 +9,7 @@ static int nb_temporaries; -void on_temporary_creation(int size) { +inline void on_temporary_creation(int size) { // here's a great place to set a breakpoint when debugging failures in this test! if(size!=0) nb_temporaries++; } @@ -77,6 +77,9 @@ template<typename MatrixType> void product_notemporary(const MatrixType& m) VERIFY_EVALUATION_COUNT( m3.noalias() -= (s1 * m1).template triangularView<Lower>() * m2, 0); VERIFY_EVALUATION_COUNT( rm3.noalias() = (s1 * m1.adjoint()).template triangularView<Upper>() * (m2+m2), 1); VERIFY_EVALUATION_COUNT( rm3.noalias() = (s1 * m1.adjoint()).template triangularView<UnitUpper>() * m2.adjoint(), 0); + + VERIFY_EVALUATION_COUNT( m3.template triangularView<Upper>() = (m1 * m2.adjoint()), 0); + VERIFY_EVALUATION_COUNT( m3.template triangularView<Upper>() -= (m1 * m2.adjoint()), 0); // NOTE this is because the blas_traits require innerstride==1 to avoid a temporary, but that doesn't seem to be actually needed for the triangular products VERIFY_EVALUATION_COUNT( rm3.col(c0).noalias() = (s1 * m1.adjoint()).template triangularView<UnitUpper>() * (s2*m2.row(c0)).adjoint(), 1); diff --git a/test/product_selfadjoint.cpp b/test/product_selfadjoint.cpp index 95693b155..374e2393b 100644 --- a/test/product_selfadjoint.cpp +++ b/test/product_selfadjoint.cpp @@ -13,7 +13,6 @@ template<typename MatrixType> void product_selfadjoint(const MatrixType& m) { typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits<Scalar>::Real RealScalar; typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; typedef Matrix<Scalar, 1, MatrixType::RowsAtCompileTime> RowVectorType; @@ -45,11 +44,11 @@ template<typename MatrixType> void product_selfadjoint(const MatrixType& m) m2 = m1.template triangularView<Upper>(); m2.template selfadjointView<Upper>().rankUpdate(-v1,s2*v2,s3); - VERIFY_IS_APPROX(m2, (m1 + (s3*(-v1)*(s2*v2).adjoint()+internal::conj(s3)*(s2*v2)*(-v1).adjoint())).template triangularView<Upper>().toDenseMatrix()); + VERIFY_IS_APPROX(m2, (m1 + (s3*(-v1)*(s2*v2).adjoint()+numext::conj(s3)*(s2*v2)*(-v1).adjoint())).template triangularView<Upper>().toDenseMatrix()); m2 = m1.template triangularView<Upper>(); m2.template selfadjointView<Upper>().rankUpdate(-s2*r1.adjoint(),r2.adjoint()*s3,s1); - VERIFY_IS_APPROX(m2, (m1 + s1*(-s2*r1.adjoint())*(r2.adjoint()*s3).adjoint() + internal::conj(s1)*(r2.adjoint()*s3) * (-s2*r1.adjoint()).adjoint()).template triangularView<Upper>().toDenseMatrix()); + VERIFY_IS_APPROX(m2, (m1 + s1*(-s2*r1.adjoint())*(r2.adjoint()*s3).adjoint() + numext::conj(s1)*(r2.adjoint()*s3) * (-s2*r1.adjoint()).adjoint()).template triangularView<Upper>().toDenseMatrix()); if (rows>1) { @@ -63,7 +62,7 @@ template<typename MatrixType> void product_selfadjoint(const MatrixType& m) void test_product_selfadjoint() { - int s; + int s = 0; for(int i = 0; i < g_repeat ; i++) { CALL_SUBTEST_1( product_selfadjoint(Matrix<float, 1, 1>()) ); CALL_SUBTEST_2( product_selfadjoint(Matrix<float, 2, 2>()) ); @@ -77,5 +76,5 @@ void test_product_selfadjoint() s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE); CALL_SUBTEST_7( product_selfadjoint(Matrix<float,Dynamic,Dynamic,RowMajor>(s,s)) ); } - EIGEN_UNUSED_VARIABLE(s) + TEST_SET_BUT_UNUSED_VARIABLE(s) } diff --git a/test/product_symm.cpp b/test/product_symm.cpp index 2f7a0d231..74d7329b1 100644 --- a/test/product_symm.cpp +++ b/test/product_symm.cpp @@ -11,8 +11,6 @@ template<typename Scalar, int Size, int OtherSize> void symm(int size = Size, int othersize = OtherSize) { - typedef typename NumTraits<Scalar>::Real RealScalar; - typedef Matrix<Scalar, Size, Size> MatrixType; typedef Matrix<Scalar, Size, OtherSize> Rhs1; typedef Matrix<Scalar, OtherSize, Size> Rhs2; diff --git a/test/product_syrk.cpp b/test/product_syrk.cpp index 5855c2181..73c95000c 100644 --- a/test/product_syrk.cpp +++ b/test/product_syrk.cpp @@ -13,7 +13,7 @@ template<typename MatrixType> void syrk(const MatrixType& m) { typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime, RowMajor> RMatrixType; typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, Dynamic> Rhs1; typedef Matrix<Scalar, Dynamic, MatrixType::RowsAtCompileTime> Rhs2; typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, Dynamic,RowMajor> Rhs3; @@ -22,10 +22,12 @@ template<typename MatrixType> void syrk(const MatrixType& m) Index cols = m.cols(); MatrixType m1 = MatrixType::Random(rows, cols), - m2 = MatrixType::Random(rows, cols); + m2 = MatrixType::Random(rows, cols), + m3 = MatrixType::Random(rows, cols); + RMatrixType rm2 = MatrixType::Random(rows, cols); - Rhs1 rhs1 = Rhs1::Random(internal::random<int>(1,320), cols); - Rhs2 rhs2 = Rhs2::Random(rows, internal::random<int>(1,320)); + Rhs1 rhs1 = Rhs1::Random(internal::random<int>(1,320), cols); Rhs1 rhs11 = Rhs1::Random(rhs1.rows(), cols); + Rhs2 rhs2 = Rhs2::Random(rows, internal::random<int>(1,320)); Rhs2 rhs22 = Rhs2::Random(rows, rhs2.cols()); Rhs3 rhs3 = Rhs3::Random(internal::random<int>(1,320), rows); Scalar s1 = internal::random<Scalar>(); @@ -35,19 +37,34 @@ template<typename MatrixType> void syrk(const MatrixType& m) m2.setZero(); VERIFY_IS_APPROX((m2.template selfadjointView<Lower>().rankUpdate(rhs2,s1)._expression()), ((s1 * rhs2 * rhs2.adjoint()).eval().template triangularView<Lower>().toDenseMatrix())); + m2.setZero(); + VERIFY_IS_APPROX(((m2.template triangularView<Lower>() += s1 * rhs2 * rhs22.adjoint()).nestedExpression()), + ((s1 * rhs2 * rhs22.adjoint()).eval().template triangularView<Lower>().toDenseMatrix())); + m2.setZero(); VERIFY_IS_APPROX(m2.template selfadjointView<Upper>().rankUpdate(rhs2,s1)._expression(), (s1 * rhs2 * rhs2.adjoint()).eval().template triangularView<Upper>().toDenseMatrix()); + m2.setZero(); + VERIFY_IS_APPROX((m2.template triangularView<Upper>() += s1 * rhs22 * rhs2.adjoint()).nestedExpression(), + (s1 * rhs22 * rhs2.adjoint()).eval().template triangularView<Upper>().toDenseMatrix()); + m2.setZero(); VERIFY_IS_APPROX(m2.template selfadjointView<Lower>().rankUpdate(rhs1.adjoint(),s1)._expression(), (s1 * rhs1.adjoint() * rhs1).eval().template triangularView<Lower>().toDenseMatrix()); - + m2.setZero(); + VERIFY_IS_APPROX((m2.template triangularView<Lower>() += s1 * rhs11.adjoint() * rhs1).nestedExpression(), + (s1 * rhs11.adjoint() * rhs1).eval().template triangularView<Lower>().toDenseMatrix()); + + m2.setZero(); VERIFY_IS_APPROX(m2.template selfadjointView<Upper>().rankUpdate(rhs1.adjoint(),s1)._expression(), (s1 * rhs1.adjoint() * rhs1).eval().template triangularView<Upper>().toDenseMatrix()); + VERIFY_IS_APPROX((m2.template triangularView<Upper>() = s1 * rhs1.adjoint() * rhs11).nestedExpression(), + (s1 * rhs1.adjoint() * rhs11).eval().template triangularView<Upper>().toDenseMatrix()); + m2.setZero(); VERIFY_IS_APPROX(m2.template selfadjointView<Lower>().rankUpdate(rhs3.adjoint(),s1)._expression(), (s1 * rhs3.adjoint() * rhs3).eval().template triangularView<Lower>().toDenseMatrix()); @@ -63,6 +80,15 @@ template<typename MatrixType> void syrk(const MatrixType& m) m2.setZero(); VERIFY_IS_APPROX((m2.template selfadjointView<Upper>().rankUpdate(m1.col(c),s1)._expression()), ((s1 * m1.col(c) * m1.col(c).adjoint()).eval().template triangularView<Upper>().toDenseMatrix())); + rm2.setZero(); + VERIFY_IS_APPROX((rm2.template selfadjointView<Upper>().rankUpdate(m1.col(c),s1)._expression()), + ((s1 * m1.col(c) * m1.col(c).adjoint()).eval().template triangularView<Upper>().toDenseMatrix())); + m2.setZero(); + VERIFY_IS_APPROX((m2.template triangularView<Upper>() += s1 * m3.col(c) * m1.col(c).adjoint()).nestedExpression(), + ((s1 * m3.col(c) * m1.col(c).adjoint()).eval().template triangularView<Upper>().toDenseMatrix())); + rm2.setZero(); + VERIFY_IS_APPROX((rm2.template triangularView<Upper>() += s1 * m1.col(c) * m3.col(c).adjoint()).nestedExpression(), + ((s1 * m1.col(c) * m3.col(c).adjoint()).eval().template triangularView<Upper>().toDenseMatrix())); m2.setZero(); VERIFY_IS_APPROX((m2.template selfadjointView<Lower>().rankUpdate(m1.col(c).conjugate(),s1)._expression()), @@ -72,9 +98,20 @@ template<typename MatrixType> void syrk(const MatrixType& m) VERIFY_IS_APPROX((m2.template selfadjointView<Upper>().rankUpdate(m1.col(c).conjugate(),s1)._expression()), ((s1 * m1.col(c).conjugate() * m1.col(c).conjugate().adjoint()).eval().template triangularView<Upper>().toDenseMatrix())); + m2.setZero(); VERIFY_IS_APPROX((m2.template selfadjointView<Lower>().rankUpdate(m1.row(c),s1)._expression()), ((s1 * m1.row(c).transpose() * m1.row(c).transpose().adjoint()).eval().template triangularView<Lower>().toDenseMatrix())); + rm2.setZero(); + VERIFY_IS_APPROX((rm2.template selfadjointView<Lower>().rankUpdate(m1.row(c),s1)._expression()), + ((s1 * m1.row(c).transpose() * m1.row(c).transpose().adjoint()).eval().template triangularView<Lower>().toDenseMatrix())); + m2.setZero(); + VERIFY_IS_APPROX((m2.template triangularView<Lower>() += s1 * m3.row(c).transpose() * m1.row(c).transpose().adjoint()).nestedExpression(), + ((s1 * m3.row(c).transpose() * m1.row(c).transpose().adjoint()).eval().template triangularView<Lower>().toDenseMatrix())); + rm2.setZero(); + VERIFY_IS_APPROX((rm2.template triangularView<Lower>() += s1 * m3.row(c).transpose() * m1.row(c).transpose().adjoint()).nestedExpression(), + ((s1 * m3.row(c).transpose() * m1.row(c).transpose().adjoint()).eval().template triangularView<Lower>().toDenseMatrix())); + m2.setZero(); VERIFY_IS_APPROX((m2.template selfadjointView<Upper>().rankUpdate(m1.row(c).adjoint(),s1)._expression()), diff --git a/test/product_trmm.cpp b/test/product_trmm.cpp index 64244c18f..d715b9a36 100644 --- a/test/product_trmm.cpp +++ b/test/product_trmm.cpp @@ -14,8 +14,6 @@ void trmm(int rows=internal::random<int>(1,EIGEN_TEST_MAX_SIZE), int cols=internal::random<int>(1,EIGEN_TEST_MAX_SIZE), int otherCols = OtherCols==Dynamic?internal::random<int>(1,EIGEN_TEST_MAX_SIZE):OtherCols) { - typedef typename NumTraits<Scalar>::Real RealScalar; - typedef Matrix<Scalar,Dynamic,Dynamic,TriOrder> TriMatrix; typedef Matrix<Scalar,Dynamic,OtherCols,OtherCols==1?ColMajor:OtherOrder> OnTheRight; typedef Matrix<Scalar,OtherCols,Dynamic,OtherCols==1?RowMajor:OtherOrder> OnTheLeft; @@ -53,10 +51,11 @@ void trmm(int rows=internal::random<int>(1,EIGEN_TEST_MAX_SIZE), ge_xs_save = ge_xs; VERIFY_IS_APPROX( (ge_xs_save + s1*triTr.conjugate() * (s2*ge_left.adjoint())).eval(), ge_xs.noalias() += (s1*mat.adjoint()).template triangularView<Mode>() * (s2*ge_left.adjoint()) ); + ge_sx.setRandom(); ge_sx_save = ge_sx; VERIFY_IS_APPROX( ge_sx_save - (ge_right.adjoint() * (-s1 * triTr).conjugate()).eval(), ge_sx.noalias() -= (ge_right.adjoint() * (-s1 * mat).adjoint().template triangularView<Mode>()).eval()); - VERIFY_IS_APPROX( ge_xs = (s1*mat).adjoint().template triangularView<Mode>() * ge_left.adjoint(), internal::conj(s1) * triTr.conjugate() * ge_left.adjoint()); + VERIFY_IS_APPROX( ge_xs = (s1*mat).adjoint().template triangularView<Mode>() * ge_left.adjoint(), numext::conj(s1) * triTr.conjugate() * ge_left.adjoint()); // TODO check with sub-matrix expressions ? } diff --git a/test/product_trmv.cpp b/test/product_trmv.cpp index 435018e8e..4c3c435c2 100644 --- a/test/product_trmv.cpp +++ b/test/product_trmv.cpp @@ -73,7 +73,7 @@ template<typename MatrixType> void trmv(const MatrixType& m) void test_product_trmv() { - int s; + int s = 0; for(int i = 0; i < g_repeat ; i++) { CALL_SUBTEST_1( trmv(Matrix<float, 1, 1>()) ); CALL_SUBTEST_2( trmv(Matrix<float, 2, 2>()) ); @@ -85,5 +85,5 @@ void test_product_trmv() s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE); CALL_SUBTEST_6( trmv(Matrix<float,Dynamic,Dynamic,RowMajor>(s, s)) ); } - EIGEN_UNUSED_VARIABLE(s); + TEST_SET_BUT_UNUSED_VARIABLE(s); } diff --git a/test/qr.cpp b/test/qr.cpp index 37fb7aa4d..a79e0dd34 100644 --- a/test/qr.cpp +++ b/test/qr.cpp @@ -19,7 +19,6 @@ template<typename MatrixType> void qr(const MatrixType& m) typedef typename MatrixType::Scalar Scalar; typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType; - typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType; MatrixType a = MatrixType::Random(rows,cols); HouseholderQR<MatrixType> qrOfA(a); @@ -53,6 +52,8 @@ template<typename MatrixType, int Cols2> void qr_fixedsize() template<typename MatrixType> void qr_invertible() { + using std::log; + using std::abs; typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; typedef typename MatrixType::Scalar Scalar; @@ -76,12 +77,12 @@ template<typename MatrixType> void qr_invertible() // now construct a matrix with prescribed determinant m1.setZero(); for(int i = 0; i < size; i++) m1(i,i) = internal::random<Scalar>(); - RealScalar absdet = internal::abs(m1.diagonal().prod()); + RealScalar absdet = abs(m1.diagonal().prod()); m3 = qr.householderQ(); // get a unitary m1 = m3 * m1 * m3; qr.compute(m1); VERIFY_IS_APPROX(absdet, qr.absDeterminant()); - VERIFY_IS_APPROX(internal::log(absdet), qr.logAbsDeterminant()); + VERIFY_IS_APPROX(log(absdet), qr.logAbsDeterminant()); } template<typename MatrixType> void qr_verify_assert() diff --git a/test/qr_colpivoting.cpp b/test/qr_colpivoting.cpp index dd0812819..eb3feac01 100644 --- a/test/qr_colpivoting.cpp +++ b/test/qr_colpivoting.cpp @@ -19,9 +19,7 @@ template<typename MatrixType> void qr() Index rank = internal::random<Index>(1, (std::min)(rows, cols)-1); typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType; - typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType; MatrixType m1; createRandomPIMatrixOfRank(rank,rows,cols,m1); ColPivHouseholderQR<MatrixType> qr(m1); @@ -72,6 +70,8 @@ template<typename MatrixType, int Cols2> void qr_fixedsize() template<typename MatrixType> void qr_invertible() { + using std::log; + using std::abs; typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; typedef typename MatrixType::Scalar Scalar; @@ -95,12 +95,12 @@ template<typename MatrixType> void qr_invertible() // now construct a matrix with prescribed determinant m1.setZero(); for(int i = 0; i < size; i++) m1(i,i) = internal::random<Scalar>(); - RealScalar absdet = internal::abs(m1.diagonal().prod()); + RealScalar absdet = abs(m1.diagonal().prod()); m3 = qr.householderQ(); // get a unitary m1 = m3 * m1 * m3; qr.compute(m1); VERIFY_IS_APPROX(absdet, qr.absDeterminant()); - VERIFY_IS_APPROX(internal::log(absdet), qr.logAbsDeterminant()); + VERIFY_IS_APPROX(log(absdet), qr.logAbsDeterminant()); } template<typename MatrixType> void qr_verify_assert() diff --git a/test/qr_fullpivoting.cpp b/test/qr_fullpivoting.cpp index e5c9790c8..511f2473f 100644 --- a/test/qr_fullpivoting.cpp +++ b/test/qr_fullpivoting.cpp @@ -20,7 +20,6 @@ template<typename MatrixType> void qr() typedef typename MatrixType::Scalar Scalar; typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType; - typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType; MatrixType m1; createRandomPIMatrixOfRank(rank,rows,cols,m1); FullPivHouseholderQR<MatrixType> qr(m1); @@ -51,6 +50,8 @@ template<typename MatrixType> void qr() template<typename MatrixType> void qr_invertible() { + using std::log; + using std::abs; typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; typedef typename MatrixType::Scalar Scalar; @@ -78,12 +79,12 @@ template<typename MatrixType> void qr_invertible() // now construct a matrix with prescribed determinant m1.setZero(); for(int i = 0; i < size; i++) m1(i,i) = internal::random<Scalar>(); - RealScalar absdet = internal::abs(m1.diagonal().prod()); + RealScalar absdet = abs(m1.diagonal().prod()); m3 = qr.matrixQ(); // get a unitary m1 = m3 * m1 * m3; qr.compute(m1); VERIFY_IS_APPROX(absdet, qr.absDeterminant()); - VERIFY_IS_APPROX(internal::log(absdet), qr.logAbsDeterminant()); + VERIFY_IS_APPROX(log(absdet), qr.logAbsDeterminant()); } template<typename MatrixType> void qr_verify_assert() @@ -129,4 +130,8 @@ void test_qr_fullpivoting() // Test problem size constructors CALL_SUBTEST_7(FullPivHouseholderQR<MatrixXf>(10, 20)); + CALL_SUBTEST_7((FullPivHouseholderQR<Matrix<float,10,20> >(10,20))); + CALL_SUBTEST_7((FullPivHouseholderQR<Matrix<float,10,20> >(Matrix<float,10,20>::Random()))); + CALL_SUBTEST_7((FullPivHouseholderQR<Matrix<float,20,10> >(20,10))); + CALL_SUBTEST_7((FullPivHouseholderQR<Matrix<float,20,10> >(Matrix<float,20,10>::Random()))); } diff --git a/test/real_qz.cpp b/test/real_qz.cpp new file mode 100644 index 000000000..7d743a734 --- /dev/null +++ b/test/real_qz.cpp @@ -0,0 +1,65 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Alexey Korepanov <kaikaikai@yandex.ru> +// +// 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/. + +#include "main.h" +#include <limits> +#include <Eigen/Eigenvalues> + +template<typename MatrixType> void real_qz(const MatrixType& m) +{ + /* this test covers the following files: + RealQZ.h + */ + using std::abs; + typedef typename MatrixType::Index Index; + typedef typename MatrixType::Scalar Scalar; + + Index dim = m.cols(); + + MatrixType A = MatrixType::Random(dim,dim), + B = MatrixType::Random(dim,dim); + + RealQZ<MatrixType> qz(A,B); + + VERIFY_IS_EQUAL(qz.info(), Success); + // check for zeros + bool all_zeros = true; + for (Index i=0; i<A.cols(); i++) + for (Index j=0; j<i; j++) { + if (abs(qz.matrixT()(i,j))!=Scalar(0.0)) + all_zeros = false; + if (j<i-1 && abs(qz.matrixS()(i,j))!=Scalar(0.0)) + all_zeros = false; + if (j==i-1 && j>0 && abs(qz.matrixS()(i,j))!=Scalar(0.0) && abs(qz.matrixS()(i-1,j-1))!=Scalar(0.0)) + all_zeros = false; + } + VERIFY_IS_EQUAL(all_zeros, true); + VERIFY_IS_APPROX(qz.matrixQ()*qz.matrixS()*qz.matrixZ(), A); + VERIFY_IS_APPROX(qz.matrixQ()*qz.matrixT()*qz.matrixZ(), B); + VERIFY_IS_APPROX(qz.matrixQ()*qz.matrixQ().adjoint(), MatrixType::Identity(dim,dim)); + VERIFY_IS_APPROX(qz.matrixZ()*qz.matrixZ().adjoint(), MatrixType::Identity(dim,dim)); +} + +void test_real_qz() +{ + int s = 0; + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( real_qz(Matrix4f()) ); + s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4); + CALL_SUBTEST_2( real_qz(MatrixXd(s,s)) ); + + // some trivial but implementation-wise tricky cases + CALL_SUBTEST_2( real_qz(MatrixXd(1,1)) ); + CALL_SUBTEST_2( real_qz(MatrixXd(2,2)) ); + CALL_SUBTEST_3( real_qz(Matrix<double,1,1>()) ); + CALL_SUBTEST_4( real_qz(Matrix2d()) ); + } + + TEST_SET_BUT_UNUSED_VARIABLE(s) +} diff --git a/test/redux.cpp b/test/redux.cpp index e07d4b1e4..0d176e500 100644 --- a/test/redux.cpp +++ b/test/redux.cpp @@ -22,26 +22,26 @@ template<typename MatrixType> void matrixRedux(const MatrixType& m) // The entries of m1 are uniformly distributed in [0,1], so m1.prod() is very small. This may lead to test // failures if we underflow into denormals. Thus, we scale so that entires are close to 1. - MatrixType m1_for_prod = MatrixType::Ones(rows, cols) + Scalar(0.2) * m1; + MatrixType m1_for_prod = MatrixType::Ones(rows, cols) + RealScalar(0.2) * m1; VERIFY_IS_MUCH_SMALLER_THAN(MatrixType::Zero(rows, cols).sum(), Scalar(1)); VERIFY_IS_APPROX(MatrixType::Ones(rows, cols).sum(), Scalar(float(rows*cols))); // the float() here to shut up excessive MSVC warning about int->complex conversion being lossy - Scalar s(0), p(1), minc(internal::real(m1.coeff(0))), maxc(internal::real(m1.coeff(0))); + Scalar s(0), p(1), minc(numext::real(m1.coeff(0))), maxc(numext::real(m1.coeff(0))); for(int j = 0; j < cols; j++) for(int i = 0; i < rows; i++) { s += m1(i,j); p *= m1_for_prod(i,j); - minc = (std::min)(internal::real(minc), internal::real(m1(i,j))); - maxc = (std::max)(internal::real(maxc), internal::real(m1(i,j))); + minc = (std::min)(numext::real(minc), numext::real(m1(i,j))); + maxc = (std::max)(numext::real(maxc), numext::real(m1(i,j))); } const Scalar mean = s/Scalar(RealScalar(rows*cols)); VERIFY_IS_APPROX(m1.sum(), s); VERIFY_IS_APPROX(m1.mean(), mean); VERIFY_IS_APPROX(m1_for_prod.prod(), p); - VERIFY_IS_APPROX(m1.real().minCoeff(), internal::real(minc)); - VERIFY_IS_APPROX(m1.real().maxCoeff(), internal::real(maxc)); + VERIFY_IS_APPROX(m1.real().minCoeff(), numext::real(minc)); + VERIFY_IS_APPROX(m1.real().maxCoeff(), numext::real(maxc)); // test slice vectorization assuming assign is ok Index r0 = internal::random<Index>(0,rows-1); @@ -61,6 +61,7 @@ template<typename MatrixType> void matrixRedux(const MatrixType& m) template<typename VectorType> void vectorRedux(const VectorType& w) { + using std::abs; typedef typename VectorType::Index Index; typedef typename VectorType::Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; @@ -72,15 +73,15 @@ template<typename VectorType> void vectorRedux(const VectorType& w) for(int i = 1; i < size; i++) { Scalar s(0), p(1); - RealScalar minc(internal::real(v.coeff(0))), maxc(internal::real(v.coeff(0))); + RealScalar minc(numext::real(v.coeff(0))), maxc(numext::real(v.coeff(0))); for(int j = 0; j < i; j++) { s += v[j]; p *= v_for_prod[j]; - minc = (std::min)(minc, internal::real(v[j])); - maxc = (std::max)(maxc, internal::real(v[j])); + minc = (std::min)(minc, numext::real(v[j])); + maxc = (std::max)(maxc, numext::real(v[j])); } - VERIFY_IS_MUCH_SMALLER_THAN(internal::abs(s - v.head(i).sum()), Scalar(1)); + VERIFY_IS_MUCH_SMALLER_THAN(abs(s - v.head(i).sum()), Scalar(1)); VERIFY_IS_APPROX(p, v_for_prod.head(i).prod()); VERIFY_IS_APPROX(minc, v.real().head(i).minCoeff()); VERIFY_IS_APPROX(maxc, v.real().head(i).maxCoeff()); @@ -89,15 +90,15 @@ template<typename VectorType> void vectorRedux(const VectorType& w) for(int i = 0; i < size-1; i++) { Scalar s(0), p(1); - RealScalar minc(internal::real(v.coeff(i))), maxc(internal::real(v.coeff(i))); + RealScalar minc(numext::real(v.coeff(i))), maxc(numext::real(v.coeff(i))); for(int j = i; j < size; j++) { s += v[j]; p *= v_for_prod[j]; - minc = (std::min)(minc, internal::real(v[j])); - maxc = (std::max)(maxc, internal::real(v[j])); + minc = (std::min)(minc, numext::real(v[j])); + maxc = (std::max)(maxc, numext::real(v[j])); } - VERIFY_IS_MUCH_SMALLER_THAN(internal::abs(s - v.tail(size-i).sum()), Scalar(1)); + VERIFY_IS_MUCH_SMALLER_THAN(abs(s - v.tail(size-i).sum()), Scalar(1)); VERIFY_IS_APPROX(p, v_for_prod.tail(size-i).prod()); VERIFY_IS_APPROX(minc, v.real().tail(size-i).minCoeff()); VERIFY_IS_APPROX(maxc, v.real().tail(size-i).maxCoeff()); @@ -106,15 +107,15 @@ template<typename VectorType> void vectorRedux(const VectorType& w) for(int i = 0; i < size/2; i++) { Scalar s(0), p(1); - RealScalar minc(internal::real(v.coeff(i))), maxc(internal::real(v.coeff(i))); + RealScalar minc(numext::real(v.coeff(i))), maxc(numext::real(v.coeff(i))); for(int j = i; j < size-i; j++) { s += v[j]; p *= v_for_prod[j]; - minc = (std::min)(minc, internal::real(v[j])); - maxc = (std::max)(maxc, internal::real(v[j])); + minc = (std::min)(minc, numext::real(v[j])); + maxc = (std::max)(maxc, numext::real(v[j])); } - VERIFY_IS_MUCH_SMALLER_THAN(internal::abs(s - v.segment(i, size-2*i).sum()), Scalar(1)); + VERIFY_IS_MUCH_SMALLER_THAN(abs(s - v.segment(i, size-2*i).sum()), Scalar(1)); VERIFY_IS_APPROX(p, v_for_prod.segment(i, size-2*i).prod()); VERIFY_IS_APPROX(minc, v.real().segment(i, size-2*i).minCoeff()); VERIFY_IS_APPROX(maxc, v.real().segment(i, size-2*i).maxCoeff()); @@ -132,7 +133,7 @@ void test_redux() { // the max size cannot be too large, otherwise reduxion operations obviously generate large errors. int maxsize = (std::min)(100,EIGEN_TEST_MAX_SIZE); - EIGEN_UNUSED_VARIABLE(maxsize); + TEST_SET_BUT_UNUSED_VARIABLE(maxsize); for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1( matrixRedux(Matrix<float, 1, 1>()) ); CALL_SUBTEST_1( matrixRedux(Array<float, 1, 1>()) ); diff --git a/test/ref.cpp b/test/ref.cpp new file mode 100644 index 000000000..19e81549c --- /dev/null +++ b/test/ref.cpp @@ -0,0 +1,252 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 20013 Gael Guennebaud <gael.guennebaud@inria.fr> +// +// 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/. + +// This unit test cannot be easily written to work with EIGEN_DEFAULT_TO_ROW_MAJOR +#ifdef EIGEN_DEFAULT_TO_ROW_MAJOR +#undef EIGEN_DEFAULT_TO_ROW_MAJOR +#endif + +static int nb_temporaries; + +inline void on_temporary_creation(int) { + // here's a great place to set a breakpoint when debugging failures in this test! + nb_temporaries++; +} + + +#define EIGEN_DENSE_STORAGE_CTOR_PLUGIN { on_temporary_creation(size); } + +#include "main.h" + +#define VERIFY_EVALUATION_COUNT(XPR,N) {\ + nb_temporaries = 0; \ + XPR; \ + if(nb_temporaries!=N) std::cerr << "nb_temporaries == " << nb_temporaries << "\n"; \ + VERIFY( (#XPR) && nb_temporaries==N ); \ + } + + +// test Ref.h + +template<typename MatrixType> void ref_matrix(const MatrixType& m) +{ + typedef typename MatrixType::Index Index; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef Matrix<Scalar,Dynamic,Dynamic,MatrixType::Options> DynMatrixType; + typedef Matrix<RealScalar,Dynamic,Dynamic,MatrixType::Options> RealDynMatrixType; + + typedef Ref<MatrixType> RefMat; + typedef Ref<DynMatrixType> RefDynMat; + typedef Ref<const DynMatrixType> ConstRefDynMat; + typedef Ref<RealDynMatrixType , 0, Stride<Dynamic,Dynamic> > RefRealMatWithStride; + + Index rows = m.rows(), cols = m.cols(); + + MatrixType m1 = MatrixType::Random(rows, cols), + m2 = m1; + + Index i = internal::random<Index>(0,rows-1); + Index j = internal::random<Index>(0,cols-1); + Index brows = internal::random<Index>(1,rows-i); + Index bcols = internal::random<Index>(1,cols-j); + + RefMat rm0 = m1; + VERIFY_IS_EQUAL(rm0, m1); + RefDynMat rm1 = m1; + VERIFY_IS_EQUAL(rm1, m1); + RefDynMat rm2 = m1.block(i,j,brows,bcols); + VERIFY_IS_EQUAL(rm2, m1.block(i,j,brows,bcols)); + rm2.setOnes(); + m2.block(i,j,brows,bcols).setOnes(); + VERIFY_IS_EQUAL(m1, m2); + + m2.block(i,j,brows,bcols).setRandom(); + rm2 = m2.block(i,j,brows,bcols); + VERIFY_IS_EQUAL(m1, m2); + + + ConstRefDynMat rm3 = m1.block(i,j,brows,bcols); + m1.block(i,j,brows,bcols) *= 2; + m2.block(i,j,brows,bcols) *= 2; + VERIFY_IS_EQUAL(rm3, m2.block(i,j,brows,bcols)); + RefRealMatWithStride rm4 = m1.real(); + VERIFY_IS_EQUAL(rm4, m2.real()); + rm4.array() += 1; + m2.real().array() += 1; + VERIFY_IS_EQUAL(m1, m2); +} + +template<typename VectorType> void ref_vector(const VectorType& m) +{ + typedef typename VectorType::Index Index; + typedef typename VectorType::Scalar Scalar; + typedef typename VectorType::RealScalar RealScalar; + typedef Matrix<Scalar,Dynamic,1,VectorType::Options> DynMatrixType; + typedef Matrix<Scalar,Dynamic,Dynamic,ColMajor> MatrixType; + typedef Matrix<RealScalar,Dynamic,1,VectorType::Options> RealDynMatrixType; + + typedef Ref<VectorType> RefMat; + typedef Ref<DynMatrixType> RefDynMat; + typedef Ref<const DynMatrixType> ConstRefDynMat; + typedef Ref<RealDynMatrixType , 0, InnerStride<> > RefRealMatWithStride; + typedef Ref<DynMatrixType , 0, InnerStride<> > RefMatWithStride; + + Index size = m.size(); + + VectorType v1 = VectorType::Random(size), + v2 = v1; + MatrixType mat1 = MatrixType::Random(size,size), + mat2 = mat1, + mat3 = MatrixType::Random(size,size); + + Index i = internal::random<Index>(0,size-1); + Index bsize = internal::random<Index>(1,size-i); + + RefMat rm0 = v1; + VERIFY_IS_EQUAL(rm0, v1); + RefDynMat rv1 = v1; + VERIFY_IS_EQUAL(rv1, v1); + RefDynMat rv2 = v1.segment(i,bsize); + VERIFY_IS_EQUAL(rv2, v1.segment(i,bsize)); + rv2.setOnes(); + v2.segment(i,bsize).setOnes(); + VERIFY_IS_EQUAL(v1, v2); + + v2.segment(i,bsize).setRandom(); + rv2 = v2.segment(i,bsize); + VERIFY_IS_EQUAL(v1, v2); + + ConstRefDynMat rm3 = v1.segment(i,bsize); + v1.segment(i,bsize) *= 2; + v2.segment(i,bsize) *= 2; + VERIFY_IS_EQUAL(rm3, v2.segment(i,bsize)); + + RefRealMatWithStride rm4 = v1.real(); + VERIFY_IS_EQUAL(rm4, v2.real()); + rm4.array() += 1; + v2.real().array() += 1; + VERIFY_IS_EQUAL(v1, v2); + + RefMatWithStride rm5 = mat1.row(i).transpose(); + VERIFY_IS_EQUAL(rm5, mat1.row(i).transpose()); + rm5.array() += 1; + mat2.row(i).array() += 1; + VERIFY_IS_EQUAL(mat1, mat2); + rm5.noalias() = rm4.transpose() * mat3; + mat2.row(i) = v2.real().transpose() * mat3; + VERIFY_IS_APPROX(mat1, mat2); +} + +template<typename PlainObjectType> void check_const_correctness(const PlainObjectType&) +{ + // verify that ref-to-const don't have LvalueBit + typedef typename internal::add_const<PlainObjectType>::type ConstPlainObjectType; + VERIFY( !(internal::traits<Ref<ConstPlainObjectType> >::Flags & LvalueBit) ); + VERIFY( !(internal::traits<Ref<ConstPlainObjectType, Aligned> >::Flags & LvalueBit) ); + VERIFY( !(Ref<ConstPlainObjectType>::Flags & LvalueBit) ); + VERIFY( !(Ref<ConstPlainObjectType, Aligned>::Flags & LvalueBit) ); +} + +template<typename B> +EIGEN_DONT_INLINE void call_ref_1(Ref<VectorXf> a, const B &b) { VERIFY_IS_EQUAL(a,b); } +template<typename B> +EIGEN_DONT_INLINE void call_ref_2(const Ref<const VectorXf>& a, const B &b) { VERIFY_IS_EQUAL(a,b); } +template<typename B> +EIGEN_DONT_INLINE void call_ref_3(Ref<VectorXf,0,InnerStride<> > a, const B &b) { VERIFY_IS_EQUAL(a,b); } +template<typename B> +EIGEN_DONT_INLINE void call_ref_4(const Ref<const VectorXf,0,InnerStride<> >& a, const B &b) { VERIFY_IS_EQUAL(a,b); } +template<typename B> +EIGEN_DONT_INLINE void call_ref_5(Ref<MatrixXf,0,OuterStride<> > a, const B &b) { VERIFY_IS_EQUAL(a,b); } +template<typename B> +EIGEN_DONT_INLINE void call_ref_6(const Ref<const MatrixXf,0,OuterStride<> >& a, const B &b) { VERIFY_IS_EQUAL(a,b); } +template<typename B> +EIGEN_DONT_INLINE void call_ref_7(Ref<Matrix<float,Dynamic,3> > a, const B &b) { VERIFY_IS_EQUAL(a,b); } + +void call_ref() +{ + VectorXcf ca = VectorXcf::Random(10); + VectorXf a = VectorXf::Random(10); + RowVectorXf b = RowVectorXf::Random(10); + MatrixXf A = MatrixXf::Random(10,10); + RowVector3f c = RowVector3f::Random(); + const VectorXf& ac(a); + VectorBlock<VectorXf> ab(a,0,3); + const VectorBlock<VectorXf> abc(a,0,3); + + + VERIFY_EVALUATION_COUNT( call_ref_1(a,a), 0); + VERIFY_EVALUATION_COUNT( call_ref_1(b,b.transpose()), 0); +// call_ref_1(ac); // does not compile because ac is const + VERIFY_EVALUATION_COUNT( call_ref_1(ab,ab), 0); + VERIFY_EVALUATION_COUNT( call_ref_1(a.head(4),a.head(4)), 0); + VERIFY_EVALUATION_COUNT( call_ref_1(abc,abc), 0); + VERIFY_EVALUATION_COUNT( call_ref_1(A.col(3),A.col(3)), 0); +// call_ref_1(A.row(3)); // does not compile because innerstride!=1 + VERIFY_EVALUATION_COUNT( call_ref_3(A.row(3),A.row(3).transpose()), 0); + VERIFY_EVALUATION_COUNT( call_ref_4(A.row(3),A.row(3).transpose()), 0); +// call_ref_1(a+a); // does not compile for obvious reason + + MatrixXf tmp = A*A.col(1); + VERIFY_EVALUATION_COUNT( call_ref_2(A*A.col(1), tmp), 1); // evaluated into a temp + VERIFY_EVALUATION_COUNT( call_ref_2(ac.head(5),ac.head(5)), 0); + VERIFY_EVALUATION_COUNT( call_ref_2(ac,ac), 0); + VERIFY_EVALUATION_COUNT( call_ref_2(a,a), 0); + VERIFY_EVALUATION_COUNT( call_ref_2(ab,ab), 0); + VERIFY_EVALUATION_COUNT( call_ref_2(a.head(4),a.head(4)), 0); + tmp = a+a; + VERIFY_EVALUATION_COUNT( call_ref_2(a+a,tmp), 1); // evaluated into a temp + VERIFY_EVALUATION_COUNT( call_ref_2(ca.imag(),ca.imag()), 1); // evaluated into a temp + + VERIFY_EVALUATION_COUNT( call_ref_4(ac.head(5),ac.head(5)), 0); + tmp = a+a; + VERIFY_EVALUATION_COUNT( call_ref_4(a+a,tmp), 1); // evaluated into a temp + VERIFY_EVALUATION_COUNT( call_ref_4(ca.imag(),ca.imag()), 0); + + VERIFY_EVALUATION_COUNT( call_ref_5(a,a), 0); + VERIFY_EVALUATION_COUNT( call_ref_5(a.head(3),a.head(3)), 0); + VERIFY_EVALUATION_COUNT( call_ref_5(A,A), 0); +// call_ref_5(A.transpose()); // does not compile + VERIFY_EVALUATION_COUNT( call_ref_5(A.block(1,1,2,2),A.block(1,1,2,2)), 0); + VERIFY_EVALUATION_COUNT( call_ref_5(b,b), 0); // storage order do not match, but this is a degenerate case that should work + VERIFY_EVALUATION_COUNT( call_ref_5(a.row(3),a.row(3)), 0); + + VERIFY_EVALUATION_COUNT( call_ref_6(a,a), 0); + VERIFY_EVALUATION_COUNT( call_ref_6(a.head(3),a.head(3)), 0); + VERIFY_EVALUATION_COUNT( call_ref_6(A.row(3),A.row(3)), 1); // evaluated into a temp thouth it could be avoided by viewing it as a 1xn matrix + tmp = A+A; + VERIFY_EVALUATION_COUNT( call_ref_6(A+A,tmp), 1); // evaluated into a temp + VERIFY_EVALUATION_COUNT( call_ref_6(A,A), 0); + VERIFY_EVALUATION_COUNT( call_ref_6(A.transpose(),A.transpose()), 1); // evaluated into a temp because the storage orders do not match + VERIFY_EVALUATION_COUNT( call_ref_6(A.block(1,1,2,2),A.block(1,1,2,2)), 0); + + VERIFY_EVALUATION_COUNT( call_ref_7(c,c), 0); +} + +void test_ref() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( ref_vector(Matrix<float, 1, 1>()) ); + CALL_SUBTEST_1( check_const_correctness(Matrix<float, 1, 1>()) ); + CALL_SUBTEST_2( ref_vector(Vector4d()) ); + CALL_SUBTEST_2( check_const_correctness(Matrix4d()) ); + CALL_SUBTEST_3( ref_vector(Vector4cf()) ); + CALL_SUBTEST_4( ref_vector(VectorXcf(8)) ); + CALL_SUBTEST_5( ref_vector(VectorXi(12)) ); + CALL_SUBTEST_5( check_const_correctness(VectorXi(12)) ); + + CALL_SUBTEST_1( ref_matrix(Matrix<float, 1, 1>()) ); + CALL_SUBTEST_2( ref_matrix(Matrix4d()) ); + CALL_SUBTEST_1( ref_matrix(Matrix<float,3,5>()) ); + CALL_SUBTEST_4( ref_matrix(MatrixXcf(internal::random<int>(1,10),internal::random<int>(1,10))) ); + CALL_SUBTEST_4( ref_matrix(Matrix<std::complex<double>,10,15>()) ); + CALL_SUBTEST_5( ref_matrix(MatrixXi(internal::random<int>(1,10),internal::random<int>(1,10))) ); + CALL_SUBTEST_6( call_ref() ); + } +} diff --git a/test/schur_complex.cpp b/test/schur_complex.cpp index a6f66ab02..5e869790f 100644 --- a/test/schur_complex.cpp +++ b/test/schur_complex.cpp @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk> +// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk> // // 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 @@ -47,6 +47,23 @@ template<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTim VERIFY_IS_EQUAL(cs1.matrixT(), cs2.matrixT()); VERIFY_IS_EQUAL(cs1.matrixU(), cs2.matrixU()); + // Test maximum number of iterations + ComplexSchur<MatrixType> cs3; + cs3.setMaxIterations(ComplexSchur<MatrixType>::m_maxIterationsPerRow * size).compute(A); + VERIFY_IS_EQUAL(cs3.info(), Success); + VERIFY_IS_EQUAL(cs3.matrixT(), cs1.matrixT()); + VERIFY_IS_EQUAL(cs3.matrixU(), cs1.matrixU()); + cs3.setMaxIterations(1).compute(A); + VERIFY_IS_EQUAL(cs3.info(), size > 1 ? NoConvergence : Success); + VERIFY_IS_EQUAL(cs3.getMaxIterations(), 1); + + MatrixType Atriangular = A; + Atriangular.template triangularView<StrictlyLower>().setZero(); + cs3.setMaxIterations(1).compute(Atriangular); // triangular matrices do not need any iterations + VERIFY_IS_EQUAL(cs3.info(), Success); + VERIFY_IS_EQUAL(cs3.matrixT(), Atriangular.template cast<ComplexScalar>()); + VERIFY_IS_EQUAL(cs3.matrixU(), ComplexMatrixType::Identity(size, size)); + // Test computation of only T, not U ComplexSchur<MatrixType> csOnlyT(A, false); VERIFY_IS_EQUAL(csOnlyT.info(), Success); diff --git a/test/schur_real.cpp b/test/schur_real.cpp index e6351d94a..36b9c24d1 100644 --- a/test/schur_real.cpp +++ b/test/schur_real.cpp @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk> +// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk> // // 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 @@ -66,6 +66,25 @@ template<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTim VERIFY_IS_EQUAL(rs1.matrixT(), rs2.matrixT()); VERIFY_IS_EQUAL(rs1.matrixU(), rs2.matrixU()); + // Test maximum number of iterations + RealSchur<MatrixType> rs3; + rs3.setMaxIterations(RealSchur<MatrixType>::m_maxIterationsPerRow * size).compute(A); + VERIFY_IS_EQUAL(rs3.info(), Success); + VERIFY_IS_EQUAL(rs3.matrixT(), rs1.matrixT()); + VERIFY_IS_EQUAL(rs3.matrixU(), rs1.matrixU()); + if (size > 2) { + rs3.setMaxIterations(1).compute(A); + VERIFY_IS_EQUAL(rs3.info(), NoConvergence); + VERIFY_IS_EQUAL(rs3.getMaxIterations(), 1); + } + + MatrixType Atriangular = A; + Atriangular.template triangularView<StrictlyLower>().setZero(); + rs3.setMaxIterations(1).compute(Atriangular); // triangular matrices do not need any iterations + VERIFY_IS_EQUAL(rs3.info(), Success); + VERIFY_IS_EQUAL(rs3.matrixT(), Atriangular); + VERIFY_IS_EQUAL(rs3.matrixU(), MatrixType::Identity(size, size)); + // Test computation of only T, not U RealSchur<MatrixType> rsOnlyT(A, false); VERIFY_IS_EQUAL(rsOnlyT.info(), Success); diff --git a/test/selfadjoint.cpp b/test/selfadjoint.cpp index 6d3ec6536..76dab6d64 100644 --- a/test/selfadjoint.cpp +++ b/test/selfadjoint.cpp @@ -16,7 +16,6 @@ template<typename MatrixType> void selfadjoint(const MatrixType& m) { typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits<Scalar>::Real RealScalar; Index rows = m.rows(); Index cols = m.cols(); @@ -47,7 +46,7 @@ void test_selfadjoint() { for(int i = 0; i < g_repeat ; i++) { - int s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE); EIGEN_UNUSED_VARIABLE(s); + int s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE); CALL_SUBTEST_1( selfadjoint(Matrix<float, 1, 1>()) ); CALL_SUBTEST_2( selfadjoint(Matrix<float, 2, 2>()) ); @@ -55,7 +54,7 @@ void test_selfadjoint() CALL_SUBTEST_4( selfadjoint(MatrixXcd(s,s)) ); CALL_SUBTEST_5( selfadjoint(Matrix<float,Dynamic,Dynamic,RowMajor>(s, s)) ); - EIGEN_UNUSED_VARIABLE(s) + TEST_SET_BUT_UNUSED_VARIABLE(s) } CALL_SUBTEST_1( bug_159() ); diff --git a/test/simplicial_cholesky.cpp b/test/simplicial_cholesky.cpp index e93a52e9c..786468421 100644 --- a/test/simplicial_cholesky.cpp +++ b/test/simplicial_cholesky.cpp @@ -11,26 +11,31 @@ template<typename T> void test_simplicial_cholesky_T() { - SimplicialCholesky<SparseMatrix<T>, Lower> chol_colmajor_lower; - SimplicialCholesky<SparseMatrix<T>, Upper> chol_colmajor_upper; - SimplicialLLT<SparseMatrix<T>, Lower> llt_colmajor_lower; - SimplicialLDLT<SparseMatrix<T>, Upper> llt_colmajor_upper; - SimplicialLDLT<SparseMatrix<T>, Lower> ldlt_colmajor_lower; - SimplicialLDLT<SparseMatrix<T>, Upper> ldlt_colmajor_upper; + SimplicialCholesky<SparseMatrix<T>, Lower> chol_colmajor_lower_amd; + SimplicialCholesky<SparseMatrix<T>, Upper> chol_colmajor_upper_amd; + SimplicialLLT<SparseMatrix<T>, Lower> llt_colmajor_lower_amd; + SimplicialLLT<SparseMatrix<T>, Upper> llt_colmajor_upper_amd; + SimplicialLDLT<SparseMatrix<T>, Lower> ldlt_colmajor_lower_amd; + SimplicialLDLT<SparseMatrix<T>, Upper> ldlt_colmajor_upper_amd; + SimplicialLDLT<SparseMatrix<T>, Lower, NaturalOrdering<int> > ldlt_colmajor_lower_nat; + SimplicialLDLT<SparseMatrix<T>, Upper, NaturalOrdering<int> > ldlt_colmajor_upper_nat; - check_sparse_spd_solving(chol_colmajor_lower); - check_sparse_spd_solving(chol_colmajor_upper); - check_sparse_spd_solving(llt_colmajor_lower); - check_sparse_spd_solving(llt_colmajor_upper); - check_sparse_spd_solving(ldlt_colmajor_lower); - check_sparse_spd_solving(ldlt_colmajor_upper); + check_sparse_spd_solving(chol_colmajor_lower_amd); + check_sparse_spd_solving(chol_colmajor_upper_amd); + check_sparse_spd_solving(llt_colmajor_lower_amd); + check_sparse_spd_solving(llt_colmajor_upper_amd); + check_sparse_spd_solving(ldlt_colmajor_lower_amd); + check_sparse_spd_solving(ldlt_colmajor_upper_amd); - check_sparse_spd_determinant(chol_colmajor_lower); - check_sparse_spd_determinant(chol_colmajor_upper); - check_sparse_spd_determinant(llt_colmajor_lower); - check_sparse_spd_determinant(llt_colmajor_upper); - check_sparse_spd_determinant(ldlt_colmajor_lower); - check_sparse_spd_determinant(ldlt_colmajor_upper); + check_sparse_spd_determinant(chol_colmajor_lower_amd); + check_sparse_spd_determinant(chol_colmajor_upper_amd); + check_sparse_spd_determinant(llt_colmajor_lower_amd); + check_sparse_spd_determinant(llt_colmajor_upper_amd); + check_sparse_spd_determinant(ldlt_colmajor_lower_amd); + check_sparse_spd_determinant(ldlt_colmajor_upper_amd); + + check_sparse_spd_solving(ldlt_colmajor_lower_nat); + check_sparse_spd_solving(ldlt_colmajor_upper_nat); } void test_simplicial_cholesky() diff --git a/test/sizeof.cpp b/test/sizeof.cpp index 68463c9b6..d9ad35620 100644 --- a/test/sizeof.cpp +++ b/test/sizeof.cpp @@ -13,7 +13,7 @@ template<typename MatrixType> void verifySizeOf(const MatrixType&) { typedef typename MatrixType::Scalar Scalar; if (MatrixType::RowsAtCompileTime!=Dynamic && MatrixType::ColsAtCompileTime!=Dynamic) - VERIFY(sizeof(MatrixType)==sizeof(Scalar)*size_t(MatrixType::SizeAtCompileTime)); + VERIFY(std::ptrdiff_t(sizeof(MatrixType))==std::ptrdiff_t(sizeof(Scalar))*std::ptrdiff_t(MatrixType::SizeAtCompileTime)); else VERIFY(sizeof(MatrixType)==sizeof(Scalar*) + 2 * sizeof(typename MatrixType::Index)); } diff --git a/test/sparse.h b/test/sparse.h index 4db0004aa..e19a76316 100644 --- a/test/sparse.h +++ b/test/sparse.h @@ -58,18 +58,18 @@ initSparse(double density, Matrix<Scalar,Dynamic,Dynamic,Opt1>& refMat, SparseMatrix<Scalar,Opt2,Index>& sparseMat, int flags = 0, - std::vector<Vector2i>* zeroCoords = 0, - std::vector<Vector2i>* nonzeroCoords = 0) + std::vector<Matrix<Index,2,1> >* zeroCoords = 0, + std::vector<Matrix<Index,2,1> >* nonzeroCoords = 0) { enum { IsRowMajor = SparseMatrix<Scalar,Opt2,Index>::IsRowMajor }; sparseMat.setZero(); //sparseMat.reserve(int(refMat.rows()*refMat.cols()*density)); sparseMat.reserve(VectorXi::Constant(IsRowMajor ? refMat.rows() : refMat.cols(), int((1.5*density)*(IsRowMajor?refMat.cols():refMat.rows())))); - for(int j=0; j<sparseMat.outerSize(); j++) + for(Index j=0; j<sparseMat.outerSize(); j++) { //sparseMat.startVec(j); - for(int i=0; i<sparseMat.innerSize(); i++) + for(Index i=0; i<sparseMat.innerSize(); i++) { int ai(i), aj(j); if(IsRowMajor) @@ -86,18 +86,18 @@ initSparse(double density, v = Scalar(0); if ((flags&ForceRealDiag) && (i==j)) - v = internal::real(v); + v = numext::real(v); if (v!=Scalar(0)) { //sparseMat.insertBackByOuterInner(j,i) = v; sparseMat.insertByOuterInner(j,i) = v; if (nonzeroCoords) - nonzeroCoords->push_back(Vector2i(ai,aj)); + nonzeroCoords->push_back(Matrix<Index,2,1> (ai,aj)); } else if (zeroCoords) { - zeroCoords->push_back(Vector2i(ai,aj)); + zeroCoords->push_back(Matrix<Index,2,1> (ai,aj)); } refMat(ai,aj) = v; } @@ -110,8 +110,8 @@ initSparse(double density, Matrix<Scalar,Dynamic,Dynamic, Opt1>& refMat, DynamicSparseMatrix<Scalar, Opt2, Index>& sparseMat, int flags = 0, - std::vector<Vector2i>* zeroCoords = 0, - std::vector<Vector2i>* nonzeroCoords = 0) + std::vector<Matrix<Index,2,1> >* zeroCoords = 0, + std::vector<Matrix<Index,2,1> >* nonzeroCoords = 0) { enum { IsRowMajor = DynamicSparseMatrix<Scalar,Opt2,Index>::IsRowMajor }; sparseMat.setZero(); @@ -136,17 +136,17 @@ initSparse(double density, v = Scalar(0); if ((flags&ForceRealDiag) && (i==j)) - v = internal::real(v); + v = numext::real(v); if (v!=Scalar(0)) { sparseMat.insertBackByOuterInner(j,i) = v; if (nonzeroCoords) - nonzeroCoords->push_back(Vector2i(ai,aj)); + nonzeroCoords->push_back(Matrix<Index,2,1> (ai,aj)); } else if (zeroCoords) { - zeroCoords->push_back(Vector2i(ai,aj)); + zeroCoords->push_back(Matrix<Index,2,1> (ai,aj)); } refMat(ai,aj) = v; } @@ -154,10 +154,34 @@ initSparse(double density, sparseMat.finalize(); } -template<typename Scalar> void +template<typename Scalar,int Options,typename Index> void initSparse(double density, Matrix<Scalar,Dynamic,1>& refVec, - SparseVector<Scalar>& sparseVec, + SparseVector<Scalar,Options,Index>& sparseVec, + std::vector<int>* zeroCoords = 0, + std::vector<int>* nonzeroCoords = 0) +{ + sparseVec.reserve(int(refVec.size()*density)); + sparseVec.setZero(); + for(Index i=0; i<refVec.size(); i++) + { + Scalar v = (internal::random<double>(0,1) < density) ? internal::random<Scalar>() : Scalar(0); + if (v!=Scalar(0)) + { + sparseVec.insertBack(i) = v; + if (nonzeroCoords) + nonzeroCoords->push_back(i); + } + else if (zeroCoords) + zeroCoords->push_back(i); + refVec[i] = v; + } +} + +template<typename Scalar,int Options,typename Index> void +initSparse(double density, + Matrix<Scalar,1,Dynamic>& refVec, + SparseVector<Scalar,Options,Index>& sparseVec, std::vector<int>* zeroCoords = 0, std::vector<int>* nonzeroCoords = 0) { @@ -178,5 +202,6 @@ initSparse(double density, } } + #include <unsupported/Eigen/SparseExtra> #endif // EIGEN_TESTSPARSE_H diff --git a/test/sparseLM.cpp b/test/sparseLM.cpp new file mode 100644 index 000000000..8e148f9bc --- /dev/null +++ b/test/sparseLM.cpp @@ -0,0 +1,176 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Desire Nuentsa <desire.nuentsa_wakam@inria.fr> +// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr> +// +// 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/. +#include <iostream> +#include <fstream> +#include <iomanip> + +#include "main.h" +#include <Eigen/LevenbergMarquardt> + +using namespace std; +using namespace Eigen; + +template <typename Scalar> +struct sparseGaussianTest : SparseFunctor<Scalar, int> +{ + typedef Matrix<Scalar,Dynamic,1> VectorType; + typedef SparseFunctor<Scalar,int> Base; + typedef typename Base::JacobianType JacobianType; + sparseGaussianTest(int inputs, int values) : SparseFunctor<Scalar,int>(inputs,values) + { } + + VectorType model(const VectorType& uv, VectorType& x) + { + VectorType y; //Change this to use expression template + int m = Base::values(); + int n = Base::inputs(); + eigen_assert(uv.size()%2 == 0); + eigen_assert(uv.size() == n); + eigen_assert(x.size() == m); + y.setZero(m); + int half = n/2; + VectorBlock<const VectorType> u(uv, 0, half); + VectorBlock<const VectorType> v(uv, half, half); + Scalar coeff; + for (int j = 0; j < m; j++) + { + for (int i = 0; i < half; i++) + { + coeff = (x(j)-i)/v(i); + coeff *= coeff; + if (coeff < 1. && coeff > 0.) + y(j) += u(i)*std::pow((1-coeff), 2); + } + } + return y; + } + void initPoints(VectorType& uv_ref, VectorType& x) + { + m_x = x; + m_y = this->model(uv_ref,x); + } + int operator()(const VectorType& uv, VectorType& fvec) + { + int m = Base::values(); + int n = Base::inputs(); + eigen_assert(uv.size()%2 == 0); + eigen_assert(uv.size() == n); + int half = n/2; + VectorBlock<const VectorType> u(uv, 0, half); + VectorBlock<const VectorType> v(uv, half, half); + fvec = m_y; + Scalar coeff; + for (int j = 0; j < m; j++) + { + for (int i = 0; i < half; i++) + { + coeff = (m_x(j)-i)/v(i); + coeff *= coeff; + if (coeff < 1. && coeff > 0.) + fvec(j) -= u(i)*std::pow((1-coeff), 2); + } + } + return 0; + } + + int df(const VectorType& uv, JacobianType& fjac) + { + int m = Base::values(); + int n = Base::inputs(); + eigen_assert(n == uv.size()); + eigen_assert(fjac.rows() == m); + eigen_assert(fjac.cols() == n); + int half = n/2; + VectorBlock<const VectorType> u(uv, 0, half); + VectorBlock<const VectorType> v(uv, half, half); + Scalar coeff; + + //Derivatives with respect to u + for (int col = 0; col < half; col++) + { + for (int row = 0; row < m; row++) + { + coeff = (m_x(row)-col)/v(col); + coeff = coeff*coeff; + if(coeff < 1. && coeff > 0.) + { + fjac.coeffRef(row,col) = -(1-coeff)*(1-coeff); + } + } + } + //Derivatives with respect to v + for (int col = 0; col < half; col++) + { + for (int row = 0; row < m; row++) + { + coeff = (m_x(row)-col)/v(col); + coeff = coeff*coeff; + if(coeff < 1. && coeff > 0.) + { + fjac.coeffRef(row,col+half) = -4 * (u(col)/v(col))*coeff*(1-coeff); + } + } + } + return 0; + } + + VectorType m_x, m_y; //Data points +}; + + +template<typename T> +void test_sparseLM_T() +{ + typedef Matrix<T,Dynamic,1> VectorType; + + int inputs = 10; + int values = 2000; + sparseGaussianTest<T> sparse_gaussian(inputs, values); + VectorType uv(inputs),uv_ref(inputs); + VectorType x(values); + // Generate the reference solution + uv_ref << -2, 1, 4 ,8, 6, 1.8, 1.2, 1.1, 1.9 , 3; + //Generate the reference data points + x.setRandom(); + x = 10*x; + x.array() += 10; + sparse_gaussian.initPoints(uv_ref, x); + + + // Generate the initial parameters + VectorBlock<VectorType> u(uv, 0, inputs/2); + VectorBlock<VectorType> v(uv, inputs/2, inputs/2); + v.setOnes(); + //Generate u or Solve for u from v + u.setOnes(); + + // Solve the optimization problem + LevenbergMarquardt<sparseGaussianTest<T> > lm(sparse_gaussian); + int info; +// info = lm.minimize(uv); + + VERIFY_IS_EQUAL(info,1); + // Do a step by step solution and save the residual + int maxiter = 200; + int iter = 0; + MatrixXd Err(values, maxiter); + MatrixXd Mod(values, maxiter); + LevenbergMarquardtSpace::Status status; + status = lm.minimizeInit(uv); + if (status==LevenbergMarquardtSpace::ImproperInputParameters) + return ; + +} +void test_sparseLM() +{ + CALL_SUBTEST_1(test_sparseLM_T<double>()); + + // CALL_SUBTEST_2(test_sparseLM_T<std::complex<double>()); +} diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp index 8897a9dca..498ecfe29 100644 --- a/test/sparse_basic.cpp +++ b/test/sparse_basic.cpp @@ -3,6 +3,7 @@ // // Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr> // Copyright (C) 2008 Daniel Gomez Ferro <dgomezferro@gmail.com> +// Copyright (C) 2013 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> // // 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 @@ -13,7 +14,8 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& ref) { typedef typename SparseMatrixType::Index Index; - + typedef Matrix<Index,2,1> Vector2; + const Index rows = ref.rows(); const Index cols = ref.cols(); typedef typename SparseMatrixType::Scalar Scalar; @@ -24,71 +26,77 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re typedef Matrix<Scalar,Dynamic,1> DenseVector; Scalar eps = 1e-6; - SparseMatrixType m(rows, cols); - DenseMatrix refMat = DenseMatrix::Zero(rows, cols); - DenseVector vec1 = DenseVector::Random(rows); Scalar s1 = internal::random<Scalar>(); - - std::vector<Vector2i> zeroCoords; - std::vector<Vector2i> nonzeroCoords; - initSparse<Scalar>(density, refMat, m, 0, &zeroCoords, &nonzeroCoords); - - if (zeroCoords.size()==0 || nonzeroCoords.size()==0) - return; - - // test coeff and coeffRef - for (int i=0; i<(int)zeroCoords.size(); ++i) { - VERIFY_IS_MUCH_SMALLER_THAN( m.coeff(zeroCoords[i].x(),zeroCoords[i].y()), eps ); - if(internal::is_same<SparseMatrixType,SparseMatrix<Scalar,Flags> >::value) - VERIFY_RAISES_ASSERT( m.coeffRef(zeroCoords[0].x(),zeroCoords[0].y()) = 5 ); - } - VERIFY_IS_APPROX(m, refMat); + SparseMatrixType m(rows, cols); + DenseMatrix refMat = DenseMatrix::Zero(rows, cols); + DenseVector vec1 = DenseVector::Random(rows); - m.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5); - refMat.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5); + std::vector<Vector2> zeroCoords; + std::vector<Vector2> nonzeroCoords; + initSparse<Scalar>(density, refMat, m, 0, &zeroCoords, &nonzeroCoords); - VERIFY_IS_APPROX(m, refMat); - /* - // test InnerIterators and Block expressions - for (int t=0; t<10; ++t) - { - int j = internal::random<int>(0,cols-1); - int i = internal::random<int>(0,rows-1); - int w = internal::random<int>(1,cols-j-1); - int h = internal::random<int>(1,rows-i-1); + if (zeroCoords.size()==0 || nonzeroCoords.size()==0) + return; -// VERIFY_IS_APPROX(m.block(i,j,h,w), refMat.block(i,j,h,w)); - for(int c=0; c<w; c++) + // test coeff and coeffRef + for (int i=0; i<(int)zeroCoords.size(); ++i) { - VERIFY_IS_APPROX(m.block(i,j,h,w).col(c), refMat.block(i,j,h,w).col(c)); - for(int r=0; r<h; r++) + VERIFY_IS_MUCH_SMALLER_THAN( m.coeff(zeroCoords[i].x(),zeroCoords[i].y()), eps ); + if(internal::is_same<SparseMatrixType,SparseMatrix<Scalar,Flags> >::value) + VERIFY_RAISES_ASSERT( m.coeffRef(zeroCoords[0].x(),zeroCoords[0].y()) = 5 ); + } + VERIFY_IS_APPROX(m, refMat); + + m.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5); + refMat.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5); + + VERIFY_IS_APPROX(m, refMat); + /* + // test InnerIterators and Block expressions + for (int t=0; t<10; ++t) { -// VERIFY_IS_APPROX(m.block(i,j,h,w).col(c).coeff(r), refMat.block(i,j,h,w).col(c).coeff(r)); + int j = internal::random<int>(0,cols-1); + int i = internal::random<int>(0,rows-1); + int w = internal::random<int>(1,cols-j-1); + int h = internal::random<int>(1,rows-i-1); + + // VERIFY_IS_APPROX(m.block(i,j,h,w), refMat.block(i,j,h,w)); + for(int c=0; c<w; c++) + { + VERIFY_IS_APPROX(m.block(i,j,h,w).col(c), refMat.block(i,j,h,w).col(c)); + for(int r=0; r<h; r++) + { + // VERIFY_IS_APPROX(m.block(i,j,h,w).col(c).coeff(r), refMat.block(i,j,h,w).col(c).coeff(r)); + } + } + // for(int r=0; r<h; r++) + // { + // VERIFY_IS_APPROX(m.block(i,j,h,w).row(r), refMat.block(i,j,h,w).row(r)); + // for(int c=0; c<w; c++) + // { + // VERIFY_IS_APPROX(m.block(i,j,h,w).row(r).coeff(c), refMat.block(i,j,h,w).row(r).coeff(c)); + // } + // } } - } -// for(int r=0; r<h; r++) -// { -// VERIFY_IS_APPROX(m.block(i,j,h,w).row(r), refMat.block(i,j,h,w).row(r)); -// for(int c=0; c<w; c++) -// { -// VERIFY_IS_APPROX(m.block(i,j,h,w).row(r).coeff(c), refMat.block(i,j,h,w).row(r).coeff(c)); -// } -// } - } - for(int c=0; c<cols; c++) - { - VERIFY_IS_APPROX(m.col(c) + m.col(c), (m + m).col(c)); - VERIFY_IS_APPROX(m.col(c) + m.col(c), refMat.col(c) + refMat.col(c)); - } + for(int c=0; c<cols; c++) + { + VERIFY_IS_APPROX(m.col(c) + m.col(c), (m + m).col(c)); + VERIFY_IS_APPROX(m.col(c) + m.col(c), refMat.col(c) + refMat.col(c)); + } - for(int r=0; r<rows; r++) - { - VERIFY_IS_APPROX(m.row(r) + m.row(r), (m + m).row(r)); - VERIFY_IS_APPROX(m.row(r) + m.row(r), refMat.row(r) + refMat.row(r)); - } - */ + for(int r=0; r<rows; r++) + { + VERIFY_IS_APPROX(m.row(r) + m.row(r), (m + m).row(r)); + VERIFY_IS_APPROX(m.row(r) + m.row(r), refMat.row(r) + refMat.row(r)); + } + */ + + // test assertion + VERIFY_RAISES_ASSERT( m.coeffRef(-1,1) = 0 ); + VERIFY_RAISES_ASSERT( m.coeffRef(0,m.cols()) = 0 ); + } // test insert (inner random) { @@ -97,11 +105,11 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re SparseMatrixType m2(rows,cols); if(internal::random<int>()%2) m2.reserve(VectorXi::Constant(m2.outerSize(), 2)); - for (int j=0; j<cols; ++j) + for (Index j=0; j<cols; ++j) { - for (int k=0; k<rows/2; ++k) + for (Index k=0; k<rows/2; ++k) { - int i = internal::random<int>(0,rows-1); + Index i = internal::random<Index>(0,rows-1); if (m1.coeff(i,j)==Scalar(0)) m2.insert(i,j) = m1(i,j) = internal::random<Scalar>(); } @@ -119,8 +127,8 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re m2.reserve(VectorXi::Constant(m2.outerSize(), 2)); for (int k=0; k<rows*cols; ++k) { - int i = internal::random<int>(0,rows-1); - int j = internal::random<int>(0,cols-1); + Index i = internal::random<Index>(0,rows-1); + Index j = internal::random<Index>(0,cols-1); if ((m1.coeff(i,j)==Scalar(0)) && (internal::random<int>()%2)) m2.insert(i,j) = m1(i,j) = internal::random<Scalar>(); else @@ -143,8 +151,8 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re m2.reserve(r); for (int k=0; k<rows*cols; ++k) { - int i = internal::random<int>(0,rows-1); - int j = internal::random<int>(0,cols-1); + Index i = internal::random<Index>(0,rows-1); + Index j = internal::random<Index>(0,cols-1); if (m1.coeff(i,j)==Scalar(0)) m2.insert(i,j) = m1(i,j) = internal::random<Scalar>(); if(mode==3) @@ -155,6 +163,80 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re VERIFY_IS_APPROX(m2,m1); } + // test innerVector() + { + DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); + SparseMatrixType m2(rows, rows); + initSparse<Scalar>(density, refMat2, m2); + Index j0 = internal::random<Index>(0,rows-1); + Index j1 = internal::random<Index>(0,rows-1); + if(SparseMatrixType::IsRowMajor) + VERIFY_IS_APPROX(m2.innerVector(j0), refMat2.row(j0)); + else + VERIFY_IS_APPROX(m2.innerVector(j0), refMat2.col(j0)); + + if(SparseMatrixType::IsRowMajor) + VERIFY_IS_APPROX(m2.innerVector(j0)+m2.innerVector(j1), refMat2.row(j0)+refMat2.row(j1)); + else + VERIFY_IS_APPROX(m2.innerVector(j0)+m2.innerVector(j1), refMat2.col(j0)+refMat2.col(j1)); + + SparseMatrixType m3(rows,rows); + m3.reserve(VectorXi::Constant(rows,rows/2)); + for(Index j=0; j<rows; ++j) + for(Index k=0; k<j; ++k) + m3.insertByOuterInner(j,k) = k+1; + for(Index j=0; j<rows; ++j) + { + VERIFY(j==numext::real(m3.innerVector(j).nonZeros())); + if(j>0) + VERIFY(j==numext::real(m3.innerVector(j).lastCoeff())); + } + m3.makeCompressed(); + for(Index j=0; j<rows; ++j) + { + VERIFY(j==numext::real(m3.innerVector(j).nonZeros())); + if(j>0) + VERIFY(j==numext::real(m3.innerVector(j).lastCoeff())); + } + + //m2.innerVector(j0) = 2*m2.innerVector(j1); + //refMat2.col(j0) = 2*refMat2.col(j1); + //VERIFY_IS_APPROX(m2, refMat2); + } + + // test innerVectors() + { + DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); + SparseMatrixType m2(rows, rows); + initSparse<Scalar>(density, refMat2, m2); + if(internal::random<float>(0,1)>0.5) m2.makeCompressed(); + + Index j0 = internal::random<Index>(0,rows-2); + Index j1 = internal::random<Index>(0,rows-2); + Index n0 = internal::random<Index>(1,rows-(std::max)(j0,j1)); + if(SparseMatrixType::IsRowMajor) + VERIFY_IS_APPROX(m2.innerVectors(j0,n0), refMat2.block(j0,0,n0,cols)); + else + VERIFY_IS_APPROX(m2.innerVectors(j0,n0), refMat2.block(0,j0,rows,n0)); + if(SparseMatrixType::IsRowMajor) + VERIFY_IS_APPROX(m2.innerVectors(j0,n0)+m2.innerVectors(j1,n0), + refMat2.middleRows(j0,n0)+refMat2.middleRows(j1,n0)); + else + VERIFY_IS_APPROX(m2.innerVectors(j0,n0)+m2.innerVectors(j1,n0), + refMat2.block(0,j0,rows,n0)+refMat2.block(0,j1,rows,n0)); + + VERIFY_IS_APPROX(m2, refMat2); + + m2.innerVectors(j0,n0) = m2.innerVectors(j0,n0) + m2.innerVectors(j1,n0); + if(SparseMatrixType::IsRowMajor) + refMat2.middleRows(j0,n0) = (refMat2.middleRows(j0,n0) + refMat2.middleRows(j1,n0)).eval(); + else + refMat2.middleCols(j0,n0) = (refMat2.middleCols(j0,n0) + refMat2.middleCols(j1,n0)).eval(); + + VERIFY_IS_APPROX(m2, refMat2); + + } + // test basic computations { DenseMatrix refM1 = DenseMatrix::Zero(rows, rows); @@ -193,6 +275,12 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re // sparse cwise* dense VERIFY_IS_APPROX(m3.cwiseProduct(refM4), refM3.cwiseProduct(refM4)); // VERIFY_IS_APPROX(m3.cwise()/refM4, refM3.cwise()/refM4); + + // test aliasing + VERIFY_IS_APPROX((m1 = -m1), (refM1 = -refM1)); + VERIFY_IS_APPROX((m1 = m1.transpose()), (refM1 = refM1.transpose().eval())); + VERIFY_IS_APPROX((m1 = -m1.transpose()), (refM1 = -refM1.transpose().eval())); + VERIFY_IS_APPROX((m1 += -m1), (refM1 += -refM1)); } // test transpose @@ -206,67 +294,38 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re VERIFY_IS_APPROX(SparseMatrixType(m2.adjoint()), refMat2.adjoint()); } - // test innerVector() - { - DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); - SparseMatrixType m2(rows, rows); - initSparse<Scalar>(density, refMat2, m2); - int j0 = internal::random<int>(0,rows-1); - int j1 = internal::random<int>(0,rows-1); - if(SparseMatrixType::IsRowMajor) - VERIFY_IS_APPROX(m2.innerVector(j0), refMat2.row(j0)); - else - VERIFY_IS_APPROX(m2.innerVector(j0), refMat2.col(j0)); - - if(SparseMatrixType::IsRowMajor) - VERIFY_IS_APPROX(m2.innerVector(j0)+m2.innerVector(j1), refMat2.row(j0)+refMat2.row(j1)); - else - VERIFY_IS_APPROX(m2.innerVector(j0)+m2.innerVector(j1), refMat2.col(j0)+refMat2.col(j1)); - - SparseMatrixType m3(rows,rows); - m3.reserve(VectorXi::Constant(rows,rows/2)); - for(int j=0; j<rows; ++j) - for(int k=0; k<j; ++k) - m3.insertByOuterInner(j,k) = k+1; - for(int j=0; j<rows; ++j) - { - VERIFY(j==internal::real(m3.innerVector(j).nonZeros())); - if(j>0) - VERIFY(j==internal::real(m3.innerVector(j).lastCoeff())); - } - m3.makeCompressed(); - for(int j=0; j<rows; ++j) - { - VERIFY(j==internal::real(m3.innerVector(j).nonZeros())); - if(j>0) - VERIFY(j==internal::real(m3.innerVector(j).lastCoeff())); - } - - //m2.innerVector(j0) = 2*m2.innerVector(j1); - //refMat2.col(j0) = 2*refMat2.col(j1); - //VERIFY_IS_APPROX(m2, refMat2); - } - - // test innerVectors() + + + // test generic blocks { DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); SparseMatrixType m2(rows, rows); initSparse<Scalar>(density, refMat2, m2); - int j0 = internal::random<int>(0,rows-2); - int j1 = internal::random<int>(0,rows-2); - int n0 = internal::random<int>(1,rows-(std::max)(j0,j1)); + Index j0 = internal::random<Index>(0,rows-2); + Index j1 = internal::random<Index>(0,rows-2); + Index n0 = internal::random<Index>(1,rows-(std::max)(j0,j1)); if(SparseMatrixType::IsRowMajor) - VERIFY_IS_APPROX(m2.innerVectors(j0,n0), refMat2.block(j0,0,n0,cols)); + VERIFY_IS_APPROX(m2.block(j0,0,n0,cols), refMat2.block(j0,0,n0,cols)); else - VERIFY_IS_APPROX(m2.innerVectors(j0,n0), refMat2.block(0,j0,rows,n0)); + VERIFY_IS_APPROX(m2.block(0,j0,rows,n0), refMat2.block(0,j0,rows,n0)); + if(SparseMatrixType::IsRowMajor) - VERIFY_IS_APPROX(m2.innerVectors(j0,n0)+m2.innerVectors(j1,n0), + VERIFY_IS_APPROX(m2.block(j0,0,n0,cols)+m2.block(j1,0,n0,cols), refMat2.block(j0,0,n0,cols)+refMat2.block(j1,0,n0,cols)); else - VERIFY_IS_APPROX(m2.innerVectors(j0,n0)+m2.innerVectors(j1,n0), + VERIFY_IS_APPROX(m2.block(0,j0,rows,n0)+m2.block(0,j1,rows,n0), refMat2.block(0,j0,rows,n0)+refMat2.block(0,j1,rows,n0)); - //m2.innerVectors(j0,n0) = m2.innerVectors(j0,n0) + m2.innerVectors(j1,n0); - //refMat2.block(0,j0,rows,n0) = refMat2.block(0,j0,rows,n0) + refMat2.block(0,j1,rows,n0); + + Index i = internal::random<Index>(0,m2.outerSize()-1); + if(SparseMatrixType::IsRowMajor) { + m2.innerVector(i) = m2.innerVector(i) * s1; + refMat2.row(i) = refMat2.row(i) * s1; + VERIFY_IS_APPROX(m2,refMat2); + } else { + m2.innerVector(i) = m2.innerVector(i) * s1; + refMat2.col(i) = refMat2.col(i) * s1; + VERIFY_IS_APPROX(m2,refMat2); + } } // test prune @@ -276,10 +335,10 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re refM2.setZero(); int countFalseNonZero = 0; int countTrueNonZero = 0; - for (int j=0; j<m2.outerSize(); ++j) + for (Index j=0; j<m2.outerSize(); ++j) { m2.startVec(j); - for (int i=0; i<m2.innerSize(); ++i) + for (Index i=0; i<m2.innerSize(); ++i) { float x = internal::random<float>(0,1); if (x<0.1) @@ -320,8 +379,8 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re refMat.setZero(); for(int i=0;i<ntriplets;++i) { - int r = internal::random<int>(0,rows-1); - int c = internal::random<int>(0,cols-1); + Index r = internal::random<Index>(0,rows-1); + Index c = internal::random<Index>(0,cols-1); Scalar v = internal::random<Scalar>(); triplets.push_back(TripletType(r,c,v)); refMat(r,c) += v; @@ -351,6 +410,14 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re refMat3 = refMat2.template triangularView<UnitLower>(); m3 = m2.template triangularView<UnitLower>(); VERIFY_IS_APPROX(m3, refMat3); + + refMat3 = refMat2.template triangularView<StrictlyUpper>(); + m3 = m2.template triangularView<StrictlyUpper>(); + VERIFY_IS_APPROX(m3, refMat3); + + refMat3 = refMat2.template triangularView<StrictlyLower>(); + m3 = m2.template triangularView<StrictlyLower>(); + VERIFY_IS_APPROX(m3, refMat3); } // test selfadjointView @@ -379,17 +446,64 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re initSparse<Scalar>(density, refMat2, m2); VERIFY_IS_APPROX(m2.diagonal(), refMat2.diagonal().eval()); } + + // test conservative resize + { + std::vector< std::pair<Index,Index> > inc; + inc.push_back(std::pair<Index,Index>(-3,-2)); + inc.push_back(std::pair<Index,Index>(0,0)); + inc.push_back(std::pair<Index,Index>(3,2)); + inc.push_back(std::pair<Index,Index>(3,0)); + inc.push_back(std::pair<Index,Index>(0,3)); + + for(size_t i = 0; i< inc.size(); i++) { + Index incRows = inc[i].first; + Index incCols = inc[i].second; + SparseMatrixType m1(rows, cols); + DenseMatrix refMat1 = DenseMatrix::Zero(rows, cols); + initSparse<Scalar>(density, refMat1, m1); + + m1.conservativeResize(rows+incRows, cols+incCols); + refMat1.conservativeResize(rows+incRows, cols+incCols); + if (incRows > 0) refMat1.bottomRows(incRows).setZero(); + if (incCols > 0) refMat1.rightCols(incCols).setZero(); + + VERIFY_IS_APPROX(m1, refMat1); + + // Insert new values + if (incRows > 0) + m1.insert(m1.rows()-1, 0) = refMat1(refMat1.rows()-1, 0) = 1; + if (incCols > 0) + m1.insert(0, m1.cols()-1) = refMat1(0, refMat1.cols()-1) = 1; + + VERIFY_IS_APPROX(m1, refMat1); + + + } + } + + // test Identity matrix + { + DenseMatrix refMat1 = DenseMatrix::Identity(rows, rows); + SparseMatrixType m1(rows, rows); + m1.setIdentity(); + VERIFY_IS_APPROX(m1, refMat1); + } } void test_sparse_basic() { for(int i = 0; i < g_repeat; i++) { int s = Eigen::internal::random<int>(1,50); + EIGEN_UNUSED_VARIABLE(s); CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double>(8, 8)) )); CALL_SUBTEST_2(( sparse_basic(SparseMatrix<std::complex<double>, ColMajor>(s, s)) )); CALL_SUBTEST_2(( sparse_basic(SparseMatrix<std::complex<double>, RowMajor>(s, s)) )); CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double>(s, s)) )); CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double,ColMajor,long int>(s, s)) )); CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double,RowMajor,long int>(s, s)) )); + + CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double,ColMajor,short int>(short(s), short(s))) )); + CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double,RowMajor,short int>(short(s), short(s))) )); } } diff --git a/test/sparse_product.cpp b/test/sparse_product.cpp index 17a955c9d..a2ea9d5b7 100644 --- a/test/sparse_product.cpp +++ b/test/sparse_product.cpp @@ -13,8 +13,9 @@ template<typename SparseMatrixType, typename DenseMatrix, bool IsRowMajor=Sparse template<typename SparseMatrixType, typename DenseMatrix> struct test_outer<SparseMatrixType,DenseMatrix,false> { static void run(SparseMatrixType& m2, SparseMatrixType& m4, DenseMatrix& refMat2, DenseMatrix& refMat4) { - int c = internal::random(0,m2.cols()-1); - int c1 = internal::random(0,m2.cols()-1); + typedef typename SparseMatrixType::Index Index; + Index c = internal::random<Index>(0,m2.cols()-1); + Index c1 = internal::random<Index>(0,m2.cols()-1); VERIFY_IS_APPROX(m4=m2.col(c)*refMat2.col(c1).transpose(), refMat4=refMat2.col(c)*refMat2.col(c1).transpose()); VERIFY_IS_APPROX(m4=refMat2.col(c1)*m2.col(c).transpose(), refMat4=refMat2.col(c1)*refMat2.col(c).transpose()); } @@ -22,8 +23,9 @@ template<typename SparseMatrixType, typename DenseMatrix> struct test_outer<Spar template<typename SparseMatrixType, typename DenseMatrix> struct test_outer<SparseMatrixType,DenseMatrix,true> { static void run(SparseMatrixType& m2, SparseMatrixType& m4, DenseMatrix& refMat2, DenseMatrix& refMat4) { - int r = internal::random(0,m2.rows()-1); - int c1 = internal::random(0,m2.cols()-1); + typedef typename SparseMatrixType::Index Index; + Index r = internal::random<Index>(0,m2.rows()-1); + Index c1 = internal::random<Index>(0,m2.cols()-1); VERIFY_IS_APPROX(m4=m2.row(r).transpose()*refMat2.col(c1).transpose(), refMat4=refMat2.row(r).transpose()*refMat2.col(c1).transpose()); VERIFY_IS_APPROX(m4=refMat2.col(c1)*m2.row(r), refMat4=refMat2.col(c1)*refMat2.row(r)); } @@ -37,15 +39,18 @@ template<typename SparseMatrixType> void sparse_product() { typedef typename SparseMatrixType::Index Index; Index n = 100; - const Index rows = internal::random<int>(1,n); - const Index cols = internal::random<int>(1,n); - const Index depth = internal::random<int>(1,n); + const Index rows = internal::random<Index>(1,n); + const Index cols = internal::random<Index>(1,n); + const Index depth = internal::random<Index>(1,n); typedef typename SparseMatrixType::Scalar Scalar; enum { Flags = SparseMatrixType::Flags }; - double density = (std::max)(8./(rows*cols), 0.01); + double density = (std::max)(8./(rows*cols), 0.1); typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; typedef Matrix<Scalar,Dynamic,1> DenseVector; + typedef Matrix<Scalar,1,Dynamic> RowDenseVector; + typedef SparseVector<Scalar,0,Index> ColSpVector; + typedef SparseVector<Scalar,RowMajor,Index> RowSpVector; Scalar s1 = internal::random<Scalar>(); Scalar s2 = internal::random<Scalar>(); @@ -117,21 +122,54 @@ template<typename SparseMatrixType> void sparse_product() test_outer<SparseMatrixType,DenseMatrix>::run(m2,m4,refMat2,refMat4); VERIFY_IS_APPROX(m6=m6*m6, refMat6=refMat6*refMat6); + + // sparse matrix * sparse vector + ColSpVector cv0(cols), cv1; + DenseVector dcv0(cols), dcv1; + initSparse(2*density,dcv0, cv0); + + RowSpVector rv0(depth), rv1; + RowDenseVector drv0(depth), drv1(rv1); + initSparse(2*density,drv0, rv0); + + VERIFY_IS_APPROX(cv1=rv0*m3, dcv1=drv0*refMat3); + VERIFY_IS_APPROX(rv1=rv0*m3, drv1=drv0*refMat3); + VERIFY_IS_APPROX(cv1=m3*cv0, dcv1=refMat3*dcv0); + VERIFY_IS_APPROX(cv1=m3t.adjoint()*cv0, dcv1=refMat3t.adjoint()*dcv0); + VERIFY_IS_APPROX(rv1=m3*cv0, drv1=refMat3*dcv0); } - + // test matrix - diagonal product { - DenseMatrix refM2 = DenseMatrix::Zero(rows, rows); - DenseMatrix refM3 = DenseMatrix::Zero(rows, rows); - DiagonalMatrix<Scalar,Dynamic> d1(DenseVector::Random(rows)); - SparseMatrixType m2(rows, rows); - SparseMatrixType m3(rows, rows); + DenseMatrix refM2 = DenseMatrix::Zero(rows, cols); + DenseMatrix refM3 = DenseMatrix::Zero(rows, cols); + DenseMatrix d3 = DenseMatrix::Zero(rows, cols); + DiagonalMatrix<Scalar,Dynamic> d1(DenseVector::Random(cols)); + DiagonalMatrix<Scalar,Dynamic> d2(DenseVector::Random(rows)); + SparseMatrixType m2(rows, cols); + SparseMatrixType m3(rows, cols); initSparse<Scalar>(density, refM2, m2); initSparse<Scalar>(density, refM3, m3); VERIFY_IS_APPROX(m3=m2*d1, refM3=refM2*d1); - VERIFY_IS_APPROX(m3=m2.transpose()*d1, refM3=refM2.transpose()*d1); - VERIFY_IS_APPROX(m3=d1*m2, refM3=d1*refM2); - VERIFY_IS_APPROX(m3=d1*m2.transpose(), refM3=d1 * refM2.transpose()); + VERIFY_IS_APPROX(m3=m2.transpose()*d2, refM3=refM2.transpose()*d2); + VERIFY_IS_APPROX(m3=d2*m2, refM3=d2*refM2); + VERIFY_IS_APPROX(m3=d1*m2.transpose(), refM3=d1*refM2.transpose()); + + // also check with a SparseWrapper: + DenseVector v1 = DenseVector::Random(cols); + DenseVector v2 = DenseVector::Random(rows); + VERIFY_IS_APPROX(m3=m2*v1.asDiagonal(), refM3=refM2*v1.asDiagonal()); + VERIFY_IS_APPROX(m3=m2.transpose()*v2.asDiagonal(), refM3=refM2.transpose()*v2.asDiagonal()); + VERIFY_IS_APPROX(m3=v2.asDiagonal()*m2, refM3=v2.asDiagonal()*refM2); + VERIFY_IS_APPROX(m3=v1.asDiagonal()*m2.transpose(), refM3=v1.asDiagonal()*refM2.transpose()); + + VERIFY_IS_APPROX(m3=v2.asDiagonal()*m2*v1.asDiagonal(), refM3=v2.asDiagonal()*refM2*v1.asDiagonal()); + + // evaluate to a dense matrix to check the .row() and .col() iterator functions + VERIFY_IS_APPROX(d3=m2*d1, refM3=refM2*d1); + VERIFY_IS_APPROX(d3=m2.transpose()*d2, refM3=refM2.transpose()*d2); + VERIFY_IS_APPROX(d3=d2*m2, refM3=d2*refM2); + VERIFY_IS_APPROX(d3=d1*m2.transpose(), refM3=d1*refM2.transpose()); } // test self adjoint products @@ -167,7 +205,16 @@ template<typename SparseMatrixType> void sparse_product() VERIFY_IS_APPROX(x=mUp.template selfadjointView<Upper>()*b, refX=refS*b); VERIFY_IS_APPROX(x=mLo.template selfadjointView<Lower>()*b, refX=refS*b); VERIFY_IS_APPROX(x=mS.template selfadjointView<Upper|Lower>()*b, refX=refS*b); + + // sparse selfadjointView * sparse + SparseMatrixType mSres(rows,rows); + VERIFY_IS_APPROX(mSres = mLo.template selfadjointView<Lower>()*mS, + refX = refLo.template selfadjointView<Lower>()*refS); + // sparse * sparse selfadjointview + VERIFY_IS_APPROX(mSres = mS * mLo.template selfadjointView<Lower>(), + refX = refS * refLo.template selfadjointView<Lower>()); } + } // New test for Bug in SparseTimeDenseProduct @@ -199,6 +246,7 @@ void test_sparse_product() CALL_SUBTEST_1( (sparse_product<SparseMatrix<double,RowMajor> >()) ); CALL_SUBTEST_2( (sparse_product<SparseMatrix<std::complex<double>, ColMajor > >()) ); CALL_SUBTEST_2( (sparse_product<SparseMatrix<std::complex<double>, RowMajor > >()) ); + CALL_SUBTEST_3( (sparse_product<SparseMatrix<float,ColMajor,long int> >()) ); CALL_SUBTEST_4( (sparse_product_regression_test<SparseMatrix<double,RowMajor>, Matrix<double, Dynamic, Dynamic, RowMajor> >()) ); } } diff --git a/test/sparse_solver.h b/test/sparse_solver.h index 75fa85082..d84aff070 100644 --- a/test/sparse_solver.h +++ b/test/sparse_solver.h @@ -17,54 +17,54 @@ void check_sparse_solving(Solver& solver, const typename Solver::MatrixType& A, typedef typename Mat::Scalar Scalar; DenseRhs refX = dA.lu().solve(db); - - Rhs x(b.rows(), b.cols()); - Rhs oldb = b; - - solver.compute(A); - if (solver.info() != Success) { - std::cerr << "sparse solver testing: factorization failed (check_sparse_solving)\n"; - exit(0); - return; - } - x = solver.solve(b); - if (solver.info() != Success) - { - std::cerr << "sparse solver testing: solving failed\n"; - return; - } - VERIFY(oldb.isApprox(b) && "sparse solver testing: the rhs should not be modified!"); + Rhs x(b.rows(), b.cols()); + Rhs oldb = b; + + solver.compute(A); + if (solver.info() != Success) + { + std::cerr << "sparse solver testing: factorization failed (check_sparse_solving)\n"; + exit(0); + return; + } + x = solver.solve(b); + if (solver.info() != Success) + { + std::cerr << "sparse solver testing: solving failed\n"; + return; + } + VERIFY(oldb.isApprox(b) && "sparse solver testing: the rhs should not be modified!"); - VERIFY(x.isApprox(refX,test_precision<Scalar>())); - - x.setZero(); - // test the analyze/factorize API - solver.analyzePattern(A); - solver.factorize(A); - if (solver.info() != Success) - { - std::cerr << "sparse solver testing: factorization failed (check_sparse_solving)\n"; - exit(0); - return; - } - x = solver.solve(b); - if (solver.info() != Success) - { - std::cerr << "sparse solver testing: solving failed\n"; - return; - } - VERIFY(oldb.isApprox(b) && "sparse solver testing: the rhs should not be modified!"); + VERIFY(x.isApprox(refX,test_precision<Scalar>())); + x.setZero(); + // test the analyze/factorize API + solver.analyzePattern(A); + solver.factorize(A); + if (solver.info() != Success) + { + std::cerr << "sparse solver testing: factorization failed (check_sparse_solving)\n"; + exit(0); + return; + } + x = solver.solve(b); + if (solver.info() != Success) + { + std::cerr << "sparse solver testing: solving failed\n"; + return; + } + VERIFY(oldb.isApprox(b) && "sparse solver testing: the rhs should not be modified!"); - VERIFY(x.isApprox(refX,test_precision<Scalar>())); + VERIFY(x.isApprox(refX,test_precision<Scalar>())); + } - // test Block as the result and rhs: + // test dense Block as the result and rhs: { DenseRhs x(db.rows(), db.cols()); - DenseRhs b(db), oldb(db); + DenseRhs oldb(db); x.setZero(); - x.block(0,0,x.rows(),x.cols()) = solver.solve(b.block(0,0,b.rows(),b.cols())); - VERIFY(oldb.isApprox(b) && "sparse solver testing: the rhs should not be modified!"); + x.block(0,0,x.rows(),x.cols()) = solver.solve(db.block(0,0,db.rows(),db.cols())); + VERIFY(oldb.isApprox(db) && "sparse solver testing: the rhs should not be modified!"); VERIFY(x.isApprox(refX,test_precision<Scalar>())); } } @@ -113,7 +113,6 @@ void check_sparse_determinant(Solver& solver, const typename Solver::MatrixType& { typedef typename Solver::MatrixType Mat; typedef typename Mat::Scalar Scalar; - typedef typename Mat::RealScalar RealScalar; solver.compute(A); if (solver.info() != Success) @@ -158,9 +157,9 @@ inline std::string get_matrixfolder() { std::string mat_folder = TEST_REAL_CASES; if( internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value ) - mat_folder = mat_folder + static_cast<string>("/complex/"); + mat_folder = mat_folder + static_cast<std::string>("/complex/"); else - mat_folder = mat_folder + static_cast<string>("/real/"); + mat_folder = mat_folder + static_cast<std::string>("/real/"); return mat_folder; } #endif @@ -169,7 +168,6 @@ template<typename Solver> void check_sparse_spd_solving(Solver& solver) { typedef typename Solver::MatrixType Mat; typedef typename Mat::Scalar Scalar; - typedef typename Mat::Index Index; typedef SparseMatrix<Scalar,ColMajor> SpMat; typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; typedef Matrix<Scalar,Dynamic,1> DenseVector; @@ -177,25 +175,32 @@ template<typename Solver> void check_sparse_spd_solving(Solver& solver) // generate the problem Mat A, halfA; DenseMatrix dA; - int size = generate_sparse_spd_problem(solver, A, halfA, dA); - - // generate the right hand sides - int rhsCols = internal::random<int>(1,16); - double density = (std::max)(8./(size*rhsCols), 0.1); - SpMat B(size,rhsCols); - DenseVector b = DenseVector::Random(size); - DenseMatrix dB(size,rhsCols); - initSparse<Scalar>(density, dB, B, ForceNonZeroDiag); - for (int i = 0; i < g_repeat; i++) { + int size = generate_sparse_spd_problem(solver, A, halfA, dA); + + // generate the right hand sides + int rhsCols = internal::random<int>(1,16); + double density = (std::max)(8./(size*rhsCols), 0.1); + SpMat B(size,rhsCols); + DenseVector b = DenseVector::Random(size); + DenseMatrix dB(size,rhsCols); + initSparse<Scalar>(density, dB, B, ForceNonZeroDiag); + check_sparse_solving(solver, A, b, dA, b); check_sparse_solving(solver, halfA, b, dA, b); check_sparse_solving(solver, A, dB, dA, dB); check_sparse_solving(solver, halfA, dB, dA, dB); check_sparse_solving(solver, A, B, dA, dB); check_sparse_solving(solver, halfA, B, dA, dB); + + // check only once + if(i==0) + { + b = DenseVector::Zero(size); + check_sparse_solving(solver, A, b, dA, b); + } } - + // First, get the folder #ifdef TEST_REAL_CASES if (internal::is_same<Scalar, float>::value @@ -241,7 +246,6 @@ int generate_sparse_square_problem(Solver&, typename Solver::MatrixType& A, Dens { typedef typename Solver::MatrixType Mat; typedef typename Mat::Scalar Scalar; - typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; int size = internal::random<int>(1,maxSize); double density = (std::max)(8./(size*size), 0.01); @@ -258,6 +262,7 @@ template<typename Solver> void check_sparse_square_solving(Solver& solver) { typedef typename Solver::MatrixType Mat; typedef typename Mat::Scalar Scalar; + typedef SparseMatrix<Scalar,ColMajor> SpMat; typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; typedef Matrix<Scalar,Dynamic,1> DenseVector; @@ -265,16 +270,28 @@ template<typename Solver> void check_sparse_square_solving(Solver& solver) Mat A; DenseMatrix dA; - int size = generate_sparse_square_problem(solver, A, dA); - - DenseVector b = DenseVector::Random(size); - DenseMatrix dB = DenseMatrix::Random(size,rhsCols); - A.makeCompressed(); for (int i = 0; i < g_repeat; i++) { + int size = generate_sparse_square_problem(solver, A, dA); + + A.makeCompressed(); + DenseVector b = DenseVector::Random(size); + DenseMatrix dB(size,rhsCols); + SpMat B(size,rhsCols); + double density = (std::max)(8./(size*rhsCols), 0.1); + initSparse<Scalar>(density, dB, B, ForceNonZeroDiag); + B.makeCompressed(); check_sparse_solving(solver, A, b, dA, b); check_sparse_solving(solver, A, dB, dA, dB); + check_sparse_solving(solver, A, B, dA, dB); + + // check only once + if(i==0) + { + b = DenseVector::Zero(size); + check_sparse_solving(solver, A, b, dA, b); + } } - + // First, get the folder #ifdef TEST_REAL_CASES if (internal::is_same<Scalar, float>::value diff --git a/test/sparse_vector.cpp b/test/sparse_vector.cpp index 7201afe5b..0c9476803 100644 --- a/test/sparse_vector.cpp +++ b/test/sparse_vector.cpp @@ -9,14 +9,14 @@ #include "sparse.h" -template<typename Scalar> void sparse_vector(int rows, int cols) +template<typename Scalar,typename Index> void sparse_vector(int rows, int cols) { double densityMat = (std::max)(8./(rows*cols), 0.01); double densityVec = (std::max)(8./float(rows), 0.1); typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; typedef Matrix<Scalar,Dynamic,1> DenseVector; - typedef SparseVector<Scalar> SparseVectorType; - typedef SparseMatrix<Scalar> SparseMatrixType; + typedef SparseVector<Scalar,0,Index> SparseVectorType; + typedef SparseMatrix<Scalar,0,Index> SparseMatrixType; Scalar eps = 1e-6; SparseMatrixType m1(rows,rows); @@ -77,15 +77,34 @@ template<typename Scalar> void sparse_vector(int rows, int cols) VERIFY_IS_APPROX(v1.squaredNorm(), refV1.squaredNorm()); + + VERIFY_IS_APPROX(v1.blueNorm(), refV1.blueNorm()); + + // test aliasing + VERIFY_IS_APPROX((v1 = -v1), (refV1 = -refV1)); + VERIFY_IS_APPROX((v1 = v1.transpose()), (refV1 = refV1.transpose().eval())); + VERIFY_IS_APPROX((v1 += -v1), (refV1 += -refV1)); + + // sparse matrix to sparse vector + SparseMatrixType mv1; + VERIFY_IS_APPROX((mv1=v1),v1); + VERIFY_IS_APPROX(mv1,(v1=mv1)); + VERIFY_IS_APPROX(mv1,(v1=mv1.transpose())); + + // check copy to dense vector with transpose + refV3.resize(0); + VERIFY_IS_APPROX(refV3 = v1.transpose(),v1.toDense()); + VERIFY_IS_APPROX(DenseVector(v1),v1.toDense()); } void test_sparse_vector() { for(int i = 0; i < g_repeat; i++) { - CALL_SUBTEST_1( sparse_vector<double>(8, 8) ); - CALL_SUBTEST_2( sparse_vector<std::complex<double> >(16, 16) ); - CALL_SUBTEST_1( sparse_vector<double>(299, 535) ); + CALL_SUBTEST_1(( sparse_vector<double,int>(8, 8) )); + CALL_SUBTEST_2(( sparse_vector<std::complex<double>, int>(16, 16) )); + CALL_SUBTEST_1(( sparse_vector<double,long int>(299, 535) )); + CALL_SUBTEST_1(( sparse_vector<double,short>(299, 535) )); } } diff --git a/test/sparselu.cpp b/test/sparselu.cpp new file mode 100644 index 000000000..37980defc --- /dev/null +++ b/test/sparselu.cpp @@ -0,0 +1,55 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + + +// SparseLU solve does not accept column major matrices for the destination. +// However, as expected, the generic check_sparse_square_solving routines produces row-major +// rhs and destination matrices when compiled with EIGEN_DEFAULT_TO_ROW_MAJOR + +#ifdef EIGEN_DEFAULT_TO_ROW_MAJOR +#undef EIGEN_DEFAULT_TO_ROW_MAJOR +#endif + +#include "sparse_solver.h" +#include <Eigen/SparseLU> +#include <unsupported/Eigen/SparseExtra> + +template<typename T> void test_sparselu_T() +{ + SparseLU<SparseMatrix<T, ColMajor> /*, COLAMDOrdering<int>*/ > sparselu_colamd; // COLAMDOrdering is the default + SparseLU<SparseMatrix<T, ColMajor>, AMDOrdering<int> > sparselu_amd; + SparseLU<SparseMatrix<T, ColMajor, long int>, NaturalOrdering<long int> > sparselu_natural; + + check_sparse_square_solving(sparselu_colamd); + check_sparse_square_solving(sparselu_amd); + check_sparse_square_solving(sparselu_natural); +} + +void test_sparselu() +{ + CALL_SUBTEST_1(test_sparselu_T<float>()); + CALL_SUBTEST_2(test_sparselu_T<double>()); + CALL_SUBTEST_3(test_sparselu_T<std::complex<float> >()); + CALL_SUBTEST_4(test_sparselu_T<std::complex<double> >()); +} diff --git a/test/sparseqr.cpp b/test/sparseqr.cpp new file mode 100644 index 000000000..1fe4a98ee --- /dev/null +++ b/test/sparseqr.cpp @@ -0,0 +1,99 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Desire Nuentsa Wakam <desire.nuentsa_wakam@inria.fr> +// Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr> +// +// 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 +#include "sparse.h" +#include <Eigen/SparseQR> + +template<typename MatrixType,typename DenseMat> +int generate_sparse_rectangular_problem(MatrixType& A, DenseMat& dA, int maxRows = 300, int maxCols = 150) +{ + eigen_assert(maxRows >= maxCols); + typedef typename MatrixType::Scalar Scalar; + int rows = internal::random<int>(1,maxRows); + int cols = internal::random<int>(1,maxCols); + double density = (std::max)(8./(rows*cols), 0.01); + + A.resize(rows,cols); + dA.resize(rows,cols); + initSparse<Scalar>(density, dA, A,ForceNonZeroDiag); + A.makeCompressed(); + int nop = internal::random<int>(0, internal::random<double>(0,1) > 0.5 ? cols/2 : 0); + for(int k=0; k<nop; ++k) + { + int j0 = internal::random<int>(0,cols-1); + int j1 = internal::random<int>(0,cols-1); + Scalar s = internal::random<Scalar>(); + A.col(j0) = s * A.col(j1); + dA.col(j0) = s * dA.col(j1); + } + +// if(rows<cols) { +// A.conservativeResize(cols,cols); +// dA.conservativeResize(cols,cols); +// dA.bottomRows(cols-rows).setZero(); +// } + + return rows; +} + +template<typename Scalar> void test_sparseqr_scalar() +{ + typedef SparseMatrix<Scalar,ColMajor> MatrixType; + typedef Matrix<Scalar,Dynamic,Dynamic> DenseMat; + typedef Matrix<Scalar,Dynamic,1> DenseVector; + MatrixType A; + DenseMat dA; + DenseVector refX,x,b; + SparseQR<MatrixType, COLAMDOrdering<int> > solver; + generate_sparse_rectangular_problem(A,dA); + + b = dA * DenseVector::Random(A.cols()); + solver.compute(A); + if (solver.info() != Success) + { + std::cerr << "sparse QR factorization failed\n"; + exit(0); + return; + } + x = solver.solve(b); + if (solver.info() != Success) + { + std::cerr << "sparse QR factorization failed\n"; + exit(0); + return; + } + + VERIFY_IS_APPROX(A * x, b); + + //Compare with a dense QR solver + ColPivHouseholderQR<DenseMat> dqr(dA); + refX = dqr.solve(b); + + VERIFY_IS_EQUAL(dqr.rank(), solver.rank()); + if(solver.rank()==A.cols()) // full rank + VERIFY_IS_APPROX(x, refX); +// else +// VERIFY((dA * refX - b).norm() * 2 > (A * x - b).norm() ); + + // Compute explicitly the matrix Q + MatrixType Q, QtQ, idM; + Q = solver.matrixQ(); + //Check ||Q' * Q - I || + QtQ = Q * Q.adjoint(); + idM.resize(Q.rows(), Q.rows()); idM.setIdentity(); + VERIFY(idM.isApprox(QtQ)); +} +void test_sparseqr() +{ + for(int i=0; i<g_repeat; ++i) + { + CALL_SUBTEST_1(test_sparseqr_scalar<double>()); + CALL_SUBTEST_2(test_sparseqr_scalar<std::complex<double> >()); + } +} + diff --git a/test/special_numbers.cpp b/test/special_numbers.cpp new file mode 100644 index 000000000..2f1b704be --- /dev/null +++ b/test/special_numbers.cpp @@ -0,0 +1,58 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2013 Gael Guennebaud <gael.guennebaud@inria.fr> +// +// 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/. + +#include "main.h" + +template<typename Scalar> void special_numbers() +{ + typedef Matrix<Scalar, Dynamic,Dynamic> MatType; + int rows = internal::random<int>(1,300); + int cols = internal::random<int>(1,300); + + Scalar nan = std::numeric_limits<Scalar>::quiet_NaN(); + Scalar inf = std::numeric_limits<Scalar>::infinity(); + Scalar s1 = internal::random<Scalar>(); + + MatType m1 = MatType::Random(rows,cols), + mnan = MatType::Random(rows,cols), + minf = MatType::Random(rows,cols), + mboth = MatType::Random(rows,cols); + + int n = internal::random<int>(1,10); + for(int k=0; k<n; ++k) + { + mnan(internal::random<int>(0,rows-1), internal::random<int>(0,cols-1)) = nan; + minf(internal::random<int>(0,rows-1), internal::random<int>(0,cols-1)) = inf; + } + mboth = mnan + minf; + + VERIFY(!m1.hasNaN()); + VERIFY(m1.allFinite()); + + VERIFY(mnan.hasNaN()); + VERIFY((s1*mnan).hasNaN()); + VERIFY(!minf.hasNaN()); + VERIFY(!(2*minf).hasNaN()); + VERIFY(mboth.hasNaN()); + VERIFY(mboth.array().hasNaN()); + + VERIFY(!mnan.allFinite()); + VERIFY(!minf.allFinite()); + VERIFY(!(minf-mboth).allFinite()); + VERIFY(!mboth.allFinite()); + VERIFY(!mboth.array().allFinite()); +} + +void test_special_numbers() +{ + for(int i = 0; i < 10*g_repeat; i++) { + CALL_SUBTEST_1( special_numbers<float>() ); + CALL_SUBTEST_1( special_numbers<double>() ); + } +} diff --git a/test/spqr_support.cpp b/test/spqr_support.cpp new file mode 100644 index 000000000..b8980e081 --- /dev/null +++ b/test/spqr_support.cpp @@ -0,0 +1,62 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Desire Nuentsa Wakam <desire.nuentsa_wakam@inria.fr> +// +// 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 +#include "sparse.h" +#include <Eigen/SPQRSupport> + + +template<typename MatrixType,typename DenseMat> +int generate_sparse_rectangular_problem(MatrixType& A, DenseMat& dA, int maxRows = 300, int maxCols = 300) +{ + eigen_assert(maxRows >= maxCols); + typedef typename MatrixType::Scalar Scalar; + int rows = internal::random<int>(1,maxRows); + int cols = internal::random<int>(1,rows); + double density = (std::max)(8./(rows*cols), 0.01); + + A.resize(rows,rows); + dA.resize(rows,rows); + initSparse<Scalar>(density, dA, A,ForceNonZeroDiag); + A.makeCompressed(); + return rows; +} + +template<typename Scalar> void test_spqr_scalar() +{ + typedef SparseMatrix<Scalar,ColMajor> MatrixType; + MatrixType A; + Matrix<Scalar,Dynamic,Dynamic> dA; + typedef Matrix<Scalar,Dynamic,1> DenseVector; + DenseVector refX,x,b; + SPQR<MatrixType> solver; + generate_sparse_rectangular_problem(A,dA); + + int m = A.rows(); + b = DenseVector::Random(m); + solver.compute(A); + if (solver.info() != Success) + { + std::cerr << "sparse QR factorization failed\n"; + exit(0); + return; + } + x = solver.solve(b); + if (solver.info() != Success) + { + std::cerr << "sparse QR factorization failed\n"; + exit(0); + return; + } + //Compare with a dense solver + refX = dA.colPivHouseholderQr().solve(b); + VERIFY(x.isApprox(refX,test_precision<Scalar>())); +} +void test_spqr_support() +{ + CALL_SUBTEST_1(test_spqr_scalar<double>()); + CALL_SUBTEST_2(test_spqr_scalar<std::complex<double> >()); +} diff --git a/test/stable_norm.cpp b/test/stable_norm.cpp index a25dbf51c..549f91fbf 100644 --- a/test/stable_norm.cpp +++ b/test/stable_norm.cpp @@ -32,6 +32,8 @@ template<typename MatrixType> void stable_norm(const MatrixType& m) /* this test covers the following files: StableNorm.h */ + using std::sqrt; + using std::abs; typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; @@ -53,8 +55,16 @@ template<typename MatrixType> void stable_norm(const MatrixType& m) Index rows = m.rows(); Index cols = m.cols(); - Scalar big = internal::random<Scalar>() * ((std::numeric_limits<RealScalar>::max)() * RealScalar(1e-4)); - Scalar small = internal::random<Scalar>() * ((std::numeric_limits<RealScalar>::min)() * RealScalar(1e4)); + // get a non-zero random factor + Scalar factor = internal::random<Scalar>(); + while(numext::abs2(factor)<RealScalar(1e-4)) + factor = internal::random<Scalar>(); + Scalar big = factor * ((std::numeric_limits<RealScalar>::max)() * RealScalar(1e-4)); + + factor = internal::random<Scalar>(); + while(numext::abs2(factor)<RealScalar(1e-4)) + factor = internal::random<Scalar>(); + Scalar small = factor * ((std::numeric_limits<RealScalar>::min)() * RealScalar(1e4)); MatrixType vzero = MatrixType::Zero(rows, cols), vrand = MatrixType::Random(rows, cols), @@ -73,23 +83,23 @@ template<typename MatrixType> void stable_norm(const MatrixType& m) // test isFinite VERIFY(!isFinite( std::numeric_limits<RealScalar>::infinity())); - VERIFY(!isFinite(internal::sqrt(-internal::abs(big)))); + VERIFY(!isFinite(sqrt(-abs(big)))); // test overflow - VERIFY(isFinite(internal::sqrt(size)*internal::abs(big))); - VERIFY_IS_NOT_APPROX(internal::sqrt(copy(vbig.squaredNorm())), internal::abs(internal::sqrt(size)*big)); // here the default norm must fail - VERIFY_IS_APPROX(vbig.stableNorm(), internal::sqrt(size)*internal::abs(big)); - VERIFY_IS_APPROX(vbig.blueNorm(), internal::sqrt(size)*internal::abs(big)); - VERIFY_IS_APPROX(vbig.hypotNorm(), internal::sqrt(size)*internal::abs(big)); + VERIFY(isFinite(sqrt(size)*abs(big))); + VERIFY_IS_NOT_APPROX(sqrt(copy(vbig.squaredNorm())), abs(sqrt(size)*big)); // here the default norm must fail + VERIFY_IS_APPROX(vbig.stableNorm(), sqrt(size)*abs(big)); + VERIFY_IS_APPROX(vbig.blueNorm(), sqrt(size)*abs(big)); + VERIFY_IS_APPROX(vbig.hypotNorm(), sqrt(size)*abs(big)); // test underflow - VERIFY(isFinite(internal::sqrt(size)*internal::abs(small))); - VERIFY_IS_NOT_APPROX(internal::sqrt(copy(vsmall.squaredNorm())), internal::abs(internal::sqrt(size)*small)); // here the default norm must fail - VERIFY_IS_APPROX(vsmall.stableNorm(), internal::sqrt(size)*internal::abs(small)); - VERIFY_IS_APPROX(vsmall.blueNorm(), internal::sqrt(size)*internal::abs(small)); - VERIFY_IS_APPROX(vsmall.hypotNorm(), internal::sqrt(size)*internal::abs(small)); + VERIFY(isFinite(sqrt(size)*abs(small))); + VERIFY_IS_NOT_APPROX(sqrt(copy(vsmall.squaredNorm())), abs(sqrt(size)*small)); // here the default norm must fail + VERIFY_IS_APPROX(vsmall.stableNorm(), sqrt(size)*abs(small)); + VERIFY_IS_APPROX(vsmall.blueNorm(), sqrt(size)*abs(small)); + VERIFY_IS_APPROX(vsmall.hypotNorm(), sqrt(size)*abs(small)); -// Test compilation of cwise() version + // Test compilation of cwise() version VERIFY_IS_APPROX(vrand.colwise().stableNorm(), vrand.colwise().norm()); VERIFY_IS_APPROX(vrand.colwise().blueNorm(), vrand.colwise().norm()); VERIFY_IS_APPROX(vrand.colwise().hypotNorm(), vrand.colwise().norm()); diff --git a/test/triangular.cpp b/test/triangular.cpp index 0e8ee5487..54320390b 100644 --- a/test/triangular.cpp +++ b/test/triangular.cpp @@ -65,7 +65,7 @@ template<typename MatrixType> void triangular_square(const MatrixType& m) m1 = MatrixType::Random(rows, cols); for (int i=0; i<rows; ++i) - while (internal::abs2(m1(i,i))<1e-1) m1(i,i) = internal::random<Scalar>(); + while (numext::abs2(m1(i,i))<1e-1) m1(i,i) = internal::random<Scalar>(); Transpose<MatrixType> trm4(m4); // test back and forward subsitution with a vector as the rhs @@ -123,9 +123,6 @@ template<typename MatrixType> void triangular_rect(const MatrixType& m) typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime }; - typedef Matrix<Scalar, Rows, 1> VectorType; - typedef Matrix<Scalar, Rows, Rows> RMatrixType; - Index rows = m.rows(); Index cols = m.cols(); @@ -214,8 +211,8 @@ void test_triangular() int maxsize = (std::min)(EIGEN_TEST_MAX_SIZE,20); for(int i = 0; i < g_repeat ; i++) { - int r = internal::random<int>(2,maxsize); EIGEN_UNUSED_VARIABLE(r); - int c = internal::random<int>(2,maxsize); EIGEN_UNUSED_VARIABLE(c); + int r = internal::random<int>(2,maxsize); TEST_SET_BUT_UNUSED_VARIABLE(r) + int c = internal::random<int>(2,maxsize); TEST_SET_BUT_UNUSED_VARIABLE(c) CALL_SUBTEST_1( triangular_square(Matrix<float, 1, 1>()) ); CALL_SUBTEST_2( triangular_square(Matrix<float, 2, 2>()) ); diff --git a/test/umeyama.cpp b/test/umeyama.cpp index b6c9be3a5..2e8092434 100644 --- a/test/umeyama.cpp +++ b/test/umeyama.cpp @@ -22,8 +22,6 @@ template <typename T> Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> randMatrixUnitary(int size) { typedef T Scalar; - typedef typename NumTraits<Scalar>::Real RealScalar; - typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> MatrixType; MatrixType Q; @@ -77,7 +75,6 @@ template <typename T> Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> randMatrixSpecialUnitary(int size) { typedef T Scalar; - typedef typename NumTraits<Scalar>::Real RealScalar; typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> MatrixType; @@ -85,7 +82,7 @@ Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> randMatrixSpecialUnitary(int si MatrixType Q = randMatrixUnitary<Scalar>(size); // tweak the first column to make the determinant be 1 - Q.col(0) *= internal::conj(Q.determinant()); + Q.col(0) *= numext::conj(Q.determinant()); return Q; } @@ -93,13 +90,14 @@ Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> randMatrixSpecialUnitary(int si template <typename MatrixType> void run_test(int dim, int num_elements) { + using std::abs; typedef typename internal::traits<MatrixType>::Scalar Scalar; typedef Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> MatrixX; typedef Matrix<Scalar, Eigen::Dynamic, 1> VectorX; // MUST be positive because in any other case det(cR_t) may become negative for // odd dimensions! - const Scalar c = internal::abs(internal::random<Scalar>()); + const Scalar c = abs(internal::random<Scalar>()); MatrixX R = randMatrixSpecialUnitary<Scalar>(dim); VectorX t = Scalar(50)*VectorX::Random(dim,1); @@ -122,6 +120,7 @@ void run_test(int dim, int num_elements) template<typename Scalar, int Dimension> void run_fixed_size_test(int num_elements) { + using std::abs; typedef Matrix<Scalar, Dimension+1, Dynamic> MatrixX; typedef Matrix<Scalar, Dimension+1, Dimension+1> HomMatrix; typedef Matrix<Scalar, Dimension, Dimension> FixedMatrix; @@ -131,10 +130,11 @@ void run_fixed_size_test(int num_elements) // MUST be positive because in any other case det(cR_t) may become negative for // odd dimensions! - const Scalar c = internal::abs(internal::random<Scalar>()); + // Also if c is to small compared to t.norm(), problem is ill-posed (cf. Bug 744) + const Scalar c = internal::random<Scalar>(0.5, 2.0); FixedMatrix R = randMatrixSpecialUnitary<Scalar>(dim); - FixedVector t = Scalar(50)*FixedVector::Random(dim,1); + FixedVector t = Scalar(32)*FixedVector::Random(dim,1); HomMatrix cR_t = HomMatrix::Identity(dim+1,dim+1); cR_t.block(0,0,dim,dim) = c*R; @@ -150,9 +150,9 @@ void run_fixed_size_test(int num_elements) HomMatrix cR_t_umeyama = umeyama(src_block, dst_block); - const Scalar error = ( cR_t_umeyama*src - dst ).array().square().sum(); + const Scalar error = ( cR_t_umeyama*src - dst ).squaredNorm(); - VERIFY(error < Scalar(10)*std::numeric_limits<Scalar>::epsilon()); + VERIFY(error < Scalar(16)*std::numeric_limits<Scalar>::epsilon()); } void test_umeyama() diff --git a/test/unalignedcount.cpp b/test/unalignedcount.cpp index 5451159e6..ca7e159f3 100644 --- a/test/unalignedcount.cpp +++ b/test/unalignedcount.cpp @@ -40,5 +40,7 @@ void test_unalignedcount() #else // The following line is to eliminate "variable not used" warnings nb_load = nb_loadu = nb_store = nb_storeu = 0; + int a(0), b(0); + VERIFY(a==b); #endif } diff --git a/test/upperbidiagonalization.cpp b/test/upperbidiagonalization.cpp index db6ce383e..d15bf588b 100644 --- a/test/upperbidiagonalization.cpp +++ b/test/upperbidiagonalization.cpp @@ -15,8 +15,8 @@ template<typename MatrixType> void upperbidiag(const MatrixType& m) const typename MatrixType::Index rows = m.rows(); const typename MatrixType::Index cols = m.cols(); - typedef typename MatrixType::Scalar Scalar; typedef Matrix<typename MatrixType::RealScalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime> RealMatrixType; + typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime> TransposeMatrixType; MatrixType a = MatrixType::Random(rows,cols); internal::UpperBidiagonalization<MatrixType> ubd(a); @@ -25,6 +25,8 @@ template<typename MatrixType> void upperbidiag(const MatrixType& m) b.block(0,0,cols,cols) = ubd.bidiagonal(); MatrixType c = ubd.householderU() * b * ubd.householderV().adjoint(); VERIFY_IS_APPROX(a,c); + TransposeMatrixType d = ubd.householderV() * b.adjoint() * ubd.householderU().adjoint(); + VERIFY_IS_APPROX(a.adjoint(),d); } void test_upperbidiagonalization() diff --git a/test/vectorwiseop.cpp b/test/vectorwiseop.cpp index b938e3957..6cd1acdda 100644 --- a/test/vectorwiseop.cpp +++ b/test/vectorwiseop.cpp @@ -15,7 +15,6 @@ template<typename ArrayType> void vectorwiseop_array(const ArrayType& m) { typedef typename ArrayType::Index Index; typedef typename ArrayType::Scalar Scalar; - typedef typename NumTraits<Scalar>::Real RealScalar; typedef Array<Scalar, ArrayType::RowsAtCompileTime, 1> ColVectorType; typedef Array<Scalar, 1, ArrayType::ColsAtCompileTime> RowVectorType; @@ -102,6 +101,16 @@ template<typename ArrayType> void vectorwiseop_array(const ArrayType& m) VERIFY_RAISES_ASSERT(m2.rowwise() /= rowvec.transpose()); VERIFY_RAISES_ASSERT(m1.rowwise() / rowvec.transpose()); + + m2 = m1; + // yes, there might be an aliasing issue there but ".rowwise() /=" + // is suppposed to evaluate " m2.colwise().sum()" into to temporary to avoid + // evaluating the reducions multiple times + if(ArrayType::RowsAtCompileTime>2 || ArrayType::RowsAtCompileTime==Dynamic) + { + m2.rowwise() /= m2.colwise().sum(); + VERIFY_IS_APPROX(m2, m1.rowwise() / m1.colwise().sum()); + } } template<typename MatrixType> void vectorwiseop_matrix(const MatrixType& m) @@ -111,6 +120,8 @@ template<typename MatrixType> void vectorwiseop_matrix(const MatrixType& m) typedef typename NumTraits<Scalar>::Real RealScalar; typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVectorType; typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType; + typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealColVectorType; + typedef Matrix<RealScalar, 1, MatrixType::ColsAtCompileTime> RealRowVectorType; Index rows = m.rows(); Index cols = m.cols(); @@ -123,6 +134,8 @@ template<typename MatrixType> void vectorwiseop_matrix(const MatrixType& m) ColVectorType colvec = ColVectorType::Random(rows); RowVectorType rowvec = RowVectorType::Random(cols); + RealColVectorType rcres; + RealRowVectorType rrres; // test addition @@ -159,6 +172,26 @@ template<typename MatrixType> void vectorwiseop_matrix(const MatrixType& m) VERIFY_RAISES_ASSERT(m2.rowwise() -= rowvec.transpose()); VERIFY_RAISES_ASSERT(m1.rowwise() - rowvec.transpose()); + + // test norm + rrres = m1.colwise().norm(); + VERIFY_IS_APPROX(rrres(c), m1.col(c).norm()); + rcres = m1.rowwise().norm(); + VERIFY_IS_APPROX(rcres(r), m1.row(r).norm()); + + // test normalized + m2 = m1.colwise().normalized(); + VERIFY_IS_APPROX(m2.col(c), m1.col(c).normalized()); + m2 = m1.rowwise().normalized(); + VERIFY_IS_APPROX(m2.row(r), m1.row(r).normalized()); + + // test normalize + m2 = m1; + m2.colwise().normalize(); + VERIFY_IS_APPROX(m2.col(c), m1.col(c).normalized()); + m2 = m1; + m2.rowwise().normalize(); + VERIFY_IS_APPROX(m2.row(r), m1.row(r).normalized()); } void test_vectorwiseop() diff --git a/test/visitor.cpp b/test/visitor.cpp index 087306258..39a5d6b5f 100644 --- a/test/visitor.cpp +++ b/test/visitor.cpp @@ -72,8 +72,8 @@ template<typename VectorType> void vectorVisitor(const VectorType& w) while(v(i) == v(i2)) // yes, == v(i) = internal::random<Scalar>(); - Scalar minc = Scalar(1000), maxc = Scalar(-1000); - Index minidx=0,maxidx=0; + Scalar minc = v(0), maxc = v(0); + Index minidx=0, maxidx=0; for(Index i = 0; i < size; i++) { if(v(i) < minc) @@ -97,6 +97,17 @@ template<typename VectorType> void vectorVisitor(const VectorType& w) VERIFY_IS_APPROX(maxc, eigen_maxc); VERIFY_IS_APPROX(minc, v.minCoeff()); VERIFY_IS_APPROX(maxc, v.maxCoeff()); + + Index idx0 = internal::random<Index>(0,size-1); + Index idx1 = eigen_minidx; + Index idx2 = eigen_maxidx; + VectorType v1(v), v2(v); + v1(idx0) = v1(idx1); + v2(idx0) = v2(idx2); + v1.minCoeff(&eigen_minidx); + v2.maxCoeff(&eigen_maxidx); + VERIFY(eigen_minidx == (std::min)(idx0,idx1)); + VERIFY(eigen_maxidx == (std::min)(idx0,idx2)); } void test_visitor() @@ -111,6 +122,7 @@ void test_visitor() } for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST_7( vectorVisitor(Vector4f()) ); + CALL_SUBTEST_7( vectorVisitor(Matrix<int,12,1>()) ); CALL_SUBTEST_8( vectorVisitor(VectorXd(10)) ); CALL_SUBTEST_9( vectorVisitor(RowVectorXd(10)) ); CALL_SUBTEST_10( vectorVisitor(VectorXf(33)) ); diff --git a/test/zerosized.cpp b/test/zerosized.cpp index 6905e584e..da7dd0481 100644 --- a/test/zerosized.cpp +++ b/test/zerosized.cpp @@ -9,12 +9,26 @@ #include "main.h" + +template<typename MatrixType> void zeroReduction(const MatrixType& m) { + // Reductions that must hold for zero sized objects + VERIFY(m.all()); + VERIFY(!m.any()); + VERIFY(m.prod()==1); + VERIFY(m.sum()==0); + VERIFY(m.count()==0); + VERIFY(m.allFinite()); + VERIFY(!m.hasNaN()); +} + + template<typename MatrixType> void zeroSizedMatrix() { MatrixType t1; - if (MatrixType::SizeAtCompileTime == Dynamic) + if (MatrixType::SizeAtCompileTime == Dynamic || MatrixType::SizeAtCompileTime == 0) { + zeroReduction(t1); if (MatrixType::RowsAtCompileTime == Dynamic) VERIFY(t1.rows() == 0); if (MatrixType::ColsAtCompileTime == Dynamic) @@ -22,9 +36,13 @@ template<typename MatrixType> void zeroSizedMatrix() if (MatrixType::RowsAtCompileTime == Dynamic && MatrixType::ColsAtCompileTime == Dynamic) { + MatrixType t2(0, 0); VERIFY(t2.rows() == 0); VERIFY(t2.cols() == 0); + + zeroReduction(t2); + VERIFY(t1==t2); } } } @@ -33,11 +51,15 @@ template<typename VectorType> void zeroSizedVector() { VectorType t1; - if (VectorType::SizeAtCompileTime == Dynamic) + if (VectorType::SizeAtCompileTime == Dynamic || VectorType::SizeAtCompileTime==0) { + zeroReduction(t1); VERIFY(t1.size() == 0); VectorType t2(DenseIndex(0)); // DenseIndex disambiguates with 0-the-null-pointer (error with gcc 4.4 and MSVC8) VERIFY(t2.size() == 0); + zeroReduction(t2); + + VERIFY(t1==t2); } } @@ -51,9 +73,12 @@ void test_zerosized() zeroSizedMatrix<Matrix<float, Dynamic, 0, 0, 0, 0> >(); zeroSizedMatrix<Matrix<float, 0, Dynamic, 0, 0, 0> >(); zeroSizedMatrix<Matrix<float, Dynamic, Dynamic, 0, 0, 0> >(); - + zeroSizedMatrix<Matrix<float, 0, 4> >(); + zeroSizedMatrix<Matrix<float, 4, 0> >(); + zeroSizedVector<Vector2d>(); zeroSizedVector<Vector3i>(); zeroSizedVector<VectorXf>(); zeroSizedVector<Matrix<float, 0, 1> >(); + zeroSizedVector<Matrix<float, 1, 0> >(); } |