From fcf821b77d06e7a0ac95104ef094121f8ef70dff Mon Sep 17 00:00:00 2001 From: Jitse Niesen Date: Mon, 28 Dec 2009 22:30:08 +0000 Subject: Rename test per naming convention. --- unsupported/test/matrix_exponential.cpp | 163 ++++++++++++++++++++++++++++++++ 1 file changed, 163 insertions(+) create mode 100644 unsupported/test/matrix_exponential.cpp (limited to 'unsupported/test/matrix_exponential.cpp') diff --git a/unsupported/test/matrix_exponential.cpp b/unsupported/test/matrix_exponential.cpp new file mode 100644 index 000000000..f155e5f98 --- /dev/null +++ b/unsupported/test/matrix_exponential.cpp @@ -0,0 +1,163 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Jitse Niesen +// +// 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 . + +#include "main.h" +#include + +double binom(int n, int k) +{ + double res = 1; + for (int i=0; i +double relerr(const MatrixBase& A, const MatrixBase& B) +{ + return std::sqrt((A - B).cwise().abs2().sum() / std::min(A.cwise().abs2().sum(), B.cwise().abs2().sum())); +} + +template +T expfn(T x, int) +{ + return std::exp(x); +} + +template +void test2dRotation(double tol) +{ + Matrix A, B, C; + T angle; + + A << 0, 1, -1, 0; + for (int i=0; i<=20; i++) + { + angle = static_cast(pow(10, i / 5. - 2)); + B << cos(angle), sin(angle), -sin(angle), cos(angle); + + ei_matrix_function(angle*A, expfn, &C); + std::cout << "test2dRotation: i = " << i << " error funm = " << relerr(C, B); + VERIFY(C.isApprox(B, static_cast(tol))); + + ei_matrix_exponential(angle*A, &C); + std::cout << " error expm = " << relerr(C, B) << "\n"; + VERIFY(C.isApprox(B, static_cast(tol))); + } +} + +template +void test2dHyperbolicRotation(double tol) +{ + Matrix,2,2> A, B, C; + std::complex imagUnit(0,1); + T angle, ch, sh; + + for (int i=0; i<=20; i++) + { + angle = static_cast((i-10) / 2.0); + ch = std::cosh(angle); + sh = std::sinh(angle); + A << 0, angle*imagUnit, -angle*imagUnit, 0; + B << ch, sh*imagUnit, -sh*imagUnit, ch; + + ei_matrix_function(A, expfn, &C); + std::cout << "test2dHyperbolicRotation: i = " << i << " error funm = " << relerr(C, B); + VERIFY(C.isApprox(B, static_cast(tol))); + + ei_matrix_exponential(A, &C); + std::cout << " error expm = " << relerr(C, B) << "\n"; + VERIFY(C.isApprox(B, static_cast(tol))); + } +} + +template +void testPascal(double tol) +{ + for (int size=1; size<20; size++) + { + Matrix A(size,size), B(size,size), C(size,size); + A.setZero(); + for (int i=0; i(i+1); + B.setZero(); + for (int i=0; i(binom(i,j)); + + ei_matrix_function(A, expfn, &C); + std::cout << "testPascal: size = " << size << " error funm = " << relerr(C, B); + VERIFY(C.isApprox(B, static_cast(tol))); + + ei_matrix_exponential(A, &C); + std::cout << " error expm = " << relerr(C, B) << "\n"; + VERIFY(C.isApprox(B, static_cast(tol))); + } +} + +template +void randomTest(const MatrixType& m, double tol) +{ + /* this test covers the following files: + Inverse.h + */ + int rows = m.rows(); + int cols = m.cols(); + MatrixType m1(rows, cols), m2(rows, cols), m3(rows, cols), + identity = MatrixType::Identity(rows, rows); + + typedef typename NumTraits::Scalar>::Real RealScalar; + + for(int i = 0; i < g_repeat; i++) { + m1 = MatrixType::Random(rows, cols); + + ei_matrix_function(m1, expfn, &m2); + ei_matrix_function(-m1, expfn, &m3); + std::cout << "randomTest: error funm = " << relerr(identity, m2 * m3); + VERIFY(identity.isApprox(m2 * m3, static_cast(tol))); + + ei_matrix_exponential(m1, &m2); + ei_matrix_exponential(-m1, &m3); + std::cout << " error expm = " << relerr(identity, m2 * m3) << "\n"; + VERIFY(identity.isApprox(m2 * m3, static_cast(tol))); + } +} + +void test_matrix_exponential() +{ + CALL_SUBTEST_2(test2dRotation(1e-13)); + CALL_SUBTEST_1(test2dRotation(1e-5)); + CALL_SUBTEST_2(test2dHyperbolicRotation(1e-14)); + CALL_SUBTEST_1(test2dHyperbolicRotation(1e-5)); + CALL_SUBTEST_6(testPascal(1e-6)); + CALL_SUBTEST_5(testPascal(1e-15)); + CALL_SUBTEST_2(randomTest(Matrix2d(), 1e-13)); + CALL_SUBTEST_7(randomTest(Matrix(), 1e-13)); + CALL_SUBTEST_3(randomTest(Matrix4cd(), 1e-13)); + CALL_SUBTEST_4(randomTest(MatrixXd(8,8), 1e-13)); + CALL_SUBTEST_1(randomTest(Matrix2f(), 1e-4)); + CALL_SUBTEST_5(randomTest(Matrix3cf(), 1e-4)); + CALL_SUBTEST_1(randomTest(Matrix4f(), 1e-4)); + CALL_SUBTEST_6(randomTest(MatrixXf(8,8), 1e-4)); +} -- cgit v1.2.3