1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> 57faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk> 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include "main.h" 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include <limits> 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include <Eigen/Eigenvalues> 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> void eigensolver(const MatrixType& m) 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* this test covers the following files: 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenSolver.h 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index rows = m.rows(); 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index cols = m.cols(); 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename NumTraits<Scalar>::Real RealScalar; 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealVectorType; 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename std::complex<typename NumTraits<typename MatrixType::Scalar>::Real> Complex; 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType a = MatrixType::Random(rows,cols); 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType a1 = MatrixType::Random(rows,cols); 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType symmA = a.adjoint() * a + a1.adjoint() * a1; 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenSolver<MatrixType> ei0(symmA); 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(ei0.info(), Success); 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(symmA * ei0.pseudoEigenvectors(), ei0.pseudoEigenvectors() * ei0.pseudoEigenvalueMatrix()); 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX((symmA.template cast<Complex>()) * (ei0.pseudoEigenvectors().template cast<Complex>()), 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath (ei0.pseudoEigenvectors().template cast<Complex>()) * (ei0.eigenvalues().asDiagonal())); 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenSolver<MatrixType> ei1(a); 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(ei1.info(), Success); 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(a * ei1.pseudoEigenvectors(), ei1.pseudoEigenvectors() * ei1.pseudoEigenvalueMatrix()); 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(a.template cast<Complex>() * ei1.eigenvectors(), 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ei1.eigenvectors() * ei1.eigenvalues().asDiagonal()); 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(ei1.eigenvectors().colwise().norm(), RealVectorType::Ones(rows).transpose()); 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(a.eigenvalues(), ei1.eigenvalues()); 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez EigenSolver<MatrixType> ei2; 487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ei2.setMaxIterations(RealSchur<MatrixType>::m_maxIterationsPerRow * rows).compute(a); 497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(ei2.info(), Success); 507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(ei2.eigenvectors(), ei1.eigenvectors()); 517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(ei2.eigenvalues(), ei1.eigenvalues()); 527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (rows > 2) { 537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ei2.setMaxIterations(1).compute(a); 547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(ei2.info(), NoConvergence); 557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(ei2.getMaxIterations(), 1); 567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenSolver<MatrixType> eiNoEivecs(a, false); 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(eiNoEivecs.info(), Success); 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(ei1.eigenvalues(), eiNoEivecs.eigenvalues()); 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(ei1.pseudoEigenvalueMatrix(), eiNoEivecs.pseudoEigenvalueMatrix()); 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType id = MatrixType::Identity(rows, cols); 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(id.operatorNorm(), RealScalar(1)); 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (rows > 2 && rows < 20) 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Test matrix with NaN 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath a(0,0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN(); 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenSolver<MatrixType> eiNaN(a); 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(eiNaN.info(), NoConvergence); 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // regression test for bug 1098 752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang EigenSolver<MatrixType> eig(a.adjoint() * a); 772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eig.compute(a.adjoint() * a); 782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // regression test for bug 478 812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang a.setZero(); 832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang EigenSolver<MatrixType> ei3(a); 842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_EQUAL(ei3.info(), Success); 852b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_MUCH_SMALLER_THAN(ei3.eigenvalues().norm(),RealScalar(1)); 862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY((ei3.eigenvectors().transpose()*ei3.eigenvectors().transpose()).eval().isIdentity()); 872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m) 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenSolver<MatrixType> eig; 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(eig.eigenvectors()); 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(eig.pseudoEigenvectors()); 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(eig.pseudoEigenvalueMatrix()); 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(eig.eigenvalues()); 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType a = MatrixType::Random(m.rows(),m.cols()); 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eig.compute(a, false); 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(eig.eigenvectors()); 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(eig.pseudoEigenvectors()); 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid test_eigensolver_generic() 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int s = 0; 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i = 0; i < g_repeat; i++) { 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_1( eigensolver(Matrix4f()) ); 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4); 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_2( eigensolver(MatrixXd(s,s)) ); 1112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang TEST_SET_BUT_UNUSED_VARIABLE(s) 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // some trivial but implementation-wise tricky cases 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_2( eigensolver(MatrixXd(1,1)) ); 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_2( eigensolver(MatrixXd(2,2)) ); 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_3( eigensolver(Matrix<double,1,1>()) ); 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_4( eigensolver(Matrix2d()) ); 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_1( eigensolver_verify_assert(Matrix4f()) ); 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4); 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_2( eigensolver_verify_assert(MatrixXd(s,s)) ); 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_3( eigensolver_verify_assert(Matrix<double,1,1>()) ); 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_4( eigensolver_verify_assert(Matrix2d()) ); 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Test problem size constructors 1277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST_5(EigenSolver<MatrixXf> tmp(s)); 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // regression test for bug 410 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_2( 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixXd A(1,1); 1332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang A(0,0) = std::sqrt(-1.); // is Not-a-Number 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Eigen::EigenSolver<MatrixXd> solver(A); 1352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_EQUAL(solver.info(), NumericalIssue); 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ); 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#ifdef EIGEN_TEST_PART_2 1402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 1412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // regression test for bug 793 1422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang MatrixXd a(3,3); 1432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang a << 0, 0, 1, 1442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1, 1, 1, 1452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1, 1e+200, 1; 1462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Eigen::EigenSolver<MatrixXd> eig(a); 1472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang double scale = 1e-200; // scale to avoid overflow during the comparisons 1482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_APPROX(a * eig.pseudoEigenvectors()*scale, eig.pseudoEigenvectors() * eig.pseudoEigenvalueMatrix()*scale); 1492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_APPROX(a * eig.eigenvectors()*scale, eig.eigenvectors() * eig.eigenvalues().asDiagonal()*scale); 1502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 1512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 1522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // check a case where all eigenvalues are null. 1532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang MatrixXd a(2,2); 1542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang a << 1, 1, 1552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang -1, -1; 1562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Eigen::EigenSolver<MatrixXd> eig(a); 1572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_APPROX(eig.pseudoEigenvectors().squaredNorm(), 2.); 1582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_APPROX((a * eig.pseudoEigenvectors()).norm()+1., 1.); 1592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_APPROX((eig.pseudoEigenvectors() * eig.pseudoEigenvalueMatrix()).norm()+1., 1.); 1602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_APPROX((a * eig.eigenvectors()).norm()+1., 1.); 1612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_APPROX((eig.eigenvectors() * eig.eigenvalues().asDiagonal()).norm()+1., 1.); 1622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 1632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#endif 1642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez TEST_SET_BUT_UNUSED_VARIABLE(s) 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 167