diff --git a/Eigen/QR b/Eigen/QR index cee480594..7e3a0bd8c 100644 --- a/Eigen/QR +++ b/Eigen/QR @@ -6,6 +6,7 @@ namespace Eigen { #include "src/QR/QR.h" +#include "src/QR/Tridiagonalization.h" #include "src/QR/EigenSolver.h" } // namespace Eigen diff --git a/Eigen/src/Core/CwiseUnaryOp.h b/Eigen/src/Core/CwiseUnaryOp.h index 893c8abf6..8715f931e 100644 --- a/Eigen/src/Core/CwiseUnaryOp.h +++ b/Eigen/src/Core/CwiseUnaryOp.h @@ -139,7 +139,7 @@ MatrixBase::cwiseAbs2() const return derived(); } -/** \returns an expression of the complex conjugate of *this. +/** \returns an expression of the complex conjugate of \c *this. * * \sa adjoint() */ template @@ -149,6 +149,16 @@ MatrixBase::conjugate() const return ConjugateReturnType(derived()); } +/** \returns an expression of the real part of \c *this. + * + * \sa adjoint() */ +template +inline const typename MatrixBase::RealReturnType +MatrixBase::real() const +{ + return derived(); +} + /** \returns an expression of *this with the \a Scalar type casted to * \a NewScalar. * diff --git a/Eigen/src/Core/Extract.h b/Eigen/src/Core/Extract.h index ea7240ada..6671e45b2 100755 --- a/Eigen/src/Core/Extract.h +++ b/Eigen/src/Core/Extract.h @@ -77,7 +77,7 @@ template class Extract inline Scalar _coeff(int row, int col) const { if(Flags & LowerTriangularBit ? col>row : row>col) - return (Scalar)0; + return (Flags & SelfAdjointBit) ? ei_conj(m_matrix.coeff(col, row)) : (Scalar)0; if(Flags & UnitDiagBit) return col==row ? (Scalar)1 : m_matrix.coeff(row, col); else if(Flags & ZeroDiagBit) diff --git a/Eigen/src/Core/Functors.h b/Eigen/src/Core/Functors.h index 6a4a1213d..f3ffe10aa 100644 --- a/Eigen/src/Core/Functors.h +++ b/Eigen/src/Core/Functors.h @@ -204,6 +204,19 @@ template struct ei_functor_traits > { enum { Cost = ei_is_same_type::ret ? 0 : NumTraits::AddCost, IsVectorizable = false }; }; +/** \internal + * \brief Template functor to extract the real part of a complex + * + * \sa class CwiseUnaryOp, MatrixBase::real() + */ +template +struct ei_scalar_real_op EIGEN_EMPTY_STRUCT { + typedef typename NumTraits::Real result_type; + inline result_type operator() (const Scalar& a) const { return ei_real(a); } +}; +template +struct ei_functor_traits > +{ enum { Cost = 0, IsVectorizable = false }; }; /** \internal * \brief Template functor to multiply a scalar by a fixed other one diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index ab35b94a6..a29504b33 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -197,6 +197,8 @@ template class MatrixBase : public ArrayBase CwiseUnaryOp, Derived>, Derived& >::ret ConjugateReturnType; + /** the return type of MatrixBase::real() */ + typedef CwiseUnaryOp, Derived> RealReturnType; /** the return type of MatrixBase::adjoint() */ typedef Transpose::type> > AdjointReturnType; @@ -477,6 +479,7 @@ template class MatrixBase : public ArrayBase /// \name Coefficient-wise operations //@{ const ConjugateReturnType conjugate() const; + const RealReturnType real() const; template const CwiseBinaryOp::Scalar>, Derived, OtherDerived> diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index e17563c9a..fd3451922 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -59,8 +59,6 @@ const unsigned int Upper = UpperTriangularBit; const unsigned int StrictlyUpper = UpperTriangularBit | ZeroDiagBit; const unsigned int Lower = LowerTriangularBit; const unsigned int StrictlyLower = LowerTriangularBit | ZeroDiagBit; - -// additional possible values for the Mode parameter of part() const unsigned int SelfAdjoint = SelfAdjointBit; // additional possible values for the Mode parameter of extract() diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index d48d4c325..a9974c38b 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -59,6 +59,7 @@ template struct ei_scalar_product_op; template struct ei_scalar_quotient_op; template struct ei_scalar_opposite_op; template struct ei_scalar_conjugate_op; +template struct ei_scalar_real_op; template struct ei_scalar_abs_op; template struct ei_scalar_abs2_op; template struct ei_scalar_sqrt_op; diff --git a/Eigen/src/QR/Tridiagonalization.h b/Eigen/src/QR/Tridiagonalization.h new file mode 100755 index 000000000..a85068e38 --- /dev/null +++ b/Eigen/src/QR/Tridiagonalization.h @@ -0,0 +1,239 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// 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 + +/** \class Tridiagonalization + * + * \brief Trigiagonal decomposition of a selfadjoint matrix + * + * \param MatrixType the type of the matrix of which we are computing the eigen decomposition + * + * 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 unitatry 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; + + enum {SizeMinusOne = MatrixType::RowsAtCompileTime==Dynamic + ? Dynamic + : MatrixType::RowsAtCompileTime-1}; + + typedef Matrix CoeffVectorType; + + typedef typename NestByValue >::RealReturnType DiagonalType; + + typedef typename NestByValue > > >::RealReturnType SubDiagonalType; + + Tridiagonalization() + {} + + Tridiagonalization(int rows, int cols) + : m_matrix(rows,cols), m_hCoeffs(rows-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); + _compute(m_matrix, m_hCoeffs); + } + + /** \returns the householder coefficients allowing to + * reconstruct the matrix Q from the packed data. + * + * \sa packedMatrix() + */ + 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 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(void) const; + const DiagonalType diagonal(void) const; + const SubDiagonalType subDiagonal(void) const; + + private: + + static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs); + + protected: + MatrixType m_matrix; + CoeffVectorType m_hCoeffs; +}; + + +/** \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: + * + * \sa packedMatrix() + */ +template +void Tridiagonalization::_compute(MatrixType& matA, CoeffVectorType& hCoeffs) +{ + assert(matA.rows()==matA.cols()); + int n = matA.rows(); + for (int i = 0; i(1))) + { + hCoeffs.coeffRef(i) = 0.; + } + else + { + Scalar v0 = matA.col(i).coeff(i+1); + RealScalar beta = ei_sqrt(ei_abs2(v0)+v1norm2); + if (ei_real(v0)>=0.) + beta = -beta; + matA.col(i).end(n-(i+2)) *= (1./(v0-beta)); + matA.col(i).coeffRef(i+1) = beta; + Scalar h = (beta - v0) / beta; + // end of the householder transformation + + // Apply similarity transformation to remaining columns, + // i.e., A = H' A H where H = I - h v v' and v = matA.col(i).end(n-i-1) + + matA.col(i).coeffRef(i+1) = 1; + // let's use the end of hCoeffs to store temporary values + hCoeffs.end(n-i-1) = h * (matA.corner(BottomRight,n-i-1,n-i-1).template extract() + * matA.col(i).end(n-i-1)); + + + hCoeffs.end(n-i-1) += (h * (-0.5) * matA.col(i).end(n-i-1).dot(hCoeffs.end(n-i-1))) + * matA.col(i).end(n-i-1); + + matA.corner(BottomRight,n-i-1,n-i-1).template part() = + matA.corner(BottomRight,n-i-1,n-i-1) - ( + (matA.col(i).end(n-i-1) * hCoeffs.end(n-i-1).adjoint()).lazy() + + (hCoeffs.end(n-i-1) * matA.col(i).end(n-i-1).adjoint()).lazy() ); + // FIXME check that the above expression does follow the lazy path (no temporary and + // only lower products are evaluated) + // FIXME can we avoid to evaluate twice the diagonal products ? + // (in a simple way otherwise it's overkill) + + matA.col(i).coeffRef(i+1) = beta; + + hCoeffs.coeffRef(i) = h; + } + } + if (NumTraits::IsComplex) + { + // householder transformation on the remaining single scalar + int i = n-2; + Scalar v0 = matA.col(i).coeff(i+1); + RealScalar beta = ei_abs(v0); + if (ei_real(v0)>=0.) + beta = -beta; + matA.col(i).coeffRef(i+1) = beta; + hCoeffs.coeffRef(i) = (beta - v0) / beta; + } +} + +/** reconstructs and returns the matrix Q */ +template +typename Tridiagonalization::MatrixType +Tridiagonalization::matrixQ(void) const +{ + int n = m_matrix.rows(); + MatrixType matQ = MatrixType::identity(n,n); + for (int i = n-2; i>=0; i--) + { + Scalar tmp = m_matrix.coeff(i+1,i); + m_matrix.const_cast_derived().coeffRef(i+1,i) = 1; + + matQ.corner(BottomRight,n-i-1,n-i-1) -= + ((m_hCoeffs[i] * m_matrix.col(i).end(n-i-1)) * + (m_matrix.col(i).end(n-i-1).adjoint() * matQ.corner(BottomRight,n-i-1,n-i-1)).lazy()).lazy(); + + m_matrix.const_cast_derived().coeffRef(i+1,i) = tmp; + } + return matQ; +} + +/** \returns an expression of the diagonal vector */ +template +const typename Tridiagonalization::DiagonalType +Tridiagonalization::diagonal(void) const +{ + return m_matrix.diagonal().nestByValue().real(); +} + +/** \returns an expression of the sub-diagonal vector */ +template +const typename Tridiagonalization::SubDiagonalType +Tridiagonalization::subDiagonal(void) const +{ + int n = m_matrix.rows(); + return Block(m_matrix, 1, 0, n-1,n-1) + .nestByValue().diagonal().nestByValue().real(); +} + +#endif // EIGEN_TRIDIAGONALIZATION_H