Home | History | Annotate | Download | only in test
      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