From 68b28f7bfb29474ad21036476618a3730fa7fffa Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 4 Sep 2009 09:23:38 +0200 Subject: rename the EigenSolver module to Eigenvalues --- Eigen/CMakeLists.txt | 2 +- Eigen/EigenSolver | 74 --- Eigen/Eigenvalues | 74 +++ Eigen/LeastSquares | 2 +- Eigen/QR | 4 +- Eigen/src/CMakeLists.txt | 2 +- Eigen/src/EigenSolver/CMakeLists.txt | 6 - Eigen/src/EigenSolver/ComplexEigenSolver.h | 139 ----- Eigen/src/EigenSolver/ComplexSchur.h | 238 -------- Eigen/src/EigenSolver/EigenSolver.h | 723 ------------------------ Eigen/src/EigenSolver/HessenbergDecomposition.h | 206 ------- Eigen/src/EigenSolver/SelfAdjointEigenSolver.h | 366 ------------ Eigen/src/EigenSolver/Tridiagonalization.h | 317 ----------- Eigen/src/Eigenvalues/CMakeLists.txt | 6 + Eigen/src/Eigenvalues/ComplexEigenSolver.h | 148 +++++ Eigen/src/Eigenvalues/ComplexSchur.h | 237 ++++++++ Eigen/src/Eigenvalues/EigenSolver.h | 723 ++++++++++++++++++++++++ Eigen/src/Eigenvalues/HessenbergDecomposition.h | 206 +++++++ Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h | 366 ++++++++++++ Eigen/src/Eigenvalues/Tridiagonalization.h | 317 +++++++++++ 20 files changed, 2082 insertions(+), 2074 deletions(-) delete mode 100644 Eigen/EigenSolver create mode 100644 Eigen/Eigenvalues delete mode 100644 Eigen/src/EigenSolver/CMakeLists.txt delete mode 100644 Eigen/src/EigenSolver/ComplexEigenSolver.h delete mode 100644 Eigen/src/EigenSolver/ComplexSchur.h delete mode 100644 Eigen/src/EigenSolver/EigenSolver.h delete mode 100644 Eigen/src/EigenSolver/HessenbergDecomposition.h delete mode 100644 Eigen/src/EigenSolver/SelfAdjointEigenSolver.h delete mode 100644 Eigen/src/EigenSolver/Tridiagonalization.h create mode 100644 Eigen/src/Eigenvalues/CMakeLists.txt create mode 100644 Eigen/src/Eigenvalues/ComplexEigenSolver.h create mode 100644 Eigen/src/Eigenvalues/ComplexSchur.h create mode 100644 Eigen/src/Eigenvalues/EigenSolver.h create mode 100644 Eigen/src/Eigenvalues/HessenbergDecomposition.h create mode 100644 Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h create mode 100644 Eigen/src/Eigenvalues/Tridiagonalization.h (limited to 'Eigen') diff --git a/Eigen/CMakeLists.txt b/Eigen/CMakeLists.txt index ebdf57812..931cc6e20 100644 --- a/Eigen/CMakeLists.txt +++ b/Eigen/CMakeLists.txt @@ -1,4 +1,4 @@ -set(Eigen_HEADERS Core LU Cholesky QR Geometry Sparse Array SVD LeastSquares QtAlignedMalloc StdVector Householder Jacobi EigenSolver) +set(Eigen_HEADERS Core LU Cholesky QR Geometry Sparse Array SVD LeastSquares QtAlignedMalloc StdVector Householder Jacobi Eigenvalues) if(EIGEN_BUILD_LIB) set(Eigen_SRCS diff --git a/Eigen/EigenSolver b/Eigen/EigenSolver deleted file mode 100644 index fd126d282..000000000 --- a/Eigen/EigenSolver +++ /dev/null @@ -1,74 +0,0 @@ -#ifndef EIGEN_EIGEN_SOLVER_MODULE_H -#define EIGEN_EIGEN_SOLVER_MODULE_H - -#include "Core" - -#include "src/Core/util/DisableMSVCWarnings.h" - -#include "Cholesky" -#include "Jacobi" -#include "Householder" - -// Note that EIGEN_HIDE_HEAVY_CODE has to be defined per module -#if (defined EIGEN_EXTERN_INSTANTIATIONS) && (EIGEN_EXTERN_INSTANTIATIONS>=2) - #ifndef EIGEN_HIDE_HEAVY_CODE - #define EIGEN_HIDE_HEAVY_CODE - #endif -#elif defined EIGEN_HIDE_HEAVY_CODE - #undef EIGEN_HIDE_HEAVY_CODE -#endif - -namespace Eigen { - -/** \defgroup EigenSolver_Module Eigen Solver module - * - * \nonstableyet - * - * This module mainly provides various eigen value solvers. - * This module also provides some MatrixBase methods, including: - * - MatrixBase::eigenvalues(), - * - MatrixBase::operatorNorm() - * - * \code - * #include - * \endcode - */ - -#include "src/EigenSolver/Tridiagonalization.h" -#include "src/EigenSolver/EigenSolver.h" -#include "src/EigenSolver/SelfAdjointEigenSolver.h" -#include "src/EigenSolver/HessenbergDecomposition.h" -#include "src/EigenSolver/ComplexSchur.h" -#include "src/EigenSolver/ComplexEigenSolver.h" - -// declare all classes for a given matrix type -#define EIGEN_EIGEN_SOLVER_MODULE_INSTANTIATE_TYPE(MATRIXTYPE,PREFIX) \ - PREFIX template class Tridiagonalization; \ - PREFIX template class HessenbergDecomposition; \ - PREFIX template class SelfAdjointEigenSolver - -// removed because it does not support complex yet -// PREFIX template class EigenSolver - -// declare all class for all types -#define EIGEN_EIGEN_SOLVER_MODULE_INSTANTIATE(PREFIX) \ - EIGEN_EIGEN_SOLVER_MODULE_INSTANTIATE_TYPE(Matrix2f,PREFIX); \ - EIGEN_EIGEN_SOLVER_MODULE_INSTANTIATE_TYPE(Matrix2d,PREFIX); \ - EIGEN_EIGEN_SOLVER_MODULE_INSTANTIATE_TYPE(Matrix3f,PREFIX); \ - EIGEN_EIGEN_SOLVER_MODULE_INSTANTIATE_TYPE(Matrix3d,PREFIX); \ - EIGEN_EIGEN_SOLVER_MODULE_INSTANTIATE_TYPE(Matrix4f,PREFIX); \ - EIGEN_EIGEN_SOLVER_MODULE_INSTANTIATE_TYPE(Matrix4d,PREFIX); \ - EIGEN_EIGEN_SOLVER_MODULE_INSTANTIATE_TYPE(MatrixXf,PREFIX); \ - EIGEN_EIGEN_SOLVER_MODULE_INSTANTIATE_TYPE(MatrixXd,PREFIX); \ - EIGEN_EIGEN_SOLVER_MODULE_INSTANTIATE_TYPE(MatrixXcf,PREFIX); \ - EIGEN_EIGEN_SOLVER_MODULE_INSTANTIATE_TYPE(MatrixXcd,PREFIX) - -#ifdef EIGEN_EXTERN_INSTANTIATIONS - EIGEN_EIGEN_SOLVER_MODULE_INSTANTIATE(extern); -#endif // EIGEN_EXTERN_INSTANTIATIONS - -} // namespace Eigen - -#include "src/Core/util/EnableMSVCWarnings.h" - -#endif // EIGEN_EIGEN_SOLVER_MODULE_H diff --git a/Eigen/Eigenvalues b/Eigen/Eigenvalues new file mode 100644 index 000000000..9a6443f39 --- /dev/null +++ b/Eigen/Eigenvalues @@ -0,0 +1,74 @@ +#ifndef EIGEN_EIGENVALUES_MODULE_H +#define EIGEN_EIGENVALUES_MODULE_H + +#include "Core" + +#include "src/Core/util/DisableMSVCWarnings.h" + +#include "Cholesky" +#include "Jacobi" +#include "Householder" + +// Note that EIGEN_HIDE_HEAVY_CODE has to be defined per module +#if (defined EIGEN_EXTERN_INSTANTIATIONS) && (EIGEN_EXTERN_INSTANTIATIONS>=2) + #ifndef EIGEN_HIDE_HEAVY_CODE + #define EIGEN_HIDE_HEAVY_CODE + #endif +#elif defined EIGEN_HIDE_HEAVY_CODE + #undef EIGEN_HIDE_HEAVY_CODE +#endif + +namespace Eigen { + +/** \defgroup Eigenvalues_Module Eigenvalues module + * + * \nonstableyet + * + * This module mainly provides various eigenvalue solvers. + * This module also provides some MatrixBase methods, including: + * - MatrixBase::eigenvalues(), + * - MatrixBase::operatorNorm() + * + * \code + * #include + * \endcode + */ + +#include "src/Eigenvalues/Tridiagonalization.h" +#include "src/Eigenvalues/EigenSolver.h" +#include "src/Eigenvalues/SelfAdjointEigenSolver.h" +#include "src/Eigenvalues/HessenbergDecomposition.h" +#include "src/Eigenvalues/ComplexSchur.h" +#include "src/Eigenvalues/ComplexEigenSolver.h" + +// declare all classes for a given matrix type +#define EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(MATRIXTYPE,PREFIX) \ + PREFIX template class Tridiagonalization; \ + PREFIX template class HessenbergDecomposition; \ + PREFIX template class SelfAdjointEigenSolver + +// removed because it does not support complex yet +// PREFIX template class EigenSolver + +// declare all class for all types +#define EIGEN_EIGENVALUES_MODULE_INSTANTIATE(PREFIX) \ + EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(Matrix2f,PREFIX); \ + EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(Matrix2d,PREFIX); \ + EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(Matrix3f,PREFIX); \ + EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(Matrix3d,PREFIX); \ + EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(Matrix4f,PREFIX); \ + EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(Matrix4d,PREFIX); \ + EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(MatrixXf,PREFIX); \ + EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(MatrixXd,PREFIX); \ + EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(MatrixXcf,PREFIX); \ + EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(MatrixXcd,PREFIX) + +#ifdef EIGEN_EXTERN_INSTANTIATIONS + EIGEN_EIGENVALUES_MODULE_INSTANTIATE(extern); +#endif // EIGEN_EXTERN_INSTANTIATIONS + +} // namespace Eigen + +#include "src/Core/util/EnableMSVCWarnings.h" + +#endif // EIGEN_EIGENVALUES_MODULE_H diff --git a/Eigen/LeastSquares b/Eigen/LeastSquares index db620f548..75620a349 100644 --- a/Eigen/LeastSquares +++ b/Eigen/LeastSquares @@ -5,7 +5,7 @@ #include "src/Core/util/DisableMSVCWarnings.h" -#include "EigenSolver" +#include "Eigenvalues" #include "Geometry" namespace Eigen { diff --git a/Eigen/QR b/Eigen/QR index a7273bc8a..f38e96c5e 100644 --- a/Eigen/QR +++ b/Eigen/QR @@ -62,7 +62,7 @@ namespace Eigen { #include "src/Core/util/EnableMSVCWarnings.h" -// FIXME for compatibility we include EigenSolver here: -#include "EigenSolver" +// FIXME for compatibility we include Eigenvalues here: +#include "Eigenvalues" #endif // EIGEN_QR_MODULE_H diff --git a/Eigen/src/CMakeLists.txt b/Eigen/src/CMakeLists.txt index 6b1333860..0df8273d1 100644 --- a/Eigen/src/CMakeLists.txt +++ b/Eigen/src/CMakeLists.txt @@ -9,4 +9,4 @@ ADD_SUBDIRECTORY(LeastSquares) ADD_SUBDIRECTORY(Sparse) ADD_SUBDIRECTORY(Jacobi) ADD_SUBDIRECTORY(Householder) -ADD_SUBDIRECTORY(EigenSolver) +ADD_SUBDIRECTORY(Eigenvalues) diff --git a/Eigen/src/EigenSolver/CMakeLists.txt b/Eigen/src/EigenSolver/CMakeLists.txt deleted file mode 100644 index 63bc75e0c..000000000 --- a/Eigen/src/EigenSolver/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_EIGENSOLVER_SRCS "*.h") - -INSTALL(FILES - ${Eigen_EIGENSOLVER_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/EigenSolver COMPONENT Devel - ) diff --git a/Eigen/src/EigenSolver/ComplexEigenSolver.h b/Eigen/src/EigenSolver/ComplexEigenSolver.h deleted file mode 100644 index 2ea7464a6..000000000 --- a/Eigen/src/EigenSolver/ComplexEigenSolver.h +++ /dev/null @@ -1,139 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2009 Claire Maurice -// Copyright (C) 2009 Gael Guennebaud -// -// Eigen is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 3 of the License, or (at your option) any later version. -// -// Alternatively, you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of -// the License, or (at your option) any later version. -// -// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License and a copy of the GNU General Public License along with -// Eigen. If not, see . - -#ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H -#define EIGEN_COMPLEX_EIGEN_SOLVER_H - -#define MAXITER 30 - -template class ComplexEigenSolver -{ - public: - typedef _MatrixType MatrixType; - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; - typedef std::complex Complex; - typedef Matrix EigenvalueType; - typedef Matrix EigenvectorType; - - /** - * \brief Default Constructor. - * - * The default constructor is useful in cases in which the user intends to - * perform decompositions via ComplexEigenSolver::compute(const MatrixType&). - */ - ComplexEigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false) - {} - - ComplexEigenSolver(const MatrixType& matrix) - : m_eivec(matrix.rows(),matrix.cols()), - m_eivalues(matrix.cols()), - m_isInitialized(false) - { - compute(matrix); - } - - EigenvectorType eigenvectors(void) const - { - ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); - return m_eivec; - } - - EigenvalueType eigenvalues() const - { - ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); - return m_eivalues; - } - - void compute(const MatrixType& matrix); - - protected: - MatrixType m_eivec; - EigenvalueType m_eivalues; - bool m_isInitialized; -}; - - -template -void ComplexEigenSolver::compute(const MatrixType& matrix) -{ - // this code is inspired from Jampack - assert(matrix.cols() == matrix.rows()); - int n = matrix.cols(); - m_eivalues.resize(n,1); - - RealScalar eps = epsilon(); - - // Reduce to complex Schur form - ComplexSchur schur(matrix); - - m_eivalues = schur.matrixT().diagonal(); - - m_eivec.setZero(); - - Scalar d2, z; - RealScalar norm = matrix.norm(); - - // compute the (normalized) eigenvectors - for(int k=n-1 ; k>=0 ; k--) - { - d2 = schur.matrixT().coeff(k,k); - m_eivec.coeffRef(k,k) = Scalar(1.0,0.0); - for(int i=k-1 ; i>=0 ; i--) - { - m_eivec.coeffRef(i,k) = -schur.matrixT().coeff(i,k); - if(k-i-1>0) - m_eivec.coeffRef(i,k) -= (schur.matrixT().row(i).segment(i+1,k-i-1) * m_eivec.col(k).segment(i+1,k-i-1)).value(); - z = schur.matrixT().coeff(i,i) - d2; - if(z==Scalar(0)) - ei_real_ref(z) = eps * norm; - m_eivec.coeffRef(i,k) = m_eivec.coeff(i,k) / z; - - } - m_eivec.col(k).normalize(); - } - - m_eivec = schur.matrixU() * m_eivec; - m_isInitialized = true; - - // sort the eigenvalues - { - for (int i=0; i -// -// Eigen is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 3 of the License, or (at your option) any later version. -// -// Alternatively, you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of -// the License, or (at your option) any later version. -// -// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License and a copy of the GNU General Public License along with -// Eigen. If not, see . - -#ifndef EIGEN_COMPLEX_SCHUR_H -#define EIGEN_COMPLEX_SCHUR_H - -#define MAXITER 30 - -/** \ingroup QR - * - * \class ComplexShur - * - * \brief Performs a complex Shur decomposition of a real or complex square matrix - * - */ -template class ComplexSchur -{ - public: - typedef _MatrixType MatrixType; - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; - typedef std::complex Complex; - typedef Matrix ComplexMatrixType; - - /** - * \brief Default Constructor. - * - * The default constructor is useful in cases in which the user intends to - * perform decompositions via ComplexSchur::compute(const MatrixType&). - */ - ComplexSchur() : m_matT(), m_matU(), m_isInitialized(false) - {} - - ComplexSchur(const MatrixType& matrix) - : m_matT(matrix.rows(),matrix.cols()), - m_matU(matrix.rows(),matrix.cols()), - m_isInitialized(false) - { - compute(matrix); - } - - ComplexMatrixType matrixU() const - { - ei_assert(m_isInitialized && "ComplexSchur is not initialized."); - return m_matU; - } - - ComplexMatrixType matrixT() const - { - ei_assert(m_isInitialized && "ComplexShur is not initialized."); - return m_matT; - } - - void compute(const MatrixType& matrix); - - protected: - ComplexMatrixType m_matT, m_matU; - bool m_isInitialized; -}; - -/** Computes the principal value of the square root of the complex \a z. */ -template -std::complex ei_sqrt(const std::complex &z) -{ - RealScalar t, tre, tim; - - t = ei_abs(z); - - if (ei_abs(ei_real(z)) <= ei_abs(ei_imag(z))) - { - // No cancellation in these formulas - tre = ei_sqrt(0.5*(t + ei_real(z))); - tim = ei_sqrt(0.5*(t - ei_real(z))); - } - else - { - // Stable computation of the above formulas - if (z.real() > 0) - { - tre = t + z.real(); - tim = ei_abs(ei_imag(z))*ei_sqrt(0.5/tre); - tre = ei_sqrt(0.5*tre); - } - else - { - tim = t - z.real(); - tre = ei_abs(ei_imag(z))*ei_sqrt(0.5/tim); - tim = ei_sqrt(0.5*tim); - } - } - if(z.imag() < 0) - tim = -tim; - - return (std::complex(tre,tim)); - -} - -template -void ComplexSchur::compute(const MatrixType& matrix) -{ - // this code is inspired from Jampack - assert(matrix.cols() == matrix.rows()); - int n = matrix.cols(); - - // Reduce to Hessenberg form - HessenbergDecomposition hess(matrix); - - m_matT = hess.matrixH(); - m_matU = hess.matrixQ(); - - int iu = m_matT.cols() - 1; - int il; - RealScalar d,sd,sf; - Complex c,b,disc,r1,r2,kappa; - - RealScalar eps = epsilon(); - - int iter = 0; - while(true) - { - //locate the range in which to iterate - while(iu > 0) - { - d = ei_norm1(m_matT.coeffRef(iu,iu)) + ei_norm1(m_matT.coeffRef(iu-1,iu-1)); - sd = ei_norm1(m_matT.coeffRef(iu,iu-1)); - - if(sd >= eps * d) break; // FIXME : precision criterion ?? - - m_matT.coeffRef(iu,iu-1) = Complex(0); - iter = 0; - --iu; - } - if(iu==0) break; - iter++; - - if(iter >= MAXITER) - { - // FIXME : what to do when iter==MAXITER ?? - std::cerr << "MAXITER" << std::endl; - return; - } - - il = iu-1; - while( il > 0 ) - { - // check if the current 2x2 block on the diagonal is upper triangular - d = ei_norm1(m_matT.coeffRef(il,il)) + ei_norm1(m_matT.coeffRef(il-1,il-1)); - sd = ei_norm1(m_matT.coeffRef(il,il-1)); - - if(sd < eps * d) break; // FIXME : precision criterion ?? - - --il; - } - - if( il != 0 ) m_matT.coeffRef(il,il-1) = Complex(0); - - // compute the shift (the normalization by sf is to avoid under/overflow) - Matrix t = m_matT.template block<2,2>(iu-1,iu-1); - sf = t.cwise().abs().sum(); - t /= sf; - - c = t.determinant(); - b = t.diagonal().sum(); - - disc = ei_sqrt(b*b - RealScalar(4)*c); - - r1 = (b+disc)/RealScalar(2); - r2 = (b-disc)/RealScalar(2); - - if(ei_norm1(r1) > ei_norm1(r2)) - r2 = c/r1; - else - r1 = c/r2; - - if(ei_norm1(r1-t.coeff(1,1)) < ei_norm1(r2-t.coeff(1,1))) - kappa = sf * r1; - else - kappa = sf * r2; - - // perform the QR step using Givens rotations - PlanarRotation rot; - rot.makeGivens(m_matT.coeff(il,il) - kappa, m_matT.coeff(il+1,il)); - - for(int i=il ; i -// -// Eigen is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 3 of the License, or (at your option) any later version. -// -// Alternatively, you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of -// the License, or (at your option) any later version. -// -// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License and a copy of the GNU General Public License along with -// Eigen. If not, see . - -#ifndef EIGEN_EIGENSOLVER_H -#define EIGEN_EIGENSOLVER_H - -/** \ingroup EigenSolver_Module - * \nonstableyet - * - * \class EigenSolver - * - * \brief Eigen values/vectors solver for non selfadjoint matrices - * - * \param MatrixType the type of the matrix of which we are computing the eigen decomposition - * - * Currently it only support real matrices. - * - * \note this code was adapted from JAMA (public domain) - * - * \sa MatrixBase::eigenvalues(), SelfAdjointEigenSolver - */ -template class EigenSolver -{ - public: - - typedef _MatrixType MatrixType; - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; - typedef std::complex Complex; - typedef Matrix EigenvalueType; - typedef Matrix EigenvectorType; - typedef Matrix RealVectorType; - typedef Matrix RealVectorTypeX; - - /** - * \brief Default Constructor. - * - * The default constructor is useful in cases in which the user intends to - * perform decompositions via EigenSolver::compute(const MatrixType&). - */ - EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false) {} - - EigenSolver(const MatrixType& matrix) - : m_eivec(matrix.rows(), matrix.cols()), - m_eivalues(matrix.cols()), - m_isInitialized(false) - { - compute(matrix); - } - - - EigenvectorType eigenvectors(void) const; - - /** \returns a real matrix V of pseudo eigenvectors. - * - * Let D be the block diagonal matrix with the real eigenvalues in 1x1 blocks, - * and any complex values u+iv in 2x2 blocks [u v ; -v u]. Then, the matrices D - * and V satisfy A*V = V*D. - * - * More precisely, if the diagonal matrix of the eigen values is:\n - * \f$ - * \left[ \begin{array}{cccccc} - * u+iv & & & & & \\ - * & u-iv & & & & \\ - * & & a+ib & & & \\ - * & & & a-ib & & \\ - * & & & & x & \\ - * & & & & & y \\ - * \end{array} \right] - * \f$ \n - * then, we have:\n - * \f$ - * D =\left[ \begin{array}{cccccc} - * u & v & & & & \\ - * -v & u & & & & \\ - * & & a & b & & \\ - * & & -b & a & & \\ - * & & & & x & \\ - * & & & & & y \\ - * \end{array} \right] - * \f$ - * - * \sa pseudoEigenvalueMatrix() - */ - const MatrixType& pseudoEigenvectors() const - { - ei_assert(m_isInitialized && "EigenSolver is not initialized."); - return m_eivec; - } - - MatrixType pseudoEigenvalueMatrix() const; - - /** \returns the eigenvalues as a column vector */ - EigenvalueType eigenvalues() const - { - ei_assert(m_isInitialized && "EigenSolver is not initialized."); - return m_eivalues; - } - - EigenSolver& compute(const MatrixType& matrix); - - private: - - void orthes(MatrixType& matH, RealVectorType& ort); - void hqr2(MatrixType& matH); - - protected: - MatrixType m_eivec; - EigenvalueType m_eivalues; - bool m_isInitialized; -}; - -/** \returns the real block diagonal matrix D of the eigenvalues. - * - * See pseudoEigenvectors() for the details. - */ -template -MatrixType EigenSolver::pseudoEigenvalueMatrix() const -{ - ei_assert(m_isInitialized && "EigenSolver is not initialized."); - int n = m_eivec.cols(); - MatrixType matD = MatrixType::Zero(n,n); - for (int i=0; i(i,i) << ei_real(m_eivalues.coeff(i)), ei_imag(m_eivalues.coeff(i)), - -ei_imag(m_eivalues.coeff(i)), ei_real(m_eivalues.coeff(i)); - ++i; - } - } - return matD; -} - -/** \returns the normalized complex eigenvectors as a matrix of column vectors. - * - * \sa eigenvalues(), pseudoEigenvectors() - */ -template -typename EigenSolver::EigenvectorType EigenSolver::eigenvectors(void) const -{ - ei_assert(m_isInitialized && "EigenSolver is not initialized."); - int n = m_eivec.cols(); - EigenvectorType matV(n,n); - for (int j=0; j(); - } - else - { - // we have a pair of complex eigen values - for (int i=0; i -EigenSolver& EigenSolver::compute(const MatrixType& matrix) -{ - assert(matrix.cols() == matrix.rows()); - int n = matrix.cols(); - m_eivalues.resize(n,1); - - MatrixType matH = matrix; - RealVectorType ort(n); - - // Reduce to Hessenberg form. - orthes(matH, ort); - - // Reduce Hessenberg to real Schur form. - hqr2(matH); - - m_isInitialized = true; - return *this; -} - -// Nonsymmetric reduction to Hessenberg form. -template -void EigenSolver::orthes(MatrixType& matH, RealVectorType& ort) -{ - // This is derived from the Algol procedures orthes and ortran, - // by Martin and Wilkinson, Handbook for Auto. Comp., - // Vol.ii-Linear Algebra, and the corresponding - // Fortran subroutines in EISPACK. - - int n = m_eivec.cols(); - int low = 0; - int high = n-1; - - for (int m = low+1; m <= high-1; ++m) - { - // Scale column. - RealScalar scale = matH.block(m, m-1, high-m+1, 1).cwise().abs().sum(); - if (scale != 0.0) - { - // Compute Householder transformation. - RealScalar h = 0.0; - // FIXME could be rewritten, but this one looks better wrt cache - for (int i = high; i >= m; i--) - { - ort.coeffRef(i) = matH.coeff(i,m-1)/scale; - h += ort.coeff(i) * ort.coeff(i); - } - RealScalar g = ei_sqrt(h); - if (ort.coeff(m) > 0) - g = -g; - h = h - ort.coeff(m) * g; - ort.coeffRef(m) = ort.coeff(m) - g; - - // Apply Householder similarity transformation - // H = (I-u*u'/h)*H*(I-u*u')/h) - int bSize = high-m+1; - matH.block(m, m, bSize, n-m).noalias() -= ((ort.segment(m, bSize)/h) - * (ort.segment(m, bSize).transpose() * matH.block(m, m, bSize, n-m))); - - matH.block(0, m, high+1, bSize).noalias() -= ((matH.block(0, m, high+1, bSize) * ort.segment(m, bSize)) - * (ort.segment(m, bSize)/h).transpose()); - - ort.coeffRef(m) = scale*ort.coeff(m); - matH.coeffRef(m,m-1) = scale*g; - } - } - - // Accumulate transformations (Algol's ortran). - m_eivec.setIdentity(); - - for (int m = high-1; m >= low+1; m--) - { - if (matH.coeff(m,m-1) != 0.0) - { - ort.segment(m+1, high-m) = matH.col(m-1).segment(m+1, high-m); - - int bSize = high-m+1; - m_eivec.block(m, m, bSize, bSize).noalias() += ( (ort.segment(m, bSize) / (matH.coeff(m,m-1) * ort.coeff(m))) - * (ort.segment(m, bSize).transpose() * m_eivec.block(m, m, bSize, bSize)) ); - } - } -} - -// Complex scalar division. -template -std::complex cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi) -{ - Scalar r,d; - if (ei_abs(yr) > ei_abs(yi)) - { - r = yi/yr; - d = yr + r*yi; - return std::complex((xr + r*xi)/d, (xi - r*xr)/d); - } - else - { - r = yr/yi; - d = yi + r*yr; - return std::complex((r*xr + xi)/d, (r*xi - xr)/d); - } -} - - -// Nonsymmetric reduction from Hessenberg to real Schur form. -template -void EigenSolver::hqr2(MatrixType& matH) -{ - // This is derived from the Algol procedure hqr2, - // by Martin and Wilkinson, Handbook for Auto. Comp., - // Vol.ii-Linear Algebra, and the corresponding - // Fortran subroutine in EISPACK. - - // Initialize - int nn = m_eivec.cols(); - int n = nn-1; - int low = 0; - int high = nn-1; - Scalar eps = ei_pow(Scalar(2),ei_is_same_type::ret ? Scalar(-23) : Scalar(-52)); - Scalar exshift = 0.0; - Scalar p=0,q=0,r=0,s=0,z=0,t,w,x,y; - - // Store roots isolated by balanc and compute matrix norm - // FIXME to be efficient the following would requires a triangular reduxion code - // Scalar norm = matH.upper().cwise().abs().sum() + matH.corner(BottomLeft,n,n).diagonal().cwise().abs().sum(); - Scalar norm = 0.0; - for (int j = 0; j < nn; ++j) - { - // FIXME what's the purpose of the following since the condition is always false - if ((j < low) || (j > high)) - { - m_eivalues.coeffRef(j) = Complex(matH.coeff(j,j), 0.0); - } - norm += matH.row(j).segment(std::max(j-1,0), nn-std::max(j-1,0)).cwise().abs().sum(); - } - - // Outer loop over eigenvalue index - int iter = 0; - while (n >= low) - { - // Look for single small sub-diagonal element - int l = n; - while (l > low) - { - s = ei_abs(matH.coeff(l-1,l-1)) + ei_abs(matH.coeff(l,l)); - if (s == 0.0) - s = norm; - if (ei_abs(matH.coeff(l,l-1)) < eps * s) - break; - l--; - } - - // Check for convergence - // One root found - if (l == n) - { - matH.coeffRef(n,n) = matH.coeff(n,n) + exshift; - m_eivalues.coeffRef(n) = Complex(matH.coeff(n,n), 0.0); - n--; - iter = 0; - } - else if (l == n-1) // Two roots found - { - w = matH.coeff(n,n-1) * matH.coeff(n-1,n); - p = (matH.coeff(n-1,n-1) - matH.coeff(n,n)) * Scalar(0.5); - q = p * p + w; - z = ei_sqrt(ei_abs(q)); - matH.coeffRef(n,n) = matH.coeff(n,n) + exshift; - matH.coeffRef(n-1,n-1) = matH.coeff(n-1,n-1) + exshift; - x = matH.coeff(n,n); - - // Scalar pair - if (q >= 0) - { - if (p >= 0) - z = p + z; - else - z = p - z; - - m_eivalues.coeffRef(n-1) = Complex(x + z, 0.0); - m_eivalues.coeffRef(n) = Complex(z!=0.0 ? x - w / z : m_eivalues.coeff(n-1).real(), 0.0); - - x = matH.coeff(n,n-1); - s = ei_abs(x) + ei_abs(z); - p = x / s; - q = z / s; - r = ei_sqrt(p * p+q * q); - p = p / r; - q = q / r; - - // Row modification - for (int j = n-1; j < nn; ++j) - { - z = matH.coeff(n-1,j); - matH.coeffRef(n-1,j) = q * z + p * matH.coeff(n,j); - matH.coeffRef(n,j) = q * matH.coeff(n,j) - p * z; - } - - // Column modification - for (int i = 0; i <= n; ++i) - { - z = matH.coeff(i,n-1); - matH.coeffRef(i,n-1) = q * z + p * matH.coeff(i,n); - matH.coeffRef(i,n) = q * matH.coeff(i,n) - p * z; - } - - // Accumulate transformations - for (int i = low; i <= high; ++i) - { - z = m_eivec.coeff(i,n-1); - m_eivec.coeffRef(i,n-1) = q * z + p * m_eivec.coeff(i,n); - m_eivec.coeffRef(i,n) = q * m_eivec.coeff(i,n) - p * z; - } - } - else // Complex pair - { - m_eivalues.coeffRef(n-1) = Complex(x + p, z); - m_eivalues.coeffRef(n) = Complex(x + p, -z); - } - n = n - 2; - iter = 0; - } - else // No convergence yet - { - // Form shift - x = matH.coeff(n,n); - y = 0.0; - w = 0.0; - if (l < n) - { - y = matH.coeff(n-1,n-1); - w = matH.coeff(n,n-1) * matH.coeff(n-1,n); - } - - // Wilkinson's original ad hoc shift - if (iter == 10) - { - exshift += x; - for (int i = low; i <= n; ++i) - matH.coeffRef(i,i) -= x; - s = ei_abs(matH.coeff(n,n-1)) + ei_abs(matH.coeff(n-1,n-2)); - x = y = Scalar(0.75) * s; - w = Scalar(-0.4375) * s * s; - } - - // MATLAB's new ad hoc shift - if (iter == 30) - { - s = Scalar((y - x) / 2.0); - s = s * s + w; - if (s > 0) - { - s = ei_sqrt(s); - if (y < x) - s = -s; - s = Scalar(x - w / ((y - x) / 2.0 + s)); - for (int i = low; i <= n; ++i) - matH.coeffRef(i,i) -= s; - exshift += s; - x = y = w = Scalar(0.964); - } - } - - iter = iter + 1; // (Could check iteration count here.) - - // Look for two consecutive small sub-diagonal elements - int m = n-2; - while (m >= l) - { - z = matH.coeff(m,m); - r = x - z; - s = y - z; - p = (r * s - w) / matH.coeff(m+1,m) + matH.coeff(m,m+1); - q = matH.coeff(m+1,m+1) - z - r - s; - r = matH.coeff(m+2,m+1); - s = ei_abs(p) + ei_abs(q) + ei_abs(r); - p = p / s; - q = q / s; - r = r / s; - if (m == l) { - break; - } - if (ei_abs(matH.coeff(m,m-1)) * (ei_abs(q) + ei_abs(r)) < - eps * (ei_abs(p) * (ei_abs(matH.coeff(m-1,m-1)) + ei_abs(z) + - ei_abs(matH.coeff(m+1,m+1))))) - { - break; - } - m--; - } - - for (int i = m+2; i <= n; ++i) - { - matH.coeffRef(i,i-2) = 0.0; - if (i > m+2) - matH.coeffRef(i,i-3) = 0.0; - } - - // Double QR step involving rows l:n and columns m:n - for (int k = m; k <= n-1; ++k) - { - int notlast = (k != n-1); - if (k != m) { - p = matH.coeff(k,k-1); - q = matH.coeff(k+1,k-1); - r = notlast ? matH.coeff(k+2,k-1) : Scalar(0); - x = ei_abs(p) + ei_abs(q) + ei_abs(r); - if (x != 0.0) - { - p = p / x; - q = q / x; - r = r / x; - } - } - - if (x == 0.0) - break; - - s = ei_sqrt(p * p + q * q + r * r); - - if (p < 0) - s = -s; - - if (s != 0) - { - if (k != m) - matH.coeffRef(k,k-1) = -s * x; - else if (l != m) - matH.coeffRef(k,k-1) = -matH.coeff(k,k-1); - - p = p + s; - x = p / s; - y = q / s; - z = r / s; - q = q / p; - r = r / p; - - // Row modification - for (int j = k; j < nn; ++j) - { - p = matH.coeff(k,j) + q * matH.coeff(k+1,j); - if (notlast) - { - p = p + r * matH.coeff(k+2,j); - matH.coeffRef(k+2,j) = matH.coeff(k+2,j) - p * z; - } - matH.coeffRef(k,j) = matH.coeff(k,j) - p * x; - matH.coeffRef(k+1,j) = matH.coeff(k+1,j) - p * y; - } - - // Column modification - for (int i = 0; i <= std::min(n,k+3); ++i) - { - p = x * matH.coeff(i,k) + y * matH.coeff(i,k+1); - if (notlast) - { - p = p + z * matH.coeff(i,k+2); - matH.coeffRef(i,k+2) = matH.coeff(i,k+2) - p * r; - } - matH.coeffRef(i,k) = matH.coeff(i,k) - p; - matH.coeffRef(i,k+1) = matH.coeff(i,k+1) - p * q; - } - - // Accumulate transformations - for (int i = low; i <= high; ++i) - { - p = x * m_eivec.coeff(i,k) + y * m_eivec.coeff(i,k+1); - if (notlast) - { - p = p + z * m_eivec.coeff(i,k+2); - m_eivec.coeffRef(i,k+2) = m_eivec.coeff(i,k+2) - p * r; - } - m_eivec.coeffRef(i,k) = m_eivec.coeff(i,k) - p; - m_eivec.coeffRef(i,k+1) = m_eivec.coeff(i,k+1) - p * q; - } - } // (s != 0) - } // k loop - } // check convergence - } // while (n >= low) - - // Backsubstitute to find vectors of upper triangular form - if (norm == 0.0) - { - return; - } - - for (n = nn-1; n >= 0; n--) - { - p = m_eivalues.coeff(n).real(); - q = m_eivalues.coeff(n).imag(); - - // Scalar vector - if (q == 0) - { - int l = n; - matH.coeffRef(n,n) = 1.0; - for (int i = n-1; i >= 0; i--) - { - w = matH.coeff(i,i) - p; - r = matH.row(i).segment(l,n-l+1).dot(matH.col(n).segment(l, n-l+1)); - - if (m_eivalues.coeff(i).imag() < 0.0) - { - z = w; - s = r; - } - else - { - l = i; - if (m_eivalues.coeff(i).imag() == 0.0) - { - if (w != 0.0) - matH.coeffRef(i,n) = -r / w; - else - matH.coeffRef(i,n) = -r / (eps * norm); - } - else // Solve real equations - { - x = matH.coeff(i,i+1); - y = matH.coeff(i+1,i); - q = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag(); - t = (x * s - z * r) / q; - matH.coeffRef(i,n) = t; - if (ei_abs(x) > ei_abs(z)) - matH.coeffRef(i+1,n) = (-r - w * t) / x; - else - matH.coeffRef(i+1,n) = (-s - y * t) / z; - } - - // Overflow control - t = ei_abs(matH.coeff(i,n)); - if ((eps * t) * t > 1) - matH.col(n).end(nn-i) /= t; - } - } - } - else if (q < 0) // Complex vector - { - std::complex cc; - int l = n-1; - - // Last vector component imaginary so matrix is triangular - if (ei_abs(matH.coeff(n,n-1)) > ei_abs(matH.coeff(n-1,n))) - { - matH.coeffRef(n-1,n-1) = q / matH.coeff(n,n-1); - matH.coeffRef(n-1,n) = -(matH.coeff(n,n) - p) / matH.coeff(n,n-1); - } - else - { - cc = cdiv(0.0,-matH.coeff(n-1,n),matH.coeff(n-1,n-1)-p,q); - matH.coeffRef(n-1,n-1) = ei_real(cc); - matH.coeffRef(n-1,n) = ei_imag(cc); - } - matH.coeffRef(n,n-1) = 0.0; - matH.coeffRef(n,n) = 1.0; - for (int i = n-2; i >= 0; i--) - { - Scalar ra,sa,vr,vi; - ra = matH.row(i).segment(l, n-l+1).dot(matH.col(n-1).segment(l, n-l+1)); - sa = matH.row(i).segment(l, n-l+1).dot(matH.col(n).segment(l, n-l+1)); - w = matH.coeff(i,i) - p; - - if (m_eivalues.coeff(i).imag() < 0.0) - { - z = w; - r = ra; - s = sa; - } - else - { - l = i; - if (m_eivalues.coeff(i).imag() == 0) - { - cc = cdiv(-ra,-sa,w,q); - matH.coeffRef(i,n-1) = ei_real(cc); - matH.coeffRef(i,n) = ei_imag(cc); - } - else - { - // Solve complex equations - x = matH.coeff(i,i+1); - y = matH.coeff(i+1,i); - vr = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag() - q * q; - vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q; - if ((vr == 0.0) && (vi == 0.0)) - vr = eps * norm * (ei_abs(w) + ei_abs(q) + ei_abs(x) + ei_abs(y) + ei_abs(z)); - - cc= cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi); - matH.coeffRef(i,n-1) = ei_real(cc); - matH.coeffRef(i,n) = ei_imag(cc); - if (ei_abs(x) > (ei_abs(z) + ei_abs(q))) - { - matH.coeffRef(i+1,n-1) = (-ra - w * matH.coeff(i,n-1) + q * matH.coeff(i,n)) / x; - matH.coeffRef(i+1,n) = (-sa - w * matH.coeff(i,n) - q * matH.coeff(i,n-1)) / x; - } - else - { - cc = cdiv(-r-y*matH.coeff(i,n-1),-s-y*matH.coeff(i,n),z,q); - matH.coeffRef(i+1,n-1) = ei_real(cc); - matH.coeffRef(i+1,n) = ei_imag(cc); - } - } - - // Overflow control - t = std::max(ei_abs(matH.coeff(i,n-1)),ei_abs(matH.coeff(i,n))); - if ((eps * t) * t > 1) - matH.block(i, n-1, nn-i, 2) /= t; - - } - } - } - } - - // Vectors of isolated roots - for (int i = 0; i < nn; ++i) - { - // FIXME again what's the purpose of this test ? - // in this algo low==0 and high==nn-1 !! - if (i < low || i > high) - { - m_eivec.row(i).end(nn-i) = matH.row(i).end(nn-i); - } - } - - // Back transformation to get eigenvectors of original matrix - int bRows = high-low+1; - for (int j = nn-1; j >= low; j--) - { - int bSize = std::min(j,high)-low+1; - m_eivec.col(j).segment(low, bRows) = (m_eivec.block(low, low, bRows, bSize) * matH.col(j).segment(low, bSize)); - } -} - -#endif // EIGEN_EIGENSOLVER_H diff --git a/Eigen/src/EigenSolver/HessenbergDecomposition.h b/Eigen/src/EigenSolver/HessenbergDecomposition.h deleted file mode 100644 index f782cef30..000000000 --- a/Eigen/src/EigenSolver/HessenbergDecomposition.h +++ /dev/null @@ -1,206 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008 Gael Guennebaud -// -// Eigen is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 3 of the License, or (at your option) any later version. -// -// Alternatively, you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of -// the License, or (at your option) any later version. -// -// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License and a copy of the GNU General Public License along with -// Eigen. If not, see . - -#ifndef EIGEN_HESSENBERGDECOMPOSITION_H -#define EIGEN_HESSENBERGDECOMPOSITION_H - -/** \ingroup EigenSolver_Module - * \nonstableyet - * - * \class HessenbergDecomposition - * - * \brief Reduces a squared matrix to an Hessemberg form - * - * \param MatrixType the type of the matrix of which we are computing the Hessenberg decomposition - * - * This class performs an Hessenberg decomposition of a matrix \f$ A \f$ such that: - * \f$ A = Q H Q^* \f$ where \f$ Q \f$ is unitary and \f$ H \f$ a Hessenberg matrix. - * - * \sa class Tridiagonalization, class Qr - */ -template class HessenbergDecomposition -{ - public: - - typedef _MatrixType MatrixType; - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; - - enum { - Size = MatrixType::RowsAtCompileTime, - SizeMinusOne = MatrixType::RowsAtCompileTime==Dynamic - ? Dynamic - : MatrixType::RowsAtCompileTime-1 - }; - - typedef Matrix CoeffVectorType; - typedef Matrix DiagonalType; - typedef Matrix SubDiagonalType; - - typedef typename NestByValue >::RealReturnType DiagonalReturnType; - - typedef typename NestByValue >,0 > >::RealReturnType SubDiagonalReturnType; - - /** This constructor initializes a HessenbergDecomposition object for - * further use with HessenbergDecomposition::compute() - */ - HessenbergDecomposition(int size = Size==Dynamic ? 2 : Size) - : m_matrix(size,size), m_hCoeffs(size-1) - {} - - HessenbergDecomposition(const MatrixType& matrix) - : m_matrix(matrix), - m_hCoeffs(matrix.cols()-1) - { - _compute(m_matrix, m_hCoeffs); - } - - /** Computes or re-compute the Hessenberg decomposition for the matrix \a matrix. - * - * This method allows to re-use the allocated data. - */ - void compute(const MatrixType& matrix) - { - m_matrix = matrix; - m_hCoeffs.resize(matrix.rows()-1,1); - _compute(m_matrix, m_hCoeffs); - } - - /** \returns the householder coefficients allowing to - * reconstruct the matrix Q from the packed data. - * - * \sa packedMatrix() - */ - CoeffVectorType householderCoefficients() const { return m_hCoeffs; } - - /** \returns the internal result of the decomposition. - * - * The returned matrix contains the following information: - * - the upper part and lower sub-diagonal represent the Hessenberg matrix H - * - the rest of the lower part contains the Householder vectors that, combined with - * Householder coefficients returned by householderCoefficients(), - * allows to reconstruct the matrix Q as follow: - * Q = H_{N-1} ... H_1 H_0 - * where the matrices H are the Householder transformation: - * H_i = (I - h_i * v_i * v_i') - * where h_i == householderCoefficients()[i] and v_i is a Householder vector: - * v_i = [ 0, ..., 0, 1, M(i+2,i), ..., M(N-1,i) ] - * - * See LAPACK for further details on this packed storage. - */ - const MatrixType& packedMatrix(void) const { return m_matrix; } - - MatrixType matrixQ() const; - MatrixType matrixH() const; - - private: - - static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs); - - protected: - MatrixType m_matrix; - CoeffVectorType m_hCoeffs; -}; - -#ifndef EIGEN_HIDE_HEAVY_CODE - -/** \internal - * Performs a tridiagonal decomposition of \a matA in place. - * - * \param matA the input selfadjoint matrix - * \param hCoeffs returned Householder coefficients - * - * The result is written in the lower triangular part of \a matA. - * - * Implemented from Golub's "Matrix Computations", algorithm 8.3.1. - * - * \sa packedMatrix() - */ -template -void HessenbergDecomposition::_compute(MatrixType& matA, CoeffVectorType& hCoeffs) -{ - assert(matA.rows()==matA.cols()); - int n = matA.rows(); - Matrix temp(n); - for (int i = 0; i -typename HessenbergDecomposition::MatrixType -HessenbergDecomposition::matrixQ() const -{ - int n = m_matrix.rows(); - MatrixType matQ = MatrixType::Identity(n,n); - Matrix temp(n); - for (int i = n-2; i>=0; i--) - { - matQ.corner(BottomRight,n-i-1,n-i-1) - .applyHouseholderOnTheLeft(m_matrix.col(i).end(n-i-2), ei_conj(m_hCoeffs.coeff(i)), &temp.coeffRef(0,0)); - } - return matQ; -} - -#endif // EIGEN_HIDE_HEAVY_CODE - -/** constructs and returns the matrix H. - * Note that the matrix H is equivalent to the upper part of the packed matrix - * (including the lower sub-diagonal). Therefore, it might be often sufficient - * to directly use the packed matrix instead of creating a new one. - */ -template -typename HessenbergDecomposition::MatrixType -HessenbergDecomposition::matrixH() const -{ - // FIXME should this function (and other similar) rather take a matrix as argument - // and fill it (to avoid temporaries) - int n = m_matrix.rows(); - MatrixType matH = m_matrix; - if (n>2) - matH.corner(BottomLeft,n-2, n-2).template triangularView().setZero(); - return matH; -} - -#endif // EIGEN_HESSENBERGDECOMPOSITION_H diff --git a/Eigen/src/EigenSolver/SelfAdjointEigenSolver.h b/Eigen/src/EigenSolver/SelfAdjointEigenSolver.h deleted file mode 100644 index 40b06df2c..000000000 --- a/Eigen/src/EigenSolver/SelfAdjointEigenSolver.h +++ /dev/null @@ -1,366 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008 Gael Guennebaud -// -// Eigen is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 3 of the License, or (at your option) any later version. -// -// Alternatively, you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of -// the License, or (at your option) any later version. -// -// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License and a copy of the GNU General Public License along with -// Eigen. If not, see . - -#ifndef EIGEN_SELFADJOINTEIGENSOLVER_H -#define EIGEN_SELFADJOINTEIGENSOLVER_H - -/** \eigensolver_module \ingroup EigenSolver_Module - * \nonstableyet - * - * \class SelfAdjointEigenSolver - * - * \brief Eigen values/vectors solver for selfadjoint matrix - * - * \param MatrixType the type of the matrix of which we are computing the eigen decomposition - * - * \note MatrixType must be an actual Matrix type, it can't be an expression type. - * - * \sa MatrixBase::eigenvalues(), class EigenSolver - */ -template class SelfAdjointEigenSolver -{ - public: - - enum {Size = _MatrixType::RowsAtCompileTime }; - typedef _MatrixType MatrixType; - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; - typedef std::complex Complex; - typedef Matrix RealVectorType; - typedef Matrix RealVectorTypeX; - typedef Tridiagonalization TridiagonalizationType; -// typedef typename TridiagonalizationType::TridiagonalMatrixType TridiagonalMatrixType; - - SelfAdjointEigenSolver() - : m_eivec(int(Size), int(Size)), - m_eivalues(int(Size)) - { - ei_assert(Size!=Dynamic); - } - - SelfAdjointEigenSolver(int size) - : m_eivec(size, size), - m_eivalues(size) - {} - - /** Constructors computing the eigenvalues of the selfadjoint matrix \a matrix, - * as well as the eigenvectors if \a computeEigenvectors is true. - * - * \sa compute(MatrixType,bool), SelfAdjointEigenSolver(MatrixType,MatrixType,bool) - */ - SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true) - : m_eivec(matrix.rows(), matrix.cols()), - m_eivalues(matrix.cols()) - { - compute(matrix, computeEigenvectors); - } - - /** Constructors computing the eigenvalues of the generalized eigen problem - * \f$ Ax = lambda B x \f$ with \a matA the selfadjoint matrix \f$ A \f$ - * and \a matB the positive definite matrix \f$ B \f$ . The eigenvectors - * are computed if \a computeEigenvectors is true. - * - * \sa compute(MatrixType,MatrixType,bool), SelfAdjointEigenSolver(MatrixType,bool) - */ - SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true) - : m_eivec(matA.rows(), matA.cols()), - m_eivalues(matA.cols()) - { - compute(matA, matB, computeEigenvectors); - } - - SelfAdjointEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true); - - SelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true); - - /** \returns the computed eigen vectors as a matrix of column vectors */ - MatrixType eigenvectors(void) const - { - #ifndef NDEBUG - ei_assert(m_eigenvectorsOk); - #endif - return m_eivec; - } - - /** \returns the computed eigen values */ - RealVectorType eigenvalues(void) const { return m_eivalues; } - - /** \returns the positive square root of the matrix - * - * \note the matrix itself must be positive in order for this to make sense. - */ - MatrixType operatorSqrt() const - { - return m_eivec * m_eivalues.cwise().sqrt().asDiagonal() * m_eivec.adjoint(); - } - - /** \returns the positive inverse square root of the matrix - * - * \note the matrix itself must be positive definite in order for this to make sense. - */ - MatrixType operatorInverseSqrt() const - { - return m_eivec * m_eivalues.cwise().inverse().cwise().sqrt().asDiagonal() * m_eivec.adjoint(); - } - - - protected: - MatrixType m_eivec; - RealVectorType m_eivalues; - #ifndef NDEBUG - bool m_eigenvectorsOk; - #endif -}; - -#ifndef EIGEN_HIDE_HEAVY_CODE - -/** \internal - * - * \eigensolver_module - * - * Performs a QR step on a tridiagonal symmetric matrix represented as a - * pair of two vectors \a diag and \a subdiag. - * - * \param matA the input selfadjoint matrix - * \param hCoeffs returned Householder coefficients - * - * For compilation efficiency reasons, this procedure does not use eigen expression - * for its arguments. - * - * Implemented from Golub's "Matrix Computations", algorithm 8.3.2: - * "implicit symmetric QR step with Wilkinson shift" - */ -template -static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int start, int end, Scalar* matrixQ, int n); - -/** Computes the eigenvalues of the selfadjoint matrix \a matrix, - * as well as the eigenvectors if \a computeEigenvectors is true. - * - * \sa SelfAdjointEigenSolver(MatrixType,bool), compute(MatrixType,MatrixType,bool) - */ -template -SelfAdjointEigenSolver& SelfAdjointEigenSolver::compute(const MatrixType& matrix, bool computeEigenvectors) -{ - #ifndef NDEBUG - m_eigenvectorsOk = computeEigenvectors; - #endif - assert(matrix.cols() == matrix.rows()); - int n = matrix.cols(); - m_eivalues.resize(n,1); - - if(n==1) - { - m_eivalues.coeffRef(0,0) = ei_real(matrix.coeff(0,0)); - m_eivec.setOnes(); - return *this; - } - - m_eivec = matrix; - - // FIXME, should tridiag be a local variable of this function or an attribute of SelfAdjointEigenSolver ? - // the latter avoids multiple memory allocation when the same SelfAdjointEigenSolver is used multiple times... - // (same for diag and subdiag) - RealVectorType& diag = m_eivalues; - typename TridiagonalizationType::SubDiagonalType subdiag(n-1); - TridiagonalizationType::decomposeInPlace(m_eivec, diag, subdiag, computeEigenvectors); - - int end = n-1; - int start = 0; - while (end>0) - { - for (int i = start; i0 && subdiag[end-1]==0) - end--; - if (end<=0) - break; - start = end - 1; - while (start>0 && subdiag[start-1]!=0) - start--; - - ei_tridiagonal_qr_step(diag.data(), subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n); - } - - // Sort eigenvalues and corresponding vectors. - // TODO make the sort optional ? - // TODO use a better sort algorithm !! - for (int i = 0; i < n-1; ++i) - { - int k; - m_eivalues.segment(i,n-i).minCoeff(&k); - if (k > 0) - { - std::swap(m_eivalues[i], m_eivalues[k+i]); - m_eivec.col(i).swap(m_eivec.col(k+i)); - } - } - return *this; -} - -/** Computes the eigenvalues of the generalized eigen problem - * \f$ Ax = lambda B x \f$ with \a matA the selfadjoint matrix \f$ A \f$ - * and \a matB the positive definite matrix \f$ B \f$ . The eigenvectors - * are computed if \a computeEigenvectors is true. - * - * \sa SelfAdjointEigenSolver(MatrixType,MatrixType,bool), compute(MatrixType,bool) - */ -template -SelfAdjointEigenSolver& SelfAdjointEigenSolver:: -compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors) -{ - ei_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows()); - - // Compute the cholesky decomposition of matB = L L' - LLT cholB(matB); - - // compute C = inv(L) A inv(L') - MatrixType matC = matA; - cholB.matrixL().solveInPlace(matC); - // FIXME since we currently do not support A * inv(L'), let's do (inv(L) A')' : - matC.adjointInPlace(); - cholB.matrixL().solveInPlace(matC); - matC.adjointInPlace(); - // this version works too: -// matC = matC.transpose(); -// cholB.matrixL().conjugate().template marked().solveTriangularInPlace(matC); -// matC = matC.transpose(); - // FIXME: this should work: (currently it only does for small matrices) -// Transpose trMatC(matC); -// cholB.matrixL().conjugate().eval().template marked().solveTriangularInPlace(trMatC); - - compute(matC, computeEigenvectors); - - if (computeEigenvectors) - { - // transform back the eigen vectors: evecs = inv(U) * evecs - cholB.matrixU().solveInPlace(m_eivec); - for (int i=0; i -inline Matrix::Scalar>::Real, ei_traits::ColsAtCompileTime, 1> -MatrixBase::eigenvalues() const -{ - ei_assert(Flags&SelfAdjointBit); - return SelfAdjointEigenSolver(eval(),false).eigenvalues(); -} - -template -struct ei_operatorNorm_selector -{ - static inline typename NumTraits::Scalar>::Real - operatorNorm(const MatrixBase& m) - { - // FIXME if it is really guaranteed that the eigenvalues are already sorted, - // then we don't need to compute a maxCoeff() here, comparing the 1st and last ones is enough. - return m.eigenvalues().cwise().abs().maxCoeff(); - } -}; - -template struct ei_operatorNorm_selector -{ - static inline typename NumTraits::Scalar>::Real - operatorNorm(const MatrixBase& m) - { - typename Derived::PlainMatrixType m_eval(m); - // FIXME if it is really guaranteed that the eigenvalues are already sorted, - // then we don't need to compute a maxCoeff() here, comparing the 1st and last ones is enough. - return ei_sqrt( - (m_eval*m_eval.adjoint()) - .template marked() - .eigenvalues() - .maxCoeff() - ); - } -}; - -/** \eigensolver_module - * - * \returns the matrix norm of this matrix. - */ -template -inline typename NumTraits::Scalar>::Real -MatrixBase::operatorNorm() const -{ - return ei_operatorNorm_selector - ::operatorNorm(derived()); -} - -#ifndef EIGEN_EXTERN_INSTANTIATIONS -template -static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int start, int end, Scalar* matrixQ, int n) -{ - RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5); - RealScalar e2 = ei_abs2(subdiag[end-1]); - RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * ei_sqrt(td*td + e2)); - RealScalar x = diag[start] - mu; - RealScalar z = subdiag[start]; - - for (int k = start; k < end; ++k) - { - PlanarRotation rot; - rot.makeGivens(x, z); - - // do T = G' T G - RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k]; - RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k+1]; - - diag[k] = rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k+1]); - diag[k+1] = rot.s() * sdk + rot.c() * dkp1; - subdiag[k] = rot.c() * sdk - rot.s() * dkp1; - - if (k > start) - subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z; - - x = subdiag[k]; - - if (k < end - 1) - { - z = -rot.s() * subdiag[k+1]; - subdiag[k + 1] = rot.c() * subdiag[k+1]; - } - - // apply the givens rotation to the unit matrix Q = Q * G - if (matrixQ) - { - Map > q(matrixQ,n,n); - q.applyOnTheRight(k,k+1,rot); - } - } -} -#endif - -#endif // EIGEN_SELFADJOINTEIGENSOLVER_H diff --git a/Eigen/src/EigenSolver/Tridiagonalization.h b/Eigen/src/EigenSolver/Tridiagonalization.h deleted file mode 100644 index e0bff17b9..000000000 --- a/Eigen/src/EigenSolver/Tridiagonalization.h +++ /dev/null @@ -1,317 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008 Gael Guennebaud -// -// Eigen is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 3 of the License, or (at your option) any later version. -// -// Alternatively, you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of -// the License, or (at your option) any later version. -// -// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License and a copy of the GNU General Public License along with -// Eigen. If not, see . - -#ifndef EIGEN_TRIDIAGONALIZATION_H -#define EIGEN_TRIDIAGONALIZATION_H - -/** \ingroup EigenSolver_Module - * \nonstableyet - * - * \class Tridiagonalization - * - * \brief Trigiagonal decomposition of a selfadjoint matrix - * - * \param MatrixType the type of the matrix of which we are performing the tridiagonalization - * - * This class performs a tridiagonal decomposition of a selfadjoint matrix \f$ A \f$ such that: - * \f$ A = Q T Q^* \f$ where \f$ Q \f$ is unitary and \f$ T \f$ a real symmetric tridiagonal matrix. - * - * \sa MatrixBase::tridiagonalize() - */ -template class Tridiagonalization -{ - public: - - typedef _MatrixType MatrixType; - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; - typedef typename ei_packet_traits::type Packet; - - enum { - Size = MatrixType::RowsAtCompileTime, - SizeMinusOne = MatrixType::RowsAtCompileTime==Dynamic - ? Dynamic - : MatrixType::RowsAtCompileTime-1, - PacketSize = ei_packet_traits::size - }; - - typedef Matrix CoeffVectorType; - typedef Matrix DiagonalType; - typedef Matrix SubDiagonalType; - - typedef typename ei_meta_if::IsComplex, - typename NestByValue >::RealReturnType, - Diagonal - >::ret DiagonalReturnType; - - typedef typename ei_meta_if::IsComplex, - typename NestByValue >,0 > >::RealReturnType, - Diagonal< - NestByValue >,0 > - >::ret SubDiagonalReturnType; - - /** This constructor initializes a Tridiagonalization object for - * further use with Tridiagonalization::compute() - */ - Tridiagonalization(int size = Size==Dynamic ? 2 : Size) - : m_matrix(size,size), m_hCoeffs(size-1) - {} - - Tridiagonalization(const MatrixType& matrix) - : m_matrix(matrix), m_hCoeffs(matrix.cols()-1) - { - _compute(m_matrix, m_hCoeffs); - } - - /** Computes or re-compute the tridiagonalization for the matrix \a matrix. - * - * This method allows to re-use the allocated data. - */ - void compute(const MatrixType& matrix) - { - m_matrix = matrix; - m_hCoeffs.resize(matrix.rows()-1, 1); - _compute(m_matrix, m_hCoeffs); - } - - /** \returns the householder coefficients allowing to - * reconstruct the matrix Q from the packed data. - * - * \sa packedMatrix() - */ - inline CoeffVectorType householderCoefficients(void) const { return m_hCoeffs; } - - /** \returns the internal result of the decomposition. - * - * The returned matrix contains the following information: - * - the strict upper part is equal to the input matrix A - * - the diagonal and lower sub-diagonal represent the tridiagonal symmetric matrix (real). - * - the rest of the lower part contains the Householder vectors that, combined with - * Householder coefficients returned by householderCoefficients(), - * allows to reconstruct the matrix Q as follow: - * Q = H_{N-1} ... H_1 H_0 - * where the matrices H are the Householder transformations: - * H_i = (I - h_i * v_i * v_i') - * where h_i == householderCoefficients()[i] and v_i is a Householder vector: - * v_i = [ 0, ..., 0, 1, M(i+2,i), ..., M(N-1,i) ] - * - * See LAPACK for further details on this packed storage. - */ - inline const MatrixType& packedMatrix(void) const { return m_matrix; } - - MatrixType matrixQ() const; - template void matrixQInPlace(MatrixBase* q) const; - MatrixType matrixT() const; - const DiagonalReturnType diagonal(void) const; - const SubDiagonalReturnType subDiagonal(void) const; - - static void decomposeInPlace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ = true); - - static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs); - - protected: - - static void _decomposeInPlace3x3(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ = true); - - MatrixType m_matrix; - CoeffVectorType m_hCoeffs; -}; - -/** \returns an expression of the diagonal vector */ -template -const typename Tridiagonalization::DiagonalReturnType -Tridiagonalization::diagonal(void) const -{ - return m_matrix.diagonal().nestByValue(); -} - -/** \returns an expression of the sub-diagonal vector */ -template -const typename Tridiagonalization::SubDiagonalReturnType -Tridiagonalization::subDiagonal(void) const -{ - int n = m_matrix.rows(); - return Block(m_matrix, 1, 0, n-1,n-1) - .nestByValue().diagonal().nestByValue(); -} - -/** constructs and returns the tridiagonal matrix T. - * Note that the matrix T is equivalent to the diagonal and sub-diagonal of the packed matrix. - * Therefore, it might be often sufficient to directly use the packed matrix, or the vector - * expressions returned by diagonal() and subDiagonal() instead of creating a new matrix. - */ -template -typename Tridiagonalization::MatrixType -Tridiagonalization::matrixT(void) const -{ - // FIXME should this function (and other similar ones) rather take a matrix as argument - // and fill it ? (to avoid temporaries) - int n = m_matrix.rows(); - MatrixType matT = m_matrix; - matT.corner(TopRight,n-1, n-1).diagonal() = subDiagonal().template cast().conjugate(); - if (n>2) - { - matT.corner(TopRight,n-2, n-2).template triangularView().setZero(); - matT.corner(BottomLeft,n-2, n-2).template triangularView().setZero(); - } - return matT; -} - -#ifndef EIGEN_HIDE_HEAVY_CODE - -/** \internal - * Performs a tridiagonal decomposition of \a matA in place. - * - * \param matA the input selfadjoint matrix - * \param hCoeffs returned Householder coefficients - * - * The result is written in the lower triangular part of \a matA. - * - * Implemented from Golub's "Matrix Computations", algorithm 8.3.1. - * - * \sa packedMatrix() - */ -template -void Tridiagonalization::_compute(MatrixType& matA, CoeffVectorType& hCoeffs) -{ - assert(matA.rows()==matA.cols()); - int n = matA.rows(); - Matrix aux(n); - for (int i = 0; i() - * (ei_conj(h) * matA.col(i).end(remainingSize))); - - hCoeffs.end(n-i-1) += (ei_conj(h)*Scalar(-0.5)*(hCoeffs.end(remainingSize).dot(matA.col(i).end(remainingSize)))) * matA.col(i).end(n-i-1); - - matA.corner(BottomRight, remainingSize, remainingSize).template selfadjointView() - .rankUpdate(matA.col(i).end(remainingSize), hCoeffs.end(remainingSize), -1); - - matA.col(i).coeffRef(i+1) = beta; - hCoeffs.coeffRef(i) = h; - } -} - -/** reconstructs and returns the matrix Q */ -template -typename Tridiagonalization::MatrixType -Tridiagonalization::matrixQ(void) const -{ - MatrixType matQ; - matrixQInPlace(&matQ); - return matQ; -} - -template -template -void Tridiagonalization::matrixQInPlace(MatrixBase* q) const -{ - QDerived& matQ = q->derived(); - int n = m_matrix.rows(); - matQ = MatrixType::Identity(n,n); - Matrix aux(n); - for (int i = n-2; i>=0; i--) - { - matQ.corner(BottomRight,n-i-1,n-i-1) - .applyHouseholderOnTheLeft(m_matrix.col(i).end(n-i-2), ei_conj(m_hCoeffs.coeff(i)), &aux.coeffRef(0,0)); - } -} - -/** Performs a full decomposition in place */ -template -void Tridiagonalization::decomposeInPlace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ) -{ - int n = mat.rows(); - ei_assert(mat.cols()==n && diag.size()==n && subdiag.size()==n-1); - if (n==3 && (!NumTraits::IsComplex) ) - { - _decomposeInPlace3x3(mat, diag, subdiag, extractQ); - } - else - { - Tridiagonalization tridiag(mat); - diag = tridiag.diagonal(); - subdiag = tridiag.subDiagonal(); - if (extractQ) - tridiag.matrixQInPlace(&mat); - } -} - -/** \internal - * Optimized path for 3x3 matrices. - * Especially useful for plane fitting. - */ -template -void Tridiagonalization::_decomposeInPlace3x3(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ) -{ - diag[0] = ei_real(mat(0,0)); - RealScalar v1norm2 = ei_abs2(mat(0,2)); - if (ei_isMuchSmallerThan(v1norm2, RealScalar(1))) - { - diag[1] = ei_real(mat(1,1)); - diag[2] = ei_real(mat(2,2)); - subdiag[0] = ei_real(mat(0,1)); - subdiag[1] = ei_real(mat(1,2)); - if (extractQ) - mat.setIdentity(); - } - else - { - RealScalar beta = ei_sqrt(ei_abs2(mat(0,1))+v1norm2); - RealScalar invBeta = RealScalar(1)/beta; - Scalar m01 = mat(0,1) * invBeta; - Scalar m02 = mat(0,2) * invBeta; - Scalar q = RealScalar(2)*m01*mat(1,2) + m02*(mat(2,2) - mat(1,1)); - diag[1] = ei_real(mat(1,1) + m02*q); - diag[2] = ei_real(mat(2,2) - m02*q); - subdiag[0] = beta; - subdiag[1] = ei_real(mat(1,2) - m01 * q); - if (extractQ) - { - mat(0,0) = 1; - mat(0,1) = 0; - mat(0,2) = 0; - mat(1,0) = 0; - mat(1,1) = m01; - mat(1,2) = m02; - mat(2,0) = 0; - mat(2,1) = m02; - mat(2,2) = -m01; - } - } -} - -#endif // EIGEN_HIDE_HEAVY_CODE - -#endif // EIGEN_TRIDIAGONALIZATION_H diff --git a/Eigen/src/Eigenvalues/CMakeLists.txt b/Eigen/src/Eigenvalues/CMakeLists.txt new file mode 100644 index 000000000..193e02685 --- /dev/null +++ b/Eigen/src/Eigenvalues/CMakeLists.txt @@ -0,0 +1,6 @@ +FILE(GLOB Eigen_EIGENVALUES_SRCS "*.h") + +INSTALL(FILES + ${Eigen_EIGENVALUES_SRCS} + DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Eigenvalues COMPONENT Devel + ) diff --git a/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/Eigen/src/Eigenvalues/ComplexEigenSolver.h new file mode 100644 index 000000000..666381949 --- /dev/null +++ b/Eigen/src/Eigenvalues/ComplexEigenSolver.h @@ -0,0 +1,148 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Claire Maurice +// Copyright (C) 2009 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H +#define EIGEN_COMPLEX_EIGEN_SOLVER_H + +/** \eigenvalues_module \ingroup Eigenvalues_Module + * \nonstableyet + * + * \class ComplexEigenSolver + * + * \brief Eigen values/vectors solver for general complex matrices + * + * \param MatrixType the type of the matrix of which we are computing the eigen decomposition + * + * \sa class EigenSolver, class SelfAdjointEigenSolver + */ +template class ComplexEigenSolver +{ + public: + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef std::complex Complex; + typedef Matrix EigenvalueType; + typedef Matrix EigenvectorType; + + /** + * \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via ComplexEigenSolver::compute(const MatrixType&). + */ + ComplexEigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false) + {} + + ComplexEigenSolver(const MatrixType& matrix) + : m_eivec(matrix.rows(),matrix.cols()), + m_eivalues(matrix.cols()), + m_isInitialized(false) + { + compute(matrix); + } + + EigenvectorType eigenvectors(void) const + { + ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); + return m_eivec; + } + + EigenvalueType eigenvalues() const + { + ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); + return m_eivalues; + } + + void compute(const MatrixType& matrix); + + protected: + MatrixType m_eivec; + EigenvalueType m_eivalues; + bool m_isInitialized; +}; + + +template +void ComplexEigenSolver::compute(const MatrixType& matrix) +{ + // this code is inspired from Jampack + assert(matrix.cols() == matrix.rows()); + int n = matrix.cols(); + m_eivalues.resize(n,1); + + RealScalar eps = epsilon(); + + // Reduce to complex Schur form + ComplexSchur schur(matrix); + + m_eivalues = schur.matrixT().diagonal(); + + m_eivec.setZero(); + + Scalar d2, z; + RealScalar norm = matrix.norm(); + + // compute the (normalized) eigenvectors + for(int k=n-1 ; k>=0 ; k--) + { + d2 = schur.matrixT().coeff(k,k); + m_eivec.coeffRef(k,k) = Scalar(1.0,0.0); + for(int i=k-1 ; i>=0 ; i--) + { + m_eivec.coeffRef(i,k) = -schur.matrixT().coeff(i,k); + if(k-i-1>0) + m_eivec.coeffRef(i,k) -= (schur.matrixT().row(i).segment(i+1,k-i-1) * m_eivec.col(k).segment(i+1,k-i-1)).value(); + z = schur.matrixT().coeff(i,i) - d2; + if(z==Scalar(0)) + ei_real_ref(z) = eps * norm; + m_eivec.coeffRef(i,k) = m_eivec.coeff(i,k) / z; + + } + m_eivec.col(k).normalize(); + } + + m_eivec = schur.matrixU() * m_eivec; + m_isInitialized = true; + + // sort the eigenvalues + { + for (int i=0; i +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_COMPLEX_SCHUR_H +#define EIGEN_COMPLEX_SCHUR_H + +/** \eigenvalues_module \ingroup Eigenvalues_Module + * \nonstableyet + * + * \class ComplexShur + * + * \brief Performs a complex Shur decomposition of a real or complex square matrix + * + */ +template class ComplexSchur +{ + public: + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef std::complex Complex; + typedef Matrix ComplexMatrixType; + + /** + * \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via ComplexSchur::compute(const MatrixType&). + */ + ComplexSchur() : m_matT(), m_matU(), m_isInitialized(false) + {} + + ComplexSchur(const MatrixType& matrix) + : m_matT(matrix.rows(),matrix.cols()), + m_matU(matrix.rows(),matrix.cols()), + m_isInitialized(false) + { + compute(matrix); + } + + ComplexMatrixType matrixU() const + { + ei_assert(m_isInitialized && "ComplexSchur is not initialized."); + return m_matU; + } + + ComplexMatrixType matrixT() const + { + ei_assert(m_isInitialized && "ComplexShur is not initialized."); + return m_matT; + } + + void compute(const MatrixType& matrix); + + protected: + ComplexMatrixType m_matT, m_matU; + bool m_isInitialized; +}; + +/** Computes the principal value of the square root of the complex \a z. */ +template +std::complex ei_sqrt(const std::complex &z) +{ + RealScalar t, tre, tim; + + t = ei_abs(z); + + if (ei_abs(ei_real(z)) <= ei_abs(ei_imag(z))) + { + // No cancellation in these formulas + tre = ei_sqrt(0.5*(t + ei_real(z))); + tim = ei_sqrt(0.5*(t - ei_real(z))); + } + else + { + // Stable computation of the above formulas + if (z.real() > 0) + { + tre = t + z.real(); + tim = ei_abs(ei_imag(z))*ei_sqrt(0.5/tre); + tre = ei_sqrt(0.5*tre); + } + else + { + tim = t - z.real(); + tre = ei_abs(ei_imag(z))*ei_sqrt(0.5/tim); + tim = ei_sqrt(0.5*tim); + } + } + if(z.imag() < 0) + tim = -tim; + + return (std::complex(tre,tim)); + +} + +template +void ComplexSchur::compute(const MatrixType& matrix) +{ + // this code is inspired from Jampack + assert(matrix.cols() == matrix.rows()); + int n = matrix.cols(); + + // Reduce to Hessenberg form + HessenbergDecomposition hess(matrix); + + m_matT = hess.matrixH(); + m_matU = hess.matrixQ(); + + int iu = m_matT.cols() - 1; + int il; + RealScalar d,sd,sf; + Complex c,b,disc,r1,r2,kappa; + + RealScalar eps = epsilon(); + + int iter = 0; + while(true) + { + //locate the range in which to iterate + while(iu > 0) + { + d = ei_norm1(m_matT.coeffRef(iu,iu)) + ei_norm1(m_matT.coeffRef(iu-1,iu-1)); + sd = ei_norm1(m_matT.coeffRef(iu,iu-1)); + + if(sd >= eps * d) break; // FIXME : precision criterion ?? + + m_matT.coeffRef(iu,iu-1) = Complex(0); + iter = 0; + --iu; + } + if(iu==0) break; + iter++; + + if(iter >= 30) + { + // FIXME : what to do when iter==MAXITER ?? + std::cerr << "MAXITER" << std::endl; + return; + } + + il = iu-1; + while( il > 0 ) + { + // check if the current 2x2 block on the diagonal is upper triangular + d = ei_norm1(m_matT.coeffRef(il,il)) + ei_norm1(m_matT.coeffRef(il-1,il-1)); + sd = ei_norm1(m_matT.coeffRef(il,il-1)); + + if(sd < eps * d) break; // FIXME : precision criterion ?? + + --il; + } + + if( il != 0 ) m_matT.coeffRef(il,il-1) = Complex(0); + + // compute the shift (the normalization by sf is to avoid under/overflow) + Matrix t = m_matT.template block<2,2>(iu-1,iu-1); + sf = t.cwise().abs().sum(); + t /= sf; + + c = t.determinant(); + b = t.diagonal().sum(); + + disc = ei_sqrt(b*b - RealScalar(4)*c); + + r1 = (b+disc)/RealScalar(2); + r2 = (b-disc)/RealScalar(2); + + if(ei_norm1(r1) > ei_norm1(r2)) + r2 = c/r1; + else + r1 = c/r2; + + if(ei_norm1(r1-t.coeff(1,1)) < ei_norm1(r2-t.coeff(1,1))) + kappa = sf * r1; + else + kappa = sf * r2; + + // perform the QR step using Givens rotations + PlanarRotation rot; + rot.makeGivens(m_matT.coeff(il,il) - kappa, m_matT.coeff(il+1,il)); + + for(int i=il ; i +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_EIGENSOLVER_H +#define EIGEN_EIGENSOLVER_H + +/** \eigenvalues_module \ingroup Eigenvalues_Module + * \nonstableyet + * + * \class EigenSolver + * + * \brief Eigen values/vectors solver for non selfadjoint matrices + * + * \param MatrixType the type of the matrix of which we are computing the eigen decomposition + * + * Currently it only support real matrices. + * + * \note this code was adapted from JAMA (public domain) + * + * \sa MatrixBase::eigenvalues(), SelfAdjointEigenSolver + */ +template class EigenSolver +{ + public: + + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef std::complex Complex; + typedef Matrix EigenvalueType; + typedef Matrix EigenvectorType; + typedef Matrix RealVectorType; + typedef Matrix RealVectorTypeX; + + /** + * \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via EigenSolver::compute(const MatrixType&). + */ + EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false) {} + + EigenSolver(const MatrixType& matrix) + : m_eivec(matrix.rows(), matrix.cols()), + m_eivalues(matrix.cols()), + m_isInitialized(false) + { + compute(matrix); + } + + + EigenvectorType eigenvectors(void) const; + + /** \returns a real matrix V of pseudo eigenvectors. + * + * Let D be the block diagonal matrix with the real eigenvalues in 1x1 blocks, + * and any complex values u+iv in 2x2 blocks [u v ; -v u]. Then, the matrices D + * and V satisfy A*V = V*D. + * + * More precisely, if the diagonal matrix of the eigen values is:\n + * \f$ + * \left[ \begin{array}{cccccc} + * u+iv & & & & & \\ + * & u-iv & & & & \\ + * & & a+ib & & & \\ + * & & & a-ib & & \\ + * & & & & x & \\ + * & & & & & y \\ + * \end{array} \right] + * \f$ \n + * then, we have:\n + * \f$ + * D =\left[ \begin{array}{cccccc} + * u & v & & & & \\ + * -v & u & & & & \\ + * & & a & b & & \\ + * & & -b & a & & \\ + * & & & & x & \\ + * & & & & & y \\ + * \end{array} \right] + * \f$ + * + * \sa pseudoEigenvalueMatrix() + */ + const MatrixType& pseudoEigenvectors() const + { + ei_assert(m_isInitialized && "EigenSolver is not initialized."); + return m_eivec; + } + + MatrixType pseudoEigenvalueMatrix() const; + + /** \returns the eigenvalues as a column vector */ + EigenvalueType eigenvalues() const + { + ei_assert(m_isInitialized && "EigenSolver is not initialized."); + return m_eivalues; + } + + EigenSolver& compute(const MatrixType& matrix); + + private: + + void orthes(MatrixType& matH, RealVectorType& ort); + void hqr2(MatrixType& matH); + + protected: + MatrixType m_eivec; + EigenvalueType m_eivalues; + bool m_isInitialized; +}; + +/** \returns the real block diagonal matrix D of the eigenvalues. + * + * See pseudoEigenvectors() for the details. + */ +template +MatrixType EigenSolver::pseudoEigenvalueMatrix() const +{ + ei_assert(m_isInitialized && "EigenSolver is not initialized."); + int n = m_eivec.cols(); + MatrixType matD = MatrixType::Zero(n,n); + for (int i=0; i(i,i) << ei_real(m_eivalues.coeff(i)), ei_imag(m_eivalues.coeff(i)), + -ei_imag(m_eivalues.coeff(i)), ei_real(m_eivalues.coeff(i)); + ++i; + } + } + return matD; +} + +/** \returns the normalized complex eigenvectors as a matrix of column vectors. + * + * \sa eigenvalues(), pseudoEigenvectors() + */ +template +typename EigenSolver::EigenvectorType EigenSolver::eigenvectors(void) const +{ + ei_assert(m_isInitialized && "EigenSolver is not initialized."); + int n = m_eivec.cols(); + EigenvectorType matV(n,n); + for (int j=0; j(); + } + else + { + // we have a pair of complex eigen values + for (int i=0; i +EigenSolver& EigenSolver::compute(const MatrixType& matrix) +{ + assert(matrix.cols() == matrix.rows()); + int n = matrix.cols(); + m_eivalues.resize(n,1); + + MatrixType matH = matrix; + RealVectorType ort(n); + + // Reduce to Hessenberg form. + orthes(matH, ort); + + // Reduce Hessenberg to real Schur form. + hqr2(matH); + + m_isInitialized = true; + return *this; +} + +// Nonsymmetric reduction to Hessenberg form. +template +void EigenSolver::orthes(MatrixType& matH, RealVectorType& ort) +{ + // This is derived from the Algol procedures orthes and ortran, + // by Martin and Wilkinson, Handbook for Auto. Comp., + // Vol.ii-Linear Algebra, and the corresponding + // Fortran subroutines in EISPACK. + + int n = m_eivec.cols(); + int low = 0; + int high = n-1; + + for (int m = low+1; m <= high-1; ++m) + { + // Scale column. + RealScalar scale = matH.block(m, m-1, high-m+1, 1).cwise().abs().sum(); + if (scale != 0.0) + { + // Compute Householder transformation. + RealScalar h = 0.0; + // FIXME could be rewritten, but this one looks better wrt cache + for (int i = high; i >= m; i--) + { + ort.coeffRef(i) = matH.coeff(i,m-1)/scale; + h += ort.coeff(i) * ort.coeff(i); + } + RealScalar g = ei_sqrt(h); + if (ort.coeff(m) > 0) + g = -g; + h = h - ort.coeff(m) * g; + ort.coeffRef(m) = ort.coeff(m) - g; + + // Apply Householder similarity transformation + // H = (I-u*u'/h)*H*(I-u*u')/h) + int bSize = high-m+1; + matH.block(m, m, bSize, n-m).noalias() -= ((ort.segment(m, bSize)/h) + * (ort.segment(m, bSize).transpose() * matH.block(m, m, bSize, n-m))); + + matH.block(0, m, high+1, bSize).noalias() -= ((matH.block(0, m, high+1, bSize) * ort.segment(m, bSize)) + * (ort.segment(m, bSize)/h).transpose()); + + ort.coeffRef(m) = scale*ort.coeff(m); + matH.coeffRef(m,m-1) = scale*g; + } + } + + // Accumulate transformations (Algol's ortran). + m_eivec.setIdentity(); + + for (int m = high-1; m >= low+1; m--) + { + if (matH.coeff(m,m-1) != 0.0) + { + ort.segment(m+1, high-m) = matH.col(m-1).segment(m+1, high-m); + + int bSize = high-m+1; + m_eivec.block(m, m, bSize, bSize).noalias() += ( (ort.segment(m, bSize) / (matH.coeff(m,m-1) * ort.coeff(m))) + * (ort.segment(m, bSize).transpose() * m_eivec.block(m, m, bSize, bSize)) ); + } + } +} + +// Complex scalar division. +template +std::complex cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi) +{ + Scalar r,d; + if (ei_abs(yr) > ei_abs(yi)) + { + r = yi/yr; + d = yr + r*yi; + return std::complex((xr + r*xi)/d, (xi - r*xr)/d); + } + else + { + r = yr/yi; + d = yi + r*yr; + return std::complex((r*xr + xi)/d, (r*xi - xr)/d); + } +} + + +// Nonsymmetric reduction from Hessenberg to real Schur form. +template +void EigenSolver::hqr2(MatrixType& matH) +{ + // This is derived from the Algol procedure hqr2, + // by Martin and Wilkinson, Handbook for Auto. Comp., + // Vol.ii-Linear Algebra, and the corresponding + // Fortran subroutine in EISPACK. + + // Initialize + int nn = m_eivec.cols(); + int n = nn-1; + int low = 0; + int high = nn-1; + Scalar eps = ei_pow(Scalar(2),ei_is_same_type::ret ? Scalar(-23) : Scalar(-52)); + Scalar exshift = 0.0; + Scalar p=0,q=0,r=0,s=0,z=0,t,w,x,y; + + // Store roots isolated by balanc and compute matrix norm + // FIXME to be efficient the following would requires a triangular reduxion code + // Scalar norm = matH.upper().cwise().abs().sum() + matH.corner(BottomLeft,n,n).diagonal().cwise().abs().sum(); + Scalar norm = 0.0; + for (int j = 0; j < nn; ++j) + { + // FIXME what's the purpose of the following since the condition is always false + if ((j < low) || (j > high)) + { + m_eivalues.coeffRef(j) = Complex(matH.coeff(j,j), 0.0); + } + norm += matH.row(j).segment(std::max(j-1,0), nn-std::max(j-1,0)).cwise().abs().sum(); + } + + // Outer loop over eigenvalue index + int iter = 0; + while (n >= low) + { + // Look for single small sub-diagonal element + int l = n; + while (l > low) + { + s = ei_abs(matH.coeff(l-1,l-1)) + ei_abs(matH.coeff(l,l)); + if (s == 0.0) + s = norm; + if (ei_abs(matH.coeff(l,l-1)) < eps * s) + break; + l--; + } + + // Check for convergence + // One root found + if (l == n) + { + matH.coeffRef(n,n) = matH.coeff(n,n) + exshift; + m_eivalues.coeffRef(n) = Complex(matH.coeff(n,n), 0.0); + n--; + iter = 0; + } + else if (l == n-1) // Two roots found + { + w = matH.coeff(n,n-1) * matH.coeff(n-1,n); + p = (matH.coeff(n-1,n-1) - matH.coeff(n,n)) * Scalar(0.5); + q = p * p + w; + z = ei_sqrt(ei_abs(q)); + matH.coeffRef(n,n) = matH.coeff(n,n) + exshift; + matH.coeffRef(n-1,n-1) = matH.coeff(n-1,n-1) + exshift; + x = matH.coeff(n,n); + + // Scalar pair + if (q >= 0) + { + if (p >= 0) + z = p + z; + else + z = p - z; + + m_eivalues.coeffRef(n-1) = Complex(x + z, 0.0); + m_eivalues.coeffRef(n) = Complex(z!=0.0 ? x - w / z : m_eivalues.coeff(n-1).real(), 0.0); + + x = matH.coeff(n,n-1); + s = ei_abs(x) + ei_abs(z); + p = x / s; + q = z / s; + r = ei_sqrt(p * p+q * q); + p = p / r; + q = q / r; + + // Row modification + for (int j = n-1; j < nn; ++j) + { + z = matH.coeff(n-1,j); + matH.coeffRef(n-1,j) = q * z + p * matH.coeff(n,j); + matH.coeffRef(n,j) = q * matH.coeff(n,j) - p * z; + } + + // Column modification + for (int i = 0; i <= n; ++i) + { + z = matH.coeff(i,n-1); + matH.coeffRef(i,n-1) = q * z + p * matH.coeff(i,n); + matH.coeffRef(i,n) = q * matH.coeff(i,n) - p * z; + } + + // Accumulate transformations + for (int i = low; i <= high; ++i) + { + z = m_eivec.coeff(i,n-1); + m_eivec.coeffRef(i,n-1) = q * z + p * m_eivec.coeff(i,n); + m_eivec.coeffRef(i,n) = q * m_eivec.coeff(i,n) - p * z; + } + } + else // Complex pair + { + m_eivalues.coeffRef(n-1) = Complex(x + p, z); + m_eivalues.coeffRef(n) = Complex(x + p, -z); + } + n = n - 2; + iter = 0; + } + else // No convergence yet + { + // Form shift + x = matH.coeff(n,n); + y = 0.0; + w = 0.0; + if (l < n) + { + y = matH.coeff(n-1,n-1); + w = matH.coeff(n,n-1) * matH.coeff(n-1,n); + } + + // Wilkinson's original ad hoc shift + if (iter == 10) + { + exshift += x; + for (int i = low; i <= n; ++i) + matH.coeffRef(i,i) -= x; + s = ei_abs(matH.coeff(n,n-1)) + ei_abs(matH.coeff(n-1,n-2)); + x = y = Scalar(0.75) * s; + w = Scalar(-0.4375) * s * s; + } + + // MATLAB's new ad hoc shift + if (iter == 30) + { + s = Scalar((y - x) / 2.0); + s = s * s + w; + if (s > 0) + { + s = ei_sqrt(s); + if (y < x) + s = -s; + s = Scalar(x - w / ((y - x) / 2.0 + s)); + for (int i = low; i <= n; ++i) + matH.coeffRef(i,i) -= s; + exshift += s; + x = y = w = Scalar(0.964); + } + } + + iter = iter + 1; // (Could check iteration count here.) + + // Look for two consecutive small sub-diagonal elements + int m = n-2; + while (m >= l) + { + z = matH.coeff(m,m); + r = x - z; + s = y - z; + p = (r * s - w) / matH.coeff(m+1,m) + matH.coeff(m,m+1); + q = matH.coeff(m+1,m+1) - z - r - s; + r = matH.coeff(m+2,m+1); + s = ei_abs(p) + ei_abs(q) + ei_abs(r); + p = p / s; + q = q / s; + r = r / s; + if (m == l) { + break; + } + if (ei_abs(matH.coeff(m,m-1)) * (ei_abs(q) + ei_abs(r)) < + eps * (ei_abs(p) * (ei_abs(matH.coeff(m-1,m-1)) + ei_abs(z) + + ei_abs(matH.coeff(m+1,m+1))))) + { + break; + } + m--; + } + + for (int i = m+2; i <= n; ++i) + { + matH.coeffRef(i,i-2) = 0.0; + if (i > m+2) + matH.coeffRef(i,i-3) = 0.0; + } + + // Double QR step involving rows l:n and columns m:n + for (int k = m; k <= n-1; ++k) + { + int notlast = (k != n-1); + if (k != m) { + p = matH.coeff(k,k-1); + q = matH.coeff(k+1,k-1); + r = notlast ? matH.coeff(k+2,k-1) : Scalar(0); + x = ei_abs(p) + ei_abs(q) + ei_abs(r); + if (x != 0.0) + { + p = p / x; + q = q / x; + r = r / x; + } + } + + if (x == 0.0) + break; + + s = ei_sqrt(p * p + q * q + r * r); + + if (p < 0) + s = -s; + + if (s != 0) + { + if (k != m) + matH.coeffRef(k,k-1) = -s * x; + else if (l != m) + matH.coeffRef(k,k-1) = -matH.coeff(k,k-1); + + p = p + s; + x = p / s; + y = q / s; + z = r / s; + q = q / p; + r = r / p; + + // Row modification + for (int j = k; j < nn; ++j) + { + p = matH.coeff(k,j) + q * matH.coeff(k+1,j); + if (notlast) + { + p = p + r * matH.coeff(k+2,j); + matH.coeffRef(k+2,j) = matH.coeff(k+2,j) - p * z; + } + matH.coeffRef(k,j) = matH.coeff(k,j) - p * x; + matH.coeffRef(k+1,j) = matH.coeff(k+1,j) - p * y; + } + + // Column modification + for (int i = 0; i <= std::min(n,k+3); ++i) + { + p = x * matH.coeff(i,k) + y * matH.coeff(i,k+1); + if (notlast) + { + p = p + z * matH.coeff(i,k+2); + matH.coeffRef(i,k+2) = matH.coeff(i,k+2) - p * r; + } + matH.coeffRef(i,k) = matH.coeff(i,k) - p; + matH.coeffRef(i,k+1) = matH.coeff(i,k+1) - p * q; + } + + // Accumulate transformations + for (int i = low; i <= high; ++i) + { + p = x * m_eivec.coeff(i,k) + y * m_eivec.coeff(i,k+1); + if (notlast) + { + p = p + z * m_eivec.coeff(i,k+2); + m_eivec.coeffRef(i,k+2) = m_eivec.coeff(i,k+2) - p * r; + } + m_eivec.coeffRef(i,k) = m_eivec.coeff(i,k) - p; + m_eivec.coeffRef(i,k+1) = m_eivec.coeff(i,k+1) - p * q; + } + } // (s != 0) + } // k loop + } // check convergence + } // while (n >= low) + + // Backsubstitute to find vectors of upper triangular form + if (norm == 0.0) + { + return; + } + + for (n = nn-1; n >= 0; n--) + { + p = m_eivalues.coeff(n).real(); + q = m_eivalues.coeff(n).imag(); + + // Scalar vector + if (q == 0) + { + int l = n; + matH.coeffRef(n,n) = 1.0; + for (int i = n-1; i >= 0; i--) + { + w = matH.coeff(i,i) - p; + r = matH.row(i).segment(l,n-l+1).dot(matH.col(n).segment(l, n-l+1)); + + if (m_eivalues.coeff(i).imag() < 0.0) + { + z = w; + s = r; + } + else + { + l = i; + if (m_eivalues.coeff(i).imag() == 0.0) + { + if (w != 0.0) + matH.coeffRef(i,n) = -r / w; + else + matH.coeffRef(i,n) = -r / (eps * norm); + } + else // Solve real equations + { + x = matH.coeff(i,i+1); + y = matH.coeff(i+1,i); + q = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag(); + t = (x * s - z * r) / q; + matH.coeffRef(i,n) = t; + if (ei_abs(x) > ei_abs(z)) + matH.coeffRef(i+1,n) = (-r - w * t) / x; + else + matH.coeffRef(i+1,n) = (-s - y * t) / z; + } + + // Overflow control + t = ei_abs(matH.coeff(i,n)); + if ((eps * t) * t > 1) + matH.col(n).end(nn-i) /= t; + } + } + } + else if (q < 0) // Complex vector + { + std::complex cc; + int l = n-1; + + // Last vector component imaginary so matrix is triangular + if (ei_abs(matH.coeff(n,n-1)) > ei_abs(matH.coeff(n-1,n))) + { + matH.coeffRef(n-1,n-1) = q / matH.coeff(n,n-1); + matH.coeffRef(n-1,n) = -(matH.coeff(n,n) - p) / matH.coeff(n,n-1); + } + else + { + cc = cdiv(0.0,-matH.coeff(n-1,n),matH.coeff(n-1,n-1)-p,q); + matH.coeffRef(n-1,n-1) = ei_real(cc); + matH.coeffRef(n-1,n) = ei_imag(cc); + } + matH.coeffRef(n,n-1) = 0.0; + matH.coeffRef(n,n) = 1.0; + for (int i = n-2; i >= 0; i--) + { + Scalar ra,sa,vr,vi; + ra = matH.row(i).segment(l, n-l+1).dot(matH.col(n-1).segment(l, n-l+1)); + sa = matH.row(i).segment(l, n-l+1).dot(matH.col(n).segment(l, n-l+1)); + w = matH.coeff(i,i) - p; + + if (m_eivalues.coeff(i).imag() < 0.0) + { + z = w; + r = ra; + s = sa; + } + else + { + l = i; + if (m_eivalues.coeff(i).imag() == 0) + { + cc = cdiv(-ra,-sa,w,q); + matH.coeffRef(i,n-1) = ei_real(cc); + matH.coeffRef(i,n) = ei_imag(cc); + } + else + { + // Solve complex equations + x = matH.coeff(i,i+1); + y = matH.coeff(i+1,i); + vr = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag() - q * q; + vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q; + if ((vr == 0.0) && (vi == 0.0)) + vr = eps * norm * (ei_abs(w) + ei_abs(q) + ei_abs(x) + ei_abs(y) + ei_abs(z)); + + cc= cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi); + matH.coeffRef(i,n-1) = ei_real(cc); + matH.coeffRef(i,n) = ei_imag(cc); + if (ei_abs(x) > (ei_abs(z) + ei_abs(q))) + { + matH.coeffRef(i+1,n-1) = (-ra - w * matH.coeff(i,n-1) + q * matH.coeff(i,n)) / x; + matH.coeffRef(i+1,n) = (-sa - w * matH.coeff(i,n) - q * matH.coeff(i,n-1)) / x; + } + else + { + cc = cdiv(-r-y*matH.coeff(i,n-1),-s-y*matH.coeff(i,n),z,q); + matH.coeffRef(i+1,n-1) = ei_real(cc); + matH.coeffRef(i+1,n) = ei_imag(cc); + } + } + + // Overflow control + t = std::max(ei_abs(matH.coeff(i,n-1)),ei_abs(matH.coeff(i,n))); + if ((eps * t) * t > 1) + matH.block(i, n-1, nn-i, 2) /= t; + + } + } + } + } + + // Vectors of isolated roots + for (int i = 0; i < nn; ++i) + { + // FIXME again what's the purpose of this test ? + // in this algo low==0 and high==nn-1 !! + if (i < low || i > high) + { + m_eivec.row(i).end(nn-i) = matH.row(i).end(nn-i); + } + } + + // Back transformation to get eigenvectors of original matrix + int bRows = high-low+1; + for (int j = nn-1; j >= low; j--) + { + int bSize = std::min(j,high)-low+1; + m_eivec.col(j).segment(low, bRows) = (m_eivec.block(low, low, bRows, bSize) * matH.col(j).segment(low, bSize)); + } +} + +#endif // EIGEN_EIGENSOLVER_H diff --git a/Eigen/src/Eigenvalues/HessenbergDecomposition.h b/Eigen/src/Eigenvalues/HessenbergDecomposition.h new file mode 100644 index 000000000..b1e21d4ee --- /dev/null +++ b/Eigen/src/Eigenvalues/HessenbergDecomposition.h @@ -0,0 +1,206 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_HESSENBERGDECOMPOSITION_H +#define EIGEN_HESSENBERGDECOMPOSITION_H + +/** \eigenvalues_module \ingroup Eigenvalues_Module + * \nonstableyet + * + * \class HessenbergDecomposition + * + * \brief Reduces a squared matrix to an Hessemberg form + * + * \param MatrixType the type of the matrix of which we are computing the Hessenberg decomposition + * + * This class performs an Hessenberg decomposition of a matrix \f$ A \f$ such that: + * \f$ A = Q H Q^* \f$ where \f$ Q \f$ is unitary and \f$ H \f$ a Hessenberg matrix. + * + * \sa class Tridiagonalization, class Qr + */ +template class HessenbergDecomposition +{ + public: + + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + + enum { + Size = MatrixType::RowsAtCompileTime, + SizeMinusOne = MatrixType::RowsAtCompileTime==Dynamic + ? Dynamic + : MatrixType::RowsAtCompileTime-1 + }; + + typedef Matrix CoeffVectorType; + typedef Matrix DiagonalType; + typedef Matrix SubDiagonalType; + + typedef typename NestByValue >::RealReturnType DiagonalReturnType; + + typedef typename NestByValue >,0 > >::RealReturnType SubDiagonalReturnType; + + /** This constructor initializes a HessenbergDecomposition object for + * further use with HessenbergDecomposition::compute() + */ + HessenbergDecomposition(int size = Size==Dynamic ? 2 : Size) + : m_matrix(size,size), m_hCoeffs(size-1) + {} + + HessenbergDecomposition(const MatrixType& matrix) + : m_matrix(matrix), + m_hCoeffs(matrix.cols()-1) + { + _compute(m_matrix, m_hCoeffs); + } + + /** Computes or re-compute the Hessenberg decomposition for the matrix \a matrix. + * + * This method allows to re-use the allocated data. + */ + void compute(const MatrixType& matrix) + { + m_matrix = matrix; + m_hCoeffs.resize(matrix.rows()-1,1); + _compute(m_matrix, m_hCoeffs); + } + + /** \returns the householder coefficients allowing to + * reconstruct the matrix Q from the packed data. + * + * \sa packedMatrix() + */ + CoeffVectorType householderCoefficients() const { return m_hCoeffs; } + + /** \returns the internal result of the decomposition. + * + * The returned matrix contains the following information: + * - the upper part and lower sub-diagonal represent the Hessenberg matrix H + * - the rest of the lower part contains the Householder vectors that, combined with + * Householder coefficients returned by householderCoefficients(), + * allows to reconstruct the matrix Q as follow: + * Q = H_{N-1} ... H_1 H_0 + * where the matrices H are the Householder transformation: + * H_i = (I - h_i * v_i * v_i') + * where h_i == householderCoefficients()[i] and v_i is a Householder vector: + * v_i = [ 0, ..., 0, 1, M(i+2,i), ..., M(N-1,i) ] + * + * See LAPACK for further details on this packed storage. + */ + const MatrixType& packedMatrix(void) const { return m_matrix; } + + MatrixType matrixQ() const; + MatrixType matrixH() const; + + private: + + static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs); + + protected: + MatrixType m_matrix; + CoeffVectorType m_hCoeffs; +}; + +#ifndef EIGEN_HIDE_HEAVY_CODE + +/** \internal + * Performs a tridiagonal decomposition of \a matA in place. + * + * \param matA the input selfadjoint matrix + * \param hCoeffs returned Householder coefficients + * + * The result is written in the lower triangular part of \a matA. + * + * Implemented from Golub's "Matrix Computations", algorithm 8.3.1. + * + * \sa packedMatrix() + */ +template +void HessenbergDecomposition::_compute(MatrixType& matA, CoeffVectorType& hCoeffs) +{ + assert(matA.rows()==matA.cols()); + int n = matA.rows(); + Matrix temp(n); + for (int i = 0; i +typename HessenbergDecomposition::MatrixType +HessenbergDecomposition::matrixQ() const +{ + int n = m_matrix.rows(); + MatrixType matQ = MatrixType::Identity(n,n); + Matrix temp(n); + for (int i = n-2; i>=0; i--) + { + matQ.corner(BottomRight,n-i-1,n-i-1) + .applyHouseholderOnTheLeft(m_matrix.col(i).end(n-i-2), ei_conj(m_hCoeffs.coeff(i)), &temp.coeffRef(0,0)); + } + return matQ; +} + +#endif // EIGEN_HIDE_HEAVY_CODE + +/** constructs and returns the matrix H. + * Note that the matrix H is equivalent to the upper part of the packed matrix + * (including the lower sub-diagonal). Therefore, it might be often sufficient + * to directly use the packed matrix instead of creating a new one. + */ +template +typename HessenbergDecomposition::MatrixType +HessenbergDecomposition::matrixH() const +{ + // FIXME should this function (and other similar) rather take a matrix as argument + // and fill it (to avoid temporaries) + int n = m_matrix.rows(); + MatrixType matH = m_matrix; + if (n>2) + matH.corner(BottomLeft,n-2, n-2).template triangularView().setZero(); + return matH; +} + +#endif // EIGEN_HESSENBERGDECOMPOSITION_H diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h new file mode 100644 index 000000000..84856aa66 --- /dev/null +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -0,0 +1,366 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_SELFADJOINTEIGENSOLVER_H +#define EIGEN_SELFADJOINTEIGENSOLVER_H + +/** \eigenvalues_module \ingroup Eigenvalues_Module + * \nonstableyet + * + * \class SelfAdjointEigenSolver + * + * \brief Eigen values/vectors solver for selfadjoint matrix + * + * \param MatrixType the type of the matrix of which we are computing the eigen decomposition + * + * \note MatrixType must be an actual Matrix type, it can't be an expression type. + * + * \sa MatrixBase::eigenvalues(), class EigenSolver + */ +template class SelfAdjointEigenSolver +{ + public: + + enum {Size = _MatrixType::RowsAtCompileTime }; + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef std::complex Complex; + typedef Matrix RealVectorType; + typedef Matrix RealVectorTypeX; + typedef Tridiagonalization TridiagonalizationType; +// typedef typename TridiagonalizationType::TridiagonalMatrixType TridiagonalMatrixType; + + SelfAdjointEigenSolver() + : m_eivec(int(Size), int(Size)), + m_eivalues(int(Size)) + { + ei_assert(Size!=Dynamic); + } + + SelfAdjointEigenSolver(int size) + : m_eivec(size, size), + m_eivalues(size) + {} + + /** Constructors computing the eigenvalues of the selfadjoint matrix \a matrix, + * as well as the eigenvectors if \a computeEigenvectors is true. + * + * \sa compute(MatrixType,bool), SelfAdjointEigenSolver(MatrixType,MatrixType,bool) + */ + SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true) + : m_eivec(matrix.rows(), matrix.cols()), + m_eivalues(matrix.cols()) + { + compute(matrix, computeEigenvectors); + } + + /** Constructors computing the eigenvalues of the generalized eigen problem + * \f$ Ax = lambda B x \f$ with \a matA the selfadjoint matrix \f$ A \f$ + * and \a matB the positive definite matrix \f$ B \f$ . The eigenvectors + * are computed if \a computeEigenvectors is true. + * + * \sa compute(MatrixType,MatrixType,bool), SelfAdjointEigenSolver(MatrixType,bool) + */ + SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true) + : m_eivec(matA.rows(), matA.cols()), + m_eivalues(matA.cols()) + { + compute(matA, matB, computeEigenvectors); + } + + SelfAdjointEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true); + + SelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true); + + /** \returns the computed eigen vectors as a matrix of column vectors */ + MatrixType eigenvectors(void) const + { + #ifndef NDEBUG + ei_assert(m_eigenvectorsOk); + #endif + return m_eivec; + } + + /** \returns the computed eigen values */ + RealVectorType eigenvalues(void) const { return m_eivalues; } + + /** \returns the positive square root of the matrix + * + * \note the matrix itself must be positive in order for this to make sense. + */ + MatrixType operatorSqrt() const + { + return m_eivec * m_eivalues.cwise().sqrt().asDiagonal() * m_eivec.adjoint(); + } + + /** \returns the positive inverse square root of the matrix + * + * \note the matrix itself must be positive definite in order for this to make sense. + */ + MatrixType operatorInverseSqrt() const + { + return m_eivec * m_eivalues.cwise().inverse().cwise().sqrt().asDiagonal() * m_eivec.adjoint(); + } + + + protected: + MatrixType m_eivec; + RealVectorType m_eivalues; + #ifndef NDEBUG + bool m_eigenvectorsOk; + #endif +}; + +#ifndef EIGEN_HIDE_HEAVY_CODE + +/** \internal + * + * \eigenvalues_module \ingroup Eigenvalues_Module + * + * Performs a QR step on a tridiagonal symmetric matrix represented as a + * pair of two vectors \a diag and \a subdiag. + * + * \param matA the input selfadjoint matrix + * \param hCoeffs returned Householder coefficients + * + * For compilation efficiency reasons, this procedure does not use eigen expression + * for its arguments. + * + * Implemented from Golub's "Matrix Computations", algorithm 8.3.2: + * "implicit symmetric QR step with Wilkinson shift" + */ +template +static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int start, int end, Scalar* matrixQ, int n); + +/** Computes the eigenvalues of the selfadjoint matrix \a matrix, + * as well as the eigenvectors if \a computeEigenvectors is true. + * + * \sa SelfAdjointEigenSolver(MatrixType,bool), compute(MatrixType,MatrixType,bool) + */ +template +SelfAdjointEigenSolver& SelfAdjointEigenSolver::compute(const MatrixType& matrix, bool computeEigenvectors) +{ + #ifndef NDEBUG + m_eigenvectorsOk = computeEigenvectors; + #endif + assert(matrix.cols() == matrix.rows()); + int n = matrix.cols(); + m_eivalues.resize(n,1); + + if(n==1) + { + m_eivalues.coeffRef(0,0) = ei_real(matrix.coeff(0,0)); + m_eivec.setOnes(); + return *this; + } + + m_eivec = matrix; + + // FIXME, should tridiag be a local variable of this function or an attribute of SelfAdjointEigenSolver ? + // the latter avoids multiple memory allocation when the same SelfAdjointEigenSolver is used multiple times... + // (same for diag and subdiag) + RealVectorType& diag = m_eivalues; + typename TridiagonalizationType::SubDiagonalType subdiag(n-1); + TridiagonalizationType::decomposeInPlace(m_eivec, diag, subdiag, computeEigenvectors); + + int end = n-1; + int start = 0; + while (end>0) + { + for (int i = start; i0 && subdiag[end-1]==0) + end--; + if (end<=0) + break; + start = end - 1; + while (start>0 && subdiag[start-1]!=0) + start--; + + ei_tridiagonal_qr_step(diag.data(), subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n); + } + + // Sort eigenvalues and corresponding vectors. + // TODO make the sort optional ? + // TODO use a better sort algorithm !! + for (int i = 0; i < n-1; ++i) + { + int k; + m_eivalues.segment(i,n-i).minCoeff(&k); + if (k > 0) + { + std::swap(m_eivalues[i], m_eivalues[k+i]); + m_eivec.col(i).swap(m_eivec.col(k+i)); + } + } + return *this; +} + +/** Computes the eigenvalues of the generalized eigen problem + * \f$ Ax = lambda B x \f$ with \a matA the selfadjoint matrix \f$ A \f$ + * and \a matB the positive definite matrix \f$ B \f$ . The eigenvectors + * are computed if \a computeEigenvectors is true. + * + * \sa SelfAdjointEigenSolver(MatrixType,MatrixType,bool), compute(MatrixType,bool) + */ +template +SelfAdjointEigenSolver& SelfAdjointEigenSolver:: +compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors) +{ + ei_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows()); + + // Compute the cholesky decomposition of matB = L L' + LLT cholB(matB); + + // compute C = inv(L) A inv(L') + MatrixType matC = matA; + cholB.matrixL().solveInPlace(matC); + // FIXME since we currently do not support A * inv(L'), let's do (inv(L) A')' : + matC.adjointInPlace(); + cholB.matrixL().solveInPlace(matC); + matC.adjointInPlace(); + // this version works too: +// matC = matC.transpose(); +// cholB.matrixL().conjugate().template marked().solveTriangularInPlace(matC); +// matC = matC.transpose(); + // FIXME: this should work: (currently it only does for small matrices) +// Transpose trMatC(matC); +// cholB.matrixL().conjugate().eval().template marked().solveTriangularInPlace(trMatC); + + compute(matC, computeEigenvectors); + + if (computeEigenvectors) + { + // transform back the eigen vectors: evecs = inv(U) * evecs + cholB.matrixU().solveInPlace(m_eivec); + for (int i=0; i +inline Matrix::Scalar>::Real, ei_traits::ColsAtCompileTime, 1> +MatrixBase::eigenvalues() const +{ + ei_assert(Flags&SelfAdjointBit); + return SelfAdjointEigenSolver(eval(),false).eigenvalues(); +} + +template +struct ei_operatorNorm_selector +{ + static inline typename NumTraits::Scalar>::Real + operatorNorm(const MatrixBase& m) + { + // FIXME if it is really guaranteed that the eigenvalues are already sorted, + // then we don't need to compute a maxCoeff() here, comparing the 1st and last ones is enough. + return m.eigenvalues().cwise().abs().maxCoeff(); + } +}; + +template struct ei_operatorNorm_selector +{ + static inline typename NumTraits::Scalar>::Real + operatorNorm(const MatrixBase& m) + { + typename Derived::PlainMatrixType m_eval(m); + // FIXME if it is really guaranteed that the eigenvalues are already sorted, + // then we don't need to compute a maxCoeff() here, comparing the 1st and last ones is enough. + return ei_sqrt( + (m_eval*m_eval.adjoint()) + .template marked() + .eigenvalues() + .maxCoeff() + ); + } +}; + +/** \eigenvalues_module + * + * \returns the matrix norm of this matrix. + */ +template +inline typename NumTraits::Scalar>::Real +MatrixBase::operatorNorm() const +{ + return ei_operatorNorm_selector + ::operatorNorm(derived()); +} + +#ifndef EIGEN_EXTERN_INSTANTIATIONS +template +static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int start, int end, Scalar* matrixQ, int n) +{ + RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5); + RealScalar e2 = ei_abs2(subdiag[end-1]); + RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * ei_sqrt(td*td + e2)); + RealScalar x = diag[start] - mu; + RealScalar z = subdiag[start]; + + for (int k = start; k < end; ++k) + { + PlanarRotation rot; + rot.makeGivens(x, z); + + // do T = G' T G + RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k]; + RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k+1]; + + diag[k] = rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k+1]); + diag[k+1] = rot.s() * sdk + rot.c() * dkp1; + subdiag[k] = rot.c() * sdk - rot.s() * dkp1; + + if (k > start) + subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z; + + x = subdiag[k]; + + if (k < end - 1) + { + z = -rot.s() * subdiag[k+1]; + subdiag[k + 1] = rot.c() * subdiag[k+1]; + } + + // apply the givens rotation to the unit matrix Q = Q * G + if (matrixQ) + { + Map > q(matrixQ,n,n); + q.applyOnTheRight(k,k+1,rot); + } + } +} +#endif + +#endif // EIGEN_SELFADJOINTEIGENSOLVER_H diff --git a/Eigen/src/Eigenvalues/Tridiagonalization.h b/Eigen/src/Eigenvalues/Tridiagonalization.h new file mode 100644 index 000000000..5f891bfa6 --- /dev/null +++ b/Eigen/src/Eigenvalues/Tridiagonalization.h @@ -0,0 +1,317 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_TRIDIAGONALIZATION_H +#define EIGEN_TRIDIAGONALIZATION_H + +/** \eigenvalues_module \ingroup Eigenvalues_Module + * \nonstableyet + * + * \class Tridiagonalization + * + * \brief Trigiagonal decomposition of a selfadjoint matrix + * + * \param MatrixType the type of the matrix of which we are performing the tridiagonalization + * + * This class performs a tridiagonal decomposition of a selfadjoint matrix \f$ A \f$ such that: + * \f$ A = Q T Q^* \f$ where \f$ Q \f$ is unitary and \f$ T \f$ a real symmetric tridiagonal matrix. + * + * \sa MatrixBase::tridiagonalize() + */ +template class Tridiagonalization +{ + public: + + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef typename ei_packet_traits::type Packet; + + enum { + Size = MatrixType::RowsAtCompileTime, + SizeMinusOne = MatrixType::RowsAtCompileTime==Dynamic + ? Dynamic + : MatrixType::RowsAtCompileTime-1, + PacketSize = ei_packet_traits::size + }; + + typedef Matrix CoeffVectorType; + typedef Matrix DiagonalType; + typedef Matrix SubDiagonalType; + + typedef typename ei_meta_if::IsComplex, + typename NestByValue >::RealReturnType, + Diagonal + >::ret DiagonalReturnType; + + typedef typename ei_meta_if::IsComplex, + typename NestByValue >,0 > >::RealReturnType, + Diagonal< + NestByValue >,0 > + >::ret SubDiagonalReturnType; + + /** This constructor initializes a Tridiagonalization object for + * further use with Tridiagonalization::compute() + */ + Tridiagonalization(int size = Size==Dynamic ? 2 : Size) + : m_matrix(size,size), m_hCoeffs(size-1) + {} + + Tridiagonalization(const MatrixType& matrix) + : m_matrix(matrix), m_hCoeffs(matrix.cols()-1) + { + _compute(m_matrix, m_hCoeffs); + } + + /** Computes or re-compute the tridiagonalization for the matrix \a matrix. + * + * This method allows to re-use the allocated data. + */ + void compute(const MatrixType& matrix) + { + m_matrix = matrix; + m_hCoeffs.resize(matrix.rows()-1, 1); + _compute(m_matrix, m_hCoeffs); + } + + /** \returns the householder coefficients allowing to + * reconstruct the matrix Q from the packed data. + * + * \sa packedMatrix() + */ + inline CoeffVectorType householderCoefficients(void) const { return m_hCoeffs; } + + /** \returns the internal result of the decomposition. + * + * The returned matrix contains the following information: + * - the strict upper part is equal to the input matrix A + * - the diagonal and lower sub-diagonal represent the tridiagonal symmetric matrix (real). + * - the rest of the lower part contains the Householder vectors that, combined with + * Householder coefficients returned by householderCoefficients(), + * allows to reconstruct the matrix Q as follow: + * Q = H_{N-1} ... H_1 H_0 + * where the matrices H are the Householder transformations: + * H_i = (I - h_i * v_i * v_i') + * where h_i == householderCoefficients()[i] and v_i is a Householder vector: + * v_i = [ 0, ..., 0, 1, M(i+2,i), ..., M(N-1,i) ] + * + * See LAPACK for further details on this packed storage. + */ + inline const MatrixType& packedMatrix(void) const { return m_matrix; } + + MatrixType matrixQ() const; + template void matrixQInPlace(MatrixBase* q) const; + MatrixType matrixT() const; + const DiagonalReturnType diagonal(void) const; + const SubDiagonalReturnType subDiagonal(void) const; + + static void decomposeInPlace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ = true); + + static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs); + + protected: + + static void _decomposeInPlace3x3(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ = true); + + MatrixType m_matrix; + CoeffVectorType m_hCoeffs; +}; + +/** \returns an expression of the diagonal vector */ +template +const typename Tridiagonalization::DiagonalReturnType +Tridiagonalization::diagonal(void) const +{ + return m_matrix.diagonal().nestByValue(); +} + +/** \returns an expression of the sub-diagonal vector */ +template +const typename Tridiagonalization::SubDiagonalReturnType +Tridiagonalization::subDiagonal(void) const +{ + int n = m_matrix.rows(); + return Block(m_matrix, 1, 0, n-1,n-1) + .nestByValue().diagonal().nestByValue(); +} + +/** constructs and returns the tridiagonal matrix T. + * Note that the matrix T is equivalent to the diagonal and sub-diagonal of the packed matrix. + * Therefore, it might be often sufficient to directly use the packed matrix, or the vector + * expressions returned by diagonal() and subDiagonal() instead of creating a new matrix. + */ +template +typename Tridiagonalization::MatrixType +Tridiagonalization::matrixT(void) const +{ + // FIXME should this function (and other similar ones) rather take a matrix as argument + // and fill it ? (to avoid temporaries) + int n = m_matrix.rows(); + MatrixType matT = m_matrix; + matT.corner(TopRight,n-1, n-1).diagonal() = subDiagonal().template cast().conjugate(); + if (n>2) + { + matT.corner(TopRight,n-2, n-2).template triangularView().setZero(); + matT.corner(BottomLeft,n-2, n-2).template triangularView().setZero(); + } + return matT; +} + +#ifndef EIGEN_HIDE_HEAVY_CODE + +/** \internal + * Performs a tridiagonal decomposition of \a matA in place. + * + * \param matA the input selfadjoint matrix + * \param hCoeffs returned Householder coefficients + * + * The result is written in the lower triangular part of \a matA. + * + * Implemented from Golub's "Matrix Computations", algorithm 8.3.1. + * + * \sa packedMatrix() + */ +template +void Tridiagonalization::_compute(MatrixType& matA, CoeffVectorType& hCoeffs) +{ + assert(matA.rows()==matA.cols()); + int n = matA.rows(); + Matrix aux(n); + for (int i = 0; i() + * (ei_conj(h) * matA.col(i).end(remainingSize))); + + hCoeffs.end(n-i-1) += (ei_conj(h)*Scalar(-0.5)*(hCoeffs.end(remainingSize).dot(matA.col(i).end(remainingSize)))) * matA.col(i).end(n-i-1); + + matA.corner(BottomRight, remainingSize, remainingSize).template selfadjointView() + .rankUpdate(matA.col(i).end(remainingSize), hCoeffs.end(remainingSize), -1); + + matA.col(i).coeffRef(i+1) = beta; + hCoeffs.coeffRef(i) = h; + } +} + +/** reconstructs and returns the matrix Q */ +template +typename Tridiagonalization::MatrixType +Tridiagonalization::matrixQ(void) const +{ + MatrixType matQ; + matrixQInPlace(&matQ); + return matQ; +} + +template +template +void Tridiagonalization::matrixQInPlace(MatrixBase* q) const +{ + QDerived& matQ = q->derived(); + int n = m_matrix.rows(); + matQ = MatrixType::Identity(n,n); + Matrix aux(n); + for (int i = n-2; i>=0; i--) + { + matQ.corner(BottomRight,n-i-1,n-i-1) + .applyHouseholderOnTheLeft(m_matrix.col(i).end(n-i-2), ei_conj(m_hCoeffs.coeff(i)), &aux.coeffRef(0,0)); + } +} + +/** Performs a full decomposition in place */ +template +void Tridiagonalization::decomposeInPlace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ) +{ + int n = mat.rows(); + ei_assert(mat.cols()==n && diag.size()==n && subdiag.size()==n-1); + if (n==3 && (!NumTraits::IsComplex) ) + { + _decomposeInPlace3x3(mat, diag, subdiag, extractQ); + } + else + { + Tridiagonalization tridiag(mat); + diag = tridiag.diagonal(); + subdiag = tridiag.subDiagonal(); + if (extractQ) + tridiag.matrixQInPlace(&mat); + } +} + +/** \internal + * Optimized path for 3x3 matrices. + * Especially useful for plane fitting. + */ +template +void Tridiagonalization::_decomposeInPlace3x3(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ) +{ + diag[0] = ei_real(mat(0,0)); + RealScalar v1norm2 = ei_abs2(mat(0,2)); + if (ei_isMuchSmallerThan(v1norm2, RealScalar(1))) + { + diag[1] = ei_real(mat(1,1)); + diag[2] = ei_real(mat(2,2)); + subdiag[0] = ei_real(mat(0,1)); + subdiag[1] = ei_real(mat(1,2)); + if (extractQ) + mat.setIdentity(); + } + else + { + RealScalar beta = ei_sqrt(ei_abs2(mat(0,1))+v1norm2); + RealScalar invBeta = RealScalar(1)/beta; + Scalar m01 = mat(0,1) * invBeta; + Scalar m02 = mat(0,2) * invBeta; + Scalar q = RealScalar(2)*m01*mat(1,2) + m02*(mat(2,2) - mat(1,1)); + diag[1] = ei_real(mat(1,1) + m02*q); + diag[2] = ei_real(mat(2,2) - m02*q); + subdiag[0] = beta; + subdiag[1] = ei_real(mat(1,2) - m01 * q); + if (extractQ) + { + mat(0,0) = 1; + mat(0,1) = 0; + mat(0,2) = 0; + mat(1,0) = 0; + mat(1,1) = m01; + mat(1,2) = m02; + mat(2,0) = 0; + mat(2,1) = m02; + mat(2,2) = -m01; + } + } +} + +#endif // EIGEN_HIDE_HEAVY_CODE + +#endif // EIGEN_TRIDIAGONALIZATION_H -- cgit v1.2.3