1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008 Gael Guennebaud <[email protected]>
5 // Copyright (C) 2010,2012 Jitse Niesen <[email protected]>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11 #include "main.h"
12 #include <limits>
13 #include <Eigen/Eigenvalues>
14
15 template<typename EigType,typename MatType>
check_eigensolver_for_given_mat(const EigType & eig,const MatType & a)16 void check_eigensolver_for_given_mat(const EigType &eig, const MatType& a)
17 {
18 typedef typename NumTraits<typename MatType::Scalar>::Real RealScalar;
19 typedef Matrix<RealScalar, MatType::RowsAtCompileTime, 1> RealVectorType;
20 typedef typename std::complex<RealScalar> Complex;
21 Index n = a.rows();
22 VERIFY_IS_EQUAL(eig.info(), Success);
23 VERIFY_IS_APPROX(a * eig.pseudoEigenvectors(), eig.pseudoEigenvectors() * eig.pseudoEigenvalueMatrix());
24 VERIFY_IS_APPROX(a.template cast<Complex>() * eig.eigenvectors(),
25 eig.eigenvectors() * eig.eigenvalues().asDiagonal());
26 VERIFY_IS_APPROX(eig.eigenvectors().colwise().norm(), RealVectorType::Ones(n).transpose());
27 VERIFY_IS_APPROX(a.eigenvalues(), eig.eigenvalues());
28 }
29
eigensolver(const MatrixType & m)30 template<typename MatrixType> void eigensolver(const MatrixType& m)
31 {
32 /* this test covers the following files:
33 EigenSolver.h
34 */
35 Index rows = m.rows();
36 Index cols = m.cols();
37
38 typedef typename MatrixType::Scalar Scalar;
39 typedef typename NumTraits<Scalar>::Real RealScalar;
40 typedef typename std::complex<RealScalar> Complex;
41
42 MatrixType a = MatrixType::Random(rows,cols);
43 MatrixType a1 = MatrixType::Random(rows,cols);
44 MatrixType symmA = a.adjoint() * a + a1.adjoint() * a1;
45
46 EigenSolver<MatrixType> ei0(symmA);
47 VERIFY_IS_EQUAL(ei0.info(), Success);
48 VERIFY_IS_APPROX(symmA * ei0.pseudoEigenvectors(), ei0.pseudoEigenvectors() * ei0.pseudoEigenvalueMatrix());
49 VERIFY_IS_APPROX((symmA.template cast<Complex>()) * (ei0.pseudoEigenvectors().template cast<Complex>()),
50 (ei0.pseudoEigenvectors().template cast<Complex>()) * (ei0.eigenvalues().asDiagonal()));
51
52 EigenSolver<MatrixType> ei1(a);
53 CALL_SUBTEST( check_eigensolver_for_given_mat(ei1,a) );
54
55 EigenSolver<MatrixType> ei2;
56 ei2.setMaxIterations(RealSchur<MatrixType>::m_maxIterationsPerRow * rows).compute(a);
57 VERIFY_IS_EQUAL(ei2.info(), Success);
58 VERIFY_IS_EQUAL(ei2.eigenvectors(), ei1.eigenvectors());
59 VERIFY_IS_EQUAL(ei2.eigenvalues(), ei1.eigenvalues());
60 if (rows > 2) {
61 ei2.setMaxIterations(1).compute(a);
62 VERIFY_IS_EQUAL(ei2.info(), NoConvergence);
63 VERIFY_IS_EQUAL(ei2.getMaxIterations(), 1);
64 }
65
66 EigenSolver<MatrixType> eiNoEivecs(a, false);
67 VERIFY_IS_EQUAL(eiNoEivecs.info(), Success);
68 VERIFY_IS_APPROX(ei1.eigenvalues(), eiNoEivecs.eigenvalues());
69 VERIFY_IS_APPROX(ei1.pseudoEigenvalueMatrix(), eiNoEivecs.pseudoEigenvalueMatrix());
70
71 MatrixType id = MatrixType::Identity(rows, cols);
72 VERIFY_IS_APPROX(id.operatorNorm(), RealScalar(1));
73
74 if (rows > 2 && rows < 20)
75 {
76 // Test matrix with NaN
77 a(0,0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN();
78 EigenSolver<MatrixType> eiNaN(a);
79 VERIFY_IS_NOT_EQUAL(eiNaN.info(), Success);
80 }
81
82 // regression test for bug 1098
83 {
84 EigenSolver<MatrixType> eig(a.adjoint() * a);
85 eig.compute(a.adjoint() * a);
86 }
87
88 // regression test for bug 478
89 {
90 a.setZero();
91 EigenSolver<MatrixType> ei3(a);
92 VERIFY_IS_EQUAL(ei3.info(), Success);
93 VERIFY_IS_MUCH_SMALLER_THAN(ei3.eigenvalues().norm(),RealScalar(1));
94 VERIFY((ei3.eigenvectors().transpose()*ei3.eigenvectors().transpose()).eval().isIdentity());
95 }
96 }
97
eigensolver_verify_assert(const MatrixType & m)98 template<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m)
99 {
100 EigenSolver<MatrixType> eig;
101 VERIFY_RAISES_ASSERT(eig.eigenvectors());
102 VERIFY_RAISES_ASSERT(eig.pseudoEigenvectors());
103 VERIFY_RAISES_ASSERT(eig.pseudoEigenvalueMatrix());
104 VERIFY_RAISES_ASSERT(eig.eigenvalues());
105
106 MatrixType a = MatrixType::Random(m.rows(),m.cols());
107 eig.compute(a, false);
108 VERIFY_RAISES_ASSERT(eig.eigenvectors());
109 VERIFY_RAISES_ASSERT(eig.pseudoEigenvectors());
110 }
111
112
113 template<typename CoeffType>
114 Matrix<typename CoeffType::Scalar,Dynamic,Dynamic>
make_companion(const CoeffType & coeffs)115 make_companion(const CoeffType& coeffs)
116 {
117 Index n = coeffs.size()-1;
118 Matrix<typename CoeffType::Scalar,Dynamic,Dynamic> res(n,n);
119 res.setZero();
120 res.row(0) = -coeffs.tail(n) / coeffs(0);
121 res.diagonal(-1).setOnes();
122 return res;
123 }
124
125 template<int>
eigensolver_generic_extra()126 void eigensolver_generic_extra()
127 {
128 {
129 // regression test for bug 793
130 MatrixXd a(3,3);
131 a << 0, 0, 1,
132 1, 1, 1,
133 1, 1e+200, 1;
134 Eigen::EigenSolver<MatrixXd> eig(a);
135 double scale = 1e-200; // scale to avoid overflow during the comparisons
136 VERIFY_IS_APPROX(a * eig.pseudoEigenvectors()*scale, eig.pseudoEigenvectors() * eig.pseudoEigenvalueMatrix()*scale);
137 VERIFY_IS_APPROX(a * eig.eigenvectors()*scale, eig.eigenvectors() * eig.eigenvalues().asDiagonal()*scale);
138 }
139 {
140 // check a case where all eigenvalues are null.
141 MatrixXd a(2,2);
142 a << 1, 1,
143 -1, -1;
144 Eigen::EigenSolver<MatrixXd> eig(a);
145 VERIFY_IS_APPROX(eig.pseudoEigenvectors().squaredNorm(), 2.);
146 VERIFY_IS_APPROX((a * eig.pseudoEigenvectors()).norm()+1., 1.);
147 VERIFY_IS_APPROX((eig.pseudoEigenvectors() * eig.pseudoEigenvalueMatrix()).norm()+1., 1.);
148 VERIFY_IS_APPROX((a * eig.eigenvectors()).norm()+1., 1.);
149 VERIFY_IS_APPROX((eig.eigenvectors() * eig.eigenvalues().asDiagonal()).norm()+1., 1.);
150 }
151
152 // regression test for bug 933
153 {
154 {
155 VectorXd coeffs(5); coeffs << 1, -3, -175, -225, 2250;
156 MatrixXd C = make_companion(coeffs);
157 EigenSolver<MatrixXd> eig(C);
158 CALL_SUBTEST( check_eigensolver_for_given_mat(eig,C) );
159 }
160 {
161 // this test is tricky because it requires high accuracy in smallest eigenvalues
162 VectorXd coeffs(5); coeffs << 6.154671e-15, -1.003870e-10, -9.819570e-01, 3.995715e+03, 2.211511e+08;
163 MatrixXd C = make_companion(coeffs);
164 EigenSolver<MatrixXd> eig(C);
165 CALL_SUBTEST( check_eigensolver_for_given_mat(eig,C) );
166 Index n = C.rows();
167 for(Index i=0;i<n;++i)
168 {
169 typedef std::complex<double> Complex;
170 MatrixXcd ac = C.cast<Complex>();
171 ac.diagonal().array() -= eig.eigenvalues()(i);
172 VectorXd sv = ac.jacobiSvd().singularValues();
173 // comparing to sv(0) is not enough here to catch the "bug",
174 // the hard-coded 1.0 is important!
175 VERIFY_IS_MUCH_SMALLER_THAN(sv(n-1), 1.0);
176 }
177 }
178 }
179 // regression test for bug 1557
180 {
181 // this test is interesting because it contains zeros on the diagonal.
182 MatrixXd A_bug1557(3,3);
183 A_bug1557 << 0, 0, 0, 1, 0, 0.5887907064808635127, 0, 1, 0;
184 EigenSolver<MatrixXd> eig(A_bug1557);
185 CALL_SUBTEST( check_eigensolver_for_given_mat(eig,A_bug1557) );
186 }
187
188 // regression test for bug 1174
189 {
190 Index n = 12;
191 MatrixXf A_bug1174(n,n);
192 A_bug1174 << 262144, 0, 0, 262144, 786432, 0, 0, 0, 0, 0, 0, 786432,
193 262144, 0, 0, 262144, 786432, 0, 0, 0, 0, 0, 0, 786432,
194 262144, 0, 0, 262144, 786432, 0, 0, 0, 0, 0, 0, 786432,
195 262144, 0, 0, 262144, 786432, 0, 0, 0, 0, 0, 0, 786432,
196 0, 262144, 262144, 0, 0, 262144, 262144, 262144, 262144, 262144, 262144, 0,
197 0, 262144, 262144, 0, 0, 262144, 262144, 262144, 262144, 262144, 262144, 0,
198 0, 262144, 262144, 0, 0, 262144, 262144, 262144, 262144, 262144, 262144, 0,
199 0, 262144, 262144, 0, 0, 262144, 262144, 262144, 262144, 262144, 262144, 0,
200 0, 262144, 262144, 0, 0, 262144, 262144, 262144, 262144, 262144, 262144, 0,
201 0, 262144, 262144, 0, 0, 262144, 262144, 262144, 262144, 262144, 262144, 0,
202 0, 262144, 262144, 0, 0, 262144, 262144, 262144, 262144, 262144, 262144, 0,
203 0, 262144, 262144, 0, 0, 262144, 262144, 262144, 262144, 262144, 262144, 0;
204 EigenSolver<MatrixXf> eig(A_bug1174);
205 CALL_SUBTEST( check_eigensolver_for_given_mat(eig,A_bug1174) );
206 }
207 }
208
EIGEN_DECLARE_TEST(eigensolver_generic)209 EIGEN_DECLARE_TEST(eigensolver_generic)
210 {
211 int s = 0;
212 for(int i = 0; i < g_repeat; i++) {
213 CALL_SUBTEST_1( eigensolver(Matrix4f()) );
214 s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4);
215 CALL_SUBTEST_2( eigensolver(MatrixXd(s,s)) );
216 TEST_SET_BUT_UNUSED_VARIABLE(s)
217
218 // some trivial but implementation-wise tricky cases
219 CALL_SUBTEST_2( eigensolver(MatrixXd(1,1)) );
220 CALL_SUBTEST_2( eigensolver(MatrixXd(2,2)) );
221 CALL_SUBTEST_3( eigensolver(Matrix<double,1,1>()) );
222 CALL_SUBTEST_4( eigensolver(Matrix2d()) );
223 }
224
225 CALL_SUBTEST_1( eigensolver_verify_assert(Matrix4f()) );
226 s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4);
227 CALL_SUBTEST_2( eigensolver_verify_assert(MatrixXd(s,s)) );
228 CALL_SUBTEST_3( eigensolver_verify_assert(Matrix<double,1,1>()) );
229 CALL_SUBTEST_4( eigensolver_verify_assert(Matrix2d()) );
230
231 // Test problem size constructors
232 CALL_SUBTEST_5(EigenSolver<MatrixXf> tmp(s));
233
234 // regression test for bug 410
235 CALL_SUBTEST_2(
236 {
237 MatrixXd A(1,1);
238 A(0,0) = std::sqrt(-1.); // is Not-a-Number
239 Eigen::EigenSolver<MatrixXd> solver(A);
240 VERIFY_IS_EQUAL(solver.info(), NumericalIssue);
241 }
242 );
243
244 CALL_SUBTEST_2( eigensolver_generic_extra<0>() );
245
246 TEST_SET_BUT_UNUSED_VARIABLE(s)
247 }
248