test/eigensolver_selfadjoint.cpp - eigen - Git at Google

 // 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)
 }
	// 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)
	}