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
17 template<typename Type1, typename Type2>
test_isApprox_abs(const Type1 & a,const Type2 & b)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>
randomMatrixWithRealEivals(const typename MatrixType::Index size)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 {
runrandomMatrixWithImagEivals52 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 {
runrandomMatrixWithImagEivals80 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>
testMatrixExponential(const MatrixType & A)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>
testMatrixLogarithm(const MatrixType & A)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>
testHyperbolicFunctions(const MatrixType & A)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>
testGonioFunctions(const MatrixType & A)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>
testMatrix(const MatrixType & A)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>
testMatrixType(const MatrixType & m)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
test_matrix_function()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