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 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (rows > 2) 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 } 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m) 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenSolver<MatrixType> eig; 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(eig.eigenvectors()); 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(eig.pseudoEigenvectors()); 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(eig.pseudoEigenvalueMatrix()); 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(eig.eigenvalues()); 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType a = MatrixType::Random(m.rows(),m.cols()); 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eig.compute(a, false); 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(eig.eigenvectors()); 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(eig.pseudoEigenvectors()); 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid test_eigensolver_generic() 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int s = 0; 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i = 0; i < g_repeat; i++) { 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_1( eigensolver(Matrix4f()) ); 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4); 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_2( eigensolver(MatrixXd(s,s)) ); 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // some trivial but implementation-wise tricky cases 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_2( eigensolver(MatrixXd(1,1)) ); 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_2( eigensolver(MatrixXd(2,2)) ); 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_3( eigensolver(Matrix<double,1,1>()) ); 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_4( eigensolver(Matrix2d()) ); 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_1( eigensolver_verify_assert(Matrix4f()) ); 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4); 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_2( eigensolver_verify_assert(MatrixXd(s,s)) ); 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_3( eigensolver_verify_assert(Matrix<double,1,1>()) ); 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_4( eigensolver_verify_assert(Matrix2d()) ); 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Test problem size constructors 1117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST_5(EigenSolver<MatrixXf> tmp(s)); 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // regression test for bug 410 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_2( 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixXd A(1,1); 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath A(0,0) = std::sqrt(-1.); 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Eigen::EigenSolver<MatrixXd> solver(A); 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixXd V(1, 1); 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath V(0,0) = solver.eigenvectors()(0,0).real(); 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ); 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez TEST_SET_BUT_UNUSED_VARIABLE(s) 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 126