1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2012 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
12 template <typename MatrixType, int IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex>
13 struct generateTriangularMatrix;
14
15 // for real matrices, make sure none of the eigenvalues are negative
16 template <typename MatrixType>
17 struct generateTriangularMatrix<MatrixType,0>
18 {
rungenerateTriangularMatrix19 static void run(MatrixType& result, typename MatrixType::Index size)
20 {
21 result.resize(size, size);
22 result.template triangularView<Upper>() = MatrixType::Random(size, size);
23 for (typename MatrixType::Index i = 0; i < size; ++i)
24 result.coeffRef(i,i) = std::abs(result.coeff(i,i));
25 }
26 };
27
28 // for complex matrices, any matrix is fine
29 template <typename MatrixType>
30 struct generateTriangularMatrix<MatrixType,1>
31 {
rungenerateTriangularMatrix32 static void run(MatrixType& result, typename MatrixType::Index size)
33 {
34 result.resize(size, size);
35 result.template triangularView<Upper>() = MatrixType::Random(size, size);
36 }
37 };
38
39 template<typename T>
test2dRotation(double tol)40 void test2dRotation(double tol)
41 {
42 Matrix<T,2,2> A, B, C;
43 T angle, c, s;
44
45 A << 0, 1, -1, 0;
46 MatrixPower<Matrix<T,2,2> > Apow(A);
47
48 for (int i=0; i<=20; ++i) {
49 angle = pow(10, (i-10) / 5.);
50 c = std::cos(angle);
51 s = std::sin(angle);
52 B << c, s, -s, c;
53
54 C = Apow(std::ldexp(angle,1) / M_PI);
55 std::cout << "test2dRotation: i = " << i << " error powerm = " << relerr(C,B) << '\n';
56 VERIFY(C.isApprox(B, static_cast<T>(tol)));
57 }
58 }
59
60 template<typename T>
test2dHyperbolicRotation(double tol)61 void test2dHyperbolicRotation(double tol)
62 {
63 Matrix<std::complex<T>,2,2> A, B, C;
64 T angle, ch = std::cosh((T)1);
65 std::complex<T> ish(0, std::sinh((T)1));
66
67 A << ch, ish, -ish, ch;
68 MatrixPower<Matrix<std::complex<T>,2,2> > Apow(A);
69
70 for (int i=0; i<=20; ++i) {
71 angle = std::ldexp(static_cast<T>(i-10), -1);
72 ch = std::cosh(angle);
73 ish = std::complex<T>(0, std::sinh(angle));
74 B << ch, ish, -ish, ch;
75
76 C = Apow(angle);
77 std::cout << "test2dHyperbolicRotation: i = " << i << " error powerm = " << relerr(C,B) << '\n';
78 VERIFY(C.isApprox(B, static_cast<T>(tol)));
79 }
80 }
81
82 template<typename MatrixType>
testExponentLaws(const MatrixType & m,double tol)83 void testExponentLaws(const MatrixType& m, double tol)
84 {
85 typedef typename MatrixType::RealScalar RealScalar;
86 MatrixType m1, m2, m3, m4, m5;
87 RealScalar x, y;
88
89 for (int i=0; i < g_repeat; ++i) {
90 generateTestMatrix<MatrixType>::run(m1, m.rows());
91 MatrixPower<MatrixType> mpow(m1);
92
93 x = internal::random<RealScalar>();
94 y = internal::random<RealScalar>();
95 m2 = mpow(x);
96 m3 = mpow(y);
97
98 m4 = mpow(x+y);
99 m5.noalias() = m2 * m3;
100 VERIFY(m4.isApprox(m5, static_cast<RealScalar>(tol)));
101
102 m4 = mpow(x*y);
103 m5 = m2.pow(y);
104 VERIFY(m4.isApprox(m5, static_cast<RealScalar>(tol)));
105
106 m4 = (std::abs(x) * m1).pow(y);
107 m5 = std::pow(std::abs(x), y) * m3;
108 VERIFY(m4.isApprox(m5, static_cast<RealScalar>(tol)));
109 }
110 }
111
112 typedef Matrix<double,3,3,RowMajor> Matrix3dRowMajor;
113 typedef Matrix<long double,Dynamic,Dynamic> MatrixXe;
114
test_matrix_power()115 void test_matrix_power()
116 {
117 CALL_SUBTEST_2(test2dRotation<double>(1e-13));
118 CALL_SUBTEST_1(test2dRotation<float>(2e-5)); // was 1e-5, relaxed for clang 2.8 / linux / x86-64
119 CALL_SUBTEST_9(test2dRotation<long double>(1e-13));
120 CALL_SUBTEST_2(test2dHyperbolicRotation<double>(1e-14));
121 CALL_SUBTEST_1(test2dHyperbolicRotation<float>(1e-5));
122 CALL_SUBTEST_9(test2dHyperbolicRotation<long double>(1e-14));
123
124 CALL_SUBTEST_2(testExponentLaws(Matrix2d(), 1e-13));
125 CALL_SUBTEST_7(testExponentLaws(Matrix3dRowMajor(), 1e-13));
126 CALL_SUBTEST_3(testExponentLaws(Matrix4cd(), 1e-13));
127 CALL_SUBTEST_4(testExponentLaws(MatrixXd(8,8), 2e-12));
128 CALL_SUBTEST_1(testExponentLaws(Matrix2f(), 1e-4));
129 CALL_SUBTEST_5(testExponentLaws(Matrix3cf(), 1e-4));
130 CALL_SUBTEST_8(testExponentLaws(Matrix4f(), 1e-4));
131 CALL_SUBTEST_6(testExponentLaws(MatrixXf(2,2), 1e-3)); // see bug 614
132 CALL_SUBTEST_9(testExponentLaws(MatrixXe(7,7), 1e-13));
133 }
134