1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
|
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#include "main.h"
#include <unsupported/Eigen/MatrixFunctions>
// Variant of VERIFY_IS_APPROX which uses absolute error instead of
// relative error.
#define VERIFY_IS_APPROX_ABS(a, b) VERIFY(test_isApprox_abs(a, b))
template<typename Type1, typename Type2>
inline bool test_isApprox_abs(const Type1& a, const Type2& b)
{
return ((a-b).array().abs() < test_precision<typename Type1::RealScalar>()).all();
}
// Returns a matrix with eigenvalues clustered around 0, 1 and 2.
template<typename MatrixType>
MatrixType randomMatrixWithRealEivals(const typename MatrixType::Index size)
{
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
MatrixType diag = MatrixType::Zero(size, size);
for (Index i = 0; i < size; ++i) {
diag(i, i) = Scalar(RealScalar(internal::random<int>(0,2)))
+ internal::random<Scalar>() * Scalar(RealScalar(0.01));
}
MatrixType A = MatrixType::Random(size, size);
HouseholderQR<MatrixType> QRofA(A);
return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
}
template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
struct randomMatrixWithImagEivals
{
// Returns a matrix with eigenvalues clustered around 0 and +/- i.
static MatrixType run(const typename MatrixType::Index size);
};
// Partial specialization for real matrices
template<typename MatrixType>
struct randomMatrixWithImagEivals<MatrixType, 0>
{
static MatrixType run(const typename MatrixType::Index size)
{
typedef typename MatrixType::Scalar Scalar;
MatrixType diag = MatrixType::Zero(size, size);
Index i = 0;
while (i < size) {
Index randomInt = internal::random<Index>(-1, 1);
if (randomInt == 0 || i == size-1) {
diag(i, i) = internal::random<Scalar>() * Scalar(0.01);
++i;
} else {
Scalar alpha = Scalar(randomInt) + internal::random<Scalar>() * Scalar(0.01);
diag(i, i+1) = alpha;
diag(i+1, i) = -alpha;
i += 2;
}
}
MatrixType A = MatrixType::Random(size, size);
HouseholderQR<MatrixType> QRofA(A);
return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
}
};
// Partial specialization for complex matrices
template<typename MatrixType>
struct randomMatrixWithImagEivals<MatrixType, 1>
{
static MatrixType run(const typename MatrixType::Index size)
{
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
const Scalar imagUnit(0, 1);
MatrixType diag = MatrixType::Zero(size, size);
for (Index i = 0; i < size; ++i) {
diag(i, i) = Scalar(RealScalar(internal::random<Index>(-1, 1))) * imagUnit
+ internal::random<Scalar>() * Scalar(RealScalar(0.01));
}
MatrixType A = MatrixType::Random(size, size);
HouseholderQR<MatrixType> QRofA(A);
return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
}
};
template<typename MatrixType>
void testMatrixExponential(const MatrixType& A)
{
typedef typename internal::traits<MatrixType>::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef std::complex<RealScalar> ComplexScalar;
VERIFY_IS_APPROX(A.exp(), A.matrixFunction(internal::stem_function_exp<ComplexScalar>));
}
template<typename MatrixType>
void testMatrixLogarithm(const MatrixType& A)
{
typedef typename internal::traits<MatrixType>::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
MatrixType scaledA;
RealScalar maxImagPartOfSpectrum = A.eigenvalues().imag().cwiseAbs().maxCoeff();
if (maxImagPartOfSpectrum >= RealScalar(0.9L * EIGEN_PI))
scaledA = A * RealScalar(0.9L * EIGEN_PI) / maxImagPartOfSpectrum;
else
scaledA = A;
// identity X.exp().log() = X only holds if Im(lambda) < pi for all eigenvalues of X
MatrixType expA = scaledA.exp();
MatrixType logExpA = expA.log();
VERIFY_IS_APPROX(logExpA, scaledA);
}
template<typename MatrixType>
void testHyperbolicFunctions(const MatrixType& A)
{
// Need to use absolute error because of possible cancellation when
// adding/subtracting expA and expmA.
VERIFY_IS_APPROX_ABS(A.sinh(), (A.exp() - (-A).exp()) / 2);
VERIFY_IS_APPROX_ABS(A.cosh(), (A.exp() + (-A).exp()) / 2);
}
template<typename MatrixType>
void testGonioFunctions(const MatrixType& A)
{
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef std::complex<RealScalar> ComplexScalar;
typedef Matrix<ComplexScalar, MatrixType::RowsAtCompileTime,
MatrixType::ColsAtCompileTime, MatrixType::Options> ComplexMatrix;
ComplexScalar imagUnit(0,1);
ComplexScalar two(2,0);
ComplexMatrix Ac = A.template cast<ComplexScalar>();
ComplexMatrix exp_iA = (imagUnit * Ac).exp();
ComplexMatrix exp_miA = (-imagUnit * Ac).exp();
ComplexMatrix sinAc = A.sin().template cast<ComplexScalar>();
VERIFY_IS_APPROX_ABS(sinAc, (exp_iA - exp_miA) / (two*imagUnit));
ComplexMatrix cosAc = A.cos().template cast<ComplexScalar>();
VERIFY_IS_APPROX_ABS(cosAc, (exp_iA + exp_miA) / 2);
}
template<typename MatrixType>
void testMatrix(const MatrixType& A)
{
testMatrixExponential(A);
testMatrixLogarithm(A);
testHyperbolicFunctions(A);
testGonioFunctions(A);
}
template<typename MatrixType>
void testMatrixType(const MatrixType& m)
{
// Matrices with clustered eigenvalue lead to different code paths
// in MatrixFunction.h and are thus useful for testing.
const Index size = m.rows();
for (int i = 0; i < g_repeat; i++) {
testMatrix(MatrixType::Random(size, size).eval());
testMatrix(randomMatrixWithRealEivals<MatrixType>(size));
testMatrix(randomMatrixWithImagEivals<MatrixType>::run(size));
}
}
void test_matrix_function()
{
CALL_SUBTEST_1(testMatrixType(Matrix<float,1,1>()));
CALL_SUBTEST_2(testMatrixType(Matrix3cf()));
CALL_SUBTEST_3(testMatrixType(MatrixXf(8,8)));
CALL_SUBTEST_4(testMatrixType(Matrix2d()));
CALL_SUBTEST_5(testMatrixType(Matrix<double,5,5,RowMajor>()));
CALL_SUBTEST_6(testMatrixType(Matrix4cd()));
CALL_SUBTEST_7(testMatrixType(MatrixXd(13,13)));
}
|
