diff options
Diffstat (limited to 'test')
| -rw-r--r-- | test/CMakeLists.txt | 5 | ||||
| -rw-r--r-- | test/adjoint.cpp | 9 | ||||
| -rw-r--r-- | test/bandmatrix.cpp | 12 | ||||
| -rw-r--r-- | test/basicstuff.cpp | 9 | ||||
| -rw-r--r-- | test/cholesky.cpp | 18 | ||||
| -rw-r--r-- | test/conservative_resize.cpp | 129 | ||||
| -rw-r--r-- | test/eigensolver_complex.cpp | 61 | ||||
| -rw-r--r-- | test/eigensolver_generic.cpp | 2 | ||||
| -rw-r--r-- | test/eigensolver_selfadjoint.cpp | 2 | ||||
| -rw-r--r-- | test/geo_orthomethods.cpp | 6 | ||||
| -rw-r--r-- | test/jacobisvd.cpp | 106 | ||||
| -rw-r--r-- | test/main.h | 12 | ||||
| -rw-r--r-- | test/mixingtypes.cpp | 127 | ||||
| -rw-r--r-- | test/product.h | 12 | ||||
| -rw-r--r-- | test/product_extra.cpp | 36 | ||||
| -rw-r--r-- | test/product_notemporary.cpp | 44 | ||||
| -rw-r--r-- | test/product_symm.cpp | 2 | ||||
| -rw-r--r-- | test/product_trsm.cpp | 4 | ||||
| -rw-r--r-- | test/stable_norm.cpp | 79 | ||||
| -rw-r--r-- | test/submatrices.cpp | 45 | ||||
| -rw-r--r-- | test/svd.cpp | 6 | ||||
| -rw-r--r-- | test/swap.cpp | 98 | ||||
| -rw-r--r-- | test/testsuite.cmake | 38 | ||||
| -rw-r--r-- | test/umeyama.cpp | 4 |
24 files changed, 745 insertions, 121 deletions
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index fc43bbb1d..4e279ea47 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -115,6 +115,7 @@ ei_add_test(product_trmv ${EI_OFLAG}) ei_add_test(product_trmm ${EI_OFLAG}) ei_add_test(product_trsm ${EI_OFLAG}) ei_add_test(product_notemporary ${EI_OFLAG}) +ei_add_test(stable_norm) ei_add_test(bandmatrix) ei_add_test(cholesky " " "${GSL_LIBRARIES}") ei_add_test(lu ${EI_OFLAG}) @@ -125,7 +126,9 @@ ei_add_test(qr_colpivoting) ei_add_test(qr_fullpivoting) ei_add_test(eigensolver_selfadjoint " " "${GSL_LIBRARIES}") ei_add_test(eigensolver_generic " " "${GSL_LIBRARIES}") +ei_add_test(eigensolver_complex) ei_add_test(svd) +ei_add_test(jacobisvd ${EI_OFLAG}) ei_add_test(geo_orthomethods) ei_add_test(geo_homogeneous) ei_add_test(geo_quaternion) @@ -146,6 +149,8 @@ ei_add_test(sparse_product) ei_add_test(sparse_solvers " " "${SPARSE_LIBS}") ei_add_test(umeyama) ei_add_test(householder) +ei_add_test(swap) +ei_add_test(conservative_resize) ei_add_property(EIGEN_TESTING_SUMMARY "CXX: ${CMAKE_CXX_COMPILER}\n") if(CMAKE_COMPILER_IS_GNUCXX) diff --git a/test/adjoint.cpp b/test/adjoint.cpp index 964658c65..bebf47ac3 100644 --- a/test/adjoint.cpp +++ b/test/adjoint.cpp @@ -72,13 +72,6 @@ template<typename MatrixType> void adjoint(const MatrixType& m) if(NumTraits<Scalar>::HasFloatingPoint) VERIFY_IS_APPROX(v1.squaredNorm(), v1.norm() * v1.norm()); VERIFY_IS_MUCH_SMALLER_THAN(ei_abs(vzero.dot(v1)), static_cast<RealScalar>(1)); - if(NumTraits<Scalar>::HasFloatingPoint) - { - VERIFY_IS_MUCH_SMALLER_THAN(vzero.norm(), static_cast<RealScalar>(1)); - VERIFY_IS_APPROX(v1.norm(), v1.stableNorm()); - VERIFY_IS_APPROX(v1.blueNorm(), v1.stableNorm()); - VERIFY_IS_APPROX(v1.hypotNorm(), v1.stableNorm()); - } // check compatibility of dot and adjoint VERIFY(ei_isApprox(v1.dot(square * v2), (square.adjoint() * v1).dot(v2), largerEps)); @@ -124,7 +117,7 @@ void test_adjoint() } // test a large matrix only once CALL_SUBTEST( adjoint(Matrix<float, 100, 100>()) ); - + { MatrixXcf a(10,10), b(10,10); VERIFY_RAISES_ASSERT(a = a.transpose()); diff --git a/test/bandmatrix.cpp b/test/bandmatrix.cpp index f677e7b85..2bdc67e28 100644 --- a/test/bandmatrix.cpp +++ b/test/bandmatrix.cpp @@ -44,21 +44,21 @@ template<typename MatrixType> void bandmatrix(const MatrixType& _m) dm1.diagonal().setConstant(123); for (int i=1; i<=m.supers();++i) { - m.diagonal(i).setConstant(i); - dm1.diagonal(i).setConstant(i); + m.diagonal(i).setConstant(static_cast<RealScalar>(i)); + dm1.diagonal(i).setConstant(static_cast<RealScalar>(i)); } for (int i=1; i<=m.subs();++i) { - m.diagonal(-i).setConstant(-i); - dm1.diagonal(-i).setConstant(-i); + m.diagonal(-i).setConstant(-static_cast<RealScalar>(i)); + dm1.diagonal(-i).setConstant(-static_cast<RealScalar>(i)); } //std::cerr << m.m_data << "\n\n" << m.toDense() << "\n\n" << dm1 << "\n\n\n\n"; VERIFY_IS_APPROX(dm1,m.toDense()); for (int i=0; i<cols; ++i) { - m.col(i).setConstant(i+1); - dm1.col(i).setConstant(i+1); + m.col(i).setConstant(static_cast<RealScalar>(i+1)); + dm1.col(i).setConstant(static_cast<RealScalar>(i+1)); } int d = std::min(rows,cols); int a = std::max(0,cols-d-supers); diff --git a/test/basicstuff.cpp b/test/basicstuff.cpp index 52aff93ee..29df99f9e 100644 --- a/test/basicstuff.cpp +++ b/test/basicstuff.cpp @@ -99,15 +99,6 @@ template<typename MatrixType> void basicStuff(const MatrixType& m) MatrixType m4; VERIFY_IS_APPROX(m4 = m1,m1); - // test swap - m3 = m1; - m1.swap(m2); - VERIFY_IS_APPROX(m3, m2); - if(rows*cols>=3) - { - VERIFY_IS_NOT_APPROX(m3, m1); - } - m3.real() = m1.real(); VERIFY_IS_APPROX(static_cast<const MatrixType&>(m3).real(), static_cast<const MatrixType&>(m1).real()); VERIFY_IS_APPROX(static_cast<const MatrixType&>(m3).real(), m1.real()); diff --git a/test/cholesky.cpp b/test/cholesky.cpp index df937fd0f..526a9f9d0 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -82,7 +82,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m) // // test gsl itself ! // VERIFY_IS_APPROX(vecB, _vecB); // VERIFY_IS_APPROX(vecX, _vecX); -// +// // Gsl::free(gMatA); // Gsl::free(gSymm); // Gsl::free(gVecB); @@ -149,16 +149,16 @@ void test_cholesky() { for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST( cholesky(Matrix<double,1,1>()) ); -// CALL_SUBTEST( cholesky(MatrixXd(1,1)) ); -// CALL_SUBTEST( cholesky(Matrix2d()) ); -// CALL_SUBTEST( cholesky(Matrix3f()) ); -// CALL_SUBTEST( cholesky(Matrix4d()) ); + CALL_SUBTEST( cholesky(MatrixXd(1,1)) ); + CALL_SUBTEST( cholesky(Matrix2d()) ); + CALL_SUBTEST( cholesky(Matrix3f()) ); + CALL_SUBTEST( cholesky(Matrix4d()) ); CALL_SUBTEST( cholesky(MatrixXd(200,200)) ); CALL_SUBTEST( cholesky(MatrixXcd(100,100)) ); } -// CALL_SUBTEST( cholesky_verify_assert<Matrix3f>() ); -// CALL_SUBTEST( cholesky_verify_assert<Matrix3d>() ); -// CALL_SUBTEST( cholesky_verify_assert<MatrixXf>() ); -// CALL_SUBTEST( cholesky_verify_assert<MatrixXd>() ); + CALL_SUBTEST( cholesky_verify_assert<Matrix3f>() ); + CALL_SUBTEST( cholesky_verify_assert<Matrix3d>() ); + CALL_SUBTEST( cholesky_verify_assert<MatrixXf>() ); + CALL_SUBTEST( cholesky_verify_assert<MatrixXd>() ); } diff --git a/test/conservative_resize.cpp b/test/conservative_resize.cpp new file mode 100644 index 000000000..b92dd5449 --- /dev/null +++ b/test/conservative_resize.cpp @@ -0,0 +1,129 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2009 Hauke Heibel <hauke.heibel@gmail.com> +// +// 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 or1 FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#include "main.h" + +#include <Eigen/Core> +#include <Eigen/Array> + +using namespace Eigen; + +template <typename Scalar, int Storage> +void run_matrix_tests() +{ + typedef Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Storage> MatrixType; + + MatrixType m, n; + + // boundary cases ... + m = n = MatrixType::Random(50,50); + m.conservativeResize(1,50); + VERIFY_IS_APPROX(m, n.block(0,0,1,50)); + + m = n = MatrixType::Random(50,50); + m.conservativeResize(50,1); + VERIFY_IS_APPROX(m, n.block(0,0,50,1)); + + m = n = MatrixType::Random(50,50); + m.conservativeResize(50,50); + VERIFY_IS_APPROX(m, n.block(0,0,50,50)); + + // random shrinking ... + for (int i=0; i<25; ++i) + { + const int rows = ei_random<int>(1,50); + const int cols = ei_random<int>(1,50); + m = n = MatrixType::Random(50,50); + m.conservativeResize(rows,cols); + VERIFY_IS_APPROX(m, n.block(0,0,rows,cols)); + } + + // random growing with zeroing ... + for (int i=0; i<25; ++i) + { + const int rows = ei_random<int>(50,75); + const int cols = ei_random<int>(50,75); + m = n = MatrixType::Random(50,50); + m.conservativeResizeLike(MatrixType::Zero(rows,cols)); + VERIFY_IS_APPROX(m.block(0,0,n.rows(),n.cols()), n); + VERIFY( rows<=50 || m.block(50,0,rows-50,cols).sum() == Scalar(0) ); + VERIFY( cols<=50 || m.block(0,50,rows,cols-50).sum() == Scalar(0) ); + } +} + +template <typename Scalar> +void run_vector_tests() +{ + typedef Matrix<Scalar, 1, Eigen::Dynamic> MatrixType; + + MatrixType m, n; + + // boundary cases ... + m = n = MatrixType::Random(50); + m.conservativeResize(1); + VERIFY_IS_APPROX(m, n.segment(0,1)); + + m = n = MatrixType::Random(50); + m.conservativeResize(50); + VERIFY_IS_APPROX(m, n.segment(0,50)); + + // random shrinking ... + for (int i=0; i<50; ++i) + { + const int size = ei_random<int>(1,50); + m = n = MatrixType::Random(50); + m.conservativeResize(size); + VERIFY_IS_APPROX(m, n.segment(0,size)); + } + + // random growing with zeroing ... + for (int i=0; i<50; ++i) + { + const int size = ei_random<int>(50,100); + m = n = MatrixType::Random(50); + m.conservativeResizeLike(MatrixType::Zero(size)); + VERIFY_IS_APPROX(m.segment(0,50), n); + VERIFY( size<=50 || m.segment(50,size-50).sum() == Scalar(0) ); + } +} + +void test_conservative_resize() +{ + run_matrix_tests<int, Eigen::RowMajor>(); + run_matrix_tests<int, Eigen::ColMajor>(); + run_matrix_tests<float, Eigen::RowMajor>(); + run_matrix_tests<float, Eigen::ColMajor>(); + run_matrix_tests<double, Eigen::RowMajor>(); + run_matrix_tests<double, Eigen::ColMajor>(); + run_matrix_tests<std::complex<float>, Eigen::RowMajor>(); + run_matrix_tests<std::complex<float>, Eigen::ColMajor>(); + run_matrix_tests<std::complex<double>, Eigen::RowMajor>(); + run_matrix_tests<std::complex<double>, Eigen::ColMajor>(); + + run_vector_tests<int>(); + run_vector_tests<float>(); + run_vector_tests<double>(); + run_vector_tests<std::complex<float> >(); + run_vector_tests<std::complex<double> >(); +} diff --git a/test/eigensolver_complex.cpp b/test/eigensolver_complex.cpp new file mode 100644 index 000000000..38ede7c4a --- /dev/null +++ b/test/eigensolver_complex.cpp @@ -0,0 +1,61 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#include "main.h" +#include <Eigen/Eigenvalues> +#include <Eigen/LU> + +template<typename MatrixType> void eigensolver(const MatrixType& m) +{ + /* this test covers the following files: + ComplexEigenSolver.h, and indirectly ComplexSchur.h + */ + int rows = m.rows(); + int cols = m.cols(); + + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; + typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealVectorType; + typedef typename std::complex<typename NumTraits<typename MatrixType::Scalar>::Real> Complex; + + MatrixType a = MatrixType::Random(rows,cols); + MatrixType symmA = a.adjoint() * a; + + ComplexEigenSolver<MatrixType> ei0(symmA); + VERIFY_IS_APPROX(symmA * ei0.eigenvectors(), ei0.eigenvectors() * ei0.eigenvalues().asDiagonal()); + + ComplexEigenSolver<MatrixType> ei1(a); + VERIFY_IS_APPROX(a * ei1.eigenvectors(), ei1.eigenvectors() * ei1.eigenvalues().asDiagonal()); + +} + +void test_eigensolver_complex() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST( eigensolver(Matrix4cf()) ); + CALL_SUBTEST( eigensolver(MatrixXcd(14,14)) ); + } +} + diff --git a/test/eigensolver_generic.cpp b/test/eigensolver_generic.cpp index 2be49faf4..e2b2055b4 100644 --- a/test/eigensolver_generic.cpp +++ b/test/eigensolver_generic.cpp @@ -23,7 +23,7 @@ // Eigen. If not, see <http://www.gnu.org/licenses/>. #include "main.h" -#include <Eigen/QR> +#include <Eigen/Eigenvalues> #ifdef HAS_GSL #include "gsl_helper.h" diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp index 2571dbf43..3836c074b 100644 --- a/test/eigensolver_selfadjoint.cpp +++ b/test/eigensolver_selfadjoint.cpp @@ -23,7 +23,7 @@ // Eigen. If not, see <http://www.gnu.org/licenses/>. #include "main.h" -#include <Eigen/QR> +#include <Eigen/Eigenvalues> #ifdef HAS_GSL #include "gsl_helper.h" diff --git a/test/geo_orthomethods.cpp b/test/geo_orthomethods.cpp index 5e1d5bdb4..540a63b82 100644 --- a/test/geo_orthomethods.cpp +++ b/test/geo_orthomethods.cpp @@ -86,10 +86,10 @@ template<typename Scalar, int Size> void orthomethods(int size=Size) VERIFY_IS_MUCH_SMALLER_THAN(v0.unitOrthogonal().dot(v0), Scalar(1)); VERIFY_IS_APPROX(v0.unitOrthogonal().norm(), RealScalar(1)); - if (size>3) + if (size>=3) { - v0.template start<3>().setZero(); - v0.end(size-3).setRandom(); + v0.template start<2>().setZero(); + v0.end(size-2).setRandom(); VERIFY_IS_MUCH_SMALLER_THAN(v0.unitOrthogonal().dot(v0), Scalar(1)); VERIFY_IS_APPROX(v0.unitOrthogonal().norm(), RealScalar(1)); diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp new file mode 100644 index 000000000..5940b8497 --- /dev/null +++ b/test/jacobisvd.cpp @@ -0,0 +1,106 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#include "main.h" +#include <Eigen/SVD> +#include <Eigen/LU> + +template<typename MatrixType, unsigned int Options> void svd(const MatrixType& m = MatrixType(), bool pickrandom = true) +{ + int rows = m.rows(); + int cols = m.cols(); + + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime + }; + + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixUType; + typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime> MatrixVType; + typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType; + typedef Matrix<Scalar, ColsAtCompileTime, 1> InputVectorType; + + MatrixType a; + if(pickrandom) a = MatrixType::Random(rows,cols); + else a = m; + + JacobiSVD<MatrixType,Options> svd(a); + MatrixType sigma = MatrixType::Zero(rows,cols); + sigma.diagonal() = svd.singularValues().template cast<Scalar>(); + MatrixUType u = svd.matrixU(); + MatrixVType v = svd.matrixV(); + + VERIFY_IS_APPROX(a, u * sigma * v.adjoint()); + VERIFY_IS_UNITARY(u); + VERIFY_IS_UNITARY(v); +} + +template<typename MatrixType> void svd_verify_assert() +{ + MatrixType tmp; + + SVD<MatrixType> svd; + //VERIFY_RAISES_ASSERT(svd.solve(tmp, &tmp)) + VERIFY_RAISES_ASSERT(svd.matrixU()) + VERIFY_RAISES_ASSERT(svd.singularValues()) + VERIFY_RAISES_ASSERT(svd.matrixV()) + /*VERIFY_RAISES_ASSERT(svd.computeUnitaryPositive(&tmp,&tmp)) + VERIFY_RAISES_ASSERT(svd.computePositiveUnitary(&tmp,&tmp)) + VERIFY_RAISES_ASSERT(svd.computeRotationScaling(&tmp,&tmp)) + VERIFY_RAISES_ASSERT(svd.computeScalingRotation(&tmp,&tmp))*/ +} + +void test_jacobisvd() +{ + for(int i = 0; i < g_repeat; i++) { + Matrix2cd m; + m << 0, 1, + 0, 1; + CALL_SUBTEST(( svd<Matrix2cd,0>(m, false) )); + m << 1, 0, + 1, 0; + CALL_SUBTEST(( svd<Matrix2cd,0>(m, false) )); + Matrix2d n; + n << 1, 1, + 1, -1; + CALL_SUBTEST(( svd<Matrix2d,0>(n, false) )); + CALL_SUBTEST(( svd<Matrix3f,0>() )); + CALL_SUBTEST(( svd<Matrix4d,Square>() )); + CALL_SUBTEST(( svd<Matrix<float,3,5> , AtLeastAsManyColsAsRows>() )); + CALL_SUBTEST(( svd<Matrix<double,Dynamic,2> , AtLeastAsManyRowsAsCols>(Matrix<double,Dynamic,2>(10,2)) )); + + CALL_SUBTEST(( svd<MatrixXf,Square>(MatrixXf(50,50)) )); + CALL_SUBTEST(( svd<MatrixXcd,AtLeastAsManyRowsAsCols>(MatrixXcd(14,7)) )); + } + CALL_SUBTEST(( svd<MatrixXf,0>(MatrixXf(300,200)) )); + CALL_SUBTEST(( svd<MatrixXcd,AtLeastAsManyColsAsRows>(MatrixXcd(100,150)) )); + + CALL_SUBTEST(( svd_verify_assert<Matrix3f>() )); + CALL_SUBTEST(( svd_verify_assert<Matrix3d>() )); + CALL_SUBTEST(( svd_verify_assert<MatrixXf>() )); + CALL_SUBTEST(( svd_verify_assert<MatrixXd>() )); +} diff --git a/test/main.h b/test/main.h index e3866c68b..619fc9e06 100644 --- a/test/main.h +++ b/test/main.h @@ -30,6 +30,10 @@ #include <vector> #include <typeinfo> +#ifdef NDEBUG +#undef NDEBUG +#endif + #ifndef EIGEN_TEST_FUNC #error EIGEN_TEST_FUNC must be defined #endif @@ -153,6 +157,8 @@ namespace Eigen #define VERIFY_IS_APPROX_OR_LESS_THAN(a, b) VERIFY(test_ei_isApproxOrLessThan(a, b)) #define VERIFY_IS_NOT_APPROX_OR_LESS_THAN(a, b) VERIFY(!test_ei_isApproxOrLessThan(a, b)) +#define VERIFY_IS_UNITARY(a) VERIFY(test_isUnitary(a)) + #define CALL_SUBTEST(FUNC) do { \ g_test_stack.push_back(EI_PP_MAKE_STRING(FUNC)); \ FUNC; \ @@ -227,6 +233,12 @@ inline bool test_ei_isMuchSmallerThan(const MatrixBase<Derived>& m, return m.isMuchSmallerThan(s, test_precision<typename ei_traits<Derived>::Scalar>()); } +template<typename Derived> +inline bool test_isUnitary(const MatrixBase<Derived>& m) +{ + return m.isUnitary(test_precision<typename ei_traits<Derived>::Scalar>()); +} + template<typename MatrixType> void createRandomMatrixOfRank(int desired_rank, int rows, int cols, MatrixType& m) { diff --git a/test/mixingtypes.cpp b/test/mixingtypes.cpp index d14232bd4..3e322c7fe 100644 --- a/test/mixingtypes.cpp +++ b/test/mixingtypes.cpp @@ -33,6 +33,7 @@ #include "main.h" +using namespace std; template<int SizeAtCompileType> void mixingtypes(int size = SizeAtCompileType) { @@ -45,6 +46,104 @@ template<int SizeAtCompileType> void mixingtypes(int size = SizeAtCompileType) typedef Matrix<std::complex<float>, SizeAtCompileType, 1> Vec_cf; typedef Matrix<std::complex<double>, SizeAtCompileType, 1> Vec_cd; + Mat_f mf = Mat_f::Random(size,size); + Mat_d md = mf.template cast<double>(); + Mat_cf mcf = Mat_cf::Random(size,size); + Mat_cd mcd = mcf.template cast<complex<double> >(); + Vec_f vf = Vec_f::Random(size,1); + Vec_d vd = vf.template cast<double>(); + Vec_cf vcf = Vec_cf::Random(size,1); + Vec_cd vcd = vcf.template cast<complex<double> >(); + float sf = ei_random<float>(); + double sd = ei_random<double>(); + complex<float> scf = ei_random<complex<float> >(); + complex<double> scd = ei_random<complex<double> >(); + + + mf+mf; + VERIFY_RAISES_ASSERT(mf+md); + VERIFY_RAISES_ASSERT(mf+mcf); + VERIFY_RAISES_ASSERT(vf=vd); + VERIFY_RAISES_ASSERT(vf+=vd); + VERIFY_RAISES_ASSERT(mcd=md); + + // check scalar products + VERIFY_IS_APPROX(vcf * sf , vcf * complex<float>(sf)); + VERIFY_IS_APPROX(sd * vcd, complex<double>(sd) * vcd); + VERIFY_IS_APPROX(vf * scf , vf.template cast<complex<float> >() * scf); + VERIFY_IS_APPROX(scd * vd, scd * vd.template cast<complex<double> >()); + + // check dot product + vf.dot(vf); + VERIFY_RAISES_ASSERT(vd.dot(vf)); + VERIFY_RAISES_ASSERT(vcf.dot(vf)); // yeah eventually we should allow this but i'm too lazy to make that change now in Dot.h + // especially as that might be rewritten as cwise product .sum() which would make that automatic. + + // check diagonal product + VERIFY_IS_APPROX(vf.asDiagonal() * mcf, vf.template cast<complex<float> >().asDiagonal() * mcf); + VERIFY_IS_APPROX(vcd.asDiagonal() * md, vcd.asDiagonal() * md.template cast<complex<double> >()); + VERIFY_IS_APPROX(mcf * vf.asDiagonal(), mcf * vf.template cast<complex<float> >().asDiagonal()); + VERIFY_IS_APPROX(md * vcd.asDiagonal(), md.template cast<complex<double> >() * vcd.asDiagonal()); +// vd.asDiagonal() * mf; // does not even compile +// vcd.asDiagonal() * mf; // does not even compile + + // check inner product + VERIFY_IS_APPROX((vf.transpose() * vcf).value(), (vf.template cast<complex<float> >().transpose() * vcf).value()); + + // check outer product + VERIFY_IS_APPROX((vf * vcf.transpose()).eval(), (vf.template cast<complex<float> >() * vcf.transpose()).eval()); +} + + +void mixingtypes_large(int size) +{ + static const int SizeAtCompileType = Dynamic; + typedef Matrix<float, SizeAtCompileType, SizeAtCompileType> Mat_f; + typedef Matrix<double, SizeAtCompileType, SizeAtCompileType> Mat_d; + typedef Matrix<std::complex<float>, SizeAtCompileType, SizeAtCompileType> Mat_cf; + typedef Matrix<std::complex<double>, SizeAtCompileType, SizeAtCompileType> Mat_cd; + typedef Matrix<float, SizeAtCompileType, 1> Vec_f; + typedef Matrix<double, SizeAtCompileType, 1> Vec_d; + typedef Matrix<std::complex<float>, SizeAtCompileType, 1> Vec_cf; + typedef Matrix<std::complex<double>, SizeAtCompileType, 1> Vec_cd; + + Mat_f mf(size,size); + Mat_d md(size,size); + Mat_cf mcf(size,size); + Mat_cd mcd(size,size); + Vec_f vf(size,1); + Vec_d vd(size,1); + Vec_cf vcf(size,1); + Vec_cd vcd(size,1); + + mf*mf; + // FIXME large products does not allow mixing types + VERIFY_RAISES_ASSERT(md*mcd); + VERIFY_RAISES_ASSERT(mcd*md); + VERIFY_RAISES_ASSERT(mf*vcf); + VERIFY_RAISES_ASSERT(mcf*vf); + VERIFY_RAISES_ASSERT(mcf *= mf); + // VERIFY_RAISES_ASSERT(vcd = md*vcd); // does not even compile (cannot convert complex to double) + VERIFY_RAISES_ASSERT(vcf = mcf*vf); + +// VERIFY_RAISES_ASSERT(mf*md); // does not even compile +// VERIFY_RAISES_ASSERT(mcf*mcd); // does not even compile +// VERIFY_RAISES_ASSERT(mcf*vcd); // does not even compile + VERIFY_RAISES_ASSERT(vcf = mf*vf); +} + +template<int SizeAtCompileType> void mixingtypes_small() +{ + int size = SizeAtCompileType; + typedef Matrix<float, SizeAtCompileType, SizeAtCompileType> Mat_f; + typedef Matrix<double, SizeAtCompileType, SizeAtCompileType> Mat_d; + typedef Matrix<std::complex<float>, SizeAtCompileType, SizeAtCompileType> Mat_cf; + typedef Matrix<std::complex<double>, SizeAtCompileType, SizeAtCompileType> Mat_cd; + typedef Matrix<float, SizeAtCompileType, 1> Vec_f; + typedef Matrix<double, SizeAtCompileType, 1> Vec_d; + typedef Matrix<std::complex<float>, SizeAtCompileType, 1> Vec_cf; + typedef Matrix<std::complex<double>, SizeAtCompileType, 1> Vec_cd; + Mat_f mf(size,size); Mat_d md(size,size); Mat_cf mcf(size,size); @@ -54,14 +153,11 @@ template<int SizeAtCompileType> void mixingtypes(int size = SizeAtCompileType) Vec_cf vcf(size,1); Vec_cd vcd(size,1); - mf+mf; - VERIFY_RAISES_ASSERT(mf+md); - VERIFY_RAISES_ASSERT(mf+mcf); - VERIFY_RAISES_ASSERT(vf=vd); - VERIFY_RAISES_ASSERT(vf+=vd); - VERIFY_RAISES_ASSERT(mcd=md); mf*mf; + // FIXME shall we discard those products ? + // 1) currently they work only if SizeAtCompileType is small enough + // 2) in case we vectorize complexes this might be difficult to still allow that md*mcd; mcd*md; mf*vcf; @@ -69,20 +165,19 @@ template<int SizeAtCompileType> void mixingtypes(int size = SizeAtCompileType) mcf *= mf; vcd = md*vcd; vcf = mcf*vf; - VERIFY_RAISES_ASSERT(mf*md); - VERIFY_RAISES_ASSERT(mcf*mcd); - VERIFY_RAISES_ASSERT(mcf*vcd); +// VERIFY_RAISES_ASSERT(mf*md); // does not even compile +// VERIFY_RAISES_ASSERT(mcf*mcd); // does not even compile +// VERIFY_RAISES_ASSERT(mcf*vcd); // does not even compile VERIFY_RAISES_ASSERT(vcf = mf*vf); - - vf.dot(vf); - VERIFY_RAISES_ASSERT(vd.dot(vf)); - VERIFY_RAISES_ASSERT(vcf.dot(vf)); // yeah eventually we should allow this but i'm too lazy to make that change now in Dot.h -} // especially as that might be rewritten as cwise product .sum() which would make that automatic. +} void test_mixingtypes() { // check that our operator new is indeed called: - CALL_SUBTEST(mixingtypes<3>(3)); - CALL_SUBTEST(mixingtypes<4>(4)); + CALL_SUBTEST(mixingtypes<3>()); + CALL_SUBTEST(mixingtypes<4>()); CALL_SUBTEST(mixingtypes<Dynamic>(20)); + + CALL_SUBTEST(mixingtypes_small<4>()); + CALL_SUBTEST(mixingtypes_large(20)); } diff --git a/test/product.h b/test/product.h index 157f6262b..40773ad90 100644 --- a/test/product.h +++ b/test/product.h @@ -81,7 +81,7 @@ template<typename MatrixType> void product(const MatrixType& m) m3 = m1; m3 *= m1.transpose() * m2; VERIFY_IS_APPROX(m3, m1 * (m1.transpose()*m2)); - VERIFY_IS_APPROX(m3, m1.lazy() * (m1.transpose()*m2)); + VERIFY_IS_APPROX(m3, m1 * (m1.transpose()*m2)); // continue testing Product.h: distributivity VERIFY_IS_APPROX(square*(m1 + m2), square*m1+square*m2); @@ -109,26 +109,26 @@ template<typename MatrixType> void product(const MatrixType& m) // test optimized operator+= path res = square; - res += (m1 * m2.transpose()).lazy(); + res.noalias() += m1 * m2.transpose(); VERIFY_IS_APPROX(res, square + m1 * m2.transpose()); if (NumTraits<Scalar>::HasFloatingPoint && std::min(rows,cols)>1) { VERIFY(areNotApprox(res,square + m2 * m1.transpose())); } vcres = vc2; - vcres += (m1.transpose() * v1).lazy(); + vcres.noalias() += m1.transpose() * v1; VERIFY_IS_APPROX(vcres, vc2 + m1.transpose() * v1); // test optimized operator-= path res = square; - res -= (m1 * m2.transpose()).lazy(); + res.noalias() -= m1 * m2.transpose(); VERIFY_IS_APPROX(res, square - (m1 * m2.transpose())); if (NumTraits<Scalar>::HasFloatingPoint && std::min(rows,cols)>1) { VERIFY(areNotApprox(res,square - m2 * m1.transpose())); } vcres = vc2; - vcres -= (m1.transpose() * v1).lazy(); + vcres.noalias() -= m1.transpose() * v1; VERIFY_IS_APPROX(vcres, vc2 - m1.transpose() * v1); tm1 = m1; @@ -145,7 +145,7 @@ template<typename MatrixType> void product(const MatrixType& m) VERIFY_IS_APPROX(res, m1 * m2.transpose()); res2 = square2; - res2 += (m1.transpose() * m2).lazy(); + res2.noalias() += m1.transpose() * m2; VERIFY_IS_APPROX(res2, square2 + m1.transpose() * m2); if (NumTraits<Scalar>::HasFloatingPoint && std::min(rows,cols)>1) { diff --git a/test/product_extra.cpp b/test/product_extra.cpp index 213dbced6..3ad99fc7a 100644 --- a/test/product_extra.cpp +++ b/test/product_extra.cpp @@ -59,16 +59,13 @@ template<typename MatrixType> void product_extra(const MatrixType& m) // r0 = ei_random<int>(0,rows/2-1), // r1 = ei_random<int>(rows/2,rows); - // all the expressions in this test should be compiled as a single matrix product - // TODO: add internal checks to verify that - - VERIFY_IS_APPROX(m3 = (m1 * m2.adjoint()).lazy(), m1 * m2.adjoint().eval()); - VERIFY_IS_APPROX(m3 = (m1.adjoint() * square.adjoint()).lazy(), m1.adjoint().eval() * square.adjoint().eval()); - VERIFY_IS_APPROX(m3 = (m1.adjoint() * m2).lazy(), m1.adjoint().eval() * m2); - VERIFY_IS_APPROX(m3 = ((s1 * m1.adjoint()) * m2).lazy(), (s1 * m1.adjoint()).eval() * m2); - VERIFY_IS_APPROX(m3 = ((- m1.adjoint() * s1) * (s3 * m2)).lazy(), (- m1.adjoint() * s1).eval() * (s3 * m2).eval()); - VERIFY_IS_APPROX(m3 = ((s2 * m1.adjoint() * s1) * m2).lazy(), (s2 * m1.adjoint() * s1).eval() * m2); - VERIFY_IS_APPROX(m3 = ((-m1*s2) * s1*m2.adjoint()).lazy(), (-m1*s2).eval() * (s1*m2.adjoint()).eval()); + VERIFY_IS_APPROX(m3.noalias() = m1 * m2.adjoint(), m1 * m2.adjoint().eval()); + 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() = (- 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()); // a very tricky case where a scale factor has to be automatically conjugated: VERIFY_IS_APPROX( m1.adjoint() * (s1*m2).conjugate(), (m1.adjoint()).eval() * ((s1*m2).conjugate()).eval()); @@ -76,7 +73,6 @@ template<typename MatrixType> void product_extra(const MatrixType& m) // test all possible conjugate combinations for the four matrix-vector product cases: -// std::cerr << "a\n"; VERIFY_IS_APPROX((-m1.conjugate() * s2) * (s1 * vc2), (-m1.conjugate()*s2).eval() * (s1 * vc2).eval()); VERIFY_IS_APPROX((-m1 * s2) * (s1 * vc2.conjugate()), @@ -84,7 +80,6 @@ template<typename MatrixType> void product_extra(const MatrixType& m) VERIFY_IS_APPROX((-m1.conjugate() * s2) * (s1 * vc2.conjugate()), (-m1.conjugate()*s2).eval() * (s1 * vc2.conjugate()).eval()); -// std::cerr << "b\n"; VERIFY_IS_APPROX((s1 * vc2.transpose()) * (-m1.adjoint() * s2), (s1 * vc2.transpose()).eval() * (-m1.adjoint()*s2).eval()); VERIFY_IS_APPROX((s1 * vc2.adjoint()) * (-m1.transpose() * s2), @@ -92,7 +87,6 @@ template<typename MatrixType> void product_extra(const MatrixType& m) VERIFY_IS_APPROX((s1 * vc2.adjoint()) * (-m1.adjoint() * s2), (s1 * vc2.adjoint()).eval() * (-m1.adjoint()*s2).eval()); -// std::cerr << "c\n"; VERIFY_IS_APPROX((-m1.adjoint() * s2) * (s1 * v1.transpose()), (-m1.adjoint()*s2).eval() * (s1 * v1.transpose()).eval()); VERIFY_IS_APPROX((-m1.transpose() * s2) * (s1 * v1.adjoint()), @@ -100,7 +94,6 @@ template<typename MatrixType> void product_extra(const MatrixType& m) VERIFY_IS_APPROX((-m1.adjoint() * s2) * (s1 * v1.adjoint()), (-m1.adjoint()*s2).eval() * (s1 * v1.adjoint()).eval()); -// std::cerr << "d\n"; VERIFY_IS_APPROX((s1 * v1) * (-m1.conjugate() * s2), (s1 * v1).eval() * (-m1.conjugate()*s2).eval()); VERIFY_IS_APPROX((s1 * v1.conjugate()) * (-m1 * s2), @@ -111,13 +104,24 @@ template<typename MatrixType> void product_extra(const MatrixType& m) VERIFY_IS_APPROX((-m1.adjoint() * s2) * (s1 * v1.adjoint()), (-m1.adjoint()*s2).eval() * (s1 * v1.adjoint()).eval()); + // test the vector-matrix product with non aligned starts + int i = ei_random<int>(0,m1.rows()-2); + int j = ei_random<int>(0,m1.cols()-2); + int r = ei_random<int>(1,m1.rows()-i); + int c = ei_random<int>(1,m1.cols()-j); + int i2 = ei_random<int>(0,m1.rows()-1); + int j2 = ei_random<int>(0,m1.cols()-1); + + VERIFY_IS_APPROX(m1.col(j2).adjoint() * m1.block(0,j,m1.rows(),c), m1.col(j2).adjoint().eval() * m1.block(0,j,m1.rows(),c).eval()); + VERIFY_IS_APPROX(m1.block(i,0,r,m1.cols()) * m1.row(i2).adjoint(), m1.block(i,0,r,m1.cols()).eval() * m1.row(i2).adjoint().eval()); + } void test_product_extra() { for(int i = 0; i < g_repeat; i++) { - CALL_SUBTEST( product_extra(MatrixXf(ei_random<int>(1,320), ei_random<int>(1,320))) ); + CALL_SUBTEST( product_extra(MatrixXf(ei_random<int>(2,320), ei_random<int>(2,320))) ); CALL_SUBTEST( product_extra(MatrixXcf(ei_random<int>(50,50), ei_random<int>(50,50))) ); - CALL_SUBTEST( product_extra(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(ei_random<int>(1,50), ei_random<int>(1,50))) ); + CALL_SUBTEST( product_extra(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(ei_random<int>(2,50), ei_random<int>(2,50))) ); } } diff --git a/test/product_notemporary.cpp b/test/product_notemporary.cpp index f6e8f4101..1a3d65291 100644 --- a/test/product_notemporary.cpp +++ b/test/product_notemporary.cpp @@ -69,53 +69,47 @@ template<typename MatrixType> void product_notemporary(const MatrixType& m) r1 = ei_random<int>(8,rows-r0); VERIFY_EVALUATION_COUNT( m3 = (m1 * m2.adjoint()), 1); - VERIFY_EVALUATION_COUNT( m3 = (m1 * m2.adjoint()).lazy(), 0); - VERIFY_EVALUATION_COUNT( m3.noalias() = m1 * m2.adjoint(), 0); - // NOTE in this case the slow product is used: - // FIXME: VERIFY_EVALUATION_COUNT( m3.noalias() = s1 * (m1 * m2.transpose()), 0); - VERIFY_EVALUATION_COUNT( m3 = (s1 * m1 * s2 * m2.adjoint()).lazy(), 0); - VERIFY_EVALUATION_COUNT( m3 = (s1 * m1 * s2 * (m1*s3+m2*s2).adjoint()).lazy(), 1); - VERIFY_EVALUATION_COUNT( m3 = ((s1 * m1).adjoint() * s2 * m2).lazy(), 0); - VERIFY_EVALUATION_COUNT( m3 += (s1 * (-m1*s3).adjoint() * (s2 * m2 * s3)).lazy(), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() = s1 * m1 * s2 * m2.adjoint(), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() = s1 * m1 * s2 * (m1*s3+m2*s2).adjoint(), 1); + VERIFY_EVALUATION_COUNT( m3.noalias() = (s1 * m1).adjoint() * s2 * m2, 0); VERIFY_EVALUATION_COUNT( m3.noalias() += s1 * (-m1*s3).adjoint() * (s2 * m2 * s3), 0); - VERIFY_EVALUATION_COUNT( m3 -= (s1 * (m1.transpose() * m2)).lazy(), 0); VERIFY_EVALUATION_COUNT( m3.noalias() -= s1 * (m1.transpose() * m2), 0); - VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1) += (-m1.block(r0,c0,r1,c1) * (s2*m2.block(r0,c0,r1,c1)).adjoint()).lazy() ), 0); - VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1) -= (s1 * m1.block(r0,c0,r1,c1) * m2.block(c0,r0,c1,r1)).lazy() ), 0); + VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1).noalias() += -m1.block(r0,c0,r1,c1) * (s2*m2.block(r0,c0,r1,c1)).adjoint() ), 0); + VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1).noalias() -= s1 * m1.block(r0,c0,r1,c1) * m2.block(c0,r0,c1,r1) ), 0); // NOTE this is because the Block expression is not handled yet by our expression analyser - VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1) = (s1 * m1.block(r0,c0,r1,c1) * (s1*m2).block(c0,r0,c1,r1)).lazy() ), 1); + VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1).noalias() = s1 * m1.block(r0,c0,r1,c1) * (s1*m2).block(c0,r0,c1,r1) ), 1); - VERIFY_EVALUATION_COUNT( m3 -= ((s1 * m1).template triangularView<LowerTriangular>() * m2).lazy(), 0); - VERIFY_EVALUATION_COUNT( rm3 = ((s1 * m1.adjoint()).template triangularView<UpperTriangular>() * (m2+m2)).lazy(), 1); - VERIFY_EVALUATION_COUNT( rm3 = ((s1 * m1.adjoint()).template triangularView<UnitUpperTriangular>() * m2.adjoint()).lazy(), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() -= (s1 * m1).template triangularView<LowerTriangular>() * m2, 0); + VERIFY_EVALUATION_COUNT( rm3.noalias() = (s1 * m1.adjoint()).template triangularView<UpperTriangular>() * (m2+m2), 1); + VERIFY_EVALUATION_COUNT( rm3.noalias() = (s1 * m1.adjoint()).template triangularView<UnitUpperTriangular>() * m2.adjoint(), 0); - VERIFY_EVALUATION_COUNT( rm3.col(c0) = ((s1 * m1.adjoint()).template triangularView<UnitUpperTriangular>() * (s2*m2.row(c0)).adjoint()).lazy(), 0); + VERIFY_EVALUATION_COUNT( rm3.col(c0).noalias() = (s1 * m1.adjoint()).template triangularView<UnitUpperTriangular>() * (s2*m2.row(c0)).adjoint(), 0); VERIFY_EVALUATION_COUNT( m1.template triangularView<LowerTriangular>().solveInPlace(m3), 0); VERIFY_EVALUATION_COUNT( m1.adjoint().template triangularView<LowerTriangular>().solveInPlace(m3.transpose()), 0); - VERIFY_EVALUATION_COUNT( m3 -= ((s1 * m1).adjoint().template selfadjointView<LowerTriangular>() * (-m2*s3).adjoint()).lazy(), 0); - VERIFY_EVALUATION_COUNT( m3 = (s2 * m2.adjoint() * (s1 * m1.adjoint()).template selfadjointView<UpperTriangular>()).lazy(), 0); - VERIFY_EVALUATION_COUNT( rm3 = ((s1 * m1.adjoint()).template selfadjointView<LowerTriangular>() * m2.adjoint()).lazy(), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() -= (s1 * m1).adjoint().template selfadjointView<LowerTriangular>() * (-m2*s3).adjoint(), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() = s2 * m2.adjoint() * (s1 * m1.adjoint()).template selfadjointView<UpperTriangular>(), 0); + VERIFY_EVALUATION_COUNT( rm3.noalias() = (s1 * m1.adjoint()).template selfadjointView<LowerTriangular>() * m2.adjoint(), 0); - VERIFY_EVALUATION_COUNT( m3.col(c0) = ((s1 * m1).adjoint().template selfadjointView<LowerTriangular>() * (-m2.row(c0)*s3).adjoint()).lazy(), 0); - VERIFY_EVALUATION_COUNT( m3.col(c0) -= ((s1 * m1).adjoint().template selfadjointView<UpperTriangular>() * (-m2.row(c0)*s3).adjoint()).lazy(), 0); + VERIFY_EVALUATION_COUNT( m3.col(c0).noalias() = (s1 * m1).adjoint().template selfadjointView<LowerTriangular>() * (-m2.row(c0)*s3).adjoint(), 0); + VERIFY_EVALUATION_COUNT( m3.col(c0).noalias() -= (s1 * m1).adjoint().template selfadjointView<UpperTriangular>() * (-m2.row(c0)*s3).adjoint(), 0); - VERIFY_EVALUATION_COUNT( m3.block(r0,c0,r1,c1) += ((m1.block(r0,r0,r1,r1).template selfadjointView<UpperTriangular>() * (s1*m2.block(r0,c0,r1,c1)) )).lazy(), 0); - VERIFY_EVALUATION_COUNT( m3.block(r0,c0,r1,c1) = ((m1.block(r0,r0,r1,r1).template selfadjointView<UpperTriangular>() * m2.block(r0,c0,r1,c1) )).lazy(), 0); + VERIFY_EVALUATION_COUNT( m3.block(r0,c0,r1,c1).noalias() += m1.block(r0,r0,r1,r1).template selfadjointView<UpperTriangular>() * (s1*m2.block(r0,c0,r1,c1)), 0); + VERIFY_EVALUATION_COUNT( m3.block(r0,c0,r1,c1).noalias() = m1.block(r0,r0,r1,r1).template selfadjointView<UpperTriangular>() * m2.block(r0,c0,r1,c1), 0); VERIFY_EVALUATION_COUNT( m3.template selfadjointView<LowerTriangular>().rankUpdate(m2.adjoint()), 0); m3.resize(1,1); - VERIFY_EVALUATION_COUNT( m3 = ((m1.block(r0,r0,r1,r1).template selfadjointView<LowerTriangular>() * m2.block(r0,c0,r1,c1) )).lazy(), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() = m1.block(r0,r0,r1,r1).template selfadjointView<LowerTriangular>() * m2.block(r0,c0,r1,c1), 0); m3.resize(1,1); - VERIFY_EVALUATION_COUNT( m3 = ((m1.block(r0,r0,r1,r1).template triangularView<UnitUpperTriangular>() * m2.block(r0,c0,r1,c1) )).lazy(), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() = m1.block(r0,r0,r1,r1).template triangularView<UnitUpperTriangular>() * m2.block(r0,c0,r1,c1), 0); } void test_product_notemporary() diff --git a/test/product_symm.cpp b/test/product_symm.cpp index 1300928a2..cf0299c64 100644 --- a/test/product_symm.cpp +++ b/test/product_symm.cpp @@ -96,7 +96,7 @@ template<typename Scalar, int Size, int OtherSize> void symm(int size = Size, in m2 = m1.template triangularView<UpperTriangular>(); rhs13 = rhs12; - VERIFY_IS_APPROX(rhs12 += (s1 * ((m2.adjoint()).template selfadjointView<LowerTriangular>() * (s2*rhs3).conjugate())).lazy(), + VERIFY_IS_APPROX(rhs12.noalias() += s1 * ((m2.adjoint()).template selfadjointView<LowerTriangular>() * (s2*rhs3).conjugate()), rhs13 += (s1*m1.adjoint()) * (s2*rhs3).conjugate()); // test matrix * selfadjoint diff --git a/test/product_trsm.cpp b/test/product_trsm.cpp index 9f372df91..756034df9 100644 --- a/test/product_trsm.cpp +++ b/test/product_trsm.cpp @@ -40,8 +40,8 @@ template<typename Scalar> void trsm(int size,int cols) Matrix<Scalar,Dynamic,Dynamic,ColMajor> cmRhs(size,cols), ref(size,cols); Matrix<Scalar,Dynamic,Dynamic,RowMajor> rmRhs(size,cols); - cmLhs.setRandom(); cmLhs *= 0.1; cmLhs.diagonal().cwise() += 1; - rmLhs.setRandom(); rmLhs *= 0.1; rmLhs.diagonal().cwise() += 1; + cmLhs.setRandom(); cmLhs *= static_cast<RealScalar>(0.1); cmLhs.diagonal().cwise() += static_cast<RealScalar>(1); + rmLhs.setRandom(); rmLhs *= static_cast<RealScalar>(0.1); rmLhs.diagonal().cwise() += static_cast<RealScalar>(1); VERIFY_TRSM(cmLhs.conjugate().template triangularView<LowerTriangular>(), cmRhs); VERIFY_TRSM(cmLhs .template triangularView<UpperTriangular>(), cmRhs); diff --git a/test/stable_norm.cpp b/test/stable_norm.cpp new file mode 100644 index 000000000..b8fbf5271 --- /dev/null +++ b/test/stable_norm.cpp @@ -0,0 +1,79 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#include "main.h" + +template<typename MatrixType> void stable_norm(const MatrixType& m) +{ + /* this test covers the following files: + StableNorm.h + */ + + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + + int rows = m.rows(); + int cols = m.cols(); + + Scalar big = ei_random<Scalar>() * std::numeric_limits<RealScalar>::max() * RealScalar(1e-4); + Scalar small = static_cast<RealScalar>(1)/big; + + MatrixType vzero = MatrixType::Zero(rows, cols), + vrand = MatrixType::Random(rows, cols), + vbig(rows, cols), + vsmall(rows,cols); + + vbig.fill(big); + vsmall.fill(small); + + VERIFY_IS_MUCH_SMALLER_THAN(vzero.norm(), static_cast<RealScalar>(1)); + VERIFY_IS_APPROX(vrand.stableNorm(), vrand.norm()); + VERIFY_IS_APPROX(vrand.blueNorm(), vrand.norm()); + VERIFY_IS_APPROX(vrand.hypotNorm(), vrand.norm()); + + RealScalar size = static_cast<RealScalar>(m.size()); + + // test overflow + VERIFY_IS_NOT_APPROX(static_cast<Scalar>(vbig.norm()), ei_sqrt(size)*big); // here the default norm must fail + VERIFY_IS_APPROX(static_cast<Scalar>(vbig.stableNorm()), ei_sqrt(size)*big); + VERIFY_IS_APPROX(static_cast<Scalar>(vbig.blueNorm()), ei_sqrt(size)*big); + VERIFY_IS_APPROX(static_cast<Scalar>(vbig.hypotNorm()), ei_sqrt(size)*big); + + // test underflow + VERIFY_IS_NOT_APPROX(static_cast<Scalar>(vsmall.norm()), ei_sqrt(size)*small); // here the default norm must fail + VERIFY_IS_APPROX(static_cast<Scalar>(vsmall.stableNorm()), ei_sqrt(size)*small); + VERIFY_IS_APPROX(static_cast<Scalar>(vsmall.blueNorm()), ei_sqrt(size)*small); + VERIFY_IS_APPROX(static_cast<Scalar>(vsmall.hypotNorm()), ei_sqrt(size)*small); +} + +void test_stable_norm() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST( stable_norm(Matrix<float, 1, 1>()) ); + CALL_SUBTEST( stable_norm(Vector4d()) ); + CALL_SUBTEST( stable_norm(VectorXd(ei_random<int>(10,2000))) ); + CALL_SUBTEST( stable_norm(VectorXf(ei_random<int>(10,2000))) ); + CALL_SUBTEST( stable_norm(VectorXcd(ei_random<int>(10,2000))) ); + } +} diff --git a/test/submatrices.cpp b/test/submatrices.cpp index a819cadc2..6fe86c281 100644 --- a/test/submatrices.cpp +++ b/test/submatrices.cpp @@ -170,6 +170,48 @@ template<typename MatrixType> void submatrices(const MatrixType& m) VERIFY(ei_real(ones.row(r1).dot(ones.row(r2))) == RealScalar(cols)); } + +template<typename MatrixType> +void compare_using_data_and_stride(const MatrixType& m) +{ + int rows = m.rows(); + int cols = m.cols(); + int size = m.size(); + int stride = m.stride(); + const typename MatrixType::Scalar* data = m.data(); + + for(int j=0;j<cols;++j) + for(int i=0;i<rows;++i) + VERIFY_IS_APPROX(m.coeff(i,j), data[(MatrixType::Flags&RowMajorBit) ? i*stride+j : j*stride + i]); + + if(MatrixType::IsVectorAtCompileTime) + { + VERIFY_IS_APPROX(stride, int((&m.coeff(1))-(&m.coeff(0)))); + for (int i=0;i<size;++i) + VERIFY_IS_APPROX(m.coeff(i), data[i*stride]); + } +} + +template<typename MatrixType> +void data_and_stride(const MatrixType& m) +{ + int rows = m.rows(); + int cols = m.cols(); + + int r1 = ei_random<int>(0,rows-1); + int r2 = ei_random<int>(r1,rows-1); + int c1 = ei_random<int>(0,cols-1); + int c2 = ei_random<int>(c1,cols-1); + + MatrixType m1 = MatrixType::Random(rows, cols); + compare_using_data_and_stride(m1.block(r1, c1, r2-r1+1, c2-c1+1)); + compare_using_data_and_stride(m1.transpose().block(c1, r1, c2-c1+1, r2-r1+1)); + compare_using_data_and_stride(m1.row(r1)); + compare_using_data_and_stride(m1.col(c1)); + compare_using_data_and_stride(m1.row(r1).transpose()); + compare_using_data_and_stride(m1.col(c1).transpose()); +} + void test_submatrices() { for(int i = 0; i < g_repeat; i++) { @@ -179,5 +221,8 @@ void test_submatrices() CALL_SUBTEST( submatrices(MatrixXi(8, 12)) ); CALL_SUBTEST( submatrices(MatrixXcd(20, 20)) ); CALL_SUBTEST( submatrices(MatrixXf(20, 20)) ); + + CALL_SUBTEST( data_and_stride(MatrixXf(ei_random(5,50), ei_random(5,50))) ); + CALL_SUBTEST( data_and_stride(Matrix<int,Dynamic,Dynamic,RowMajor>(ei_random(5,50), ei_random(5,50))) ); } } diff --git a/test/svd.cpp b/test/svd.cpp index 2ccd94764..e6a32bd3f 100644 --- a/test/svd.cpp +++ b/test/svd.cpp @@ -41,15 +41,11 @@ template<typename MatrixType> void svd(const MatrixType& m) Matrix<Scalar, MatrixType::RowsAtCompileTime, 1>::Random(rows,1); Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> x(cols,1), x2(cols,1); - RealScalar largerEps = test_precision<RealScalar>(); - if (ei_is_same_type<RealScalar,float>::ret) - largerEps = 1e-3f; - { SVD<MatrixType> svd(a); MatrixType sigma = MatrixType::Zero(rows,cols); MatrixType matU = MatrixType::Zero(rows,rows); - sigma.block(0,0,cols,cols) = svd.singularValues().asDiagonal(); + sigma.diagonal() = svd.singularValues(); matU = svd.matrixU(); VERIFY_IS_APPROX(a, matU * sigma * svd.matrixV().transpose()); } diff --git a/test/swap.cpp b/test/swap.cpp new file mode 100644 index 000000000..8b325992c --- /dev/null +++ b/test/swap.cpp @@ -0,0 +1,98 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#define EIGEN_NO_STATIC_ASSERT +#include "main.h" + +template<typename T> +struct other_matrix_type +{ + typedef int type; +}; + +template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols> +struct other_matrix_type<Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> > +{ + typedef Matrix<_Scalar, _Rows, _Cols, _Options^RowMajor, _MaxRows, _MaxCols> type; +}; + +template<typename MatrixType> void swap(const MatrixType& m) +{ + typedef typename other_matrix_type<MatrixType>::type OtherMatrixType; + typedef typename MatrixType::Scalar Scalar; + + ei_assert((!ei_is_same_type<MatrixType,OtherMatrixType>::ret)); + int rows = m.rows(); + int cols = m.cols(); + + // construct 3 matrix guaranteed to be distinct + MatrixType m1 = MatrixType::Random(rows,cols); + MatrixType m2 = MatrixType::Random(rows,cols) + Scalar(100) * MatrixType::Identity(rows,cols); + OtherMatrixType m3 = OtherMatrixType::Random(rows,cols) + Scalar(200) * OtherMatrixType::Identity(rows,cols); + + MatrixType m1_copy = m1; + MatrixType m2_copy = m2; + OtherMatrixType m3_copy = m3; + + // test swapping 2 matrices of same type + m1.swap(m2); + VERIFY_IS_APPROX(m1,m2_copy); + VERIFY_IS_APPROX(m2,m1_copy); + m1 = m1_copy; + m2 = m2_copy; + + // test swapping 2 matrices of different types + m1.swap(m3); + VERIFY_IS_APPROX(m1,m3_copy); + VERIFY_IS_APPROX(m3,m1_copy); + m1 = m1_copy; + m3 = m3_copy; + + // test swapping matrix with expression + m1.swap(m2.block(0,0,rows,cols)); + VERIFY_IS_APPROX(m1,m2_copy); + VERIFY_IS_APPROX(m2,m1_copy); + m1 = m1_copy; + m2 = m2_copy; + + // test swapping two expressions of different types + m1.transpose().swap(m3.transpose()); + VERIFY_IS_APPROX(m1,m3_copy); + VERIFY_IS_APPROX(m3,m1_copy); + m1 = m1_copy; + m3 = m3_copy; + + // test assertion on mismatching size -- matrix case + VERIFY_RAISES_ASSERT(m1.swap(m1.row(0))); + // test assertion on mismatching size -- xpr case + VERIFY_RAISES_ASSERT(m1.row(0).swap(m1)); +} + +void test_swap() +{ + CALL_SUBTEST( swap(Matrix3f()) ); // fixed size, no vectorization + CALL_SUBTEST( swap(Matrix4d()) ); // fixed size, possible vectorization + CALL_SUBTEST( swap(MatrixXd(3,3)) ); // dyn size, no vectorization + CALL_SUBTEST( swap(MatrixXf(30,30)) ); // dyn size, possible vectorization +} diff --git a/test/testsuite.cmake b/test/testsuite.cmake index 6a44ce239..37ee87565 100644 --- a/test/testsuite.cmake +++ b/test/testsuite.cmake @@ -27,6 +27,13 @@ # - EIGEN_WORK_DIR: directory used to download the source files and make the builds # default: folder which contains this script # - EIGEN_CMAKE_ARGS: additional arguments passed to cmake +# - EIGEN_GENERATOR_TYPE: allows to overwrite the generator type +# default: nmake (windows +# See http://www.cmake.org/cmake/help/cmake2.6docs.html#section_Generators for a complete +# list of supported generators. +# - EIGEN_NO_UPDATE: allows to submit dash boards from local repositories +# This might be interesting in case you want to submit dashboards +# including local changes. # - CTEST_SOURCE_DIRECTORY: path to eigen's src (use a new and empty folder, not the one you are working on) # default: <EIGEN_WORK_DIR>/src # - CTEST_BINARY_DIRECTORY: build directory @@ -132,11 +139,11 @@ endif(NOT EIGEN_MODE) ## mandatory variables (the default should be ok in most cases): -if(NOT IGNORE_CVS) +if(NOT EIGEN_NO_UPDATE) SET (CTEST_CVS_COMMAND "hg") SET (CTEST_CVS_CHECKOUT "${CTEST_CVS_COMMAND} clone http://bitbucket.org/eigen/eigen2 \"${CTEST_SOURCE_DIRECTORY}\"") SET(CTEST_BACKUP_AND_RESTORE TRUE) # the backup is CVS related ... -endif(NOT IGNORE_CVS) +endif(NOT EIGEN_NO_UPDATE) # which ctest command to use for running the dashboard SET (CTEST_COMMAND "${EIGEN_CMAKE_DIR}ctest -D ${EIGEN_MODE}") @@ -158,15 +165,24 @@ SET(CTEST_START_WITH_EMPTY_BINARY_DIRECTORY TRUE) # any quotes inside of this string if you use it if(WIN32 AND NOT UNIX) #message(SEND_ERROR "win32") - set(CTEST_CMAKE_COMMAND "${CTEST_CMAKE_COMMAND} -G \"NMake Makefiles\" -DCMAKE_MAKE_PROGRAM=nmake") - SET (CTEST_INITIAL_CACHE " - MAKECOMMAND:STRING=nmake -i - CMAKE_MAKE_PROGRAM:FILEPATH=nmake - CMAKE_GENERATOR:INTERNAL=NMake Makefiles - CMAKE_BUILD_TYPE:STRING=Release - BUILDNAME:STRING=${EIGEN_BUILD_STRING} - SITE:STRING=${EIGEN_SITE} - ") + if(EIGEN_GENERATOR_TYPE) + set(CTEST_CMAKE_COMMAND "${CTEST_CMAKE_COMMAND} -G \"${EIGEN_GENERATOR_TYPE}\"") + SET (CTEST_INITIAL_CACHE " + CMAKE_BUILD_TYPE:STRING=Release + BUILDNAME:STRING=${EIGEN_BUILD_STRING} + SITE:STRING=${EIGEN_SITE} + ") + else(EIGEN_GENERATOR_TYPE) + set(CTEST_CMAKE_COMMAND "${CTEST_CMAKE_COMMAND} -G \"NMake Makefiles\" -DCMAKE_MAKE_PROGRAM=nmake") + SET (CTEST_INITIAL_CACHE " + MAKECOMMAND:STRING=nmake -i + CMAKE_MAKE_PROGRAM:FILEPATH=nmake + CMAKE_GENERATOR:INTERNAL=NMake Makefiles + CMAKE_BUILD_TYPE:STRING=Release + BUILDNAME:STRING=${EIGEN_BUILD_STRING} + SITE:STRING=${EIGEN_SITE} + ") + endif(EIGEN_GENERATOR_TYPE) else(WIN32 AND NOT UNIX) SET (CTEST_INITIAL_CACHE " BUILDNAME:STRING=${EIGEN_BUILD_STRING} diff --git a/test/umeyama.cpp b/test/umeyama.cpp index b6d394883..0999c59c9 100644 --- a/test/umeyama.cpp +++ b/test/umeyama.cpp @@ -127,7 +127,7 @@ void run_test(int dim, int num_elements) MatrixX src = MatrixX::Random(dim+1, num_elements); src.row(dim) = Matrix<Scalar, 1, Dynamic>::Constant(num_elements, Scalar(1)); - MatrixX dst = (cR_t*src).lazy(); + MatrixX dst = cR_t*src; MatrixX cR_t_umeyama = umeyama(src.block(0,0,dim,num_elements), dst.block(0,0,dim,num_elements)); @@ -160,7 +160,7 @@ void run_fixed_size_test(int num_elements) MatrixX src = MatrixX::Random(dim+1, num_elements); src.row(dim) = Matrix<Scalar, 1, Dynamic>::Constant(num_elements, Scalar(1)); - MatrixX dst = (cR_t*src).lazy(); + MatrixX dst = cR_t*src; Block<MatrixX, Dimension, Dynamic> src_block(src,0,0,dim,num_elements); Block<MatrixX, Dimension, Dynamic> dst_block(dst,0,0,dim,num_elements); |