UpperBidiagonalization.h

00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2010 Benoit Jacob <jacob.benoit.1@gmail.com>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_BIDIAGONALIZATION_H
00011 #define EIGEN_BIDIAGONALIZATION_H
00012 
00013 namespace Eigen { 
00014 
00015 namespace internal {
00016 // UpperBidiagonalization will probably be replaced by a Bidiagonalization class, don't want to make it stable API.
00017 // At the same time, it's useful to keep for now as it's about the only thing that is testing the BandMatrix class.
00018 
00019 template<typename _MatrixType> class UpperBidiagonalization
00020 {
00021   public:
00022 
00023     typedef _MatrixType MatrixType;
00024     enum {
00025       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00026       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00027       ColsAtCompileTimeMinusOne = internal::decrement_size<ColsAtCompileTime>::ret
00028     };
00029     typedef typename MatrixType::Scalar Scalar;
00030     typedef typename MatrixType::RealScalar RealScalar;
00031     typedef typename MatrixType::Index Index;
00032     typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
00033     typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType;
00034     typedef BandMatrix<RealScalar, ColsAtCompileTime, ColsAtCompileTime, 1, 0> BidiagonalType;
00035     typedef Matrix<Scalar, ColsAtCompileTime, 1> DiagVectorType;
00036     typedef Matrix<Scalar, ColsAtCompileTimeMinusOne, 1> SuperDiagVectorType;
00037     typedef HouseholderSequence<
00038               const MatrixType,
00039               CwiseUnaryOp<internal::scalar_conjugate_op<Scalar>, const Diagonal<const MatrixType,0> >
00040             > HouseholderUSequenceType;
00041     typedef HouseholderSequence<
00042               const typename internal::remove_all<typename MatrixType::ConjugateReturnType>::type,
00043               Diagonal<const MatrixType,1>,
00044               OnTheRight
00045             > HouseholderVSequenceType;
00046     
00047     /**
00048     * \brief Default Constructor.
00049     *
00050     * The default constructor is useful in cases in which the user intends to
00051     * perform decompositions via Bidiagonalization::compute(const MatrixType&).
00052     */
00053     UpperBidiagonalization() : m_householder(), m_bidiagonal(), m_isInitialized(false) {}
00054 
00055     UpperBidiagonalization(const MatrixType& matrix)
00056       : m_householder(matrix.rows(), matrix.cols()),
00057         m_bidiagonal(matrix.cols(), matrix.cols()),
00058         m_isInitialized(false)
00059     {
00060       compute(matrix);
00061     }
00062     
00063     UpperBidiagonalization& compute(const MatrixType& matrix);
00064     
00065     const MatrixType& householder() const { return m_householder; }
00066     const BidiagonalType& bidiagonal() const { return m_bidiagonal; }
00067     
00068     const HouseholderUSequenceType householderU() const
00069     {
00070       eigen_assert(m_isInitialized && "UpperBidiagonalization is not initialized.");
00071       return HouseholderUSequenceType(m_householder, m_householder.diagonal().conjugate());
00072     }
00073 
00074     const HouseholderVSequenceType householderV() // const here gives nasty errors and i'm lazy
00075     {
00076       eigen_assert(m_isInitialized && "UpperBidiagonalization is not initialized.");
00077       return HouseholderVSequenceType(m_householder.conjugate(), m_householder.const_derived().template diagonal<1>())
00078              .setLength(m_householder.cols()-1)
00079              .setShift(1);
00080     }
00081     
00082   protected:
00083     MatrixType m_householder;
00084     BidiagonalType m_bidiagonal;
00085     bool m_isInitialized;
00086 };
00087 
00088 template<typename _MatrixType>
00089 UpperBidiagonalization<_MatrixType>& UpperBidiagonalization<_MatrixType>::compute(const _MatrixType& matrix)
00090 {
00091   Index rows = matrix.rows();
00092   Index cols = matrix.cols();
00093   
00094   eigen_assert(rows >= cols && "UpperBidiagonalization is only for matrices satisfying rows>=cols.");
00095   
00096   m_householder = matrix;
00097 
00098   ColVectorType temp(rows);
00099 
00100   for (Index k = 0; /* breaks at k==cols-1 below */ ; ++k)
00101   {
00102     Index remainingRows = rows - k;
00103     Index remainingCols = cols - k - 1;
00104 
00105     // construct left householder transform in-place in m_householder
00106     m_householder.col(k).tail(remainingRows)
00107                  .makeHouseholderInPlace(m_householder.coeffRef(k,k),
00108                                          m_bidiagonal.template diagonal<0>().coeffRef(k));
00109     // apply householder transform to remaining part of m_householder on the left
00110     m_householder.bottomRightCorner(remainingRows, remainingCols)
00111                  .applyHouseholderOnTheLeft(m_householder.col(k).tail(remainingRows-1),
00112                                             m_householder.coeff(k,k),
00113                                             temp.data());
00114 
00115     if(k == cols-1) break;
00116     
00117     // construct right householder transform in-place in m_householder
00118     m_householder.row(k).tail(remainingCols)
00119                  .makeHouseholderInPlace(m_householder.coeffRef(k,k+1),
00120                                          m_bidiagonal.template diagonal<1>().coeffRef(k));
00121     // apply householder transform to remaining part of m_householder on the left
00122     m_householder.bottomRightCorner(remainingRows-1, remainingCols)
00123                  .applyHouseholderOnTheRight(m_householder.row(k).tail(remainingCols-1).transpose(),
00124                                              m_householder.coeff(k,k+1),
00125                                              temp.data());
00126   }
00127   m_isInitialized = true;
00128   return *this;
00129 }
00130 
00131 #if 0
00132 /** \return the Householder QR decomposition of \c *this.
00133   *
00134   * \sa class Bidiagonalization
00135   */
00136 template<typename Derived>
00137 const UpperBidiagonalization<typename MatrixBase<Derived>::PlainObject>
00138 MatrixBase<Derived>::bidiagonalization() const
00139 {
00140   return UpperBidiagonalization<PlainObject>(eval());
00141 }
00142 #endif
00143 
00144 } // end namespace internal
00145 
00146 } // end namespace Eigen
00147 
00148 #endif // EIGEN_BIDIAGONALIZATION_H