1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2012, 2013 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 T> 13 void test2dRotation(const T& tol) 14 { 15 Matrix<T,2,2> A, B, C; 16 T angle, c, s; 17 18 A << 0, 1, -1, 0; 19 MatrixPower<Matrix<T,2,2> > Apow(A); 20 21 for (int i=0; i<=20; ++i) { 22 angle = std::pow(T(10), (i-10) / T(5.)); 23 c = std::cos(angle); 24 s = std::sin(angle); 25 B << c, s, -s, c; 26 27 C = Apow(std::ldexp(angle,1) / T(EIGEN_PI)); 28 std::cout << "test2dRotation: i = " << i << " error powerm = " << relerr(C,B) << '\n'; 29 VERIFY(C.isApprox(B, tol)); 30 } 31 } 32 33 template<typename T> 34 void test2dHyperbolicRotation(const T& tol) 35 { 36 Matrix<std::complex<T>,2,2> A, B, C; 37 T angle, ch = std::cosh((T)1); 38 std::complex<T> ish(0, std::sinh((T)1)); 39 40 A << ch, ish, -ish, ch; 41 MatrixPower<Matrix<std::complex<T>,2,2> > Apow(A); 42 43 for (int i=0; i<=20; ++i) { 44 angle = std::ldexp(static_cast<T>(i-10), -1); 45 ch = std::cosh(angle); 46 ish = std::complex<T>(0, std::sinh(angle)); 47 B << ch, ish, -ish, ch; 48 49 C = Apow(angle); 50 std::cout << "test2dHyperbolicRotation: i = " << i << " error powerm = " << relerr(C,B) << '\n'; 51 VERIFY(C.isApprox(B, tol)); 52 } 53 } 54 55 template<typename T> 56 void test3dRotation(const T& tol) 57 { 58 Matrix<T,3,1> v; 59 T angle; 60 61 for (int i=0; i<=20; ++i) { 62 v = Matrix<T,3,1>::Random(); 63 v.normalize(); 64 angle = std::pow(T(10), (i-10) / T(5.)); 65 VERIFY(AngleAxis<T>(angle, v).matrix().isApprox(AngleAxis<T>(1,v).matrix().pow(angle), tol)); 66 } 67 } 68 69 template<typename MatrixType> 70 void testGeneral(const MatrixType& m, const typename MatrixType::RealScalar& tol) 71 { 72 typedef typename MatrixType::RealScalar RealScalar; 73 MatrixType m1, m2, m3, m4, m5; 74 RealScalar x, y; 75 76 for (int i=0; i < g_repeat; ++i) { 77 generateTestMatrix<MatrixType>::run(m1, m.rows()); 78 MatrixPower<MatrixType> mpow(m1); 79 80 x = internal::random<RealScalar>(); 81 y = internal::random<RealScalar>(); 82 m2 = mpow(x); 83 m3 = mpow(y); 84 85 m4 = mpow(x+y); 86 m5.noalias() = m2 * m3; 87 VERIFY(m4.isApprox(m5, tol)); 88 89 m4 = mpow(x*y); 90 m5 = m2.pow(y); 91 VERIFY(m4.isApprox(m5, tol)); 92 93 m4 = (std::abs(x) * m1).pow(y); 94 m5 = std::pow(std::abs(x), y) * m3; 95 VERIFY(m4.isApprox(m5, tol)); 96 } 97 } 98 99 template<typename MatrixType> 100 void testSingular(const MatrixType& m_const, const typename MatrixType::RealScalar& tol) 101 { 102 // we need to pass by reference in order to prevent errors with 103 // MSVC for aligned data types ... 104 MatrixType& m = const_cast<MatrixType&>(m_const); 105 106 const int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex; 107 typedef typename internal::conditional<IsComplex, TriangularView<MatrixType,Upper>, const MatrixType&>::type TriangularType; 108 typename internal::conditional< IsComplex, ComplexSchur<MatrixType>, RealSchur<MatrixType> >::type schur; 109 MatrixType T; 110 111 for (int i=0; i < g_repeat; ++i) { 112 m.setRandom(); 113 m.col(0).fill(0); 114 115 schur.compute(m); 116 T = schur.matrixT(); 117 const MatrixType& U = schur.matrixU(); 118 processTriangularMatrix<MatrixType>::run(m, T, U); 119 MatrixPower<MatrixType> mpow(m); 120 121 T = T.sqrt(); 122 VERIFY(mpow(0.5L).isApprox(U * (TriangularType(T) * U.adjoint()), tol)); 123 124 T = T.sqrt(); 125 VERIFY(mpow(0.25L).isApprox(U * (TriangularType(T) * U.adjoint()), tol)); 126 127 T = T.sqrt(); 128 VERIFY(mpow(0.125L).isApprox(U * (TriangularType(T) * U.adjoint()), tol)); 129 } 130 } 131 132 template<typename MatrixType> 133 void testLogThenExp(const MatrixType& m_const, const typename MatrixType::RealScalar& tol) 134 { 135 // we need to pass by reference in order to prevent errors with 136 // MSVC for aligned data types ... 137 MatrixType& m = const_cast<MatrixType&>(m_const); 138 139 typedef typename MatrixType::Scalar Scalar; 140 Scalar x; 141 142 for (int i=0; i < g_repeat; ++i) { 143 generateTestMatrix<MatrixType>::run(m, m.rows()); 144 x = internal::random<Scalar>(); 145 VERIFY(m.pow(x).isApprox((x * m.log()).exp(), tol)); 146 } 147 } 148 149 typedef Matrix<double,3,3,RowMajor> Matrix3dRowMajor; 150 typedef Matrix<long double,3,3> Matrix3e; 151 typedef Matrix<long double,Dynamic,Dynamic> MatrixXe; 152 153 void test_matrix_power() 154 { 155 CALL_SUBTEST_2(test2dRotation<double>(1e-13)); 156 CALL_SUBTEST_1(test2dRotation<float>(2e-5)); // was 1e-5, relaxed for clang 2.8 / linux / x86-64 157 CALL_SUBTEST_9(test2dRotation<long double>(1e-13L)); 158 CALL_SUBTEST_2(test2dHyperbolicRotation<double>(1e-14)); 159 CALL_SUBTEST_1(test2dHyperbolicRotation<float>(1e-5)); 160 CALL_SUBTEST_9(test2dHyperbolicRotation<long double>(1e-14L)); 161 162 CALL_SUBTEST_10(test3dRotation<double>(1e-13)); 163 CALL_SUBTEST_11(test3dRotation<float>(1e-5)); 164 CALL_SUBTEST_12(test3dRotation<long double>(1e-13L)); 165 166 CALL_SUBTEST_2(testGeneral(Matrix2d(), 1e-13)); 167 CALL_SUBTEST_7(testGeneral(Matrix3dRowMajor(), 1e-13)); 168 CALL_SUBTEST_3(testGeneral(Matrix4cd(), 1e-13)); 169 CALL_SUBTEST_4(testGeneral(MatrixXd(8,8), 2e-12)); 170 CALL_SUBTEST_1(testGeneral(Matrix2f(), 1e-4)); 171 CALL_SUBTEST_5(testGeneral(Matrix3cf(), 1e-4)); 172 CALL_SUBTEST_8(testGeneral(Matrix4f(), 1e-4)); 173 CALL_SUBTEST_6(testGeneral(MatrixXf(2,2), 1e-3)); // see bug 614 174 CALL_SUBTEST_9(testGeneral(MatrixXe(7,7), 1e-13L)); 175 CALL_SUBTEST_10(testGeneral(Matrix3d(), 1e-13)); 176 CALL_SUBTEST_11(testGeneral(Matrix3f(), 1e-4)); 177 CALL_SUBTEST_12(testGeneral(Matrix3e(), 1e-13L)); 178 179 CALL_SUBTEST_2(testSingular(Matrix2d(), 1e-13)); 180 CALL_SUBTEST_7(testSingular(Matrix3dRowMajor(), 1e-13)); 181 CALL_SUBTEST_3(testSingular(Matrix4cd(), 1e-13)); 182 CALL_SUBTEST_4(testSingular(MatrixXd(8,8), 2e-12)); 183 CALL_SUBTEST_1(testSingular(Matrix2f(), 1e-4)); 184 CALL_SUBTEST_5(testSingular(Matrix3cf(), 1e-4)); 185 CALL_SUBTEST_8(testSingular(Matrix4f(), 1e-4)); 186 CALL_SUBTEST_6(testSingular(MatrixXf(2,2), 1e-3)); 187 CALL_SUBTEST_9(testSingular(MatrixXe(7,7), 1e-13L)); 188 CALL_SUBTEST_10(testSingular(Matrix3d(), 1e-13)); 189 CALL_SUBTEST_11(testSingular(Matrix3f(), 1e-4)); 190 CALL_SUBTEST_12(testSingular(Matrix3e(), 1e-13L)); 191 192 CALL_SUBTEST_2(testLogThenExp(Matrix2d(), 1e-13)); 193 CALL_SUBTEST_7(testLogThenExp(Matrix3dRowMajor(), 1e-13)); 194 CALL_SUBTEST_3(testLogThenExp(Matrix4cd(), 1e-13)); 195 CALL_SUBTEST_4(testLogThenExp(MatrixXd(8,8), 2e-12)); 196 CALL_SUBTEST_1(testLogThenExp(Matrix2f(), 1e-4)); 197 CALL_SUBTEST_5(testLogThenExp(Matrix3cf(), 1e-4)); 198 CALL_SUBTEST_8(testLogThenExp(Matrix4f(), 1e-4)); 199 CALL_SUBTEST_6(testLogThenExp(MatrixXf(2,2), 1e-3)); 200 CALL_SUBTEST_9(testLogThenExp(MatrixXe(7,7), 1e-13L)); 201 CALL_SUBTEST_10(testLogThenExp(Matrix3d(), 1e-13)); 202 CALL_SUBTEST_11(testLogThenExp(Matrix3f(), 1e-4)); 203 CALL_SUBTEST_12(testLogThenExp(Matrix3e(), 1e-13L)); 204 } 205
