1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2012, 2013 Chen-Pang He <jdh8@ms63.hinet.net>
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
12template<typename T>
13void test2dRotation(const T& tol)
14{
15  Matrix<T,2,2> A, B, C;
16  T angle, c, s;
17
18  A << 0, 1, -1, 0;
19  MatrixPower<Matrix<T,2,2> > Apow(A);
20
21  for (int i=0; i<=20; ++i) {
22    angle = std::pow(T(10), (i-10) / T(5.));
23    c = std::cos(angle);
24    s = std::sin(angle);
25    B << c, s, -s, c;
26
27    C = Apow(std::ldexp(angle,1) / T(EIGEN_PI));
28    std::cout << "test2dRotation: i = " << i << "   error powerm = " << relerr(C,B) << '\n';
29    VERIFY(C.isApprox(B, tol));
30  }
31}
32
33template<typename T>
34void test2dHyperbolicRotation(const T& tol)
35{
36  Matrix<std::complex<T>,2,2> A, B, C;
37  T angle, ch = std::cosh((T)1);
38  std::complex<T> ish(0, std::sinh((T)1));
39
40  A << ch, ish, -ish, ch;
41  MatrixPower<Matrix<std::complex<T>,2,2> > Apow(A);
42
43  for (int i=0; i<=20; ++i) {
44    angle = std::ldexp(static_cast<T>(i-10), -1);
45    ch = std::cosh(angle);
46    ish = std::complex<T>(0, std::sinh(angle));
47    B << ch, ish, -ish, ch;
48
49    C = Apow(angle);
50    std::cout << "test2dHyperbolicRotation: i = " << i << "   error powerm = " << relerr(C,B) << '\n';
51    VERIFY(C.isApprox(B, tol));
52  }
53}
54
55template<typename T>
56void test3dRotation(const T& tol)
57{
58  Matrix<T,3,1> v;
59  T angle;
60
61  for (int i=0; i<=20; ++i) {
62    v = Matrix<T,3,1>::Random();
63    v.normalize();
64    angle = std::pow(T(10), (i-10) / T(5.));
65    VERIFY(AngleAxis<T>(angle, v).matrix().isApprox(AngleAxis<T>(1,v).matrix().pow(angle), tol));
66  }
67}
68
69template<typename MatrixType>
70void testGeneral(const MatrixType& m, const typename MatrixType::RealScalar& tol)
71{
72  typedef typename MatrixType::RealScalar RealScalar;
73  MatrixType m1, m2, m3, m4, m5;
74  RealScalar x, y;
75
76  for (int i=0; i < g_repeat; ++i) {
77    generateTestMatrix<MatrixType>::run(m1, m.rows());
78    MatrixPower<MatrixType> mpow(m1);
79
80    x = internal::random<RealScalar>();
81    y = internal::random<RealScalar>();
82    m2 = mpow(x);
83    m3 = mpow(y);
84
85    m4 = mpow(x+y);
86    m5.noalias() = m2 * m3;
87    VERIFY(m4.isApprox(m5, tol));
88
89    m4 = mpow(x*y);
90    m5 = m2.pow(y);
91    VERIFY(m4.isApprox(m5, tol));
92
93    m4 = (std::abs(x) * m1).pow(y);
94    m5 = std::pow(std::abs(x), y) * m3;
95    VERIFY(m4.isApprox(m5, tol));
96  }
97}
98
99template<typename MatrixType>
100void testSingular(const MatrixType& m_const, const typename MatrixType::RealScalar& tol)
101{
102  // we need to pass by reference in order to prevent errors with
103  // MSVC for aligned data types ...
104  MatrixType& m = const_cast<MatrixType&>(m_const);
105
106  const int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex;
107  typedef typename internal::conditional<IsComplex, TriangularView<MatrixType,Upper>, const MatrixType&>::type TriangularType;
108  typename internal::conditional< IsComplex, ComplexSchur<MatrixType>, RealSchur<MatrixType> >::type schur;
109  MatrixType T;
110
111  for (int i=0; i < g_repeat; ++i) {
112    m.setRandom();
113    m.col(0).fill(0);
114
115    schur.compute(m);
116    T = schur.matrixT();
117    const MatrixType& U = schur.matrixU();
118    processTriangularMatrix<MatrixType>::run(m, T, U);
119    MatrixPower<MatrixType> mpow(m);
120
121    T = T.sqrt();
122    VERIFY(mpow(0.5L).isApprox(U * (TriangularType(T) * U.adjoint()), tol));
123
124    T = T.sqrt();
125    VERIFY(mpow(0.25L).isApprox(U * (TriangularType(T) * U.adjoint()), tol));
126
127    T = T.sqrt();
128    VERIFY(mpow(0.125L).isApprox(U * (TriangularType(T) * U.adjoint()), tol));
129  }
130}
131
132template<typename MatrixType>
133void testLogThenExp(const MatrixType& m_const, const typename MatrixType::RealScalar& tol)
134{
135  // we need to pass by reference in order to prevent errors with
136  // MSVC for aligned data types ...
137  MatrixType& m = const_cast<MatrixType&>(m_const);
138
139  typedef typename MatrixType::Scalar Scalar;
140  Scalar x;
141
142  for (int i=0; i < g_repeat; ++i) {
143    generateTestMatrix<MatrixType>::run(m, m.rows());
144    x = internal::random<Scalar>();
145    VERIFY(m.pow(x).isApprox((x * m.log()).exp(), tol));
146  }
147}
148
149typedef Matrix<double,3,3,RowMajor>         Matrix3dRowMajor;
150typedef Matrix<long double,3,3>             Matrix3e;
151typedef Matrix<long double,Dynamic,Dynamic> MatrixXe;
152
153void test_matrix_power()
154{
155  CALL_SUBTEST_2(test2dRotation<double>(1e-13));
156  CALL_SUBTEST_1(test2dRotation<float>(2e-5));  // was 1e-5, relaxed for clang 2.8 / linux / x86-64
157  CALL_SUBTEST_9(test2dRotation<long double>(1e-13L));
158  CALL_SUBTEST_2(test2dHyperbolicRotation<double>(1e-14));
159  CALL_SUBTEST_1(test2dHyperbolicRotation<float>(1e-5));
160  CALL_SUBTEST_9(test2dHyperbolicRotation<long double>(1e-14L));
161
162  CALL_SUBTEST_10(test3dRotation<double>(1e-13));
163  CALL_SUBTEST_11(test3dRotation<float>(1e-5));
164  CALL_SUBTEST_12(test3dRotation<long double>(1e-13L));
165
166  CALL_SUBTEST_2(testGeneral(Matrix2d(),         1e-13));
167  CALL_SUBTEST_7(testGeneral(Matrix3dRowMajor(), 1e-13));
168  CALL_SUBTEST_3(testGeneral(Matrix4cd(),        1e-13));
169  CALL_SUBTEST_4(testGeneral(MatrixXd(8,8),      2e-12));
170  CALL_SUBTEST_1(testGeneral(Matrix2f(),         1e-4));
171  CALL_SUBTEST_5(testGeneral(Matrix3cf(),        1e-4));
172  CALL_SUBTEST_8(testGeneral(Matrix4f(),         1e-4));
173  CALL_SUBTEST_6(testGeneral(MatrixXf(2,2),      1e-3)); // see bug 614
174  CALL_SUBTEST_9(testGeneral(MatrixXe(7,7),      1e-13L));
175  CALL_SUBTEST_10(testGeneral(Matrix3d(),        1e-13));
176  CALL_SUBTEST_11(testGeneral(Matrix3f(),        1e-4));
177  CALL_SUBTEST_12(testGeneral(Matrix3e(),        1e-13L));
178
179  CALL_SUBTEST_2(testSingular(Matrix2d(),         1e-13));
180  CALL_SUBTEST_7(testSingular(Matrix3dRowMajor(), 1e-13));
181  CALL_SUBTEST_3(testSingular(Matrix4cd(),        1e-13));
182  CALL_SUBTEST_4(testSingular(MatrixXd(8,8),      2e-12));
183  CALL_SUBTEST_1(testSingular(Matrix2f(),         1e-4));
184  CALL_SUBTEST_5(testSingular(Matrix3cf(),        1e-4));
185  CALL_SUBTEST_8(testSingular(Matrix4f(),         1e-4));
186  CALL_SUBTEST_6(testSingular(MatrixXf(2,2),      1e-3));
187  CALL_SUBTEST_9(testSingular(MatrixXe(7,7),      1e-13L));
188  CALL_SUBTEST_10(testSingular(Matrix3d(),        1e-13));
189  CALL_SUBTEST_11(testSingular(Matrix3f(),        1e-4));
190  CALL_SUBTEST_12(testSingular(Matrix3e(),        1e-13L));
191
192  CALL_SUBTEST_2(testLogThenExp(Matrix2d(),         1e-13));
193  CALL_SUBTEST_7(testLogThenExp(Matrix3dRowMajor(), 1e-13));
194  CALL_SUBTEST_3(testLogThenExp(Matrix4cd(),        1e-13));
195  CALL_SUBTEST_4(testLogThenExp(MatrixXd(8,8),      2e-12));
196  CALL_SUBTEST_1(testLogThenExp(Matrix2f(),         1e-4));
197  CALL_SUBTEST_5(testLogThenExp(Matrix3cf(),        1e-4));
198  CALL_SUBTEST_8(testLogThenExp(Matrix4f(),         1e-4));
199  CALL_SUBTEST_6(testLogThenExp(MatrixXf(2,2),      1e-3));
200  CALL_SUBTEST_9(testLogThenExp(MatrixXe(7,7),      1e-13L));
201  CALL_SUBTEST_10(testLogThenExp(Matrix3d(),        1e-13));
202  CALL_SUBTEST_11(testLogThenExp(Matrix3f(),        1e-4));
203  CALL_SUBTEST_12(testLogThenExp(Matrix3e(),        1e-13L));
204}
205