1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2010 Jitse Niesen <[email protected]>
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 Index size)26 MatrixType randomMatrixWithRealEivals(const Index size)
27 {
28 typedef typename MatrixType::Scalar Scalar;
29 typedef typename MatrixType::RealScalar RealScalar;
30 MatrixType diag = MatrixType::Zero(size, size);
31 for (Index i = 0; i < size; ++i) {
32 diag(i, i) = Scalar(RealScalar(internal::random<int>(0,2)))
33 + internal::random<Scalar>() * Scalar(RealScalar(0.01));
34 }
35 MatrixType A = MatrixType::Random(size, size);
36 HouseholderQR<MatrixType> QRofA(A);
37 return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
38 }
39
40 template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
41 struct randomMatrixWithImagEivals
42 {
43 // Returns a matrix with eigenvalues clustered around 0 and +/- i.
44 static MatrixType run(const Index size);
45 };
46
47 // Partial specialization for real matrices
48 template<typename MatrixType>
49 struct randomMatrixWithImagEivals<MatrixType, 0>
50 {
runrandomMatrixWithImagEivals51 static MatrixType run(const Index size)
52 {
53 typedef typename MatrixType::Scalar Scalar;
54 MatrixType diag = MatrixType::Zero(size, size);
55 Index i = 0;
56 while (i < size) {
57 Index randomInt = internal::random<Index>(-1, 1);
58 if (randomInt == 0 || i == size-1) {
59 diag(i, i) = internal::random<Scalar>() * Scalar(0.01);
60 ++i;
61 } else {
62 Scalar alpha = Scalar(randomInt) + internal::random<Scalar>() * Scalar(0.01);
63 diag(i, i+1) = alpha;
64 diag(i+1, i) = -alpha;
65 i += 2;
66 }
67 }
68 MatrixType A = MatrixType::Random(size, size);
69 HouseholderQR<MatrixType> QRofA(A);
70 return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
71 }
72 };
73
74 // Partial specialization for complex matrices
75 template<typename MatrixType>
76 struct randomMatrixWithImagEivals<MatrixType, 1>
77 {
runrandomMatrixWithImagEivals78 static MatrixType run(const Index size)
79 {
80 typedef typename MatrixType::Scalar Scalar;
81 typedef typename MatrixType::RealScalar RealScalar;
82 const Scalar imagUnit(0, 1);
83 MatrixType diag = MatrixType::Zero(size, size);
84 for (Index i = 0; i < size; ++i) {
85 diag(i, i) = Scalar(RealScalar(internal::random<Index>(-1, 1))) * imagUnit
86 + internal::random<Scalar>() * Scalar(RealScalar(0.01));
87 }
88 MatrixType A = MatrixType::Random(size, size);
89 HouseholderQR<MatrixType> QRofA(A);
90 return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
91 }
92 };
93
94
95 template<typename MatrixType>
testMatrixExponential(const MatrixType & A)96 void testMatrixExponential(const MatrixType& A)
97 {
98 typedef typename internal::traits<MatrixType>::Scalar Scalar;
99 typedef typename NumTraits<Scalar>::Real RealScalar;
100 typedef std::complex<RealScalar> ComplexScalar;
101
102 VERIFY_IS_APPROX(A.exp(), A.matrixFunction(internal::stem_function_exp<ComplexScalar>));
103 }
104
105 template<typename MatrixType>
testMatrixLogarithm(const MatrixType & A)106 void testMatrixLogarithm(const MatrixType& A)
107 {
108 typedef typename internal::traits<MatrixType>::Scalar Scalar;
109 typedef typename NumTraits<Scalar>::Real RealScalar;
110
111 MatrixType scaledA;
112 RealScalar maxImagPartOfSpectrum = A.eigenvalues().imag().cwiseAbs().maxCoeff();
113 if (maxImagPartOfSpectrum >= RealScalar(0.9L * EIGEN_PI))
114 scaledA = A * RealScalar(0.9L * EIGEN_PI) / maxImagPartOfSpectrum;
115 else
116 scaledA = A;
117
118 // identity X.exp().log() = X only holds if Im(lambda) < pi for all eigenvalues of X
119 MatrixType expA = scaledA.exp();
120 MatrixType logExpA = expA.log();
121 VERIFY_IS_APPROX(logExpA, scaledA);
122 }
123
124 template<typename MatrixType>
testHyperbolicFunctions(const MatrixType & A)125 void testHyperbolicFunctions(const MatrixType& A)
126 {
127 // Need to use absolute error because of possible cancellation when
128 // adding/subtracting expA and expmA.
129 VERIFY_IS_APPROX_ABS(A.sinh(), (A.exp() - (-A).exp()) / 2);
130 VERIFY_IS_APPROX_ABS(A.cosh(), (A.exp() + (-A).exp()) / 2);
131 }
132
133 template<typename MatrixType>
testGonioFunctions(const MatrixType & A)134 void testGonioFunctions(const MatrixType& A)
135 {
136 typedef typename MatrixType::Scalar Scalar;
137 typedef typename NumTraits<Scalar>::Real RealScalar;
138 typedef std::complex<RealScalar> ComplexScalar;
139 typedef Matrix<ComplexScalar, MatrixType::RowsAtCompileTime,
140 MatrixType::ColsAtCompileTime, MatrixType::Options> ComplexMatrix;
141
142 ComplexScalar imagUnit(0,1);
143 ComplexScalar two(2,0);
144
145 ComplexMatrix Ac = A.template cast<ComplexScalar>();
146
147 ComplexMatrix exp_iA = (imagUnit * Ac).exp();
148 ComplexMatrix exp_miA = (-imagUnit * Ac).exp();
149
150 ComplexMatrix sinAc = A.sin().template cast<ComplexScalar>();
151 VERIFY_IS_APPROX_ABS(sinAc, (exp_iA - exp_miA) / (two*imagUnit));
152
153 ComplexMatrix cosAc = A.cos().template cast<ComplexScalar>();
154 VERIFY_IS_APPROX_ABS(cosAc, (exp_iA + exp_miA) / 2);
155 }
156
157 template<typename MatrixType>
testMatrix(const MatrixType & A)158 void testMatrix(const MatrixType& A)
159 {
160 testMatrixExponential(A);
161 testMatrixLogarithm(A);
162 testHyperbolicFunctions(A);
163 testGonioFunctions(A);
164 }
165
166 template<typename MatrixType>
testMatrixType(const MatrixType & m)167 void testMatrixType(const MatrixType& m)
168 {
169 // Matrices with clustered eigenvalue lead to different code paths
170 // in MatrixFunction.h and are thus useful for testing.
171
172 const Index size = m.rows();
173 for (int i = 0; i < g_repeat; i++) {
174 testMatrix(MatrixType::Random(size, size).eval());
175 testMatrix(randomMatrixWithRealEivals<MatrixType>(size));
176 testMatrix(randomMatrixWithImagEivals<MatrixType>::run(size));
177 }
178 }
179
180 template<typename MatrixType>
testMapRef(const MatrixType & A)181 void testMapRef(const MatrixType& A)
182 {
183 // Test if passing Ref and Map objects is possible
184 // (Regression test for Bug #1796)
185 Index size = A.rows();
186 MatrixType X; X.setRandom(size, size);
187 MatrixType Y(size,size);
188 Ref< MatrixType> R(Y);
189 Ref<const MatrixType> Rc(X);
190 Map< MatrixType> M(Y.data(), size, size);
191 Map<const MatrixType> Mc(X.data(), size, size);
192
193 X = X*X; // make sure sqrt is possible
194 Y = X.sqrt();
195 R = Rc.sqrt();
196 M = Mc.sqrt();
197 Y = X.exp();
198 R = Rc.exp();
199 M = Mc.exp();
200 X = Y; // make sure log is possible
201 Y = X.log();
202 R = Rc.log();
203 M = Mc.log();
204
205 Y = X.cos() + Rc.cos() + Mc.cos();
206 Y = X.sin() + Rc.sin() + Mc.sin();
207
208 Y = X.cosh() + Rc.cosh() + Mc.cosh();
209 Y = X.sinh() + Rc.sinh() + Mc.sinh();
210 }
211
212
EIGEN_DECLARE_TEST(matrix_function)213 EIGEN_DECLARE_TEST(matrix_function)
214 {
215 CALL_SUBTEST_1(testMatrixType(Matrix<float,1,1>()));
216 CALL_SUBTEST_2(testMatrixType(Matrix3cf()));
217 CALL_SUBTEST_3(testMatrixType(MatrixXf(8,8)));
218 CALL_SUBTEST_4(testMatrixType(Matrix2d()));
219 CALL_SUBTEST_5(testMatrixType(Matrix<double,5,5,RowMajor>()));
220 CALL_SUBTEST_6(testMatrixType(Matrix4cd()));
221 CALL_SUBTEST_7(testMatrixType(MatrixXd(13,13)));
222
223 CALL_SUBTEST_1(testMapRef(Matrix<float,1,1>()));
224 CALL_SUBTEST_2(testMapRef(Matrix3cf()));
225 CALL_SUBTEST_3(testMapRef(MatrixXf(8,8)));
226 CALL_SUBTEST_7(testMapRef(MatrixXd(13,13)));
227 }
228