diff options
| -rw-r--r-- | test/CMakeLists.txt | 1 | ||||
| -rw-r--r-- | test/jacobi.cpp | 93 | ||||
| -rw-r--r-- | test/jacobisvd.cpp | 1 |
3 files changed, 94 insertions, 1 deletions
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 95115e9d5..052d4ec13 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -91,6 +91,7 @@ ei_add_test(eigensolver_selfadjoint " " "${GSL_LIBRARIES}") ei_add_test(eigensolver_generic " " "${GSL_LIBRARIES}") ei_add_test(eigensolver_complex) ei_add_test(svd) +ei_add_test(jacobi) ei_add_test(jacobisvd) ei_add_test(geo_orthomethods) ei_add_test(geo_homogeneous) diff --git a/test/jacobi.cpp b/test/jacobi.cpp new file mode 100644 index 000000000..eb839fa77 --- /dev/null +++ b/test/jacobi.cpp @@ -0,0 +1,93 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#include "main.h" +#include <Eigen/SVD> + +template<typename MatrixType, typename JacobiScalar> +void jacobi(const MatrixType& m = MatrixType()) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::Index Index; + Index rows = m.rows(); + Index cols = m.cols(); + + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime + }; + + typedef Matrix<JacobiScalar, 2, 1> JacobiVector; + + const MatrixType a(MatrixType::Random(rows, cols)); + + JacobiVector v = JacobiVector::Random().normalized(); + JacobiScalar c = v.x(), s = v.y(); + PlanarRotation<JacobiScalar> rot(c, s); + + { + Index p = ei_random<Index>(0, rows-1); + Index q; + do { + q = ei_random<Index>(0, rows-1); + } while (q == p); + + MatrixType b = a; + b.applyOnTheLeft(p, q, rot); + VERIFY_IS_APPROX(b.row(p), c * a.row(p) + ei_conj(s) * a.row(q)); + VERIFY_IS_APPROX(b.row(q), -s * a.row(p) + ei_conj(c) * a.row(q)); + } + + { + Index p = ei_random<Index>(0, cols-1); + Index q; + do { + q = ei_random<Index>(0, cols-1); + } while (q == p); + + MatrixType b = a; + b.applyOnTheRight(p, q, rot); + VERIFY_IS_APPROX(b.col(p), c * a.col(p) - s * a.col(q)); + VERIFY_IS_APPROX(b.col(q), ei_conj(s) * a.col(p) + ei_conj(c) * a.col(q)); + } +} + +void test_jacobi() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1(( jacobi<Matrix3f, float>() )); + CALL_SUBTEST_2(( jacobi<Matrix4d, double>() )); + CALL_SUBTEST_3(( jacobi<Matrix4cf, float>() )); + CALL_SUBTEST_3(( jacobi<Matrix4cf, std::complex<float> >() )); + + int r = ei_random<int>(2, 20), + c = ei_random<int>(2, 20); + CALL_SUBTEST_4(( jacobi<MatrixXf, float>(MatrixXf(r,c)) )); + CALL_SUBTEST_5(( jacobi<MatrixXcd, double>(MatrixXcd(r,c)) )); + CALL_SUBTEST_5(( jacobi<MatrixXcd, std::complex<double> >(MatrixXcd(r,c)) )); + (void) r; + (void) c; + } +} diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp index 0039be388..18ca86c2d 100644 --- a/test/jacobisvd.cpp +++ b/test/jacobisvd.cpp @@ -25,7 +25,6 @@ #include "main.h" #include <Eigen/SVD> -#include <Eigen/LU> template<typename MatrixType, int QRPreconditioner> void jacobisvd_check_full(const MatrixType& m, const JacobiSVD<MatrixType, QRPreconditioner>& svd) |