From 7faaa9f3f0df9d23790277834d426c3d992ac3ba Mon Sep 17 00:00:00 2001 From: Carlos Hernandez Date: Tue, 5 Aug 2014 17:53:32 -0700 Subject: Update Eigen to the latest stable release, 3.2.2 ./Eigen/src/Core/util/NonMPL2.h is left untouched, so that usage of non MPL2 code is disabled. Change-Id: I86fc9257b3c30d0ca15b268d4ef07bf038bba7ca --- test/CMakeLists.txt | 51 ++++- test/adjoint.cpp | 83 ++++---- test/array.cpp | 81 ++++---- test/array_for_matrix.cpp | 74 +++++-- test/array_replicate.cpp | 1 - test/basicstuff.cpp | 11 +- test/block.cpp | 39 +++- test/cholesky.cpp | 132 +++++++++++-- test/conservative_resize.cpp | 66 ++++--- test/corners.cpp | 20 ++ test/cwiseop.cpp | 63 +++--- test/denseLM.cpp | 190 ++++++++++++++++++ test/determinant.cpp | 6 +- test/diagonal.cpp | 16 +- test/diagonalmatrices.cpp | 10 +- test/dynalloc.cpp | 3 +- test/eigen2support.cpp | 4 +- test/eigensolver_complex.cpp | 21 +- test/eigensolver_generalized_real.cpp | 59 ++++++ test/eigensolver_generic.cpp | 20 +- test/eigensolver_selfadjoint.cpp | 11 +- test/exceptions.cpp | 4 + test/geo_alignedbox.cpp | 11 +- test/geo_eulerangles.cpp | 100 ++++++++-- test/geo_hyperplane.cpp | 4 +- test/geo_parametrizedline.cpp | 5 +- test/geo_quaternion.cpp | 63 ++++-- test/geo_transformations.cpp | 34 ++-- test/householder.cpp | 29 ++- test/inverse.cpp | 14 +- test/jacobi.cpp | 12 +- test/jacobisvd.cpp | 154 ++++++++++++--- test/linearstructure.cpp | 3 +- test/lu.cpp | 4 - test/main.h | 73 +++++-- test/map.cpp | 144 -------------- test/mapped_matrix.cpp | 141 ++++++++++++++ test/mapstride.cpp | 4 +- test/meta.cpp | 5 +- test/metis_support.cpp | 39 ++++ test/miscmatrices.cpp | 1 - test/nesting_ops.cpp | 12 +- test/nomalloc.cpp | 7 +- test/nullary.cpp | 6 + test/packetmath.cpp | 54 ++++-- test/pastix_support.cpp | 2 +- test/permutationmatrices.cpp | 1 - test/prec_inverse_4x4.cpp | 4 +- test/product.h | 1 - test/product_extra.cpp | 35 +++- test/product_mmtr.cpp | 2 - test/product_notemporary.cpp | 5 +- test/product_selfadjoint.cpp | 9 +- test/product_symm.cpp | 2 - test/product_syrk.cpp | 47 ++++- test/product_trmm.cpp | 5 +- test/product_trmv.cpp | 4 +- test/qr.cpp | 7 +- test/qr_colpivoting.cpp | 8 +- test/qr_fullpivoting.cpp | 11 +- test/real_qz.cpp | 65 +++++++ test/redux.cpp | 39 ++-- test/ref.cpp | 252 ++++++++++++++++++++++++ test/schur_complex.cpp | 19 +- test/schur_real.cpp | 21 +- test/selfadjoint.cpp | 5 +- test/simplicial_cholesky.cpp | 41 ++-- test/sizeof.cpp | 2 +- test/sparse.h | 53 +++-- test/sparseLM.cpp | 176 +++++++++++++++++ test/sparse_basic.cpp | 350 ++++++++++++++++++++++------------ test/sparse_product.cpp | 82 ++++++-- test/sparse_solver.h | 143 ++++++++------ test/sparse_vector.cpp | 31 ++- test/sparselu.cpp | 55 ++++++ test/sparseqr.cpp | 99 ++++++++++ test/special_numbers.cpp | 58 ++++++ test/spqr_support.cpp | 62 ++++++ test/stable_norm.cpp | 38 ++-- test/triangular.cpp | 9 +- test/umeyama.cpp | 18 +- test/unalignedcount.cpp | 2 + test/upperbidiagonalization.cpp | 4 +- test/vectorwiseop.cpp | 35 +++- test/visitor.cpp | 16 +- test/zerosized.cpp | 31 ++- 86 files changed, 2885 insertions(+), 848 deletions(-) create mode 100644 test/denseLM.cpp create mode 100644 test/eigensolver_generalized_real.cpp delete mode 100644 test/map.cpp create mode 100644 test/mapped_matrix.cpp create mode 100644 test/metis_support.cpp create mode 100644 test/real_qz.cpp create mode 100644 test/ref.cpp create mode 100644 test/sparseLM.cpp create mode 100644 test/sparselu.cpp create mode 100644 test/sparseqr.cpp create mode 100644 test/special_numbers.cpp create mode 100644 test/spqr_support.cpp (limited to 'test') 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 struct adjoint_specific; + +template<> struct adjoint_specific { + template + 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 { + template + static void run(const Vec& v1, const Vec& v2, Vec& v3, const Mat& square, Scalar s1, Scalar s2) { + typedef typename NumTraits::Real RealScalar; + using std::abs; + + RealScalar ref = NumTraits::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::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())); + + // 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 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::Real RealScalar; @@ -43,45 +86,21 @@ template 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::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::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::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(1)); + adjoint_specific::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(1)); - ref = NumTraits::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(0, rows-1), c = internal::random(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::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 void array(const ArrayType& m) { typedef typename ArrayType::Index Index; typedef typename ArrayType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; typedef Array ColVectorType; typedef Array RowVectorType; @@ -64,9 +63,12 @@ template 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())) + 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())) VERIFY_IS_NOT_APPROX(((m1+m2).rowwise().sum()).sum(), m1.sum()); VERIFY_IS_APPROX(m1.colwise().sum(), m1.colwise().redux(internal::scalar_sum_op())); @@ -83,10 +85,10 @@ template void array(const ArrayType& m) template void comparisons(const ArrayType& m) { + using std::abs; typedef typename ArrayType::Index Index; typedef typename ArrayType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; - typedef Array VectorType; Index rows = m.rows(); Index cols = m.cols(); @@ -120,7 +122,7 @@ template void comparisons(const ArrayType& m) Scalar mid = (m1.cwiseAbs().minCoeff() + m1.cwiseAbs().maxCoeff())/Scalar(2); for (int j=0; j void comparisons(const ArrayType& m) template 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::Real RealScalar; @@ -157,58 +161,58 @@ template 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(); // 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::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::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::epsilon()); + const RealScalar tiny = sqrt(std::numeric_limits::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 void array_complex(const ArrayType& m) @@ -223,11 +227,10 @@ template 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 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 void array_for_matrix(const MatrixType& m) { typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; typedef Matrix ColVectorType; typedef Matrix RowVectorType; @@ -26,10 +25,10 @@ template void array_for_matrix(const MatrixType& m) ColVectorType cv1 = ColVectorType::Random(rows); RowVectorType rv1 = RowVectorType::Random(cols); - + Scalar s1 = internal::random(), s2 = internal::random(); - + // 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 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())); // vector-wise ops @@ -73,10 +72,10 @@ template void array_for_matrix(const MatrixType& m) template void comparisons(const MatrixType& m) { + using std::abs; typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; - typedef Matrix VectorType; Index rows = m.rows(); Index cols = m.cols(); @@ -110,7 +109,7 @@ template void comparisons(const MatrixType& m) Scalar mid = (m1.cwiseAbs().minCoeff() + m1.cwiseAbs().maxCoeff())/Scalar(2); for (int j=0; j void comparisons(const MatrixType& m) template void lpNorm(const VectorType& v) { + using std::sqrt; VectorType u = VectorType::Random(v.size()); VERIFY_IS_APPROX(u.template lpNorm(), 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 void cwise_min_max(const MatrixType& m) @@ -163,11 +163,53 @@ template 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 void resize(const MatrixTraits& t) +{ + typedef typename MatrixTraits::Index Index; + typedef typename MatrixTraits::Scalar Scalar; + typedef Matrix MatrixType; + typedef Array Array2DType; + typedef Matrix VectorType; + typedef Array 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(3); } void test_array_for_matrix() @@ -202,4 +244,10 @@ void test_array_for_matrix() CALL_SUBTEST_5( lpNorm(VectorXf(internal::random(1,EIGEN_TEST_MAX_SIZE))) ); CALL_SUBTEST_4( lpNorm(VectorXcf(internal::random(1,EIGEN_TEST_MAX_SIZE))) ); } + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_4( resize(MatrixXcf(internal::random(1,EIGEN_TEST_MAX_SIZE), internal::random(1,EIGEN_TEST_MAX_SIZE))) ); + CALL_SUBTEST_5( resize(MatrixXf(internal::random(1,EIGEN_TEST_MAX_SIZE), internal::random(1,EIGEN_TEST_MAX_SIZE))) ); + CALL_SUBTEST_6( resize(MatrixXi(internal::random(1,EIGEN_TEST_MAX_SIZE), internal::random(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 void replicate(const MatrixType& m) */ typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; typedef Matrix VectorType; typedef Matrix MatrixX; typedef Matrix 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 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::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 void basicStuffComplex(const MatrixType& m) Scalar s1 = internal::random(), s2 = internal::random(); - 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::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 Eigen::internal::enable_if::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 Eigen::internal::enable_if::IsComplex,typename MatrixType::Scalar>::type +block_real_only(const MatrixType &, Index, Index, Index, Index, const Scalar&) { + return Scalar(0); +} + + template void block(const MatrixType& m) { typedef typename MatrixType::Index Index; @@ -37,6 +57,8 @@ template void block(const MatrixType& m) Index c1 = internal::random(0,cols-1); Index c2 = internal::random(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 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 b1(1,1); b1(0,0) = m1(r1,c1); @@ -77,6 +100,12 @@ template void block(const MatrixType& m) // check that fixed block() and block() agree Matrix b = m1.template block(3,3); VERIFY_IS_EQUAL(b, m1.block(3,3,BlockRows,BlockCols)); + + // same tests with mixed fixed/dynamic size + m1.template block(1,1,BlockRows,BlockCols) *= s1; + m1.template block(1,1,BlockRows,BlockCols)(0,3) = m1.template block<2,5>(1,1)(1,2); + Matrix b2 = m1.template block(3,3,2,5); + VERIFY_IS_EQUAL(b2, m1.block(3,3,BlockRows,BlockCols)); } if (rows>2) @@ -96,11 +125,11 @@ template 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 void cholesky(const MatrixType& m) symm += a1 * a1.adjoint(); } - SquareMatrixType symmUp = symm.template triangularView(); - SquareMatrixType symmLo = symm.template triangularView(); - // to test if really Cholesky only uses the upper triangular part, uncomment the following // FIXME: currently that fails !! //symm.template part().setZero(); { + SquareMatrixType symmUp = symm.template triangularView(); + SquareMatrixType symmLo = symm.template triangularView(); + LLT chollo(symmLo); VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix()); vecX = chollo.solve(vecB); @@ -114,6 +114,21 @@ template 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().llt().solve(matB); + VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView().llt().solve(matB)); + m2 = m1; + m2 -= symmLo.template selfadjointView().llt().solve(matB); + VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView().llt().solve(matB)); + m2 = m1; + m2.noalias() += symmLo.template selfadjointView().llt().solve(matB); + VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView().llt().solve(matB)); + m2 = m1; + m2.noalias() -= symmLo.template selfadjointView().llt().solve(matB); + VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView().llt().solve(matB)); } // LDLT @@ -165,22 +180,58 @@ template 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().llt().solve(matB); - VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView().llt().solve(matB)); - m2 = m1; - m2 -= symmLo.template selfadjointView().llt().solve(matB); - VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView().llt().solve(matB)); - m2 = m1; - m2.noalias() += symmLo.template selfadjointView().llt().solve(matB); - VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView().llt().solve(matB)); - m2 = m1; - m2.noalias() -= symmLo.template selfadjointView().llt().solve(matB); - VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView().llt().solve(matB)); + // check matrices coming from linear constraints with Lagrange multipliers + if(rows>=3) + { + SquareMatrixType A = symm; + int c = internal::random(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(1,rows-1); + Matrix a = Matrix::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::max_exponent10/8); + Matrix a = Matrix::Random(rows,rows); + Matrix d = Matrix::Random(rows); + for(int k=0; k(-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(symm) )); @@ -262,6 +313,45 @@ template 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 void cholesky_definiteness(const MatrixType& m) +{ + eigen_assert(m.rows() == 2 && m.cols() == 2); + MatrixType mat; + { + mat << 1, 0, 0, -1; + LDLT ldlt(mat); + VERIFY(!ldlt.isNegative()); + VERIFY(!ldlt.isPositive()); + } + { + mat << 1, 2, 2, 1; + LDLT ldlt(mat); + VERIFY(!ldlt.isNegative()); + VERIFY(!ldlt.isPositive()); + } + { + mat << 0, 0, 0, 0; + LDLT ldlt(mat); + VERIFY(ldlt.isNegative()); + VERIFY(ldlt.isPositive()); + } + { + mat << 0, 0, 0, 1; + LDLT ldlt(mat); + VERIFY(!ldlt.isNegative()); + VERIFY(ldlt.isPositive()); + } + { + mat << -1, 0, 0, 0; + LDLT ldlt(mat); + VERIFY(ldlt.isNegative()); + VERIFY(!ldlt.isPositive()); + } +} + template void cholesky_verify_assert() { MatrixType tmp; @@ -284,11 +374,12 @@ template 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()) ); 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(1,EIGEN_TEST_MAX_SIZE); @@ -306,5 +397,6 @@ void test_cholesky() CALL_SUBTEST_9( LLT(10) ); CALL_SUBTEST_9( LDLT(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 void run_vector_tests() { - typedef Matrix MatrixType; + typedef Matrix 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(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(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::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())); - CALL_SUBTEST_1((run_matrix_tests())); - CALL_SUBTEST_2((run_matrix_tests())); - CALL_SUBTEST_2((run_matrix_tests())); - CALL_SUBTEST_3((run_matrix_tests())); - CALL_SUBTEST_3((run_matrix_tests())); - CALL_SUBTEST_4((run_matrix_tests, Eigen::RowMajor>())); - CALL_SUBTEST_4((run_matrix_tests, Eigen::ColMajor>())); - CALL_SUBTEST_5((run_matrix_tests, Eigen::RowMajor>())); - CALL_SUBTEST_6((run_matrix_tests, Eigen::ColMajor>())); - - CALL_SUBTEST_1((run_vector_tests())); - CALL_SUBTEST_2((run_vector_tests())); - CALL_SUBTEST_3((run_vector_tests())); - CALL_SUBTEST_4((run_vector_tests >())); - CALL_SUBTEST_5((run_vector_tests >())); + for(int i=0; i())); + CALL_SUBTEST_1((run_matrix_tests())); + CALL_SUBTEST_2((run_matrix_tests())); + CALL_SUBTEST_2((run_matrix_tests())); + CALL_SUBTEST_3((run_matrix_tests())); + CALL_SUBTEST_3((run_matrix_tests())); + CALL_SUBTEST_4((run_matrix_tests, Eigen::RowMajor>())); + CALL_SUBTEST_4((run_matrix_tests, Eigen::ColMajor>())); + CALL_SUBTEST_5((run_matrix_tests, Eigen::RowMajor>())); + CALL_SUBTEST_6((run_matrix_tests, Eigen::ColMajor>())); + + CALL_SUBTEST_1((run_vector_tests())); + CALL_SUBTEST_2((run_vector_tests())); + CALL_SUBTEST_3((run_vector_tests())); + CALL_SUBTEST_4((run_vector_tests >())); + CALL_SUBTEST_5((run_vector_tests >())); + } } 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 void c VERIFY_IS_EQUAL((matrix.template bottomLeftCorner()), (matrix.template block(rows-r,0))); VERIFY_IS_EQUAL((matrix.template bottomRightCorner()), (matrix.template block(rows-r,cols-c))); + VERIFY_IS_EQUAL((matrix.template topLeftCorner()), (matrix.template topLeftCorner(r,c))); + VERIFY_IS_EQUAL((matrix.template topRightCorner()), (matrix.template topRightCorner(r,c))); + VERIFY_IS_EQUAL((matrix.template bottomLeftCorner()), (matrix.template bottomLeftCorner(r,c))); + VERIFY_IS_EQUAL((matrix.template bottomRightCorner()), (matrix.template bottomRightCorner(r,c))); + + VERIFY_IS_EQUAL((matrix.template topLeftCorner()), (matrix.template topLeftCorner(r,c))); + VERIFY_IS_EQUAL((matrix.template topRightCorner()), (matrix.template topRightCorner(r,c))); + VERIFY_IS_EQUAL((matrix.template bottomLeftCorner()), (matrix.template bottomLeftCorner(r,c))); + VERIFY_IS_EQUAL((matrix.template bottomRightCorner()), (matrix.template bottomRightCorner(r,c))); + VERIFY_IS_EQUAL((matrix.template topRows()), (matrix.template block(0,0))); VERIFY_IS_EQUAL((matrix.template middleRows(sr)), (matrix.template block(sr,0))); VERIFY_IS_EQUAL((matrix.template bottomRows()), (matrix.template block(rows-r,0))); @@ -74,6 +84,16 @@ template void c VERIFY_IS_EQUAL((const_matrix.template bottomLeftCorner()), (const_matrix.template block(rows-r,0))); VERIFY_IS_EQUAL((const_matrix.template bottomRightCorner()), (const_matrix.template block(rows-r,cols-c))); + VERIFY_IS_EQUAL((const_matrix.template topLeftCorner()), (const_matrix.template topLeftCorner(r,c))); + VERIFY_IS_EQUAL((const_matrix.template topRightCorner()), (const_matrix.template topRightCorner(r,c))); + VERIFY_IS_EQUAL((const_matrix.template bottomLeftCorner()), (const_matrix.template bottomLeftCorner(r,c))); + VERIFY_IS_EQUAL((const_matrix.template bottomRightCorner()), (const_matrix.template bottomRightCorner(r,c))); + + VERIFY_IS_EQUAL((const_matrix.template topLeftCorner()), (const_matrix.template topLeftCorner(r,c))); + VERIFY_IS_EQUAL((const_matrix.template topRightCorner()), (const_matrix.template topRightCorner(r,c))); + VERIFY_IS_EQUAL((const_matrix.template bottomLeftCorner()), (const_matrix.template bottomLeftCorner(r,c))); + VERIFY_IS_EQUAL((const_matrix.template bottomRightCorner()), (const_matrix.template bottomRightCorner(r,c))); + VERIFY_IS_EQUAL((const_matrix.template topRows()), (const_matrix.template block(0,0))); VERIFY_IS_EQUAL((const_matrix.template middleRows(sr)), (const_matrix.template block(sr,0))); VERIFY_IS_EQUAL((const_matrix.template bottomRows()), (const_matrix.template block(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 struct AddIfNull { enum { Cost = NumTraits::AddCost }; }; +template +typename Eigen::internal::enable_if::IsInteger,typename MatrixType::Scalar>::type +cwiseops_real_only(MatrixType& m1, MatrixType& m2, MatrixType& m3, MatrixType& mones) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::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 Eigen::internal::enable_if::IsInteger,typename MatrixType::Scalar>::type +cwiseops_real_only(MatrixType& , MatrixType& , MatrixType& , MatrixType& ) +{ + return 0; +} + template void cwiseops(const MatrixType& m) { typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; typedef Matrix 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 void cwiseops(const MatrixType& m) VERIFY_IS_APPROX(m3, m1.cwise() * m2); VERIFY_IS_APPROX(mones, m2.cwise()/m2); - if(!NumTraits::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 void cwiseops(const MatrixType& m) VERIFY( (m1.cwise().max(m2).cwise() > (m1-mones)).all() ); VERIFY( (m1.cwise()(), Scalar(1)))).all() ); - VERIFY( !(m1.cwise()(), Scalar(1)))).all() ); - VERIFY( !(m1.cwise()>m1.unaryExpr(bind2nd(plus(), Scalar(1)))).any() ); + VERIFY( !(m1.cwise()(), Scalar(1)))).all() ); + VERIFY( !(m1.cwise()>m1bis.unaryExpr(bind2nd(plus(), 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 +// Copyright (C) 2012 Gael Guennebaud +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include +#include +#include + +#include "main.h" +#include +using namespace std; +using namespace Eigen; + +template +struct DenseLM : DenseFunctor +{ + typedef DenseFunctor Base; + typedef typename Base::JacobianType JacobianType; + typedef Matrix VectorType; + + DenseLM(int n, int m) : DenseFunctor(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 u(uv, 0, half); + VectorBlock 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 u(uv, 0, half); + VectorBlock 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 u(uv, 0, half); + VectorBlock 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 +int test_minimizeLM(FunctorType& functor, VectorType& uv) +{ + LevenbergMarquardt lm(functor); + LevenbergMarquardtSpace::Status info; + + info = lm.minimize(uv); + + VERIFY_IS_EQUAL(info, 1); + //FIXME Check other parameters + return info; +} + +template +int test_lmder(FunctorType& functor, VectorType& uv) +{ + typedef typename VectorType::Scalar Scalar; + LevenbergMarquardtSpace::Status info; + LevenbergMarquardt lm(functor); + info = lm.lmder1(uv); + + VERIFY_IS_EQUAL(info, 1); + //FIXME Check other parameters + return info; +} + +template +int test_minimizeSteps(FunctorType& functor, VectorType& uv) +{ + LevenbergMarquardtSpace::Status info; + LevenbergMarquardt 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 +void test_denseLM_T() +{ + typedef Matrix VectorType; + + int inputs = 10; + int values = 1000; + DenseLM 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 u(uv, 0, inputs/2); + VectorBlock 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()); + + // CALL_SUBTEST_2(test_sparseLM_T()); +} 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 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 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()) ); CALL_SUBTEST_2( determinant(Matrix()) ); CALL_SUBTEST_3( determinant(Matrix()) ); @@ -62,6 +62,6 @@ void test_determinant() CALL_SUBTEST_5( determinant(Matrix, 10, 10>()) ); s = internal::random(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 void diagonal(const MatrixType& m) { typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; - typedef Matrix VectorType; - typedef Matrix RowVectorType; Index rows = m.rows(); Index cols = m.cols(); @@ -31,17 +28,14 @@ template 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().RowsAtCompileTime!=Dynamic) + if(MatrixType::SizeAtCompileTime!=Dynamic) { VERIFY(m1.template diagonal().RowsAtCompileTime == m1.diagonal(N1).size()); - } - if(m1.template diagonal().RowsAtCompileTime!=Dynamic) - { VERIFY(m1.template diagonal().RowsAtCompileTime == m1.diagonal(N2).size()); } @@ -56,12 +50,12 @@ template void diagonal(const MatrixType& m) VERIFY_IS_APPROX(m2.template diagonal()[0], static_cast(6) * m1.template diagonal()[0]); m2.diagonal(N1) = 2 * m1.diagonal(N1); - VERIFY_IS_APPROX(m2.diagonal(), static_cast(2) * m1.diagonal(N1)); + VERIFY_IS_APPROX(m2.template diagonal(), static_cast(2) * m1.diagonal(N1)); m2.diagonal(N1)[0] *= 3; VERIFY_IS_APPROX(m2.diagonal(N1)[0], static_cast(6) * m1.diagonal(N1)[0]); m2.diagonal(N2) = 2 * m1.diagonal(N2); - VERIFY_IS_APPROX(m2.diagonal(), static_cast(2) * m1.diagonal(N2)); + VERIFY_IS_APPROX(m2.template diagonal(), static_cast(2) * m1.diagonal(N2)); m2.diagonal(N2)[0] *= 3; VERIFY_IS_APPROX(m2.diagonal(N2)[0], static_cast(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 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 VectorType; typedef Matrix RowVectorType; @@ -32,6 +31,8 @@ template void diagonalmatrices(const MatrixType& m) rv2 = RowVectorType::Random(cols); LeftDiagonalMatrix ldm1(v1), ldm2(v2); RightDiagonalMatrix rdm1(rv1), rdm2(rv2); + + Scalar s1 = internal::random(); SquareMatrixType sq_m1 (v1.asDiagonal()); VERIFY_IS_APPROX(sq_m1, v1.asDiagonal().toDenseMatrix()); @@ -76,6 +77,13 @@ template 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 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 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 void eigensolver(const MatrixType& m) typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; - typedef Matrix VectorType; - typedef Matrix RealVectorType; - typedef typename std::complex::Real> Complex; MatrixType a = MatrixType::Random(rows,cols); MatrixType symmA = a.adjoint() * a; @@ -59,6 +56,17 @@ template void eigensolver(const MatrixType& m) // another algorithm so results may differ slightly verify_is_approx_upto_permutation(a.eigenvalues(), ei1.eigenvalues()); + ComplexEigenSolver ei2; + ei2.setMaxIterations(ComplexSchur::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 eiNoEivecs(a, false); VERIFY_IS_EQUAL(eiNoEivecs.info(), Success); VERIFY_IS_APPROX(ei1.eigenvalues(), eiNoEivecs.eigenvalues()); @@ -93,7 +101,7 @@ template 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(1,EIGEN_TEST_MAX_SIZE/4); @@ -101,7 +109,6 @@ void test_eigensolver_complex() CALL_SUBTEST_3( eigensolver(Matrix, 1, 1>()) ); CALL_SUBTEST_4( eigensolver(Matrix3f()) ); } - CALL_SUBTEST_1( eigensolver_verify_assert(Matrix4cf()) ); s = internal::random(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(s)); + CALL_SUBTEST_5(ComplexEigenSolver 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 +// +// 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 +#include + +template 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 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 symmEig(spdA, spdB); + GeneralizedEigenSolver 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(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()) ); + 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 -// Copyright (C) 2010 Jitse Niesen +// Copyright (C) 2010,2012 Jitse Niesen // // 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 void eigensolver(const MatrixType& m) typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; - typedef Matrix VectorType; typedef Matrix RealVectorType; typedef typename std::complex::Real> Complex; @@ -45,6 +44,17 @@ template void eigensolver(const MatrixType& m) VERIFY_IS_APPROX(ei1.eigenvectors().colwise().norm(), RealVectorType::Ones(rows).transpose()); VERIFY_IS_APPROX(a.eigenvalues(), ei1.eigenvalues()); + EigenSolver ei2; + ei2.setMaxIterations(RealSchur::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 eiNoEivecs(a, false); VERIFY_IS_EQUAL(eiNoEivecs.info(), Success); VERIFY_IS_APPROX(ei1.eigenvalues(), eiNoEivecs.eigenvalues()); @@ -78,7 +88,7 @@ template 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(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(s)); + CALL_SUBTEST_5(EigenSolver 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 void selfadjointeigensolver(const MatrixType& m) typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; - typedef Matrix VectorType; - typedef Matrix RealVectorType; - typedef typename std::complex::Real> Complex; RealScalar largerEps = 10*test_precision(); @@ -113,7 +110,7 @@ template 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(1,EIGEN_TEST_MAX_SIZE/4); - CALL_SUBTEST_8(SelfAdjointEigenSolver(s)); - CALL_SUBTEST_8(Tridiagonalization(s)); + CALL_SUBTEST_8(SelfAdjointEigenSolver tmp1(s)); + CALL_SUBTEST_8(Tridiagonalization 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 using namespace std; +template EIGEN_DONT_INLINE +void kill_extra_precision(T& x) { eigen_assert(&x != 0); } + + template void alignedbox(const BoxType& _box) { /* this test covers the following files: @@ -36,6 +40,10 @@ template 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::Real RealScalar; typedef Matrix 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 +// Copyright (C) 2008-2012 Gael Guennebaud // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed @@ -12,10 +12,55 @@ #include #include -template void eulerangles(void) + +template +void verify_euler(const Matrix& ea, int i, int j, int k) { typedef Matrix Matrix3; typedef Matrix Vector3; + typedef AngleAxis 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())) ) + VERIFY((ea-eabis).norm() <= test_precision()); + + // 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 void check_all_var(const Matrix& 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 void eulerangles() +{ + typedef Matrix Matrix3; + typedef Matrix Vector3; + typedef Array Array3; typedef Quaternion Quaternionx; typedef AngleAxis AngleAxisx; @@ -24,25 +69,38 @@ template 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(0,Scalar(M_PI)); + check_all_var(ea); + + ea[0] = ea[1] = internal::random(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 void hyperplane(const HyperplaneType& _plane) const Index dim = _plane.dim(); enum { Options = HyperplaneType::Options }; typedef typename HyperplaneType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; typedef Matrix VectorType; typedef Matrix MatrixType; @@ -79,6 +78,7 @@ template void hyperplane(const HyperplaneType& _plane) template void lines() { + using std::abs; typedef Hyperplane HLine; typedef ParametrizedLine PLine; typedef Matrix Vector; @@ -90,7 +90,7 @@ template void lines() Vector u = Vector::Random(); Vector v = Vector::Random(); Scalar a = internal::random(); - while (internal::abs(a-1) < 1e-4) a = internal::random(); + while (abs(a-1) < 1e-4) a = internal::random(); 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 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::Real RealScalar; typedef Matrix VectorType; - typedef Matrix MatrixType; typedef Hyperplane HyperplaneType; VectorType p0 = VectorType::Random(dim); @@ -35,7 +34,7 @@ template void parametrizedline(const LineType& _line) LineType l0(p0, d0); Scalar s0 = internal::random(); - Scalar s1 = internal::abs(internal::random()); + Scalar s1 = abs(internal::random()); 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 T bounded_acos(T v) template void check_slerp(const QuatType& q0, const QuatType& q1) { + using std::abs; typedef typename QuatType::Scalar Scalar; - typedef Matrix VectorType; typedef AngleAxis AA; Scalar largeEps = test_precision(); 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 void quaternion(void) /* this test covers the following files: Quaternion.h */ - - typedef Matrix Matrix3; + using std::abs; typedef Matrix Vector3; typedef Matrix Vector4; typedef Quaternion Quaternionx; @@ -82,13 +81,13 @@ template 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 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() + if (abs(aa.angle()) > 5*test_precision() && (aa.axis() - v1.normalized()).norm() < 1.99 && (aa.axis() + v1.normalized()).norm() < 1.99) { @@ -157,7 +156,7 @@ template 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 void quaternion(void) template void mapQuaternion(void){ typedef Map, Aligned> MQuaternionA; + typedef Map, Aligned> MCQuaternionA; typedef Map > MQuaternionUA; typedef Map > MCQuaternionUA; typedef Quaternion Quaternionx; + typedef Matrix Vector3; + typedef AngleAxis AngleAxisx; + + Vector3 v0 = Vector3::Random(), + v1 = Vector3::Random(); + Scalar a = internal::random(-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 void mapQuaternion(void){ if(internal::packet_traits::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 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 void non_projective_only() /* this test covers the following files: Cross.h Quaternion.h, Transform.cpp */ - typedef Matrix Matrix2; - typedef Matrix Matrix3; - typedef Matrix Matrix4; - typedef Matrix Vector2; typedef Matrix Vector3; - typedef Matrix Vector4; typedef Quaternion Quaternionx; typedef AngleAxis AngleAxisx; - typedef Transform Transform2; typedef Transform Transform3; - typedef Transform Isometry2; - typedef Transform Isometry3; - typedef typename Transform3::MatrixType MatrixType; - typedef DiagonalMatrix AlignedScaling2; typedef DiagonalMatrix AlignedScaling3; - typedef Translation Translation2; typedef Translation Translation3; Vector3 v0 = Vector3::Random(), @@ -88,7 +77,8 @@ template void transformations() /* this test covers the following files: Cross.h Quaternion.h, Transform.cpp */ - typedef Matrix Matrix2; + using std::cos; + using std::abs; typedef Matrix Matrix3; typedef Matrix Matrix4; typedef Matrix Vector2; @@ -98,10 +88,7 @@ template void transformations() typedef AngleAxis AngleAxisx; typedef Transform Transform2; typedef Transform Transform3; - typedef Transform Isometry2; - typedef Transform Isometry3; typedef typename Transform3::MatrixType MatrixType; - typedef DiagonalMatrix AlignedScaling2; typedef DiagonalMatrix AlignedScaling3; typedef Translation Translation2; typedef Translation Translation3; @@ -115,7 +102,7 @@ template 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 void transformations() // Transform // TODO complete the tests ! a = 0; - while (internal::abs(a)(-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 void transformations() Vector2 v20 = Vector2::Random(); Vector2 v21 = Vector2::Random(); for (int k=0; k<2; ++k) - if (internal::abs(v21[k])(a).toRotationMatrix(); VERIFY_IS_APPROX(t20.fromPositionOrientationScale(v20,a,v21).matrix(), @@ -292,6 +279,13 @@ template 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 void transformations() Rotation2D r2d1d = r2d1.template cast(); VERIFY_IS_APPROX(r2d1d.template cast(),r2d1); - t20 = Translation2(v20) * (Rotation2D(s0) * Scaling(s0)); - t21 = Translation2(v20) * Rotation2D(s0) * Scaling(s0); + t20 = Translation2(v20) * (Rotation2D(s0) * Eigen::Scaling(s0)); + t21 = Translation2(v20) * Rotation2D(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 void householder(const MatrixType& m) typedef Matrix HBlockMatrixType; typedef Matrix HCoeffsVectorType; - typedef Matrix RightSquareMatrixType; - typedef Matrix VBlockMatrixType; typedef Matrix TMatrixType; Matrix _tmp((std::max)(rows,cols)); @@ -62,8 +60,8 @@ template 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 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 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().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 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 void inverse(const MatrixType& m) Index cols = m.cols(); typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; - typedef Matrix VectorType; MatrixType m1(rows, cols), m2(rows, cols), @@ -42,6 +41,9 @@ template 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::Real RealScalar; + typedef Matrix VectorType; + //computeInverseAndDetWithCheck tests //First: an invertible matrix bool invertible; @@ -63,7 +65,7 @@ template 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 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()) ); CALL_SUBTEST_2( inverse(Matrix2d()) ); CALL_SUBTEST_3( inverse(Matrix3f()) ); CALL_SUBTEST_4( inverse(Matrix4f()) ); CALL_SUBTEST_4( inverse(Matrix()) ); - s = internal::random(50,320); + s = internal::random(50,320); CALL_SUBTEST_5( inverse(MatrixXf(s,s)) ); s = internal::random(25,100); CALL_SUBTEST_6( inverse(MatrixXcd(s,s)) ); CALL_SUBTEST_7( inverse(Matrix4d()) ); CALL_SUBTEST_7( inverse(Matrix()) ); } - 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 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 is really important to test as it is the only way to cover conjugation issues in certain unaligned paths CALL_SUBTEST_6(( jacobi(MatrixXcf(r,c)) )); CALL_SUBTEST_6(( jacobi >(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::Real RealScalar; typedef Matrix MatrixUType; typedef Matrix MatrixVType; - typedef Matrix ColVectorType; - typedef Matrix InputVectorType; MatrixType sigma = MatrixType::Zero(rows,cols); sigma.diagonal() = svd.singularValues().template cast(); @@ -70,6 +67,7 @@ template 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(1, cols)); JacobiSVD svd(m, computationOptions); + + if(internal::is_same::value) svd.setThreshold(1e-8); + else if(internal::is_same::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::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::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::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 MatrixType2; + typedef Matrix RhsType2; + typedef Matrix MatrixType2T; + Index rank = RankAtCompileTime2==Dynamic ? internal::random(1,cols) : Index(RankAtCompileTime2); + MatrixType2 m2(rank,cols); + int guard = 0; + do { + m2.setRandom(); + } while(m2.jacobiSvd().setThreshold(test_precision()).rank()!=rank && (++guard)<10); + VERIFY(guard<10); + RhsType2 rhs2 = RhsType2::Random(rank); + // use QR to find a reference minimal norm solution + HouseholderQR qr(m2.adjoint()); + Matrix tmp = qr.matrixQR().topLeftCorner(rank,rank).template triangularView().adjoint().solve(rhs2); + tmp.conservativeResize(cols); + tmp.tail(cols-rank).setZero(); + SolutionType x21 = qr.householderQ() * tmp; + // now check with SVD + JacobiSVD 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 MatrixType3; + typedef Matrix RhsType3; + Index rows3 = RowsAtCompileTime3==Dynamic ? internal::random(rank+1,2*cols) : Index(RowsAtCompileTime3); + Matrix C = Matrix::Random(rows3,rank); + MatrixType3 m3 = C * m2; + RhsType3 rhs3 = C * rhs2; + JacobiSVD 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 @@ -95,27 +174,30 @@ void jacobisvd_test_all_computation_options(const MatrixType& m) if (QRPreconditioner == NoQRPreconditioner && m.rows() != m.cols()) return; JacobiSVD fullSvd(m, ComputeFullU|ComputeFullV); - - jacobisvd_check_full(m, fullSvd); - jacobisvd_solve(m, ComputeFullU | ComputeFullV); - + CALL_SUBTEST(( jacobisvd_check_full(m, fullSvd) )); + CALL_SUBTEST(( jacobisvd_solve(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(m, ComputeFullU | ComputeThinV); - jacobisvd_solve(m, ComputeThinU | ComputeFullV); - jacobisvd_solve(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(m, ComputeFullU | ComputeThinV) )); + CALL_SUBTEST(( jacobisvd_solve(m, ComputeThinU | ComputeFullV) )); + CALL_SUBTEST(( jacobisvd_solve(m, ComputeThinU | ComputeThinV) )); // test reconstruction typedef typename MatrixType::Index Index; @@ -128,12 +210,29 @@ void jacobisvd_test_all_computation_options(const MatrixType& m) template 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::max_exponent10/4; + s = internal::random(1,s); + Matrix d = Matrix::Random(diagSize); + for(Index k=0; k(-s,s)); + m = Matrix::Random(a.rows(),diagSize) * d.asDiagonal() * Matrix::Random(diagSize,a.cols()); + // cancel some coeffs + Index n = internal::random(0,m.size()-1); + for(Index i=0; i(0,m.rows()-1), internal::random(0,m.cols()-1)) = Scalar(0); + } - jacobisvd_test_all_computation_options(m); - jacobisvd_test_all_computation_options(m); - jacobisvd_test_all_computation_options(m); - jacobisvd_test_all_computation_options(m); + CALL_SUBTEST(( jacobisvd_test_all_computation_options(m) )); + CALL_SUBTEST(( jacobisvd_test_all_computation_options(m) )); + CALL_SUBTEST(( jacobisvd_test_all_computation_options(m) )); + CALL_SUBTEST(( jacobisvd_test_all_computation_options(m) )); } template 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 svd; internal::set_is_malloc_allowed(true); svd.compute(m); @@ -323,6 +422,11 @@ void test_jacobisvd() int r = internal::random(1, 30), c = internal::random(1, 30); + + TEST_SET_BUT_UNUSED_VARIABLE(r) + TEST_SET_BUT_UNUSED_VARIABLE(c) + + CALL_SUBTEST_10(( jacobisvd(MatrixXd(r,c)) )); CALL_SUBTEST_7(( jacobisvd(MatrixXf(r,c)) )); CALL_SUBTEST_8(( jacobisvd(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 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 void linearStructure(const MatrixType& m) m3(rows, cols); Scalar s1 = internal::random(); - while (internal::abs(s1)<1e-3) s1 = internal::random(); + while (abs(s1)<1e-3) s1 = internal::random(); Index r = internal::random(0, rows-1), c = internal::random(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 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 void lu_invertible() /* this test covers the following files: LU.h */ - typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; int size = internal::random(1,EIGEN_TEST_MAX_SIZE); @@ -132,8 +130,6 @@ template void lu_partial_piv() PartialPivLU.h */ typedef typename MatrixType::Index Index; - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; Index rows = internal::random(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 #include #include +#include #include #include #include @@ -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 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 g_test_stack; @@ -169,12 +174,17 @@ namespace Eigen #define EIGEN_INTERNAL_DEBUGGING #include // 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(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& m) return m.isUnitary(test_precision::Scalar>()); } +// Forward declaration to avoid ICC warning +template +bool test_is_equal(const T& actual, const U& expected); + template 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 +void createRandomPIMatrixOfRank(typename MatrixType::Index desired_rank, typename MatrixType::Index rows, typename MatrixType::Index cols, MatrixType& m); template 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 +void randomPermutationVector(PermutationVectorType& v, typename PermutationVectorType::Index size); template void randomPermutationVector(PermutationVectorType& v, typename PermutationVectorType::Index size) { @@ -371,20 +392,22 @@ template<> struct GetDifferentType { typedef float type; }; template struct GetDifferentType > { typedef std::complex::type> type; }; -template std::string type_name() { return "other"; } -template<> std::string type_name() { return "float"; } -template<> std::string type_name() { return "double"; } -template<> std::string type_name() { return "int"; } -template<> std::string type_name >() { return "complex"; } +// Forward declaration to avoid ICC warning +template std::string type_name(); +template std::string type_name() { return "other"; } +template<> std::string type_name() { return "float"; } +template<> std::string type_name() { return "double"; } +template<> std::string type_name() { return "int"; } +template<> std::string type_name >() { return "complex"; } template<> std::string type_name >() { return "complex"; } -template<> std::string type_name >() { return "complex"; } +template<> std::string type_name >() { return "complex"; } // 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 + // 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/map.cpp deleted file mode 100644 index fe983e802..000000000 --- a/test/map.cpp +++ /dev/null @@ -1,144 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2006-2010 Benoit Jacob -// -// 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/. - -#ifndef EIGEN_NO_STATIC_ASSERT -#define EIGEN_NO_STATIC_ASSERT // turn static asserts into runtime asserts in order to check them -#endif - -#include "main.h" - -template void map_class_vector(const VectorType& m) -{ - typedef typename VectorType::Index Index; - typedef typename VectorType::Scalar Scalar; - - Index size = m.size(); - - // test Map.h - Scalar* array1 = internal::aligned_new(size); - Scalar* array2 = internal::aligned_new(size); - Scalar* array3 = new Scalar[size+1]; - Scalar* array3unaligned = size_t(array3)%16 == 0 ? array3+1 : array3; - - Map(array1, size) = VectorType::Random(size); - Map(array2, size) = Map(array1, size); - Map(array3unaligned, size) = Map(array1, size); - VectorType ma1 = Map(array1, size); - VectorType ma2 = Map(array2, size); - VectorType ma3 = Map(array3unaligned, size); - VERIFY_IS_EQUAL(ma1, ma2); - VERIFY_IS_EQUAL(ma1, ma3); - #ifdef EIGEN_VECTORIZE - if(internal::packet_traits::Vectorizable) - VERIFY_RAISES_ASSERT((Map(array3unaligned, size))) - #endif - - internal::aligned_delete(array1, size); - internal::aligned_delete(array2, size); - delete[] array3; -} - -template void map_class_matrix(const MatrixType& m) -{ - typedef typename MatrixType::Index Index; - typedef typename MatrixType::Scalar Scalar; - - Index rows = m.rows(), cols = m.cols(), size = rows*cols; - - // test Map.h - Scalar* array1 = internal::aligned_new(size); - for(int i = 0; i < size; i++) array1[i] = Scalar(1); - Scalar* array2 = internal::aligned_new(size); - for(int i = 0; i < size; i++) array2[i] = Scalar(1); - Scalar* array3 = new Scalar[size+1]; - for(int i = 0; i < size+1; i++) array3[i] = Scalar(1); - Scalar* array3unaligned = size_t(array3)%16 == 0 ? array3+1 : array3; - Map(array1, rows, cols) = MatrixType::Ones(rows,cols); - Map(array2, rows, cols) = Map(array1, rows, cols); - Map(array3unaligned, rows, cols) = Map(array1, rows, cols); - MatrixType ma1 = Map(array1, rows, cols); - MatrixType ma2 = Map(array2, rows, cols); - VERIFY_IS_EQUAL(ma1, ma2); - MatrixType ma3 = Map(array3unaligned, rows, cols); - VERIFY_IS_EQUAL(ma1, ma3); - - internal::aligned_delete(array1, size); - internal::aligned_delete(array2, size); - delete[] array3; -} - -template void map_static_methods(const VectorType& m) -{ - typedef typename VectorType::Index Index; - typedef typename VectorType::Scalar Scalar; - - Index size = m.size(); - - // test Map.h - Scalar* array1 = internal::aligned_new(size); - Scalar* array2 = internal::aligned_new(size); - Scalar* array3 = new Scalar[size+1]; - Scalar* array3unaligned = size_t(array3)%16 == 0 ? array3+1 : array3; - - VectorType::MapAligned(array1, size) = VectorType::Random(size); - VectorType::Map(array2, size) = VectorType::Map(array1, size); - VectorType::Map(array3unaligned, size) = VectorType::Map(array1, size); - VectorType ma1 = VectorType::Map(array1, size); - VectorType ma2 = VectorType::MapAligned(array2, size); - VectorType ma3 = VectorType::Map(array3unaligned, size); - VERIFY_IS_EQUAL(ma1, ma2); - VERIFY_IS_EQUAL(ma1, ma3); - - internal::aligned_delete(array1, size); - internal::aligned_delete(array2, size); - delete[] array3; -} - -template 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. - - // verify that map-to-const don't have LvalueBit - typedef typename internal::add_const::type ConstPlainObjectType; - VERIFY( !(internal::traits >::Flags & LvalueBit) ); - VERIFY( !(internal::traits >::Flags & LvalueBit) ); - VERIFY( !(Map::Flags & LvalueBit) ); - VERIFY( !(Map::Flags & LvalueBit) ); -} - -void test_map() -{ - for(int i = 0; i < g_repeat; i++) { - CALL_SUBTEST_1( map_class_vector(Matrix()) ); - CALL_SUBTEST_1( check_const_correctness(Matrix()) ); - CALL_SUBTEST_2( map_class_vector(Vector4d()) ); - CALL_SUBTEST_2( check_const_correctness(Matrix4d()) ); - CALL_SUBTEST_3( map_class_vector(RowVector4f()) ); - CALL_SUBTEST_4( map_class_vector(VectorXcf(8)) ); - CALL_SUBTEST_5( map_class_vector(VectorXi(12)) ); - CALL_SUBTEST_5( check_const_correctness(VectorXi(12)) ); - - CALL_SUBTEST_1( map_class_matrix(Matrix()) ); - CALL_SUBTEST_2( map_class_matrix(Matrix4d()) ); - CALL_SUBTEST_11( map_class_matrix(Matrix()) ); - CALL_SUBTEST_4( map_class_matrix(MatrixXcf(internal::random(1,10),internal::random(1,10))) ); - CALL_SUBTEST_5( map_class_matrix(MatrixXi(internal::random(1,10),internal::random(1,10))) ); - - CALL_SUBTEST_6( map_static_methods(Matrix()) ); - CALL_SUBTEST_7( map_static_methods(Vector3f()) ); - CALL_SUBTEST_8( map_static_methods(RowVector3d()) ); - CALL_SUBTEST_9( map_static_methods(VectorXcd(8)) ); - CALL_SUBTEST_10( map_static_methods(VectorXf(12)) ); - } -} diff --git a/test/mapped_matrix.cpp b/test/mapped_matrix.cpp new file mode 100644 index 000000000..de9dbbde3 --- /dev/null +++ b/test/mapped_matrix.cpp @@ -0,0 +1,141 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2006-2010 Benoit Jacob +// +// 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/. + +#ifndef EIGEN_NO_STATIC_ASSERT +#define EIGEN_NO_STATIC_ASSERT // turn static asserts into runtime asserts in order to check them +#endif + +#include "main.h" + +template void map_class_vector(const VectorType& m) +{ + typedef typename VectorType::Index Index; + typedef typename VectorType::Scalar Scalar; + + Index size = m.size(); + + // test Map.h + Scalar* array1 = internal::aligned_new(size); + Scalar* array2 = internal::aligned_new(size); + Scalar* array3 = new Scalar[size+1]; + Scalar* array3unaligned = size_t(array3)%16 == 0 ? array3+1 : array3; + + Map(array1, size) = VectorType::Random(size); + Map(array2, size) = Map(array1, size); + Map(array3unaligned, size) = Map(array1, size); + VectorType ma1 = Map(array1, size); + VectorType ma2 = Map(array2, size); + VectorType ma3 = Map(array3unaligned, size); + VERIFY_IS_EQUAL(ma1, ma2); + VERIFY_IS_EQUAL(ma1, ma3); + #ifdef EIGEN_VECTORIZE + if(internal::packet_traits::Vectorizable) + VERIFY_RAISES_ASSERT((Map(array3unaligned, size))) + #endif + + internal::aligned_delete(array1, size); + internal::aligned_delete(array2, size); + delete[] array3; +} + +template void map_class_matrix(const MatrixType& m) +{ + typedef typename MatrixType::Index Index; + typedef typename MatrixType::Scalar Scalar; + + Index rows = m.rows(), cols = m.cols(), size = rows*cols; + + // test Map.h + Scalar* array1 = internal::aligned_new(size); + for(int i = 0; i < size; i++) array1[i] = Scalar(1); + Scalar* array2 = internal::aligned_new(size); + for(int i = 0; i < size; i++) array2[i] = Scalar(1); + Scalar* array3 = new Scalar[size+1]; + for(int i = 0; i < size+1; i++) array3[i] = Scalar(1); + Scalar* array3unaligned = size_t(array3)%16 == 0 ? array3+1 : array3; + Map(array1, rows, cols) = MatrixType::Ones(rows,cols); + Map(array2, rows, cols) = Map(array1, rows, cols); + Map(array3unaligned, rows, cols) = Map(array1, rows, cols); + MatrixType ma1 = Map(array1, rows, cols); + MatrixType ma2 = Map(array2, rows, cols); + VERIFY_IS_EQUAL(ma1, ma2); + MatrixType ma3 = Map(array3unaligned, rows, cols); + VERIFY_IS_EQUAL(ma1, ma3); + + internal::aligned_delete(array1, size); + internal::aligned_delete(array2, size); + delete[] array3; +} + +template void map_static_methods(const VectorType& m) +{ + typedef typename VectorType::Index Index; + typedef typename VectorType::Scalar Scalar; + + Index size = m.size(); + + // test Map.h + Scalar* array1 = internal::aligned_new(size); + Scalar* array2 = internal::aligned_new(size); + Scalar* array3 = new Scalar[size+1]; + Scalar* array3unaligned = size_t(array3)%16 == 0 ? array3+1 : array3; + + VectorType::MapAligned(array1, size) = VectorType::Random(size); + VectorType::Map(array2, size) = VectorType::Map(array1, size); + VectorType::Map(array3unaligned, size) = VectorType::Map(array1, size); + VectorType ma1 = VectorType::Map(array1, size); + VectorType ma2 = VectorType::MapAligned(array2, size); + VectorType ma3 = VectorType::Map(array3unaligned, size); + VERIFY_IS_EQUAL(ma1, ma2); + VERIFY_IS_EQUAL(ma1, ma3); + + internal::aligned_delete(array1, size); + internal::aligned_delete(array2, size); + delete[] array3; +} + +template void check_const_correctness(const PlainObjectType&) +{ + // 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. + + // verify that map-to-const don't have LvalueBit + typedef typename internal::add_const::type ConstPlainObjectType; + VERIFY( !(internal::traits >::Flags & LvalueBit) ); + VERIFY( !(internal::traits >::Flags & LvalueBit) ); + VERIFY( !(Map::Flags & LvalueBit) ); + VERIFY( !(Map::Flags & LvalueBit) ); +} + +void test_mapped_matrix() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( map_class_vector(Matrix()) ); + CALL_SUBTEST_1( check_const_correctness(Matrix()) ); + CALL_SUBTEST_2( map_class_vector(Vector4d()) ); + CALL_SUBTEST_2( check_const_correctness(Matrix4d()) ); + CALL_SUBTEST_3( map_class_vector(RowVector4f()) ); + CALL_SUBTEST_4( map_class_vector(VectorXcf(8)) ); + CALL_SUBTEST_5( map_class_vector(VectorXi(12)) ); + CALL_SUBTEST_5( check_const_correctness(VectorXi(12)) ); + + CALL_SUBTEST_1( map_class_matrix(Matrix()) ); + CALL_SUBTEST_2( map_class_matrix(Matrix4d()) ); + CALL_SUBTEST_11( map_class_matrix(Matrix()) ); + CALL_SUBTEST_4( map_class_matrix(MatrixXcf(internal::random(1,10),internal::random(1,10))) ); + CALL_SUBTEST_5( map_class_matrix(MatrixXi(internal::random(1,10),internal::random(1,10))) ); + + CALL_SUBTEST_6( map_static_methods(Matrix()) ); + CALL_SUBTEST_7( map_static_methods(Vector3f()) ); + CALL_SUBTEST_8( map_static_methods(RowVector3d()) ); + CALL_SUBTEST_9( map_static_methods(VectorXcd(8)) ); + CALL_SUBTEST_10( map_static_methods(VectorXf(12)) ); + } +} 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 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(Matrix()) ); CALL_SUBTEST_1( map_class_vector(Matrix()) ); CALL_SUBTEST_2( map_class_vector(Vector4d()) ); @@ -142,5 +142,7 @@ void test_mapstride() CALL_SUBTEST_5( map_class_matrix(MatrixXi(internal::random(1,maxn),internal::random(1,maxn))) ); CALL_SUBTEST_6( map_class_matrix(MatrixXcd(internal::random(1,maxn),internal::random(1,maxn))) ); CALL_SUBTEST_6( map_class_matrix(MatrixXcd(internal::random(1,maxn),internal::random(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::value)); VERIFY((!internal::is_same::value)); @@ -56,7 +53,7 @@ void test_meta() VERIFY(( internal::is_same::type >::value)); VERIFY(internal::meta_sqrt<1>::ret == 1); - #define VERIFY_META_SQRT(X) VERIFY(internal::meta_sqrt::ret == int(internal::sqrt(double(X)))) + #define VERIFY_META_SQRT(X) VERIFY(internal::meta_sqrt::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 +// +// 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 . +#include "sparse_solver.h" +#include +#include +#include + +template void test_metis_T() +{ + SparseLU, MetisOrdering > sparselu_metis; + + check_sparse_square_solving(sparselu_metis); +} + +void test_metis_support() +{ + CALL_SUBTEST_1(test_metis_T()); +} 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 void miscMatrices(const MatrixType& m) typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; typedef Matrix VectorType; - typedef Matrix 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 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 void nomalloc(const MatrixType& m) */ typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; - typedef Matrix 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 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 bool areApprox(const Scalar* a, const Scalar* b, int s { for (int i=0; i >(a,size) << "]" << " != " << Map >(b,size) << "\n"; return false; @@ -99,6 +99,7 @@ struct packet_helper template void packetmath() { + using std::abs; typedef typename internal::packet_traits::type Packet; const int PacketSize = internal::packet_traits::size; typedef typename NumTraits::Real RealScalar; @@ -113,7 +114,7 @@ template void packetmath() { data1[i] = internal::random()/RealScalar(PacketSize); data2[i] = internal::random()/RealScalar(PacketSize); - refvalue = (std::max)(refvalue,internal::abs(data1[i])); + refvalue = (std::max)(refvalue,abs(data1[i])); } internal::pstore(data2, internal::pload(data1)); @@ -144,7 +145,6 @@ template void packetmath() for (int i=0; i Vector; VERIFY(areApprox(ref, data2, PacketSize) && "internal::palign"); } @@ -156,7 +156,7 @@ template 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 void packetmath() template void packetmath_real() { + using std::abs; typedef typename internal::packet_traits::type Packet; const int PacketSize = internal::packet_traits::size; @@ -217,35 +218,50 @@ template void packetmath_real() for (int i=0; i(-1e3,1e3); - data2[i] = internal::random(-1e3,1e3); + data1[i] = internal::random(-1,1) * std::pow(Scalar(10), internal::random(-3,3)); + data2[i] = internal::random(-1,1) * std::pow(Scalar(10), internal::random(-3,3)); } - CHECK_CWISE1_IF(internal::packet_traits::HasSin, internal::sin, internal::psin); - CHECK_CWISE1_IF(internal::packet_traits::HasCos, internal::cos, internal::pcos); - CHECK_CWISE1_IF(internal::packet_traits::HasTan, internal::tan, internal::ptan); + CHECK_CWISE1_IF(internal::packet_traits::HasSin, std::sin, internal::psin); + CHECK_CWISE1_IF(internal::packet_traits::HasCos, std::cos, internal::pcos); + CHECK_CWISE1_IF(internal::packet_traits::HasTan, std::tan, internal::ptan); for (int i=0; i(-1,1); data2[i] = internal::random(-1,1); } - CHECK_CWISE1_IF(internal::packet_traits::HasASin, internal::asin, internal::pasin); - CHECK_CWISE1_IF(internal::packet_traits::HasACos, internal::acos, internal::pacos); + CHECK_CWISE1_IF(internal::packet_traits::HasASin, std::asin, internal::pasin); + CHECK_CWISE1_IF(internal::packet_traits::HasACos, std::acos, internal::pacos); for (int i=0; i(-87,88); data2[i] = internal::random(-87,88); } - CHECK_CWISE1_IF(internal::packet_traits::HasExp, internal::exp, internal::pexp); + CHECK_CWISE1_IF(internal::packet_traits::HasExp, std::exp, internal::pexp); for (int i=0; i(0,1e6); - data2[i] = internal::random(0,1e6); + data1[i] = internal::random(0,1) * std::pow(Scalar(10), internal::random(-6,6)); + data2[i] = internal::random(0,1) * std::pow(Scalar(10), internal::random(-6,6)); } - CHECK_CWISE1_IF(internal::packet_traits::HasLog, internal::log, internal::plog); - CHECK_CWISE1_IF(internal::packet_traits::HasSqrt, internal::sqrt, internal::psqrt); + if(internal::random(0,1)<0.1) + data1[internal::random(0, PacketSize)] = 0; + CHECK_CWISE1_IF(internal::packet_traits::HasLog, std::log, internal::plog); + CHECK_CWISE1_IF(internal::packet_traits::HasSqrt, std::sqrt, internal::psqrt); +} + +template void packetmath_notcomplex() +{ + using std::abs; + typedef typename internal::packet_traits::type Packet; + const int PacketSize = internal::packet_traits::size; + + EIGEN_ALIGN16 Scalar data1[internal::packet_traits::size*4]; + EIGEN_ALIGN16 Scalar data2[internal::packet_traits::size*4]; + EIGEN_ALIGN16 Scalar ref[internal::packet_traits::size*4]; + + Array::Map(data1, internal::packet_traits::size*4).setRandom(); ref[0] = data1[0]; for (int i=0; i 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 >() ); CALL_SUBTEST_2( packetmath >() ); + CALL_SUBTEST_1( packetmath_notcomplex() ); + CALL_SUBTEST_2( packetmath_notcomplex() ); + CALL_SUBTEST_3( packetmath_notcomplex() ); + CALL_SUBTEST_1( packetmath_real() ); CALL_SUBTEST_2( packetmath_real() ); 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()); CALL_SUBTEST_3( (test_pastix_T_LU >()) ); CALL_SUBTEST_4(test_pastix_T_LU >()); -} \ 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 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 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 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 void inverse_permutation_4x4() template 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 void inverse_general_4x4(int repeat) RealScalar absdet; do { m = MatrixType::Random(); - absdet = internal::abs(m.determinant()); + absdet = abs(m.determinant()); } while(absdet < NumTraits::epsilon()); MatrixType inv = m.inverse(); double error = double( (m*inv-MatrixType::Identity()).norm() * absdet / NumTraits::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 void product(const MatrixType& m) */ typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits::NonInteger NonInteger; typedef Matrix RowVectorType; typedef Matrix ColVectorType; typedef Matrix 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 void product_extra(const MatrixType& m) { typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits::NonInteger NonInteger; typedef Matrix RowVectorType; typedef Matrix ColVectorType; typedef Matrix 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(1,EIGEN_TEST_MAX_SIZE/2), internal::random(1,EIGEN_TEST_MAX_SIZE/2))) ); CALL_SUBTEST_4( product_extra(MatrixXcd(internal::random(1,EIGEN_TEST_MAX_SIZE/2), internal::random(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 void mmtr(int size) { - typedef typename NumTraits::Real RealScalar; - typedef Matrix MatrixColMaj; typedef Matrix 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 void product_notemporary(const MatrixType& m) VERIFY_EVALUATION_COUNT( m3.noalias() -= (s1 * m1).template triangularView() * m2, 0); VERIFY_EVALUATION_COUNT( rm3.noalias() = (s1 * m1.adjoint()).template triangularView() * (m2+m2), 1); VERIFY_EVALUATION_COUNT( rm3.noalias() = (s1 * m1.adjoint()).template triangularView() * m2.adjoint(), 0); + + VERIFY_EVALUATION_COUNT( m3.template triangularView() = (m1 * m2.adjoint()), 0); + VERIFY_EVALUATION_COUNT( m3.template triangularView() -= (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() * (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 void product_selfadjoint(const MatrixType& m) { typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; typedef Matrix VectorType; typedef Matrix RowVectorType; @@ -45,11 +44,11 @@ template void product_selfadjoint(const MatrixType& m) m2 = m1.template triangularView(); m2.template selfadjointView().rankUpdate(-v1,s2*v2,s3); - VERIFY_IS_APPROX(m2, (m1 + (s3*(-v1)*(s2*v2).adjoint()+internal::conj(s3)*(s2*v2)*(-v1).adjoint())).template triangularView().toDenseMatrix()); + VERIFY_IS_APPROX(m2, (m1 + (s3*(-v1)*(s2*v2).adjoint()+numext::conj(s3)*(s2*v2)*(-v1).adjoint())).template triangularView().toDenseMatrix()); m2 = m1.template triangularView(); m2.template selfadjointView().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().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().toDenseMatrix()); if (rows>1) { @@ -63,7 +62,7 @@ template 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()) ); CALL_SUBTEST_2( product_selfadjoint(Matrix()) ); @@ -77,5 +76,5 @@ void test_product_selfadjoint() s = internal::random(1,EIGEN_TEST_MAX_SIZE); CALL_SUBTEST_7( product_selfadjoint(Matrix(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 void symm(int size = Size, int othersize = OtherSize) { - typedef typename NumTraits::Real RealScalar; - typedef Matrix MatrixType; typedef Matrix Rhs1; typedef Matrix 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 void syrk(const MatrixType& m) { typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; + typedef Matrix RMatrixType; typedef Matrix Rhs1; typedef Matrix Rhs2; typedef Matrix Rhs3; @@ -22,10 +22,12 @@ template 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(1,320), cols); - Rhs2 rhs2 = Rhs2::Random(rows, internal::random(1,320)); + Rhs1 rhs1 = Rhs1::Random(internal::random(1,320), cols); Rhs1 rhs11 = Rhs1::Random(rhs1.rows(), cols); + Rhs2 rhs2 = Rhs2::Random(rows, internal::random(1,320)); Rhs2 rhs22 = Rhs2::Random(rows, rhs2.cols()); Rhs3 rhs3 = Rhs3::Random(internal::random(1,320), rows); Scalar s1 = internal::random(); @@ -35,19 +37,34 @@ template void syrk(const MatrixType& m) m2.setZero(); VERIFY_IS_APPROX((m2.template selfadjointView().rankUpdate(rhs2,s1)._expression()), ((s1 * rhs2 * rhs2.adjoint()).eval().template triangularView().toDenseMatrix())); + m2.setZero(); + VERIFY_IS_APPROX(((m2.template triangularView() += s1 * rhs2 * rhs22.adjoint()).nestedExpression()), + ((s1 * rhs2 * rhs22.adjoint()).eval().template triangularView().toDenseMatrix())); + m2.setZero(); VERIFY_IS_APPROX(m2.template selfadjointView().rankUpdate(rhs2,s1)._expression(), (s1 * rhs2 * rhs2.adjoint()).eval().template triangularView().toDenseMatrix()); + m2.setZero(); + VERIFY_IS_APPROX((m2.template triangularView() += s1 * rhs22 * rhs2.adjoint()).nestedExpression(), + (s1 * rhs22 * rhs2.adjoint()).eval().template triangularView().toDenseMatrix()); + m2.setZero(); VERIFY_IS_APPROX(m2.template selfadjointView().rankUpdate(rhs1.adjoint(),s1)._expression(), (s1 * rhs1.adjoint() * rhs1).eval().template triangularView().toDenseMatrix()); - + m2.setZero(); + VERIFY_IS_APPROX((m2.template triangularView() += s1 * rhs11.adjoint() * rhs1).nestedExpression(), + (s1 * rhs11.adjoint() * rhs1).eval().template triangularView().toDenseMatrix()); + + m2.setZero(); VERIFY_IS_APPROX(m2.template selfadjointView().rankUpdate(rhs1.adjoint(),s1)._expression(), (s1 * rhs1.adjoint() * rhs1).eval().template triangularView().toDenseMatrix()); + VERIFY_IS_APPROX((m2.template triangularView() = s1 * rhs1.adjoint() * rhs11).nestedExpression(), + (s1 * rhs1.adjoint() * rhs11).eval().template triangularView().toDenseMatrix()); + m2.setZero(); VERIFY_IS_APPROX(m2.template selfadjointView().rankUpdate(rhs3.adjoint(),s1)._expression(), (s1 * rhs3.adjoint() * rhs3).eval().template triangularView().toDenseMatrix()); @@ -63,6 +80,15 @@ template void syrk(const MatrixType& m) m2.setZero(); VERIFY_IS_APPROX((m2.template selfadjointView().rankUpdate(m1.col(c),s1)._expression()), ((s1 * m1.col(c) * m1.col(c).adjoint()).eval().template triangularView().toDenseMatrix())); + rm2.setZero(); + VERIFY_IS_APPROX((rm2.template selfadjointView().rankUpdate(m1.col(c),s1)._expression()), + ((s1 * m1.col(c) * m1.col(c).adjoint()).eval().template triangularView().toDenseMatrix())); + m2.setZero(); + VERIFY_IS_APPROX((m2.template triangularView() += s1 * m3.col(c) * m1.col(c).adjoint()).nestedExpression(), + ((s1 * m3.col(c) * m1.col(c).adjoint()).eval().template triangularView().toDenseMatrix())); + rm2.setZero(); + VERIFY_IS_APPROX((rm2.template triangularView() += s1 * m1.col(c) * m3.col(c).adjoint()).nestedExpression(), + ((s1 * m1.col(c) * m3.col(c).adjoint()).eval().template triangularView().toDenseMatrix())); m2.setZero(); VERIFY_IS_APPROX((m2.template selfadjointView().rankUpdate(m1.col(c).conjugate(),s1)._expression()), @@ -72,9 +98,20 @@ template void syrk(const MatrixType& m) VERIFY_IS_APPROX((m2.template selfadjointView().rankUpdate(m1.col(c).conjugate(),s1)._expression()), ((s1 * m1.col(c).conjugate() * m1.col(c).conjugate().adjoint()).eval().template triangularView().toDenseMatrix())); + m2.setZero(); VERIFY_IS_APPROX((m2.template selfadjointView().rankUpdate(m1.row(c),s1)._expression()), ((s1 * m1.row(c).transpose() * m1.row(c).transpose().adjoint()).eval().template triangularView().toDenseMatrix())); + rm2.setZero(); + VERIFY_IS_APPROX((rm2.template selfadjointView().rankUpdate(m1.row(c),s1)._expression()), + ((s1 * m1.row(c).transpose() * m1.row(c).transpose().adjoint()).eval().template triangularView().toDenseMatrix())); + m2.setZero(); + VERIFY_IS_APPROX((m2.template triangularView() += s1 * m3.row(c).transpose() * m1.row(c).transpose().adjoint()).nestedExpression(), + ((s1 * m3.row(c).transpose() * m1.row(c).transpose().adjoint()).eval().template triangularView().toDenseMatrix())); + rm2.setZero(); + VERIFY_IS_APPROX((rm2.template triangularView() += s1 * m3.row(c).transpose() * m1.row(c).transpose().adjoint()).nestedExpression(), + ((s1 * m3.row(c).transpose() * m1.row(c).transpose().adjoint()).eval().template triangularView().toDenseMatrix())); + m2.setZero(); VERIFY_IS_APPROX((m2.template selfadjointView().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(1,EIGEN_TEST_MAX_SIZE), int cols=internal::random(1,EIGEN_TEST_MAX_SIZE), int otherCols = OtherCols==Dynamic?internal::random(1,EIGEN_TEST_MAX_SIZE):OtherCols) { - typedef typename NumTraits::Real RealScalar; - typedef Matrix TriMatrix; typedef Matrix OnTheRight; typedef Matrix OnTheLeft; @@ -53,10 +51,11 @@ void trmm(int rows=internal::random(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() * (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()).eval()); - VERIFY_IS_APPROX( ge_xs = (s1*mat).adjoint().template triangularView() * ge_left.adjoint(), internal::conj(s1) * triTr.conjugate() * ge_left.adjoint()); + VERIFY_IS_APPROX( ge_xs = (s1*mat).adjoint().template triangularView() * 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 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()) ); CALL_SUBTEST_2( trmv(Matrix()) ); @@ -85,5 +85,5 @@ void test_product_trmv() s = internal::random(1,EIGEN_TEST_MAX_SIZE); CALL_SUBTEST_6( trmv(Matrix(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 void qr(const MatrixType& m) typedef typename MatrixType::Scalar Scalar; typedef Matrix MatrixQType; - typedef Matrix VectorType; MatrixType a = MatrixType::Random(rows,cols); HouseholderQR qrOfA(a); @@ -53,6 +52,8 @@ template void qr_fixedsize() template void qr_invertible() { + using std::log; + using std::abs; typedef typename NumTraits::Real RealScalar; typedef typename MatrixType::Scalar Scalar; @@ -76,12 +77,12 @@ template void qr_invertible() // now construct a matrix with prescribed determinant m1.setZero(); for(int i = 0; i < size; i++) m1(i,i) = internal::random(); - 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 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 void qr() Index rank = internal::random(1, (std::min)(rows, cols)-1); typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; typedef Matrix MatrixQType; - typedef Matrix VectorType; MatrixType m1; createRandomPIMatrixOfRank(rank,rows,cols,m1); ColPivHouseholderQR qr(m1); @@ -72,6 +70,8 @@ template void qr_fixedsize() template void qr_invertible() { + using std::log; + using std::abs; typedef typename NumTraits::Real RealScalar; typedef typename MatrixType::Scalar Scalar; @@ -95,12 +95,12 @@ template void qr_invertible() // now construct a matrix with prescribed determinant m1.setZero(); for(int i = 0; i < size; i++) m1(i,i) = internal::random(); - 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 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 void qr() typedef typename MatrixType::Scalar Scalar; typedef Matrix MatrixQType; - typedef Matrix VectorType; MatrixType m1; createRandomPIMatrixOfRank(rank,rows,cols,m1); FullPivHouseholderQR qr(m1); @@ -51,6 +50,8 @@ template void qr() template void qr_invertible() { + using std::log; + using std::abs; typedef typename NumTraits::Real RealScalar; typedef typename MatrixType::Scalar Scalar; @@ -78,12 +79,12 @@ template void qr_invertible() // now construct a matrix with prescribed determinant m1.setZero(); for(int i = 0; i < size; i++) m1(i,i) = internal::random(); - 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 void qr_verify_assert() @@ -129,4 +130,8 @@ void test_qr_fullpivoting() // Test problem size constructors CALL_SUBTEST_7(FullPivHouseholderQR(10, 20)); + CALL_SUBTEST_7((FullPivHouseholderQR >(10,20))); + CALL_SUBTEST_7((FullPivHouseholderQR >(Matrix::Random()))); + CALL_SUBTEST_7((FullPivHouseholderQR >(20,10))); + CALL_SUBTEST_7((FullPivHouseholderQR >(Matrix::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 +// +// 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 +#include + +template 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 qz(A,B); + + VERIFY_IS_EQUAL(qz.info(), Success); + // check for zeros + bool all_zeros = true; + for (Index i=0; i0 && 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(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()) ); + 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 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(0,rows-1); @@ -61,6 +61,7 @@ template void matrixRedux(const MatrixType& m) template void vectorRedux(const VectorType& w) { + using std::abs; typedef typename VectorType::Index Index; typedef typename VectorType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; @@ -72,15 +73,15 @@ template 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 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 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()) ); CALL_SUBTEST_1( matrixRedux(Array()) ); 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 +// +// 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 void ref_matrix(const MatrixType& m) +{ + typedef typename MatrixType::Index Index; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef Matrix DynMatrixType; + typedef Matrix RealDynMatrixType; + + typedef Ref RefMat; + typedef Ref RefDynMat; + typedef Ref ConstRefDynMat; + typedef Ref > RefRealMatWithStride; + + Index rows = m.rows(), cols = m.cols(); + + MatrixType m1 = MatrixType::Random(rows, cols), + m2 = m1; + + Index i = internal::random(0,rows-1); + Index j = internal::random(0,cols-1); + Index brows = internal::random(1,rows-i); + Index bcols = internal::random(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 void ref_vector(const VectorType& m) +{ + typedef typename VectorType::Index Index; + typedef typename VectorType::Scalar Scalar; + typedef typename VectorType::RealScalar RealScalar; + typedef Matrix DynMatrixType; + typedef Matrix MatrixType; + typedef Matrix RealDynMatrixType; + + typedef Ref RefMat; + typedef Ref RefDynMat; + typedef Ref ConstRefDynMat; + typedef Ref > RefRealMatWithStride; + typedef Ref > 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(0,size-1); + Index bsize = internal::random(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 void check_const_correctness(const PlainObjectType&) +{ + // verify that ref-to-const don't have LvalueBit + typedef typename internal::add_const::type ConstPlainObjectType; + VERIFY( !(internal::traits >::Flags & LvalueBit) ); + VERIFY( !(internal::traits >::Flags & LvalueBit) ); + VERIFY( !(Ref::Flags & LvalueBit) ); + VERIFY( !(Ref::Flags & LvalueBit) ); +} + +template +EIGEN_DONT_INLINE void call_ref_1(Ref a, const B &b) { VERIFY_IS_EQUAL(a,b); } +template +EIGEN_DONT_INLINE void call_ref_2(const Ref& a, const B &b) { VERIFY_IS_EQUAL(a,b); } +template +EIGEN_DONT_INLINE void call_ref_3(Ref > a, const B &b) { VERIFY_IS_EQUAL(a,b); } +template +EIGEN_DONT_INLINE void call_ref_4(const Ref >& a, const B &b) { VERIFY_IS_EQUAL(a,b); } +template +EIGEN_DONT_INLINE void call_ref_5(Ref > a, const B &b) { VERIFY_IS_EQUAL(a,b); } +template +EIGEN_DONT_INLINE void call_ref_6(const Ref >& a, const B &b) { VERIFY_IS_EQUAL(a,b); } +template +EIGEN_DONT_INLINE void call_ref_7(Ref > 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 ab(a,0,3); + const VectorBlock 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()) ); + CALL_SUBTEST_1( check_const_correctness(Matrix()) ); + 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()) ); + CALL_SUBTEST_2( ref_matrix(Matrix4d()) ); + CALL_SUBTEST_1( ref_matrix(Matrix()) ); + CALL_SUBTEST_4( ref_matrix(MatrixXcf(internal::random(1,10),internal::random(1,10))) ); + CALL_SUBTEST_4( ref_matrix(Matrix,10,15>()) ); + CALL_SUBTEST_5( ref_matrix(MatrixXi(internal::random(1,10),internal::random(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 +// Copyright (C) 2010,2012 Jitse Niesen // // 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 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 cs3; + cs3.setMaxIterations(ComplexSchur::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().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()); + VERIFY_IS_EQUAL(cs3.matrixU(), ComplexMatrixType::Identity(size, size)); + // Test computation of only T, not U ComplexSchur 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 +// Copyright (C) 2010,2012 Jitse Niesen // // 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 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 rs3; + rs3.setMaxIterations(RealSchur::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().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 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 void selfadjoint(const MatrixType& m) { typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits::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(1,EIGEN_TEST_MAX_SIZE); EIGEN_UNUSED_VARIABLE(s); + int s = internal::random(1,EIGEN_TEST_MAX_SIZE); CALL_SUBTEST_1( selfadjoint(Matrix()) ); CALL_SUBTEST_2( selfadjoint(Matrix()) ); @@ -55,7 +54,7 @@ void test_selfadjoint() CALL_SUBTEST_4( selfadjoint(MatrixXcd(s,s)) ); CALL_SUBTEST_5( selfadjoint(Matrix(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 void test_simplicial_cholesky_T() { - SimplicialCholesky, Lower> chol_colmajor_lower; - SimplicialCholesky, Upper> chol_colmajor_upper; - SimplicialLLT, Lower> llt_colmajor_lower; - SimplicialLDLT, Upper> llt_colmajor_upper; - SimplicialLDLT, Lower> ldlt_colmajor_lower; - SimplicialLDLT, Upper> ldlt_colmajor_upper; + SimplicialCholesky, Lower> chol_colmajor_lower_amd; + SimplicialCholesky, Upper> chol_colmajor_upper_amd; + SimplicialLLT, Lower> llt_colmajor_lower_amd; + SimplicialLLT, Upper> llt_colmajor_upper_amd; + SimplicialLDLT, Lower> ldlt_colmajor_lower_amd; + SimplicialLDLT, Upper> ldlt_colmajor_upper_amd; + SimplicialLDLT, Lower, NaturalOrdering > ldlt_colmajor_lower_nat; + SimplicialLDLT, Upper, NaturalOrdering > 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 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& refMat, SparseMatrix& sparseMat, int flags = 0, - std::vector* zeroCoords = 0, - std::vector* nonzeroCoords = 0) + std::vector >* zeroCoords = 0, + std::vector >* nonzeroCoords = 0) { enum { IsRowMajor = SparseMatrix::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; jpush_back(Vector2i(ai,aj)); + nonzeroCoords->push_back(Matrix (ai,aj)); } else if (zeroCoords) { - zeroCoords->push_back(Vector2i(ai,aj)); + zeroCoords->push_back(Matrix (ai,aj)); } refMat(ai,aj) = v; } @@ -110,8 +110,8 @@ initSparse(double density, Matrix& refMat, DynamicSparseMatrix& sparseMat, int flags = 0, - std::vector* zeroCoords = 0, - std::vector* nonzeroCoords = 0) + std::vector >* zeroCoords = 0, + std::vector >* nonzeroCoords = 0) { enum { IsRowMajor = DynamicSparseMatrix::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 (ai,aj)); } else if (zeroCoords) { - zeroCoords->push_back(Vector2i(ai,aj)); + zeroCoords->push_back(Matrix (ai,aj)); } refMat(ai,aj) = v; } @@ -154,10 +154,34 @@ initSparse(double density, sparseMat.finalize(); } -template void +template void initSparse(double density, Matrix& refVec, - SparseVector& sparseVec, + SparseVector& sparseVec, + std::vector* zeroCoords = 0, + std::vector* nonzeroCoords = 0) +{ + sparseVec.reserve(int(refVec.size()*density)); + sparseVec.setZero(); + for(Index i=0; i(0,1) < density) ? internal::random() : 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 void +initSparse(double density, + Matrix& refVec, + SparseVector& sparseVec, std::vector* zeroCoords = 0, std::vector* nonzeroCoords = 0) { @@ -178,5 +202,6 @@ initSparse(double density, } } + #include #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 +// Copyright (C) 2012 Gael Guennebaud +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. +#include +#include +#include + +#include "main.h" +#include + +using namespace std; +using namespace Eigen; + +template +struct sparseGaussianTest : SparseFunctor +{ + typedef Matrix VectorType; + typedef SparseFunctor Base; + typedef typename Base::JacobianType JacobianType; + sparseGaussianTest(int inputs, int values) : SparseFunctor(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 u(uv, 0, half); + VectorBlock 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 u(uv, 0, half); + VectorBlock 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 u(uv, 0, half); + VectorBlock 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 +void test_sparseLM_T() +{ + typedef Matrix VectorType; + + int inputs = 10; + int values = 2000; + sparseGaussianTest 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 u(uv, 0, inputs/2); + VectorBlock v(uv, inputs/2, inputs/2); + v.setOnes(); + //Generate u or Solve for u from v + u.setOnes(); + + // Solve the optimization problem + LevenbergMarquardt > 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()); + + // CALL_SUBTEST_2(test_sparseLM_T()); +} 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 // Copyright (C) 2008 Daniel Gomez Ferro +// Copyright (C) 2013 Désiré Nuentsa-Wakam // // 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 void sparse_basic(const SparseMatrixType& ref) { typedef typename SparseMatrixType::Index Index; - + typedef Matrix Vector2; + const Index rows = ref.rows(); const Index cols = ref.cols(); typedef typename SparseMatrixType::Scalar Scalar; @@ -24,71 +26,77 @@ template void sparse_basic(const SparseMatrixType& re typedef Matrix DenseVector; Scalar eps = 1e-6; - SparseMatrixType m(rows, cols); - DenseMatrix refMat = DenseMatrix::Zero(rows, cols); - DenseVector vec1 = DenseVector::Random(rows); Scalar s1 = internal::random(); - - std::vector zeroCoords; - std::vector nonzeroCoords; - initSparse(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 >::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 zeroCoords; + std::vector nonzeroCoords; + initSparse(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(0,cols-1); - int i = internal::random(0,rows-1); - int w = internal::random(1,cols-j-1); - int h = internal::random(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 >::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(0,cols-1); + int i = internal::random(0,rows-1); + int w = internal::random(1,cols-j-1); + int h = internal::random(1,rows-i-1); + + // VERIFY_IS_APPROX(m.block(i,j,h,w), refMat.block(i,j,h,w)); + for(int c=0; c void sparse_basic(const SparseMatrixType& re SparseMatrixType m2(rows,cols); if(internal::random()%2) m2.reserve(VectorXi::Constant(m2.outerSize(), 2)); - for (int j=0; j(0,rows-1); + Index i = internal::random(0,rows-1); if (m1.coeff(i,j)==Scalar(0)) m2.insert(i,j) = m1(i,j) = internal::random(); } @@ -119,8 +127,8 @@ template void sparse_basic(const SparseMatrixType& re m2.reserve(VectorXi::Constant(m2.outerSize(), 2)); for (int k=0; k(0,rows-1); - int j = internal::random(0,cols-1); + Index i = internal::random(0,rows-1); + Index j = internal::random(0,cols-1); if ((m1.coeff(i,j)==Scalar(0)) && (internal::random()%2)) m2.insert(i,j) = m1(i,j) = internal::random(); else @@ -143,8 +151,8 @@ template void sparse_basic(const SparseMatrixType& re m2.reserve(r); for (int k=0; k(0,rows-1); - int j = internal::random(0,cols-1); + Index i = internal::random(0,rows-1); + Index j = internal::random(0,cols-1); if (m1.coeff(i,j)==Scalar(0)) m2.insert(i,j) = m1(i,j) = internal::random(); if(mode==3) @@ -155,6 +163,80 @@ template void sparse_basic(const SparseMatrixType& re VERIFY_IS_APPROX(m2,m1); } + // test innerVector() + { + DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); + SparseMatrixType m2(rows, rows); + initSparse(density, refMat2, m2); + Index j0 = internal::random(0,rows-1); + Index j1 = internal::random(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; j0) + VERIFY(j==numext::real(m3.innerVector(j).lastCoeff())); + } + m3.makeCompressed(); + for(Index j=0; j0) + 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(density, refMat2, m2); + if(internal::random(0,1)>0.5) m2.makeCompressed(); + + Index j0 = internal::random(0,rows-2); + Index j1 = internal::random(0,rows-2); + Index n0 = internal::random(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 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 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(density, refMat2, m2); - int j0 = internal::random(0,rows-1); - int j1 = internal::random(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; j0) - VERIFY(j==internal::real(m3.innerVector(j).lastCoeff())); - } - m3.makeCompressed(); - for(int j=0; j0) - 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(density, refMat2, m2); - int j0 = internal::random(0,rows-2); - int j1 = internal::random(0,rows-2); - int n0 = internal::random(1,rows-(std::max)(j0,j1)); + Index j0 = internal::random(0,rows-2); + Index j1 = internal::random(0,rows-2); + Index n0 = internal::random(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(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 void sparse_basic(const SparseMatrixType& re refM2.setZero(); int countFalseNonZero = 0; int countTrueNonZero = 0; - for (int j=0; j(0,1); if (x<0.1) @@ -320,8 +379,8 @@ template void sparse_basic(const SparseMatrixType& re refMat.setZero(); for(int i=0;i(0,rows-1); - int c = internal::random(0,cols-1); + Index r = internal::random(0,rows-1); + Index c = internal::random(0,cols-1); Scalar v = internal::random(); triplets.push_back(TripletType(r,c,v)); refMat(r,c) += v; @@ -351,6 +410,14 @@ template void sparse_basic(const SparseMatrixType& re refMat3 = refMat2.template triangularView(); m3 = m2.template triangularView(); VERIFY_IS_APPROX(m3, refMat3); + + refMat3 = refMat2.template triangularView(); + m3 = m2.template triangularView(); + VERIFY_IS_APPROX(m3, refMat3); + + refMat3 = refMat2.template triangularView(); + m3 = m2.template triangularView(); + VERIFY_IS_APPROX(m3, refMat3); } // test selfadjointView @@ -379,17 +446,64 @@ template void sparse_basic(const SparseMatrixType& re initSparse(density, refMat2, m2); VERIFY_IS_APPROX(m2.diagonal(), refMat2.diagonal().eval()); } + + // test conservative resize + { + std::vector< std::pair > inc; + inc.push_back(std::pair(-3,-2)); + inc.push_back(std::pair(0,0)); + inc.push_back(std::pair(3,2)); + inc.push_back(std::pair(3,0)); + inc.push_back(std::pair(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(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(1,50); + EIGEN_UNUSED_VARIABLE(s); CALL_SUBTEST_1(( sparse_basic(SparseMatrix(8, 8)) )); CALL_SUBTEST_2(( sparse_basic(SparseMatrix, ColMajor>(s, s)) )); CALL_SUBTEST_2(( sparse_basic(SparseMatrix, RowMajor>(s, s)) )); CALL_SUBTEST_1(( sparse_basic(SparseMatrix(s, s)) )); CALL_SUBTEST_1(( sparse_basic(SparseMatrix(s, s)) )); CALL_SUBTEST_1(( sparse_basic(SparseMatrix(s, s)) )); + + CALL_SUBTEST_1(( sparse_basic(SparseMatrix(short(s), short(s))) )); + CALL_SUBTEST_1(( sparse_basic(SparseMatrix(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 struct test_outer { 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(0,m2.cols()-1); + Index c1 = internal::random(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 struct test_outer struct test_outer { 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(0,m2.rows()-1); + Index c1 = internal::random(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 void sparse_product() { typedef typename SparseMatrixType::Index Index; Index n = 100; - const Index rows = internal::random(1,n); - const Index cols = internal::random(1,n); - const Index depth = internal::random(1,n); + const Index rows = internal::random(1,n); + const Index cols = internal::random(1,n); + const Index depth = internal::random(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 DenseMatrix; typedef Matrix DenseVector; + typedef Matrix RowDenseVector; + typedef SparseVector ColSpVector; + typedef SparseVector RowSpVector; Scalar s1 = internal::random(); Scalar s2 = internal::random(); @@ -117,21 +122,54 @@ template void sparse_product() test_outer::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 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 d1(DenseVector::Random(cols)); + DiagonalMatrix d2(DenseVector::Random(rows)); + SparseMatrixType m2(rows, cols); + SparseMatrixType m3(rows, cols); initSparse(density, refM2, m2); initSparse(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 void sparse_product() VERIFY_IS_APPROX(x=mUp.template selfadjointView()*b, refX=refS*b); VERIFY_IS_APPROX(x=mLo.template selfadjointView()*b, refX=refS*b); VERIFY_IS_APPROX(x=mS.template selfadjointView()*b, refX=refS*b); + + // sparse selfadjointView * sparse + SparseMatrixType mSres(rows,rows); + VERIFY_IS_APPROX(mSres = mLo.template selfadjointView()*mS, + refX = refLo.template selfadjointView()*refS); + // sparse * sparse selfadjointview + VERIFY_IS_APPROX(mSres = mS * mLo.template selfadjointView(), + refX = refS * refLo.template selfadjointView()); } + } // New test for Bug in SparseTimeDenseProduct @@ -199,6 +246,7 @@ void test_sparse_product() CALL_SUBTEST_1( (sparse_product >()) ); CALL_SUBTEST_2( (sparse_product, ColMajor > >()) ); CALL_SUBTEST_2( (sparse_product, RowMajor > >()) ); + CALL_SUBTEST_3( (sparse_product >()) ); CALL_SUBTEST_4( (sparse_product_regression_test, Matrix >()) ); } } 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())); - - 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())); + 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())); + VERIFY(x.isApprox(refX,test_precision())); + } - // 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())); } } @@ -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 >::value || internal::is_same >::value ) - mat_folder = mat_folder + static_cast("/complex/"); + mat_folder = mat_folder + static_cast("/complex/"); else - mat_folder = mat_folder + static_cast("/real/"); + mat_folder = mat_folder + static_cast("/real/"); return mat_folder; } #endif @@ -169,7 +168,6 @@ template void check_sparse_spd_solving(Solver& solver) { typedef typename Solver::MatrixType Mat; typedef typename Mat::Scalar Scalar; - typedef typename Mat::Index Index; typedef SparseMatrix SpMat; typedef Matrix DenseMatrix; typedef Matrix DenseVector; @@ -177,25 +175,32 @@ template 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(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(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(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(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::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 DenseMatrix; int size = internal::random(1,maxSize); double density = (std::max)(8./(size*size), 0.01); @@ -258,6 +262,7 @@ template void check_sparse_square_solving(Solver& solver) { typedef typename Solver::MatrixType Mat; typedef typename Mat::Scalar Scalar; + typedef SparseMatrix SpMat; typedef Matrix DenseMatrix; typedef Matrix DenseVector; @@ -265,16 +270,28 @@ template 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(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::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 void sparse_vector(int rows, int cols) +template 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 DenseMatrix; typedef Matrix DenseVector; - typedef SparseVector SparseVectorType; - typedef SparseMatrix SparseMatrixType; + typedef SparseVector SparseVectorType; + typedef SparseMatrix SparseMatrixType; Scalar eps = 1e-6; SparseMatrixType m1(rows,rows); @@ -77,15 +77,34 @@ template 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(8, 8) ); - CALL_SUBTEST_2( sparse_vector >(16, 16) ); - CALL_SUBTEST_1( sparse_vector(299, 535) ); + CALL_SUBTEST_1(( sparse_vector(8, 8) )); + CALL_SUBTEST_2(( sparse_vector, int>(16, 16) )); + CALL_SUBTEST_1(( sparse_vector(299, 535) )); + CALL_SUBTEST_1(( sparse_vector(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 +// +// 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 . + + +// 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 +#include + +template void test_sparselu_T() +{ + SparseLU /*, COLAMDOrdering*/ > sparselu_colamd; // COLAMDOrdering is the default + SparseLU, AMDOrdering > sparselu_amd; + SparseLU, NaturalOrdering > 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()); + CALL_SUBTEST_2(test_sparselu_T()); + CALL_SUBTEST_3(test_sparselu_T >()); + CALL_SUBTEST_4(test_sparselu_T >()); +} 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 +// Copyright (C) 2014 Gael Guennebaud +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +#include "sparse.h" +#include + +template +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(1,maxRows); + int cols = internal::random(1,maxCols); + double density = (std::max)(8./(rows*cols), 0.01); + + A.resize(rows,cols); + dA.resize(rows,cols); + initSparse(density, dA, A,ForceNonZeroDiag); + A.makeCompressed(); + int nop = internal::random(0, internal::random(0,1) > 0.5 ? cols/2 : 0); + for(int k=0; k(0,cols-1); + int j1 = internal::random(0,cols-1); + Scalar s = internal::random(); + A.col(j0) = s * A.col(j1); + dA.col(j0) = s * dA.col(j1); + } + +// if(rows void test_sparseqr_scalar() +{ + typedef SparseMatrix MatrixType; + typedef Matrix DenseMat; + typedef Matrix DenseVector; + MatrixType A; + DenseMat dA; + DenseVector refX,x,b; + SparseQR > 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 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()); + CALL_SUBTEST_2(test_sparseqr_scalar >()); + } +} + 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 +// +// 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 void special_numbers() +{ + typedef Matrix MatType; + int rows = internal::random(1,300); + int cols = internal::random(1,300); + + Scalar nan = std::numeric_limits::quiet_NaN(); + Scalar inf = std::numeric_limits::infinity(); + Scalar s1 = internal::random(); + + 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(1,10); + for(int k=0; k(0,rows-1), internal::random(0,cols-1)) = nan; + minf(internal::random(0,rows-1), internal::random(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() ); + CALL_SUBTEST_1( special_numbers() ); + } +} 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 +// +// 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 + + +template +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(1,maxRows); + int cols = internal::random(1,rows); + double density = (std::max)(8./(rows*cols), 0.01); + + A.resize(rows,rows); + dA.resize(rows,rows); + initSparse(density, dA, A,ForceNonZeroDiag); + A.makeCompressed(); + return rows; +} + +template void test_spqr_scalar() +{ + typedef SparseMatrix MatrixType; + MatrixType A; + Matrix dA; + typedef Matrix DenseVector; + DenseVector refX,x,b; + SPQR 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())); +} +void test_spqr_support() +{ + CALL_SUBTEST_1(test_spqr_scalar()); + CALL_SUBTEST_2(test_spqr_scalar >()); +} 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 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::Real RealScalar; @@ -53,8 +55,16 @@ template void stable_norm(const MatrixType& m) Index rows = m.rows(); Index cols = m.cols(); - Scalar big = internal::random() * ((std::numeric_limits::max)() * RealScalar(1e-4)); - Scalar small = internal::random() * ((std::numeric_limits::min)() * RealScalar(1e4)); + // get a non-zero random factor + Scalar factor = internal::random(); + while(numext::abs2(factor)(); + Scalar big = factor * ((std::numeric_limits::max)() * RealScalar(1e-4)); + + factor = internal::random(); + while(numext::abs2(factor)(); + Scalar small = factor * ((std::numeric_limits::min)() * RealScalar(1e4)); MatrixType vzero = MatrixType::Zero(rows, cols), vrand = MatrixType::Random(rows, cols), @@ -73,23 +83,23 @@ template void stable_norm(const MatrixType& m) // test isFinite VERIFY(!isFinite( std::numeric_limits::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 void triangular_square(const MatrixType& m) m1 = MatrixType::Random(rows, cols); for (int i=0; i(); + while (numext::abs2(m1(i,i))<1e-1) m1(i,i) = internal::random(); Transpose trm4(m4); // test back and forward subsitution with a vector as the rhs @@ -123,9 +123,6 @@ template void triangular_rect(const MatrixType& m) typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime }; - typedef Matrix VectorType; - typedef Matrix 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(2,maxsize); EIGEN_UNUSED_VARIABLE(r); - int c = internal::random(2,maxsize); EIGEN_UNUSED_VARIABLE(c); + int r = internal::random(2,maxsize); TEST_SET_BUT_UNUSED_VARIABLE(r) + int c = internal::random(2,maxsize); TEST_SET_BUT_UNUSED_VARIABLE(c) CALL_SUBTEST_1( triangular_square(Matrix()) ); CALL_SUBTEST_2( triangular_square(Matrix()) ); 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 Eigen::Matrix randMatrixUnitary(int size) { typedef T Scalar; - typedef typename NumTraits::Real RealScalar; - typedef Eigen::Matrix MatrixType; MatrixType Q; @@ -77,7 +75,6 @@ template Eigen::Matrix randMatrixSpecialUnitary(int size) { typedef T Scalar; - typedef typename NumTraits::Real RealScalar; typedef Eigen::Matrix MatrixType; @@ -85,7 +82,7 @@ Eigen::Matrix randMatrixSpecialUnitary(int si MatrixType Q = randMatrixUnitary(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 randMatrixSpecialUnitary(int si template void run_test(int dim, int num_elements) { + using std::abs; typedef typename internal::traits::Scalar Scalar; typedef Matrix MatrixX; typedef Matrix 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()); + const Scalar c = abs(internal::random()); MatrixX R = randMatrixSpecialUnitary(dim); VectorX t = Scalar(50)*VectorX::Random(dim,1); @@ -122,6 +120,7 @@ void run_test(int dim, int num_elements) template void run_fixed_size_test(int num_elements) { + using std::abs; typedef Matrix MatrixX; typedef Matrix HomMatrix; typedef Matrix 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()); + // Also if c is to small compared to t.norm(), problem is ill-posed (cf. Bug 744) + const Scalar c = internal::random(0.5, 2.0); FixedMatrix R = randMatrixSpecialUnitary(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::epsilon()); + VERIFY(error < Scalar(16)*std::numeric_limits::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 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 RealMatrixType; + typedef Matrix TransposeMatrixType; MatrixType a = MatrixType::Random(rows,cols); internal::UpperBidiagonalization ubd(a); @@ -25,6 +25,8 @@ template 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 void vectorwiseop_array(const ArrayType& m) { typedef typename ArrayType::Index Index; typedef typename ArrayType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; typedef Array ColVectorType; typedef Array RowVectorType; @@ -102,6 +101,16 @@ template 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 void vectorwiseop_matrix(const MatrixType& m) @@ -111,6 +120,8 @@ template void vectorwiseop_matrix(const MatrixType& m) typedef typename NumTraits::Real RealScalar; typedef Matrix ColVectorType; typedef Matrix RowVectorType; + typedef Matrix RealColVectorType; + typedef Matrix RealRowVectorType; Index rows = m.rows(); Index cols = m.cols(); @@ -123,6 +134,8 @@ template 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 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 void vectorVisitor(const VectorType& w) while(v(i) == v(i2)) // yes, == v(i) = internal::random(); - 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 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(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()) ); 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 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 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 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 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 >(); zeroSizedMatrix >(); zeroSizedMatrix >(); - + zeroSizedMatrix >(); + zeroSizedMatrix >(); + zeroSizedVector(); zeroSizedVector(); zeroSizedVector(); zeroSizedVector >(); + zeroSizedVector >(); } -- cgit v1.2.3