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 (at) 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 } 134
