Squashed 'third_party/eigen/' content from commit 61d72f6
Change-Id: Iccc90fa0b55ab44037f018046d2fcffd90d9d025
git-subtree-dir: third_party/eigen
git-subtree-split: 61d72f6383cfa842868c53e30e087b0258177257
diff --git a/test/eigen2/CMakeLists.txt b/test/eigen2/CMakeLists.txt
new file mode 100644
index 0000000..9615a60
--- /dev/null
+++ b/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/test/eigen2/eigen2_adjoint.cpp b/test/eigen2/eigen2_adjoint.cpp
new file mode 100644
index 0000000..c0f8114
--- /dev/null
+++ b/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/test/eigen2/eigen2_alignedbox.cpp b/test/eigen2/eigen2_alignedbox.cpp
new file mode 100644
index 0000000..35043b9
--- /dev/null
+++ b/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/test/eigen2/eigen2_array.cpp b/test/eigen2/eigen2_array.cpp
new file mode 100644
index 0000000..c1ff40c
--- /dev/null
+++ b/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/test/eigen2/eigen2_basicstuff.cpp b/test/eigen2/eigen2_basicstuff.cpp
new file mode 100644
index 0000000..dd2dec1
--- /dev/null
+++ b/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/test/eigen2/eigen2_bug_132.cpp b/test/eigen2/eigen2_bug_132.cpp
new file mode 100644
index 0000000..7fe3610
--- /dev/null
+++ b/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/test/eigen2/eigen2_cholesky.cpp b/test/eigen2/eigen2_cholesky.cpp
new file mode 100644
index 0000000..9c4b6f5
--- /dev/null
+++ b/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/test/eigen2/eigen2_commainitializer.cpp b/test/eigen2/eigen2_commainitializer.cpp
new file mode 100644
index 0000000..e0f901e
--- /dev/null
+++ b/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/test/eigen2/eigen2_cwiseop.cpp b/test/eigen2/eigen2_cwiseop.cpp
new file mode 100644
index 0000000..22e1cc3
--- /dev/null
+++ b/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/test/eigen2/eigen2_determinant.cpp b/test/eigen2/eigen2_determinant.cpp
new file mode 100644
index 0000000..c7b4ad0
--- /dev/null
+++ b/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/test/eigen2/eigen2_dynalloc.cpp b/test/eigen2/eigen2_dynalloc.cpp
new file mode 100644
index 0000000..1891a9e
--- /dev/null
+++ b/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/test/eigen2/eigen2_eigensolver.cpp b/test/eigen2/eigen2_eigensolver.cpp
new file mode 100644
index 0000000..48b4ace
--- /dev/null
+++ b/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/test/eigen2/eigen2_first_aligned.cpp b/test/eigen2/eigen2_first_aligned.cpp
new file mode 100644
index 0000000..51bb3ca
--- /dev/null
+++ b/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/test/eigen2/eigen2_geometry.cpp b/test/eigen2/eigen2_geometry.cpp
new file mode 100644
index 0000000..5140407
--- /dev/null
+++ b/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/test/eigen2/eigen2_geometry_with_eigen2_prefix.cpp b/test/eigen2/eigen2_geometry_with_eigen2_prefix.cpp
new file mode 100644
index 0000000..12d4a71
--- /dev/null
+++ b/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/test/eigen2/eigen2_hyperplane.cpp b/test/eigen2/eigen2_hyperplane.cpp
new file mode 100644
index 0000000..f3f85e1
--- /dev/null
+++ b/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/test/eigen2/eigen2_inverse.cpp b/test/eigen2/eigen2_inverse.cpp
new file mode 100644
index 0000000..ccd24a1
--- /dev/null
+++ b/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/test/eigen2/eigen2_linearstructure.cpp b/test/eigen2/eigen2_linearstructure.cpp
new file mode 100644
index 0000000..488f4c4
--- /dev/null
+++ b/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/test/eigen2/eigen2_lu.cpp b/test/eigen2/eigen2_lu.cpp
new file mode 100644
index 0000000..e993b1c
--- /dev/null
+++ b/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/test/eigen2/eigen2_map.cpp b/test/eigen2/eigen2_map.cpp
new file mode 100644
index 0000000..4a1c4e1
--- /dev/null
+++ b/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/test/eigen2/eigen2_meta.cpp b/test/eigen2/eigen2_meta.cpp
new file mode 100644
index 0000000..1d01bd8
--- /dev/null
+++ b/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/test/eigen2/eigen2_miscmatrices.cpp b/test/eigen2/eigen2_miscmatrices.cpp
new file mode 100644
index 0000000..8bbb20c
--- /dev/null
+++ b/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/test/eigen2/eigen2_mixingtypes.cpp b/test/eigen2/eigen2_mixingtypes.cpp
new file mode 100644
index 0000000..fb5ac5d
--- /dev/null
+++ b/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/test/eigen2/eigen2_newstdvector.cpp b/test/eigen2/eigen2_newstdvector.cpp
new file mode 100644
index 0000000..5f90099
--- /dev/null
+++ b/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/test/eigen2/eigen2_nomalloc.cpp b/test/eigen2/eigen2_nomalloc.cpp
new file mode 100644
index 0000000..d34a699
--- /dev/null
+++ b/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/test/eigen2/eigen2_packetmath.cpp b/test/eigen2/eigen2_packetmath.cpp
new file mode 100644
index 0000000..b1f325f
--- /dev/null
+++ b/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/test/eigen2/eigen2_parametrizedline.cpp b/test/eigen2/eigen2_parametrizedline.cpp
new file mode 100644
index 0000000..8147288
--- /dev/null
+++ b/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/test/eigen2/eigen2_prec_inverse_4x4.cpp b/test/eigen2/eigen2_prec_inverse_4x4.cpp
new file mode 100644
index 0000000..8bfa556
--- /dev/null
+++ b/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/test/eigen2/eigen2_product_large.cpp b/test/eigen2/eigen2_product_large.cpp
new file mode 100644
index 0000000..5149ef7
--- /dev/null
+++ b/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/test/eigen2/eigen2_product_small.cpp b/test/eigen2/eigen2_product_small.cpp
new file mode 100644
index 0000000..4cd8c10
--- /dev/null
+++ b/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/test/eigen2/eigen2_qr.cpp b/test/eigen2/eigen2_qr.cpp
new file mode 100644
index 0000000..76977e4
--- /dev/null
+++ b/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/test/eigen2/eigen2_qtvector.cpp b/test/eigen2/eigen2_qtvector.cpp
new file mode 100644
index 0000000..6cfb58a
--- /dev/null
+++ b/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/test/eigen2/eigen2_regression.cpp b/test/eigen2/eigen2_regression.cpp
new file mode 100644
index 0000000..c68b58d
--- /dev/null
+++ b/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/test/eigen2/eigen2_sizeof.cpp b/test/eigen2/eigen2_sizeof.cpp
new file mode 100644
index 0000000..ec1af5a
--- /dev/null
+++ b/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/test/eigen2/eigen2_smallvectors.cpp b/test/eigen2/eigen2_smallvectors.cpp
new file mode 100644
index 0000000..03962b1
--- /dev/null
+++ b/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/test/eigen2/eigen2_sparse_basic.cpp b/test/eigen2/eigen2_sparse_basic.cpp
new file mode 100644
index 0000000..0490776
--- /dev/null
+++ b/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/test/eigen2/eigen2_sparse_product.cpp b/test/eigen2/eigen2_sparse_product.cpp
new file mode 100644
index 0000000..d28e76d
--- /dev/null
+++ b/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/test/eigen2/eigen2_sparse_solvers.cpp b/test/eigen2/eigen2_sparse_solvers.cpp
new file mode 100644
index 0000000..3aef27a
--- /dev/null
+++ b/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/test/eigen2/eigen2_sparse_vector.cpp b/test/eigen2/eigen2_sparse_vector.cpp
new file mode 100644
index 0000000..e6d2d77
--- /dev/null
+++ b/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/test/eigen2/eigen2_stdvector.cpp b/test/eigen2/eigen2_stdvector.cpp
new file mode 100644
index 0000000..6ab05c2
--- /dev/null
+++ b/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/test/eigen2/eigen2_submatrices.cpp b/test/eigen2/eigen2_submatrices.cpp
new file mode 100644
index 0000000..dee970b
--- /dev/null
+++ b/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/test/eigen2/eigen2_sum.cpp b/test/eigen2/eigen2_sum.cpp
new file mode 100644
index 0000000..b47057c
--- /dev/null
+++ b/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/test/eigen2/eigen2_svd.cpp b/test/eigen2/eigen2_svd.cpp
new file mode 100644
index 0000000..d4689a5
--- /dev/null
+++ b/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/test/eigen2/eigen2_swap.cpp b/test/eigen2/eigen2_swap.cpp
new file mode 100644
index 0000000..f3a8846
--- /dev/null
+++ b/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/test/eigen2/eigen2_triangular.cpp b/test/eigen2/eigen2_triangular.cpp
new file mode 100644
index 0000000..6f17b7d
--- /dev/null
+++ b/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/test/eigen2/eigen2_unalignedassert.cpp b/test/eigen2/eigen2_unalignedassert.cpp
new file mode 100644
index 0000000..d10b6f5
--- /dev/null
+++ b/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/test/eigen2/eigen2_visitor.cpp b/test/eigen2/eigen2_visitor.cpp
new file mode 100644
index 0000000..4781991
--- /dev/null
+++ b/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/test/eigen2/gsl_helper.h b/test/eigen2/gsl_helper.h
new file mode 100644
index 0000000..d1d8545
--- /dev/null
+++ b/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/test/eigen2/main.h b/test/eigen2/main.h
new file mode 100644
index 0000000..ad2ba19
--- /dev/null
+++ b/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/test/eigen2/product.h b/test/eigen2/product.h
new file mode 100644
index 0000000..ae1b4ba
--- /dev/null
+++ b/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/test/eigen2/runtest.sh b/test/eigen2/runtest.sh
new file mode 100755
index 0000000..bc693af
--- /dev/null
+++ b/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/test/eigen2/sparse.h b/test/eigen2/sparse.h
new file mode 100644
index 0000000..e12f899
--- /dev/null
+++ b/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/test/eigen2/testsuite.cmake b/test/eigen2/testsuite.cmake
new file mode 100644
index 0000000..12b6bfa
--- /dev/null
+++ b/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)