1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2009 Jitse Niesen <jitse (at) maths.leeds.ac.uk> 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 double binom(int n, int k) 13 { 14 double res = 1; 15 for (int i=0; i<k; i++) 16 res = res * (n-k+i+1) / (i+1); 17 return res; 18 } 19 20 template <typename T> 21 T expfn(T x, int) 22 { 23 return std::exp(x); 24 } 25 26 template <typename T> 27 void test2dRotation(double tol) 28 { 29 Matrix<T,2,2> A, B, C; 30 T angle; 31 32 A << 0, 1, -1, 0; 33 for (int i=0; i<=20; i++) 34 { 35 angle = static_cast<T>(pow(10, i / 5. - 2)); 36 B << std::cos(angle), std::sin(angle), -std::sin(angle), std::cos(angle); 37 38 C = (angle*A).matrixFunction(expfn); 39 std::cout << "test2dRotation: i = " << i << " error funm = " << relerr(C, B); 40 VERIFY(C.isApprox(B, static_cast<T>(tol))); 41 42 C = (angle*A).exp(); 43 std::cout << " error expm = " << relerr(C, B) << "\n"; 44 VERIFY(C.isApprox(B, static_cast<T>(tol))); 45 } 46 } 47 48 template <typename T> 49 void test2dHyperbolicRotation(double tol) 50 { 51 Matrix<std::complex<T>,2,2> A, B, C; 52 std::complex<T> imagUnit(0,1); 53 T angle, ch, sh; 54 55 for (int i=0; i<=20; i++) 56 { 57 angle = static_cast<T>((i-10) / 2.0); 58 ch = std::cosh(angle); 59 sh = std::sinh(angle); 60 A << 0, angle*imagUnit, -angle*imagUnit, 0; 61 B << ch, sh*imagUnit, -sh*imagUnit, ch; 62 63 C = A.matrixFunction(expfn); 64 std::cout << "test2dHyperbolicRotation: i = " << i << " error funm = " << relerr(C, B); 65 VERIFY(C.isApprox(B, static_cast<T>(tol))); 66 67 C = A.exp(); 68 std::cout << " error expm = " << relerr(C, B) << "\n"; 69 VERIFY(C.isApprox(B, static_cast<T>(tol))); 70 } 71 } 72 73 template <typename T> 74 void testPascal(double tol) 75 { 76 for (int size=1; size<20; size++) 77 { 78 Matrix<T,Dynamic,Dynamic> A(size,size), B(size,size), C(size,size); 79 A.setZero(); 80 for (int i=0; i<size-1; i++) 81 A(i+1,i) = static_cast<T>(i+1); 82 B.setZero(); 83 for (int i=0; i<size; i++) 84 for (int j=0; j<=i; j++) 85 B(i,j) = static_cast<T>(binom(i,j)); 86 87 C = A.matrixFunction(expfn); 88 std::cout << "testPascal: size = " << size << " error funm = " << relerr(C, B); 89 VERIFY(C.isApprox(B, static_cast<T>(tol))); 90 91 C = A.exp(); 92 std::cout << " error expm = " << relerr(C, B) << "\n"; 93 VERIFY(C.isApprox(B, static_cast<T>(tol))); 94 } 95 } 96 97 template<typename MatrixType> 98 void randomTest(const MatrixType& m, double tol) 99 { 100 /* this test covers the following files: 101 Inverse.h 102 */ 103 typename MatrixType::Index rows = m.rows(); 104 typename MatrixType::Index cols = m.cols(); 105 MatrixType m1(rows, cols), m2(rows, cols), identity = MatrixType::Identity(rows, cols); 106 107 typedef typename NumTraits<typename internal::traits<MatrixType>::Scalar>::Real RealScalar; 108 109 for(int i = 0; i < g_repeat; i++) { 110 m1 = MatrixType::Random(rows, cols); 111 112 m2 = m1.matrixFunction(expfn) * (-m1).matrixFunction(expfn); 113 std::cout << "randomTest: error funm = " << relerr(identity, m2); 114 VERIFY(identity.isApprox(m2, static_cast<RealScalar>(tol))); 115 116 m2 = m1.exp() * (-m1).exp(); 117 std::cout << " error expm = " << relerr(identity, m2) << "\n"; 118 VERIFY(identity.isApprox(m2, static_cast<RealScalar>(tol))); 119 } 120 } 121 122 void test_matrix_exponential() 123 { 124 CALL_SUBTEST_2(test2dRotation<double>(1e-13)); 125 CALL_SUBTEST_1(test2dRotation<float>(2e-5)); // was 1e-5, relaxed for clang 2.8 / linux / x86-64 126 CALL_SUBTEST_8(test2dRotation<long double>(1e-13)); 127 CALL_SUBTEST_2(test2dHyperbolicRotation<double>(1e-14)); 128 CALL_SUBTEST_1(test2dHyperbolicRotation<float>(1e-5)); 129 CALL_SUBTEST_8(test2dHyperbolicRotation<long double>(1e-14)); 130 CALL_SUBTEST_6(testPascal<float>(1e-6)); 131 CALL_SUBTEST_5(testPascal<double>(1e-15)); 132 CALL_SUBTEST_2(randomTest(Matrix2d(), 1e-13)); 133 CALL_SUBTEST_7(randomTest(Matrix<double,3,3,RowMajor>(), 1e-13)); 134 CALL_SUBTEST_3(randomTest(Matrix4cd(), 1e-13)); 135 CALL_SUBTEST_4(randomTest(MatrixXd(8,8), 1e-13)); 136 CALL_SUBTEST_1(randomTest(Matrix2f(), 1e-4)); 137 CALL_SUBTEST_5(randomTest(Matrix3cf(), 1e-4)); 138 CALL_SUBTEST_1(randomTest(Matrix4f(), 1e-4)); 139 CALL_SUBTEST_6(randomTest(MatrixXf(8,8), 1e-4)); 140 CALL_SUBTEST_9(randomTest(Matrix<long double,Dynamic,Dynamic>(7,7), 1e-13)); 141 } 142
