diff options
| -rw-r--r-- | Eigen/QR | 1 | ||||
| -rw-r--r-- | Eigen/src/Core/CwiseUnaryOp.h | 12 | ||||
| -rwxr-xr-x | Eigen/src/Core/Extract.h | 2 | ||||
| -rw-r--r-- | Eigen/src/Core/Functors.h | 13 | ||||
| -rw-r--r-- | Eigen/src/Core/MatrixBase.h | 3 | ||||
| -rw-r--r-- | Eigen/src/Core/util/Constants.h | 2 | ||||
| -rw-r--r-- | Eigen/src/Core/util/ForwardDeclarations.h | 1 | ||||
| -rwxr-xr-x | Eigen/src/QR/Tridiagonalization.h | 239 |
8 files changed, 269 insertions, 4 deletions
@@ -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<Derived>::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<typename Derived> @@ -149,6 +149,16 @@ MatrixBase<Derived>::conjugate() const return ConjugateReturnType(derived()); } +/** \returns an expression of the real part of \c *this. + * + * \sa adjoint() */ +template<typename Derived> +inline const typename MatrixBase<Derived>::RealReturnType +MatrixBase<Derived>::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<typename MatrixType, unsigned int Mode> 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<typename Scalar, typename NewType> struct ei_functor_traits<ei_scalar_cast_op<Scalar,NewType> > { enum { Cost = ei_is_same_type<Scalar, NewType>::ret ? 0 : NumTraits<NewType>::AddCost, IsVectorizable = false }; }; +/** \internal + * \brief Template functor to extract the real part of a complex + * + * \sa class CwiseUnaryOp, MatrixBase::real() + */ +template<typename Scalar> +struct ei_scalar_real_op EIGEN_EMPTY_STRUCT { + typedef typename NumTraits<Scalar>::Real result_type; + inline result_type operator() (const Scalar& a) const { return ei_real(a); } +}; +template<typename Scalar> +struct ei_functor_traits<ei_scalar_real_op<Scalar> > +{ 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<typename Derived> class MatrixBase : public ArrayBase<Derived> CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, Derived>, Derived& >::ret ConjugateReturnType; + /** the return type of MatrixBase::real() */ + typedef CwiseUnaryOp<ei_scalar_real_op<Scalar>, Derived> RealReturnType; /** the return type of MatrixBase::adjoint() */ typedef Transpose<NestByValue<typename ei_unref<ConjugateReturnType>::type> > AdjointReturnType; @@ -477,6 +479,7 @@ template<typename Derived> class MatrixBase : public ArrayBase<Derived> /// \name Coefficient-wise operations //@{ const ConjugateReturnType conjugate() const; + const RealReturnType real() const; template<typename OtherDerived> const CwiseBinaryOp<ei_scalar_product_op<typename ei_traits<Derived>::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<typename Scalar> struct ei_scalar_product_op; template<typename Scalar> struct ei_scalar_quotient_op; template<typename Scalar> struct ei_scalar_opposite_op; template<typename Scalar> struct ei_scalar_conjugate_op; +template<typename Scalar> struct ei_scalar_real_op; template<typename Scalar> struct ei_scalar_abs_op; template<typename Scalar> struct ei_scalar_abs2_op; template<typename Scalar> 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 <g.gael@free.fr> +// +// 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 <http://www.gnu.org/licenses/>. + +#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<typename _MatrixType> class Tridiagonalization +{ + public: + + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + + enum {SizeMinusOne = MatrixType::RowsAtCompileTime==Dynamic + ? Dynamic + : MatrixType::RowsAtCompileTime-1}; + + typedef Matrix<Scalar, SizeMinusOne, 1> CoeffVectorType; + + typedef typename NestByValue<DiagonalCoeffs<MatrixType> >::RealReturnType DiagonalType; + + typedef typename NestByValue<DiagonalCoeffs< + NestByValue<Block< + MatrixType,SizeMinusOne,SizeMinusOne> > > >::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<typename MatrixType> +void Tridiagonalization<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs) +{ + assert(matA.rows()==matA.cols()); + int n = matA.rows(); + for (int i = 0; i<n-2; ++i) + { + // let's consider the vector v = i-th column starting at position i+1 + + // start of the householder transformation + // squared norm of the vector v skipping the first element + RealScalar v1norm2 = matA.col(i).end(n-(i+2)).norm2(); + + if (ei_isMuchSmallerThan(v1norm2,static_cast<Scalar>(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<Lower|SelfAdjoint>() + * 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<Lower>() = + 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<Scalar>::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 MatrixType> +typename Tridiagonalization<MatrixType>::MatrixType +Tridiagonalization<MatrixType>::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<typename MatrixType> +const typename Tridiagonalization<MatrixType>::DiagonalType +Tridiagonalization<MatrixType>::diagonal(void) const +{ + return m_matrix.diagonal().nestByValue().real(); +} + +/** \returns an expression of the sub-diagonal vector */ +template<typename MatrixType> +const typename Tridiagonalization<MatrixType>::SubDiagonalType +Tridiagonalization<MatrixType>::subDiagonal(void) const +{ + int n = m_matrix.rows(); + return Block<MatrixType,SizeMinusOne,SizeMinusOne>(m_matrix, 1, 0, n-1,n-1) + .nestByValue().diagonal().nestByValue().real(); +} + +#endif // EIGEN_TRIDIAGONALIZATION_H |