eigensolver_generic.cpp
Go to the documentation of this file.
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 
15 template<typename MatrixType> void eigensolver(const MatrixType& m)
16 {
17  /* this test covers the following files:
18  EigenSolver.h
19  */
20  Index rows = m.rows();
21  Index cols = m.cols();
22 
23  typedef typename MatrixType::Scalar Scalar;
24  typedef typename NumTraits<Scalar>::Real RealScalar;
27 
28  MatrixType a = MatrixType::Random(rows,cols);
29  MatrixType a1 = MatrixType::Random(rows,cols);
30  MatrixType symmA = a.adjoint() * a + a1.adjoint() * a1;
31 
32  EigenSolver<MatrixType> ei0(symmA);
35  VERIFY_IS_APPROX((symmA.template cast<Complex>()) * (ei0.pseudoEigenvectors().template cast<Complex>()),
36  (ei0.pseudoEigenvectors().template cast<Complex>()) * (ei0.eigenvalues().asDiagonal()));
37 
41  VERIFY_IS_APPROX(a.template cast<Complex>() * ei1.eigenvectors(),
42  ei1.eigenvectors() * ei1.eigenvalues().asDiagonal());
43  VERIFY_IS_APPROX(ei1.eigenvectors().colwise().norm(), RealVectorType::Ones(rows).transpose());
44  VERIFY_IS_APPROX(a.eigenvalues(), ei1.eigenvalues());
45 
48  VERIFY_IS_EQUAL(ei2.info(), Success);
49  VERIFY_IS_EQUAL(ei2.eigenvectors(), ei1.eigenvectors());
50  VERIFY_IS_EQUAL(ei2.eigenvalues(), ei1.eigenvalues());
51  if (rows > 2) {
52  ei2.setMaxIterations(1).compute(a);
53  VERIFY_IS_EQUAL(ei2.info(), NoConvergence);
54  VERIFY_IS_EQUAL(ei2.getMaxIterations(), 1);
55  }
56 
57  EigenSolver<MatrixType> eiNoEivecs(a, false);
58  VERIFY_IS_EQUAL(eiNoEivecs.info(), Success);
59  VERIFY_IS_APPROX(ei1.eigenvalues(), eiNoEivecs.eigenvalues());
60  VERIFY_IS_APPROX(ei1.pseudoEigenvalueMatrix(), eiNoEivecs.pseudoEigenvalueMatrix());
61 
62  MatrixType id = MatrixType::Identity(rows, cols);
63  VERIFY_IS_APPROX(id.operatorNorm(), RealScalar(1));
64 
65  if (rows > 2 && rows < 20)
66  {
67  // Test matrix with NaN
68  a(0,0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN();
69  EigenSolver<MatrixType> eiNaN(a);
71  }
72 
73  // regression test for bug 1098
74  {
75  EigenSolver<MatrixType> eig(a.adjoint() * a);
76  eig.compute(a.adjoint() * a);
77  }
78 
79  // regression test for bug 478
80  {
81  a.setZero();
85  VERIFY((ei3.eigenvectors().transpose()*ei3.eigenvectors().transpose()).eval().isIdentity());
86  }
87 }
88 
89 template<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m)
90 {
96 
97  MatrixType a = MatrixType::Random(m.rows(),m.cols());
98  eig.compute(a, false);
101 }
102 
104 {
105  int s = 0;
106  for(int i = 0; i < g_repeat; i++) {
107  CALL_SUBTEST_1( eigensolver(Matrix4f()) );
108  s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4);
109  CALL_SUBTEST_2( eigensolver(MatrixXd(s,s)) );
111 
112  // some trivial but implementation-wise tricky cases
113  CALL_SUBTEST_2( eigensolver(MatrixXd(1,1)) );
114  CALL_SUBTEST_2( eigensolver(MatrixXd(2,2)) );
115  CALL_SUBTEST_3( eigensolver(Matrix<double,1,1>()) );
116  CALL_SUBTEST_4( eigensolver(Matrix2d()) );
117  }
118 
119  CALL_SUBTEST_1( eigensolver_verify_assert(Matrix4f()) );
120  s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4);
121  CALL_SUBTEST_2( eigensolver_verify_assert(MatrixXd(s,s)) );
122  CALL_SUBTEST_3( eigensolver_verify_assert(Matrix<double,1,1>()) );
123  CALL_SUBTEST_4( eigensolver_verify_assert(Matrix2d()) );
124 
125  // Test problem size constructors
126  CALL_SUBTEST_5(EigenSolver<MatrixXf> tmp(s));
127 
128  // regression test for bug 410
129  CALL_SUBTEST_2(
130  {
131  MatrixXd A(1,1);
132  A(0,0) = std::sqrt(-1.); // is Not-a-Number
135  }
136  );
137 
138 #ifdef EIGEN_TEST_PART_2
139  {
140  // regression test for bug 793
141  MatrixXd a(3,3);
142  a << 0, 0, 1,
143  1, 1, 1,
144  1, 1e+200, 1;
146  double scale = 1e-200; // scale to avoid overflow during the comparisons
148  VERIFY_IS_APPROX(a * eig.eigenvectors()*scale, eig.eigenvectors() * eig.eigenvalues().asDiagonal()*scale);
149  }
150  {
151  // check a case where all eigenvalues are null.
152  MatrixXd a(2,2);
153  a << 1, 1,
154  -1, -1;
156  VERIFY_IS_APPROX(eig.pseudoEigenvectors().squaredNorm(), 2.);
157  VERIFY_IS_APPROX((a * eig.pseudoEigenvectors()).norm()+1., 1.);
158  VERIFY_IS_APPROX((eig.pseudoEigenvectors() * eig.pseudoEigenvalueMatrix()).norm()+1., 1.);
159  VERIFY_IS_APPROX((a * eig.eigenvectors()).norm()+1., 1.);
160  VERIFY_IS_APPROX((eig.eigenvectors() * eig.eigenvalues().asDiagonal()).norm()+1., 1.);
161  }
162 #endif
163 
165 }
Matrix3f m
SCALAR Scalar
Definition: bench_gemm.cpp:33
#define VERIFY_RAISES_ASSERT(a)
Definition: main.h:285
EigenSolver & compute(const EigenBase< InputType > &matrix, bool computeEigenvectors=true)
Computes eigendecomposition of given matrix.
EigenvectorsType eigenvectors() const
Returns the eigenvectors of given matrix.
Definition: EigenSolver.h:345
void eigensolver_verify_assert(const MatrixType &m)
void test_eigensolver_generic()
Performs a real Schur decomposition of a square matrix.
Definition: RealSchur.h:54
EIGEN_DEVICE_FUNC const SqrtReturnType sqrt() const
MatrixXf MatrixType
BiCGSTAB< SparseMatrix< double > > solver
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:150
Matrix< SCALARA, Dynamic, Dynamic > A
Definition: bench_gemm.cpp:35
ComputationInfo info() const
Definition: EigenSolver.h:281
std::complex< RealScalar > Complex
Array33i a
MatrixType pseudoEigenvalueMatrix() const
Returns the block-diagonal matrix in the pseudo-eigendecomposition.
Definition: EigenSolver.h:324
#define VERIFY_IS_APPROX(a, b)
#define VERIFY_IS_EQUAL(a, b)
Definition: main.h:331
const MatrixType & pseudoEigenvectors() const
Returns the pseudo-eigenvectors of given matrix.
Definition: EigenSolver.h:199
static int g_repeat
Definition: main.h:144
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
SelfAdjointEigenSolver< PlainMatrixType > eig(mat, computeVectors?ComputeEigenvectors:EigenvaluesOnly)
Array< double, 1, 3 > e(1./3., 0.5, 2.)
RealScalar s
EigenSolver & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition: EigenSolver.h:288
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:34
#define VERIFY_IS_MUCH_SMALLER_THAN(a, b)
Definition: main.h:335
mp::number< mp::cpp_dec_float< 100 >, mp::et_on > Real
#define TEST_SET_BUT_UNUSED_VARIABLE(X)
Definition: main.h:91
#define VERIFY(a)
Definition: main.h:325
#define EIGEN_TEST_MAX_SIZE
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy y set format x g set format y g set format x2 g set format y2 g set format z g set angles radians set nogrid set key title set key left top Right noreverse box linetype linewidth samplen spacing width set nolabel set noarrow set nologscale set logscale x set set pointsize set encoding default set nopolar set noparametric set set set set surface set nocontour set clabel set mapping cartesian set nohidden3d set cntrparam order set cntrparam linear set cntrparam levels auto set cntrparam points set size set set xzeroaxis lt lw set x2zeroaxis lt lw set yzeroaxis lt lw set y2zeroaxis lt lw set tics in set ticslevel set tics scale
internal::nested_eval< T, 1 >::type eval(const T &xpr)
Computes eigenvalues and eigenvectors of general matrices.
Definition: EigenSolver.h:64
The matrix class, also used for vectors and row-vectors.
void eigensolver(const MatrixType &m)
const EigenvalueType & eigenvalues() const
Returns the eigenvalues of given matrix.
Definition: EigenSolver.h:244


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:42:01