1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2010 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 "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 typename MatrixType::Index size) 27 { 28 typedef typename MatrixType::Index Index; 29 typedef typename MatrixType::Scalar Scalar; 30 typedef typename MatrixType::RealScalar RealScalar; 31 MatrixType diag = MatrixType::Zero(size, size); 32 for (Index i = 0; i < size; ++i) { 33 diag(i, i) = Scalar(RealScalar(internal::random<int>(0,2))) 34 + internal::random<Scalar>() * Scalar(RealScalar(0.01)); 35 } 36 MatrixType A = MatrixType::Random(size, size); 37 HouseholderQR<MatrixType> QRofA(A); 38 return QRofA.householderQ().inverse() * diag * QRofA.householderQ(); 39 } 40 41 template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex> 42 struct randomMatrixWithImagEivals 43 { 44 // Returns a matrix with eigenvalues clustered around 0 and +/- i. 45 static MatrixType run(const typename MatrixType::Index size); 46 }; 47 48 // Partial specialization for real matrices 49 template<typename MatrixType> 50 struct randomMatrixWithImagEivals<MatrixType, 0> 51 { 52 static MatrixType run(const typename MatrixType::Index size) 53 { 54 typedef typename MatrixType::Index Index; 55 typedef typename MatrixType::Scalar Scalar; 56 MatrixType diag = MatrixType::Zero(size, size); 57 Index i = 0; 58 while (i < size) { 59 Index randomInt = internal::random<Index>(-1, 1); 60 if (randomInt == 0 || i == size-1) { 61 diag(i, i) = internal::random<Scalar>() * Scalar(0.01); 62 ++i; 63 } else { 64 Scalar alpha = Scalar(randomInt) + internal::random<Scalar>() * Scalar(0.01); 65 diag(i, i+1) = alpha; 66 diag(i+1, i) = -alpha; 67 i += 2; 68 } 69 } 70 MatrixType A = MatrixType::Random(size, size); 71 HouseholderQR<MatrixType> QRofA(A); 72 return QRofA.householderQ().inverse() * diag * QRofA.householderQ(); 73 } 74 }; 75 76 // Partial specialization for complex matrices 77 template<typename MatrixType> 78 struct randomMatrixWithImagEivals<MatrixType, 1> 79 { 80 static MatrixType run(const typename MatrixType::Index size) 81 { 82 typedef typename MatrixType::Index Index; 83 typedef typename MatrixType::Scalar Scalar; 84 typedef typename MatrixType::RealScalar RealScalar; 85 const Scalar imagUnit(0, 1); 86 MatrixType diag = MatrixType::Zero(size, size); 87 for (Index i = 0; i < size; ++i) { 88 diag(i, i) = Scalar(RealScalar(internal::random<Index>(-1, 1))) * imagUnit 89 + internal::random<Scalar>() * Scalar(RealScalar(0.01)); 90 } 91 MatrixType A = MatrixType::Random(size, size); 92 HouseholderQR<MatrixType> QRofA(A); 93 return QRofA.householderQ().inverse() * diag * QRofA.householderQ(); 94 } 95 }; 96 97 98 template<typename MatrixType> 99 void testMatrixExponential(const MatrixType& A) 100 { 101 typedef typename internal::traits<MatrixType>::Scalar Scalar; 102 typedef typename NumTraits<Scalar>::Real RealScalar; 103 typedef std::complex<RealScalar> ComplexScalar; 104 105 VERIFY_IS_APPROX(A.exp(), A.matrixFunction(StdStemFunctions<ComplexScalar>::exp)); 106 } 107 108 template<typename MatrixType> 109 void testMatrixLogarithm(const MatrixType& A) 110 { 111 typedef typename internal::traits<MatrixType>::Scalar Scalar; 112 typedef typename NumTraits<Scalar>::Real RealScalar; 113 114 MatrixType scaledA; 115 RealScalar maxImagPartOfSpectrum = A.eigenvalues().imag().cwiseAbs().maxCoeff(); 116 if (maxImagPartOfSpectrum >= 0.9 * M_PI) 117 scaledA = A * 0.9 * M_PI / maxImagPartOfSpectrum; 118 else 119 scaledA = A; 120 121 // identity X.exp().log() = X only holds if Im(lambda) < pi for all eigenvalues of X 122 MatrixType expA = scaledA.exp(); 123 MatrixType logExpA = expA.log(); 124 VERIFY_IS_APPROX(logExpA, scaledA); 125 } 126 127 template<typename MatrixType> 128 void testHyperbolicFunctions(const MatrixType& A) 129 { 130 // Need to use absolute error because of possible cancellation when 131 // adding/subtracting expA and expmA. 132 VERIFY_IS_APPROX_ABS(A.sinh(), (A.exp() - (-A).exp()) / 2); 133 VERIFY_IS_APPROX_ABS(A.cosh(), (A.exp() + (-A).exp()) / 2); 134 } 135 136 template<typename MatrixType> 137 void testGonioFunctions(const MatrixType& A) 138 { 139 typedef typename MatrixType::Scalar Scalar; 140 typedef typename NumTraits<Scalar>::Real RealScalar; 141 typedef std::complex<RealScalar> ComplexScalar; 142 typedef Matrix<ComplexScalar, MatrixType::RowsAtCompileTime, 143 MatrixType::ColsAtCompileTime, MatrixType::Options> ComplexMatrix; 144 145 ComplexScalar imagUnit(0,1); 146 ComplexScalar two(2,0); 147 148 ComplexMatrix Ac = A.template cast<ComplexScalar>(); 149 150 ComplexMatrix exp_iA = (imagUnit * Ac).exp(); 151 ComplexMatrix exp_miA = (-imagUnit * Ac).exp(); 152 153 ComplexMatrix sinAc = A.sin().template cast<ComplexScalar>(); 154 VERIFY_IS_APPROX_ABS(sinAc, (exp_iA - exp_miA) / (two*imagUnit)); 155 156 ComplexMatrix cosAc = A.cos().template cast<ComplexScalar>(); 157 VERIFY_IS_APPROX_ABS(cosAc, (exp_iA + exp_miA) / 2); 158 } 159 160 template<typename MatrixType> 161 void testMatrix(const MatrixType& A) 162 { 163 testMatrixExponential(A); 164 testMatrixLogarithm(A); 165 testHyperbolicFunctions(A); 166 testGonioFunctions(A); 167 } 168 169 template<typename MatrixType> 170 void testMatrixType(const MatrixType& m) 171 { 172 // Matrices with clustered eigenvalue lead to different code paths 173 // in MatrixFunction.h and are thus useful for testing. 174 typedef typename MatrixType::Index Index; 175 176 const Index size = m.rows(); 177 for (int i = 0; i < g_repeat; i++) { 178 testMatrix(MatrixType::Random(size, size).eval()); 179 testMatrix(randomMatrixWithRealEivals<MatrixType>(size)); 180 testMatrix(randomMatrixWithImagEivals<MatrixType>::run(size)); 181 } 182 } 183 184 void test_matrix_function() 185 { 186 CALL_SUBTEST_1(testMatrixType(Matrix<float,1,1>())); 187 CALL_SUBTEST_2(testMatrixType(Matrix3cf())); 188 CALL_SUBTEST_3(testMatrixType(MatrixXf(8,8))); 189 CALL_SUBTEST_4(testMatrixType(Matrix2d())); 190 CALL_SUBTEST_5(testMatrixType(Matrix<double,5,5,RowMajor>())); 191 CALL_SUBTEST_6(testMatrixType(Matrix4cd())); 192 CALL_SUBTEST_7(testMatrixType(MatrixXd(13,13))); 193 } 194
