diff options
Diffstat (limited to 'eigen/Eigen/src/SVD/UpperBidiagonalization.h')
| -rw-r--r-- | eigen/Eigen/src/SVD/UpperBidiagonalization.h | 148 |
1 files changed, 148 insertions, 0 deletions
diff --git a/eigen/Eigen/src/SVD/UpperBidiagonalization.h b/eigen/Eigen/src/SVD/UpperBidiagonalization.h new file mode 100644 index 0000000..587de37 --- /dev/null +++ b/eigen/Eigen/src/SVD/UpperBidiagonalization.h @@ -0,0 +1,148 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2010 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_BIDIAGONALIZATION_H +#define EIGEN_BIDIAGONALIZATION_H + +namespace Eigen { + +namespace internal { +// UpperBidiagonalization will probably be replaced by a Bidiagonalization class, don't want to make it stable API. +// At the same time, it's useful to keep for now as it's about the only thing that is testing the BandMatrix class. + +template<typename _MatrixType> class UpperBidiagonalization +{ + public: + + typedef _MatrixType MatrixType; + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + ColsAtCompileTimeMinusOne = internal::decrement_size<ColsAtCompileTime>::ret + }; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef typename MatrixType::Index Index; + typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType; + typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType; + typedef BandMatrix<RealScalar, ColsAtCompileTime, ColsAtCompileTime, 1, 0> BidiagonalType; + typedef Matrix<Scalar, ColsAtCompileTime, 1> DiagVectorType; + typedef Matrix<Scalar, ColsAtCompileTimeMinusOne, 1> SuperDiagVectorType; + typedef HouseholderSequence< + const MatrixType, + CwiseUnaryOp<internal::scalar_conjugate_op<Scalar>, const Diagonal<const MatrixType,0> > + > HouseholderUSequenceType; + typedef HouseholderSequence< + const typename internal::remove_all<typename MatrixType::ConjugateReturnType>::type, + Diagonal<const MatrixType,1>, + OnTheRight + > HouseholderVSequenceType; + + /** + * \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via Bidiagonalization::compute(const MatrixType&). + */ + UpperBidiagonalization() : m_householder(), m_bidiagonal(), m_isInitialized(false) {} + + UpperBidiagonalization(const MatrixType& matrix) + : m_householder(matrix.rows(), matrix.cols()), + m_bidiagonal(matrix.cols(), matrix.cols()), + m_isInitialized(false) + { + compute(matrix); + } + + UpperBidiagonalization& compute(const MatrixType& matrix); + + const MatrixType& householder() const { return m_householder; } + const BidiagonalType& bidiagonal() const { return m_bidiagonal; } + + const HouseholderUSequenceType householderU() const + { + eigen_assert(m_isInitialized && "UpperBidiagonalization is not initialized."); + return HouseholderUSequenceType(m_householder, m_householder.diagonal().conjugate()); + } + + const HouseholderVSequenceType householderV() // const here gives nasty errors and i'm lazy + { + eigen_assert(m_isInitialized && "UpperBidiagonalization is not initialized."); + return HouseholderVSequenceType(m_householder.conjugate(), m_householder.const_derived().template diagonal<1>()) + .setLength(m_householder.cols()-1) + .setShift(1); + } + + protected: + MatrixType m_householder; + BidiagonalType m_bidiagonal; + bool m_isInitialized; +}; + +template<typename _MatrixType> +UpperBidiagonalization<_MatrixType>& UpperBidiagonalization<_MatrixType>::compute(const _MatrixType& matrix) +{ + Index rows = matrix.rows(); + Index cols = matrix.cols(); + + eigen_assert(rows >= cols && "UpperBidiagonalization is only for matrices satisfying rows>=cols."); + + m_householder = matrix; + + ColVectorType temp(rows); + + for (Index k = 0; /* breaks at k==cols-1 below */ ; ++k) + { + Index remainingRows = rows - k; + Index remainingCols = cols - k - 1; + + // construct left householder transform in-place in m_householder + m_householder.col(k).tail(remainingRows) + .makeHouseholderInPlace(m_householder.coeffRef(k,k), + m_bidiagonal.template diagonal<0>().coeffRef(k)); + // apply householder transform to remaining part of m_householder on the left + m_householder.bottomRightCorner(remainingRows, remainingCols) + .applyHouseholderOnTheLeft(m_householder.col(k).tail(remainingRows-1), + m_householder.coeff(k,k), + temp.data()); + + if(k == cols-1) break; + + // construct right householder transform in-place in m_householder + m_householder.row(k).tail(remainingCols) + .makeHouseholderInPlace(m_householder.coeffRef(k,k+1), + m_bidiagonal.template diagonal<1>().coeffRef(k)); + // apply householder transform to remaining part of m_householder on the left + m_householder.bottomRightCorner(remainingRows-1, remainingCols) + .applyHouseholderOnTheRight(m_householder.row(k).tail(remainingCols-1).transpose(), + m_householder.coeff(k,k+1), + temp.data()); + } + m_isInitialized = true; + return *this; +} + +#if 0 +/** \return the Householder QR decomposition of \c *this. + * + * \sa class Bidiagonalization + */ +template<typename Derived> +const UpperBidiagonalization<typename MatrixBase<Derived>::PlainObject> +MatrixBase<Derived>::bidiagonalization() const +{ + return UpperBidiagonalization<PlainObject>(eval()); +} +#endif + +} // end namespace internal + +} // end namespace Eigen + +#endif // EIGEN_BIDIAGONALIZATION_H |