eigensolver_generic.cpp (9459B)
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 EigType,typename MatType> 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 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 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> 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> 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 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 }