diff -r dfb0ec8e0d86 Eigen/src/Eigenvalues/GeneralizedEigenSolver.h --- a/Eigen/src/Eigenvalues/GeneralizedEigenSolver.h Tue Aug 16 16:00:30 2016 -0700 +++ b/Eigen/src/Eigenvalues/GeneralizedEigenSolver.h Thu Aug 18 16:11:14 2016 +0100 @@ -89,7 +89,7 @@ */ typedef Matrix VectorType; - /** \brief Type for vector of complex scalar values eigenvalues as returned by betas(). + /** \brief Type for vector of complex scalar values eigenvalues as returned by alphas(). * * This is a column vector with entries of type #ComplexScalar. * The length of the vector is the size of #MatrixType. @@ -114,7 +114,14 @@ * * \sa compute() for an example. */ - GeneralizedEigenSolver() : m_eivec(), m_alphas(), m_betas(), m_isInitialized(false), m_realQZ(), m_matS(), m_tmp() {} + GeneralizedEigenSolver() + : m_eivec(), + m_alphas(), + m_betas(), + m_valuesOkay(false), + m_vectorsOkay(false), + m_realQZ() + {} /** \brief Default constructor with memory preallocation * @@ -126,11 +133,9 @@ : m_eivec(size, size), m_alphas(size), m_betas(size), - m_isInitialized(false), - m_eigenvectorsOk(false), - m_realQZ(size), - m_matS(size, size), - m_tmp(size) + m_valuesOkay(false), + m_vectorsOkay(false), + m_realQZ(size) {} /** \brief Constructor; computes the generalized eigendecomposition of given matrix pair. @@ -149,33 +154,29 @@ : m_eivec(A.rows(), A.cols()), m_alphas(A.cols()), m_betas(A.cols()), - m_isInitialized(false), - m_eigenvectorsOk(false), - m_realQZ(A.cols()), - m_matS(A.rows(), A.cols()), - m_tmp(A.cols()) + m_valuesOkay(false), + m_vectorsOkay(false), + m_realQZ(A.cols()) { compute(A, B, computeEigenvectors); } /* \brief Returns the computed generalized eigenvectors. * - * \returns %Matrix whose columns are the (possibly complex) eigenvectors. + * \returns %Matrix whose columns are the (possibly complex) right eigenvectors. + * i.e. the eigenvectors that solve (A - l*B)x = 0. The ordering matches the eigenvalues. * * \pre Either the constructor * GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function * compute(const MatrixType&, const MatrixType& bool) has been called before, and * \p computeEigenvectors was set to true (the default). * - * Column \f\$ k \f\$ of the returned matrix is an eigenvector corresponding - * to eigenvalue number \f\$ k \f\$ as returned by eigenvalues(). The - * eigenvectors are normalized to have (Euclidean) norm equal to one. The - * matrix returned by this function is the matrix \f\$ V \f\$ in the - * generalized eigendecomposition \f\$ A = B V D V^{-1} \f\$, if it exists. - * * \sa eigenvalues() */ -// EigenvectorsType eigenvectors() const; + EigenvectorsType eigenvectors() const { + eigen_assert(m_vectorsOkay && "Eigenvectors for GeneralizedEigenSolver were not calculated.") + return m_eivec; + } /** \brief Returns an expression of the computed generalized eigenvalues. * @@ -197,7 +198,7 @@ */ EigenvalueType eigenvalues() const { - eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized."); + eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized."); return EigenvalueType(m_alphas,m_betas); } @@ -208,7 +209,7 @@ * \sa betas(), eigenvalues() */ ComplexVectorType alphas() const { - eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized."); + eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized."); return m_alphas; } @@ -219,7 +220,7 @@ * \sa alphas(), eigenvalues() */ VectorType betas() const { - eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized."); + eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized."); return m_betas; } @@ -250,7 +251,7 @@ ComputationInfo info() const { - eigen_assert(m_isInitialized && "EigenSolver is not initialized."); + eigen_assert(m_valuesOkay && "EigenSolver is not initialized."); return m_realQZ.info(); } @@ -270,29 +271,13 @@ EIGEN_STATIC_ASSERT(!NumTraits::IsComplex, NUMERIC_TYPE_MUST_BE_REAL); } - MatrixType m_eivec; + EigenvectorsType m_eivec; ComplexVectorType m_alphas; VectorType m_betas; - bool m_isInitialized; - bool m_eigenvectorsOk; + bool m_valuesOkay, m_vectorsOkay; RealQZ m_realQZ; - MatrixType m_matS; - - typedef Matrix ColumnVectorType; - ColumnVectorType m_tmp; }; -//template -//typename GeneralizedEigenSolver::EigenvectorsType GeneralizedEigenSolver::eigenvectors() const -//{ -// eigen_assert(m_isInitialized && "EigenSolver is not initialized."); -// eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); -// Index n = m_eivec.cols(); -// EigenvectorsType matV(n,n); -// // TODO -// return matV; -//} - template GeneralizedEigenSolver& GeneralizedEigenSolver::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors) @@ -302,58 +287,111 @@ using std::sqrt; using std::abs; eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows()); - + m_valuesOkay = false; + m_vectorsOkay = false; // Reduce to generalized real Schur form: // A = Q S Z and B = Q T Z m_realQZ.compute(A, B, computeEigenvectors); - - if (m_realQZ.info() == Success) - { - m_matS = m_realQZ.matrixS(); - const MatrixType &matT = m_realQZ.matrixT(); - if (computeEigenvectors) - m_eivec = m_realQZ.matrixZ().transpose(); - - // Compute eigenvalues from matS + if (m_realQZ.info() == Success) { + // Temp space for the untransformed eigenvectors + VectorType v; + ComplexVectorType cv; + // Resize storage m_alphas.resize(A.cols()); m_betas.resize(A.cols()); + if (computeEigenvectors) { + m_eivec.resize(A.cols(), A.cols()); + v.resize(A.cols()); + cv.resize(A.cols()); + } + // Grab some references + const MatrixType &mZT = m_realQZ.matrixZ().transpose(); + const MatrixType &mS = m_realQZ.matrixS(); + const MatrixType &mT = m_realQZ.matrixT(); Index i = 0; - while (i < A.cols()) - { - if (i == A.cols() - 1 || m_matS.coeff(i+1, i) == Scalar(0)) - { - m_alphas.coeffRef(i) = m_matS.coeff(i, i); - m_betas.coeffRef(i) = matT.coeff(i,i); + while (i < A.cols()) { + if (i == A.cols() - 1 || mS.coeff(i+1, i) == Scalar(0)) { + // Real eigenvalue + m_alphas.coeffRef(i) = mS.diagonal().coeff(i); + m_betas.coeffRef(i) = mT.diagonal().coeff(i); + if (computeEigenvectors) { + v.setConstant(Scalar(0.0)); + v.coeffRef(i) = Scalar(1.0); + // For singular eigenvalues do nothing more + if (!internal::isApprox(fabs(m_betas.coeffRef(i)), 0.0)) { + //Non-singular eigenvalue + const Scalar alpha = real(m_alphas.coeffRef(i)); + const Scalar beta = m_betas.coeffRef(i); + for (Index j = i-1; j >= 0; j--) { + const Index st = j+1; + const Index sz = i-j; + if (j > 0 && mS.coeff(j, j-1) != Scalar(0)) { // 2x2 block + Matrix RHS; + RHS[0] = -v.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j-1,st,1,sz) - alpha*mT.block(j-1,st,1,sz)).sum(); + RHS[1] = -v.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum(); + Matrix LHS; + LHS << beta*mS.coeffRef(j-1,j-1) - alpha*mT.coeffRef(j-1,j-1), beta*mS.coeffRef(j-1,j) - alpha*mT.coeffRef(j-1,j), + beta*mS.coeffRef(j,j-1) - alpha*mT.coeffRef(j,j-1), beta*mS.coeffRef(j,j) - alpha*mT.coeffRef(j,j); + v.segment(j-1,2) = LHS.partialPivLu().solve(RHS); + j--; + } else { + v.coeffRef(j) = -v.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum() / (beta*mS.coeffRef(j,j) - alpha*mT.coeffRef(j,j)); + } + } + } + m_eivec.col(i).real() = (mZT * v).normalized(); + m_eivec.col(i).imag().setConstant(0); + } ++i; - } - else - { + } else { // We need to extract the generalized eigenvalues of the pair of a general 2x2 block S and a positive diagonal 2x2 block T // Then taking beta=T_00*T_11, we can avoid any division, and alpha is the eigenvalues of A = (U^-1 * S * U) * diag(T_11,T_00): // T = [a 0] // [0 b] - RealScalar a = matT.diagonal().coeff(i), - b = matT.diagonal().coeff(i+1); + RealScalar a = mT.diagonal().coeff(i), + b = mT.diagonal().coeff(i+1); + const RealScalar beta = m_betas.coeffRef(i) = m_betas.coeffRef(i+1) = a*b; + // ^^ NOTE: using diagonal()(i) instead of coeff(i,i) workarounds a MSVC bug. - Matrix S2 = m_matS.template block<2,2>(i,i) * Matrix(b,a).asDiagonal(); + Matrix S2 = mS.template block<2,2>(i,i) * Matrix(b,a).asDiagonal(); Scalar p = Scalar(0.5) * (S2.coeff(0,0) - S2.coeff(1,1)); Scalar z = sqrt(abs(p * p + S2.coeff(1,0) * S2.coeff(0,1))); - m_alphas.coeffRef(i) = ComplexScalar(S2.coeff(1,1) + p, z); - m_alphas.coeffRef(i+1) = ComplexScalar(S2.coeff(1,1) + p, -z); + const ComplexScalar alpha = ComplexScalar(S2.coeff(1,1) + p, (beta > 0) ? z : -z); + m_alphas.coeffRef(i) = conj(alpha); + m_alphas.coeffRef(i+1) = alpha; - m_betas.coeffRef(i) = - m_betas.coeffRef(i+1) = a*b; - + if (computeEigenvectors) { + // Compute eigenvector in position (i+1) and then position (i) is just the conjugate + cv.setConstant(ComplexScalar(0.0, 0.0)); + cv.coeffRef(i+1) = ComplexScalar(1.0, 0.0); + cv.coeffRef(i) = -(beta*mS.coeffRef(i,i+1) - alpha*mT.coeffRef(i,i+1)) / (beta*mS.coeffRef(i,i) - alpha*mT.coeffRef(i,i)); + for (Index j = i-1; j >= 0; j--) { + const Index st = j+1; + const Index sz = i+1-j; + if (j > 0 && mS.coeff(j, j-1) != Scalar(0)) { // 2x2 block + Matrix RHS; + RHS[0] = -cv.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j-1,st,1,sz) - alpha*mT.block(j-1,st,1,sz)).sum(); + RHS[1] = -cv.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum(); + Matrix LHS; + LHS << beta*mS.coeffRef(j-1,j-1) - alpha*mT.coeffRef(j-1,j-1), beta*mS.coeffRef(j-1,j) - alpha*mT.coeffRef(j-1,j), + beta*mS.coeffRef(j,j-1) - alpha*mT.coeffRef(j,j-1), beta*mS.coeffRef(j,j) - alpha*mT.coeffRef(j,j); + cv.segment(j-1,2) = LHS.partialPivLu().solve(RHS); + j--; + } else { + cv.coeffRef(j) = -cv.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum() / (beta*mS.coeffRef(j,j) - alpha*mT.coeffRef(j,j)); + } + } + m_eivec.col(i+1) = (mZT * cv).normalized(); + m_eivec.col(i) = m_eivec.col(i+1).conjugate(); + } i += 2; } } + m_valuesOkay = true; + m_vectorsOkay = computeEigenvectors; } - - m_isInitialized = true; - m_eigenvectorsOk = false;//computeEigenvectors; - return *this; } diff -r dfb0ec8e0d86 test/CMakeLists.txt --- a/test/CMakeLists.txt Tue Aug 16 16:00:30 2016 -0700 +++ b/test/CMakeLists.txt Thu Aug 18 16:11:14 2016 +0100 @@ -211,6 +211,7 @@ ei_add_test(eigensolver_complex) ei_add_test(real_qz) ei_add_test(eigensolver_generalized_real) +ei_add_test(eigensolver_generalized_complex) ei_add_test(jacobi) ei_add_test(jacobisvd) ei_add_test(bdcsvd) diff -r dfb0ec8e0d86 test/eigensolver_generalized_complex.cpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/eigensolver_generalized_complex.cpp Thu Aug 18 16:11:14 2016 +0100 @@ -0,0 +1,45 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Gael Guennebaud +// +// 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 +#include + +template void generalized_eigensolver_complex(const MatrixType& m) +{ + typedef typename MatrixType::Index Index; + /* this test covers the following files: + GeneralizedEigenSolver.h + */ + Index rows = m.rows(); + Index cols = m.cols(); + MatrixType a = MatrixType::Random(rows,cols); + MatrixType b = MatrixType::Random(rows,cols); + GeneralizedEigenSolver eig(a, b); + typename GeneralizedEigenSolver::EigenvectorsType D = eig.eigenvalues().asDiagonal(); + typename GeneralizedEigenSolver::EigenvectorsType V = eig.eigenvectors(); + VERIFY_IS_APPROX(a*V, b*V*D); +} + +void test_eigensolver_generalized_complex() +{ + for(int i = 0; i < g_repeat; i++) { + int s = 0; + CALL_SUBTEST_1( generalized_eigensolver_complex(Matrix4f()) ); + s = internal::random(1,EIGEN_TEST_MAX_SIZE/4); + CALL_SUBTEST_2( generalized_eigensolver_complex(MatrixXd(s,s)) ); + + // some trivial but implementation-wise tricky cases + CALL_SUBTEST_3( generalized_eigensolver_complex(MatrixXd(1,1)) ); + CALL_SUBTEST_4( generalized_eigensolver_complex(MatrixXd(2,2)) ); + CALL_SUBTEST_5( generalized_eigensolver_complex(Matrix()) ); + CALL_SUBTEST_6( generalized_eigensolver_complex(Matrix2d()) ); + TEST_SET_BUT_UNUSED_VARIABLE(s) + } +}