1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
5 // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
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 "svd_fill.h"
13 #include <limits>
14 #include <Eigen/Eigenvalues>
15 #include <Eigen/SparseCore>
16
17
selfadjointeigensolver_essential_check(const MatrixType & m)18 template<typename MatrixType> void selfadjointeigensolver_essential_check(const MatrixType& m)
19 {
20 typedef typename MatrixType::Scalar Scalar;
21 typedef typename NumTraits<Scalar>::Real RealScalar;
22 RealScalar eival_eps = numext::mini<RealScalar>(test_precision<RealScalar>(), NumTraits<Scalar>::dummy_precision()*20000);
23
24 SelfAdjointEigenSolver<MatrixType> eiSymm(m);
25 VERIFY_IS_EQUAL(eiSymm.info(), Success);
26
27 RealScalar scaling = m.cwiseAbs().maxCoeff();
28
29 if(scaling<(std::numeric_limits<RealScalar>::min)())
30 {
31 VERIFY(eiSymm.eigenvalues().cwiseAbs().maxCoeff() <= (std::numeric_limits<RealScalar>::min)());
32 }
33 else
34 {
35 VERIFY_IS_APPROX((m.template selfadjointView<Lower>() * eiSymm.eigenvectors())/scaling,
36 (eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal())/scaling);
37 }
38 VERIFY_IS_APPROX(m.template selfadjointView<Lower>().eigenvalues(), eiSymm.eigenvalues());
39 VERIFY_IS_UNITARY(eiSymm.eigenvectors());
40
41 if(m.cols()<=4)
42 {
43 SelfAdjointEigenSolver<MatrixType> eiDirect;
44 eiDirect.computeDirect(m);
45 VERIFY_IS_EQUAL(eiDirect.info(), Success);
46 if(! eiSymm.eigenvalues().isApprox(eiDirect.eigenvalues(), eival_eps) )
47 {
48 std::cerr << "reference eigenvalues: " << eiSymm.eigenvalues().transpose() << "\n"
49 << "obtained eigenvalues: " << eiDirect.eigenvalues().transpose() << "\n"
50 << "diff: " << (eiSymm.eigenvalues()-eiDirect.eigenvalues()).transpose() << "\n"
51 << "error (eps): " << (eiSymm.eigenvalues()-eiDirect.eigenvalues()).norm() / eiSymm.eigenvalues().norm() << " (" << eival_eps << ")\n";
52 }
53 if(scaling<(std::numeric_limits<RealScalar>::min)())
54 {
55 VERIFY(eiDirect.eigenvalues().cwiseAbs().maxCoeff() <= (std::numeric_limits<RealScalar>::min)());
56 }
57 else
58 {
59 VERIFY_IS_APPROX(eiSymm.eigenvalues()/scaling, eiDirect.eigenvalues()/scaling);
60 VERIFY_IS_APPROX((m.template selfadjointView<Lower>() * eiDirect.eigenvectors())/scaling,
61 (eiDirect.eigenvectors() * eiDirect.eigenvalues().asDiagonal())/scaling);
62 VERIFY_IS_APPROX(m.template selfadjointView<Lower>().eigenvalues()/scaling, eiDirect.eigenvalues()/scaling);
63 }
64
65 VERIFY_IS_UNITARY(eiDirect.eigenvectors());
66 }
67 }
68
selfadjointeigensolver(const MatrixType & m)69 template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
70 {
71 typedef typename MatrixType::Index Index;
72 /* this test covers the following files:
73 EigenSolver.h, SelfAdjointEigenSolver.h (and indirectly: Tridiagonalization.h)
74 */
75 Index rows = m.rows();
76 Index cols = m.cols();
77
78 typedef typename MatrixType::Scalar Scalar;
79 typedef typename NumTraits<Scalar>::Real RealScalar;
80
81 RealScalar largerEps = 10*test_precision<RealScalar>();
82
83 MatrixType a = MatrixType::Random(rows,cols);
84 MatrixType a1 = MatrixType::Random(rows,cols);
85 MatrixType symmA = a.adjoint() * a + a1.adjoint() * a1;
86 MatrixType symmC = symmA;
87
88 svd_fill_random(symmA,Symmetric);
89
90 symmA.template triangularView<StrictlyUpper>().setZero();
91 symmC.template triangularView<StrictlyUpper>().setZero();
92
93 MatrixType b = MatrixType::Random(rows,cols);
94 MatrixType b1 = MatrixType::Random(rows,cols);
95 MatrixType symmB = b.adjoint() * b + b1.adjoint() * b1;
96 symmB.template triangularView<StrictlyUpper>().setZero();
97
98 CALL_SUBTEST( selfadjointeigensolver_essential_check(symmA) );
99
100 SelfAdjointEigenSolver<MatrixType> eiSymm(symmA);
101 // generalized eigen pb
102 GeneralizedSelfAdjointEigenSolver<MatrixType> eiSymmGen(symmC, symmB);
103
104 SelfAdjointEigenSolver<MatrixType> eiSymmNoEivecs(symmA, false);
105 VERIFY_IS_EQUAL(eiSymmNoEivecs.info(), Success);
106 VERIFY_IS_APPROX(eiSymm.eigenvalues(), eiSymmNoEivecs.eigenvalues());
107
108 // generalized eigen problem Ax = lBx
109 eiSymmGen.compute(symmC, symmB,Ax_lBx);
110 VERIFY_IS_EQUAL(eiSymmGen.info(), Success);
111 VERIFY((symmC.template selfadjointView<Lower>() * eiSymmGen.eigenvectors()).isApprox(
112 symmB.template selfadjointView<Lower>() * (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
113
114 // generalized eigen problem BAx = lx
115 eiSymmGen.compute(symmC, symmB,BAx_lx);
116 VERIFY_IS_EQUAL(eiSymmGen.info(), Success);
117 VERIFY((symmB.template selfadjointView<Lower>() * (symmC.template selfadjointView<Lower>() * eiSymmGen.eigenvectors())).isApprox(
118 (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
119
120 // generalized eigen problem ABx = lx
121 eiSymmGen.compute(symmC, symmB,ABx_lx);
122 VERIFY_IS_EQUAL(eiSymmGen.info(), Success);
123 VERIFY((symmC.template selfadjointView<Lower>() * (symmB.template selfadjointView<Lower>() * eiSymmGen.eigenvectors())).isApprox(
124 (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
125
126
127 eiSymm.compute(symmC);
128 MatrixType sqrtSymmA = eiSymm.operatorSqrt();
129 VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView<Lower>()), sqrtSymmA*sqrtSymmA);
130 VERIFY_IS_APPROX(sqrtSymmA, symmC.template selfadjointView<Lower>()*eiSymm.operatorInverseSqrt());
131
132 MatrixType id = MatrixType::Identity(rows, cols);
133 VERIFY_IS_APPROX(id.template selfadjointView<Lower>().operatorNorm(), RealScalar(1));
134
135 SelfAdjointEigenSolver<MatrixType> eiSymmUninitialized;
136 VERIFY_RAISES_ASSERT(eiSymmUninitialized.info());
137 VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvalues());
138 VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvectors());
139 VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorSqrt());
140 VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorInverseSqrt());
141
142 eiSymmUninitialized.compute(symmA, false);
143 VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvectors());
144 VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorSqrt());
145 VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorInverseSqrt());
146
147 // test Tridiagonalization's methods
148 Tridiagonalization<MatrixType> tridiag(symmC);
149 VERIFY_IS_APPROX(tridiag.diagonal(), tridiag.matrixT().diagonal());
150 VERIFY_IS_APPROX(tridiag.subDiagonal(), tridiag.matrixT().template diagonal<-1>());
151 Matrix<RealScalar,Dynamic,Dynamic> T = tridiag.matrixT();
152 if(rows>1 && cols>1) {
153 // FIXME check that upper and lower part are 0:
154 //VERIFY(T.topRightCorner(rows-2, cols-2).template triangularView<Upper>().isZero());
155 }
156 VERIFY_IS_APPROX(tridiag.diagonal(), T.diagonal());
157 VERIFY_IS_APPROX(tridiag.subDiagonal(), T.template diagonal<1>());
158 VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView<Lower>()), tridiag.matrixQ() * tridiag.matrixT().eval() * MatrixType(tridiag.matrixQ()).adjoint());
159 VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView<Lower>()), tridiag.matrixQ() * tridiag.matrixT() * tridiag.matrixQ().adjoint());
160
161 // Test computation of eigenvalues from tridiagonal matrix
162 if(rows > 1)
163 {
164 SelfAdjointEigenSolver<MatrixType> eiSymmTridiag;
165 eiSymmTridiag.computeFromTridiagonal(tridiag.matrixT().diagonal(), tridiag.matrixT().diagonal(-1), ComputeEigenvectors);
166 VERIFY_IS_APPROX(eiSymm.eigenvalues(), eiSymmTridiag.eigenvalues());
167 VERIFY_IS_APPROX(tridiag.matrixT(), eiSymmTridiag.eigenvectors().real() * eiSymmTridiag.eigenvalues().asDiagonal() * eiSymmTridiag.eigenvectors().real().transpose());
168 }
169
170 if (rows > 1 && rows < 20)
171 {
172 // Test matrix with NaN
173 symmC(0,0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN();
174 SelfAdjointEigenSolver<MatrixType> eiSymmNaN(symmC);
175 VERIFY_IS_EQUAL(eiSymmNaN.info(), NoConvergence);
176 }
177
178 // regression test for bug 1098
179 {
180 SelfAdjointEigenSolver<MatrixType> eig(a.adjoint() * a);
181 eig.compute(a.adjoint() * a);
182 }
183
184 // regression test for bug 478
185 {
186 a.setZero();
187 SelfAdjointEigenSolver<MatrixType> ei3(a);
188 VERIFY_IS_EQUAL(ei3.info(), Success);
189 VERIFY_IS_MUCH_SMALLER_THAN(ei3.eigenvalues().norm(),RealScalar(1));
190 VERIFY((ei3.eigenvectors().transpose()*ei3.eigenvectors().transpose()).eval().isIdentity());
191 }
192 }
193
194 template<int>
bug_854()195 void bug_854()
196 {
197 Matrix3d m;
198 m << 850.961, 51.966, 0,
199 51.966, 254.841, 0,
200 0, 0, 0;
201 selfadjointeigensolver_essential_check(m);
202 }
203
204 template<int>
bug_1014()205 void bug_1014()
206 {
207 Matrix3d m;
208 m << 0.11111111111111114658, 0, 0,
209 0, 0.11111111111111109107, 0,
210 0, 0, 0.11111111111111107719;
211 selfadjointeigensolver_essential_check(m);
212 }
213
214 template<int>
bug_1225()215 void bug_1225()
216 {
217 Matrix3d m1, m2;
218 m1.setRandom();
219 m1 = m1*m1.transpose();
220 m2 = m1.triangularView<Upper>();
221 SelfAdjointEigenSolver<Matrix3d> eig1(m1);
222 SelfAdjointEigenSolver<Matrix3d> eig2(m2.selfadjointView<Upper>());
223 VERIFY_IS_APPROX(eig1.eigenvalues(), eig2.eigenvalues());
224 }
225
226 template<int>
bug_1204()227 void bug_1204()
228 {
229 SparseMatrix<double> A(2,2);
230 A.setIdentity();
231 SelfAdjointEigenSolver<Eigen::SparseMatrix<double> > eig(A);
232 }
233
test_eigensolver_selfadjoint()234 void test_eigensolver_selfadjoint()
235 {
236 int s = 0;
237 for(int i = 0; i < g_repeat; i++) {
238 // trivial test for 1x1 matrices:
239 CALL_SUBTEST_1( selfadjointeigensolver(Matrix<float, 1, 1>()));
240 CALL_SUBTEST_1( selfadjointeigensolver(Matrix<double, 1, 1>()));
241 // very important to test 3x3 and 2x2 matrices since we provide special paths for them
242 CALL_SUBTEST_12( selfadjointeigensolver(Matrix2f()) );
243 CALL_SUBTEST_12( selfadjointeigensolver(Matrix2d()) );
244 CALL_SUBTEST_13( selfadjointeigensolver(Matrix3f()) );
245 CALL_SUBTEST_13( selfadjointeigensolver(Matrix3d()) );
246 CALL_SUBTEST_2( selfadjointeigensolver(Matrix4d()) );
247
248 s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4);
249 CALL_SUBTEST_3( selfadjointeigensolver(MatrixXf(s,s)) );
250 CALL_SUBTEST_4( selfadjointeigensolver(MatrixXd(s,s)) );
251 CALL_SUBTEST_5( selfadjointeigensolver(MatrixXcd(s,s)) );
252 CALL_SUBTEST_9( selfadjointeigensolver(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(s,s)) );
253 TEST_SET_BUT_UNUSED_VARIABLE(s)
254
255 // some trivial but implementation-wise tricky cases
256 CALL_SUBTEST_4( selfadjointeigensolver(MatrixXd(1,1)) );
257 CALL_SUBTEST_4( selfadjointeigensolver(MatrixXd(2,2)) );
258 CALL_SUBTEST_6( selfadjointeigensolver(Matrix<double,1,1>()) );
259 CALL_SUBTEST_7( selfadjointeigensolver(Matrix<double,2,2>()) );
260 }
261
262 CALL_SUBTEST_13( bug_854<0>() );
263 CALL_SUBTEST_13( bug_1014<0>() );
264 CALL_SUBTEST_13( bug_1204<0>() );
265 CALL_SUBTEST_13( bug_1225<0>() );
266
267 // Test problem size constructors
268 s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4);
269 CALL_SUBTEST_8(SelfAdjointEigenSolver<MatrixXf> tmp1(s));
270 CALL_SUBTEST_8(Tridiagonalization<MatrixXf> tmp2(s));
271
272 TEST_SET_BUT_UNUSED_VARIABLE(s)
273 }
274
275