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,2012 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 <limits>
13 #include <Eigen/Eigenvalues>
14 
eigensolver(const MatrixType & m)15 template<typename MatrixType> void eigensolver(const MatrixType& m)
16 {
17   typedef typename MatrixType::Index Index;
18   /* this test covers the following files:
19      EigenSolver.h
20   */
21   Index rows = m.rows();
22   Index cols = m.cols();
23 
24   typedef typename MatrixType::Scalar Scalar;
25   typedef typename NumTraits<Scalar>::Real RealScalar;
26   typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealVectorType;
27   typedef typename std::complex<typename NumTraits<typename MatrixType::Scalar>::Real> Complex;
28 
29   MatrixType a = MatrixType::Random(rows,cols);
30   MatrixType a1 = MatrixType::Random(rows,cols);
31   MatrixType symmA =  a.adjoint() * a + a1.adjoint() * a1;
32 
33   EigenSolver<MatrixType> ei0(symmA);
34   VERIFY_IS_EQUAL(ei0.info(), Success);
35   VERIFY_IS_APPROX(symmA * ei0.pseudoEigenvectors(), ei0.pseudoEigenvectors() * ei0.pseudoEigenvalueMatrix());
36   VERIFY_IS_APPROX((symmA.template cast<Complex>()) * (ei0.pseudoEigenvectors().template cast<Complex>()),
37     (ei0.pseudoEigenvectors().template cast<Complex>()) * (ei0.eigenvalues().asDiagonal()));
38 
39   EigenSolver<MatrixType> ei1(a);
40   VERIFY_IS_EQUAL(ei1.info(), Success);
41   VERIFY_IS_APPROX(a * ei1.pseudoEigenvectors(), ei1.pseudoEigenvectors() * ei1.pseudoEigenvalueMatrix());
42   VERIFY_IS_APPROX(a.template cast<Complex>() * ei1.eigenvectors(),
43                    ei1.eigenvectors() * ei1.eigenvalues().asDiagonal());
44   VERIFY_IS_APPROX(ei1.eigenvectors().colwise().norm(), RealVectorType::Ones(rows).transpose());
45   VERIFY_IS_APPROX(a.eigenvalues(), ei1.eigenvalues());
46 
47   EigenSolver<MatrixType> ei2;
48   ei2.setMaxIterations(RealSchur<MatrixType>::m_maxIterationsPerRow * rows).compute(a);
49   VERIFY_IS_EQUAL(ei2.info(), Success);
50   VERIFY_IS_EQUAL(ei2.eigenvectors(), ei1.eigenvectors());
51   VERIFY_IS_EQUAL(ei2.eigenvalues(), ei1.eigenvalues());
52   if (rows > 2) {
53     ei2.setMaxIterations(1).compute(a);
54     VERIFY_IS_EQUAL(ei2.info(), NoConvergence);
55     VERIFY_IS_EQUAL(ei2.getMaxIterations(), 1);
56   }
57 
58   EigenSolver<MatrixType> eiNoEivecs(a, false);
59   VERIFY_IS_EQUAL(eiNoEivecs.info(), Success);
60   VERIFY_IS_APPROX(ei1.eigenvalues(), eiNoEivecs.eigenvalues());
61   VERIFY_IS_APPROX(ei1.pseudoEigenvalueMatrix(), eiNoEivecs.pseudoEigenvalueMatrix());
62 
63   MatrixType id = MatrixType::Identity(rows, cols);
64   VERIFY_IS_APPROX(id.operatorNorm(), RealScalar(1));
65 
66   if (rows > 2 && rows < 20)
67   {
68     // Test matrix with NaN
69     a(0,0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN();
70     EigenSolver<MatrixType> eiNaN(a);
71     VERIFY_IS_EQUAL(eiNaN.info(), NoConvergence);
72   }
73 
74   // regression test for bug 1098
75   {
76     EigenSolver<MatrixType> eig(a.adjoint() * a);
77     eig.compute(a.adjoint() * a);
78   }
79 
80   // regression test for bug 478
81   {
82     a.setZero();
83     EigenSolver<MatrixType> ei3(a);
84     VERIFY_IS_EQUAL(ei3.info(), Success);
85     VERIFY_IS_MUCH_SMALLER_THAN(ei3.eigenvalues().norm(),RealScalar(1));
86     VERIFY((ei3.eigenvectors().transpose()*ei3.eigenvectors().transpose()).eval().isIdentity());
87   }
88 }
89 
eigensolver_verify_assert(const MatrixType & m)90 template<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m)
91 {
92   EigenSolver<MatrixType> eig;
93   VERIFY_RAISES_ASSERT(eig.eigenvectors());
94   VERIFY_RAISES_ASSERT(eig.pseudoEigenvectors());
95   VERIFY_RAISES_ASSERT(eig.pseudoEigenvalueMatrix());
96   VERIFY_RAISES_ASSERT(eig.eigenvalues());
97 
98   MatrixType a = MatrixType::Random(m.rows(),m.cols());
99   eig.compute(a, false);
100   VERIFY_RAISES_ASSERT(eig.eigenvectors());
101   VERIFY_RAISES_ASSERT(eig.pseudoEigenvectors());
102 }
103 
test_eigensolver_generic()104 void test_eigensolver_generic()
105 {
106   int s = 0;
107   for(int i = 0; i < g_repeat; i++) {
108     CALL_SUBTEST_1( eigensolver(Matrix4f()) );
109     s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4);
110     CALL_SUBTEST_2( eigensolver(MatrixXd(s,s)) );
111     TEST_SET_BUT_UNUSED_VARIABLE(s)
112 
113     // some trivial but implementation-wise tricky cases
114     CALL_SUBTEST_2( eigensolver(MatrixXd(1,1)) );
115     CALL_SUBTEST_2( eigensolver(MatrixXd(2,2)) );
116     CALL_SUBTEST_3( eigensolver(Matrix<double,1,1>()) );
117     CALL_SUBTEST_4( eigensolver(Matrix2d()) );
118   }
119 
120   CALL_SUBTEST_1( eigensolver_verify_assert(Matrix4f()) );
121   s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4);
122   CALL_SUBTEST_2( eigensolver_verify_assert(MatrixXd(s,s)) );
123   CALL_SUBTEST_3( eigensolver_verify_assert(Matrix<double,1,1>()) );
124   CALL_SUBTEST_4( eigensolver_verify_assert(Matrix2d()) );
125 
126   // Test problem size constructors
127   CALL_SUBTEST_5(EigenSolver<MatrixXf> tmp(s));
128 
129   // regression test for bug 410
130   CALL_SUBTEST_2(
131   {
132      MatrixXd A(1,1);
133      A(0,0) = std::sqrt(-1.); // is Not-a-Number
134      Eigen::EigenSolver<MatrixXd> solver(A);
135      VERIFY_IS_EQUAL(solver.info(), NumericalIssue);
136   }
137   );
138 
139 #ifdef EIGEN_TEST_PART_2
140   {
141     // regression test for bug 793
142     MatrixXd a(3,3);
143     a << 0,  0,  1,
144         1,  1, 1,
145         1, 1e+200,  1;
146     Eigen::EigenSolver<MatrixXd> eig(a);
147     double scale = 1e-200; // scale to avoid overflow during the comparisons
148     VERIFY_IS_APPROX(a * eig.pseudoEigenvectors()*scale, eig.pseudoEigenvectors() * eig.pseudoEigenvalueMatrix()*scale);
149     VERIFY_IS_APPROX(a * eig.eigenvectors()*scale, eig.eigenvectors() * eig.eigenvalues().asDiagonal()*scale);
150   }
151   {
152     // check a case where all eigenvalues are null.
153     MatrixXd a(2,2);
154     a << 1,  1,
155         -1, -1;
156     Eigen::EigenSolver<MatrixXd> eig(a);
157     VERIFY_IS_APPROX(eig.pseudoEigenvectors().squaredNorm(), 2.);
158     VERIFY_IS_APPROX((a * eig.pseudoEigenvectors()).norm()+1., 1.);
159     VERIFY_IS_APPROX((eig.pseudoEigenvectors() * eig.pseudoEigenvalueMatrix()).norm()+1., 1.);
160     VERIFY_IS_APPROX((a * eig.eigenvectors()).norm()+1., 1.);
161     VERIFY_IS_APPROX((eig.eigenvectors() * eig.eigenvalues().asDiagonal()).norm()+1., 1.);
162   }
163 #endif
164 
165   TEST_SET_BUT_UNUSED_VARIABLE(s)
166 }
167