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