diff options
Diffstat (limited to 'eigen/test/eigen2')
53 files changed, 6445 insertions, 0 deletions
diff --git a/eigen/test/eigen2/CMakeLists.txt b/eigen/test/eigen2/CMakeLists.txt new file mode 100644 index 0000000..9615a60 --- /dev/null +++ b/eigen/test/eigen2/CMakeLists.txt @@ -0,0 +1,61 @@ +add_custom_target(eigen2_buildtests) +add_custom_target(eigen2_check COMMAND "ctest -R eigen2") +add_dependencies(eigen2_check eigen2_buildtests) +add_dependencies(buildtests eigen2_buildtests) + +add_definitions("-DEIGEN2_SUPPORT_STAGE10_FULL_EIGEN2_API") +add_definitions("-DEIGEN_NO_EIGEN2_DEPRECATED_WARNING") + +ei_add_test(eigen2_meta) +ei_add_test(eigen2_sizeof) +ei_add_test(eigen2_dynalloc) +ei_add_test(eigen2_nomalloc) +#ei_add_test(eigen2_first_aligned) +ei_add_test(eigen2_mixingtypes) +#ei_add_test(eigen2_packetmath) +ei_add_test(eigen2_unalignedassert) +#ei_add_test(eigen2_vectorization_logic) +ei_add_test(eigen2_basicstuff) +ei_add_test(eigen2_linearstructure) +ei_add_test(eigen2_cwiseop) +ei_add_test(eigen2_sum) +ei_add_test(eigen2_product_small) +ei_add_test(eigen2_product_large ${EI_OFLAG}) +ei_add_test(eigen2_adjoint) +ei_add_test(eigen2_submatrices) +ei_add_test(eigen2_miscmatrices) +ei_add_test(eigen2_commainitializer) +ei_add_test(eigen2_smallvectors) +ei_add_test(eigen2_map) +ei_add_test(eigen2_array) +ei_add_test(eigen2_triangular) +ei_add_test(eigen2_cholesky " " "${GSL_LIBRARIES}") +ei_add_test(eigen2_lu ${EI_OFLAG}) +ei_add_test(eigen2_determinant ${EI_OFLAG}) +ei_add_test(eigen2_inverse) +ei_add_test(eigen2_qr) +ei_add_test(eigen2_eigensolver " " "${GSL_LIBRARIES}") +ei_add_test(eigen2_svd) +ei_add_test(eigen2_geometry) +ei_add_test(eigen2_geometry_with_eigen2_prefix) +ei_add_test(eigen2_hyperplane) +ei_add_test(eigen2_parametrizedline) +ei_add_test(eigen2_alignedbox) +ei_add_test(eigen2_regression) +ei_add_test(eigen2_stdvector) +ei_add_test(eigen2_newstdvector) +if(QT4_FOUND) + ei_add_test(eigen2_qtvector " " "${QT_QTCORE_LIBRARY}") +endif(QT4_FOUND) +# no support for eigen2 sparse module +# if(NOT EIGEN_DEFAULT_TO_ROW_MAJOR) +# ei_add_test(eigen2_sparse_vector) +# ei_add_test(eigen2_sparse_basic) +# ei_add_test(eigen2_sparse_solvers " " "${SPARSE_LIBS}") +# ei_add_test(eigen2_sparse_product) +# endif() +ei_add_test(eigen2_swap) +ei_add_test(eigen2_visitor) +ei_add_test(eigen2_bug_132) + +ei_add_test(eigen2_prec_inverse_4x4 ${EI_OFLAG}) diff --git a/eigen/test/eigen2/eigen2_adjoint.cpp b/eigen/test/eigen2/eigen2_adjoint.cpp new file mode 100644 index 0000000..c0f8114 --- /dev/null +++ b/eigen/test/eigen2/eigen2_adjoint.cpp @@ -0,0 +1,99 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" + +template<typename MatrixType> void adjoint(const MatrixType& m) +{ + /* this test covers the following files: + Transpose.h Conjugate.h Dot.h + */ + + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; + typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType; + int rows = m.rows(); + int cols = m.cols(); + + RealScalar largerEps = test_precision<RealScalar>(); + if (ei_is_same_type<RealScalar,float>::ret) + largerEps = RealScalar(1e-3f); + + MatrixType m1 = MatrixType::Random(rows, cols), + m2 = MatrixType::Random(rows, cols), + m3(rows, cols), + square = SquareMatrixType::Random(rows, rows); + VectorType v1 = VectorType::Random(rows), + v2 = VectorType::Random(rows), + v3 = VectorType::Random(rows), + vzero = VectorType::Zero(rows); + + Scalar s1 = ei_random<Scalar>(), + s2 = ei_random<Scalar>(); + + // check basic compatibility of adjoint, transpose, conjugate + VERIFY_IS_APPROX(m1.transpose().conjugate().adjoint(), m1); + VERIFY_IS_APPROX(m1.adjoint().conjugate().transpose(), m1); + + // check multiplicative behavior + VERIFY_IS_APPROX((m1.adjoint() * m2).adjoint(), m2.adjoint() * m1); + VERIFY_IS_APPROX((s1 * m1).adjoint(), ei_conj(s1) * m1.adjoint()); + + // check basic properties of dot, norm, norm2 + typedef typename NumTraits<Scalar>::Real RealScalar; + VERIFY(ei_isApprox((s1 * v1 + s2 * v2).eigen2_dot(v3), s1 * v1.eigen2_dot(v3) + s2 * v2.eigen2_dot(v3), largerEps)); + VERIFY(ei_isApprox(v3.eigen2_dot(s1 * v1 + s2 * v2), ei_conj(s1)*v3.eigen2_dot(v1)+ei_conj(s2)*v3.eigen2_dot(v2), largerEps)); + VERIFY_IS_APPROX(ei_conj(v1.eigen2_dot(v2)), v2.eigen2_dot(v1)); + VERIFY_IS_APPROX(ei_real(v1.eigen2_dot(v1)), v1.squaredNorm()); + if(NumTraits<Scalar>::HasFloatingPoint) + VERIFY_IS_APPROX(v1.squaredNorm(), v1.norm() * v1.norm()); + VERIFY_IS_MUCH_SMALLER_THAN(ei_abs(vzero.eigen2_dot(v1)), static_cast<RealScalar>(1)); + if(NumTraits<Scalar>::HasFloatingPoint) + VERIFY_IS_MUCH_SMALLER_THAN(vzero.norm(), static_cast<RealScalar>(1)); + + // check compatibility of dot and adjoint + VERIFY(ei_isApprox(v1.eigen2_dot(square * v2), (square.adjoint() * v1).eigen2_dot(v2), largerEps)); + + // like in testBasicStuff, test operator() to check const-qualification + int r = ei_random<int>(0, rows-1), + c = ei_random<int>(0, cols-1); + VERIFY_IS_APPROX(m1.conjugate()(r,c), ei_conj(m1(r,c))); + VERIFY_IS_APPROX(m1.adjoint()(c,r), ei_conj(m1(r,c))); + + if(NumTraits<Scalar>::HasFloatingPoint) + { + // 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)); + } + + // check inplace transpose + m3 = m1; + m3.transposeInPlace(); + VERIFY_IS_APPROX(m3,m1.transpose()); + m3.transposeInPlace(); + VERIFY_IS_APPROX(m3,m1); + +} + +void test_eigen2_adjoint() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( adjoint(Matrix<float, 1, 1>()) ); + CALL_SUBTEST_2( adjoint(Matrix3d()) ); + CALL_SUBTEST_3( adjoint(Matrix4f()) ); + CALL_SUBTEST_4( adjoint(MatrixXcf(4, 4)) ); + CALL_SUBTEST_5( adjoint(MatrixXi(8, 12)) ); + CALL_SUBTEST_6( adjoint(MatrixXf(21, 21)) ); + } + // test a large matrix only once + CALL_SUBTEST_7( adjoint(Matrix<float, 100, 100>()) ); +} + diff --git a/eigen/test/eigen2/eigen2_alignedbox.cpp b/eigen/test/eigen2/eigen2_alignedbox.cpp new file mode 100644 index 0000000..35043b9 --- /dev/null +++ b/eigen/test/eigen2/eigen2_alignedbox.cpp @@ -0,0 +1,60 @@ +// 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 Gael Guennebaud <g.gael@free.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" +#include <Eigen/Geometry> +#include <Eigen/LU> +#include <Eigen/QR> + +template<typename BoxType> void alignedbox(const BoxType& _box) +{ + /* this test covers the following files: + AlignedBox.h + */ + + const int dim = _box.dim(); + typedef typename BoxType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, BoxType::AmbientDimAtCompileTime, 1> VectorType; + + VectorType p0 = VectorType::Random(dim); + VectorType p1 = VectorType::Random(dim); + RealScalar s1 = ei_random<RealScalar>(0,1); + + BoxType b0(dim); + BoxType b1(VectorType::Random(dim),VectorType::Random(dim)); + BoxType b2; + + b0.extend(p0); + b0.extend(p1); + VERIFY(b0.contains(p0*s1+(Scalar(1)-s1)*p1)); + VERIFY(!b0.contains(p0 + (1+s1)*(p1-p0))); + + (b2 = b0).extend(b1); + VERIFY(b2.contains(b0)); + VERIFY(b2.contains(b1)); + VERIFY_IS_APPROX(b2.clamp(b0), b0); + + // casting + const int Dim = BoxType::AmbientDimAtCompileTime; + typedef typename GetDifferentType<Scalar>::type OtherScalar; + AlignedBox<OtherScalar,Dim> hp1f = b0.template cast<OtherScalar>(); + VERIFY_IS_APPROX(hp1f.template cast<Scalar>(),b0); + AlignedBox<Scalar,Dim> hp1d = b0.template cast<Scalar>(); + VERIFY_IS_APPROX(hp1d.template cast<Scalar>(),b0); +} + +void test_eigen2_alignedbox() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( alignedbox(AlignedBox<float,2>()) ); + CALL_SUBTEST_2( alignedbox(AlignedBox<float,3>()) ); + CALL_SUBTEST_3( alignedbox(AlignedBox<double,4>()) ); + } +} diff --git a/eigen/test/eigen2/eigen2_array.cpp b/eigen/test/eigen2/eigen2_array.cpp new file mode 100644 index 0000000..c1ff40c --- /dev/null +++ b/eigen/test/eigen2/eigen2_array.cpp @@ -0,0 +1,142 @@ +// 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 Gael Guennebaud <g.gael@free.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" +#include <Eigen/Array> + +template<typename MatrixType> void array(const MatrixType& m) +{ + /* this test covers the following files: + Array.cpp + */ + + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; + + int rows = m.rows(); + int cols = m.cols(); + + MatrixType m1 = MatrixType::Random(rows, cols), + m2 = MatrixType::Random(rows, cols), + m3(rows, cols); + + Scalar s1 = ei_random<Scalar>(), + s2 = ei_random<Scalar>(); + + // scalar addition + VERIFY_IS_APPROX(m1.cwise() + s1, s1 + m1.cwise()); + VERIFY_IS_APPROX(m1.cwise() + s1, MatrixType::Constant(rows,cols,s1) + m1); + VERIFY_IS_APPROX((m1*Scalar(2)).cwise() - s2, (m1+m1) - MatrixType::Constant(rows,cols,s2) ); + m3 = m1; + m3.cwise() += s2; + VERIFY_IS_APPROX(m3, m1.cwise() + s2); + m3 = m1; + m3.cwise() -= s1; + VERIFY_IS_APPROX(m3, m1.cwise() - s1); + + // reductions + VERIFY_IS_APPROX(m1.colwise().sum().sum(), m1.sum()); + VERIFY_IS_APPROX(m1.rowwise().sum().sum(), m1.sum()); + if (!ei_isApprox(m1.sum(), (m1+m2).sum())) + VERIFY_IS_NOT_APPROX(((m1+m2).rowwise().sum()).sum(), m1.sum()); + VERIFY_IS_APPROX(m1.colwise().sum(), m1.colwise().redux(internal::scalar_sum_op<Scalar>())); +} + +template<typename MatrixType> void comparisons(const MatrixType& m) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; + + int rows = m.rows(); + int cols = m.cols(); + + int r = ei_random<int>(0, rows-1), + c = ei_random<int>(0, cols-1); + + MatrixType m1 = MatrixType::Random(rows, cols), + m2 = MatrixType::Random(rows, cols), + m3(rows, cols); + + VERIFY(((m1.cwise() + Scalar(1)).cwise() > m1).all()); + VERIFY(((m1.cwise() - Scalar(1)).cwise() < m1).all()); + if (rows*cols>1) + { + m3 = m1; + m3(r,c) += 1; + VERIFY(! (m1.cwise() < m3).all() ); + VERIFY(! (m1.cwise() > m3).all() ); + } + + // comparisons to scalar + VERIFY( (m1.cwise() != (m1(r,c)+1) ).any() ); + VERIFY( (m1.cwise() > (m1(r,c)-1) ).any() ); + VERIFY( (m1.cwise() < (m1(r,c)+1) ).any() ); + VERIFY( (m1.cwise() == m1(r,c) ).any() ); + + // test Select + VERIFY_IS_APPROX( (m1.cwise()<m2).select(m1,m2), m1.cwise().min(m2) ); + VERIFY_IS_APPROX( (m1.cwise()>m2).select(m1,m2), m1.cwise().max(m2) ); + Scalar mid = (m1.cwise().abs().minCoeff() + m1.cwise().abs().maxCoeff())/Scalar(2); + for (int j=0; j<cols; ++j) + for (int i=0; i<rows; ++i) + m3(i,j) = ei_abs(m1(i,j))<mid ? 0 : m1(i,j); + VERIFY_IS_APPROX( (m1.cwise().abs().cwise()<MatrixType::Constant(rows,cols,mid)) + .select(MatrixType::Zero(rows,cols),m1), m3); + // shorter versions: + VERIFY_IS_APPROX( (m1.cwise().abs().cwise()<MatrixType::Constant(rows,cols,mid)) + .select(0,m1), m3); + VERIFY_IS_APPROX( (m1.cwise().abs().cwise()>=MatrixType::Constant(rows,cols,mid)) + .select(m1,0), m3); + // even shorter version: + VERIFY_IS_APPROX( (m1.cwise().abs().cwise()<mid).select(0,m1), m3); + + // count + VERIFY(((m1.cwise().abs().cwise()+1).cwise()>RealScalar(0.1)).count() == rows*cols); + VERIFY_IS_APPROX(((m1.cwise().abs().cwise()+1).cwise()>RealScalar(0.1)).colwise().count().template cast<int>(), RowVectorXi::Constant(cols,rows)); + VERIFY_IS_APPROX(((m1.cwise().abs().cwise()+1).cwise()>RealScalar(0.1)).rowwise().count().template cast<int>(), VectorXi::Constant(rows, cols)); +} + +template<typename VectorType> void lpNorm(const VectorType& v) +{ + VectorType u = VectorType::Random(v.size()); + + VERIFY_IS_APPROX(u.template lpNorm<Infinity>(), u.cwise().abs().maxCoeff()); + VERIFY_IS_APPROX(u.template lpNorm<1>(), u.cwise().abs().sum()); + VERIFY_IS_APPROX(u.template lpNorm<2>(), ei_sqrt(u.cwise().abs().cwise().square().sum())); + VERIFY_IS_APPROX(ei_pow(u.template lpNorm<5>(), typename VectorType::RealScalar(5)), u.cwise().abs().cwise().pow(5).sum()); +} + +void test_eigen2_array() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( array(Matrix<float, 1, 1>()) ); + CALL_SUBTEST_2( array(Matrix2f()) ); + CALL_SUBTEST_3( array(Matrix4d()) ); + CALL_SUBTEST_4( array(MatrixXcf(3, 3)) ); + CALL_SUBTEST_5( array(MatrixXf(8, 12)) ); + CALL_SUBTEST_6( array(MatrixXi(8, 12)) ); + } + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( comparisons(Matrix<float, 1, 1>()) ); + CALL_SUBTEST_2( comparisons(Matrix2f()) ); + CALL_SUBTEST_3( comparisons(Matrix4d()) ); + CALL_SUBTEST_5( comparisons(MatrixXf(8, 12)) ); + CALL_SUBTEST_6( comparisons(MatrixXi(8, 12)) ); + } + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( lpNorm(Matrix<float, 1, 1>()) ); + CALL_SUBTEST_2( lpNorm(Vector2f()) ); + CALL_SUBTEST_3( lpNorm(Vector3d()) ); + CALL_SUBTEST_4( lpNorm(Vector4f()) ); + CALL_SUBTEST_5( lpNorm(VectorXf(16)) ); + CALL_SUBTEST_7( lpNorm(VectorXcd(10)) ); + } +} diff --git a/eigen/test/eigen2/eigen2_basicstuff.cpp b/eigen/test/eigen2/eigen2_basicstuff.cpp new file mode 100644 index 0000000..dd2dec1 --- /dev/null +++ b/eigen/test/eigen2/eigen2_basicstuff.cpp @@ -0,0 +1,105 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" + +template<typename MatrixType> void basicStuff(const MatrixType& m) +{ + typedef typename MatrixType::Scalar Scalar; + typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; + + int rows = m.rows(); + int cols = m.cols(); + + // this test relies a lot on Random.h, and there's not much more that we can do + // to test it, hence I consider that we will have tested Random.h + MatrixType m1 = MatrixType::Random(rows, cols), + m2 = MatrixType::Random(rows, cols), + m3(rows, cols), + mzero = MatrixType::Zero(rows, cols), + square = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>::Random(rows, rows); + VectorType v1 = VectorType::Random(rows), + vzero = VectorType::Zero(rows); + + Scalar x = ei_random<Scalar>(); + + int r = ei_random<int>(0, rows-1), + c = ei_random<int>(0, cols-1); + + m1.coeffRef(r,c) = x; + VERIFY_IS_APPROX(x, m1.coeff(r,c)); + m1(r,c) = x; + VERIFY_IS_APPROX(x, m1(r,c)); + v1.coeffRef(r) = x; + VERIFY_IS_APPROX(x, v1.coeff(r)); + v1(r) = x; + VERIFY_IS_APPROX(x, v1(r)); + v1[r] = x; + VERIFY_IS_APPROX(x, v1[r]); + + VERIFY_IS_APPROX( v1, v1); + VERIFY_IS_NOT_APPROX( v1, 2*v1); + VERIFY_IS_MUCH_SMALLER_THAN( vzero, v1); + if(NumTraits<Scalar>::HasFloatingPoint) + VERIFY_IS_MUCH_SMALLER_THAN( vzero, v1.norm()); + VERIFY_IS_NOT_MUCH_SMALLER_THAN(v1, v1); + VERIFY_IS_APPROX( vzero, v1-v1); + VERIFY_IS_APPROX( m1, m1); + VERIFY_IS_NOT_APPROX( m1, 2*m1); + VERIFY_IS_MUCH_SMALLER_THAN( mzero, m1); + VERIFY_IS_NOT_MUCH_SMALLER_THAN(m1, m1); + VERIFY_IS_APPROX( mzero, m1-m1); + + // always test operator() on each read-only expression class, + // in order to check const-qualifiers. + // indeed, if an expression class (here Zero) is meant to be read-only, + // hence has no _write() method, the corresponding MatrixBase method (here zero()) + // should return a const-qualified object so that it is the const-qualified + // operator() that gets called, which in turn calls _read(). + VERIFY_IS_MUCH_SMALLER_THAN(MatrixType::Zero(rows,cols)(r,c), static_cast<Scalar>(1)); + + // now test copying a row-vector into a (column-)vector and conversely. + square.col(r) = square.row(r).eval(); + Matrix<Scalar, 1, MatrixType::RowsAtCompileTime> rv(rows); + Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> cv(rows); + rv = square.row(r); + cv = square.col(r); + VERIFY_IS_APPROX(rv, cv.transpose()); + + if(cols!=1 && rows!=1 && MatrixType::SizeAtCompileTime!=Dynamic) + { + VERIFY_RAISES_ASSERT(m1 = (m2.block(0,0, rows-1, cols-1))); + } + + VERIFY_IS_APPROX(m3 = m1,m1); + 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); + } +} + +void test_eigen2_basicstuff() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( basicStuff(Matrix<float, 1, 1>()) ); + CALL_SUBTEST_2( basicStuff(Matrix4d()) ); + CALL_SUBTEST_3( basicStuff(MatrixXcf(3, 3)) ); + CALL_SUBTEST_4( basicStuff(MatrixXi(8, 12)) ); + CALL_SUBTEST_5( basicStuff(MatrixXcd(20, 20)) ); + CALL_SUBTEST_6( basicStuff(Matrix<float, 100, 100>()) ); + CALL_SUBTEST_7( basicStuff(Matrix<long double,Dynamic,Dynamic>(10,10)) ); + } +} diff --git a/eigen/test/eigen2/eigen2_bug_132.cpp b/eigen/test/eigen2/eigen2_bug_132.cpp new file mode 100644 index 0000000..7fe3610 --- /dev/null +++ b/eigen/test/eigen2/eigen2_bug_132.cpp @@ -0,0 +1,26 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2010 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// 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" + +void test_eigen2_bug_132() { + int size = 100; + MatrixXd A(size, size); + VectorXd b(size), c(size); + { + VectorXd y = A.transpose() * (b-c); // bug 132: infinite recursion in coeffRef + VectorXd z = (b-c).transpose() * A; // bug 132: infinite recursion in coeffRef + } + + // the following ones weren't failing, but let's include them for completeness: + { + VectorXd y = A * (b-c); + VectorXd z = (b-c).transpose() * A.transpose(); + } +} diff --git a/eigen/test/eigen2/eigen2_cholesky.cpp b/eigen/test/eigen2/eigen2_cholesky.cpp new file mode 100644 index 0000000..9c4b6f5 --- /dev/null +++ b/eigen/test/eigen2/eigen2_cholesky.cpp @@ -0,0 +1,113 @@ +// 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 Gael Guennebaud <g.gael@free.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#define EIGEN_NO_ASSERTION_CHECKING +#include "main.h" +#include <Eigen/Cholesky> +#include <Eigen/LU> + +#ifdef HAS_GSL +#include "gsl_helper.h" +#endif + +template<typename MatrixType> void cholesky(const MatrixType& m) +{ + /* this test covers the following files: + LLT.h LDLT.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, MatrixType::RowsAtCompileTime> SquareMatrixType; + typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; + + MatrixType a0 = MatrixType::Random(rows,cols); + VectorType vecB = VectorType::Random(rows), vecX(rows); + MatrixType matB = MatrixType::Random(rows,cols), matX(rows,cols); + SquareMatrixType symm = a0 * a0.adjoint(); + // let's make sure the matrix is not singular or near singular + MatrixType a1 = MatrixType::Random(rows,cols); + symm += a1 * a1.adjoint(); + + #ifdef HAS_GSL + if (ei_is_same_type<RealScalar,double>::ret) + { + typedef GslTraits<Scalar> Gsl; + typename Gsl::Matrix gMatA=0, gSymm=0; + typename Gsl::Vector gVecB=0, gVecX=0; + convert<MatrixType>(symm, gSymm); + convert<MatrixType>(symm, gMatA); + convert<VectorType>(vecB, gVecB); + convert<VectorType>(vecB, gVecX); + Gsl::cholesky(gMatA); + Gsl::cholesky_solve(gMatA, gVecB, gVecX); + VectorType vecX(rows), _vecX, _vecB; + convert(gVecX, _vecX); + symm.llt().solve(vecB, &vecX); + Gsl::prod(gSymm, gVecX, gVecB); + convert(gVecB, _vecB); + // test gsl itself ! + VERIFY_IS_APPROX(vecB, _vecB); + VERIFY_IS_APPROX(vecX, _vecX); + + Gsl::free(gMatA); + Gsl::free(gSymm); + Gsl::free(gVecB); + Gsl::free(gVecX); + } + #endif + + { + LDLT<SquareMatrixType> ldlt(symm); + VERIFY(ldlt.isPositiveDefinite()); + // in eigen3, LDLT is pivoting + //VERIFY_IS_APPROX(symm, ldlt.matrixL() * ldlt.vectorD().asDiagonal() * ldlt.matrixL().adjoint()); + ldlt.solve(vecB, &vecX); + VERIFY_IS_APPROX(symm * vecX, vecB); + ldlt.solve(matB, &matX); + VERIFY_IS_APPROX(symm * matX, matB); + } + + { + LLT<SquareMatrixType> chol(symm); + VERIFY(chol.isPositiveDefinite()); + VERIFY_IS_APPROX(symm, chol.matrixL() * chol.matrixL().adjoint()); + chol.solve(vecB, &vecX); + VERIFY_IS_APPROX(symm * vecX, vecB); + chol.solve(matB, &matX); + VERIFY_IS_APPROX(symm * matX, matB); + } + +#if 0 // cholesky is not rank-revealing anyway + // test isPositiveDefinite on non definite matrix + if (rows>4) + { + SquareMatrixType symm = a0.block(0,0,rows,cols-4) * a0.block(0,0,rows,cols-4).adjoint(); + LLT<SquareMatrixType> chol(symm); + VERIFY(!chol.isPositiveDefinite()); + LDLT<SquareMatrixType> cholnosqrt(symm); + VERIFY(!cholnosqrt.isPositiveDefinite()); + } +#endif +} + +void test_eigen2_cholesky() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( cholesky(Matrix<double,1,1>()) ); + CALL_SUBTEST_2( cholesky(Matrix2d()) ); + CALL_SUBTEST_3( cholesky(Matrix3f()) ); + CALL_SUBTEST_4( cholesky(Matrix4d()) ); + CALL_SUBTEST_5( cholesky(MatrixXcd(7,7)) ); + CALL_SUBTEST_6( cholesky(MatrixXf(17,17)) ); + CALL_SUBTEST_7( cholesky(MatrixXd(33,33)) ); + } +} diff --git a/eigen/test/eigen2/eigen2_commainitializer.cpp b/eigen/test/eigen2/eigen2_commainitializer.cpp new file mode 100644 index 0000000..e0f901e --- /dev/null +++ b/eigen/test/eigen2/eigen2_commainitializer.cpp @@ -0,0 +1,46 @@ +// 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 Gael Guennebaud <g.gael@free.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" + +void test_eigen2_commainitializer() +{ + Matrix3d m3; + Matrix4d m4; + + VERIFY_RAISES_ASSERT( (m3 << 1, 2, 3, 4, 5, 6, 7, 8) ); + + #ifndef _MSC_VER + VERIFY_RAISES_ASSERT( (m3 << 1, 2, 3, 4, 5, 6, 7, 8, 9, 10) ); + #endif + + double data[] = {1, 2, 3, 4, 5, 6, 7, 8, 9}; + Matrix3d ref = Map<Matrix<double,3,3,RowMajor> >(data); + + m3 = Matrix3d::Random(); + m3 << 1, 2, 3, 4, 5, 6, 7, 8, 9; + VERIFY_IS_APPROX(m3, ref ); + + Vector3d vec[3]; + vec[0] << 1, 4, 7; + vec[1] << 2, 5, 8; + vec[2] << 3, 6, 9; + m3 = Matrix3d::Random(); + m3 << vec[0], vec[1], vec[2]; + VERIFY_IS_APPROX(m3, ref); + + vec[0] << 1, 2, 3; + vec[1] << 4, 5, 6; + vec[2] << 7, 8, 9; + m3 = Matrix3d::Random(); + m3 << vec[0].transpose(), + 4, 5, 6, + vec[2].transpose(); + VERIFY_IS_APPROX(m3, ref); +} diff --git a/eigen/test/eigen2/eigen2_cwiseop.cpp b/eigen/test/eigen2/eigen2_cwiseop.cpp new file mode 100644 index 0000000..22e1cc3 --- /dev/null +++ b/eigen/test/eigen2/eigen2_cwiseop.cpp @@ -0,0 +1,155 @@ +// 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 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// 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 <functional> +#include <Eigen/Array> + +using namespace std; + +template<typename Scalar> struct AddIfNull { + const Scalar operator() (const Scalar a, const Scalar b) const {return a<=1e-3 ? b : a;} + enum { Cost = NumTraits<Scalar>::AddCost }; +}; + +template<typename MatrixType> void cwiseops(const MatrixType& m) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; + + int rows = m.rows(); + int cols = m.cols(); + + MatrixType m1 = MatrixType::Random(rows, cols), + m2 = MatrixType::Random(rows, cols), + m3(rows, cols), + m4(rows, cols), + mzero = MatrixType::Zero(rows, cols), + mones = MatrixType::Ones(rows, cols), + identity = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> + ::Identity(rows, rows); + VectorType vzero = VectorType::Zero(rows), + vones = VectorType::Ones(rows), + v3(rows); + + int r = ei_random<int>(0, rows-1), + c = ei_random<int>(0, cols-1); + + Scalar s1 = ei_random<Scalar>(); + + // test Zero, Ones, Constant, and the set* variants + m3 = MatrixType::Constant(rows, cols, s1); + for (int j=0; j<cols; ++j) + for (int i=0; i<rows; ++i) + { + VERIFY_IS_APPROX(mzero(i,j), Scalar(0)); + VERIFY_IS_APPROX(mones(i,j), Scalar(1)); + VERIFY_IS_APPROX(m3(i,j), s1); + } + VERIFY(mzero.isZero()); + VERIFY(mones.isOnes()); + VERIFY(m3.isConstant(s1)); + VERIFY(identity.isIdentity()); + VERIFY_IS_APPROX(m4.setConstant(s1), m3); + VERIFY_IS_APPROX(m4.setConstant(rows,cols,s1), m3); + VERIFY_IS_APPROX(m4.setZero(), mzero); + VERIFY_IS_APPROX(m4.setZero(rows,cols), mzero); + VERIFY_IS_APPROX(m4.setOnes(), mones); + VERIFY_IS_APPROX(m4.setOnes(rows,cols), mones); + m4.fill(s1); + VERIFY_IS_APPROX(m4, m3); + + VERIFY_IS_APPROX(v3.setConstant(rows, s1), VectorType::Constant(rows,s1)); + VERIFY_IS_APPROX(v3.setZero(rows), vzero); + VERIFY_IS_APPROX(v3.setOnes(rows), vones); + + m2 = m2.template binaryExpr<AddIfNull<Scalar> >(mones); + + VERIFY_IS_APPROX(m1.cwise().pow(2), m1.cwise().abs2()); + VERIFY_IS_APPROX(m1.cwise().pow(2), m1.cwise().square()); + VERIFY_IS_APPROX(m1.cwise().pow(3), m1.cwise().cube()); + + VERIFY_IS_APPROX(m1 + mones, m1.cwise()+Scalar(1)); + VERIFY_IS_APPROX(m1 - mones, m1.cwise()-Scalar(1)); + m3 = m1; m3.cwise() += 1; + VERIFY_IS_APPROX(m1 + mones, m3); + m3 = m1; m3.cwise() -= 1; + VERIFY_IS_APPROX(m1 - mones, m3); + + VERIFY_IS_APPROX(m2, m2.cwise() * mones); + VERIFY_IS_APPROX(m1.cwise() * m2, m2.cwise() * m1); + m3 = m1; + m3.cwise() *= m2; + VERIFY_IS_APPROX(m3, m1.cwise() * m2); + + VERIFY_IS_APPROX(mones, m2.cwise()/m2); + if(NumTraits<Scalar>::HasFloatingPoint) + { + 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) ); + VERIFY_IS_APPROX( m1.cwise().min(m1+mones), m1 ); + VERIFY_IS_APPROX( m1.cwise().min(m1-mones), m1-mones ); + + // check max + VERIFY_IS_APPROX( m1.cwise().max(m2), m2.cwise().max(m1) ); + VERIFY_IS_APPROX( m1.cwise().max(m1-mones), m1 ); + VERIFY_IS_APPROX( m1.cwise().max(m1+mones), m1+mones ); + + VERIFY( (m1.cwise() == m1).all() ); + VERIFY( (m1.cwise() != m2).any() ); + VERIFY(!(m1.cwise() == (m1+mones)).any() ); + if (rows*cols>1) + { + m3 = m1; + m3(r,c) += 1; + VERIFY( (m1.cwise() == m3).any() ); + VERIFY( !(m1.cwise() == m3).all() ); + } + VERIFY( (m1.cwise().min(m2).cwise() <= m2).all() ); + VERIFY( (m1.cwise().max(m2).cwise() >= m2).all() ); + VERIFY( (m1.cwise().min(m2).cwise() < (m1+mones)).all() ); + VERIFY( (m1.cwise().max(m2).cwise() > (m1-mones)).all() ); + + VERIFY( (m1.cwise()<m1.unaryExpr(bind2nd(plus<Scalar>(), Scalar(1)))).all() ); + VERIFY( !(m1.cwise()<m1.unaryExpr(bind2nd(minus<Scalar>(), Scalar(1)))).all() ); + VERIFY( !(m1.cwise()>m1.unaryExpr(bind2nd(plus<Scalar>(), Scalar(1)))).any() ); +} + +void test_eigen2_cwiseop() +{ + for(int i = 0; i < g_repeat ; i++) { + CALL_SUBTEST_1( cwiseops(Matrix<float, 1, 1>()) ); + CALL_SUBTEST_2( cwiseops(Matrix4d()) ); + CALL_SUBTEST_3( cwiseops(MatrixXf(3, 3)) ); + CALL_SUBTEST_3( cwiseops(MatrixXf(22, 22)) ); + CALL_SUBTEST_4( cwiseops(MatrixXi(8, 12)) ); + CALL_SUBTEST_5( cwiseops(MatrixXd(20, 20)) ); + } +} diff --git a/eigen/test/eigen2/eigen2_determinant.cpp b/eigen/test/eigen2/eigen2_determinant.cpp new file mode 100644 index 0000000..c7b4ad0 --- /dev/null +++ b/eigen/test/eigen2/eigen2_determinant.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 Benoit Jacob <jacob.benoit.1@gmail.com> +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" +#include <Eigen/LU> + +template<typename MatrixType> void determinant(const MatrixType& m) +{ + /* this test covers the following files: + Determinant.h + */ + int size = m.rows(); + + MatrixType m1(size, size), m2(size, size); + m1.setRandom(); + m2.setRandom(); + typedef typename MatrixType::Scalar Scalar; + Scalar x = ei_random<Scalar>(); + VERIFY_IS_APPROX(MatrixType::Identity(size, size).determinant(), Scalar(1)); + VERIFY_IS_APPROX((m1*m2).determinant(), m1.determinant() * m2.determinant()); + if(size==1) return; + int i = ei_random<int>(0, size-1); + int j; + do { + j = ei_random<int>(0, size-1); + } while(j==i); + m2 = m1; + m2.row(i).swap(m2.row(j)); + VERIFY_IS_APPROX(m2.determinant(), -m1.determinant()); + m2 = m1; + 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(ei_conj(m2.determinant()), m2.adjoint().determinant()); + m2 = m1; + m2.row(i) += x*m2.row(j); + VERIFY_IS_APPROX(m2.determinant(), m1.determinant()); + m2 = m1; + m2.row(i) *= x; + VERIFY_IS_APPROX(m2.determinant(), m1.determinant() * x); +} + +void test_eigen2_determinant() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( determinant(Matrix<float, 1, 1>()) ); + CALL_SUBTEST_2( determinant(Matrix<double, 2, 2>()) ); + CALL_SUBTEST_3( determinant(Matrix<double, 3, 3>()) ); + CALL_SUBTEST_4( determinant(Matrix<double, 4, 4>()) ); + CALL_SUBTEST_5( determinant(Matrix<std::complex<double>, 10, 10>()) ); + CALL_SUBTEST_6( determinant(MatrixXd(20, 20)) ); + } + CALL_SUBTEST_6( determinant(MatrixXd(200, 200)) ); +} diff --git a/eigen/test/eigen2/eigen2_dynalloc.cpp b/eigen/test/eigen2/eigen2_dynalloc.cpp new file mode 100644 index 0000000..1891a9e --- /dev/null +++ b/eigen/test/eigen2/eigen2_dynalloc.cpp @@ -0,0 +1,131 @@ +// 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 Gael Guennebaud <g.gael@free.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" + +#if EIGEN_ARCH_WANTS_ALIGNMENT +#define ALIGNMENT 16 +#else +#define ALIGNMENT 1 +#endif + +void check_handmade_aligned_malloc() +{ + for(int i = 1; i < 1000; i++) + { + char *p = (char*)ei_handmade_aligned_malloc(i); + VERIFY(std::size_t(p)%ALIGNMENT==0); + // if the buffer is wrongly allocated this will give a bad write --> check with valgrind + for(int j = 0; j < i; j++) p[j]=0; + ei_handmade_aligned_free(p); + } +} + +void check_aligned_malloc() +{ + for(int i = 1; i < 1000; i++) + { + char *p = (char*)ei_aligned_malloc(i); + VERIFY(std::size_t(p)%ALIGNMENT==0); + // if the buffer is wrongly allocated this will give a bad write --> check with valgrind + for(int j = 0; j < i; j++) p[j]=0; + ei_aligned_free(p); + } +} + +void check_aligned_new() +{ + for(int i = 1; i < 1000; i++) + { + float *p = ei_aligned_new<float>(i); + VERIFY(std::size_t(p)%ALIGNMENT==0); + // if the buffer is wrongly allocated this will give a bad write --> check with valgrind + for(int j = 0; j < i; j++) p[j]=0; + ei_aligned_delete(p,i); + } +} + +void check_aligned_stack_alloc() +{ + for(int i = 1; i < 1000; i++) + { + ei_declare_aligned_stack_constructed_variable(float, p, i, 0); + VERIFY(std::size_t(p)%ALIGNMENT==0); + // if the buffer is wrongly allocated this will give a bad write --> check with valgrind + for(int j = 0; j < i; j++) p[j]=0; + } +} + + +// test compilation with both a struct and a class... +struct MyStruct +{ + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + char dummychar; + Vector4f avec; +}; + +class MyClassA +{ + public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + char dummychar; + Vector4f avec; +}; + +template<typename T> void check_dynaligned() +{ + T* obj = new T; + VERIFY(std::size_t(obj)%ALIGNMENT==0); + delete obj; +} + +void test_eigen2_dynalloc() +{ + // low level dynamic memory allocation + CALL_SUBTEST(check_handmade_aligned_malloc()); + CALL_SUBTEST(check_aligned_malloc()); + CALL_SUBTEST(check_aligned_new()); + CALL_SUBTEST(check_aligned_stack_alloc()); + + for (int i=0; i<g_repeat*100; ++i) + { + CALL_SUBTEST( check_dynaligned<Vector4f>() ); + CALL_SUBTEST( check_dynaligned<Vector2d>() ); + CALL_SUBTEST( check_dynaligned<Matrix4f>() ); + CALL_SUBTEST( check_dynaligned<Vector4d>() ); + CALL_SUBTEST( check_dynaligned<Vector4i>() ); + } + + // check static allocation, who knows ? + { + MyStruct foo0; VERIFY(std::size_t(foo0.avec.data())%ALIGNMENT==0); + MyClassA fooA; VERIFY(std::size_t(fooA.avec.data())%ALIGNMENT==0); + } + + // dynamic allocation, single object + for (int i=0; i<g_repeat*100; ++i) + { + MyStruct *foo0 = new MyStruct(); VERIFY(std::size_t(foo0->avec.data())%ALIGNMENT==0); + MyClassA *fooA = new MyClassA(); VERIFY(std::size_t(fooA->avec.data())%ALIGNMENT==0); + delete foo0; + delete fooA; + } + + // dynamic allocation, array + const int N = 10; + for (int i=0; i<g_repeat*100; ++i) + { + MyStruct *foo0 = new MyStruct[N]; VERIFY(std::size_t(foo0->avec.data())%ALIGNMENT==0); + MyClassA *fooA = new MyClassA[N]; VERIFY(std::size_t(fooA->avec.data())%ALIGNMENT==0); + delete[] foo0; + delete[] fooA; + } + +} diff --git a/eigen/test/eigen2/eigen2_eigensolver.cpp b/eigen/test/eigen2/eigen2_eigensolver.cpp new file mode 100644 index 0000000..48b4ace --- /dev/null +++ b/eigen/test/eigen2/eigen2_eigensolver.cpp @@ -0,0 +1,146 @@ +// 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 Gael Guennebaud <g.gael@free.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" +#include <Eigen/QR> + +#ifdef HAS_GSL +#include "gsl_helper.h" +#endif + +template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m) +{ + /* this test covers the following files: + EigenSolver.h, SelfAdjointEigenSolver.h (and indirectly: Tridiagonalization.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; + + RealScalar largerEps = 10*test_precision<RealScalar>(); + + MatrixType a = MatrixType::Random(rows,cols); + MatrixType a1 = MatrixType::Random(rows,cols); + MatrixType symmA = a.adjoint() * a + a1.adjoint() * a1; + + MatrixType b = MatrixType::Random(rows,cols); + MatrixType b1 = MatrixType::Random(rows,cols); + MatrixType symmB = b.adjoint() * b + b1.adjoint() * b1; + + SelfAdjointEigenSolver<MatrixType> eiSymm(symmA); + // generalized eigen pb + SelfAdjointEigenSolver<MatrixType> eiSymmGen(symmA, symmB); + + #ifdef HAS_GSL + if (ei_is_same_type<RealScalar,double>::ret) + { + typedef GslTraits<Scalar> Gsl; + typename Gsl::Matrix gEvec=0, gSymmA=0, gSymmB=0; + typename GslTraits<RealScalar>::Vector gEval=0; + RealVectorType _eval; + MatrixType _evec; + convert<MatrixType>(symmA, gSymmA); + convert<MatrixType>(symmB, gSymmB); + convert<MatrixType>(symmA, gEvec); + gEval = GslTraits<RealScalar>::createVector(rows); + + Gsl::eigen_symm(gSymmA, gEval, gEvec); + convert(gEval, _eval); + convert(gEvec, _evec); + + // test gsl itself ! + VERIFY((symmA * _evec).isApprox(_evec * _eval.asDiagonal(), largerEps)); + + // compare with eigen + VERIFY_IS_APPROX(_eval, eiSymm.eigenvalues()); + VERIFY_IS_APPROX(_evec.cwise().abs(), eiSymm.eigenvectors().cwise().abs()); + + // generalized pb + Gsl::eigen_symm_gen(gSymmA, gSymmB, gEval, gEvec); + convert(gEval, _eval); + convert(gEvec, _evec); + // test GSL itself: + VERIFY((symmA * _evec).isApprox(symmB * (_evec * _eval.asDiagonal()), largerEps)); + + // compare with eigen + MatrixType normalized_eivec = eiSymmGen.eigenvectors()*eiSymmGen.eigenvectors().colwise().norm().asDiagonal().inverse(); + VERIFY_IS_APPROX(_eval, eiSymmGen.eigenvalues()); + VERIFY_IS_APPROX(_evec.cwiseAbs(), normalized_eivec.cwiseAbs()); + + Gsl::free(gSymmA); + Gsl::free(gSymmB); + GslTraits<RealScalar>::free(gEval); + Gsl::free(gEvec); + } + #endif + + VERIFY((symmA * eiSymm.eigenvectors()).isApprox( + eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal(), largerEps)); + + // generalized eigen problem Ax = lBx + VERIFY((symmA * eiSymmGen.eigenvectors()).isApprox( + symmB * (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps)); + + MatrixType sqrtSymmA = eiSymm.operatorSqrt(); + VERIFY_IS_APPROX(symmA, sqrtSymmA*sqrtSymmA); + VERIFY_IS_APPROX(sqrtSymmA, symmA*eiSymm.operatorInverseSqrt()); +} + +template<typename MatrixType> void eigensolver(const MatrixType& m) +{ + /* this test covers the following files: + EigenSolver.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; + + // RealScalar largerEps = 10*test_precision<RealScalar>(); + + MatrixType a = MatrixType::Random(rows,cols); + MatrixType a1 = MatrixType::Random(rows,cols); + MatrixType symmA = a.adjoint() * a + a1.adjoint() * a1; + + EigenSolver<MatrixType> ei0(symmA); + VERIFY_IS_APPROX(symmA * ei0.pseudoEigenvectors(), ei0.pseudoEigenvectors() * ei0.pseudoEigenvalueMatrix()); + VERIFY_IS_APPROX((symmA.template cast<Complex>()) * (ei0.pseudoEigenvectors().template cast<Complex>()), + (ei0.pseudoEigenvectors().template cast<Complex>()) * (ei0.eigenvalues().asDiagonal())); + + EigenSolver<MatrixType> ei1(a); + VERIFY_IS_APPROX(a * ei1.pseudoEigenvectors(), ei1.pseudoEigenvectors() * ei1.pseudoEigenvalueMatrix()); + VERIFY_IS_APPROX(a.template cast<Complex>() * ei1.eigenvectors(), + ei1.eigenvectors() * ei1.eigenvalues().asDiagonal()); + +} + +void test_eigen2_eigensolver() +{ + for(int i = 0; i < g_repeat; i++) { + // very important to test a 3x3 matrix since we provide a special path for it + CALL_SUBTEST_1( selfadjointeigensolver(Matrix3f()) ); + CALL_SUBTEST_2( selfadjointeigensolver(Matrix4d()) ); + CALL_SUBTEST_3( selfadjointeigensolver(MatrixXf(7,7)) ); + CALL_SUBTEST_4( selfadjointeigensolver(MatrixXcd(5,5)) ); + CALL_SUBTEST_5( selfadjointeigensolver(MatrixXd(19,19)) ); + + CALL_SUBTEST_6( eigensolver(Matrix4f()) ); + CALL_SUBTEST_5( eigensolver(MatrixXd(17,17)) ); + } +} + diff --git a/eigen/test/eigen2/eigen2_first_aligned.cpp b/eigen/test/eigen2/eigen2_first_aligned.cpp new file mode 100644 index 0000000..51bb3ca --- /dev/null +++ b/eigen/test/eigen2/eigen2_first_aligned.cpp @@ -0,0 +1,49 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" + +template<typename Scalar> +void test_eigen2_first_aligned_helper(Scalar *array, int size) +{ + const int packet_size = sizeof(Scalar) * ei_packet_traits<Scalar>::size; + VERIFY(((std::size_t(array) + sizeof(Scalar) * ei_alignmentOffset(array, size)) % packet_size) == 0); +} + +template<typename Scalar> +void test_eigen2_none_aligned_helper(Scalar *array, int size) +{ + VERIFY(ei_packet_traits<Scalar>::size == 1 || ei_alignmentOffset(array, size) == size); +} + +struct some_non_vectorizable_type { float x; }; + +void test_eigen2_first_aligned() +{ + EIGEN_ALIGN_128 float array_float[100]; + test_first_aligned_helper(array_float, 50); + test_first_aligned_helper(array_float+1, 50); + test_first_aligned_helper(array_float+2, 50); + test_first_aligned_helper(array_float+3, 50); + test_first_aligned_helper(array_float+4, 50); + test_first_aligned_helper(array_float+5, 50); + + EIGEN_ALIGN_128 double array_double[100]; + test_first_aligned_helper(array_double, 50); + test_first_aligned_helper(array_double+1, 50); + test_first_aligned_helper(array_double+2, 50); + + double *array_double_plus_4_bytes = (double*)(std::size_t(array_double)+4); + test_none_aligned_helper(array_double_plus_4_bytes, 50); + test_none_aligned_helper(array_double_plus_4_bytes+1, 50); + + some_non_vectorizable_type array_nonvec[100]; + test_first_aligned_helper(array_nonvec, 100); + test_none_aligned_helper(array_nonvec, 100); +} diff --git a/eigen/test/eigen2/eigen2_geometry.cpp b/eigen/test/eigen2/eigen2_geometry.cpp new file mode 100644 index 0000000..5140407 --- /dev/null +++ b/eigen/test/eigen2/eigen2_geometry.cpp @@ -0,0 +1,432 @@ +// 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 Gael Guennebaud <g.gael@free.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" +#include <Eigen/Geometry> +#include <Eigen/LU> +#include <Eigen/SVD> + +template<typename Scalar> void geometry(void) +{ + /* this test covers the following files: + Cross.h Quaternion.h, Transform.cpp + */ + + typedef Matrix<Scalar,2,2> Matrix2; + typedef Matrix<Scalar,3,3> Matrix3; + typedef Matrix<Scalar,4,4> Matrix4; + typedef Matrix<Scalar,2,1> Vector2; + typedef Matrix<Scalar,3,1> Vector3; + typedef Matrix<Scalar,4,1> Vector4; + typedef Quaternion<Scalar> Quaternionx; + typedef AngleAxis<Scalar> AngleAxisx; + typedef Transform<Scalar,2> Transform2; + typedef Transform<Scalar,3> Transform3; + typedef Scaling<Scalar,2> Scaling2; + typedef Scaling<Scalar,3> Scaling3; + typedef Translation<Scalar,2> Translation2; + typedef Translation<Scalar,3> Translation3; + + Scalar largeEps = test_precision<Scalar>(); + if (ei_is_same_type<Scalar,float>::ret) + largeEps = 1e-2f; + + Vector3 v0 = Vector3::Random(), + v1 = Vector3::Random(), + v2 = Vector3::Random(); + Vector2 u0 = Vector2::Random(); + Matrix3 matrot1; + + Scalar a = ei_random<Scalar>(-Scalar(M_PI), Scalar(M_PI)); + + // cross product + VERIFY_IS_MUCH_SMALLER_THAN(v1.cross(v2).eigen2_dot(v1), Scalar(1)); + Matrix3 m; + m << v0.normalized(), + (v0.cross(v1)).normalized(), + (v0.cross(v1).cross(v0)).normalized(); + VERIFY(m.isUnitary()); + + // Quaternion: Identity(), setIdentity(); + Quaternionx q1, q2; + q2.setIdentity(); + VERIFY_IS_APPROX(Quaternionx(Quaternionx::Identity()).coeffs(), q2.coeffs()); + q1.coeffs().setRandom(); + VERIFY_IS_APPROX(q1.coeffs(), (q1*q2).coeffs()); + + // unitOrthogonal + VERIFY_IS_MUCH_SMALLER_THAN(u0.unitOrthogonal().eigen2_dot(u0), Scalar(1)); + VERIFY_IS_MUCH_SMALLER_THAN(v0.unitOrthogonal().eigen2_dot(v0), Scalar(1)); + VERIFY_IS_APPROX(u0.unitOrthogonal().norm(), Scalar(1)); + VERIFY_IS_APPROX(v0.unitOrthogonal().norm(), Scalar(1)); + + + VERIFY_IS_APPROX(v0, AngleAxisx(a, v0.normalized()) * v0); + VERIFY_IS_APPROX(-v0, AngleAxisx(Scalar(M_PI), v0.unitOrthogonal()) * v0); + VERIFY_IS_APPROX(ei_cos(a)*v0.squaredNorm(), v0.eigen2_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); + + q1 = AngleAxisx(a, v0.normalized()); + q2 = AngleAxisx(a, v1.normalized()); + + // angular distance + Scalar refangle = ei_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(ei_isApprox(q1.angularDistance(q2), refangle, largeEps)); + } + + // rotation matrix conversion + VERIFY_IS_APPROX(q1 * v2, q1.toRotationMatrix() * v2); + VERIFY_IS_APPROX(q1 * q2 * v2, + q1.toRotationMatrix() * q2.toRotationMatrix() * v2); + + VERIFY( (q2*q1).isApprox(q1*q2, largeEps) || !(q2 * q1 * v2).isApprox( + q1.toRotationMatrix() * q2.toRotationMatrix() * v2)); + + q2 = q1.toRotationMatrix(); + VERIFY_IS_APPROX(q1*v1,q2*v1); + + matrot1 = AngleAxisx(Scalar(0.1), Vector3::UnitX()) + * AngleAxisx(Scalar(0.2), Vector3::UnitY()) + * AngleAxisx(Scalar(0.3), Vector3::UnitZ()); + VERIFY_IS_APPROX(matrot1 * v1, + AngleAxisx(Scalar(0.1), Vector3(1,0,0)).toRotationMatrix() + * (AngleAxisx(Scalar(0.2), Vector3(0,1,0)).toRotationMatrix() + * (AngleAxisx(Scalar(0.3), Vector3(0,0,1)).toRotationMatrix() * v1))); + + // angle-axis conversion + AngleAxisx aa = q1; + VERIFY_IS_APPROX(q1 * v1, Quaternionx(aa) * v1); + VERIFY_IS_NOT_APPROX(q1 * v1, Quaternionx(AngleAxisx(aa.angle()*2,aa.axis())) * v1); + + // from two vector creation + VERIFY_IS_APPROX(v2.normalized(),(q2.setFromTwoVectors(v1,v2)*v1).normalized()); + VERIFY_IS_APPROX(v2.normalized(),(q2.setFromTwoVectors(v1,v2)*v1).normalized()); + + // inverse and conjugate + VERIFY_IS_APPROX(q1 * (q1.inverse() * v1), v1); + VERIFY_IS_APPROX(q1 * (q1.conjugate() * v1), v1); + + // AngleAxis + VERIFY_IS_APPROX(AngleAxisx(a,v1.normalized()).toRotationMatrix(), + Quaternionx(AngleAxisx(a,v1.normalized())).toRotationMatrix()); + + AngleAxisx aa1; + m = q1.toRotationMatrix(); + aa1 = m; + VERIFY_IS_APPROX(AngleAxisx(m).toRotationMatrix(), + Quaternionx(m).toRotationMatrix()); + + // Transform + // TODO complete the tests ! + a = 0; + while (ei_abs(a)<Scalar(0.1)) + a = ei_random<Scalar>(-Scalar(0.4)*Scalar(M_PI), Scalar(0.4)*Scalar(M_PI)); + q1 = AngleAxisx(a, v0.normalized()); + Transform3 t0, t1, t2; + // first test setIdentity() and Identity() + t0.setIdentity(); + VERIFY_IS_APPROX(t0.matrix(), Transform3::MatrixType::Identity()); + t0.matrix().setZero(); + t0 = Transform3::Identity(); + VERIFY_IS_APPROX(t0.matrix(), Transform3::MatrixType::Identity()); + + t0.linear() = q1.toRotationMatrix(); + t1.setIdentity(); + t1.linear() = q1.toRotationMatrix(); + + v0 << 50, 2, 1;//= ei_random_matrix<Vector3>().cwiseProduct(Vector3(10,2,0.5)); + t0.scale(v0); + t1.prescale(v0); + + VERIFY_IS_APPROX( (t0 * Vector3(1,0,0)).norm(), v0.x()); + //VERIFY(!ei_isApprox((t1 * Vector3(1,0,0)).norm(), v0.x())); + + t0.setIdentity(); + t1.setIdentity(); + v1 << 1, 2, 3; + t0.linear() = q1.toRotationMatrix(); + t0.pretranslate(v0); + t0.scale(v1); + t1.linear() = q1.conjugate().toRotationMatrix(); + t1.prescale(v1.cwise().inverse()); + t1.translate(-v0); + + VERIFY((t0.matrix() * t1.matrix()).isIdentity(test_precision<Scalar>())); + + t1.fromPositionOrientationScale(v0, q1, v1); + VERIFY_IS_APPROX(t1.matrix(), t0.matrix()); + VERIFY_IS_APPROX(t1*v1, t0*v1); + + t0.setIdentity(); t0.scale(v0).rotate(q1.toRotationMatrix()); + t1.setIdentity(); t1.scale(v0).rotate(q1); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + + t0.setIdentity(); t0.scale(v0).rotate(AngleAxisx(q1)); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + + VERIFY_IS_APPROX(t0.scale(a).matrix(), t1.scale(Vector3::Constant(a)).matrix()); + VERIFY_IS_APPROX(t0.prescale(a).matrix(), t1.prescale(Vector3::Constant(a)).matrix()); + + // More transform constructors, operator=, operator*= + + Matrix3 mat3 = Matrix3::Random(); + Matrix4 mat4; + mat4 << mat3 , Vector3::Zero() , Vector4::Zero().transpose(); + Transform3 tmat3(mat3), tmat4(mat4); + tmat4.matrix()(3,3) = Scalar(1); + VERIFY_IS_APPROX(tmat3.matrix(), tmat4.matrix()); + + Scalar a3 = ei_random<Scalar>(-Scalar(M_PI), Scalar(M_PI)); + Vector3 v3 = Vector3::Random().normalized(); + AngleAxisx aa3(a3, v3); + Transform3 t3(aa3); + Transform3 t4; + t4 = aa3; + VERIFY_IS_APPROX(t3.matrix(), t4.matrix()); + t4.rotate(AngleAxisx(-a3,v3)); + VERIFY_IS_APPROX(t4.matrix(), Matrix4::Identity()); + t4 *= aa3; + VERIFY_IS_APPROX(t3.matrix(), t4.matrix()); + + v3 = Vector3::Random(); + Translation3 tv3(v3); + Transform3 t5(tv3); + t4 = tv3; + VERIFY_IS_APPROX(t5.matrix(), t4.matrix()); + t4.translate(-v3); + VERIFY_IS_APPROX(t4.matrix(), Matrix4::Identity()); + t4 *= tv3; + VERIFY_IS_APPROX(t5.matrix(), t4.matrix()); + + Scaling3 sv3(v3); + Transform3 t6(sv3); + t4 = sv3; + VERIFY_IS_APPROX(t6.matrix(), t4.matrix()); + t4.scale(v3.cwise().inverse()); + VERIFY_IS_APPROX(t4.matrix(), Matrix4::Identity()); + t4 *= sv3; + VERIFY_IS_APPROX(t6.matrix(), t4.matrix()); + + // matrix * transform + VERIFY_IS_APPROX(Transform3(t3.matrix()*t4).matrix(), Transform3(t3*t4).matrix()); + + // chained Transform product + VERIFY_IS_APPROX(((t3*t4)*t5).matrix(), (t3*(t4*t5)).matrix()); + + // check that Transform product doesn't have aliasing problems + t5 = t4; + t5 = t5*t5; + VERIFY_IS_APPROX(t5, t4*t4); + + // 2D transformation + Transform2 t20, t21; + Vector2 v20 = Vector2::Random(); + Vector2 v21 = Vector2::Random(); + for (int k=0; k<2; ++k) + if (ei_abs(v21[k])<Scalar(1e-3)) v21[k] = Scalar(1e-3); + t21.setIdentity(); + t21.linear() = Rotation2D<Scalar>(a).toRotationMatrix(); + VERIFY_IS_APPROX(t20.fromPositionOrientationScale(v20,a,v21).matrix(), + t21.pretranslate(v20).scale(v21).matrix()); + + t21.setIdentity(); + t21.linear() = Rotation2D<Scalar>(-a).toRotationMatrix(); + VERIFY( (t20.fromPositionOrientationScale(v20,a,v21) + * (t21.prescale(v21.cwise().inverse()).translate(-v20))).matrix().isIdentity(test_precision<Scalar>()) ); + + // Transform - new API + // 3D + t0.setIdentity(); + t0.rotate(q1).scale(v0).translate(v0); + // mat * scaling and mat * translation + t1 = (Matrix3(q1) * Scaling3(v0)) * Translation3(v0); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + // mat * transformation and scaling * translation + t1 = Matrix3(q1) * (Scaling3(v0) * Translation3(v0)); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + + t0.setIdentity(); + t0.prerotate(q1).prescale(v0).pretranslate(v0); + // translation * scaling and transformation * mat + t1 = (Translation3(v0) * Scaling3(v0)) * Matrix3(q1); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + // scaling * mat and translation * mat + t1 = Translation3(v0) * (Scaling3(v0) * Matrix3(q1)); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + + t0.setIdentity(); + t0.scale(v0).translate(v0).rotate(q1); + // translation * mat and scaling * transformation + t1 = Scaling3(v0) * (Translation3(v0) * Matrix3(q1)); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + // transformation * scaling + t0.scale(v0); + t1 = t1 * Scaling3(v0); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + // transformation * translation + t0.translate(v0); + t1 = t1 * Translation3(v0); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + // translation * transformation + t0.pretranslate(v0); + t1 = Translation3(v0) * t1; + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + + // transform * quaternion + t0.rotate(q1); + t1 = t1 * q1; + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + + // translation * quaternion + t0.translate(v1).rotate(q1); + t1 = t1 * (Translation3(v1) * q1); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + + // scaling * quaternion + t0.scale(v1).rotate(q1); + t1 = t1 * (Scaling3(v1) * q1); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + + // quaternion * transform + t0.prerotate(q1); + t1 = q1 * t1; + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + + // quaternion * translation + t0.rotate(q1).translate(v1); + t1 = t1 * (q1 * Translation3(v1)); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + + // quaternion * scaling + t0.rotate(q1).scale(v1); + t1 = t1 * (q1 * Scaling3(v1)); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + + // translation * vector + t0.setIdentity(); + t0.translate(v0); + VERIFY_IS_APPROX(t0 * v1, Translation3(v0) * v1); + + // scaling * vector + t0.setIdentity(); + t0.scale(v0); + VERIFY_IS_APPROX(t0 * v1, Scaling3(v0) * v1); + + // test transform inversion + t0.setIdentity(); + t0.translate(v0); + t0.linear().setRandom(); + VERIFY_IS_APPROX(t0.inverse(Affine), t0.matrix().inverse()); + t0.setIdentity(); + t0.translate(v0).rotate(q1); + VERIFY_IS_APPROX(t0.inverse(Isometry), t0.matrix().inverse()); + + // test extract rotation and scaling + t0.setIdentity(); + t0.translate(v0).rotate(q1).scale(v1); + VERIFY_IS_APPROX(t0.rotation() * v1, Matrix3(q1) * v1); + + Matrix3 mat_rotation, mat_scaling; + t0.setIdentity(); + t0.translate(v0).rotate(q1).scale(v1); + t0.computeRotationScaling(&mat_rotation, &mat_scaling); + VERIFY_IS_APPROX(t0.linear(), mat_rotation * mat_scaling); + VERIFY_IS_APPROX(mat_rotation*mat_rotation.adjoint(), Matrix3::Identity()); + VERIFY_IS_APPROX(mat_rotation.determinant(), Scalar(1)); + t0.computeScalingRotation(&mat_scaling, &mat_rotation); + VERIFY_IS_APPROX(t0.linear(), mat_scaling * mat_rotation); + VERIFY_IS_APPROX(mat_rotation*mat_rotation.adjoint(), Matrix3::Identity()); + VERIFY_IS_APPROX(mat_rotation.determinant(), Scalar(1)); + + // test casting + Transform<float,3> t1f = t1.template cast<float>(); + VERIFY_IS_APPROX(t1f.template cast<Scalar>(),t1); + Transform<double,3> t1d = t1.template cast<double>(); + VERIFY_IS_APPROX(t1d.template cast<Scalar>(),t1); + + Translation3 tr1(v0); + Translation<float,3> tr1f = tr1.template cast<float>(); + VERIFY_IS_APPROX(tr1f.template cast<Scalar>(),tr1); + Translation<double,3> tr1d = tr1.template cast<double>(); + VERIFY_IS_APPROX(tr1d.template cast<Scalar>(),tr1); + + Scaling3 sc1(v0); + Scaling<float,3> sc1f = sc1.template cast<float>(); + VERIFY_IS_APPROX(sc1f.template cast<Scalar>(),sc1); + Scaling<double,3> sc1d = sc1.template cast<double>(); + VERIFY_IS_APPROX(sc1d.template cast<Scalar>(),sc1); + + Quaternion<float> q1f = q1.template cast<float>(); + VERIFY_IS_APPROX(q1f.template cast<Scalar>(),q1); + Quaternion<double> q1d = q1.template cast<double>(); + VERIFY_IS_APPROX(q1d.template cast<Scalar>(),q1); + + AngleAxis<float> aa1f = aa1.template cast<float>(); + VERIFY_IS_APPROX(aa1f.template cast<Scalar>(),aa1); + AngleAxis<double> aa1d = aa1.template cast<double>(); + VERIFY_IS_APPROX(aa1d.template cast<Scalar>(),aa1); + + Rotation2D<Scalar> r2d1(ei_random<Scalar>()); + Rotation2D<float> r2d1f = r2d1.template cast<float>(); + VERIFY_IS_APPROX(r2d1f.template cast<Scalar>(),r2d1); + Rotation2D<double> r2d1d = r2d1.template cast<double>(); + VERIFY_IS_APPROX(r2d1d.template cast<Scalar>(),r2d1); + + m = q1; +// m.col(1) = Vector3(0,ei_random<Scalar>(),ei_random<Scalar>()).normalized(); +// m.col(0) = Vector3(-1,0,0).normalized(); +// m.col(2) = m.col(0).cross(m.col(1)); + #define VERIFY_EULER(I,J,K, X,Y,Z) { \ + Vector3 ea = m.eulerAngles(I,J,K); \ + Matrix3 m1 = Matrix3(AngleAxisx(ea[0], Vector3::Unit##X()) * AngleAxisx(ea[1], Vector3::Unit##Y()) * AngleAxisx(ea[2], Vector3::Unit##Z())); \ + VERIFY_IS_APPROX(m, m1); \ + 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); + + // colwise/rowwise cross product + mat3.setRandom(); + Vector3 vec3 = Vector3::Random(); + Matrix3 mcross; + int i = ei_random<int>(0,2); + mcross = mat3.colwise().cross(vec3); + VERIFY_IS_APPROX(mcross.col(i), mat3.col(i).cross(vec3)); + mcross = mat3.rowwise().cross(vec3); + VERIFY_IS_APPROX(mcross.row(i), mat3.row(i).cross(vec3)); + + +} + +void test_eigen2_geometry() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( geometry<float>() ); + CALL_SUBTEST_2( geometry<double>() ); + } +} diff --git a/eigen/test/eigen2/eigen2_geometry_with_eigen2_prefix.cpp b/eigen/test/eigen2/eigen2_geometry_with_eigen2_prefix.cpp new file mode 100644 index 0000000..12d4a71 --- /dev/null +++ b/eigen/test/eigen2/eigen2_geometry_with_eigen2_prefix.cpp @@ -0,0 +1,435 @@ +// 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 Gael Guennebaud <g.gael@free.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#define EIGEN2_SUPPORT_STAGE15_RESOLVE_API_CONFLICTS_WARN + +#include "main.h" +#include <Eigen/Geometry> +#include <Eigen/LU> +#include <Eigen/SVD> + +template<typename Scalar> void geometry(void) +{ + /* this test covers the following files: + Cross.h Quaternion.h, Transform.cpp + */ + + typedef Matrix<Scalar,2,2> Matrix2; + typedef Matrix<Scalar,3,3> Matrix3; + typedef Matrix<Scalar,4,4> Matrix4; + typedef Matrix<Scalar,2,1> Vector2; + typedef Matrix<Scalar,3,1> Vector3; + typedef Matrix<Scalar,4,1> Vector4; + typedef eigen2_Quaternion<Scalar> Quaternionx; + typedef eigen2_AngleAxis<Scalar> AngleAxisx; + typedef eigen2_Transform<Scalar,2> Transform2; + typedef eigen2_Transform<Scalar,3> Transform3; + typedef eigen2_Scaling<Scalar,2> Scaling2; + typedef eigen2_Scaling<Scalar,3> Scaling3; + typedef eigen2_Translation<Scalar,2> Translation2; + typedef eigen2_Translation<Scalar,3> Translation3; + + Scalar largeEps = test_precision<Scalar>(); + if (ei_is_same_type<Scalar,float>::ret) + largeEps = 1e-2f; + + Vector3 v0 = Vector3::Random(), + v1 = Vector3::Random(), + v2 = Vector3::Random(); + Vector2 u0 = Vector2::Random(); + Matrix3 matrot1; + + Scalar a = ei_random<Scalar>(-Scalar(M_PI), Scalar(M_PI)); + + // cross product + VERIFY_IS_MUCH_SMALLER_THAN(v1.cross(v2).eigen2_dot(v1), Scalar(1)); + Matrix3 m; + m << v0.normalized(), + (v0.cross(v1)).normalized(), + (v0.cross(v1).cross(v0)).normalized(); + VERIFY(m.isUnitary()); + + // Quaternion: Identity(), setIdentity(); + Quaternionx q1, q2; + q2.setIdentity(); + VERIFY_IS_APPROX(Quaternionx(Quaternionx::Identity()).coeffs(), q2.coeffs()); + q1.coeffs().setRandom(); + VERIFY_IS_APPROX(q1.coeffs(), (q1*q2).coeffs()); + + // unitOrthogonal + VERIFY_IS_MUCH_SMALLER_THAN(u0.unitOrthogonal().eigen2_dot(u0), Scalar(1)); + VERIFY_IS_MUCH_SMALLER_THAN(v0.unitOrthogonal().eigen2_dot(v0), Scalar(1)); + VERIFY_IS_APPROX(u0.unitOrthogonal().norm(), Scalar(1)); + VERIFY_IS_APPROX(v0.unitOrthogonal().norm(), Scalar(1)); + + + VERIFY_IS_APPROX(v0, AngleAxisx(a, v0.normalized()) * v0); + VERIFY_IS_APPROX(-v0, AngleAxisx(Scalar(M_PI), v0.unitOrthogonal()) * v0); + VERIFY_IS_APPROX(ei_cos(a)*v0.squaredNorm(), v0.eigen2_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); + + q1 = AngleAxisx(a, v0.normalized()); + q2 = AngleAxisx(a, v1.normalized()); + + // angular distance + Scalar refangle = ei_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(ei_isApprox(q1.angularDistance(q2), refangle, largeEps)); + } + + // rotation matrix conversion + VERIFY_IS_APPROX(q1 * v2, q1.toRotationMatrix() * v2); + VERIFY_IS_APPROX(q1 * q2 * v2, + q1.toRotationMatrix() * q2.toRotationMatrix() * v2); + + VERIFY( (q2*q1).isApprox(q1*q2, largeEps) || !(q2 * q1 * v2).isApprox( + q1.toRotationMatrix() * q2.toRotationMatrix() * v2)); + + q2 = q1.toRotationMatrix(); + VERIFY_IS_APPROX(q1*v1,q2*v1); + + matrot1 = AngleAxisx(Scalar(0.1), Vector3::UnitX()) + * AngleAxisx(Scalar(0.2), Vector3::UnitY()) + * AngleAxisx(Scalar(0.3), Vector3::UnitZ()); + VERIFY_IS_APPROX(matrot1 * v1, + AngleAxisx(Scalar(0.1), Vector3(1,0,0)).toRotationMatrix() + * (AngleAxisx(Scalar(0.2), Vector3(0,1,0)).toRotationMatrix() + * (AngleAxisx(Scalar(0.3), Vector3(0,0,1)).toRotationMatrix() * v1))); + + // angle-axis conversion + AngleAxisx aa = q1; + VERIFY_IS_APPROX(q1 * v1, Quaternionx(aa) * v1); + VERIFY_IS_NOT_APPROX(q1 * v1, Quaternionx(AngleAxisx(aa.angle()*2,aa.axis())) * v1); + + // from two vector creation + VERIFY_IS_APPROX(v2.normalized(),(q2.setFromTwoVectors(v1,v2)*v1).normalized()); + VERIFY_IS_APPROX(v2.normalized(),(q2.setFromTwoVectors(v1,v2)*v1).normalized()); + + // inverse and conjugate + VERIFY_IS_APPROX(q1 * (q1.inverse() * v1), v1); + VERIFY_IS_APPROX(q1 * (q1.conjugate() * v1), v1); + + // AngleAxis + VERIFY_IS_APPROX(AngleAxisx(a,v1.normalized()).toRotationMatrix(), + Quaternionx(AngleAxisx(a,v1.normalized())).toRotationMatrix()); + + AngleAxisx aa1; + m = q1.toRotationMatrix(); + aa1 = m; + VERIFY_IS_APPROX(AngleAxisx(m).toRotationMatrix(), + Quaternionx(m).toRotationMatrix()); + + // Transform + // TODO complete the tests ! + a = 0; + while (ei_abs(a)<Scalar(0.1)) + a = ei_random<Scalar>(-Scalar(0.4)*Scalar(M_PI), Scalar(0.4)*Scalar(M_PI)); + q1 = AngleAxisx(a, v0.normalized()); + Transform3 t0, t1, t2; + // first test setIdentity() and Identity() + t0.setIdentity(); + VERIFY_IS_APPROX(t0.matrix(), Transform3::MatrixType::Identity()); + t0.matrix().setZero(); + t0 = Transform3::Identity(); + VERIFY_IS_APPROX(t0.matrix(), Transform3::MatrixType::Identity()); + + t0.linear() = q1.toRotationMatrix(); + t1.setIdentity(); + t1.linear() = q1.toRotationMatrix(); + + v0 << 50, 2, 1;//= ei_random_matrix<Vector3>().cwiseProduct(Vector3(10,2,0.5)); + t0.scale(v0); + t1.prescale(v0); + + VERIFY_IS_APPROX( (t0 * Vector3(1,0,0)).norm(), v0.x()); + //VERIFY(!ei_isApprox((t1 * Vector3(1,0,0)).norm(), v0.x())); + + t0.setIdentity(); + t1.setIdentity(); + v1 << 1, 2, 3; + t0.linear() = q1.toRotationMatrix(); + t0.pretranslate(v0); + t0.scale(v1); + t1.linear() = q1.conjugate().toRotationMatrix(); + t1.prescale(v1.cwise().inverse()); + t1.translate(-v0); + + VERIFY((t0.matrix() * t1.matrix()).isIdentity(test_precision<Scalar>())); + + t1.fromPositionOrientationScale(v0, q1, v1); + VERIFY_IS_APPROX(t1.matrix(), t0.matrix()); + VERIFY_IS_APPROX(t1*v1, t0*v1); + + t0.setIdentity(); t0.scale(v0).rotate(q1.toRotationMatrix()); + t1.setIdentity(); t1.scale(v0).rotate(q1); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + + t0.setIdentity(); t0.scale(v0).rotate(AngleAxisx(q1)); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + + VERIFY_IS_APPROX(t0.scale(a).matrix(), t1.scale(Vector3::Constant(a)).matrix()); + VERIFY_IS_APPROX(t0.prescale(a).matrix(), t1.prescale(Vector3::Constant(a)).matrix()); + + // More transform constructors, operator=, operator*= + + Matrix3 mat3 = Matrix3::Random(); + Matrix4 mat4; + mat4 << mat3 , Vector3::Zero() , Vector4::Zero().transpose(); + Transform3 tmat3(mat3), tmat4(mat4); + tmat4.matrix()(3,3) = Scalar(1); + VERIFY_IS_APPROX(tmat3.matrix(), tmat4.matrix()); + + Scalar a3 = ei_random<Scalar>(-Scalar(M_PI), Scalar(M_PI)); + Vector3 v3 = Vector3::Random().normalized(); + AngleAxisx aa3(a3, v3); + Transform3 t3(aa3); + Transform3 t4; + t4 = aa3; + VERIFY_IS_APPROX(t3.matrix(), t4.matrix()); + t4.rotate(AngleAxisx(-a3,v3)); + VERIFY_IS_APPROX(t4.matrix(), Matrix4::Identity()); + t4 *= aa3; + VERIFY_IS_APPROX(t3.matrix(), t4.matrix()); + + v3 = Vector3::Random(); + Translation3 tv3(v3); + Transform3 t5(tv3); + t4 = tv3; + VERIFY_IS_APPROX(t5.matrix(), t4.matrix()); + t4.translate(-v3); + VERIFY_IS_APPROX(t4.matrix(), Matrix4::Identity()); + t4 *= tv3; + VERIFY_IS_APPROX(t5.matrix(), t4.matrix()); + + Scaling3 sv3(v3); + Transform3 t6(sv3); + t4 = sv3; + VERIFY_IS_APPROX(t6.matrix(), t4.matrix()); + t4.scale(v3.cwise().inverse()); + VERIFY_IS_APPROX(t4.matrix(), Matrix4::Identity()); + t4 *= sv3; + VERIFY_IS_APPROX(t6.matrix(), t4.matrix()); + + // matrix * transform + VERIFY_IS_APPROX(Transform3(t3.matrix()*t4).matrix(), Transform3(t3*t4).matrix()); + + // chained Transform product + VERIFY_IS_APPROX(((t3*t4)*t5).matrix(), (t3*(t4*t5)).matrix()); + + // check that Transform product doesn't have aliasing problems + t5 = t4; + t5 = t5*t5; + VERIFY_IS_APPROX(t5, t4*t4); + + // 2D transformation + Transform2 t20, t21; + Vector2 v20 = Vector2::Random(); + Vector2 v21 = Vector2::Random(); + for (int k=0; k<2; ++k) + if (ei_abs(v21[k])<Scalar(1e-3)) v21[k] = Scalar(1e-3); + t21.setIdentity(); + t21.linear() = Rotation2D<Scalar>(a).toRotationMatrix(); + VERIFY_IS_APPROX(t20.fromPositionOrientationScale(v20,a,v21).matrix(), + t21.pretranslate(v20).scale(v21).matrix()); + + t21.setIdentity(); + t21.linear() = Rotation2D<Scalar>(-a).toRotationMatrix(); + VERIFY( (t20.fromPositionOrientationScale(v20,a,v21) + * (t21.prescale(v21.cwise().inverse()).translate(-v20))).matrix().isIdentity(test_precision<Scalar>()) ); + + // Transform - new API + // 3D + t0.setIdentity(); + t0.rotate(q1).scale(v0).translate(v0); + // mat * scaling and mat * translation + t1 = (Matrix3(q1) * Scaling3(v0)) * Translation3(v0); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + // mat * transformation and scaling * translation + t1 = Matrix3(q1) * (Scaling3(v0) * Translation3(v0)); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + + t0.setIdentity(); + t0.prerotate(q1).prescale(v0).pretranslate(v0); + // translation * scaling and transformation * mat + t1 = (Translation3(v0) * Scaling3(v0)) * Matrix3(q1); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + // scaling * mat and translation * mat + t1 = Translation3(v0) * (Scaling3(v0) * Matrix3(q1)); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + + t0.setIdentity(); + t0.scale(v0).translate(v0).rotate(q1); + // translation * mat and scaling * transformation + t1 = Scaling3(v0) * (Translation3(v0) * Matrix3(q1)); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + // transformation * scaling + t0.scale(v0); + t1 = t1 * Scaling3(v0); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + // transformation * translation + t0.translate(v0); + t1 = t1 * Translation3(v0); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + // translation * transformation + t0.pretranslate(v0); + t1 = Translation3(v0) * t1; + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + + // transform * quaternion + t0.rotate(q1); + t1 = t1 * q1; + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + + // translation * quaternion + t0.translate(v1).rotate(q1); + t1 = t1 * (Translation3(v1) * q1); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + + // scaling * quaternion + t0.scale(v1).rotate(q1); + t1 = t1 * (Scaling3(v1) * q1); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + + // quaternion * transform + t0.prerotate(q1); + t1 = q1 * t1; + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + + // quaternion * translation + t0.rotate(q1).translate(v1); + t1 = t1 * (q1 * Translation3(v1)); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + + // quaternion * scaling + t0.rotate(q1).scale(v1); + t1 = t1 * (q1 * Scaling3(v1)); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + + // translation * vector + t0.setIdentity(); + t0.translate(v0); + VERIFY_IS_APPROX(t0 * v1, Translation3(v0) * v1); + + // scaling * vector + t0.setIdentity(); + t0.scale(v0); + VERIFY_IS_APPROX(t0 * v1, Scaling3(v0) * v1); + + // test transform inversion + t0.setIdentity(); + t0.translate(v0); + t0.linear().setRandom(); + VERIFY_IS_APPROX(t0.inverse(Affine), t0.matrix().inverse()); + t0.setIdentity(); + t0.translate(v0).rotate(q1); + VERIFY_IS_APPROX(t0.inverse(Isometry), t0.matrix().inverse()); + + // test extract rotation and scaling + t0.setIdentity(); + t0.translate(v0).rotate(q1).scale(v1); + VERIFY_IS_APPROX(t0.rotation() * v1, Matrix3(q1) * v1); + + Matrix3 mat_rotation, mat_scaling; + t0.setIdentity(); + t0.translate(v0).rotate(q1).scale(v1); + t0.computeRotationScaling(&mat_rotation, &mat_scaling); + VERIFY_IS_APPROX(t0.linear(), mat_rotation * mat_scaling); + VERIFY_IS_APPROX(mat_rotation*mat_rotation.adjoint(), Matrix3::Identity()); + VERIFY_IS_APPROX(mat_rotation.determinant(), Scalar(1)); + t0.computeScalingRotation(&mat_scaling, &mat_rotation); + VERIFY_IS_APPROX(t0.linear(), mat_scaling * mat_rotation); + VERIFY_IS_APPROX(mat_rotation*mat_rotation.adjoint(), Matrix3::Identity()); + VERIFY_IS_APPROX(mat_rotation.determinant(), Scalar(1)); + + // test casting + eigen2_Transform<float,3> t1f = t1.template cast<float>(); + VERIFY_IS_APPROX(t1f.template cast<Scalar>(),t1); + eigen2_Transform<double,3> t1d = t1.template cast<double>(); + VERIFY_IS_APPROX(t1d.template cast<Scalar>(),t1); + + Translation3 tr1(v0); + eigen2_Translation<float,3> tr1f = tr1.template cast<float>(); + VERIFY_IS_APPROX(tr1f.template cast<Scalar>(),tr1); + eigen2_Translation<double,3> tr1d = tr1.template cast<double>(); + VERIFY_IS_APPROX(tr1d.template cast<Scalar>(),tr1); + + Scaling3 sc1(v0); + eigen2_Scaling<float,3> sc1f = sc1.template cast<float>(); + VERIFY_IS_APPROX(sc1f.template cast<Scalar>(),sc1); + eigen2_Scaling<double,3> sc1d = sc1.template cast<double>(); + VERIFY_IS_APPROX(sc1d.template cast<Scalar>(),sc1); + + eigen2_Quaternion<float> q1f = q1.template cast<float>(); + VERIFY_IS_APPROX(q1f.template cast<Scalar>(),q1); + eigen2_Quaternion<double> q1d = q1.template cast<double>(); + VERIFY_IS_APPROX(q1d.template cast<Scalar>(),q1); + + eigen2_AngleAxis<float> aa1f = aa1.template cast<float>(); + VERIFY_IS_APPROX(aa1f.template cast<Scalar>(),aa1); + eigen2_AngleAxis<double> aa1d = aa1.template cast<double>(); + VERIFY_IS_APPROX(aa1d.template cast<Scalar>(),aa1); + + eigen2_Rotation2D<Scalar> r2d1(ei_random<Scalar>()); + eigen2_Rotation2D<float> r2d1f = r2d1.template cast<float>(); + VERIFY_IS_APPROX(r2d1f.template cast<Scalar>(),r2d1); + eigen2_Rotation2D<double> r2d1d = r2d1.template cast<double>(); + VERIFY_IS_APPROX(r2d1d.template cast<Scalar>(),r2d1); + + m = q1; +// m.col(1) = Vector3(0,ei_random<Scalar>(),ei_random<Scalar>()).normalized(); +// m.col(0) = Vector3(-1,0,0).normalized(); +// m.col(2) = m.col(0).cross(m.col(1)); + #define VERIFY_EULER(I,J,K, X,Y,Z) { \ + Vector3 ea = m.eulerAngles(I,J,K); \ + Matrix3 m1 = Matrix3(AngleAxisx(ea[0], Vector3::Unit##X()) * AngleAxisx(ea[1], Vector3::Unit##Y()) * AngleAxisx(ea[2], Vector3::Unit##Z())); \ + VERIFY_IS_APPROX(m, m1); \ + 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); + + // colwise/rowwise cross product + mat3.setRandom(); + Vector3 vec3 = Vector3::Random(); + Matrix3 mcross; + int i = ei_random<int>(0,2); + mcross = mat3.colwise().cross(vec3); + VERIFY_IS_APPROX(mcross.col(i), mat3.col(i).cross(vec3)); + mcross = mat3.rowwise().cross(vec3); + VERIFY_IS_APPROX(mcross.row(i), mat3.row(i).cross(vec3)); + + +} + +void test_eigen2_geometry_with_eigen2_prefix() +{ + std::cout << "eigen2 support: " << EIGEN2_SUPPORT_STAGE << std::endl; + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( geometry<float>() ); + CALL_SUBTEST_2( geometry<double>() ); + } +} diff --git a/eigen/test/eigen2/eigen2_hyperplane.cpp b/eigen/test/eigen2/eigen2_hyperplane.cpp new file mode 100644 index 0000000..f3f85e1 --- /dev/null +++ b/eigen/test/eigen2/eigen2_hyperplane.cpp @@ -0,0 +1,126 @@ +// 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 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// 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 <Eigen/Geometry> +#include <Eigen/LU> +#include <Eigen/QR> + +template<typename HyperplaneType> void hyperplane(const HyperplaneType& _plane) +{ + /* this test covers the following files: + Hyperplane.h + */ + + const int dim = _plane.dim(); + typedef typename HyperplaneType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, HyperplaneType::AmbientDimAtCompileTime, 1> VectorType; + typedef Matrix<Scalar, HyperplaneType::AmbientDimAtCompileTime, + HyperplaneType::AmbientDimAtCompileTime> MatrixType; + + VectorType p0 = VectorType::Random(dim); + VectorType p1 = VectorType::Random(dim); + + VectorType n0 = VectorType::Random(dim).normalized(); + VectorType n1 = VectorType::Random(dim).normalized(); + + HyperplaneType pl0(n0, p0); + HyperplaneType pl1(n1, p1); + HyperplaneType pl2 = pl1; + + Scalar s0 = ei_random<Scalar>(); + Scalar s1 = ei_random<Scalar>(); + + VERIFY_IS_APPROX( n1.eigen2_dot(n1), Scalar(1) ); + + VERIFY_IS_MUCH_SMALLER_THAN( pl0.absDistance(p0), Scalar(1) ); + VERIFY_IS_APPROX( pl1.signedDistance(p1 + n1 * s0), s0 ); + VERIFY_IS_MUCH_SMALLER_THAN( pl1.signedDistance(pl1.projection(p0)), Scalar(1) ); + VERIFY_IS_MUCH_SMALLER_THAN( pl1.absDistance(p1 + pl1.normal().unitOrthogonal() * s1), Scalar(1) ); + + // transform + if (!NumTraits<Scalar>::IsComplex) + { + MatrixType rot = MatrixType::Random(dim,dim).qr().matrixQ(); + Scaling<Scalar,HyperplaneType::AmbientDimAtCompileTime> scaling(VectorType::Random()); + Translation<Scalar,HyperplaneType::AmbientDimAtCompileTime> translation(VectorType::Random()); + + pl2 = pl1; + VERIFY_IS_MUCH_SMALLER_THAN( pl2.transform(rot).absDistance(rot * p1), Scalar(1) ); + pl2 = pl1; + VERIFY_IS_MUCH_SMALLER_THAN( pl2.transform(rot,Isometry).absDistance(rot * p1), Scalar(1) ); + pl2 = pl1; + VERIFY_IS_MUCH_SMALLER_THAN( pl2.transform(rot*scaling).absDistance((rot*scaling) * p1), Scalar(1) ); + pl2 = pl1; + VERIFY_IS_MUCH_SMALLER_THAN( pl2.transform(rot*scaling*translation) + .absDistance((rot*scaling*translation) * p1), Scalar(1) ); + pl2 = pl1; + VERIFY_IS_MUCH_SMALLER_THAN( pl2.transform(rot*translation,Isometry) + .absDistance((rot*translation) * p1), Scalar(1) ); + } + + // casting + const int Dim = HyperplaneType::AmbientDimAtCompileTime; + typedef typename GetDifferentType<Scalar>::type OtherScalar; + Hyperplane<OtherScalar,Dim> hp1f = pl1.template cast<OtherScalar>(); + VERIFY_IS_APPROX(hp1f.template cast<Scalar>(),pl1); + Hyperplane<Scalar,Dim> hp1d = pl1.template cast<Scalar>(); + VERIFY_IS_APPROX(hp1d.template cast<Scalar>(),pl1); +} + +template<typename Scalar> void lines() +{ + typedef Hyperplane<Scalar, 2> HLine; + typedef ParametrizedLine<Scalar, 2> PLine; + typedef Matrix<Scalar,2,1> Vector; + typedef Matrix<Scalar,3,1> CoeffsType; + + for(int i = 0; i < 10; i++) + { + Vector center = Vector::Random(); + Vector u = Vector::Random(); + Vector v = Vector::Random(); + Scalar a = ei_random<Scalar>(); + while (ei_abs(a-1) < 1e-4) a = ei_random<Scalar>(); + while (u.norm() < 1e-4) u = Vector::Random(); + while (v.norm() < 1e-4) v = Vector::Random(); + + HLine line_u = HLine::Through(center + u, center + a*u); + HLine line_v = HLine::Through(center + v, center + a*v); + + // the line equations should be normalized so that a^2+b^2=1 + VERIFY_IS_APPROX(line_u.normal().norm(), Scalar(1)); + VERIFY_IS_APPROX(line_v.normal().norm(), Scalar(1)); + + Vector result = line_u.intersection(line_v); + + // the lines should intersect at the point we called "center" + VERIFY_IS_APPROX(result, center); + + // check conversions between two types of lines + PLine pl(line_u); // gcc 3.3 will commit suicide if we don't name this variable + CoeffsType converted_coeffs(HLine(pl).coeffs()); + converted_coeffs *= line_u.coeffs()(0)/converted_coeffs(0); + VERIFY(line_u.coeffs().isApprox(converted_coeffs)); + } +} + +void test_eigen2_hyperplane() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( hyperplane(Hyperplane<float,2>()) ); + CALL_SUBTEST_2( hyperplane(Hyperplane<float,3>()) ); + CALL_SUBTEST_3( hyperplane(Hyperplane<double,4>()) ); + CALL_SUBTEST_4( hyperplane(Hyperplane<std::complex<double>,5>()) ); + CALL_SUBTEST_5( lines<float>() ); + CALL_SUBTEST_6( lines<double>() ); + } +} diff --git a/eigen/test/eigen2/eigen2_inverse.cpp b/eigen/test/eigen2/eigen2_inverse.cpp new file mode 100644 index 0000000..ccd24a1 --- /dev/null +++ b/eigen/test/eigen2/eigen2_inverse.cpp @@ -0,0 +1,62 @@ +// 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 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// 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 <Eigen/LU> + +template<typename MatrixType> void inverse(const MatrixType& m) +{ + /* this test covers the following files: + Inverse.h + */ + int rows = m.rows(); + int cols = m.cols(); + + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType; + + MatrixType m1 = MatrixType::Random(rows, cols), + m2(rows, cols), + identity = MatrixType::Identity(rows, rows); + + while(ei_abs(m1.determinant()) < RealScalar(0.1) && rows <= 8) + { + m1 = MatrixType::Random(rows, cols); + } + + m2 = m1.inverse(); + VERIFY_IS_APPROX(m1, m2.inverse() ); + + m1.computeInverse(&m2); + VERIFY_IS_APPROX(m1, m2.inverse() ); + + VERIFY_IS_APPROX((Scalar(2)*m2).inverse(), m2.inverse()*Scalar(0.5)); + + VERIFY_IS_APPROX(identity, m1.inverse() * m1 ); + VERIFY_IS_APPROX(identity, m1 * m1.inverse() ); + + VERIFY_IS_APPROX(m1, m1.inverse().inverse() ); + + // since for the general case we implement separately row-major and col-major, test that + VERIFY_IS_APPROX(m1.transpose().inverse(), m1.inverse().transpose()); +} + +void test_eigen2_inverse() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( inverse(Matrix<double,1,1>()) ); + CALL_SUBTEST_2( inverse(Matrix2d()) ); + CALL_SUBTEST_3( inverse(Matrix3f()) ); + CALL_SUBTEST_4( inverse(Matrix4f()) ); + CALL_SUBTEST_5( inverse(MatrixXf(8,8)) ); + CALL_SUBTEST_6( inverse(MatrixXcd(7,7)) ); + } +} diff --git a/eigen/test/eigen2/eigen2_linearstructure.cpp b/eigen/test/eigen2/eigen2_linearstructure.cpp new file mode 100644 index 0000000..488f4c4 --- /dev/null +++ b/eigen/test/eigen2/eigen2_linearstructure.cpp @@ -0,0 +1,83 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" + +template<typename MatrixType> void linearStructure(const MatrixType& m) +{ + /* this test covers the following files: + Sum.h Difference.h Opposite.h ScalarMultiple.h + */ + + typedef typename MatrixType::Scalar Scalar; + typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; + + int rows = m.rows(); + int cols = m.cols(); + + // this test relies a lot on Random.h, and there's not much more that we can do + // to test it, hence I consider that we will have tested Random.h + MatrixType m1 = MatrixType::Random(rows, cols), + m2 = MatrixType::Random(rows, cols), + m3(rows, cols); + + Scalar s1 = ei_random<Scalar>(); + while (ei_abs(s1)<1e-3) s1 = ei_random<Scalar>(); + + int r = ei_random<int>(0, rows-1), + c = ei_random<int>(0, cols-1); + + VERIFY_IS_APPROX(-(-m1), m1); + VERIFY_IS_APPROX(m1+m1, 2*m1); + VERIFY_IS_APPROX(m1+m2-m1, m2); + VERIFY_IS_APPROX(-m2+m1+m2, m1); + VERIFY_IS_APPROX(m1*s1, s1*m1); + VERIFY_IS_APPROX((m1+m2)*s1, s1*m1+s1*m2); + VERIFY_IS_APPROX((-m1+m2)*s1, -s1*m1+s1*m2); + m3 = m2; m3 += m1; + VERIFY_IS_APPROX(m3, m1+m2); + m3 = m2; m3 -= m1; + VERIFY_IS_APPROX(m3, m2-m1); + m3 = m2; m3 *= s1; + VERIFY_IS_APPROX(m3, s1*m2); + if(NumTraits<Scalar>::HasFloatingPoint) + { + m3 = m2; m3 /= s1; + VERIFY_IS_APPROX(m3, m2/s1); + } + + // again, test operator() to check const-qualification + VERIFY_IS_APPROX((-m1)(r,c), -(m1(r,c))); + VERIFY_IS_APPROX((m1-m2)(r,c), (m1(r,c))-(m2(r,c))); + VERIFY_IS_APPROX((m1+m2)(r,c), (m1(r,c))+(m2(r,c))); + VERIFY_IS_APPROX((s1*m1)(r,c), s1*(m1(r,c))); + VERIFY_IS_APPROX((m1*s1)(r,c), (m1(r,c))*s1); + if(NumTraits<Scalar>::HasFloatingPoint) + VERIFY_IS_APPROX((m1/s1)(r,c), (m1(r,c))/s1); + + // use .block to disable vectorization and compare to the vectorized version + VERIFY_IS_APPROX(m1+m1.block(0,0,rows,cols), m1+m1); + VERIFY_IS_APPROX(m1.cwise() * m1.block(0,0,rows,cols), m1.cwise() * m1); + VERIFY_IS_APPROX(m1 - m1.block(0,0,rows,cols), m1 - m1); + VERIFY_IS_APPROX(m1.block(0,0,rows,cols) * s1, m1 * s1); +} + +void test_eigen2_linearstructure() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( linearStructure(Matrix<float, 1, 1>()) ); + CALL_SUBTEST_2( linearStructure(Matrix2f()) ); + CALL_SUBTEST_3( linearStructure(Vector3d()) ); + CALL_SUBTEST_4( linearStructure(Matrix4d()) ); + CALL_SUBTEST_5( linearStructure(MatrixXcf(3, 3)) ); + CALL_SUBTEST_6( linearStructure(MatrixXf(8, 12)) ); + CALL_SUBTEST_7( linearStructure(MatrixXi(8, 12)) ); + CALL_SUBTEST_8( linearStructure(MatrixXcd(20, 20)) ); + } +} diff --git a/eigen/test/eigen2/eigen2_lu.cpp b/eigen/test/eigen2/eigen2_lu.cpp new file mode 100644 index 0000000..e993b1c --- /dev/null +++ b/eigen/test/eigen2/eigen2_lu.cpp @@ -0,0 +1,122 @@ +// 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 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// 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 <Eigen/LU> + +template<typename Derived> +void doSomeRankPreservingOperations(Eigen::MatrixBase<Derived>& m) +{ + typedef typename Derived::RealScalar RealScalar; + for(int a = 0; a < 3*(m.rows()+m.cols()); a++) + { + RealScalar d = Eigen::ei_random<RealScalar>(-1,1); + int i = Eigen::ei_random<int>(0,m.rows()-1); // i is a random row number + int j; + do { + j = Eigen::ei_random<int>(0,m.rows()-1); + } while (i==j); // j is another one (must be different) + m.row(i) += d * m.row(j); + + i = Eigen::ei_random<int>(0,m.cols()-1); // i is a random column number + do { + j = Eigen::ei_random<int>(0,m.cols()-1); + } while (i==j); // j is another one (must be different) + m.col(i) += d * m.col(j); + } +} + +template<typename MatrixType> void lu_non_invertible() +{ + /* this test covers the following files: + LU.h + */ + // NOTE there seems to be a problem with too small sizes -- could easily lie in the doSomeRankPreservingOperations function + int rows = ei_random<int>(20,200), cols = ei_random<int>(20,200), cols2 = ei_random<int>(20,200); + int rank = ei_random<int>(1, std::min(rows, cols)-1); + + MatrixType m1(rows, cols), m2(cols, cols2), m3(rows, cols2), k(1,1); + m1 = MatrixType::Random(rows,cols); + if(rows <= cols) + for(int i = rank; i < rows; i++) m1.row(i).setZero(); + else + for(int i = rank; i < cols; i++) m1.col(i).setZero(); + doSomeRankPreservingOperations(m1); + + LU<MatrixType> lu(m1); + typename LU<MatrixType>::KernelResultType m1kernel = lu.kernel(); + typename LU<MatrixType>::ImageResultType m1image = lu.image(); + + VERIFY(rank == lu.rank()); + VERIFY(cols - lu.rank() == lu.dimensionOfKernel()); + VERIFY(!lu.isInjective()); + VERIFY(!lu.isInvertible()); + VERIFY(lu.isSurjective() == (lu.rank() == rows)); + VERIFY((m1 * m1kernel).isMuchSmallerThan(m1)); + VERIFY(m1image.lu().rank() == rank); + MatrixType sidebyside(m1.rows(), m1.cols() + m1image.cols()); + sidebyside << m1, m1image; + VERIFY(sidebyside.lu().rank() == rank); + m2 = MatrixType::Random(cols,cols2); + m3 = m1*m2; + m2 = MatrixType::Random(cols,cols2); + lu.solve(m3, &m2); + VERIFY_IS_APPROX(m3, m1*m2); + /* solve now always returns true + m3 = MatrixType::Random(rows,cols2); + VERIFY(!lu.solve(m3, &m2)); + */ +} + +template<typename MatrixType> void lu_invertible() +{ + /* this test covers the following files: + LU.h + */ + typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; + int size = ei_random<int>(10,200); + + MatrixType m1(size, size), m2(size, size), m3(size, size); + m1 = MatrixType::Random(size,size); + + if (ei_is_same_type<RealScalar,float>::ret) + { + // let's build a matrix more stable to inverse + MatrixType a = MatrixType::Random(size,size*2); + m1 += a * a.adjoint(); + } + + LU<MatrixType> lu(m1); + VERIFY(0 == lu.dimensionOfKernel()); + VERIFY(size == lu.rank()); + VERIFY(lu.isInjective()); + VERIFY(lu.isSurjective()); + VERIFY(lu.isInvertible()); + VERIFY(lu.image().lu().isInvertible()); + m3 = MatrixType::Random(size,size); + lu.solve(m3, &m2); + VERIFY_IS_APPROX(m3, m1*m2); + VERIFY_IS_APPROX(m2, lu.inverse()*m3); + m3 = MatrixType::Random(size,size); + VERIFY(lu.solve(m3, &m2)); +} + +void test_eigen2_lu() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( lu_non_invertible<MatrixXf>() ); + CALL_SUBTEST_2( lu_non_invertible<MatrixXd>() ); + CALL_SUBTEST_3( lu_non_invertible<MatrixXcf>() ); + CALL_SUBTEST_4( lu_non_invertible<MatrixXcd>() ); + CALL_SUBTEST_1( lu_invertible<MatrixXf>() ); + CALL_SUBTEST_2( lu_invertible<MatrixXd>() ); + CALL_SUBTEST_3( lu_invertible<MatrixXcf>() ); + CALL_SUBTEST_4( lu_invertible<MatrixXcd>() ); + } +} diff --git a/eigen/test/eigen2/eigen2_map.cpp b/eigen/test/eigen2/eigen2_map.cpp new file mode 100644 index 0000000..4a1c4e1 --- /dev/null +++ b/eigen/test/eigen2/eigen2_map.cpp @@ -0,0 +1,114 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2007-2010 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" + +template<typename VectorType> void map_class_vector(const VectorType& m) +{ + typedef typename VectorType::Scalar Scalar; + + int size = m.size(); + + // test Map.h + Scalar* array1 = ei_aligned_new<Scalar>(size); + Scalar* array2 = ei_aligned_new<Scalar>(size); + Scalar* array3 = new Scalar[size+1]; + Scalar* array3unaligned = std::size_t(array3)%16 == 0 ? array3+1 : array3; + + Map<VectorType, Aligned>(array1, size) = VectorType::Random(size); + Map<VectorType>(array2, size) = Map<VectorType>(array1, size); + Map<VectorType>(array3unaligned, size) = Map<VectorType>((const Scalar*)array1, size); // test non-const-correctness support in eigen2 + VectorType ma1 = Map<VectorType>(array1, size); + VectorType ma2 = Map<VectorType, Aligned>(array2, size); + VectorType ma3 = Map<VectorType>(array3unaligned, size); + VERIFY_IS_APPROX(ma1, ma2); + VERIFY_IS_APPROX(ma1, ma3); + + ei_aligned_delete(array1, size); + ei_aligned_delete(array2, size); + delete[] array3; +} + +template<typename MatrixType> void map_class_matrix(const MatrixType& m) +{ + typedef typename MatrixType::Scalar Scalar; + + int rows = m.rows(), cols = m.cols(), size = rows*cols; + + // test Map.h + Scalar* array1 = ei_aligned_new<Scalar>(size); + for(int i = 0; i < size; i++) array1[i] = Scalar(1); + Scalar* array2 = ei_aligned_new<Scalar>(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 = std::size_t(array3)%16 == 0 ? array3+1 : array3; + Map<MatrixType, Aligned>(array1, rows, cols) = MatrixType::Ones(rows,cols); + Map<MatrixType>(array2, rows, cols) = Map<MatrixType>((const Scalar*)array1, rows, cols); // test non-const-correctness support in eigen2 + Map<MatrixType>(array3unaligned, rows, cols) = Map<MatrixType>(array1, rows, cols); + MatrixType ma1 = Map<MatrixType>(array1, rows, cols); + MatrixType ma2 = Map<MatrixType, Aligned>(array2, rows, cols); + VERIFY_IS_APPROX(ma1, ma2); + MatrixType ma3 = Map<MatrixType>(array3unaligned, rows, cols); + VERIFY_IS_APPROX(ma1, ma3); + + ei_aligned_delete(array1, size); + ei_aligned_delete(array2, size); + delete[] array3; +} + +template<typename VectorType> void map_static_methods(const VectorType& m) +{ + typedef typename VectorType::Scalar Scalar; + + int size = m.size(); + + // test Map.h + Scalar* array1 = ei_aligned_new<Scalar>(size); + Scalar* array2 = ei_aligned_new<Scalar>(size); + Scalar* array3 = new Scalar[size+1]; + Scalar* array3unaligned = std::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_APPROX(ma1, ma2); + VERIFY_IS_APPROX(ma1, ma3); + + ei_aligned_delete(array1, size); + ei_aligned_delete(array2, size); + delete[] array3; +} + + +void test_eigen2_map() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( map_class_vector(Matrix<float, 1, 1>()) ); + CALL_SUBTEST_2( map_class_vector(Vector4d()) ); + 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_1( map_class_matrix(Matrix<float, 1, 1>()) ); + CALL_SUBTEST_2( map_class_matrix(Matrix4d()) ); + CALL_SUBTEST_6( map_class_matrix(Matrix<float,3,5>()) ); + CALL_SUBTEST_4( map_class_matrix(MatrixXcf(ei_random<int>(1,10),ei_random<int>(1,10))) ); + CALL_SUBTEST_5( map_class_matrix(MatrixXi(ei_random<int>(1,10),ei_random<int>(1,10))) ); + + CALL_SUBTEST_1( map_static_methods(Matrix<double, 1, 1>()) ); + CALL_SUBTEST_2( map_static_methods(Vector3f()) ); + CALL_SUBTEST_7( map_static_methods(RowVector3d()) ); + CALL_SUBTEST_4( map_static_methods(VectorXcd(8)) ); + CALL_SUBTEST_5( map_static_methods(VectorXf(12)) ); + } +} diff --git a/eigen/test/eigen2/eigen2_meta.cpp b/eigen/test/eigen2/eigen2_meta.cpp new file mode 100644 index 0000000..1d01bd8 --- /dev/null +++ b/eigen/test/eigen2/eigen2_meta.cpp @@ -0,0 +1,60 @@ +// 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 Gael Guennebaud <g.gael@free.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" + +void test_eigen2_meta() +{ + typedef float & FloatRef; + typedef const float & ConstFloatRef; + + VERIFY((ei_meta_if<(3<4),ei_meta_true, ei_meta_false>::ret::ret)); + VERIFY(( ei_is_same_type<float,float>::ret)); + VERIFY((!ei_is_same_type<float,double>::ret)); + VERIFY((!ei_is_same_type<float,float&>::ret)); + VERIFY((!ei_is_same_type<float,const float&>::ret)); + + VERIFY(( ei_is_same_type<float,ei_cleantype<const float&>::type >::ret)); + VERIFY(( ei_is_same_type<float,ei_cleantype<const float*>::type >::ret)); + VERIFY(( ei_is_same_type<float,ei_cleantype<const float*&>::type >::ret)); + VERIFY(( ei_is_same_type<float,ei_cleantype<float**>::type >::ret)); + VERIFY(( ei_is_same_type<float,ei_cleantype<float**&>::type >::ret)); + VERIFY(( ei_is_same_type<float,ei_cleantype<float* const *&>::type >::ret)); + VERIFY(( ei_is_same_type<float,ei_cleantype<float* const>::type >::ret)); + + VERIFY(( ei_is_same_type<float*,ei_unconst<const float*>::type >::ret)); + VERIFY(( ei_is_same_type<float&,ei_unconst<const float&>::type >::ret)); + VERIFY(( ei_is_same_type<float&,ei_unconst<ConstFloatRef>::type >::ret)); + + VERIFY(( ei_is_same_type<float&,ei_unconst<float&>::type >::ret)); + VERIFY(( ei_is_same_type<float,ei_unref<float&>::type >::ret)); + VERIFY(( ei_is_same_type<const float,ei_unref<const float&>::type >::ret)); + VERIFY(( ei_is_same_type<float,ei_unpointer<float*>::type >::ret)); + VERIFY(( ei_is_same_type<const float,ei_unpointer<const float*>::type >::ret)); + VERIFY(( ei_is_same_type<float,ei_unpointer<float* const >::type >::ret)); + + VERIFY(ei_meta_sqrt<1>::ret == 1); + #define VERIFY_META_SQRT(X) VERIFY(ei_meta_sqrt<X>::ret == int(ei_sqrt(double(X)))) + VERIFY_META_SQRT(2); + VERIFY_META_SQRT(3); + VERIFY_META_SQRT(4); + VERIFY_META_SQRT(5); + VERIFY_META_SQRT(6); + VERIFY_META_SQRT(8); + VERIFY_META_SQRT(9); + VERIFY_META_SQRT(15); + VERIFY_META_SQRT(16); + VERIFY_META_SQRT(17); + VERIFY_META_SQRT(255); + VERIFY_META_SQRT(256); + VERIFY_META_SQRT(257); + VERIFY_META_SQRT(1023); + VERIFY_META_SQRT(1024); + VERIFY_META_SQRT(1025); +} diff --git a/eigen/test/eigen2/eigen2_miscmatrices.cpp b/eigen/test/eigen2/eigen2_miscmatrices.cpp new file mode 100644 index 0000000..8bbb20c --- /dev/null +++ b/eigen/test/eigen2/eigen2_miscmatrices.cpp @@ -0,0 +1,48 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" + +template<typename MatrixType> void miscMatrices(const MatrixType& m) +{ + /* this test covers the following files: + DiagonalMatrix.h Ones.h + */ + + typedef typename MatrixType::Scalar Scalar; + typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; + typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType; + int rows = m.rows(); + int cols = m.cols(); + + int r = ei_random<int>(0, rows-1), r2 = ei_random<int>(0, rows-1), c = ei_random<int>(0, cols-1); + VERIFY_IS_APPROX(MatrixType::Ones(rows,cols)(r,c), static_cast<Scalar>(1)); + MatrixType m1 = MatrixType::Ones(rows,cols); + VERIFY_IS_APPROX(m1(r,c), static_cast<Scalar>(1)); + VectorType v1 = VectorType::Random(rows); + v1[0]; + Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> + square = v1.asDiagonal(); + if(r==r2) VERIFY_IS_APPROX(square(r,r2), v1[r]); + else VERIFY_IS_MUCH_SMALLER_THAN(square(r,r2), static_cast<Scalar>(1)); + square = MatrixType::Zero(rows, rows); + square.diagonal() = VectorType::Ones(rows); + VERIFY_IS_APPROX(square, MatrixType::Identity(rows, rows)); +} + +void test_eigen2_miscmatrices() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( miscMatrices(Matrix<float, 1, 1>()) ); + CALL_SUBTEST_2( miscMatrices(Matrix4d()) ); + CALL_SUBTEST_3( miscMatrices(MatrixXcf(3, 3)) ); + CALL_SUBTEST_4( miscMatrices(MatrixXi(8, 12)) ); + CALL_SUBTEST_5( miscMatrices(MatrixXcd(20, 20)) ); + } +} diff --git a/eigen/test/eigen2/eigen2_mixingtypes.cpp b/eigen/test/eigen2/eigen2_mixingtypes.cpp new file mode 100644 index 0000000..fb5ac5d --- /dev/null +++ b/eigen/test/eigen2/eigen2_mixingtypes.cpp @@ -0,0 +1,77 @@ +// 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 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// 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 + +#ifndef EIGEN_DONT_VECTORIZE +#define EIGEN_DONT_VECTORIZE // SSE intrinsics aren't designed to allow mixing types +#endif + +#include "main.h" + + +template<int SizeAtCompileType> void mixingtypes(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); + 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; + 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; + md*mcd; + mcd*md; + mf*vcf; + mcf*vf; + mcf *= mf; + vcd = md*vcd; + vcf = mcf*vf; +#if 0 + // these are know generating hard build errors in eigen3 + VERIFY_RAISES_ASSERT(mf*md); + VERIFY_RAISES_ASSERT(mcf*mcd); + VERIFY_RAISES_ASSERT(mcf*vcd); + VERIFY_RAISES_ASSERT(vcf = mf*vf); + + vf.eigen2_dot(vf); + VERIFY_RAISES_ASSERT(vd.eigen2_dot(vf)); + VERIFY_RAISES_ASSERT(vcf.eigen2_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. +#endif +} + +void test_eigen2_mixingtypes() +{ + // check that our operator new is indeed called: + CALL_SUBTEST_1(mixingtypes<3>()); + CALL_SUBTEST_2(mixingtypes<4>()); + CALL_SUBTEST_3(mixingtypes<Dynamic>(20)); +} diff --git a/eigen/test/eigen2/eigen2_newstdvector.cpp b/eigen/test/eigen2/eigen2_newstdvector.cpp new file mode 100644 index 0000000..5f90099 --- /dev/null +++ b/eigen/test/eigen2/eigen2_newstdvector.cpp @@ -0,0 +1,149 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// 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/. + +#define EIGEN_USE_NEW_STDVECTOR +#include "main.h" +#include <Eigen/StdVector> +#include <Eigen/Geometry> + +template<typename MatrixType> +void check_stdvector_matrix(const MatrixType& m) +{ + int rows = m.rows(); + int cols = m.cols(); + MatrixType x = MatrixType::Random(rows,cols), y = MatrixType::Random(rows,cols); + std::vector<MatrixType,Eigen::aligned_allocator<MatrixType> > v(10, MatrixType(rows,cols)), w(20, y); + v[5] = x; + w[6] = v[5]; + VERIFY_IS_APPROX(w[6], v[5]); + v = w; + for(int i = 0; i < 20; i++) + { + VERIFY_IS_APPROX(w[i], v[i]); + } + + v.resize(21); + v[20] = x; + VERIFY_IS_APPROX(v[20], x); + v.resize(22,y); + VERIFY_IS_APPROX(v[21], y); + v.push_back(x); + VERIFY_IS_APPROX(v[22], x); + VERIFY((std::size_t)&(v[22]) == (std::size_t)&(v[21]) + sizeof(MatrixType)); + + // do a lot of push_back such that the vector gets internally resized + // (with memory reallocation) + MatrixType* ref = &w[0]; + for(int i=0; i<30 || ((ref==&w[0]) && i<300); ++i) + v.push_back(w[i%w.size()]); + for(unsigned int i=23; i<v.size(); ++i) + { + VERIFY(v[i]==w[(i-23)%w.size()]); + } +} + +template<typename TransformType> +void check_stdvector_transform(const TransformType&) +{ + typedef typename TransformType::MatrixType MatrixType; + TransformType x(MatrixType::Random()), y(MatrixType::Random()); + std::vector<TransformType,Eigen::aligned_allocator<TransformType> > v(10), w(20, y); + v[5] = x; + w[6] = v[5]; + VERIFY_IS_APPROX(w[6], v[5]); + v = w; + for(int i = 0; i < 20; i++) + { + VERIFY_IS_APPROX(w[i], v[i]); + } + + v.resize(21); + v[20] = x; + VERIFY_IS_APPROX(v[20], x); + v.resize(22,y); + VERIFY_IS_APPROX(v[21], y); + v.push_back(x); + VERIFY_IS_APPROX(v[22], x); + VERIFY((std::size_t)&(v[22]) == (std::size_t)&(v[21]) + sizeof(TransformType)); + + // do a lot of push_back such that the vector gets internally resized + // (with memory reallocation) + TransformType* ref = &w[0]; + for(int i=0; i<30 || ((ref==&w[0]) && i<300); ++i) + v.push_back(w[i%w.size()]); + for(unsigned int i=23; i<v.size(); ++i) + { + VERIFY(v[i].matrix()==w[(i-23)%w.size()].matrix()); + } +} + +template<typename QuaternionType> +void check_stdvector_quaternion(const QuaternionType&) +{ + typedef typename QuaternionType::Coefficients Coefficients; + QuaternionType x(Coefficients::Random()), y(Coefficients::Random()); + std::vector<QuaternionType,Eigen::aligned_allocator<QuaternionType> > v(10), w(20, y); + v[5] = x; + w[6] = v[5]; + VERIFY_IS_APPROX(w[6], v[5]); + v = w; + for(int i = 0; i < 20; i++) + { + VERIFY_IS_APPROX(w[i], v[i]); + } + + v.resize(21); + v[20] = x; + VERIFY_IS_APPROX(v[20], x); + v.resize(22,y); + VERIFY_IS_APPROX(v[21], y); + v.push_back(x); + VERIFY_IS_APPROX(v[22], x); + VERIFY((std::size_t)&(v[22]) == (std::size_t)&(v[21]) + sizeof(QuaternionType)); + + // do a lot of push_back such that the vector gets internally resized + // (with memory reallocation) + QuaternionType* ref = &w[0]; + for(int i=0; i<30 || ((ref==&w[0]) && i<300); ++i) + v.push_back(w[i%w.size()]); + for(unsigned int i=23; i<v.size(); ++i) + { + VERIFY(v[i].coeffs()==w[(i-23)%w.size()].coeffs()); + } +} + +void test_eigen2_newstdvector() +{ + // some non vectorizable fixed sizes + CALL_SUBTEST_1(check_stdvector_matrix(Vector2f())); + CALL_SUBTEST_1(check_stdvector_matrix(Matrix3f())); + CALL_SUBTEST_1(check_stdvector_matrix(Matrix3d())); + + // some vectorizable fixed sizes + CALL_SUBTEST_2(check_stdvector_matrix(Matrix2f())); + CALL_SUBTEST_2(check_stdvector_matrix(Vector4f())); + CALL_SUBTEST_2(check_stdvector_matrix(Matrix4f())); + CALL_SUBTEST_2(check_stdvector_matrix(Matrix4d())); + + // some dynamic sizes + CALL_SUBTEST_3(check_stdvector_matrix(MatrixXd(1,1))); + CALL_SUBTEST_3(check_stdvector_matrix(VectorXd(20))); + CALL_SUBTEST_3(check_stdvector_matrix(RowVectorXf(20))); + CALL_SUBTEST_3(check_stdvector_matrix(MatrixXcf(10,10))); + + // some Transform + CALL_SUBTEST_4(check_stdvector_transform(Transform2f())); + CALL_SUBTEST_4(check_stdvector_transform(Transform3f())); + CALL_SUBTEST_4(check_stdvector_transform(Transform3d())); + //CALL_SUBTEST(check_stdvector_transform(Transform4d())); + + // some Quaternion + CALL_SUBTEST_5(check_stdvector_quaternion(Quaternionf())); + CALL_SUBTEST_5(check_stdvector_quaternion(Quaterniond())); +} diff --git a/eigen/test/eigen2/eigen2_nomalloc.cpp b/eigen/test/eigen2/eigen2_nomalloc.cpp new file mode 100644 index 0000000..d34a699 --- /dev/null +++ b/eigen/test/eigen2/eigen2_nomalloc.cpp @@ -0,0 +1,53 @@ +// 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 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// 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 hack is needed to make this file compiles with -pedantic (gcc) +#ifdef __GNUC__ +#define throw(X) +#endif +// discard stack allocation as that too bypasses malloc +#define EIGEN_STACK_ALLOCATION_LIMIT 0 +// any heap allocation will raise an assert +#define EIGEN_NO_MALLOC + +#include "main.h" + +template<typename MatrixType> void nomalloc(const MatrixType& m) +{ + /* this test check no dynamic memory allocation are issued with fixed-size matrices + */ + + typedef typename MatrixType::Scalar Scalar; + + int rows = m.rows(); + int cols = m.cols(); + + MatrixType m1 = MatrixType::Random(rows, cols), + m2 = MatrixType::Random(rows, cols); + + Scalar s1 = ei_random<Scalar>(); + + int r = ei_random<int>(0, rows-1), + c = ei_random<int>(0, cols-1); + + VERIFY_IS_APPROX((m1+m2)*s1, s1*m1+s1*m2); + VERIFY_IS_APPROX((m1+m2)(r,c), (m1(r,c))+(m2(r,c))); + VERIFY_IS_APPROX(m1.cwise() * m1.block(0,0,rows,cols), m1.cwise() * m1); + VERIFY_IS_APPROX((m1*m1.transpose())*m2, m1*(m1.transpose()*m2)); +} + +void test_eigen2_nomalloc() +{ + // check that our operator new is indeed called: + VERIFY_RAISES_ASSERT(MatrixXd dummy = MatrixXd::Random(3,3)); + CALL_SUBTEST_1( nomalloc(Matrix<float, 1, 1>()) ); + CALL_SUBTEST_2( nomalloc(Matrix4d()) ); + CALL_SUBTEST_3( nomalloc(Matrix<float,32,32>()) ); +} diff --git a/eigen/test/eigen2/eigen2_packetmath.cpp b/eigen/test/eigen2/eigen2_packetmath.cpp new file mode 100644 index 0000000..b1f325f --- /dev/null +++ b/eigen/test/eigen2/eigen2_packetmath.cpp @@ -0,0 +1,132 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// 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" + +// using namespace Eigen; + +template<typename Scalar> bool areApprox(const Scalar* a, const Scalar* b, int size) +{ + for (int i=0; i<size; ++i) + if (!ei_isApprox(a[i],b[i])) return false; + return true; +} + +#define CHECK_CWISE(REFOP, POP) { \ + for (int i=0; i<PacketSize; ++i) \ + ref[i] = REFOP(data1[i], data1[i+PacketSize]); \ + ei_pstore(data2, POP(ei_pload(data1), ei_pload(data1+PacketSize))); \ + VERIFY(areApprox(ref, data2, PacketSize) && #POP); \ +} + +#define REF_ADD(a,b) ((a)+(b)) +#define REF_SUB(a,b) ((a)-(b)) +#define REF_MUL(a,b) ((a)*(b)) +#define REF_DIV(a,b) ((a)/(b)) + +namespace std { + +template<> const complex<float>& min(const complex<float>& a, const complex<float>& b) +{ return a.real() < b.real() ? a : b; } + +template<> const complex<float>& max(const complex<float>& a, const complex<float>& b) +{ return a.real() < b.real() ? b : a; } + +} + +template<typename Scalar> void packetmath() +{ + typedef typename ei_packet_traits<Scalar>::type Packet; + const int PacketSize = ei_packet_traits<Scalar>::size; + + const int size = PacketSize*4; + EIGEN_ALIGN_128 Scalar data1[ei_packet_traits<Scalar>::size*4]; + EIGEN_ALIGN_128 Scalar data2[ei_packet_traits<Scalar>::size*4]; + EIGEN_ALIGN_128 Packet packets[PacketSize*2]; + EIGEN_ALIGN_128 Scalar ref[ei_packet_traits<Scalar>::size*4]; + for (int i=0; i<size; ++i) + { + data1[i] = ei_random<Scalar>(); + data2[i] = ei_random<Scalar>(); + } + + ei_pstore(data2, ei_pload(data1)); + VERIFY(areApprox(data1, data2, PacketSize) && "aligned load/store"); + + for (int offset=0; offset<PacketSize; ++offset) + { + ei_pstore(data2, ei_ploadu(data1+offset)); + VERIFY(areApprox(data1+offset, data2, PacketSize) && "ei_ploadu"); + } + + for (int offset=0; offset<PacketSize; ++offset) + { + ei_pstoreu(data2+offset, ei_pload(data1)); + VERIFY(areApprox(data1, data2+offset, PacketSize) && "ei_pstoreu"); + } + + for (int offset=0; offset<PacketSize; ++offset) + { + packets[0] = ei_pload(data1); + packets[1] = ei_pload(data1+PacketSize); + if (offset==0) ei_palign<0>(packets[0], packets[1]); + else if (offset==1) ei_palign<1>(packets[0], packets[1]); + else if (offset==2) ei_palign<2>(packets[0], packets[1]); + else if (offset==3) ei_palign<3>(packets[0], packets[1]); + ei_pstore(data2, packets[0]); + + for (int i=0; i<PacketSize; ++i) + ref[i] = data1[i+offset]; + + typedef Matrix<Scalar, PacketSize, 1> Vector; + VERIFY(areApprox(ref, data2, PacketSize) && "ei_palign"); + } + + CHECK_CWISE(REF_ADD, ei_padd); + CHECK_CWISE(REF_SUB, ei_psub); + CHECK_CWISE(REF_MUL, ei_pmul); + #ifndef EIGEN_VECTORIZE_ALTIVEC + if (!ei_is_same_type<Scalar,int>::ret) + CHECK_CWISE(REF_DIV, ei_pdiv); + #endif + CHECK_CWISE(std::min, ei_pmin); + CHECK_CWISE(std::max, ei_pmax); + + for (int i=0; i<PacketSize; ++i) + ref[i] = data1[0]; + ei_pstore(data2, ei_pset1(data1[0])); + VERIFY(areApprox(ref, data2, PacketSize) && "ei_pset1"); + + VERIFY(ei_isApprox(data1[0], ei_pfirst(ei_pload(data1))) && "ei_pfirst"); + + ref[0] = 0; + for (int i=0; i<PacketSize; ++i) + ref[0] += data1[i]; + VERIFY(ei_isApprox(ref[0], ei_predux(ei_pload(data1))) && "ei_predux"); + + for (int j=0; j<PacketSize; ++j) + { + ref[j] = 0; + for (int i=0; i<PacketSize; ++i) + ref[j] += data1[i+j*PacketSize]; + packets[j] = ei_pload(data1+j*PacketSize); + } + ei_pstore(data2, ei_preduxp(packets)); + VERIFY(areApprox(ref, data2, PacketSize) && "ei_preduxp"); +} + +void test_eigen2_packetmath() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( packetmath<float>() ); + CALL_SUBTEST_2( packetmath<double>() ); + CALL_SUBTEST_3( packetmath<int>() ); + CALL_SUBTEST_4( packetmath<std::complex<float> >() ); + } +} diff --git a/eigen/test/eigen2/eigen2_parametrizedline.cpp b/eigen/test/eigen2/eigen2_parametrizedline.cpp new file mode 100644 index 0000000..8147288 --- /dev/null +++ b/eigen/test/eigen2/eigen2_parametrizedline.cpp @@ -0,0 +1,62 @@ +// 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 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// 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 <Eigen/Geometry> +#include <Eigen/LU> +#include <Eigen/QR> + +template<typename LineType> void parametrizedline(const LineType& _line) +{ + /* this test covers the following files: + ParametrizedLine.h + */ + + const int dim = _line.dim(); + typedef typename LineType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, LineType::AmbientDimAtCompileTime, 1> VectorType; + typedef Matrix<Scalar, LineType::AmbientDimAtCompileTime, + LineType::AmbientDimAtCompileTime> MatrixType; + + VectorType p0 = VectorType::Random(dim); + VectorType p1 = VectorType::Random(dim); + + VectorType d0 = VectorType::Random(dim).normalized(); + + LineType l0(p0, d0); + + Scalar s0 = ei_random<Scalar>(); + Scalar s1 = ei_abs(ei_random<Scalar>()); + + VERIFY_IS_MUCH_SMALLER_THAN( l0.distance(p0), RealScalar(1) ); + VERIFY_IS_MUCH_SMALLER_THAN( l0.distance(p0+s0*d0), RealScalar(1) ); + VERIFY_IS_APPROX( (l0.projection(p1)-p1).norm(), l0.distance(p1) ); + VERIFY_IS_MUCH_SMALLER_THAN( l0.distance(l0.projection(p1)), RealScalar(1) ); + VERIFY_IS_APPROX( Scalar(l0.distance((p0+s0*d0) + d0.unitOrthogonal() * s1)), s1 ); + + // casting + const int Dim = LineType::AmbientDimAtCompileTime; + typedef typename GetDifferentType<Scalar>::type OtherScalar; + ParametrizedLine<OtherScalar,Dim> hp1f = l0.template cast<OtherScalar>(); + VERIFY_IS_APPROX(hp1f.template cast<Scalar>(),l0); + ParametrizedLine<Scalar,Dim> hp1d = l0.template cast<Scalar>(); + VERIFY_IS_APPROX(hp1d.template cast<Scalar>(),l0); +} + +void test_eigen2_parametrizedline() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( parametrizedline(ParametrizedLine<float,2>()) ); + CALL_SUBTEST_2( parametrizedline(ParametrizedLine<float,3>()) ); + CALL_SUBTEST_3( parametrizedline(ParametrizedLine<double,4>()) ); + CALL_SUBTEST_4( parametrizedline(ParametrizedLine<std::complex<double>,5>()) ); + } +} diff --git a/eigen/test/eigen2/eigen2_prec_inverse_4x4.cpp b/eigen/test/eigen2/eigen2_prec_inverse_4x4.cpp new file mode 100644 index 0000000..8bfa556 --- /dev/null +++ b/eigen/test/eigen2/eigen2_prec_inverse_4x4.cpp @@ -0,0 +1,84 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// 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 <Eigen/LU> +#include <algorithm> + +template<typename T> std::string type_name() { return "other"; } +template<> std::string type_name<float>() { return "float"; } +template<> std::string type_name<double>() { return "double"; } +template<> std::string type_name<int>() { return "int"; } +template<> std::string type_name<std::complex<float> >() { return "complex<float>"; } +template<> std::string type_name<std::complex<double> >() { return "complex<double>"; } +template<> std::string type_name<std::complex<int> >() { return "complex<int>"; } + +#define EIGEN_DEBUG_VAR(x) std::cerr << #x << " = " << x << std::endl; + +template<typename T> inline typename NumTraits<T>::Real epsilon() +{ + return std::numeric_limits<typename NumTraits<T>::Real>::epsilon(); +} + +template<typename MatrixType> void inverse_permutation_4x4() +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + Vector4i indices(0,1,2,3); + for(int i = 0; i < 24; ++i) + { + MatrixType m = MatrixType::Zero(); + m(indices(0),0) = 1; + m(indices(1),1) = 1; + m(indices(2),2) = 1; + m(indices(3),3) = 1; + MatrixType inv = m.inverse(); + double error = double( (m*inv-MatrixType::Identity()).norm() / epsilon<Scalar>() ); + VERIFY(error == 0.0); + std::next_permutation(indices.data(),indices.data()+4); + } +} + +template<typename MatrixType> void inverse_general_4x4(int repeat) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + double error_sum = 0., error_max = 0.; + for(int i = 0; i < repeat; ++i) + { + MatrixType m; + RealScalar absdet; + do { + m = MatrixType::Random(); + absdet = ei_abs(m.determinant()); + } while(absdet < 10 * epsilon<Scalar>()); + MatrixType inv = m.inverse(); + double error = double( (m*inv-MatrixType::Identity()).norm() * absdet / epsilon<Scalar>() ); + error_sum += error; + error_max = std::max(error_max, error); + } + std::cerr << "inverse_general_4x4, Scalar = " << type_name<Scalar>() << std::endl; + double error_avg = error_sum / repeat; + EIGEN_DEBUG_VAR(error_avg); + EIGEN_DEBUG_VAR(error_max); + VERIFY(error_avg < (NumTraits<Scalar>::IsComplex ? 8.0 : 1.25)); + VERIFY(error_max < (NumTraits<Scalar>::IsComplex ? 64.0 : 20.0)); +} + +void test_eigen2_prec_inverse_4x4() +{ + CALL_SUBTEST_1((inverse_permutation_4x4<Matrix4f>())); + CALL_SUBTEST_1(( inverse_general_4x4<Matrix4f>(200000 * g_repeat) )); + + CALL_SUBTEST_2((inverse_permutation_4x4<Matrix<double,4,4,RowMajor> >())); + CALL_SUBTEST_2(( inverse_general_4x4<Matrix<double,4,4,RowMajor> >(200000 * g_repeat) )); + + CALL_SUBTEST_3((inverse_permutation_4x4<Matrix4cf>())); + CALL_SUBTEST_3((inverse_general_4x4<Matrix4cf>(50000 * g_repeat))); +} diff --git a/eigen/test/eigen2/eigen2_product_large.cpp b/eigen/test/eigen2/eigen2_product_large.cpp new file mode 100644 index 0000000..5149ef7 --- /dev/null +++ b/eigen/test/eigen2/eigen2_product_large.cpp @@ -0,0 +1,45 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// 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 "product.h" + +void test_eigen2_product_large() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( product(MatrixXf(ei_random<int>(1,320), ei_random<int>(1,320))) ); + CALL_SUBTEST_2( product(MatrixXd(ei_random<int>(1,320), ei_random<int>(1,320))) ); + CALL_SUBTEST_3( product(MatrixXi(ei_random<int>(1,320), ei_random<int>(1,320))) ); + CALL_SUBTEST_4( product(MatrixXcf(ei_random<int>(1,50), ei_random<int>(1,50))) ); + CALL_SUBTEST_5( product(Matrix<float,Dynamic,Dynamic,RowMajor>(ei_random<int>(1,320), ei_random<int>(1,320))) ); + } + +#ifdef EIGEN_TEST_PART_6 + { + // test a specific issue in DiagonalProduct + int N = 1000000; + VectorXf v = VectorXf::Ones(N); + MatrixXf m = MatrixXf::Ones(N,3); + m = (v+v).asDiagonal() * m; + VERIFY_IS_APPROX(m, MatrixXf::Constant(N,3,2)); + } + + { + // test deferred resizing in Matrix::operator= + MatrixXf a = MatrixXf::Random(10,4), b = MatrixXf::Random(4,10), c = a; + VERIFY_IS_APPROX((a = a * b), (c * b).eval()); + } + + { + MatrixXf mat1(10,10); mat1.setRandom(); + MatrixXf mat2(32,10); mat2.setRandom(); + MatrixXf result = mat1.row(2)*mat2.transpose(); + VERIFY_IS_APPROX(result, (mat1.row(2)*mat2.transpose()).eval()); + } +#endif +} diff --git a/eigen/test/eigen2/eigen2_product_small.cpp b/eigen/test/eigen2/eigen2_product_small.cpp new file mode 100644 index 0000000..4cd8c10 --- /dev/null +++ b/eigen/test/eigen2/eigen2_product_small.cpp @@ -0,0 +1,22 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// 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/. + +#define EIGEN_NO_STATIC_ASSERT +#include "product.h" + +void test_eigen2_product_small() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( product(Matrix<float, 3, 2>()) ); + CALL_SUBTEST_2( product(Matrix<int, 3, 5>()) ); + CALL_SUBTEST_3( product(Matrix3d()) ); + CALL_SUBTEST_4( product(Matrix4d()) ); + CALL_SUBTEST_5( product(Matrix4f()) ); + } +} diff --git a/eigen/test/eigen2/eigen2_qr.cpp b/eigen/test/eigen2/eigen2_qr.cpp new file mode 100644 index 0000000..76977e4 --- /dev/null +++ b/eigen/test/eigen2/eigen2_qr.cpp @@ -0,0 +1,69 @@ +// 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 Gael Guennebaud <g.gael@free.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" +#include <Eigen/QR> + +template<typename MatrixType> void qr(const MatrixType& m) +{ + /* this test covers the following files: + QR.h + */ + int rows = m.rows(); + int cols = m.cols(); + + typedef typename MatrixType::Scalar Scalar; + typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> SquareMatrixType; + typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType; + + MatrixType a = MatrixType::Random(rows,cols); + QR<MatrixType> qrOfA(a); + VERIFY_IS_APPROX(a, qrOfA.matrixQ() * qrOfA.matrixR()); + VERIFY_IS_NOT_APPROX(a+MatrixType::Identity(rows, cols), qrOfA.matrixQ() * qrOfA.matrixR()); + + #if 0 // eigenvalues module not yet ready + SquareMatrixType b = a.adjoint() * a; + + // check tridiagonalization + Tridiagonalization<SquareMatrixType> tridiag(b); + VERIFY_IS_APPROX(b, tridiag.matrixQ() * tridiag.matrixT() * tridiag.matrixQ().adjoint()); + + // check hessenberg decomposition + HessenbergDecomposition<SquareMatrixType> hess(b); + VERIFY_IS_APPROX(b, hess.matrixQ() * hess.matrixH() * hess.matrixQ().adjoint()); + VERIFY_IS_APPROX(tridiag.matrixT(), hess.matrixH()); + b = SquareMatrixType::Random(cols,cols); + hess.compute(b); + VERIFY_IS_APPROX(b, hess.matrixQ() * hess.matrixH() * hess.matrixQ().adjoint()); + #endif +} + +void test_eigen2_qr() +{ + for(int i = 0; i < 1; i++) { + CALL_SUBTEST_1( qr(Matrix2f()) ); + CALL_SUBTEST_2( qr(Matrix4d()) ); + CALL_SUBTEST_3( qr(MatrixXf(12,8)) ); + CALL_SUBTEST_4( qr(MatrixXcd(5,5)) ); + CALL_SUBTEST_4( qr(MatrixXcd(7,3)) ); + } + +#ifdef EIGEN_TEST_PART_5 + // small isFullRank test + { + Matrix3d mat; + mat << 1, 45, 1, 2, 2, 2, 1, 2, 3; + VERIFY(mat.qr().isFullRank()); + mat << 1, 1, 1, 2, 2, 2, 1, 2, 3; + //always returns true in eigen2support + //VERIFY(!mat.qr().isFullRank()); + } + +#endif +} diff --git a/eigen/test/eigen2/eigen2_qtvector.cpp b/eigen/test/eigen2/eigen2_qtvector.cpp new file mode 100644 index 0000000..6cfb58a --- /dev/null +++ b/eigen/test/eigen2/eigen2_qtvector.cpp @@ -0,0 +1,158 @@ +// 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 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// 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/. + +#define EIGEN_WORK_AROUND_QT_BUG_CALLING_WRONG_OPERATOR_NEW_FIXED_IN_QT_4_5 + +#include "main.h" + +#include <Eigen/Geometry> +#include <Eigen/QtAlignedMalloc> + +#include <QtCore/QVector> + +template<typename MatrixType> +void check_qtvector_matrix(const MatrixType& m) +{ + int rows = m.rows(); + int cols = m.cols(); + MatrixType x = MatrixType::Random(rows,cols), y = MatrixType::Random(rows,cols); + QVector<MatrixType> v(10, MatrixType(rows,cols)), w(20, y); + for(int i = 0; i < 20; i++) + { + VERIFY_IS_APPROX(w[i], y); + } + v[5] = x; + w[6] = v[5]; + VERIFY_IS_APPROX(w[6], v[5]); + v = w; + for(int i = 0; i < 20; i++) + { + VERIFY_IS_APPROX(w[i], v[i]); + } + + v.resize(21); + v[20] = x; + VERIFY_IS_APPROX(v[20], x); + v.fill(y,22); + VERIFY_IS_APPROX(v[21], y); + v.push_back(x); + VERIFY_IS_APPROX(v[22], x); + VERIFY((size_t)&(v[22]) == (size_t)&(v[21]) + sizeof(MatrixType)); + + // do a lot of push_back such that the vector gets internally resized + // (with memory reallocation) + MatrixType* ref = &w[0]; + for(int i=0; i<30 || ((ref==&w[0]) && i<300); ++i) + v.push_back(w[i%w.size()]); + for(int i=23; i<v.size(); ++i) + { + VERIFY(v[i]==w[(i-23)%w.size()]); + } +} + +template<typename TransformType> +void check_qtvector_transform(const TransformType&) +{ + typedef typename TransformType::MatrixType MatrixType; + TransformType x(MatrixType::Random()), y(MatrixType::Random()); + QVector<TransformType> v(10), w(20, y); + v[5] = x; + w[6] = v[5]; + VERIFY_IS_APPROX(w[6], v[5]); + v = w; + for(int i = 0; i < 20; i++) + { + VERIFY_IS_APPROX(w[i], v[i]); + } + + v.resize(21); + v[20] = x; + VERIFY_IS_APPROX(v[20], x); + v.fill(y,22); + VERIFY_IS_APPROX(v[21], y); + v.push_back(x); + VERIFY_IS_APPROX(v[22], x); + VERIFY((size_t)&(v[22]) == (size_t)&(v[21]) + sizeof(TransformType)); + + // do a lot of push_back such that the vector gets internally resized + // (with memory reallocation) + TransformType* ref = &w[0]; + for(int i=0; i<30 || ((ref==&w[0]) && i<300); ++i) + v.push_back(w[i%w.size()]); + for(unsigned int i=23; int(i)<v.size(); ++i) + { + VERIFY(v[i].matrix()==w[(i-23)%w.size()].matrix()); + } +} + +template<typename QuaternionType> +void check_qtvector_quaternion(const QuaternionType&) +{ + typedef typename QuaternionType::Coefficients Coefficients; + QuaternionType x(Coefficients::Random()), y(Coefficients::Random()); + QVector<QuaternionType> v(10), w(20, y); + v[5] = x; + w[6] = v[5]; + VERIFY_IS_APPROX(w[6], v[5]); + v = w; + for(int i = 0; i < 20; i++) + { + VERIFY_IS_APPROX(w[i], v[i]); + } + + v.resize(21); + v[20] = x; + VERIFY_IS_APPROX(v[20], x); + v.fill(y,22); + VERIFY_IS_APPROX(v[21], y); + v.push_back(x); + VERIFY_IS_APPROX(v[22], x); + VERIFY((size_t)&(v[22]) == (size_t)&(v[21]) + sizeof(QuaternionType)); + + // do a lot of push_back such that the vector gets internally resized + // (with memory reallocation) + QuaternionType* ref = &w[0]; + for(int i=0; i<30 || ((ref==&w[0]) && i<300); ++i) + v.push_back(w[i%w.size()]); + for(unsigned int i=23; int(i)<v.size(); ++i) + { + VERIFY(v[i].coeffs()==w[(i-23)%w.size()].coeffs()); + } +} + +void test_eigen2_qtvector() +{ + // some non vectorizable fixed sizes + CALL_SUBTEST_1(check_qtvector_matrix(Vector2f())); + CALL_SUBTEST_1(check_qtvector_matrix(Matrix3f())); + CALL_SUBTEST_1(check_qtvector_matrix(Matrix3d())); + + // some vectorizable fixed sizes + CALL_SUBTEST_2(check_qtvector_matrix(Matrix2f())); + CALL_SUBTEST_2(check_qtvector_matrix(Vector4f())); + CALL_SUBTEST_2(check_qtvector_matrix(Matrix4f())); + CALL_SUBTEST_2(check_qtvector_matrix(Matrix4d())); + + // some dynamic sizes + CALL_SUBTEST_3(check_qtvector_matrix(MatrixXd(1,1))); + CALL_SUBTEST_3(check_qtvector_matrix(VectorXd(20))); + CALL_SUBTEST_3(check_qtvector_matrix(RowVectorXf(20))); + CALL_SUBTEST_3(check_qtvector_matrix(MatrixXcf(10,10))); + + // some Transform + CALL_SUBTEST_4(check_qtvector_transform(Transform2f())); + CALL_SUBTEST_4(check_qtvector_transform(Transform3f())); + CALL_SUBTEST_4(check_qtvector_transform(Transform3d())); + //CALL_SUBTEST_4(check_qtvector_transform(Transform4d())); + + // some Quaternion + CALL_SUBTEST_5(check_qtvector_quaternion(Quaternionf())); + CALL_SUBTEST_5(check_qtvector_quaternion(Quaternionf())); +} diff --git a/eigen/test/eigen2/eigen2_regression.cpp b/eigen/test/eigen2/eigen2_regression.cpp new file mode 100644 index 0000000..c68b58d --- /dev/null +++ b/eigen/test/eigen2/eigen2_regression.cpp @@ -0,0 +1,136 @@ +// 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 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// 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 <Eigen/LeastSquares> + +template<typename VectorType, + typename HyperplaneType> +void makeNoisyCohyperplanarPoints(int numPoints, + VectorType **points, + HyperplaneType *hyperplane, + typename VectorType::Scalar noiseAmplitude) +{ + typedef typename VectorType::Scalar Scalar; + const int size = points[0]->size(); + // pick a random hyperplane, store the coefficients of its equation + hyperplane->coeffs().resize(size + 1); + for(int j = 0; j < size + 1; j++) + { + do { + hyperplane->coeffs().coeffRef(j) = ei_random<Scalar>(); + } while(ei_abs(hyperplane->coeffs().coeff(j)) < 0.5); + } + + // now pick numPoints random points on this hyperplane + for(int i = 0; i < numPoints; i++) + { + VectorType& cur_point = *(points[i]); + do + { + cur_point = VectorType::Random(size)/*.normalized()*/; + // project cur_point onto the hyperplane + Scalar x = - (hyperplane->coeffs().start(size).cwise()*cur_point).sum(); + cur_point *= hyperplane->coeffs().coeff(size) / x; + } while( cur_point.norm() < 0.5 + || cur_point.norm() > 2.0 ); + } + + // add some noise to these points + for(int i = 0; i < numPoints; i++ ) + *(points[i]) += noiseAmplitude * VectorType::Random(size); +} + +template<typename VectorType> +void check_linearRegression(int numPoints, + VectorType **points, + const VectorType& original, + typename VectorType::Scalar tolerance) +{ + int size = points[0]->size(); + assert(size==2); + VectorType result(size); + linearRegression(numPoints, points, &result, 1); + typename VectorType::Scalar error = (result - original).norm() / original.norm(); + VERIFY(ei_abs(error) < ei_abs(tolerance)); +} + +template<typename VectorType, + typename HyperplaneType> +void check_fitHyperplane(int numPoints, + VectorType **points, + const HyperplaneType& original, + typename VectorType::Scalar tolerance) +{ + int size = points[0]->size(); + HyperplaneType result(size); + fitHyperplane(numPoints, points, &result); + result.coeffs() *= original.coeffs().coeff(size)/result.coeffs().coeff(size); + typename VectorType::Scalar error = (result.coeffs() - original.coeffs()).norm() / original.coeffs().norm(); + std::cout << ei_abs(error) << " xxx " << ei_abs(tolerance) << std::endl; + VERIFY(ei_abs(error) < ei_abs(tolerance)); +} + +void test_eigen2_regression() +{ + for(int i = 0; i < g_repeat; i++) + { +#ifdef EIGEN_TEST_PART_1 + { + Vector2f points2f [1000]; + Vector2f *points2f_ptrs [1000]; + for(int i = 0; i < 1000; i++) points2f_ptrs[i] = &(points2f[i]); + Vector2f coeffs2f; + Hyperplane<float,2> coeffs3f; + makeNoisyCohyperplanarPoints(1000, points2f_ptrs, &coeffs3f, 0.01f); + coeffs2f[0] = -coeffs3f.coeffs()[0]/coeffs3f.coeffs()[1]; + coeffs2f[1] = -coeffs3f.coeffs()[2]/coeffs3f.coeffs()[1]; + CALL_SUBTEST(check_linearRegression(10, points2f_ptrs, coeffs2f, 0.05f)); + CALL_SUBTEST(check_linearRegression(100, points2f_ptrs, coeffs2f, 0.01f)); + CALL_SUBTEST(check_linearRegression(1000, points2f_ptrs, coeffs2f, 0.002f)); + } +#endif +#ifdef EIGEN_TEST_PART_2 + { + Vector2f points2f [1000]; + Vector2f *points2f_ptrs [1000]; + for(int i = 0; i < 1000; i++) points2f_ptrs[i] = &(points2f[i]); + Hyperplane<float,2> coeffs3f; + makeNoisyCohyperplanarPoints(1000, points2f_ptrs, &coeffs3f, 0.01f); + CALL_SUBTEST(check_fitHyperplane(10, points2f_ptrs, coeffs3f, 0.05f)); + CALL_SUBTEST(check_fitHyperplane(100, points2f_ptrs, coeffs3f, 0.01f)); + CALL_SUBTEST(check_fitHyperplane(1000, points2f_ptrs, coeffs3f, 0.002f)); + } +#endif +#ifdef EIGEN_TEST_PART_3 + { + Vector4d points4d [1000]; + Vector4d *points4d_ptrs [1000]; + for(int i = 0; i < 1000; i++) points4d_ptrs[i] = &(points4d[i]); + Hyperplane<double,4> coeffs5d; + makeNoisyCohyperplanarPoints(1000, points4d_ptrs, &coeffs5d, 0.01); + CALL_SUBTEST(check_fitHyperplane(10, points4d_ptrs, coeffs5d, 0.05)); + CALL_SUBTEST(check_fitHyperplane(100, points4d_ptrs, coeffs5d, 0.01)); + CALL_SUBTEST(check_fitHyperplane(1000, points4d_ptrs, coeffs5d, 0.002)); + } +#endif +#ifdef EIGEN_TEST_PART_4 + { + VectorXcd *points11cd_ptrs[1000]; + for(int i = 0; i < 1000; i++) points11cd_ptrs[i] = new VectorXcd(11); + Hyperplane<std::complex<double>,Dynamic> *coeffs12cd = new Hyperplane<std::complex<double>,Dynamic>(11); + makeNoisyCohyperplanarPoints(1000, points11cd_ptrs, coeffs12cd, 0.01); + CALL_SUBTEST(check_fitHyperplane(100, points11cd_ptrs, *coeffs12cd, 0.025)); + CALL_SUBTEST(check_fitHyperplane(1000, points11cd_ptrs, *coeffs12cd, 0.006)); + delete coeffs12cd; + for(int i = 0; i < 1000; i++) delete points11cd_ptrs[i]; + } +#endif + } +} diff --git a/eigen/test/eigen2/eigen2_sizeof.cpp b/eigen/test/eigen2/eigen2_sizeof.cpp new file mode 100644 index 0000000..ec1af5a --- /dev/null +++ b/eigen/test/eigen2/eigen2_sizeof.cpp @@ -0,0 +1,31 @@ +// 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 Gael Guennebaud <g.gael@free.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" + +template<typename MatrixType> void verifySizeOf(const MatrixType&) +{ + typedef typename MatrixType::Scalar Scalar; + if (MatrixType::RowsAtCompileTime!=Dynamic && MatrixType::ColsAtCompileTime!=Dynamic) + VERIFY(sizeof(MatrixType)==sizeof(Scalar)*MatrixType::SizeAtCompileTime); + else + VERIFY(sizeof(MatrixType)==sizeof(Scalar*) + 2 * sizeof(typename MatrixType::Index)); +} + +void test_eigen2_sizeof() +{ + CALL_SUBTEST( verifySizeOf(Matrix<float, 1, 1>()) ); + CALL_SUBTEST( verifySizeOf(Matrix4d()) ); + CALL_SUBTEST( verifySizeOf(Matrix<double, 4, 2>()) ); + CALL_SUBTEST( verifySizeOf(Matrix<bool, 7, 5>()) ); + CALL_SUBTEST( verifySizeOf(MatrixXcf(3, 3)) ); + CALL_SUBTEST( verifySizeOf(MatrixXi(8, 12)) ); + CALL_SUBTEST( verifySizeOf(MatrixXcd(20, 20)) ); + CALL_SUBTEST( verifySizeOf(Matrix<float, 100, 100>()) ); +} diff --git a/eigen/test/eigen2/eigen2_smallvectors.cpp b/eigen/test/eigen2/eigen2_smallvectors.cpp new file mode 100644 index 0000000..03962b1 --- /dev/null +++ b/eigen/test/eigen2/eigen2_smallvectors.cpp @@ -0,0 +1,42 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" + +template<typename Scalar> void smallVectors() +{ + typedef Matrix<Scalar, 1, 2> V2; + typedef Matrix<Scalar, 3, 1> V3; + typedef Matrix<Scalar, 1, 4> V4; + Scalar x1 = ei_random<Scalar>(), + x2 = ei_random<Scalar>(), + x3 = ei_random<Scalar>(), + x4 = ei_random<Scalar>(); + V2 v2(x1, x2); + V3 v3(x1, x2, x3); + V4 v4(x1, x2, x3, x4); + VERIFY_IS_APPROX(x1, v2.x()); + VERIFY_IS_APPROX(x1, v3.x()); + VERIFY_IS_APPROX(x1, v4.x()); + VERIFY_IS_APPROX(x2, v2.y()); + VERIFY_IS_APPROX(x2, v3.y()); + VERIFY_IS_APPROX(x2, v4.y()); + VERIFY_IS_APPROX(x3, v3.z()); + VERIFY_IS_APPROX(x3, v4.z()); + VERIFY_IS_APPROX(x4, v4.w()); +} + +void test_eigen2_smallvectors() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST( smallVectors<int>() ); + CALL_SUBTEST( smallVectors<float>() ); + CALL_SUBTEST( smallVectors<double>() ); + } +} diff --git a/eigen/test/eigen2/eigen2_sparse_basic.cpp b/eigen/test/eigen2/eigen2_sparse_basic.cpp new file mode 100644 index 0000000..0490776 --- /dev/null +++ b/eigen/test/eigen2/eigen2_sparse_basic.cpp @@ -0,0 +1,317 @@ +// 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 Daniel Gomez Ferro <dgomezferro@gmail.com> +// +// 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 "sparse.h" + +template<typename SetterType,typename DenseType, typename Scalar, int Options> +bool test_random_setter(SparseMatrix<Scalar,Options>& sm, const DenseType& ref, const std::vector<Vector2i>& nonzeroCoords) +{ + typedef SparseMatrix<Scalar,Options> SparseType; + { + sm.setZero(); + SetterType w(sm); + std::vector<Vector2i> remaining = nonzeroCoords; + while(!remaining.empty()) + { + int i = ei_random<int>(0,remaining.size()-1); + w(remaining[i].x(),remaining[i].y()) = ref.coeff(remaining[i].x(),remaining[i].y()); + remaining[i] = remaining.back(); + remaining.pop_back(); + } + } + return sm.isApprox(ref); +} + +template<typename SetterType,typename DenseType, typename T> +bool test_random_setter(DynamicSparseMatrix<T>& sm, const DenseType& ref, const std::vector<Vector2i>& nonzeroCoords) +{ + sm.setZero(); + std::vector<Vector2i> remaining = nonzeroCoords; + while(!remaining.empty()) + { + int i = ei_random<int>(0,remaining.size()-1); + sm.coeffRef(remaining[i].x(),remaining[i].y()) = ref.coeff(remaining[i].x(),remaining[i].y()); + remaining[i] = remaining.back(); + remaining.pop_back(); + } + return sm.isApprox(ref); +} + +template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& ref) +{ + const int rows = ref.rows(); + const int cols = ref.cols(); + typedef typename SparseMatrixType::Scalar Scalar; + enum { Flags = SparseMatrixType::Flags }; + + double density = std::max(8./(rows*cols), 0.01); + typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; + typedef Matrix<Scalar,Dynamic,1> DenseVector; + Scalar eps = 1e-6; + + SparseMatrixType m(rows, cols); + DenseMatrix refMat = DenseMatrix::Zero(rows, cols); + DenseVector vec1 = DenseVector::Random(rows); + Scalar s1 = ei_random<Scalar>(); + + std::vector<Vector2i> zeroCoords; + std::vector<Vector2i> nonzeroCoords; + initSparse<Scalar>(density, refMat, m, 0, &zeroCoords, &nonzeroCoords); + + if (zeroCoords.size()==0 || nonzeroCoords.size()==0) + return; + + // test coeff and coeffRef + for (int i=0; i<(int)zeroCoords.size(); ++i) + { + VERIFY_IS_MUCH_SMALLER_THAN( m.coeff(zeroCoords[i].x(),zeroCoords[i].y()), eps ); + if(ei_is_same_type<SparseMatrixType,SparseMatrix<Scalar,Flags> >::ret) + 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) + { + int j = ei_random<int>(0,cols-1); + int i = ei_random<int>(0,rows-1); + int w = ei_random<int>(1,cols-j-1); + int h = ei_random<int>(1,rows-i-1); + +// VERIFY_IS_APPROX(m.block(i,j,h,w), refMat.block(i,j,h,w)); + for(int c=0; c<w; c++) + { + VERIFY_IS_APPROX(m.block(i,j,h,w).col(c), refMat.block(i,j,h,w).col(c)); + for(int r=0; r<h; r++) + { +// VERIFY_IS_APPROX(m.block(i,j,h,w).col(c).coeff(r), refMat.block(i,j,h,w).col(c).coeff(r)); + } + } +// for(int r=0; r<h; r++) +// { +// VERIFY_IS_APPROX(m.block(i,j,h,w).row(r), refMat.block(i,j,h,w).row(r)); +// for(int c=0; c<w; c++) +// { +// VERIFY_IS_APPROX(m.block(i,j,h,w).row(r).coeff(c), refMat.block(i,j,h,w).row(r).coeff(c)); +// } +// } + } + + for(int c=0; c<cols; c++) + { + VERIFY_IS_APPROX(m.col(c) + m.col(c), (m + m).col(c)); + VERIFY_IS_APPROX(m.col(c) + m.col(c), refMat.col(c) + refMat.col(c)); + } + + for(int r=0; r<rows; r++) + { + VERIFY_IS_APPROX(m.row(r) + m.row(r), (m + m).row(r)); + VERIFY_IS_APPROX(m.row(r) + m.row(r), refMat.row(r) + refMat.row(r)); + } + */ + + // test SparseSetters + // coherent setter + // TODO extend the MatrixSetter +// { +// m.setZero(); +// VERIFY_IS_NOT_APPROX(m, refMat); +// SparseSetter<SparseMatrixType, FullyCoherentAccessPattern> w(m); +// for (int i=0; i<nonzeroCoords.size(); ++i) +// { +// w->coeffRef(nonzeroCoords[i].x(),nonzeroCoords[i].y()) = refMat.coeff(nonzeroCoords[i].x(),nonzeroCoords[i].y()); +// } +// } +// VERIFY_IS_APPROX(m, refMat); + + // random setter +// { +// m.setZero(); +// VERIFY_IS_NOT_APPROX(m, refMat); +// SparseSetter<SparseMatrixType, RandomAccessPattern> w(m); +// std::vector<Vector2i> remaining = nonzeroCoords; +// while(!remaining.empty()) +// { +// int i = ei_random<int>(0,remaining.size()-1); +// w->coeffRef(remaining[i].x(),remaining[i].y()) = refMat.coeff(remaining[i].x(),remaining[i].y()); +// remaining[i] = remaining.back(); +// remaining.pop_back(); +// } +// } +// VERIFY_IS_APPROX(m, refMat); + + VERIFY(( test_random_setter<RandomSetter<SparseMatrixType, StdMapTraits> >(m,refMat,nonzeroCoords) )); + #ifdef EIGEN_UNORDERED_MAP_SUPPORT + VERIFY(( test_random_setter<RandomSetter<SparseMatrixType, StdUnorderedMapTraits> >(m,refMat,nonzeroCoords) )); + #endif + #ifdef _DENSE_HASH_MAP_H_ + VERIFY(( test_random_setter<RandomSetter<SparseMatrixType, GoogleDenseHashMapTraits> >(m,refMat,nonzeroCoords) )); + #endif + #ifdef _SPARSE_HASH_MAP_H_ + VERIFY(( test_random_setter<RandomSetter<SparseMatrixType, GoogleSparseHashMapTraits> >(m,refMat,nonzeroCoords) )); + #endif + + // test fillrand + { + DenseMatrix m1(rows,cols); + m1.setZero(); + SparseMatrixType m2(rows,cols); + m2.startFill(); + for (int j=0; j<cols; ++j) + { + for (int k=0; k<rows/2; ++k) + { + int i = ei_random<int>(0,rows-1); + if (m1.coeff(i,j)==Scalar(0)) + m2.fillrand(i,j) = m1(i,j) = ei_random<Scalar>(); + } + } + m2.endFill(); + VERIFY_IS_APPROX(m2,m1); + } + + // test RandomSetter + /*{ + SparseMatrixType m1(rows,cols), m2(rows,cols); + DenseMatrix refM1 = DenseMatrix::Zero(rows, rows); + initSparse<Scalar>(density, refM1, m1); + { + Eigen::RandomSetter<SparseMatrixType > setter(m2); + for (int j=0; j<m1.outerSize(); ++j) + for (typename SparseMatrixType::InnerIterator i(m1,j); i; ++i) + setter(i.index(), j) = i.value(); + } + VERIFY_IS_APPROX(m1, m2); + }*/ +// std::cerr << m.transpose() << "\n\n" << refMat.transpose() << "\n\n"; +// VERIFY_IS_APPROX(m, refMat); + + // test basic computations + { + DenseMatrix refM1 = DenseMatrix::Zero(rows, rows); + DenseMatrix refM2 = DenseMatrix::Zero(rows, rows); + DenseMatrix refM3 = DenseMatrix::Zero(rows, rows); + DenseMatrix refM4 = DenseMatrix::Zero(rows, rows); + SparseMatrixType m1(rows, rows); + SparseMatrixType m2(rows, rows); + SparseMatrixType m3(rows, rows); + SparseMatrixType m4(rows, rows); + initSparse<Scalar>(density, refM1, m1); + initSparse<Scalar>(density, refM2, m2); + initSparse<Scalar>(density, refM3, m3); + initSparse<Scalar>(density, refM4, m4); + + VERIFY_IS_APPROX(m1+m2, refM1+refM2); + VERIFY_IS_APPROX(m1+m2+m3, refM1+refM2+refM3); + VERIFY_IS_APPROX(m3.cwise()*(m1+m2), refM3.cwise()*(refM1+refM2)); + VERIFY_IS_APPROX(m1*s1-m2, refM1*s1-refM2); + + VERIFY_IS_APPROX(m1*=s1, refM1*=s1); + VERIFY_IS_APPROX(m1/=s1, refM1/=s1); + + VERIFY_IS_APPROX(m1+=m2, refM1+=refM2); + VERIFY_IS_APPROX(m1-=m2, refM1-=refM2); + + VERIFY_IS_APPROX(m1.col(0).eigen2_dot(refM2.row(0)), refM1.col(0).eigen2_dot(refM2.row(0))); + + refM4.setRandom(); + // sparse cwise* dense + VERIFY_IS_APPROX(m3.cwise()*refM4, refM3.cwise()*refM4); +// VERIFY_IS_APPROX(m3.cwise()/refM4, refM3.cwise()/refM4); + } + + // test innerVector() + { + DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); + SparseMatrixType m2(rows, rows); + initSparse<Scalar>(density, refMat2, m2); + int j0 = ei_random(0,rows-1); + int j1 = ei_random(0,rows-1); + VERIFY_IS_APPROX(m2.innerVector(j0), refMat2.col(j0)); + VERIFY_IS_APPROX(m2.innerVector(j0)+m2.innerVector(j1), refMat2.col(j0)+refMat2.col(j1)); + //m2.innerVector(j0) = 2*m2.innerVector(j1); + //refMat2.col(j0) = 2*refMat2.col(j1); + //VERIFY_IS_APPROX(m2, refMat2); + } + + // test innerVectors() + { + DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); + SparseMatrixType m2(rows, rows); + initSparse<Scalar>(density, refMat2, m2); + int j0 = ei_random(0,rows-2); + int j1 = ei_random(0,rows-2); + int n0 = ei_random<int>(1,rows-std::max(j0,j1)); + VERIFY_IS_APPROX(m2.innerVectors(j0,n0), refMat2.block(0,j0,rows,n0)); + VERIFY_IS_APPROX(m2.innerVectors(j0,n0)+m2.innerVectors(j1,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); + } + + // test transpose + { + DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); + SparseMatrixType m2(rows, rows); + initSparse<Scalar>(density, refMat2, m2); + VERIFY_IS_APPROX(m2.transpose().eval(), refMat2.transpose().eval()); + VERIFY_IS_APPROX(m2.transpose(), refMat2.transpose()); + } + + // test prune + { + SparseMatrixType m2(rows, rows); + DenseMatrix refM2(rows, rows); + refM2.setZero(); + int countFalseNonZero = 0; + int countTrueNonZero = 0; + m2.startFill(); + for (int j=0; j<m2.outerSize(); ++j) + for (int i=0; i<m2.innerSize(); ++i) + { + float x = ei_random<float>(0,1); + if (x<0.1) + { + // do nothing + } + else if (x<0.5) + { + countFalseNonZero++; + m2.fill(i,j) = Scalar(0); + } + else + { + countTrueNonZero++; + m2.fill(i,j) = refM2(i,j) = Scalar(1); + } + } + m2.endFill(); + VERIFY(countFalseNonZero+countTrueNonZero == m2.nonZeros()); + VERIFY_IS_APPROX(m2, refM2); + m2.prune(1); + VERIFY(countTrueNonZero==m2.nonZeros()); + VERIFY_IS_APPROX(m2, refM2); + } +} + +void test_eigen2_sparse_basic() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( sparse_basic(SparseMatrix<double>(8, 8)) ); + CALL_SUBTEST_2( sparse_basic(SparseMatrix<std::complex<double> >(16, 16)) ); + CALL_SUBTEST_1( sparse_basic(SparseMatrix<double>(33, 33)) ); + + CALL_SUBTEST_3( sparse_basic(DynamicSparseMatrix<double>(8, 8)) ); + } +} diff --git a/eigen/test/eigen2/eigen2_sparse_product.cpp b/eigen/test/eigen2/eigen2_sparse_product.cpp new file mode 100644 index 0000000..d28e76d --- /dev/null +++ b/eigen/test/eigen2/eigen2_sparse_product.cpp @@ -0,0 +1,115 @@ +// 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 Daniel Gomez Ferro <dgomezferro@gmail.com> +// +// 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 "sparse.h" + +template<typename SparseMatrixType> void sparse_product(const SparseMatrixType& ref) +{ + const int rows = ref.rows(); + const int cols = ref.cols(); + typedef typename SparseMatrixType::Scalar Scalar; + enum { Flags = SparseMatrixType::Flags }; + + double density = std::max(8./(rows*cols), 0.01); + typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; + typedef Matrix<Scalar,Dynamic,1> DenseVector; + + // test matrix-matrix product + { + DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); + DenseMatrix refMat3 = DenseMatrix::Zero(rows, rows); + DenseMatrix refMat4 = DenseMatrix::Zero(rows, rows); + DenseMatrix dm4 = DenseMatrix::Zero(rows, rows); + SparseMatrixType m2(rows, rows); + SparseMatrixType m3(rows, rows); + SparseMatrixType m4(rows, rows); + initSparse<Scalar>(density, refMat2, m2); + initSparse<Scalar>(density, refMat3, m3); + initSparse<Scalar>(density, refMat4, m4); + VERIFY_IS_APPROX(m4=m2*m3, refMat4=refMat2*refMat3); + VERIFY_IS_APPROX(m4=m2.transpose()*m3, refMat4=refMat2.transpose()*refMat3); + VERIFY_IS_APPROX(m4=m2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose()); + VERIFY_IS_APPROX(m4=m2*m3.transpose(), refMat4=refMat2*refMat3.transpose()); + + // sparse * dense + VERIFY_IS_APPROX(dm4=m2*refMat3, refMat4=refMat2*refMat3); + VERIFY_IS_APPROX(dm4=m2*refMat3.transpose(), refMat4=refMat2*refMat3.transpose()); + VERIFY_IS_APPROX(dm4=m2.transpose()*refMat3, refMat4=refMat2.transpose()*refMat3); + VERIFY_IS_APPROX(dm4=m2.transpose()*refMat3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose()); + + // dense * sparse + VERIFY_IS_APPROX(dm4=refMat2*m3, refMat4=refMat2*refMat3); + VERIFY_IS_APPROX(dm4=refMat2*m3.transpose(), refMat4=refMat2*refMat3.transpose()); + VERIFY_IS_APPROX(dm4=refMat2.transpose()*m3, refMat4=refMat2.transpose()*refMat3); + VERIFY_IS_APPROX(dm4=refMat2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose()); + + VERIFY_IS_APPROX(m3=m3*m3, refMat3=refMat3*refMat3); + } + + // test matrix - diagonal product + if(false) // it compiles, but the precision is terrible. probably doesn't matter in this branch.... + { + DenseMatrix refM2 = DenseMatrix::Zero(rows, rows); + DenseMatrix refM3 = DenseMatrix::Zero(rows, rows); + DiagonalMatrix<DenseVector> d1(DenseVector::Random(rows)); + SparseMatrixType m2(rows, rows); + SparseMatrixType m3(rows, rows); + initSparse<Scalar>(density, refM2, m2); + initSparse<Scalar>(density, refM3, m3); + VERIFY_IS_APPROX(m3=m2*d1, refM3=refM2*d1); + VERIFY_IS_APPROX(m3=m2.transpose()*d1, refM3=refM2.transpose()*d1); + VERIFY_IS_APPROX(m3=d1*m2, refM3=d1*refM2); + VERIFY_IS_APPROX(m3=d1*m2.transpose(), refM3=d1 * refM2.transpose()); + } + + // test self adjoint products + { + DenseMatrix b = DenseMatrix::Random(rows, rows); + DenseMatrix x = DenseMatrix::Random(rows, rows); + DenseMatrix refX = DenseMatrix::Random(rows, rows); + DenseMatrix refUp = DenseMatrix::Zero(rows, rows); + DenseMatrix refLo = DenseMatrix::Zero(rows, rows); + DenseMatrix refS = DenseMatrix::Zero(rows, rows); + SparseMatrixType mUp(rows, rows); + SparseMatrixType mLo(rows, rows); + SparseMatrixType mS(rows, rows); + do { + initSparse<Scalar>(density, refUp, mUp, ForceRealDiag|/*ForceNonZeroDiag|*/MakeUpperTriangular); + } while (refUp.isZero()); + refLo = refUp.transpose().conjugate(); + mLo = mUp.transpose().conjugate(); + refS = refUp + refLo; + refS.diagonal() *= 0.5; + mS = mUp + mLo; + for (int k=0; k<mS.outerSize(); ++k) + for (typename SparseMatrixType::InnerIterator it(mS,k); it; ++it) + if (it.index() == k) + it.valueRef() *= 0.5; + + VERIFY_IS_APPROX(refS.adjoint(), refS); + VERIFY_IS_APPROX(mS.transpose().conjugate(), mS); + VERIFY_IS_APPROX(mS, refS); + VERIFY_IS_APPROX(x=mS*b, refX=refS*b); + VERIFY_IS_APPROX(x=mUp.template marked<UpperTriangular|SelfAdjoint>()*b, refX=refS*b); + VERIFY_IS_APPROX(x=mLo.template marked<LowerTriangular|SelfAdjoint>()*b, refX=refS*b); + VERIFY_IS_APPROX(x=mS.template marked<SelfAdjoint>()*b, refX=refS*b); + } + +} + +void test_eigen2_sparse_product() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( sparse_product(SparseMatrix<double>(8, 8)) ); + CALL_SUBTEST_2( sparse_product(SparseMatrix<std::complex<double> >(16, 16)) ); + CALL_SUBTEST_1( sparse_product(SparseMatrix<double>(33, 33)) ); + + CALL_SUBTEST_3( sparse_product(DynamicSparseMatrix<double>(8, 8)) ); + } +} diff --git a/eigen/test/eigen2/eigen2_sparse_solvers.cpp b/eigen/test/eigen2/eigen2_sparse_solvers.cpp new file mode 100644 index 0000000..3aef27a --- /dev/null +++ b/eigen/test/eigen2/eigen2_sparse_solvers.cpp @@ -0,0 +1,200 @@ +// 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 Daniel Gomez Ferro <dgomezferro@gmail.com> +// +// 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 "sparse.h" + +template<typename Scalar> void +initSPD(double density, + Matrix<Scalar,Dynamic,Dynamic>& refMat, + SparseMatrix<Scalar>& sparseMat) +{ + Matrix<Scalar,Dynamic,Dynamic> aux(refMat.rows(),refMat.cols()); + initSparse(density,refMat,sparseMat); + refMat = refMat * refMat.adjoint(); + for (int k=0; k<2; ++k) + { + initSparse(density,aux,sparseMat,ForceNonZeroDiag); + refMat += aux * aux.adjoint(); + } + sparseMat.startFill(); + for (int j=0 ; j<sparseMat.cols(); ++j) + for (int i=j ; i<sparseMat.rows(); ++i) + if (refMat(i,j)!=Scalar(0)) + sparseMat.fill(i,j) = refMat(i,j); + sparseMat.endFill(); +} + +template<typename Scalar> void sparse_solvers(int rows, int cols) +{ + double density = std::max(8./(rows*cols), 0.01); + typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; + typedef Matrix<Scalar,Dynamic,1> DenseVector; + // Scalar eps = 1e-6; + + DenseVector vec1 = DenseVector::Random(rows); + + std::vector<Vector2i> zeroCoords; + std::vector<Vector2i> nonzeroCoords; + + // test triangular solver + { + DenseVector vec2 = vec1, vec3 = vec1; + SparseMatrix<Scalar> m2(rows, cols); + DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols); + + // lower + initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); + VERIFY_IS_APPROX(refMat2.template marked<LowerTriangular>().solveTriangular(vec2), + m2.template marked<LowerTriangular>().solveTriangular(vec3)); + + // lower - transpose + initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); + VERIFY_IS_APPROX(refMat2.template marked<LowerTriangular>().transpose().solveTriangular(vec2), + m2.template marked<LowerTriangular>().transpose().solveTriangular(vec3)); + + // upper + initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords); + VERIFY_IS_APPROX(refMat2.template marked<UpperTriangular>().solveTriangular(vec2), + m2.template marked<UpperTriangular>().solveTriangular(vec3)); + + // upper - transpose + initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords); + VERIFY_IS_APPROX(refMat2.template marked<UpperTriangular>().transpose().solveTriangular(vec2), + m2.template marked<UpperTriangular>().transpose().solveTriangular(vec3)); + } + + // test LLT + { + // TODO fix the issue with complex (see SparseLLT::solveInPlace) + SparseMatrix<Scalar> m2(rows, cols); + DenseMatrix refMat2(rows, cols); + + DenseVector b = DenseVector::Random(cols); + DenseVector refX(cols), x(cols); + + initSPD(density, refMat2, m2); + + refMat2.llt().solve(b, &refX); + typedef SparseMatrix<Scalar,LowerTriangular|SelfAdjoint> SparseSelfAdjointMatrix; + if (!NumTraits<Scalar>::IsComplex) + { + x = b; + SparseLLT<SparseSelfAdjointMatrix> (m2).solveInPlace(x); + VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: default"); + } + #ifdef EIGEN_CHOLMOD_SUPPORT + x = b; + SparseLLT<SparseSelfAdjointMatrix,Cholmod>(m2).solveInPlace(x); + VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod"); + #endif + if (!NumTraits<Scalar>::IsComplex) + { + #ifdef EIGEN_TAUCS_SUPPORT + x = b; + SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,IncompleteFactorization).solveInPlace(x); + VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (IncompleteFactorization)"); + x = b; + SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,SupernodalMultifrontal).solveInPlace(x); + VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalMultifrontal)"); + x = b; + SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,SupernodalLeftLooking).solveInPlace(x); + VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalLeftLooking)"); + #endif + } + } + + // test LDLT + if (!NumTraits<Scalar>::IsComplex) + { + // TODO fix the issue with complex (see SparseLDLT::solveInPlace) + SparseMatrix<Scalar> m2(rows, cols); + DenseMatrix refMat2(rows, cols); + + DenseVector b = DenseVector::Random(cols); + DenseVector refX(cols), x(cols); + + //initSPD(density, refMat2, m2); + initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, 0, 0); + refMat2 += refMat2.adjoint(); + refMat2.diagonal() *= 0.5; + + refMat2.ldlt().solve(b, &refX); + typedef SparseMatrix<Scalar,UpperTriangular|SelfAdjoint> SparseSelfAdjointMatrix; + x = b; + SparseLDLT<SparseSelfAdjointMatrix> ldlt(m2); + if (ldlt.succeeded()) + ldlt.solveInPlace(x); + VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: default"); + } + + // test LU + { + static int count = 0; + SparseMatrix<Scalar> m2(rows, cols); + DenseMatrix refMat2(rows, cols); + + DenseVector b = DenseVector::Random(cols); + DenseVector refX(cols), x(cols); + + initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag, &zeroCoords, &nonzeroCoords); + + LU<DenseMatrix> refLu(refMat2); + refLu.solve(b, &refX); + #if defined(EIGEN_SUPERLU_SUPPORT) || defined(EIGEN_UMFPACK_SUPPORT) + Scalar refDet = refLu.determinant(); + #endif + x.setZero(); + // // SparseLU<SparseMatrix<Scalar> > (m2).solve(b,&x); + // // VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LU: default"); + #ifdef EIGEN_SUPERLU_SUPPORT + { + x.setZero(); + SparseLU<SparseMatrix<Scalar>,SuperLU> slu(m2); + if (slu.succeeded()) + { + if (slu.solve(b,&x)) { + VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LU: SuperLU"); + } + // std::cerr << refDet << " == " << slu.determinant() << "\n"; + if (count==0) { + VERIFY_IS_APPROX(refDet,slu.determinant()); // FIXME det is not very stable for complex + } + } + } + #endif + #ifdef EIGEN_UMFPACK_SUPPORT + { + // check solve + x.setZero(); + SparseLU<SparseMatrix<Scalar>,UmfPack> slu(m2); + if (slu.succeeded()) { + if (slu.solve(b,&x)) { + if (count==0) { + VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LU: umfpack"); // FIXME solve is not very stable for complex + } + } + VERIFY_IS_APPROX(refDet,slu.determinant()); + // TODO check the extracted data + //std::cerr << slu.matrixL() << "\n"; + } + } + #endif + count++; + } + +} + +void test_eigen2_sparse_solvers() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( sparse_solvers<double>(8, 8) ); + CALL_SUBTEST_2( sparse_solvers<std::complex<double> >(16, 16) ); + CALL_SUBTEST_1( sparse_solvers<double>(101, 101) ); + } +} diff --git a/eigen/test/eigen2/eigen2_sparse_vector.cpp b/eigen/test/eigen2/eigen2_sparse_vector.cpp new file mode 100644 index 0000000..e6d2d77 --- /dev/null +++ b/eigen/test/eigen2/eigen2_sparse_vector.cpp @@ -0,0 +1,84 @@ +// 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 Daniel Gomez Ferro <dgomezferro@gmail.com> +// +// 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 "sparse.h" + +template<typename Scalar> void sparse_vector(int rows, int cols) +{ + double densityMat = std::max(8./(rows*cols), 0.01); + double densityVec = std::max(8./float(rows), 0.1); + typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; + typedef Matrix<Scalar,Dynamic,1> DenseVector; + typedef SparseVector<Scalar> SparseVectorType; + typedef SparseMatrix<Scalar> SparseMatrixType; + Scalar eps = 1e-6; + + SparseMatrixType m1(rows,cols); + SparseVectorType v1(rows), v2(rows), v3(rows); + DenseMatrix refM1 = DenseMatrix::Zero(rows, cols); + DenseVector refV1 = DenseVector::Random(rows), + refV2 = DenseVector::Random(rows), + refV3 = DenseVector::Random(rows); + + std::vector<int> zerocoords, nonzerocoords; + initSparse<Scalar>(densityVec, refV1, v1, &zerocoords, &nonzerocoords); + initSparse<Scalar>(densityMat, refM1, m1); + + initSparse<Scalar>(densityVec, refV2, v2); + initSparse<Scalar>(densityVec, refV3, v3); + + Scalar s1 = ei_random<Scalar>(); + + // test coeff and coeffRef + for (unsigned int i=0; i<zerocoords.size(); ++i) + { + VERIFY_IS_MUCH_SMALLER_THAN( v1.coeff(zerocoords[i]), eps ); + //VERIFY_RAISES_ASSERT( v1.coeffRef(zerocoords[i]) = 5 ); + } + { + VERIFY(int(nonzerocoords.size()) == v1.nonZeros()); + int j=0; + for (typename SparseVectorType::InnerIterator it(v1); it; ++it,++j) + { + VERIFY(nonzerocoords[j]==it.index()); + VERIFY(it.value()==v1.coeff(it.index())); + VERIFY(it.value()==refV1.coeff(it.index())); + } + } + VERIFY_IS_APPROX(v1, refV1); + + v1.coeffRef(nonzerocoords[0]) = Scalar(5); + refV1.coeffRef(nonzerocoords[0]) = Scalar(5); + VERIFY_IS_APPROX(v1, refV1); + + VERIFY_IS_APPROX(v1+v2, refV1+refV2); + VERIFY_IS_APPROX(v1+v2+v3, refV1+refV2+refV3); + + VERIFY_IS_APPROX(v1*s1-v2, refV1*s1-refV2); + + VERIFY_IS_APPROX(v1*=s1, refV1*=s1); + VERIFY_IS_APPROX(v1/=s1, refV1/=s1); + + VERIFY_IS_APPROX(v1+=v2, refV1+=refV2); + VERIFY_IS_APPROX(v1-=v2, refV1-=refV2); + + VERIFY_IS_APPROX(v1.eigen2_dot(v2), refV1.eigen2_dot(refV2)); + VERIFY_IS_APPROX(v1.eigen2_dot(refV2), refV1.eigen2_dot(refV2)); + +} + +void test_eigen2_sparse_vector() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( sparse_vector<double>(8, 8) ); + CALL_SUBTEST_2( sparse_vector<std::complex<double> >(16, 16) ); + CALL_SUBTEST_1( sparse_vector<double>(299, 535) ); + } +} + diff --git a/eigen/test/eigen2/eigen2_stdvector.cpp b/eigen/test/eigen2/eigen2_stdvector.cpp new file mode 100644 index 0000000..6ab05c2 --- /dev/null +++ b/eigen/test/eigen2/eigen2_stdvector.cpp @@ -0,0 +1,148 @@ +// 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 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// 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 <Eigen/StdVector> +#include "main.h" +#include <Eigen/Geometry> + +template<typename MatrixType> +void check_stdvector_matrix(const MatrixType& m) +{ + int rows = m.rows(); + int cols = m.cols(); + MatrixType x = MatrixType::Random(rows,cols), y = MatrixType::Random(rows,cols); + std::vector<MatrixType, aligned_allocator<MatrixType> > v(10, MatrixType(rows,cols)), w(20, y); + v[5] = x; + w[6] = v[5]; + VERIFY_IS_APPROX(w[6], v[5]); + v = w; + for(int i = 0; i < 20; i++) + { + VERIFY_IS_APPROX(w[i], v[i]); + } + + v.resize(21); + v[20] = x; + VERIFY_IS_APPROX(v[20], x); + v.resize(22,y); + VERIFY_IS_APPROX(v[21], y); + v.push_back(x); + VERIFY_IS_APPROX(v[22], x); + VERIFY((std::size_t)&(v[22]) == (std::size_t)&(v[21]) + sizeof(MatrixType)); + + // do a lot of push_back such that the vector gets internally resized + // (with memory reallocation) + MatrixType* ref = &w[0]; + for(int i=0; i<30 || ((ref==&w[0]) && i<300); ++i) + v.push_back(w[i%w.size()]); + for(unsigned int i=23; i<v.size(); ++i) + { + VERIFY(v[i]==w[(i-23)%w.size()]); + } +} + +template<typename TransformType> +void check_stdvector_transform(const TransformType&) +{ + typedef typename TransformType::MatrixType MatrixType; + TransformType x(MatrixType::Random()), y(MatrixType::Random()); + std::vector<TransformType, aligned_allocator<TransformType> > v(10), w(20, y); + v[5] = x; + w[6] = v[5]; + VERIFY_IS_APPROX(w[6], v[5]); + v = w; + for(int i = 0; i < 20; i++) + { + VERIFY_IS_APPROX(w[i], v[i]); + } + + v.resize(21); + v[20] = x; + VERIFY_IS_APPROX(v[20], x); + v.resize(22,y); + VERIFY_IS_APPROX(v[21], y); + v.push_back(x); + VERIFY_IS_APPROX(v[22], x); + VERIFY((std::size_t)&(v[22]) == (std::size_t)&(v[21]) + sizeof(TransformType)); + + // do a lot of push_back such that the vector gets internally resized + // (with memory reallocation) + TransformType* ref = &w[0]; + for(int i=0; i<30 || ((ref==&w[0]) && i<300); ++i) + v.push_back(w[i%w.size()]); + for(unsigned int i=23; i<v.size(); ++i) + { + VERIFY(v[i].matrix()==w[(i-23)%w.size()].matrix()); + } +} + +template<typename QuaternionType> +void check_stdvector_quaternion(const QuaternionType&) +{ + typedef typename QuaternionType::Coefficients Coefficients; + QuaternionType x(Coefficients::Random()), y(Coefficients::Random()); + std::vector<QuaternionType, aligned_allocator<QuaternionType> > v(10), w(20, y); + v[5] = x; + w[6] = v[5]; + VERIFY_IS_APPROX(w[6], v[5]); + v = w; + for(int i = 0; i < 20; i++) + { + VERIFY_IS_APPROX(w[i], v[i]); + } + + v.resize(21); + v[20] = x; + VERIFY_IS_APPROX(v[20], x); + v.resize(22,y); + VERIFY_IS_APPROX(v[21], y); + v.push_back(x); + VERIFY_IS_APPROX(v[22], x); + VERIFY((std::size_t)&(v[22]) == (std::size_t)&(v[21]) + sizeof(QuaternionType)); + + // do a lot of push_back such that the vector gets internally resized + // (with memory reallocation) + QuaternionType* ref = &w[0]; + for(int i=0; i<30 || ((ref==&w[0]) && i<300); ++i) + v.push_back(w[i%w.size()]); + for(unsigned int i=23; i<v.size(); ++i) + { + VERIFY(v[i].coeffs()==w[(i-23)%w.size()].coeffs()); + } +} + +void test_eigen2_stdvector() +{ + // some non vectorizable fixed sizes + CALL_SUBTEST_1(check_stdvector_matrix(Vector2f())); + CALL_SUBTEST_1(check_stdvector_matrix(Matrix3f())); + CALL_SUBTEST_1(check_stdvector_matrix(Matrix3d())); + + // some vectorizable fixed sizes + CALL_SUBTEST_2(check_stdvector_matrix(Matrix2f())); + CALL_SUBTEST_2(check_stdvector_matrix(Vector4f())); + CALL_SUBTEST_2(check_stdvector_matrix(Matrix4f())); + CALL_SUBTEST_2(check_stdvector_matrix(Matrix4d())); + + // some dynamic sizes + CALL_SUBTEST_3(check_stdvector_matrix(MatrixXd(1,1))); + CALL_SUBTEST_3(check_stdvector_matrix(VectorXd(20))); + CALL_SUBTEST_3(check_stdvector_matrix(RowVectorXf(20))); + CALL_SUBTEST_3(check_stdvector_matrix(MatrixXcf(10,10))); + + // some Transform + CALL_SUBTEST_4(check_stdvector_transform(Transform2f())); + CALL_SUBTEST_4(check_stdvector_transform(Transform3f())); + CALL_SUBTEST_4(check_stdvector_transform(Transform3d())); + //CALL_SUBTEST_4(check_stdvector_transform(Transform4d())); + + // some Quaternion + CALL_SUBTEST_5(check_stdvector_quaternion(Quaternionf())); + CALL_SUBTEST_5(check_stdvector_quaternion(Quaternionf())); +} diff --git a/eigen/test/eigen2/eigen2_submatrices.cpp b/eigen/test/eigen2/eigen2_submatrices.cpp new file mode 100644 index 0000000..dee970b --- /dev/null +++ b/eigen/test/eigen2/eigen2_submatrices.cpp @@ -0,0 +1,142 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// 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" + +// check minor separately in order to avoid the possible creation of a zero-sized +// array. Comes from a compilation error with gcc-3.4 or gcc-4 with -ansi -pedantic. +// Another solution would be to declare the array like this: T m_data[Size==0?1:Size]; in ei_matrix_storage +// but this is probably not bad to raise such an error at compile time... +template<typename Scalar, int _Rows, int _Cols> struct CheckMinor +{ + typedef Matrix<Scalar, _Rows, _Cols> MatrixType; + CheckMinor(MatrixType& m1, int r1, int c1) + { + int rows = m1.rows(); + int cols = m1.cols(); + + Matrix<Scalar, Dynamic, Dynamic> mi = m1.minor(0,0).eval(); + VERIFY_IS_APPROX(mi, m1.block(1,1,rows-1,cols-1)); + mi = m1.minor(r1,c1); + VERIFY_IS_APPROX(mi.transpose(), m1.transpose().minor(c1,r1)); + //check operator(), both constant and non-constant, on minor() + m1.minor(r1,c1)(0,0) = m1.minor(0,0)(0,0); + } +}; + +template<typename Scalar> struct CheckMinor<Scalar,1,1> +{ + typedef Matrix<Scalar, 1, 1> MatrixType; + CheckMinor(MatrixType&, int, int) {} +}; + +template<typename MatrixType> void submatrices(const MatrixType& m) +{ + /* this test covers the following files: + Row.h Column.h Block.h Minor.h DiagonalCoeffs.h + */ + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; + typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType; + int rows = m.rows(); + int cols = m.cols(); + + MatrixType m1 = MatrixType::Random(rows, cols), + m2 = MatrixType::Random(rows, cols), + m3(rows, cols), + ones = MatrixType::Ones(rows, cols), + square = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> + ::Random(rows, rows); + VectorType v1 = VectorType::Random(rows); + + Scalar s1 = ei_random<Scalar>(); + + 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); + + //check row() and col() + VERIFY_IS_APPROX(m1.col(c1).transpose(), m1.transpose().row(c1)); + VERIFY_IS_APPROX(square.row(r1).eigen2_dot(m1.col(c1)), (square.lazy() * m1.conjugate())(r1,c1)); + //check operator(), both constant and non-constant, on row() and col() + m1.row(r1) += s1 * m1.row(r2); + m1.col(c1) += s1 * m1.col(c2); + + //check block() + Matrix<Scalar,Dynamic,Dynamic> b1(1,1); b1(0,0) = m1(r1,c1); + RowVectorType br1(m1.block(r1,0,1,cols)); + VectorType bc1(m1.block(0,c1,rows,1)); + VERIFY_IS_APPROX(b1, m1.block(r1,c1,1,1)); + VERIFY_IS_APPROX(m1.row(r1), br1); + VERIFY_IS_APPROX(m1.col(c1), bc1); + //check operator(), both constant and non-constant, on block() + m1.block(r1,c1,r2-r1+1,c2-c1+1) = s1 * m2.block(0, 0, r2-r1+1,c2-c1+1); + m1.block(r1,c1,r2-r1+1,c2-c1+1)(r2-r1,c2-c1) = m2.block(0, 0, r2-r1+1,c2-c1+1)(0,0); + + //check minor() + CheckMinor<Scalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime> checkminor(m1,r1,c1); + + //check diagonal() + VERIFY_IS_APPROX(m1.diagonal(), m1.transpose().diagonal()); + m2.diagonal() = 2 * m1.diagonal(); + m2.diagonal()[0] *= 3; + VERIFY_IS_APPROX(m2.diagonal()[0], static_cast<Scalar>(6) * m1.diagonal()[0]); + + enum { + BlockRows = EIGEN_SIZE_MIN_PREFER_FIXED(MatrixType::RowsAtCompileTime,2), + BlockCols = EIGEN_SIZE_MIN_PREFER_FIXED(MatrixType::ColsAtCompileTime,5) + }; + if (rows>=5 && cols>=8) + { + // test fixed block() as lvalue + m1.template block<BlockRows,BlockCols>(1,1) *= s1; + // test operator() on fixed block() both as constant and non-constant + m1.template block<BlockRows,BlockCols>(1,1)(0, 3) = m1.template block<2,5>(1,1)(1,2); + // check that fixed block() and block() agree + Matrix<Scalar,Dynamic,Dynamic> b = m1.template block<BlockRows,BlockCols>(3,3); + VERIFY_IS_APPROX(b, m1.block(3,3,BlockRows,BlockCols)); + } + + if (rows>2) + { + // test sub vectors + VERIFY_IS_APPROX(v1.template start<2>(), v1.block(0,0,2,1)); + VERIFY_IS_APPROX(v1.template start<2>(), v1.start(2)); + VERIFY_IS_APPROX(v1.template start<2>(), v1.segment(0,2)); + VERIFY_IS_APPROX(v1.template start<2>(), v1.template segment<2>(0)); + int i = rows-2; + VERIFY_IS_APPROX(v1.template end<2>(), v1.block(i,0,2,1)); + VERIFY_IS_APPROX(v1.template end<2>(), v1.end(2)); + VERIFY_IS_APPROX(v1.template end<2>(), v1.segment(i,2)); + VERIFY_IS_APPROX(v1.template end<2>(), v1.template segment<2>(i)); + i = ei_random(0,rows-2); + VERIFY_IS_APPROX(v1.segment(i,2), v1.template segment<2>(i)); + } + + // stress some basic stuffs with block matrices + VERIFY(ei_real(ones.col(c1).sum()) == RealScalar(rows)); + VERIFY(ei_real(ones.row(r1).sum()) == RealScalar(cols)); + + VERIFY(ei_real(ones.col(c1).eigen2_dot(ones.col(c2))) == RealScalar(rows)); + VERIFY(ei_real(ones.row(r1).eigen2_dot(ones.row(r2))) == RealScalar(cols)); +} + +void test_eigen2_submatrices() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( submatrices(Matrix<float, 1, 1>()) ); + CALL_SUBTEST_2( submatrices(Matrix4d()) ); + CALL_SUBTEST_3( submatrices(MatrixXcf(3, 3)) ); + CALL_SUBTEST_4( submatrices(MatrixXi(8, 12)) ); + CALL_SUBTEST_5( submatrices(MatrixXcd(20, 20)) ); + CALL_SUBTEST_6( submatrices(MatrixXf(20, 20)) ); + } +} diff --git a/eigen/test/eigen2/eigen2_sum.cpp b/eigen/test/eigen2/eigen2_sum.cpp new file mode 100644 index 0000000..b47057c --- /dev/null +++ b/eigen/test/eigen2/eigen2_sum.cpp @@ -0,0 +1,71 @@ +// 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 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" + +template<typename MatrixType> void matrixSum(const MatrixType& m) +{ + typedef typename MatrixType::Scalar Scalar; + + int rows = m.rows(); + int cols = m.cols(); + + MatrixType m1 = MatrixType::Random(rows, cols); + + 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 x = Scalar(0); + for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) x += m1(i,j); + VERIFY_IS_APPROX(m1.sum(), x); +} + +template<typename VectorType> void vectorSum(const VectorType& w) +{ + typedef typename VectorType::Scalar Scalar; + int size = w.size(); + + VectorType v = VectorType::Random(size); + for(int i = 1; i < size; i++) + { + Scalar s = Scalar(0); + for(int j = 0; j < i; j++) s += v[j]; + VERIFY_IS_APPROX(s, v.start(i).sum()); + } + + for(int i = 0; i < size-1; i++) + { + Scalar s = Scalar(0); + for(int j = i; j < size; j++) s += v[j]; + VERIFY_IS_APPROX(s, v.end(size-i).sum()); + } + + for(int i = 0; i < size/2; i++) + { + Scalar s = Scalar(0); + for(int j = i; j < size-i; j++) s += v[j]; + VERIFY_IS_APPROX(s, v.segment(i, size-2*i).sum()); + } +} + +void test_eigen2_sum() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( matrixSum(Matrix<float, 1, 1>()) ); + CALL_SUBTEST_2( matrixSum(Matrix2f()) ); + CALL_SUBTEST_3( matrixSum(Matrix4d()) ); + CALL_SUBTEST_4( matrixSum(MatrixXcf(3, 3)) ); + CALL_SUBTEST_5( matrixSum(MatrixXf(8, 12)) ); + CALL_SUBTEST_6( matrixSum(MatrixXi(8, 12)) ); + } + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_5( vectorSum(VectorXf(5)) ); + CALL_SUBTEST_7( vectorSum(VectorXd(10)) ); + CALL_SUBTEST_5( vectorSum(VectorXf(33)) ); + } +} diff --git a/eigen/test/eigen2/eigen2_svd.cpp b/eigen/test/eigen2/eigen2_svd.cpp new file mode 100644 index 0000000..d4689a5 --- /dev/null +++ b/eigen/test/eigen2/eigen2_svd.cpp @@ -0,0 +1,87 @@ +// 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 Gael Guennebaud <g.gael@free.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" +#include <Eigen/SVD> + +template<typename MatrixType> void svd(const MatrixType& m) +{ + /* this test covers the following files: + SVD.h + */ + int rows = m.rows(); + int cols = m.cols(); + + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + MatrixType a = MatrixType::Random(rows,cols); + Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> b = + 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(); + matU.block(0,0,rows,cols) = svd.matrixU(); + VERIFY_IS_APPROX(a, matU * sigma * svd.matrixV().transpose()); + } + + + if (rows==cols) + { + if (ei_is_same_type<RealScalar,float>::ret) + { + MatrixType a1 = MatrixType::Random(rows,cols); + a += a * a.adjoint() + a1 * a1.adjoint(); + } + SVD<MatrixType> svd(a); + svd.solve(b, &x); + VERIFY_IS_APPROX(a * x,b); + } + + + if(rows==cols) + { + SVD<MatrixType> svd(a); + MatrixType unitary, positive; + svd.computeUnitaryPositive(&unitary, &positive); + VERIFY_IS_APPROX(unitary * unitary.adjoint(), MatrixType::Identity(unitary.rows(),unitary.rows())); + VERIFY_IS_APPROX(positive, positive.adjoint()); + for(int i = 0; i < rows; i++) VERIFY(positive.diagonal()[i] >= 0); // cheap necessary (not sufficient) condition for positivity + VERIFY_IS_APPROX(unitary*positive, a); + + svd.computePositiveUnitary(&positive, &unitary); + VERIFY_IS_APPROX(unitary * unitary.adjoint(), MatrixType::Identity(unitary.rows(),unitary.rows())); + VERIFY_IS_APPROX(positive, positive.adjoint()); + for(int i = 0; i < rows; i++) VERIFY(positive.diagonal()[i] >= 0); // cheap necessary (not sufficient) condition for positivity + VERIFY_IS_APPROX(positive*unitary, a); + } +} + +void test_eigen2_svd() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( svd(Matrix3f()) ); + CALL_SUBTEST_2( svd(Matrix4d()) ); + CALL_SUBTEST_3( svd(MatrixXf(7,7)) ); + CALL_SUBTEST_4( svd(MatrixXd(14,7)) ); + // complex are not implemented yet +// CALL_SUBTEST( svd(MatrixXcd(6,6)) ); +// CALL_SUBTEST( svd(MatrixXcf(3,3)) ); + SVD<MatrixXf> s; + MatrixXf m = MatrixXf::Random(10,1); + s.compute(m); + } +} diff --git a/eigen/test/eigen2/eigen2_swap.cpp b/eigen/test/eigen2/eigen2_swap.cpp new file mode 100644 index 0000000..f3a8846 --- /dev/null +++ b/eigen/test/eigen2/eigen2_swap.cpp @@ -0,0 +1,83 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// 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/. + +#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_eigen2_swap() +{ + CALL_SUBTEST_1( swap(Matrix3f()) ); // fixed size, no vectorization + CALL_SUBTEST_1( swap(Matrix4d()) ); // fixed size, possible vectorization + CALL_SUBTEST_1( swap(MatrixXd(3,3)) ); // dyn size, no vectorization + CALL_SUBTEST_1( swap(MatrixXf(30,30)) ); // dyn size, possible vectorization +} diff --git a/eigen/test/eigen2/eigen2_triangular.cpp b/eigen/test/eigen2/eigen2_triangular.cpp new file mode 100644 index 0000000..6f17b7d --- /dev/null +++ b/eigen/test/eigen2/eigen2_triangular.cpp @@ -0,0 +1,148 @@ +// 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 Gael Guennebaud <gael.guennebaud@gmail.com> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" + +template<typename MatrixType> void triangular(const MatrixType& m) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + + RealScalar largerEps = 10*test_precision<RealScalar>(); + + int rows = m.rows(); + int cols = m.cols(); + + MatrixType m1 = MatrixType::Random(rows, cols), + m2 = MatrixType::Random(rows, cols), + m3(rows, cols), + m4(rows, cols), + r1(rows, cols), + r2(rows, cols); + + MatrixType m1up = m1.template part<Eigen::UpperTriangular>(); + MatrixType m2up = m2.template part<Eigen::UpperTriangular>(); + + if (rows*cols>1) + { + VERIFY(m1up.isUpperTriangular()); + VERIFY(m2up.transpose().isLowerTriangular()); + VERIFY(!m2.isLowerTriangular()); + } + +// VERIFY_IS_APPROX(m1up.transpose() * m2, m1.upper().transpose().lower() * m2); + + // test overloaded operator+= + r1.setZero(); + r2.setZero(); + r1.template part<Eigen::UpperTriangular>() += m1; + r2 += m1up; + VERIFY_IS_APPROX(r1,r2); + + // test overloaded operator= + m1.setZero(); + m1.template part<Eigen::UpperTriangular>() = (m2.transpose() * m2).lazy(); + m3 = m2.transpose() * m2; + VERIFY_IS_APPROX(m3.template part<Eigen::LowerTriangular>().transpose(), m1); + + // test overloaded operator= + m1.setZero(); + m1.template part<Eigen::LowerTriangular>() = (m2.transpose() * m2).lazy(); + VERIFY_IS_APPROX(m3.template part<Eigen::LowerTriangular>(), m1); + + VERIFY_IS_APPROX(m3.template part<Diagonal>(), m3.diagonal().asDiagonal()); + + m1 = MatrixType::Random(rows, cols); + for (int i=0; i<rows; ++i) + while (ei_abs2(m1(i,i))<1e-3) m1(i,i) = ei_random<Scalar>(); + + Transpose<MatrixType> trm4(m4); + // test back and forward subsitution + m3 = m1.template part<Eigen::LowerTriangular>(); + VERIFY(m3.template marked<Eigen::LowerTriangular>().solveTriangular(m3).cwise().abs().isIdentity(test_precision<RealScalar>())); + VERIFY(m3.transpose().template marked<Eigen::UpperTriangular>() + .solveTriangular(m3.transpose()).cwise().abs().isIdentity(test_precision<RealScalar>())); + // check M * inv(L) using in place API + m4 = m3; + m3.transpose().template marked<Eigen::UpperTriangular>().solveTriangularInPlace(trm4); + VERIFY(m4.cwise().abs().isIdentity(test_precision<RealScalar>())); + + m3 = m1.template part<Eigen::UpperTriangular>(); + VERIFY(m3.template marked<Eigen::UpperTriangular>().solveTriangular(m3).cwise().abs().isIdentity(test_precision<RealScalar>())); + VERIFY(m3.transpose().template marked<Eigen::LowerTriangular>() + .solveTriangular(m3.transpose()).cwise().abs().isIdentity(test_precision<RealScalar>())); + // check M * inv(U) using in place API + m4 = m3; + m3.transpose().template marked<Eigen::LowerTriangular>().solveTriangularInPlace(trm4); + VERIFY(m4.cwise().abs().isIdentity(test_precision<RealScalar>())); + + m3 = m1.template part<Eigen::UpperTriangular>(); + VERIFY(m2.isApprox(m3 * (m3.template marked<Eigen::UpperTriangular>().solveTriangular(m2)), largerEps)); + m3 = m1.template part<Eigen::LowerTriangular>(); + VERIFY(m2.isApprox(m3 * (m3.template marked<Eigen::LowerTriangular>().solveTriangular(m2)), largerEps)); + + VERIFY((m1.template part<Eigen::UpperTriangular>() * m2.template part<Eigen::UpperTriangular>()).isUpperTriangular()); + + // test swap + m1.setOnes(); + m2.setZero(); + m2.template part<Eigen::UpperTriangular>().swap(m1); + m3.setZero(); + m3.template part<Eigen::UpperTriangular>().setOnes(); + VERIFY_IS_APPROX(m2,m3); + +} + +void selfadjoint() +{ + Matrix2i m; + m << 1, 2, + 3, 4; + + Matrix2i m1 = Matrix2i::Zero(); + m1.part<SelfAdjoint>() = m; + Matrix2i ref1; + ref1 << 1, 2, + 2, 4; + VERIFY(m1 == ref1); + + Matrix2i m2 = Matrix2i::Zero(); + m2.part<SelfAdjoint>() = m.part<UpperTriangular>(); + Matrix2i ref2; + ref2 << 1, 2, + 2, 4; + VERIFY(m2 == ref2); + + Matrix2i m3 = Matrix2i::Zero(); + m3.part<SelfAdjoint>() = m.part<LowerTriangular>(); + Matrix2i ref3; + ref3 << 1, 0, + 0, 4; + VERIFY(m3 == ref3); + + // example inspired from bug 159 + int array[] = {1, 2, 3, 4}; + Matrix2i::Map(array).part<SelfAdjoint>() = Matrix2i::Random().part<LowerTriangular>(); + + std::cout << "hello\n" << array << std::endl; +} + +void test_eigen2_triangular() +{ + CALL_SUBTEST_8( selfadjoint() ); + for(int i = 0; i < g_repeat ; i++) { + CALL_SUBTEST_1( triangular(Matrix<float, 1, 1>()) ); + CALL_SUBTEST_2( triangular(Matrix<float, 2, 2>()) ); + CALL_SUBTEST_3( triangular(Matrix3d()) ); + CALL_SUBTEST_4( triangular(MatrixXcf(4, 4)) ); + CALL_SUBTEST_5( triangular(Matrix<std::complex<float>,8, 8>()) ); + CALL_SUBTEST_6( triangular(MatrixXd(17,17)) ); + CALL_SUBTEST_7( triangular(Matrix<float,Dynamic,Dynamic,RowMajor>(5, 5)) ); + } +} diff --git a/eigen/test/eigen2/eigen2_unalignedassert.cpp b/eigen/test/eigen2/eigen2_unalignedassert.cpp new file mode 100644 index 0000000..d10b6f5 --- /dev/null +++ b/eigen/test/eigen2/eigen2_unalignedassert.cpp @@ -0,0 +1,116 @@ +// 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 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// 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" + +struct Good1 +{ + MatrixXd m; // good: m will allocate its own array, taking care of alignment. + Good1() : m(20,20) {} +}; + +struct Good2 +{ + Matrix3d m; // good: m's size isn't a multiple of 16 bytes, so m doesn't have to be aligned +}; + +struct Good3 +{ + Vector2f m; // good: same reason +}; + +struct Bad4 +{ + Vector2d m; // bad: sizeof(m)%16==0 so alignment is required +}; + +struct Bad5 +{ + Matrix<float, 2, 6> m; // bad: same reason +}; + +struct Bad6 +{ + Matrix<double, 3, 4> m; // bad: same reason +}; + +struct Good7 +{ + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + Vector2d m; + float f; // make the struct have sizeof%16!=0 to make it a little more tricky when we allow an array of 2 such objects +}; + +struct Good8 +{ + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + float f; // try the f at first -- the EIGEN_ALIGN_128 attribute of m should make that still work + Matrix4f m; +}; + +struct Good9 +{ + Matrix<float,2,2,DontAlign> m; // good: no alignment requested + float f; +}; + +template<bool Align> struct Depends +{ + EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(Align) + Vector2d m; + float f; +}; + +template<typename T> +void check_unalignedassert_good() +{ + T *x, *y; + x = new T; + delete x; + y = new T[2]; + delete[] y; +} + +#if EIGEN_ARCH_WANTS_ALIGNMENT +template<typename T> +void check_unalignedassert_bad() +{ + float buf[sizeof(T)+16]; + float *unaligned = buf; + while((reinterpret_cast<std::size_t>(unaligned)&0xf)==0) ++unaligned; // make sure unaligned is really unaligned + T *x = ::new(static_cast<void*>(unaligned)) T; + x->~T(); +} +#endif + +void unalignedassert() +{ + check_unalignedassert_good<Good1>(); + check_unalignedassert_good<Good2>(); + check_unalignedassert_good<Good3>(); +#if EIGEN_ARCH_WANTS_ALIGNMENT + VERIFY_RAISES_ASSERT(check_unalignedassert_bad<Bad4>()); + VERIFY_RAISES_ASSERT(check_unalignedassert_bad<Bad5>()); + VERIFY_RAISES_ASSERT(check_unalignedassert_bad<Bad6>()); +#endif + + check_unalignedassert_good<Good7>(); + check_unalignedassert_good<Good8>(); + check_unalignedassert_good<Good9>(); + check_unalignedassert_good<Depends<true> >(); + +#if EIGEN_ARCH_WANTS_ALIGNMENT + VERIFY_RAISES_ASSERT(check_unalignedassert_bad<Depends<false> >()); +#endif +} + +void test_eigen2_unalignedassert() +{ + CALL_SUBTEST(unalignedassert()); +} diff --git a/eigen/test/eigen2/eigen2_visitor.cpp b/eigen/test/eigen2/eigen2_visitor.cpp new file mode 100644 index 0000000..4781991 --- /dev/null +++ b/eigen/test/eigen2/eigen2_visitor.cpp @@ -0,0 +1,116 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" + +template<typename MatrixType> void matrixVisitor(const MatrixType& p) +{ + typedef typename MatrixType::Scalar Scalar; + + int rows = p.rows(); + int cols = p.cols(); + + // construct a random matrix where all coefficients are different + MatrixType m; + m = MatrixType::Random(rows, cols); + for(int i = 0; i < m.size(); i++) + for(int i2 = 0; i2 < i; i2++) + while(m(i) == m(i2)) // yes, == + m(i) = ei_random<Scalar>(); + + Scalar minc = Scalar(1000), maxc = Scalar(-1000); + int minrow=0,mincol=0,maxrow=0,maxcol=0; + for(int j = 0; j < cols; j++) + for(int i = 0; i < rows; i++) + { + if(m(i,j) < minc) + { + minc = m(i,j); + minrow = i; + mincol = j; + } + if(m(i,j) > maxc) + { + maxc = m(i,j); + maxrow = i; + maxcol = j; + } + } + int eigen_minrow, eigen_mincol, eigen_maxrow, eigen_maxcol; + Scalar eigen_minc, eigen_maxc; + eigen_minc = m.minCoeff(&eigen_minrow,&eigen_mincol); + eigen_maxc = m.maxCoeff(&eigen_maxrow,&eigen_maxcol); + VERIFY(minrow == eigen_minrow); + VERIFY(maxrow == eigen_maxrow); + VERIFY(mincol == eigen_mincol); + VERIFY(maxcol == eigen_maxcol); + VERIFY_IS_APPROX(minc, eigen_minc); + VERIFY_IS_APPROX(maxc, eigen_maxc); + VERIFY_IS_APPROX(minc, m.minCoeff()); + VERIFY_IS_APPROX(maxc, m.maxCoeff()); +} + +template<typename VectorType> void vectorVisitor(const VectorType& w) +{ + typedef typename VectorType::Scalar Scalar; + + int size = w.size(); + + // construct a random vector where all coefficients are different + VectorType v; + v = VectorType::Random(size); + for(int i = 0; i < size; i++) + for(int i2 = 0; i2 < i; i2++) + while(v(i) == v(i2)) // yes, == + v(i) = ei_random<Scalar>(); + + Scalar minc = Scalar(1000), maxc = Scalar(-1000); + int minidx=0,maxidx=0; + for(int i = 0; i < size; i++) + { + if(v(i) < minc) + { + minc = v(i); + minidx = i; + } + if(v(i) > maxc) + { + maxc = v(i); + maxidx = i; + } + } + int eigen_minidx, eigen_maxidx; + Scalar eigen_minc, eigen_maxc; + eigen_minc = v.minCoeff(&eigen_minidx); + eigen_maxc = v.maxCoeff(&eigen_maxidx); + VERIFY(minidx == eigen_minidx); + VERIFY(maxidx == eigen_maxidx); + VERIFY_IS_APPROX(minc, eigen_minc); + VERIFY_IS_APPROX(maxc, eigen_maxc); + VERIFY_IS_APPROX(minc, v.minCoeff()); + VERIFY_IS_APPROX(maxc, v.maxCoeff()); +} + +void test_eigen2_visitor() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1( matrixVisitor(Matrix<float, 1, 1>()) ); + CALL_SUBTEST_2( matrixVisitor(Matrix2f()) ); + CALL_SUBTEST_3( matrixVisitor(Matrix4d()) ); + CALL_SUBTEST_4( matrixVisitor(MatrixXd(8, 12)) ); + CALL_SUBTEST_5( matrixVisitor(Matrix<double,Dynamic,Dynamic,RowMajor>(20, 20)) ); + CALL_SUBTEST_6( matrixVisitor(MatrixXi(8, 12)) ); + } + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_7( vectorVisitor(Vector4f()) ); + CALL_SUBTEST_4( vectorVisitor(VectorXd(10)) ); + CALL_SUBTEST_4( vectorVisitor(RowVectorXd(10)) ); + CALL_SUBTEST_8( vectorVisitor(VectorXf(33)) ); + } +} diff --git a/eigen/test/eigen2/gsl_helper.h b/eigen/test/eigen2/gsl_helper.h new file mode 100644 index 0000000..d1d8545 --- /dev/null +++ b/eigen/test/eigen2/gsl_helper.h @@ -0,0 +1,175 @@ +// 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 Gael Guennebaud <g.gael@free.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_GSL_HELPER +#define EIGEN_GSL_HELPER + +#include <Eigen/Core> + +#include <gsl/gsl_blas.h> +#include <gsl/gsl_multifit.h> +#include <gsl/gsl_eigen.h> +#include <gsl/gsl_linalg.h> +#include <gsl/gsl_complex.h> +#include <gsl/gsl_complex_math.h> + +namespace Eigen { + +template<typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex> struct GslTraits +{ + typedef gsl_matrix* Matrix; + typedef gsl_vector* Vector; + static Matrix createMatrix(int rows, int cols) { return gsl_matrix_alloc(rows,cols); } + static Vector createVector(int size) { return gsl_vector_alloc(size); } + static void free(Matrix& m) { gsl_matrix_free(m); m=0; } + static void free(Vector& m) { gsl_vector_free(m); m=0; } + static void prod(const Matrix& m, const Vector& v, Vector& x) { gsl_blas_dgemv(CblasNoTrans,1,m,v,0,x); } + static void cholesky(Matrix& m) { gsl_linalg_cholesky_decomp(m); } + static void cholesky_solve(const Matrix& m, const Vector& b, Vector& x) { gsl_linalg_cholesky_solve(m,b,x); } + static void eigen_symm(const Matrix& m, Vector& eval, Matrix& evec) + { + gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc(m->size1); + Matrix a = createMatrix(m->size1, m->size2); + gsl_matrix_memcpy(a, m); + gsl_eigen_symmv(a,eval,evec,w); + gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_VAL_ASC); + gsl_eigen_symmv_free(w); + free(a); + } + static void eigen_symm_gen(const Matrix& m, const Matrix& _b, Vector& eval, Matrix& evec) + { + gsl_eigen_gensymmv_workspace * w = gsl_eigen_gensymmv_alloc(m->size1); + Matrix a = createMatrix(m->size1, m->size2); + Matrix b = createMatrix(_b->size1, _b->size2); + gsl_matrix_memcpy(a, m); + gsl_matrix_memcpy(b, _b); + gsl_eigen_gensymmv(a,b,eval,evec,w); + gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_VAL_ASC); + gsl_eigen_gensymmv_free(w); + free(a); + } +}; + +template<typename Scalar> struct GslTraits<Scalar,true> +{ + typedef gsl_matrix_complex* Matrix; + typedef gsl_vector_complex* Vector; + static Matrix createMatrix(int rows, int cols) { return gsl_matrix_complex_alloc(rows,cols); } + static Vector createVector(int size) { return gsl_vector_complex_alloc(size); } + static void free(Matrix& m) { gsl_matrix_complex_free(m); m=0; } + static void free(Vector& m) { gsl_vector_complex_free(m); m=0; } + static void cholesky(Matrix& m) { gsl_linalg_complex_cholesky_decomp(m); } + static void cholesky_solve(const Matrix& m, const Vector& b, Vector& x) { gsl_linalg_complex_cholesky_solve(m,b,x); } + static void prod(const Matrix& m, const Vector& v, Vector& x) + { gsl_blas_zgemv(CblasNoTrans,gsl_complex_rect(1,0),m,v,gsl_complex_rect(0,0),x); } + static void eigen_symm(const Matrix& m, gsl_vector* &eval, Matrix& evec) + { + gsl_eigen_hermv_workspace * w = gsl_eigen_hermv_alloc(m->size1); + Matrix a = createMatrix(m->size1, m->size2); + gsl_matrix_complex_memcpy(a, m); + gsl_eigen_hermv(a,eval,evec,w); + gsl_eigen_hermv_sort(eval, evec, GSL_EIGEN_SORT_VAL_ASC); + gsl_eigen_hermv_free(w); + free(a); + } + static void eigen_symm_gen(const Matrix& m, const Matrix& _b, gsl_vector* &eval, Matrix& evec) + { + gsl_eigen_genhermv_workspace * w = gsl_eigen_genhermv_alloc(m->size1); + Matrix a = createMatrix(m->size1, m->size2); + Matrix b = createMatrix(_b->size1, _b->size2); + gsl_matrix_complex_memcpy(a, m); + gsl_matrix_complex_memcpy(b, _b); + gsl_eigen_genhermv(a,b,eval,evec,w); + gsl_eigen_hermv_sort(eval, evec, GSL_EIGEN_SORT_VAL_ASC); + gsl_eigen_genhermv_free(w); + free(a); + } +}; + +template<typename MatrixType> +void convert(const MatrixType& m, gsl_matrix* &res) +{ +// if (res) +// gsl_matrix_free(res); + res = gsl_matrix_alloc(m.rows(), m.cols()); + for (int i=0 ; i<m.rows() ; ++i) + for (int j=0 ; j<m.cols(); ++j) + gsl_matrix_set(res, i, j, m(i,j)); +} + +template<typename MatrixType> +void convert(const gsl_matrix* m, MatrixType& res) +{ + res.resize(int(m->size1), int(m->size2)); + for (int i=0 ; i<res.rows() ; ++i) + for (int j=0 ; j<res.cols(); ++j) + res(i,j) = gsl_matrix_get(m,i,j); +} + +template<typename VectorType> +void convert(const VectorType& m, gsl_vector* &res) +{ + if (res) gsl_vector_free(res); + res = gsl_vector_alloc(m.size()); + for (int i=0 ; i<m.size() ; ++i) + gsl_vector_set(res, i, m[i]); +} + +template<typename VectorType> +void convert(const gsl_vector* m, VectorType& res) +{ + res.resize (m->size); + for (int i=0 ; i<res.rows() ; ++i) + res[i] = gsl_vector_get(m, i); +} + +template<typename MatrixType> +void convert(const MatrixType& m, gsl_matrix_complex* &res) +{ + res = gsl_matrix_complex_alloc(m.rows(), m.cols()); + for (int i=0 ; i<m.rows() ; ++i) + for (int j=0 ; j<m.cols(); ++j) + { + gsl_matrix_complex_set(res, i, j, + gsl_complex_rect(m(i,j).real(), m(i,j).imag())); + } +} + +template<typename MatrixType> +void convert(const gsl_matrix_complex* m, MatrixType& res) +{ + res.resize(int(m->size1), int(m->size2)); + for (int i=0 ; i<res.rows() ; ++i) + for (int j=0 ; j<res.cols(); ++j) + res(i,j) = typename MatrixType::Scalar( + GSL_REAL(gsl_matrix_complex_get(m,i,j)), + GSL_IMAG(gsl_matrix_complex_get(m,i,j))); +} + +template<typename VectorType> +void convert(const VectorType& m, gsl_vector_complex* &res) +{ + res = gsl_vector_complex_alloc(m.size()); + for (int i=0 ; i<m.size() ; ++i) + gsl_vector_complex_set(res, i, gsl_complex_rect(m[i].real(), m[i].imag())); +} + +template<typename VectorType> +void convert(const gsl_vector_complex* m, VectorType& res) +{ + res.resize(m->size); + for (int i=0 ; i<res.rows() ; ++i) + res[i] = typename VectorType::Scalar( + GSL_REAL(gsl_vector_complex_get(m, i)), + GSL_IMAG(gsl_vector_complex_get(m, i))); +} + +} + +#endif // EIGEN_GSL_HELPER diff --git a/eigen/test/eigen2/main.h b/eigen/test/eigen2/main.h new file mode 100644 index 0000000..ad2ba19 --- /dev/null +++ b/eigen/test/eigen2/main.h @@ -0,0 +1,399 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include <cstdlib> +#include <ctime> +#include <iostream> +#include <string> +#include <vector> + +#ifndef EIGEN_TEST_FUNC +#error EIGEN_TEST_FUNC must be defined +#endif + +#define DEFAULT_REPEAT 10 + +namespace Eigen +{ + static std::vector<std::string> g_test_stack; + static int g_repeat; +} + +#define EI_PP_MAKE_STRING2(S) #S +#define EI_PP_MAKE_STRING(S) EI_PP_MAKE_STRING2(S) + +#define EI_PP_CAT2(a,b) a ## b +#define EI_PP_CAT(a,b) EI_PP_CAT2(a,b) + +#ifndef EIGEN_NO_ASSERTION_CHECKING + + namespace Eigen + { + static const bool should_raise_an_assert = false; + + // Used to avoid to raise two exceptions at a time in which + // case the exception is not properly caught. + // This may happen when a second exceptions is raise in a destructor. + static bool no_more_assert = false; + + struct eigen_assert_exception + { + eigen_assert_exception(void) {} + ~eigen_assert_exception() { Eigen::no_more_assert = false; } + }; + } + + // If EIGEN_DEBUG_ASSERTS is defined and if no assertion is raised while + // one should have been, then the list of excecuted assertions is printed out. + // + // EIGEN_DEBUG_ASSERTS is not enabled by default as it + // significantly increases the compilation time + // and might even introduce side effects that would hide + // some memory errors. + #ifdef EIGEN_DEBUG_ASSERTS + + namespace Eigen + { + static bool ei_push_assert = false; + static std::vector<std::string> eigen_assert_list; + } + + #define eigen_assert(a) \ + if( (!(a)) && (!no_more_assert) ) \ + { \ + Eigen::no_more_assert = true; \ + throw Eigen::eigen_assert_exception(); \ + } \ + else if (Eigen::ei_push_assert) \ + { \ + eigen_assert_list.push_back(std::string(EI_PP_MAKE_STRING(__FILE__)" ("EI_PP_MAKE_STRING(__LINE__)") : "#a) ); \ + } + + #define VERIFY_RAISES_ASSERT(a) \ + { \ + Eigen::no_more_assert = false; \ + try { \ + Eigen::eigen_assert_list.clear(); \ + Eigen::ei_push_assert = true; \ + a; \ + Eigen::ei_push_assert = false; \ + std::cerr << "One of the following asserts should have been raised:\n"; \ + for (uint ai=0 ; ai<eigen_assert_list.size() ; ++ai) \ + std::cerr << " " << eigen_assert_list[ai] << "\n"; \ + VERIFY(Eigen::should_raise_an_assert && # a); \ + } catch (Eigen::eigen_assert_exception e) { \ + Eigen::ei_push_assert = false; VERIFY(true); \ + } \ + } + + #else // EIGEN_DEBUG_ASSERTS + + #undef eigen_assert + + // see bug 89. The copy_bool here is working around a bug in gcc <= 4.3 + #define eigen_assert(a) \ + if( (!Eigen::internal::copy_bool(a)) && (!no_more_assert) ) \ + { \ + Eigen::no_more_assert = true; \ + throw Eigen::eigen_assert_exception(); \ + } + + #define VERIFY_RAISES_ASSERT(a) { \ + Eigen::no_more_assert = false; \ + try { a; VERIFY(Eigen::should_raise_an_assert && # a); } \ + catch (Eigen::eigen_assert_exception e) { VERIFY(true); } \ + } + + #endif // EIGEN_DEBUG_ASSERTS + + #define EIGEN_USE_CUSTOM_ASSERT + +#else // EIGEN_NO_ASSERTION_CHECKING + + #define VERIFY_RAISES_ASSERT(a) {} + +#endif // EIGEN_NO_ASSERTION_CHECKING + + +#define EIGEN_INTERNAL_DEBUGGING +#define EIGEN_NICE_RANDOM +#include <Eigen/Array> + + +#define VERIFY(a) do { if (!(a)) { \ + std::cerr << "Test " << g_test_stack.back() << " failed in "EI_PP_MAKE_STRING(__FILE__) << " (" << EI_PP_MAKE_STRING(__LINE__) << ")" \ + << std::endl << " " << EI_PP_MAKE_STRING(a) << std::endl << std::endl; \ + abort(); \ + } } while (0) + +#define VERIFY_IS_APPROX(a, b) VERIFY(test_ei_isApprox(a, b)) +#define VERIFY_IS_NOT_APPROX(a, b) VERIFY(!test_ei_isApprox(a, b)) +#define VERIFY_IS_MUCH_SMALLER_THAN(a, b) VERIFY(test_ei_isMuchSmallerThan(a, b)) +#define VERIFY_IS_NOT_MUCH_SMALLER_THAN(a, b) VERIFY(!test_ei_isMuchSmallerThan(a, b)) +#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 CALL_SUBTEST(FUNC) do { \ + g_test_stack.push_back(EI_PP_MAKE_STRING(FUNC)); \ + FUNC; \ + g_test_stack.pop_back(); \ + } while (0) + +namespace Eigen { + +template<typename T> inline typename NumTraits<T>::Real test_precision(); +template<> inline int test_precision<int>() { return 0; } +template<> inline float test_precision<float>() { return 1e-3f; } +template<> inline double test_precision<double>() { return 1e-6; } +template<> inline float test_precision<std::complex<float> >() { return test_precision<float>(); } +template<> inline double test_precision<std::complex<double> >() { return test_precision<double>(); } +template<> inline long double test_precision<long double>() { return 1e-6; } + +inline bool test_ei_isApprox(const int& a, const int& b) +{ return ei_isApprox(a, b, test_precision<int>()); } +inline bool test_ei_isMuchSmallerThan(const int& a, const int& b) +{ return ei_isMuchSmallerThan(a, b, test_precision<int>()); } +inline bool test_ei_isApproxOrLessThan(const int& a, const int& b) +{ return ei_isApproxOrLessThan(a, b, test_precision<int>()); } + +inline bool test_ei_isApprox(const float& a, const float& b) +{ return ei_isApprox(a, b, test_precision<float>()); } +inline bool test_ei_isMuchSmallerThan(const float& a, const float& b) +{ return ei_isMuchSmallerThan(a, b, test_precision<float>()); } +inline bool test_ei_isApproxOrLessThan(const float& a, const float& b) +{ return ei_isApproxOrLessThan(a, b, test_precision<float>()); } + +inline bool test_ei_isApprox(const double& a, const double& b) +{ return ei_isApprox(a, b, test_precision<double>()); } +inline bool test_ei_isMuchSmallerThan(const double& a, const double& b) +{ return ei_isMuchSmallerThan(a, b, test_precision<double>()); } +inline bool test_ei_isApproxOrLessThan(const double& a, const double& b) +{ return ei_isApproxOrLessThan(a, b, test_precision<double>()); } + +inline bool test_ei_isApprox(const std::complex<float>& a, const std::complex<float>& b) +{ return ei_isApprox(a, b, test_precision<std::complex<float> >()); } +inline bool test_ei_isMuchSmallerThan(const std::complex<float>& a, const std::complex<float>& b) +{ return ei_isMuchSmallerThan(a, b, test_precision<std::complex<float> >()); } + +inline bool test_ei_isApprox(const std::complex<double>& a, const std::complex<double>& b) +{ return ei_isApprox(a, b, test_precision<std::complex<double> >()); } +inline bool test_ei_isMuchSmallerThan(const std::complex<double>& a, const std::complex<double>& b) +{ return ei_isMuchSmallerThan(a, b, test_precision<std::complex<double> >()); } + +inline bool test_ei_isApprox(const long double& a, const long double& b) +{ return ei_isApprox(a, b, test_precision<long double>()); } +inline bool test_ei_isMuchSmallerThan(const long double& a, const long double& b) +{ return ei_isMuchSmallerThan(a, b, test_precision<long double>()); } +inline bool test_ei_isApproxOrLessThan(const long double& a, const long double& b) +{ return ei_isApproxOrLessThan(a, b, test_precision<long double>()); } + +template<typename Type1, typename Type2> +inline bool test_ei_isApprox(const Type1& a, const Type2& b) +{ + return a.isApprox(b, test_precision<typename Type1::Scalar>()); +} + +template<typename Derived1, typename Derived2> +inline bool test_ei_isMuchSmallerThan(const MatrixBase<Derived1>& m1, + const MatrixBase<Derived2>& m2) +{ + return m1.isMuchSmallerThan(m2, test_precision<typename ei_traits<Derived1>::Scalar>()); +} + +template<typename Derived> +inline bool test_ei_isMuchSmallerThan(const MatrixBase<Derived>& m, + const typename NumTraits<typename ei_traits<Derived>::Scalar>::Real& s) +{ + return m.isMuchSmallerThan(s, test_precision<typename ei_traits<Derived>::Scalar>()); +} + +} // end namespace Eigen + +template<typename T> struct GetDifferentType; + +template<> struct GetDifferentType<float> { typedef double type; }; +template<> struct GetDifferentType<double> { typedef float type; }; +template<typename T> struct GetDifferentType<std::complex<T> > +{ typedef std::complex<typename GetDifferentType<T>::type> type; }; + +// forward declaration of the main test function +void EI_PP_CAT(test_,EIGEN_TEST_FUNC)(); + +using namespace Eigen; + +#ifdef EIGEN_TEST_PART_1 +#define CALL_SUBTEST_1(FUNC) CALL_SUBTEST(FUNC) +#else +#define CALL_SUBTEST_1(FUNC) +#endif + +#ifdef EIGEN_TEST_PART_2 +#define CALL_SUBTEST_2(FUNC) CALL_SUBTEST(FUNC) +#else +#define CALL_SUBTEST_2(FUNC) +#endif + +#ifdef EIGEN_TEST_PART_3 +#define CALL_SUBTEST_3(FUNC) CALL_SUBTEST(FUNC) +#else +#define CALL_SUBTEST_3(FUNC) +#endif + +#ifdef EIGEN_TEST_PART_4 +#define CALL_SUBTEST_4(FUNC) CALL_SUBTEST(FUNC) +#else +#define CALL_SUBTEST_4(FUNC) +#endif + +#ifdef EIGEN_TEST_PART_5 +#define CALL_SUBTEST_5(FUNC) CALL_SUBTEST(FUNC) +#else +#define CALL_SUBTEST_5(FUNC) +#endif + +#ifdef EIGEN_TEST_PART_6 +#define CALL_SUBTEST_6(FUNC) CALL_SUBTEST(FUNC) +#else +#define CALL_SUBTEST_6(FUNC) +#endif + +#ifdef EIGEN_TEST_PART_7 +#define CALL_SUBTEST_7(FUNC) CALL_SUBTEST(FUNC) +#else +#define CALL_SUBTEST_7(FUNC) +#endif + +#ifdef EIGEN_TEST_PART_8 +#define CALL_SUBTEST_8(FUNC) CALL_SUBTEST(FUNC) +#else +#define CALL_SUBTEST_8(FUNC) +#endif + +#ifdef EIGEN_TEST_PART_9 +#define CALL_SUBTEST_9(FUNC) CALL_SUBTEST(FUNC) +#else +#define CALL_SUBTEST_9(FUNC) +#endif + +#ifdef EIGEN_TEST_PART_10 +#define CALL_SUBTEST_10(FUNC) CALL_SUBTEST(FUNC) +#else +#define CALL_SUBTEST_10(FUNC) +#endif + +#ifdef EIGEN_TEST_PART_11 +#define CALL_SUBTEST_11(FUNC) CALL_SUBTEST(FUNC) +#else +#define CALL_SUBTEST_11(FUNC) +#endif + +#ifdef EIGEN_TEST_PART_12 +#define CALL_SUBTEST_12(FUNC) CALL_SUBTEST(FUNC) +#else +#define CALL_SUBTEST_12(FUNC) +#endif + +#ifdef EIGEN_TEST_PART_13 +#define CALL_SUBTEST_13(FUNC) CALL_SUBTEST(FUNC) +#else +#define CALL_SUBTEST_13(FUNC) +#endif + +#ifdef EIGEN_TEST_PART_14 +#define CALL_SUBTEST_14(FUNC) CALL_SUBTEST(FUNC) +#else +#define CALL_SUBTEST_14(FUNC) +#endif + +#ifdef EIGEN_TEST_PART_15 +#define CALL_SUBTEST_15(FUNC) CALL_SUBTEST(FUNC) +#else +#define CALL_SUBTEST_15(FUNC) +#endif + +#ifdef EIGEN_TEST_PART_16 +#define CALL_SUBTEST_16(FUNC) CALL_SUBTEST(FUNC) +#else +#define CALL_SUBTEST_16(FUNC) +#endif + + + +int main(int argc, char *argv[]) +{ + bool has_set_repeat = false; + bool has_set_seed = false; + bool need_help = false; + unsigned int seed = 0; + int repeat = DEFAULT_REPEAT; + + for(int i = 1; i < argc; i++) + { + if(argv[i][0] == 'r') + { + if(has_set_repeat) + { + std::cout << "Argument " << argv[i] << " conflicting with a former argument" << std::endl; + return 1; + } + repeat = std::atoi(argv[i]+1); + has_set_repeat = true; + if(repeat <= 0) + { + std::cout << "Invalid \'repeat\' value " << argv[i]+1 << std::endl; + return 1; + } + } + else if(argv[i][0] == 's') + { + if(has_set_seed) + { + std::cout << "Argument " << argv[i] << " conflicting with a former argument" << std::endl; + return 1; + } + seed = int(std::strtoul(argv[i]+1, 0, 10)); + has_set_seed = true; + bool ok = seed!=0; + if(!ok) + { + std::cout << "Invalid \'seed\' value " << argv[i]+1 << std::endl; + return 1; + } + } + else + { + need_help = true; + } + } + + if(need_help) + { + std::cout << "This test application takes the following optional arguments:" << std::endl; + std::cout << " rN Repeat each test N times (default: " << DEFAULT_REPEAT << ")" << std::endl; + std::cout << " sN Use N as seed for random numbers (default: based on current time)" << std::endl; + return 1; + } + + if(!has_set_seed) seed = (unsigned int) std::time(NULL); + if(!has_set_repeat) repeat = DEFAULT_REPEAT; + + std::cout << "Initializing random number generator with seed " << seed << std::endl; + std::srand(seed); + std::cout << "Repeating each test " << repeat << " times" << std::endl; + + Eigen::g_repeat = repeat; + Eigen::g_test_stack.push_back(EI_PP_MAKE_STRING(EIGEN_TEST_FUNC)); + + EI_PP_CAT(test_,EIGEN_TEST_FUNC)(); + return 0; +} + + + diff --git a/eigen/test/eigen2/product.h b/eigen/test/eigen2/product.h new file mode 100644 index 0000000..ae1b4ba --- /dev/null +++ b/eigen/test/eigen2/product.h @@ -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) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// 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 <Eigen/Array> +#include <Eigen/QR> + +template<typename Derived1, typename Derived2> +bool areNotApprox(const MatrixBase<Derived1>& m1, const MatrixBase<Derived2>& m2, typename Derived1::RealScalar epsilon = precision<typename Derived1::RealScalar>()) +{ + return !((m1-m2).cwise().abs2().maxCoeff() < epsilon * epsilon + * std::max(m1.cwise().abs2().maxCoeff(), m2.cwise().abs2().maxCoeff())); +} + +template<typename MatrixType> void product(const MatrixType& m) +{ + /* this test covers the following files: + Identity.h Product.h + */ + + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::FloatingPoint FloatingPoint; + typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> RowVectorType; + typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> ColVectorType; + typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> RowSquareMatrixType; + typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> ColSquareMatrixType; + typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime, + MatrixType::Options^RowMajor> OtherMajorMatrixType; + + int rows = m.rows(); + int cols = m.cols(); + + // this test relies a lot on Random.h, and there's not much more that we can do + // to test it, hence I consider that we will have tested Random.h + MatrixType m1 = MatrixType::Random(rows, cols), + m2 = MatrixType::Random(rows, cols), + m3(rows, cols); + RowSquareMatrixType + identity = RowSquareMatrixType::Identity(rows, rows), + square = RowSquareMatrixType::Random(rows, rows), + res = RowSquareMatrixType::Random(rows, rows); + ColSquareMatrixType + square2 = ColSquareMatrixType::Random(cols, cols), + res2 = ColSquareMatrixType::Random(cols, cols); + RowVectorType v1 = RowVectorType::Random(rows); + ColVectorType vc2 = ColVectorType::Random(cols), vcres(cols); + OtherMajorMatrixType tm1 = m1; + + Scalar s1 = ei_random<Scalar>(); + + int r = ei_random<int>(0, rows-1), + c = ei_random<int>(0, cols-1); + + // begin testing Product.h: only associativity for now + // (we use Transpose.h but this doesn't count as a test for it) + + VERIFY_IS_APPROX((m1*m1.transpose())*m2, m1*(m1.transpose()*m2)); + m3 = m1; + m3 *= m1.transpose() * m2; + VERIFY_IS_APPROX(m3, m1 * (m1.transpose()*m2)); + VERIFY_IS_APPROX(m3, m1.lazy() * (m1.transpose()*m2)); + + // continue testing Product.h: distributivity + VERIFY_IS_APPROX(square*(m1 + m2), square*m1+square*m2); + VERIFY_IS_APPROX(square*(m1 - m2), square*m1-square*m2); + + // continue testing Product.h: compatibility with ScalarMultiple.h + VERIFY_IS_APPROX(s1*(square*m1), (s1*square)*m1); + VERIFY_IS_APPROX(s1*(square*m1), square*(m1*s1)); + + // again, test operator() to check const-qualification + s1 += (square.lazy() * m1)(r,c); + + // test Product.h together with Identity.h + VERIFY_IS_APPROX(v1, identity*v1); + VERIFY_IS_APPROX(v1.transpose(), v1.transpose() * identity); + // again, test operator() to check const-qualification + VERIFY_IS_APPROX(MatrixType::Identity(rows, cols)(r,c), static_cast<Scalar>(r==c)); + + if (rows!=cols) + VERIFY_RAISES_ASSERT(m3 = m1*m1); + + // test the previous tests were not screwed up because operator* returns 0 + // (we use the more accurate default epsilon) + if (NumTraits<Scalar>::HasFloatingPoint && std::min(rows,cols)>1) + { + VERIFY(areNotApprox(m1.transpose()*m2,m2.transpose()*m1)); + } + + // test optimized operator+= path + res = square; + res += (m1 * m2.transpose()).lazy(); + 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(); + VERIFY_IS_APPROX(vcres, vc2 + m1.transpose() * v1); + tm1 = m1; + VERIFY_IS_APPROX(tm1.transpose() * v1, m1.transpose() * v1); + VERIFY_IS_APPROX(v1.transpose() * tm1, v1.transpose() * m1); + + // test submatrix and matrix/vector product + for (int i=0; i<rows; ++i) + res.row(i) = m1.row(i) * m2.transpose(); + VERIFY_IS_APPROX(res, m1 * m2.transpose()); + // the other way round: + for (int i=0; i<rows; ++i) + res.col(i) = m1 * m2.transpose().col(i); + VERIFY_IS_APPROX(res, m1 * m2.transpose()); + + res2 = square2; + res2 += (m1.transpose() * m2).lazy(); + VERIFY_IS_APPROX(res2, square2 + m1.transpose() * m2); + + if (NumTraits<Scalar>::HasFloatingPoint && std::min(rows,cols)>1) + { + VERIFY(areNotApprox(res2,square2 + m2.transpose() * m1)); + } +} + diff --git a/eigen/test/eigen2/runtest.sh b/eigen/test/eigen2/runtest.sh new file mode 100644 index 0000000..bc693af --- /dev/null +++ b/eigen/test/eigen2/runtest.sh @@ -0,0 +1,28 @@ +#!/bin/bash + +black='\E[30m' +red='\E[31m' +green='\E[32m' +yellow='\E[33m' +blue='\E[34m' +magenta='\E[35m' +cyan='\E[36m' +white='\E[37m' + +if make test_$1 > /dev/null 2> .runtest.log ; then + if ! ./test_$1 r20 > /dev/null 2> .runtest.log ; then + echo -e $red Test $1 failed: $black + echo -e $blue + cat .runtest.log + echo -e $black + exit 1 + else + echo -e $green Test $1 passed$black + fi +else + echo -e $red Build of target $1 failed: $black + echo -e $blue + cat .runtest.log + echo -e $black + exit 1 +fi diff --git a/eigen/test/eigen2/sparse.h b/eigen/test/eigen2/sparse.h new file mode 100644 index 0000000..e12f899 --- /dev/null +++ b/eigen/test/eigen2/sparse.h @@ -0,0 +1,154 @@ +// 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 Daniel Gomez Ferro <dgomezferro@gmail.com> +// +// 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_TESTSPARSE_H + +#include "main.h" + +#if EIGEN_GNUC_AT_LEAST(4,0) && !defined __ICC +#include <tr1/unordered_map> +#define EIGEN_UNORDERED_MAP_SUPPORT +namespace std { + using std::tr1::unordered_map; +} +#endif + +#ifdef EIGEN_GOOGLEHASH_SUPPORT + #include <google/sparse_hash_map> +#endif + +#include <Eigen/Cholesky> +#include <Eigen/LU> +#include <Eigen/Sparse> + +enum { + ForceNonZeroDiag = 1, + MakeLowerTriangular = 2, + MakeUpperTriangular = 4, + ForceRealDiag = 8 +}; + +/* Initializes both a sparse and dense matrix with same random values, + * and a ratio of \a density non zero entries. + * \param flags is a union of ForceNonZeroDiag, MakeLowerTriangular and MakeUpperTriangular + * allowing to control the shape of the matrix. + * \param zeroCoords and nonzeroCoords allows to get the coordinate lists of the non zero, + * and zero coefficients respectively. + */ +template<typename Scalar> void +initSparse(double density, + Matrix<Scalar,Dynamic,Dynamic>& refMat, + SparseMatrix<Scalar>& sparseMat, + int flags = 0, + std::vector<Vector2i>* zeroCoords = 0, + std::vector<Vector2i>* nonzeroCoords = 0) +{ + sparseMat.startFill(int(refMat.rows()*refMat.cols()*density)); + for(int j=0; j<refMat.cols(); j++) + { + for(int i=0; i<refMat.rows(); i++) + { + Scalar v = (ei_random<double>(0,1) < density) ? ei_random<Scalar>() : Scalar(0); + if ((flags&ForceNonZeroDiag) && (i==j)) + { + v = ei_random<Scalar>()*Scalar(3.); + v = v*v + Scalar(5.); + } + if ((flags & MakeLowerTriangular) && j>i) + v = Scalar(0); + else if ((flags & MakeUpperTriangular) && j<i) + v = Scalar(0); + + if ((flags&ForceRealDiag) && (i==j)) + v = ei_real(v); + + if (v!=Scalar(0)) + { + sparseMat.fill(i,j) = v; + if (nonzeroCoords) + nonzeroCoords->push_back(Vector2i(i,j)); + } + else if (zeroCoords) + { + zeroCoords->push_back(Vector2i(i,j)); + } + refMat(i,j) = v; + } + } + sparseMat.endFill(); +} + +template<typename Scalar> void +initSparse(double density, + Matrix<Scalar,Dynamic,Dynamic>& refMat, + DynamicSparseMatrix<Scalar>& sparseMat, + int flags = 0, + std::vector<Vector2i>* zeroCoords = 0, + std::vector<Vector2i>* nonzeroCoords = 0) +{ + sparseMat.startFill(int(refMat.rows()*refMat.cols()*density)); + for(int j=0; j<refMat.cols(); j++) + { + for(int i=0; i<refMat.rows(); i++) + { + Scalar v = (ei_random<double>(0,1) < density) ? ei_random<Scalar>() : Scalar(0); + if ((flags&ForceNonZeroDiag) && (i==j)) + { + v = ei_random<Scalar>()*Scalar(3.); + v = v*v + Scalar(5.); + } + if ((flags & MakeLowerTriangular) && j>i) + v = Scalar(0); + else if ((flags & MakeUpperTriangular) && j<i) + v = Scalar(0); + + if ((flags&ForceRealDiag) && (i==j)) + v = ei_real(v); + + if (v!=Scalar(0)) + { + sparseMat.fill(i,j) = v; + if (nonzeroCoords) + nonzeroCoords->push_back(Vector2i(i,j)); + } + else if (zeroCoords) + { + zeroCoords->push_back(Vector2i(i,j)); + } + refMat(i,j) = v; + } + } + sparseMat.endFill(); +} + +template<typename Scalar> void +initSparse(double density, + Matrix<Scalar,Dynamic,1>& refVec, + SparseVector<Scalar>& sparseVec, + std::vector<int>* zeroCoords = 0, + std::vector<int>* nonzeroCoords = 0) +{ + sparseVec.reserve(int(refVec.size()*density)); + sparseVec.setZero(); + for(int i=0; i<refVec.size(); i++) + { + Scalar v = (ei_random<double>(0,1) < density) ? ei_random<Scalar>() : Scalar(0); + if (v!=Scalar(0)) + { + sparseVec.fill(i) = v; + if (nonzeroCoords) + nonzeroCoords->push_back(i); + } + else if (zeroCoords) + zeroCoords->push_back(i); + refVec[i] = v; + } +} + +#endif // EIGEN_TESTSPARSE_H diff --git a/eigen/test/eigen2/testsuite.cmake b/eigen/test/eigen2/testsuite.cmake new file mode 100644 index 0000000..12b6bfa --- /dev/null +++ b/eigen/test/eigen2/testsuite.cmake @@ -0,0 +1,197 @@ + +#################################################################### +# +# Usage: +# - create a new folder, let's call it cdash +# - in that folder, do: +# ctest -S path/to/eigen2/test/testsuite.cmake[,option1=value1[,option2=value2]] +# +# Options: +# - EIGEN_CXX: compiler, eg.: g++-4.2 +# default: default c++ compiler +# - EIGEN_SITE: eg, INRIA-Bdx_pc-gael, or the name of the contributor, etc. +# default: hostname +# - EIGEN_BUILD_STRING: a string which identify the system/compiler. It should be formed like that: +# <OS_name>-<OS_version>-<arch>-<compiler-version> +# with: +# <OS_name> = opensuse, debian, osx, windows, cygwin, freebsd, solaris, etc. +# <OS_version> = 11.1, XP, vista, leopard, etc. +# <arch> = i386, x86_64, ia64, powerpc, etc. +# <compiler-version> = gcc-4.3.2, icc-11.0, MSVC-2008, etc. +# - EIGEN_EXPLICIT_VECTORIZATION: novec, SSE2, Altivec +# default: SSE2 for x86_64 systems, novec otherwise +# Its value is automatically appended to EIGEN_BUILD_STRING +# - EIGEN_CMAKE_DIR: path to cmake executable +# - EIGEN_MODE: dashboard model, can be Experimental, Nightly, or Continuous +# default: Nightly +# - 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 +# - 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 +# default: <EIGEN_WORK_DIR>/nightly-<EIGEN_CXX> +# +# Here is an example running several compilers on a linux system: +# #!/bin/bash +# ARCH=`uname -m` +# SITE=`hostname` +# VERSION=opensuse-11.1 +# WORK_DIR=/home/gael/Coding/eigen2/cdash +# # get the last version of the script +# wget http://bitbucket.org/eigen/eigen/raw/tip/test/testsuite.cmake -o $WORK_DIR/testsuite.cmake +# COMMON="ctest -S $WORK_DIR/testsuite.cmake,EIGEN_WORK_DIR=$WORK_DIR,EIGEN_SITE=$SITE,EIGEN_MODE=$1,EIGEN_BUILD_STRING=$OS_VERSION-$ARCH" +# $COMMON-gcc-3.4.6,EIGEN_CXX=g++-3.4 +# $COMMON-gcc-4.0.1,EIGEN_CXX=g++-4.0.1 +# $COMMON-gcc-4.3.2,EIGEN_CXX=g++-4.3,EIGEN_EXPLICIT_VECTORIZATION=novec +# $COMMON-gcc-4.3.2,EIGEN_CXX=g++-4.3,EIGEN_EXPLICIT_VECTORIZATION=SSE2 +# $COMMON-icc-11.0,EIGEN_CXX=icpc +# +#################################################################### + +# process the arguments + +set(ARGLIST ${CTEST_SCRIPT_ARG}) +while(${ARGLIST} MATCHES ".+.*") + + # pick first + string(REGEX MATCH "([^,]*)(,.*)?" DUMMY ${ARGLIST}) + SET(TOP ${CMAKE_MATCH_1}) + + # remove first + string(REGEX MATCHALL "[^,]*,(.*)" DUMMY ${ARGLIST}) + SET(ARGLIST ${CMAKE_MATCH_1}) + + # decompose as a pair key=value + string(REGEX MATCH "([^=]*)(=.*)?" DUMMY ${TOP}) + SET(KEY ${CMAKE_MATCH_1}) + + string(REGEX MATCH "[^=]*=(.*)" DUMMY ${TOP}) + SET(VALUE ${CMAKE_MATCH_1}) + + # set the variable to the specified value + if(VALUE) + SET(${KEY} ${VALUE}) + else(VALUE) + SET(${KEY} ON) + endif(VALUE) + +endwhile(${ARGLIST} MATCHES ".+.*") + +#################################################################### +# Automatically set some user variables if they have not been defined manually +#################################################################### +cmake_minimum_required(VERSION 2.6 FATAL_ERROR) + +if(NOT EIGEN_SITE) + site_name(EIGEN_SITE) +endif(NOT EIGEN_SITE) + +if(NOT EIGEN_CMAKE_DIR) + SET(EIGEN_CMAKE_DIR "") +endif(NOT EIGEN_CMAKE_DIR) + +if(NOT EIGEN_BUILD_STRING) + + # let's try to find all information we need to make the build string ourself + + # OS + build_name(EIGEN_OS_VERSION) + + # arch + set(EIGEN_ARCH ${CMAKE_SYSTEM_PROCESSOR}) + if(WIN32) + set(EIGEN_ARCH $ENV{PROCESSOR_ARCHITECTURE}) + else(WIN32) + execute_process(COMMAND uname -m OUTPUT_VARIABLE EIGEN_ARCH OUTPUT_STRIP_TRAILING_WHITESPACE) + endif(WIN32) + + set(EIGEN_BUILD_STRING ${EIGEN_OS_VERSION}${EIGEN_ARCH}-${EIGEN_CXX}) + +endif(NOT EIGEN_BUILD_STRING) + +if(DEFINED EIGEN_EXPLICIT_VECTORIZATION) + set(EIGEN_BUILD_STRING ${EIGEN_BUILD_STRING}-${EIGEN_EXPLICIT_VECTORIZATION}) +endif(DEFINED EIGEN_EXPLICIT_VECTORIZATION) + +if(NOT EIGEN_WORK_DIR) + set(EIGEN_WORK_DIR ${CTEST_SCRIPT_DIRECTORY}) +endif(NOT EIGEN_WORK_DIR) + +if(NOT CTEST_SOURCE_DIRECTORY) + SET (CTEST_SOURCE_DIRECTORY "${EIGEN_WORK_DIR}/src") +endif(NOT CTEST_SOURCE_DIRECTORY) + +if(NOT CTEST_BINARY_DIRECTORY) + SET (CTEST_BINARY_DIRECTORY "${EIGEN_WORK_DIR}/nightly_${EIGEN_CXX}") +endif(NOT CTEST_BINARY_DIRECTORY) + +if(NOT EIGEN_MODE) + set(EIGEN_MODE Nightly) +endif(NOT EIGEN_MODE) + +## mandatory variables (the default should be ok in most cases): + +SET (CTEST_CVS_COMMAND "hg") +SET (CTEST_CVS_CHECKOUT "${CTEST_CVS_COMMAND} clone -r 2.0 http://bitbucket.org/eigen/eigen \"${CTEST_SOURCE_DIRECTORY}\"") + +# which ctest command to use for running the dashboard +SET (CTEST_COMMAND "${EIGEN_CMAKE_DIR}ctest -D ${EIGEN_MODE}") + +# what cmake command to use for configuring this dashboard +SET (CTEST_CMAKE_COMMAND "${EIGEN_CMAKE_DIR}cmake -DEIGEN_BUILD_TESTS=on ") + +#################################################################### +# The values in this section are optional you can either +# have them or leave them commented out +#################################################################### + +# this make sure we get consistent outputs +SET($ENV{LC_MESSAGES} "en_EN") + +# should ctest wipe the binary tree before running +SET(CTEST_START_WITH_EMPTY_BINARY_DIRECTORY TRUE) +SET(CTEST_BACKUP_AND_RESTORE TRUE) + +# this is the initial cache to use for the binary tree, be careful to escape +# 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 + BUILDNAME:STRING=${EIGEN_BUILD_STRING} + SITE:STRING=${EIGEN_SITE} + ") +else(WIN32 AND NOT UNIX) + SET (CTEST_INITIAL_CACHE " + BUILDNAME:STRING=${EIGEN_BUILD_STRING} + SITE:STRING=${EIGEN_SITE} + ") +endif(WIN32 AND NOT UNIX) + +# set any extra environment variables to use during the execution of the script here: + +if(EIGEN_CXX) + set(CTEST_ENVIRONMENT "CXX=${EIGEN_CXX}") +endif(EIGEN_CXX) + +if(DEFINED EIGEN_EXPLICIT_VECTORIZATION) + if(EIGEN_EXPLICIT_VECTORIZATION MATCHES SSE2) + set(CTEST_CMAKE_COMMAND "${CTEST_CMAKE_COMMAND} -DEIGEN_TEST_SSE2=ON") + elseif(EIGEN_EXPLICIT_VECTORIZATION MATCHES SSE3) + set(CTEST_CMAKE_COMMAND "${CTEST_CMAKE_COMMAND} -DEIGEN_TEST_SSE2=ON -DEIGEN_TEST_SSE3=ON") + elseif(EIGEN_EXPLICIT_VECTORIZATION MATCHES Altivec) + set(CTEST_CMAKE_COMMAND "${CTEST_CMAKE_COMMAND} -DEIGEN_TEST_ALTIVEC=ON") + elseif(EIGEN_EXPLICIT_VECTORIZATION MATCHES novec) + set(CTEST_CMAKE_COMMAND "${CTEST_CMAKE_COMMAND} -DEIGEN_TEST_NO_EXPLICIT_VECTORIZATION=ON") + else(EIGEN_EXPLICIT_VECTORIZATION MATCHES SSE2) + message(FATAL_ERROR "Invalid value for EIGEN_EXPLICIT_VECTORIZATION (${EIGEN_EXPLICIT_VECTORIZATION}), must be: novec, SSE2, SSE3, Altivec") + endif(EIGEN_EXPLICIT_VECTORIZATION MATCHES SSE2) +endif(DEFINED EIGEN_EXPLICIT_VECTORIZATION) + +if(DEFINED EIGEN_CMAKE_ARGS) + set(CTEST_CMAKE_COMMAND "${CTEST_CMAKE_COMMAND} ${EIGEN_CMAKE_ARGS}") +endif(DEFINED EIGEN_CMAKE_ARGS) |