From 4a04f7206914a49f5f95adc5eb786237f1a9f547 Mon Sep 17 00:00:00 2001 From: Campbell Barton Date: Sun, 23 Oct 2011 17:52:20 +0000 Subject: remove $Id: tags after discussion on the mailign list: http://markmail.org/message/fp7ozcywxum3ar7n --- .../Eigen3/Eigen/src/SVD/UpperBidiagonalization.h | 159 +++++++++++++++++++++ 1 file changed, 159 insertions(+) create mode 100644 extern/Eigen3/Eigen/src/SVD/UpperBidiagonalization.h (limited to 'extern/Eigen3/Eigen/src/SVD/UpperBidiagonalization.h') diff --git a/extern/Eigen3/Eigen/src/SVD/UpperBidiagonalization.h b/extern/Eigen3/Eigen/src/SVD/UpperBidiagonalization.h new file mode 100644 index 00000000000..2de197da953 --- /dev/null +++ b/extern/Eigen3/Eigen/src/SVD/UpperBidiagonalization.h @@ -0,0 +1,159 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2010 Benoit Jacob +// +// 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_BIDIAGONALIZATION_H +#define EIGEN_BIDIAGONALIZATION_H + +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 class UpperBidiagonalization +{ + public: + + typedef _MatrixType MatrixType; + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + ColsAtCompileTimeMinusOne = internal::decrement_size::ret + }; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef typename MatrixType::Index Index; + typedef Matrix RowVectorType; + typedef Matrix ColVectorType; + typedef BandMatrix BidiagonalType; + typedef Matrix DiagVectorType; + typedef Matrix SuperDiagVectorType; + typedef HouseholderSequence< + const MatrixType, + CwiseUnaryOp, const Diagonal > + > HouseholderUSequenceType; + typedef HouseholderSequence< + const MatrixType, + Diagonal, + 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, 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 +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 +const UpperBidiagonalization::PlainObject> +MatrixBase::bidiagonalization() const +{ + return UpperBidiagonalization(eval()); +} +#endif + +} // end namespace internal + +#endif // EIGEN_BIDIAGONALIZATION_H -- cgit v1.2.3