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 (at) inria.fr> 5 // Copyright (C) 2010 Jitse Niesen <jitse (at) 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 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 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> 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> 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> 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> 227 void bug_1204() 228 { 229 SparseMatrix<double> A(2,2); 230 A.setIdentity(); 231 SelfAdjointEigenSolver<Eigen::SparseMatrix<double> > eig(A); 232 } 233 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
