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 | {
|
---|
19 | 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 | {
|
---|
32 | 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>
|
---|
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>
|
---|
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>
|
---|
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 |
|
---|
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 | }
|
---|