diff options
author | Chen-Pang He <jdh8@ms63.hinet.net> | 2012-08-27 22:48:37 +0100 |
---|---|---|
committer | Chen-Pang He <jdh8@ms63.hinet.net> | 2012-08-27 22:48:37 +0100 |
commit | bfaa7f4ffeb55f91278d70cd56659ce866e6ef88 (patch) | |
tree | ba5e0e225845e2e4ad9d401dceb6bdb03be56dce /unsupported | |
parent | b55d260adaadaece9ed92973792c4cc846061881 (diff) |
Add test for matrix power.
Use Christoph Hertzberg's suggestion to use exponent laws.
Diffstat (limited to 'unsupported')
-rw-r--r-- | unsupported/test/CMakeLists.txt | 1 | ||||
-rw-r--r-- | unsupported/test/matrix_exponential.cpp | 12 | ||||
-rw-r--r-- | unsupported/test/matrix_functions.h | 47 | ||||
-rw-r--r-- | unsupported/test/matrix_power.cpp | 104 | ||||
-rw-r--r-- | unsupported/test/matrix_square_root.cpp | 33 |
5 files changed, 155 insertions, 42 deletions
diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt index b34b151b1..ff0137ec6 100644 --- a/unsupported/test/CMakeLists.txt +++ b/unsupported/test/CMakeLists.txt @@ -33,6 +33,7 @@ endif() ei_add_test(matrix_exponential) ei_add_test(matrix_function) +ei_add_test(matrix_power) ei_add_test(matrix_square_root) ei_add_test(alignedvector3) ei_add_test(FFT) diff --git a/unsupported/test/matrix_exponential.cpp b/unsupported/test/matrix_exponential.cpp index 695472f91..50dec083d 100644 --- a/unsupported/test/matrix_exponential.cpp +++ b/unsupported/test/matrix_exponential.cpp @@ -7,8 +7,7 @@ // 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 <unsupported/Eigen/MatrixFunctions> +#include "matrix_functions.h" double binom(int n, int k) { @@ -18,12 +17,6 @@ double binom(int n, int k) return res; } -template <typename Derived, typename OtherDerived> -double relerr(const MatrixBase<Derived>& A, const MatrixBase<OtherDerived>& B) -{ - return std::sqrt((A - B).cwiseAbs2().sum() / (std::min)(A.cwiseAbs2().sum(), B.cwiseAbs2().sum())); -} - template <typename T> T expfn(T x, int) { @@ -109,8 +102,7 @@ void randomTest(const MatrixType& m, double tol) */ typename MatrixType::Index rows = m.rows(); typename MatrixType::Index cols = m.cols(); - MatrixType m1(rows, cols), m2(rows, cols), m3(rows, cols), - identity = MatrixType::Identity(rows, rows); + MatrixType m1(rows, cols), m2(rows, cols), identity = MatrixType::Identity(rows, cols); typedef typename NumTraits<typename internal::traits<MatrixType>::Scalar>::Real RealScalar; diff --git a/unsupported/test/matrix_functions.h b/unsupported/test/matrix_functions.h new file mode 100644 index 000000000..5817caef6 --- /dev/null +++ b/unsupported/test/matrix_functions.h @@ -0,0 +1,47 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009-2011 Jitse Niesen <jitse@maths.leeds.ac.uk> +// +// 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 <unsupported/Eigen/MatrixFunctions> + +template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex> +struct generateTestMatrix; + +// for real matrices, make sure none of the eigenvalues are negative +template <typename MatrixType> +struct generateTestMatrix<MatrixType,0> +{ + static void run(MatrixType& result, typename MatrixType::Index size) + { + MatrixType mat = MatrixType::Random(size, size); + EigenSolver<MatrixType> es(mat); + typename EigenSolver<MatrixType>::EigenvalueType eivals = es.eigenvalues(); + for (typename MatrixType::Index i = 0; i < size; ++i) { + if (eivals(i).imag() == 0 && eivals(i).real() < 0) + eivals(i) = -eivals(i); + } + result = (es.eigenvectors() * eivals.asDiagonal() * es.eigenvectors().inverse()).real(); + } +}; + +// for complex matrices, any matrix is fine +template <typename MatrixType> +struct generateTestMatrix<MatrixType,1> +{ + static void run(MatrixType& result, typename MatrixType::Index size) + { + result = MatrixType::Random(size, size); + } +}; + +template <typename Derived, typename OtherDerived> +double relerr(const MatrixBase<Derived>& A, const MatrixBase<OtherDerived>& B) +{ + return std::sqrt((A - B).cwiseAbs2().sum() / (std::min)(A.cwiseAbs2().sum(), B.cwiseAbs2().sum())); +} diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp new file mode 100644 index 000000000..f3ef57157 --- /dev/null +++ b/unsupported/test/matrix_power.cpp @@ -0,0 +1,104 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net> +// +// 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 "matrix_functions.h" + +template <typename T> +void test2dRotation(double tol) +{ + Matrix<T,2,2> A, B, C; + T angle, c, s; + + A << 0, 1, -1, 0; + for (int i = 0; i <= 20; i++) { + angle = pow(10, (i-10) / 5.); + c = std::cos(angle); + s = std::sin(angle); + B << c, s, -s, c; + + C = A.pow(std::ldexp(angle, 1) / M_PI); + std::cout << "test2dRotation: i = " << i << " error powerm = " << relerr(C, B) << "\n"; + VERIFY(C.isApprox(B, T(tol))); + } +} + +template <typename T> +void test2dHyperbolicRotation(double tol) +{ + Matrix<std::complex<T>,2,2> A, B, C; + T angle, ch = std::cosh(1); + std::complex<T> ish(0, std::sinh(1)); + + A << ch, ish, -ish, ch; + for (int i = 0; i <= 20; i++) { + angle = std::ldexp(T(i-10), -1); + ch = std::cosh(angle); + ish = std::complex<T>(0, std::sinh(angle)); + B << ch, ish, -ish, ch; + + C = A.pow(angle); + std::cout << "test2dHyperbolicRotation: i = " << i << " error powerm = " << relerr(C, B) << "\n"; + VERIFY(C.isApprox(B, T(tol))); + } +} + +template <typename MatrixType> +void testExponentLaws(const MatrixType& m, double tol) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + + typename MatrixType::Index rows = m.rows(); + typename MatrixType::Index cols = m.cols(); + MatrixType m1, m1x, m1y, m2, m3; + RealScalar x = internal::random<RealScalar>(), y = internal::random<RealScalar>(); + double err[3]; + + for(int i = 0; i < g_repeat; i++) { + generateTestMatrix<MatrixType>::run(m1, m.rows()); + m1x = m1.pow(x); + m1y = m1.pow(y); + + m2 = m1.pow(x + y); + m3 = m1x * m1y; + err[0] = relerr(m2, m3); + VERIFY(m2.isApprox(m3, static_cast<RealScalar>(tol))); + + m2 = m1.pow(x * y); + m3 = m1x.pow(y); + err[1] = relerr(m2, m3); + VERIFY(m2.isApprox(m3, static_cast<RealScalar>(tol))); + + m2 = (std::abs(x) * m1).pow(y); + m3 = std::pow(std::abs(x), y) * m1y; + err[2] = relerr(m2, m3); + VERIFY(m2.isApprox(m3, static_cast<RealScalar>(tol))); + + std::cout << "testExponentLaws: error powerm = " << err[0] << " " << err[1] << " " << err[2] << "\n"; + } +} + +void test_matrix_power() +{ + CALL_SUBTEST_2(test2dRotation<double>(1e-13)); + CALL_SUBTEST_1(test2dRotation<float>(2e-5)); // was 1e-5, relaxed for clang 2.8 / linux / x86-64 + CALL_SUBTEST_8(test2dRotation<long double>(1e-13)); + CALL_SUBTEST_2(test2dHyperbolicRotation<double>(1e-14)); + CALL_SUBTEST_1(test2dHyperbolicRotation<float>(1e-5)); + CALL_SUBTEST_8(test2dHyperbolicRotation<long double>(1e-14)); + CALL_SUBTEST_2(testExponentLaws(Matrix2d(), 1e-13)); + CALL_SUBTEST_7(testExponentLaws(Matrix<double,3,3,RowMajor>(), 1e-13)); + CALL_SUBTEST_3(testExponentLaws(Matrix4cd(), 1e-13)); + CALL_SUBTEST_4(testExponentLaws(MatrixXd(8,8), 1e-13)); + CALL_SUBTEST_1(testExponentLaws(Matrix2f(), 1e-4)); + CALL_SUBTEST_5(testExponentLaws(Matrix3cf(), 1e-4)); + CALL_SUBTEST_1(testExponentLaws(Matrix4f(), 1e-4)); + CALL_SUBTEST_6(testExponentLaws(MatrixXf(8,8), 1e-4)); + CALL_SUBTEST_9(testExponentLaws(Matrix<long double,Dynamic,Dynamic>(7,7), 1e-13)); +} diff --git a/unsupported/test/matrix_square_root.cpp b/unsupported/test/matrix_square_root.cpp index 508619a7a..ea541e1ea 100644 --- a/unsupported/test/matrix_square_root.cpp +++ b/unsupported/test/matrix_square_root.cpp @@ -7,38 +7,7 @@ // 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 <unsupported/Eigen/MatrixFunctions> - -template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex> -struct generateTestMatrix; - -// for real matrices, make sure none of the eigenvalues are negative -template <typename MatrixType> -struct generateTestMatrix<MatrixType,0> -{ - static void run(MatrixType& result, typename MatrixType::Index size) - { - MatrixType mat = MatrixType::Random(size, size); - EigenSolver<MatrixType> es(mat); - typename EigenSolver<MatrixType>::EigenvalueType eivals = es.eigenvalues(); - for (typename MatrixType::Index i = 0; i < size; ++i) { - if (eivals(i).imag() == 0 && eivals(i).real() < 0) - eivals(i) = -eivals(i); - } - result = (es.eigenvectors() * eivals.asDiagonal() * es.eigenvectors().inverse()).real(); - } -}; - -// for complex matrices, any matrix is fine -template <typename MatrixType> -struct generateTestMatrix<MatrixType,1> -{ - static void run(MatrixType& result, typename MatrixType::Index size) - { - result = MatrixType::Random(size, size); - } -}; +#include "matrix_functions.h" template<typename MatrixType> void testMatrixSqrt(const MatrixType& m) |