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
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>
testHyperbolicFunctions(const MatrixType & A)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>
testGonioFunctions(const MatrixType & A)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>
testMatrix(const MatrixType & A)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>
testMatrixType(const MatrixType & m)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
test_matrix_function()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