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 typedef std::complex<RealScalar> ComplexScalar; 114 115 MatrixType scaledA; 116 RealScalar maxImagPartOfSpectrum = A.eigenvalues().imag().cwiseAbs().maxCoeff(); 117 if (maxImagPartOfSpectrum >= 0.9 * M_PI) 118 scaledA = A * 0.9 * M_PI / maxImagPartOfSpectrum; 119 else 120 scaledA = A; 121 122 // identity X.exp().log() = X only holds if Im(lambda) < pi for all eigenvalues of X 123 MatrixType expA = scaledA.exp(); 124 MatrixType logExpA = expA.log(); 125 VERIFY_IS_APPROX(logExpA, scaledA); 126 } 127 128 template<typename MatrixType> 129 void testHyperbolicFunctions(const MatrixType& A) 130 { 131 // Need to use absolute error because of possible cancellation when 132 // adding/subtracting expA and expmA. 133 VERIFY_IS_APPROX_ABS(A.sinh(), (A.exp() - (-A).exp()) / 2); 134 VERIFY_IS_APPROX_ABS(A.cosh(), (A.exp() + (-A).exp()) / 2); 135 } 136 137 template<typename MatrixType> 138 void testGonioFunctions(const MatrixType& A) 139 { 140 typedef typename MatrixType::Scalar Scalar; 141 typedef typename NumTraits<Scalar>::Real RealScalar; 142 typedef std::complex<RealScalar> ComplexScalar; 143 typedef Matrix<ComplexScalar, MatrixType::RowsAtCompileTime, 144 MatrixType::ColsAtCompileTime, MatrixType::Options> ComplexMatrix; 145 146 ComplexScalar imagUnit(0,1); 147 ComplexScalar two(2,0); 148 149 ComplexMatrix Ac = A.template cast<ComplexScalar>(); 150 151 ComplexMatrix exp_iA = (imagUnit * Ac).exp(); 152 ComplexMatrix exp_miA = (-imagUnit * Ac).exp(); 153 154 ComplexMatrix sinAc = A.sin().template cast<ComplexScalar>(); 155 VERIFY_IS_APPROX_ABS(sinAc, (exp_iA - exp_miA) / (two*imagUnit)); 156 157 ComplexMatrix cosAc = A.cos().template cast<ComplexScalar>(); 158 VERIFY_IS_APPROX_ABS(cosAc, (exp_iA + exp_miA) / 2); 159 } 160 161 template<typename MatrixType> 162 void testMatrix(const MatrixType& A) 163 { 164 testMatrixExponential(A); 165 testMatrixLogarithm(A); 166 testHyperbolicFunctions(A); 167 testGonioFunctions(A); 168 } 169 170 template<typename MatrixType> 171 void testMatrixType(const MatrixType& m) 172 { 173 // Matrices with clustered eigenvalue lead to different code paths 174 // in MatrixFunction.h and are thus useful for testing. 175 typedef typename MatrixType::Index Index; 176 177 const Index size = m.rows(); 178 for (int i = 0; i < g_repeat; i++) { 179 testMatrix(MatrixType::Random(size, size).eval()); 180 testMatrix(randomMatrixWithRealEivals<MatrixType>(size)); 181 testMatrix(randomMatrixWithImagEivals<MatrixType>::run(size)); 182 } 183 } 184 185 void test_matrix_function() 186 { 187 CALL_SUBTEST_1(testMatrixType(Matrix<float,1,1>())); 188 CALL_SUBTEST_2(testMatrixType(Matrix3cf())); 189 CALL_SUBTEST_3(testMatrixType(MatrixXf(8,8))); 190 CALL_SUBTEST_4(testMatrixType(Matrix2d())); 191 CALL_SUBTEST_5(testMatrixType(Matrix<double,5,5,RowMajor>())); 192 CALL_SUBTEST_6(testMatrixType(Matrix4cd())); 193 CALL_SUBTEST_7(testMatrixType(MatrixXd(13,13))); 194 } 195
