mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-12-21 07:19:46 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add unit test for non symmetric generalized eigenvalues

Browse Source
This commit is contained in:
Gael Guennebaud 2016-06-09 16:17:27 +02:00
parent a20d2ec1c0
commit 15890c304e
1 changed files with 24 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
test/eigensolver_generalized_real.cpp
View File

@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library // This file is part of Eigen, a lightweight C++ template library
// for linear algebra. // for linear algebra.
// //
// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr> // Copyright (C) 2012-2016 Gael Guennebaud <gael.guennebaud@inria.fr>
// //
// This Source Code Form is subject to the terms of the Mozilla // 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 // Public License v. 2.0. If a copy of the MPL was not distributed
@ -10,6 +10,7 @@
#include "main.h" #include "main.h"
#include <limits> #include <limits>
#include <Eigen/Eigenvalues> #include <Eigen/Eigenvalues>
#include <Eigen/LU>
template<typename MatrixType> void generalized_eigensolver_real(const MatrixType& m) template<typename MatrixType> void generalized_eigensolver_real(const MatrixType& m)
{ {
@ -21,6 +22,7 @@ template<typename MatrixType> void generalized_eigensolver_real(const MatrixType
Index cols = m.cols(); Index cols = m.cols();
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef std::complex<Scalar> ComplexScalar;
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
MatrixType a = MatrixType::Random(rows,cols); MatrixType a = MatrixType::Random(rows,cols);
@ -31,14 +33,28 @@ template<typename MatrixType> void generalized_eigensolver_real(const MatrixType
MatrixType spdB = b.adjoint() * b + b1.adjoint() * b1; MatrixType spdB = b.adjoint() * b + b1.adjoint() * b1;
// lets compare to GeneralizedSelfAdjointEigenSolver // lets compare to GeneralizedSelfAdjointEigenSolver
GeneralizedSelfAdjointEigenSolver<MatrixType> symmEig(spdA, spdB); {
GeneralizedEigenSolver<MatrixType> eig(spdA, spdB); GeneralizedSelfAdjointEigenSolver<MatrixType> symmEig(spdA, spdB);
GeneralizedEigenSolver<MatrixType> eig(spdA, spdB);
VERIFY_IS_EQUAL(eig.eigenvalues().imag().cwiseAbs().maxCoeff(), 0); VERIFY_IS_EQUAL(eig.eigenvalues().imag().cwiseAbs().maxCoeff(), 0);
VectorType realEigenvalues = eig.eigenvalues().real(); VectorType realEigenvalues = eig.eigenvalues().real();
std::sort(realEigenvalues.data(), realEigenvalues.data()+realEigenvalues.size()); std::sort(realEigenvalues.data(), realEigenvalues.data()+realEigenvalues.size());
VERIFY_IS_APPROX(realEigenvalues, symmEig.eigenvalues()); VERIFY_IS_APPROX(realEigenvalues, symmEig.eigenvalues());
}
// non symmetric case:
{
GeneralizedEigenSolver<MatrixType> eig(a,b);
for(Index k=0; k<cols; ++k)
{
Matrix<ComplexScalar,Dynamic,Dynamic> tmp = (eig.betas()(k)*a).template cast<ComplexScalar>() - eig.alphas()(k)*b;
if(tmp.norm()>(std::numeric_limits<Scalar>::min)())
tmp /= tmp.norm();
VERIFY_IS_MUCH_SMALLER_THAN( std::abs(tmp.determinant()), Scalar(1) );
}
}
// regression test for bug 1098 // regression test for bug 1098
{ {
@ -57,7 +73,7 @@ void test_eigensolver_generalized_real()
s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4); s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4);
CALL_SUBTEST_2( generalized_eigensolver_real(MatrixXd(s,s)) ); CALL_SUBTEST_2( generalized_eigensolver_real(MatrixXd(s,s)) );
// some trivial but implementation-wise tricky cases // some trivial but implementation-wise special cases
CALL_SUBTEST_2( generalized_eigensolver_real(MatrixXd(1,1)) ); CALL_SUBTEST_2( generalized_eigensolver_real(MatrixXd(1,1)) );
CALL_SUBTEST_2( generalized_eigensolver_real(MatrixXd(2,2)) ); CALL_SUBTEST_2( generalized_eigensolver_real(MatrixXd(2,2)) );
CALL_SUBTEST_3( generalized_eigensolver_real(Matrix<double,1,1>()) ); CALL_SUBTEST_3( generalized_eigensolver_real(Matrix<double,1,1>()) );