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 17template<typename Type1, typename Type2> 18inline 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. 25template<typename MatrixType> 26MatrixType 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 41template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex> 42struct 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 49template<typename MatrixType> 50struct 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 77template<typename MatrixType> 78struct 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 98template<typename MatrixType> 99void 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 108template<typename MatrixType> 109void 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 127template<typename MatrixType> 128void 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 136template<typename MatrixType> 137void 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 160template<typename MatrixType> 161void testMatrix(const MatrixType& A) 162{ 163 testMatrixExponential(A); 164 testMatrixLogarithm(A); 165 testHyperbolicFunctions(A); 166 testGonioFunctions(A); 167} 168 169template<typename MatrixType> 170void 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 184void 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