1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2009 Jitse Niesen <jitse@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
binom(int n,int k)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>
relerr(const MatrixBase<Derived> & A,const MatrixBase<OtherDerived> & B)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>
expfn(T x,int)28 T expfn(T x, int)
29 {
30 return std::exp(x);
31 }
32
33 template <typename T>
test2dRotation(double tol)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>
test2dHyperbolicRotation(double tol)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>
testPascal(double tol)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>
randomTest(const MatrixType & m,double tol)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
test_matrix_exponential()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