123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281 |
- // This file is part of Eigen, a lightweight C++ template library
- // for linear algebra.
- //
- // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
- // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
- //
- // This Source Code Form is subject to the terms of the Mozilla
- // Public License v. 2.0. If a copy of the MPL was not distributed
- // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
- #include "main.h"
- #include "svd_fill.h"
- #include <limits>
- #include <Eigen/Eigenvalues>
- #include <Eigen/SparseCore>
- template<typename MatrixType> void selfadjointeigensolver_essential_check(const MatrixType& m)
- {
- typedef typename MatrixType::Scalar Scalar;
- typedef typename NumTraits<Scalar>::Real RealScalar;
- RealScalar eival_eps = numext::mini<RealScalar>(test_precision<RealScalar>(), NumTraits<Scalar>::dummy_precision()*20000);
-
- SelfAdjointEigenSolver<MatrixType> eiSymm(m);
- VERIFY_IS_EQUAL(eiSymm.info(), Success);
- RealScalar scaling = m.cwiseAbs().maxCoeff();
- if(scaling<(std::numeric_limits<RealScalar>::min)())
- {
- VERIFY(eiSymm.eigenvalues().cwiseAbs().maxCoeff() <= (std::numeric_limits<RealScalar>::min)());
- }
- else
- {
- VERIFY_IS_APPROX((m.template selfadjointView<Lower>() * eiSymm.eigenvectors())/scaling,
- (eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal())/scaling);
- }
- VERIFY_IS_APPROX(m.template selfadjointView<Lower>().eigenvalues(), eiSymm.eigenvalues());
- VERIFY_IS_UNITARY(eiSymm.eigenvectors());
- if(m.cols()<=4)
- {
- SelfAdjointEigenSolver<MatrixType> eiDirect;
- eiDirect.computeDirect(m);
- VERIFY_IS_EQUAL(eiDirect.info(), Success);
- if(! eiSymm.eigenvalues().isApprox(eiDirect.eigenvalues(), eival_eps) )
- {
- std::cerr << "reference eigenvalues: " << eiSymm.eigenvalues().transpose() << "\n"
- << "obtained eigenvalues: " << eiDirect.eigenvalues().transpose() << "\n"
- << "diff: " << (eiSymm.eigenvalues()-eiDirect.eigenvalues()).transpose() << "\n"
- << "error (eps): " << (eiSymm.eigenvalues()-eiDirect.eigenvalues()).norm() / eiSymm.eigenvalues().norm() << " (" << eival_eps << ")\n";
- }
- if(scaling<(std::numeric_limits<RealScalar>::min)())
- {
- VERIFY(eiDirect.eigenvalues().cwiseAbs().maxCoeff() <= (std::numeric_limits<RealScalar>::min)());
- }
- else
- {
- VERIFY_IS_APPROX(eiSymm.eigenvalues()/scaling, eiDirect.eigenvalues()/scaling);
- VERIFY_IS_APPROX((m.template selfadjointView<Lower>() * eiDirect.eigenvectors())/scaling,
- (eiDirect.eigenvectors() * eiDirect.eigenvalues().asDiagonal())/scaling);
- VERIFY_IS_APPROX(m.template selfadjointView<Lower>().eigenvalues()/scaling, eiDirect.eigenvalues()/scaling);
- }
- VERIFY_IS_UNITARY(eiDirect.eigenvectors());
- }
- }
- template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
- {
- /* this test covers the following files:
- EigenSolver.h, SelfAdjointEigenSolver.h (and indirectly: Tridiagonalization.h)
- */
- Index rows = m.rows();
- Index cols = m.cols();
- typedef typename MatrixType::Scalar Scalar;
- typedef typename NumTraits<Scalar>::Real RealScalar;
- RealScalar largerEps = 10*test_precision<RealScalar>();
- MatrixType a = MatrixType::Random(rows,cols);
- MatrixType a1 = MatrixType::Random(rows,cols);
- MatrixType symmA = a.adjoint() * a + a1.adjoint() * a1;
- MatrixType symmC = symmA;
-
- svd_fill_random(symmA,Symmetric);
- symmA.template triangularView<StrictlyUpper>().setZero();
- symmC.template triangularView<StrictlyUpper>().setZero();
- MatrixType b = MatrixType::Random(rows,cols);
- MatrixType b1 = MatrixType::Random(rows,cols);
- MatrixType symmB = b.adjoint() * b + b1.adjoint() * b1;
- symmB.template triangularView<StrictlyUpper>().setZero();
-
- CALL_SUBTEST( selfadjointeigensolver_essential_check(symmA) );
- SelfAdjointEigenSolver<MatrixType> eiSymm(symmA);
- // generalized eigen pb
- GeneralizedSelfAdjointEigenSolver<MatrixType> eiSymmGen(symmC, symmB);
- SelfAdjointEigenSolver<MatrixType> eiSymmNoEivecs(symmA, false);
- VERIFY_IS_EQUAL(eiSymmNoEivecs.info(), Success);
- VERIFY_IS_APPROX(eiSymm.eigenvalues(), eiSymmNoEivecs.eigenvalues());
-
- // generalized eigen problem Ax = lBx
- eiSymmGen.compute(symmC, symmB,Ax_lBx);
- VERIFY_IS_EQUAL(eiSymmGen.info(), Success);
- VERIFY((symmC.template selfadjointView<Lower>() * eiSymmGen.eigenvectors()).isApprox(
- symmB.template selfadjointView<Lower>() * (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
- // generalized eigen problem BAx = lx
- eiSymmGen.compute(symmC, symmB,BAx_lx);
- VERIFY_IS_EQUAL(eiSymmGen.info(), Success);
- VERIFY((symmB.template selfadjointView<Lower>() * (symmC.template selfadjointView<Lower>() * eiSymmGen.eigenvectors())).isApprox(
- (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
- // generalized eigen problem ABx = lx
- eiSymmGen.compute(symmC, symmB,ABx_lx);
- VERIFY_IS_EQUAL(eiSymmGen.info(), Success);
- VERIFY((symmC.template selfadjointView<Lower>() * (symmB.template selfadjointView<Lower>() * eiSymmGen.eigenvectors())).isApprox(
- (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
- eiSymm.compute(symmC);
- MatrixType sqrtSymmA = eiSymm.operatorSqrt();
- VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView<Lower>()), sqrtSymmA*sqrtSymmA);
- VERIFY_IS_APPROX(sqrtSymmA, symmC.template selfadjointView<Lower>()*eiSymm.operatorInverseSqrt());
- MatrixType id = MatrixType::Identity(rows, cols);
- VERIFY_IS_APPROX(id.template selfadjointView<Lower>().operatorNorm(), RealScalar(1));
- SelfAdjointEigenSolver<MatrixType> eiSymmUninitialized;
- VERIFY_RAISES_ASSERT(eiSymmUninitialized.info());
- VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvalues());
- VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvectors());
- VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorSqrt());
- VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorInverseSqrt());
- eiSymmUninitialized.compute(symmA, false);
- VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvectors());
- VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorSqrt());
- VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorInverseSqrt());
- // test Tridiagonalization's methods
- Tridiagonalization<MatrixType> tridiag(symmC);
- VERIFY_IS_APPROX(tridiag.diagonal(), tridiag.matrixT().diagonal());
- VERIFY_IS_APPROX(tridiag.subDiagonal(), tridiag.matrixT().template diagonal<-1>());
- Matrix<RealScalar,Dynamic,Dynamic> T = tridiag.matrixT();
- if(rows>1 && cols>1) {
- // FIXME check that upper and lower part are 0:
- //VERIFY(T.topRightCorner(rows-2, cols-2).template triangularView<Upper>().isZero());
- }
- VERIFY_IS_APPROX(tridiag.diagonal(), T.diagonal());
- VERIFY_IS_APPROX(tridiag.subDiagonal(), T.template diagonal<1>());
- VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView<Lower>()), tridiag.matrixQ() * tridiag.matrixT().eval() * MatrixType(tridiag.matrixQ()).adjoint());
- VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView<Lower>()), tridiag.matrixQ() * tridiag.matrixT() * tridiag.matrixQ().adjoint());
-
- // Test computation of eigenvalues from tridiagonal matrix
- if(rows > 1)
- {
- SelfAdjointEigenSolver<MatrixType> eiSymmTridiag;
- eiSymmTridiag.computeFromTridiagonal(tridiag.matrixT().diagonal(), tridiag.matrixT().diagonal(-1), ComputeEigenvectors);
- VERIFY_IS_APPROX(eiSymm.eigenvalues(), eiSymmTridiag.eigenvalues());
- VERIFY_IS_APPROX(tridiag.matrixT(), eiSymmTridiag.eigenvectors().real() * eiSymmTridiag.eigenvalues().asDiagonal() * eiSymmTridiag.eigenvectors().real().transpose());
- }
- if (rows > 1 && rows < 20)
- {
- // Test matrix with NaN
- symmC(0,0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN();
- SelfAdjointEigenSolver<MatrixType> eiSymmNaN(symmC);
- VERIFY_IS_EQUAL(eiSymmNaN.info(), NoConvergence);
- }
- // regression test for bug 1098
- {
- SelfAdjointEigenSolver<MatrixType> eig(a.adjoint() * a);
- eig.compute(a.adjoint() * a);
- }
- // regression test for bug 478
- {
- a.setZero();
- SelfAdjointEigenSolver<MatrixType> ei3(a);
- VERIFY_IS_EQUAL(ei3.info(), Success);
- VERIFY_IS_MUCH_SMALLER_THAN(ei3.eigenvalues().norm(),RealScalar(1));
- VERIFY((ei3.eigenvectors().transpose()*ei3.eigenvectors().transpose()).eval().isIdentity());
- }
- }
- template<int>
- void bug_854()
- {
- Matrix3d m;
- m << 850.961, 51.966, 0,
- 51.966, 254.841, 0,
- 0, 0, 0;
- selfadjointeigensolver_essential_check(m);
- }
- template<int>
- void bug_1014()
- {
- Matrix3d m;
- m << 0.11111111111111114658, 0, 0,
- 0, 0.11111111111111109107, 0,
- 0, 0, 0.11111111111111107719;
- selfadjointeigensolver_essential_check(m);
- }
- template<int>
- void bug_1225()
- {
- Matrix3d m1, m2;
- m1.setRandom();
- m1 = m1*m1.transpose();
- m2 = m1.triangularView<Upper>();
- SelfAdjointEigenSolver<Matrix3d> eig1(m1);
- SelfAdjointEigenSolver<Matrix3d> eig2(m2.selfadjointView<Upper>());
- VERIFY_IS_APPROX(eig1.eigenvalues(), eig2.eigenvalues());
- }
- template<int>
- void bug_1204()
- {
- SparseMatrix<double> A(2,2);
- A.setIdentity();
- SelfAdjointEigenSolver<Eigen::SparseMatrix<double> > eig(A);
- }
- EIGEN_DECLARE_TEST(eigensolver_selfadjoint)
- {
- int s = 0;
- for(int i = 0; i < g_repeat; i++) {
- // trivial test for 1x1 matrices:
- CALL_SUBTEST_1( selfadjointeigensolver(Matrix<float, 1, 1>()));
- CALL_SUBTEST_1( selfadjointeigensolver(Matrix<double, 1, 1>()));
- CALL_SUBTEST_1( selfadjointeigensolver(Matrix<std::complex<double>, 1, 1>()));
- // very important to test 3x3 and 2x2 matrices since we provide special paths for them
- CALL_SUBTEST_12( selfadjointeigensolver(Matrix2f()) );
- CALL_SUBTEST_12( selfadjointeigensolver(Matrix2d()) );
- CALL_SUBTEST_12( selfadjointeigensolver(Matrix2cd()) );
- CALL_SUBTEST_13( selfadjointeigensolver(Matrix3f()) );
- CALL_SUBTEST_13( selfadjointeigensolver(Matrix3d()) );
- CALL_SUBTEST_13( selfadjointeigensolver(Matrix3cd()) );
- CALL_SUBTEST_2( selfadjointeigensolver(Matrix4d()) );
- CALL_SUBTEST_2( selfadjointeigensolver(Matrix4cd()) );
-
- s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4);
- CALL_SUBTEST_3( selfadjointeigensolver(MatrixXf(s,s)) );
- CALL_SUBTEST_4( selfadjointeigensolver(MatrixXd(s,s)) );
- CALL_SUBTEST_5( selfadjointeigensolver(MatrixXcd(s,s)) );
- CALL_SUBTEST_9( selfadjointeigensolver(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(s,s)) );
- TEST_SET_BUT_UNUSED_VARIABLE(s)
- // some trivial but implementation-wise tricky cases
- CALL_SUBTEST_4( selfadjointeigensolver(MatrixXd(1,1)) );
- CALL_SUBTEST_4( selfadjointeigensolver(MatrixXd(2,2)) );
- CALL_SUBTEST_5( selfadjointeigensolver(MatrixXcd(1,1)) );
- CALL_SUBTEST_5( selfadjointeigensolver(MatrixXcd(2,2)) );
- CALL_SUBTEST_6( selfadjointeigensolver(Matrix<double,1,1>()) );
- CALL_SUBTEST_7( selfadjointeigensolver(Matrix<double,2,2>()) );
- }
-
- CALL_SUBTEST_13( bug_854<0>() );
- CALL_SUBTEST_13( bug_1014<0>() );
- CALL_SUBTEST_13( bug_1204<0>() );
- CALL_SUBTEST_13( bug_1225<0>() );
- // Test problem size constructors
- s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4);
- CALL_SUBTEST_8(SelfAdjointEigenSolver<MatrixXf> tmp1(s));
- CALL_SUBTEST_8(Tridiagonalization<MatrixXf> tmp2(s));
-
- TEST_SET_BUT_UNUSED_VARIABLE(s)
- }
|