• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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