Home | History | Annotate | Download | only in test
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2009 Jitse Niesen <jitse (at) maths.leeds.ac.uk>
      5 //
      6 // This Source Code Form is subject to the terms of the Mozilla
      7 // Public License v. 2.0. If a copy of the MPL was not distributed
      8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
      9 
     10 #include "main.h"
     11 #include <unsupported/Eigen/MatrixFunctions>
     12 
     13 double binom(int n, int k)
     14 {
     15   double res = 1;
     16   for (int i=0; i<k; i++)
     17     res = res * (n-k+i+1) / (i+1);
     18   return res;
     19 }
     20 
     21 template <typename Derived, typename OtherDerived>
     22 double relerr(const MatrixBase<Derived>& A, const MatrixBase<OtherDerived>& B)
     23 {
     24   return std::sqrt((A - B).cwiseAbs2().sum() / (std::min)(A.cwiseAbs2().sum(), B.cwiseAbs2().sum()));
     25 }
     26 
     27 template <typename T>
     28 T expfn(T x, int)
     29 {
     30   return std::exp(x);
     31 }
     32 
     33 template <typename T>
     34 void test2dRotation(double tol)
     35 {
     36   Matrix<T,2,2> A, B, C;
     37   T angle;
     38 
     39   A << 0, 1, -1, 0;
     40   for (int i=0; i<=20; i++)
     41   {
     42     angle = static_cast<T>(pow(10, i / 5. - 2));
     43     B << std::cos(angle), std::sin(angle), -std::sin(angle), std::cos(angle);
     44 
     45     C = (angle*A).matrixFunction(expfn);
     46     std::cout << "test2dRotation: i = " << i << "   error funm = " << relerr(C, B);
     47     VERIFY(C.isApprox(B, static_cast<T>(tol)));
     48 
     49     C = (angle*A).exp();
     50     std::cout << "   error expm = " << relerr(C, B) << "\n";
     51     VERIFY(C.isApprox(B, static_cast<T>(tol)));
     52   }
     53 }
     54 
     55 template <typename T>
     56 void test2dHyperbolicRotation(double tol)
     57 {
     58   Matrix<std::complex<T>,2,2> A, B, C;
     59   std::complex<T> imagUnit(0,1);
     60   T angle, ch, sh;
     61 
     62   for (int i=0; i<=20; i++)
     63   {
     64     angle = static_cast<T>((i-10) / 2.0);
     65     ch = std::cosh(angle);
     66     sh = std::sinh(angle);
     67     A << 0, angle*imagUnit, -angle*imagUnit, 0;
     68     B << ch, sh*imagUnit, -sh*imagUnit, ch;
     69 
     70     C = A.matrixFunction(expfn);
     71     std::cout << "test2dHyperbolicRotation: i = " << i << "   error funm = " << relerr(C, B);
     72     VERIFY(C.isApprox(B, static_cast<T>(tol)));
     73 
     74     C = A.exp();
     75     std::cout << "   error expm = " << relerr(C, B) << "\n";
     76     VERIFY(C.isApprox(B, static_cast<T>(tol)));
     77   }
     78 }
     79 
     80 template <typename T>
     81 void testPascal(double tol)
     82 {
     83   for (int size=1; size<20; size++)
     84   {
     85     Matrix<T,Dynamic,Dynamic> A(size,size), B(size,size), C(size,size);
     86     A.setZero();
     87     for (int i=0; i<size-1; i++)
     88       A(i+1,i) = static_cast<T>(i+1);
     89     B.setZero();
     90     for (int i=0; i<size; i++)
     91       for (int j=0; j<=i; j++)
     92     B(i,j) = static_cast<T>(binom(i,j));
     93 
     94     C = A.matrixFunction(expfn);
     95     std::cout << "testPascal: size = " << size << "   error funm = " << relerr(C, B);
     96     VERIFY(C.isApprox(B, static_cast<T>(tol)));
     97 
     98     C = A.exp();
     99     std::cout << "   error expm = " << relerr(C, B) << "\n";
    100     VERIFY(C.isApprox(B, static_cast<T>(tol)));
    101   }
    102 }
    103 
    104 template<typename MatrixType>
    105 void randomTest(const MatrixType& m, double tol)
    106 {
    107   /* this test covers the following files:
    108      Inverse.h
    109   */
    110   typename MatrixType::Index rows = m.rows();
    111   typename MatrixType::Index cols = m.cols();
    112   MatrixType m1(rows, cols), m2(rows, cols), m3(rows, cols),
    113              identity = MatrixType::Identity(rows, rows);
    114 
    115   typedef typename NumTraits<typename internal::traits<MatrixType>::Scalar>::Real RealScalar;
    116 
    117   for(int i = 0; i < g_repeat; i++) {
    118     m1 = MatrixType::Random(rows, cols);
    119 
    120     m2 = m1.matrixFunction(expfn) * (-m1).matrixFunction(expfn);
    121     std::cout << "randomTest: error funm = " << relerr(identity, m2);
    122     VERIFY(identity.isApprox(m2, static_cast<RealScalar>(tol)));
    123 
    124     m2 = m1.exp() * (-m1).exp();
    125     std::cout << "   error expm = " << relerr(identity, m2) << "\n";
    126     VERIFY(identity.isApprox(m2, static_cast<RealScalar>(tol)));
    127   }
    128 }
    129 
    130 void test_matrix_exponential()
    131 {
    132   CALL_SUBTEST_2(test2dRotation<double>(1e-13));
    133   CALL_SUBTEST_1(test2dRotation<float>(2e-5));  // was 1e-5, relaxed for clang 2.8 / linux / x86-64
    134   CALL_SUBTEST_8(test2dRotation<long double>(1e-13));
    135   CALL_SUBTEST_2(test2dHyperbolicRotation<double>(1e-14));
    136   CALL_SUBTEST_1(test2dHyperbolicRotation<float>(1e-5));
    137   CALL_SUBTEST_8(test2dHyperbolicRotation<long double>(1e-14));
    138   CALL_SUBTEST_6(testPascal<float>(1e-6));
    139   CALL_SUBTEST_5(testPascal<double>(1e-15));
    140   CALL_SUBTEST_2(randomTest(Matrix2d(), 1e-13));
    141   CALL_SUBTEST_7(randomTest(Matrix<double,3,3,RowMajor>(), 1e-13));
    142   CALL_SUBTEST_3(randomTest(Matrix4cd(), 1e-13));
    143   CALL_SUBTEST_4(randomTest(MatrixXd(8,8), 1e-13));
    144   CALL_SUBTEST_1(randomTest(Matrix2f(), 1e-4));
    145   CALL_SUBTEST_5(randomTest(Matrix3cf(), 1e-4));
    146   CALL_SUBTEST_1(randomTest(Matrix4f(), 1e-4));
    147   CALL_SUBTEST_6(randomTest(MatrixXf(8,8), 1e-4));
    148   CALL_SUBTEST_9(randomTest(Matrix<long double,Dynamic,Dynamic>(7,7), 1e-13));
    149 }
    150