matrix_exponential.cpp
Go to the documentation of this file.
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 "matrix_functions.h"
11 
12 double binom(int n, int k)
13 {
14  double res = 1;
15  for (int i=0; i<k; i++)
16  res = res * (n-k+i+1) / (i+1);
17  return res;
18 }
19 
20 template <typename T>
21 T expfn(T x, int)
22 {
23  return std::exp(x);
24 }
25 
26 template <typename T>
27 void test2dRotation(double tol)
28 {
29  Matrix<T,2,2> A, B, C;
30  T angle;
31 
32  A << 0, 1, -1, 0;
33  for (int i=0; i<=20; i++)
34  {
35  angle = static_cast<T>(pow(10, i / 5. - 2));
36  B << std::cos(angle), std::sin(angle), -std::sin(angle), std::cos(angle);
37 
38  C = (angle*A).matrixFunction(expfn);
39  std::cout << "test2dRotation: i = " << i << " error funm = " << relerr(C, B);
40  VERIFY(C.isApprox(B, static_cast<T>(tol)));
41 
42  C = (angle*A).exp();
43  std::cout << " error expm = " << relerr(C, B) << "\n";
44  VERIFY(C.isApprox(B, static_cast<T>(tol)));
45  }
46 }
47 
48 template <typename T>
50 {
51  Matrix<std::complex<T>,2,2> A, B, C;
52  std::complex<T> imagUnit(0,1);
53  T angle, ch, sh;
54 
55  for (int i=0; i<=20; i++)
56  {
57  angle = static_cast<T>((i-10) / 2.0);
58  ch = std::cosh(angle);
59  sh = std::sinh(angle);
60  A << 0, angle*imagUnit, -angle*imagUnit, 0;
61  B << ch, sh*imagUnit, -sh*imagUnit, ch;
62 
63  C = A.matrixFunction(expfn);
64  std::cout << "test2dHyperbolicRotation: i = " << i << " error funm = " << relerr(C, B);
65  VERIFY(C.isApprox(B, static_cast<T>(tol)));
66 
67  C = A.exp();
68  std::cout << " error expm = " << relerr(C, B) << "\n";
69  VERIFY(C.isApprox(B, static_cast<T>(tol)));
70  }
71 }
72 
73 template <typename T>
74 void testPascal(double tol)
75 {
76  for (int size=1; size<20; size++)
77  {
79  A.setZero();
80  for (int i=0; i<size-1; i++)
81  A(i+1,i) = static_cast<T>(i+1);
82  B.setZero();
83  for (int i=0; i<size; i++)
84  for (int j=0; j<=i; j++)
85  B(i,j) = static_cast<T>(binom(i,j));
86 
87  C = A.matrixFunction(expfn);
88  std::cout << "testPascal: size = " << size << " error funm = " << relerr(C, B);
89  VERIFY(C.isApprox(B, static_cast<T>(tol)));
90 
91  C = A.exp();
92  std::cout << " error expm = " << relerr(C, B) << "\n";
93  VERIFY(C.isApprox(B, static_cast<T>(tol)));
94  }
95 }
96 
97 template<typename MatrixType>
98 void randomTest(const MatrixType& m, double tol)
99 {
100  /* this test covers the following files:
101  Inverse.h
102  */
103  typename MatrixType::Index rows = m.rows();
104  typename MatrixType::Index cols = m.cols();
105  MatrixType m1(rows, cols), m2(rows, cols), identity = MatrixType::Identity(rows, cols);
106 
108 
109  for(int i = 0; i < g_repeat; i++) {
110  m1 = MatrixType::Random(rows, cols);
111 
112  m2 = m1.matrixFunction(expfn) * (-m1).matrixFunction(expfn);
113  std::cout << "randomTest: error funm = " << relerr(identity, m2);
114  VERIFY(identity.isApprox(m2, static_cast<RealScalar>(tol)));
115 
116  m2 = m1.exp() * (-m1).exp();
117  std::cout << " error expm = " << relerr(identity, m2) << "\n";
118  VERIFY(identity.isApprox(m2, static_cast<RealScalar>(tol)));
119  }
120 }
121 
123 {
124  CALL_SUBTEST_2(test2dRotation<double>(1e-13));
125  CALL_SUBTEST_1(test2dRotation<float>(2e-5)); // was 1e-5, relaxed for clang 2.8 / linux / x86-64
126  CALL_SUBTEST_8(test2dRotation<long double>(1e-13));
127  CALL_SUBTEST_2(test2dHyperbolicRotation<double>(1e-14));
128  CALL_SUBTEST_1(test2dHyperbolicRotation<float>(1e-5));
129  CALL_SUBTEST_8(test2dHyperbolicRotation<long double>(1e-14));
130  CALL_SUBTEST_6(testPascal<float>(1e-6));
131  CALL_SUBTEST_5(testPascal<double>(1e-15));
132  CALL_SUBTEST_2(randomTest(Matrix2d(), 1e-13));
133  CALL_SUBTEST_7(randomTest(Matrix<double,3,3,RowMajor>(), 1e-13));
134  CALL_SUBTEST_3(randomTest(Matrix4cd(), 1e-13));
135  CALL_SUBTEST_4(randomTest(MatrixXd(8,8), 1e-13));
136  CALL_SUBTEST_1(randomTest(Matrix2f(), 1e-4));
137  CALL_SUBTEST_5(randomTest(Matrix3cf(), 1e-4));
138  CALL_SUBTEST_1(randomTest(Matrix4f(), 1e-4));
139  CALL_SUBTEST_6(randomTest(MatrixXf(8,8), 1e-4));
140  CALL_SUBTEST_9(randomTest(Matrix<long double,Dynamic,Dynamic>(7,7), 1e-13));
141 }
Matrix3f m
EIGEN_DEVICE_FUNC const ExpReturnType exp() const
void test2dRotation(double tol)
MatrixType m2(n_dims)
int n
void testPascal(double tol)
EIGEN_DEVICE_FUNC const CoshReturnType cosh() const
double binom(int n, int k)
MatrixXf MatrixType
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:150
Matrix< SCALARA, Dynamic, Dynamic > A
Definition: bench_gemm.cpp:35
Derived::RealScalar relerr(const MatrixBase< Derived > &A, const MatrixBase< OtherDerived > &B)
Matrix< SCALARB, Dynamic, Dynamic > B
Definition: bench_gemm.cpp:36
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
EIGEN_DEVICE_FUNC const CosReturnType cos() const
void test2dHyperbolicRotation(double tol)
internal::enable_if< !(internal::is_same< typename Derived::Scalar, ScalarExponent >::value)&&EIGEN_SCALAR_BINARY_SUPPORTED(pow, typename Derived::Scalar, ScalarExponent), const EIGEN_EXPR_BINARYOP_SCALAR_RETURN_TYPE(Derived, ScalarExponent, pow) >::type pow(const Eigen::ArrayBase< Derived > &x, const ScalarExponent &exponent)
Matrix3d m1
Definition: IOFormat.cpp:2
EIGEN_DEVICE_FUNC const SinhReturnType sinh() const
static int g_repeat
Definition: main.h:144
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
Array< double, 1, 3 > e(1./3., 0.5, 2.)
void randomTest(const MatrixType &m, double tol)
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:34
Matrix< Scalar, Dynamic, Dynamic > C
Definition: bench_gemm.cpp:37
mp::number< mp::cpp_dec_float< 100 >, mp::et_on > Real
#define VERIFY(a)
Definition: main.h:325
T expfn(T x, int)
EIGEN_DEVICE_FUNC const SinReturnType sin() const
const G double tol
Definition: Group.h:83
The matrix class, also used for vectors and row-vectors.
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
void test_matrix_exponential()
std::ptrdiff_t j


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:42:48