matrix_function.cpp (7447B)
1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2010 Jitse Niesen <jitse@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 "main.h" 11 #include <unsupported/Eigen/MatrixFunctions> 12 13 // Variant of VERIFY_IS_APPROX which uses absolute error instead of 14 // relative error. 15 #define VERIFY_IS_APPROX_ABS(a, b) VERIFY(test_isApprox_abs(a, b)) 16 17 template<typename Type1, typename Type2> 18 inline bool test_isApprox_abs(const Type1& a, const Type2& b) 19 { 20 return ((a-b).array().abs() < test_precision<typename Type1::RealScalar>()).all(); 21 } 22 23 24 // Returns a matrix with eigenvalues clustered around 0, 1 and 2. 25 template<typename MatrixType> 26 MatrixType randomMatrixWithRealEivals(const Index size) 27 { 28 typedef typename MatrixType::Scalar Scalar; 29 typedef typename MatrixType::RealScalar RealScalar; 30 MatrixType diag = MatrixType::Zero(size, size); 31 for (Index i = 0; i < size; ++i) { 32 diag(i, i) = Scalar(RealScalar(internal::random<int>(0,2))) 33 + internal::random<Scalar>() * Scalar(RealScalar(0.01)); 34 } 35 MatrixType A = MatrixType::Random(size, size); 36 HouseholderQR<MatrixType> QRofA(A); 37 return QRofA.householderQ().inverse() * diag * QRofA.householderQ(); 38 } 39 40 template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex> 41 struct randomMatrixWithImagEivals 42 { 43 // Returns a matrix with eigenvalues clustered around 0 and +/- i. 44 static MatrixType run(const Index size); 45 }; 46 47 // Partial specialization for real matrices 48 template<typename MatrixType> 49 struct randomMatrixWithImagEivals<MatrixType, 0> 50 { 51 static MatrixType run(const Index size) 52 { 53 typedef typename MatrixType::Scalar Scalar; 54 MatrixType diag = MatrixType::Zero(size, size); 55 Index i = 0; 56 while (i < size) { 57 Index randomInt = internal::random<Index>(-1, 1); 58 if (randomInt == 0 || i == size-1) { 59 diag(i, i) = internal::random<Scalar>() * Scalar(0.01); 60 ++i; 61 } else { 62 Scalar alpha = Scalar(randomInt) + internal::random<Scalar>() * Scalar(0.01); 63 diag(i, i+1) = alpha; 64 diag(i+1, i) = -alpha; 65 i += 2; 66 } 67 } 68 MatrixType A = MatrixType::Random(size, size); 69 HouseholderQR<MatrixType> QRofA(A); 70 return QRofA.householderQ().inverse() * diag * QRofA.householderQ(); 71 } 72 }; 73 74 // Partial specialization for complex matrices 75 template<typename MatrixType> 76 struct randomMatrixWithImagEivals<MatrixType, 1> 77 { 78 static MatrixType run(const Index size) 79 { 80 typedef typename MatrixType::Scalar Scalar; 81 typedef typename MatrixType::RealScalar RealScalar; 82 const Scalar imagUnit(0, 1); 83 MatrixType diag = MatrixType::Zero(size, size); 84 for (Index i = 0; i < size; ++i) { 85 diag(i, i) = Scalar(RealScalar(internal::random<Index>(-1, 1))) * imagUnit 86 + internal::random<Scalar>() * Scalar(RealScalar(0.01)); 87 } 88 MatrixType A = MatrixType::Random(size, size); 89 HouseholderQR<MatrixType> QRofA(A); 90 return QRofA.householderQ().inverse() * diag * QRofA.householderQ(); 91 } 92 }; 93 94 95 template<typename MatrixType> 96 void testMatrixExponential(const MatrixType& A) 97 { 98 typedef typename internal::traits<MatrixType>::Scalar Scalar; 99 typedef typename NumTraits<Scalar>::Real RealScalar; 100 typedef std::complex<RealScalar> ComplexScalar; 101 102 VERIFY_IS_APPROX(A.exp(), A.matrixFunction(internal::stem_function_exp<ComplexScalar>)); 103 } 104 105 template<typename MatrixType> 106 void testMatrixLogarithm(const MatrixType& A) 107 { 108 typedef typename internal::traits<MatrixType>::Scalar Scalar; 109 typedef typename NumTraits<Scalar>::Real RealScalar; 110 111 MatrixType scaledA; 112 RealScalar maxImagPartOfSpectrum = A.eigenvalues().imag().cwiseAbs().maxCoeff(); 113 if (maxImagPartOfSpectrum >= RealScalar(0.9L * EIGEN_PI)) 114 scaledA = A * RealScalar(0.9L * EIGEN_PI) / maxImagPartOfSpectrum; 115 else 116 scaledA = A; 117 118 // identity X.exp().log() = X only holds if Im(lambda) < pi for all eigenvalues of X 119 MatrixType expA = scaledA.exp(); 120 MatrixType logExpA = expA.log(); 121 VERIFY_IS_APPROX(logExpA, scaledA); 122 } 123 124 template<typename MatrixType> 125 void testHyperbolicFunctions(const MatrixType& A) 126 { 127 // Need to use absolute error because of possible cancellation when 128 // adding/subtracting expA and expmA. 129 VERIFY_IS_APPROX_ABS(A.sinh(), (A.exp() - (-A).exp()) / 2); 130 VERIFY_IS_APPROX_ABS(A.cosh(), (A.exp() + (-A).exp()) / 2); 131 } 132 133 template<typename MatrixType> 134 void testGonioFunctions(const MatrixType& A) 135 { 136 typedef typename MatrixType::Scalar Scalar; 137 typedef typename NumTraits<Scalar>::Real RealScalar; 138 typedef std::complex<RealScalar> ComplexScalar; 139 typedef Matrix<ComplexScalar, MatrixType::RowsAtCompileTime, 140 MatrixType::ColsAtCompileTime, MatrixType::Options> ComplexMatrix; 141 142 ComplexScalar imagUnit(0,1); 143 ComplexScalar two(2,0); 144 145 ComplexMatrix Ac = A.template cast<ComplexScalar>(); 146 147 ComplexMatrix exp_iA = (imagUnit * Ac).exp(); 148 ComplexMatrix exp_miA = (-imagUnit * Ac).exp(); 149 150 ComplexMatrix sinAc = A.sin().template cast<ComplexScalar>(); 151 VERIFY_IS_APPROX_ABS(sinAc, (exp_iA - exp_miA) / (two*imagUnit)); 152 153 ComplexMatrix cosAc = A.cos().template cast<ComplexScalar>(); 154 VERIFY_IS_APPROX_ABS(cosAc, (exp_iA + exp_miA) / 2); 155 } 156 157 template<typename MatrixType> 158 void testMatrix(const MatrixType& A) 159 { 160 testMatrixExponential(A); 161 testMatrixLogarithm(A); 162 testHyperbolicFunctions(A); 163 testGonioFunctions(A); 164 } 165 166 template<typename MatrixType> 167 void testMatrixType(const MatrixType& m) 168 { 169 // Matrices with clustered eigenvalue lead to different code paths 170 // in MatrixFunction.h and are thus useful for testing. 171 172 const Index size = m.rows(); 173 for (int i = 0; i < g_repeat; i++) { 174 testMatrix(MatrixType::Random(size, size).eval()); 175 testMatrix(randomMatrixWithRealEivals<MatrixType>(size)); 176 testMatrix(randomMatrixWithImagEivals<MatrixType>::run(size)); 177 } 178 } 179 180 template<typename MatrixType> 181 void testMapRef(const MatrixType& A) 182 { 183 // Test if passing Ref and Map objects is possible 184 // (Regression test for Bug #1796) 185 Index size = A.rows(); 186 MatrixType X; X.setRandom(size, size); 187 MatrixType Y(size,size); 188 Ref< MatrixType> R(Y); 189 Ref<const MatrixType> Rc(X); 190 Map< MatrixType> M(Y.data(), size, size); 191 Map<const MatrixType> Mc(X.data(), size, size); 192 193 X = X*X; // make sure sqrt is possible 194 Y = X.sqrt(); 195 R = Rc.sqrt(); 196 M = Mc.sqrt(); 197 Y = X.exp(); 198 R = Rc.exp(); 199 M = Mc.exp(); 200 X = Y; // make sure log is possible 201 Y = X.log(); 202 R = Rc.log(); 203 M = Mc.log(); 204 205 Y = X.cos() + Rc.cos() + Mc.cos(); 206 Y = X.sin() + Rc.sin() + Mc.sin(); 207 208 Y = X.cosh() + Rc.cosh() + Mc.cosh(); 209 Y = X.sinh() + Rc.sinh() + Mc.sinh(); 210 } 211 212 213 EIGEN_DECLARE_TEST(matrix_function) 214 { 215 CALL_SUBTEST_1(testMatrixType(Matrix<float,1,1>())); 216 CALL_SUBTEST_2(testMatrixType(Matrix3cf())); 217 CALL_SUBTEST_3(testMatrixType(MatrixXf(8,8))); 218 CALL_SUBTEST_4(testMatrixType(Matrix2d())); 219 CALL_SUBTEST_5(testMatrixType(Matrix<double,5,5,RowMajor>())); 220 CALL_SUBTEST_6(testMatrixType(Matrix4cd())); 221 CALL_SUBTEST_7(testMatrixType(MatrixXd(13,13))); 222 223 CALL_SUBTEST_1(testMapRef(Matrix<float,1,1>())); 224 CALL_SUBTEST_2(testMapRef(Matrix3cf())); 225 CALL_SUBTEST_3(testMapRef(MatrixXf(8,8))); 226 CALL_SUBTEST_7(testMapRef(MatrixXd(13,13))); 227 }