1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
|
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#include "matrix_functions.h"
template<typename T>
void test2dRotation(double tol)
{
Matrix<T,2,2> A, B, C;
T angle, c, s;
A << 0, 1, -1, 0;
MatrixPower<Matrix<T,2,2> > Apow(A);
for (int i=0; i<=20; ++i) {
angle = pow(10, (i-10) / 5.);
c = std::cos(angle);
s = std::sin(angle);
B << c, s, -s, c;
C = Apow(std::ldexp(angle,1) / M_PI);
std::cout << "test2dRotation: i = " << i << " error powerm = " << relerr(C,B) << '\n';
VERIFY(C.isApprox(B, static_cast<T>(tol)));
}
}
template<typename T>
void test2dHyperbolicRotation(double tol)
{
Matrix<std::complex<T>,2,2> A, B, C;
T angle, ch = std::cosh(1);
std::complex<T> ish(0, std::sinh(1));
A << ch, ish, -ish, ch;
MatrixPower<Matrix<std::complex<T>,2,2> > Apow(A);
for (int i=0; i<=20; ++i) {
angle = std::ldexp(static_cast<T>(i-10), -1);
ch = std::cosh(angle);
ish = std::complex<T>(0, std::sinh(angle));
B << ch, ish, -ish, ch;
C = Apow(angle);
std::cout << "test2dHyperbolicRotation: i = " << i << " error powerm = " << relerr(C,B) << '\n';
VERIFY(C.isApprox(B, static_cast<T>(tol)));
}
}
template<typename MatrixType>
void testExponentLaws(const MatrixType& m, double tol)
{
typedef typename MatrixType::RealScalar RealScalar;
MatrixType m1, m2, m3, m4, m5;
RealScalar x, y;
for (int i=0; i<g_repeat; ++i) {
generateTestMatrix<MatrixType>::run(m1, m.rows());
MatrixPower<MatrixType> mpow(m1);
x = internal::random<RealScalar>();
y = internal::random<RealScalar>();
m2 = mpow(x);
m3 = mpow(y);
m4 = mpow(x+y);
m5.noalias() = m2 * m3;
std::cout << "testExponentLaws: error powerm = " << relerr(m4, m5);
VERIFY(m4.isApprox(m5, static_cast<RealScalar>(tol)));
m4 = mpow(x*y);
m5 = m2.pow(y);
std::cout << " " << relerr(m4, m5);
VERIFY(m4.isApprox(m5, static_cast<RealScalar>(tol)));
m4 = (std::abs(x) * m1).pow(y);
m5 = std::pow(std::abs(x), y) * m3;
std::cout << " " << relerr(m4, m5) << '\n';
VERIFY(m4.isApprox(m5, static_cast<RealScalar>(tol)));
}
}
void test_matrix_power()
{
typedef Matrix<long double,Dynamic,Dynamic> MatrixXe;
CALL_SUBTEST_2(test2dRotation<double>(1e-13));
CALL_SUBTEST_1(test2dRotation<float>(2e-5)); // was 1e-5, relaxed for clang 2.8 / linux / x86-64
CALL_SUBTEST_9(test2dRotation<long double>(1e-13));
CALL_SUBTEST_2(test2dHyperbolicRotation<double>(1e-14));
CALL_SUBTEST_1(test2dHyperbolicRotation<float>(1e-5));
CALL_SUBTEST_9(test2dHyperbolicRotation<long double>(1e-14));
CALL_SUBTEST_2(testExponentLaws(Matrix2d(), 1e-13));
CALL_SUBTEST_7(testExponentLaws(Matrix<double,3,3,RowMajor>(), 1e-13));
CALL_SUBTEST_3(testExponentLaws(Matrix4cd(), 1e-13));
CALL_SUBTEST_4(testExponentLaws(MatrixXd(8,8), 1e-13));
CALL_SUBTEST_1(testExponentLaws(Matrix2f(), 1e-4));
CALL_SUBTEST_5(testExponentLaws(Matrix3cf(), 1e-4));
CALL_SUBTEST_8(testExponentLaws(Matrix4f(), 1e-4));
CALL_SUBTEST_6(testExponentLaws(MatrixXf(8,8), 1e-4));
}
|
