UpperBidiagonalization.h
Go to the documentation of this file.
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 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_BIDIAGONALIZATION_H
00026 #define EIGEN_BIDIAGONALIZATION_H
00027 
00028 namespace Eigen { 
00029 
00030 namespace internal {
00031 // UpperBidiagonalization will probably be replaced by a Bidiagonalization class, don't want to make it stable API.
00032 // At the same time, it's useful to keep for now as it's about the only thing that is testing the BandMatrix class.
00033 
00034 template<typename _MatrixType> class UpperBidiagonalization
00035 {
00036   public:
00037 
00038     typedef _MatrixType MatrixType;
00039     enum {
00040       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00041       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00042       ColsAtCompileTimeMinusOne = internal::decrement_size<ColsAtCompileTime>::ret
00043     };
00044     typedef typename MatrixType::Scalar Scalar;
00045     typedef typename MatrixType::RealScalar RealScalar;
00046     typedef typename MatrixType::Index Index;
00047     typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
00048     typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType;
00049     typedef BandMatrix<RealScalar, ColsAtCompileTime, ColsAtCompileTime, 1, 0> BidiagonalType;
00050     typedef Matrix<Scalar, ColsAtCompileTime, 1> DiagVectorType;
00051     typedef Matrix<Scalar, ColsAtCompileTimeMinusOne, 1> SuperDiagVectorType;
00052     typedef HouseholderSequence<
00053               const MatrixType,
00054               CwiseUnaryOp<internal::scalar_conjugate_op<Scalar>, const Diagonal<const MatrixType,0> >
00055             > HouseholderUSequenceType;
00056     typedef HouseholderSequence<
00057               const MatrixType,
00058               Diagonal<const MatrixType,1>,
00059               OnTheRight
00060             > HouseholderVSequenceType;
00061     
00068     UpperBidiagonalization() : m_householder(), m_bidiagonal(), m_isInitialized(false) {}
00069 
00070     UpperBidiagonalization(const MatrixType& matrix)
00071       : m_householder(matrix.rows(), matrix.cols()),
00072         m_bidiagonal(matrix.cols(), matrix.cols()),
00073         m_isInitialized(false)
00074     {
00075       compute(matrix);
00076     }
00077     
00078     UpperBidiagonalization& compute(const MatrixType& matrix);
00079     
00080     const MatrixType& householder() const { return m_householder; }
00081     const BidiagonalType& bidiagonal() const { return m_bidiagonal; }
00082     
00083     const HouseholderUSequenceType householderU() const
00084     {
00085       eigen_assert(m_isInitialized && "UpperBidiagonalization is not initialized.");
00086       return HouseholderUSequenceType(m_householder, m_householder.diagonal().conjugate());
00087     }
00088 
00089     const HouseholderVSequenceType householderV() // const here gives nasty errors and i'm lazy
00090     {
00091       eigen_assert(m_isInitialized && "UpperBidiagonalization is not initialized.");
00092       return HouseholderVSequenceType(m_householder, m_householder.const_derived().template diagonal<1>())
00093              .setLength(m_householder.cols()-1)
00094              .setShift(1);
00095     }
00096     
00097   protected:
00098     MatrixType m_householder;
00099     BidiagonalType m_bidiagonal;
00100     bool m_isInitialized;
00101 };
00102 
00103 template<typename _MatrixType>
00104 UpperBidiagonalization<_MatrixType>& UpperBidiagonalization<_MatrixType>::compute(const _MatrixType& matrix)
00105 {
00106   Index rows = matrix.rows();
00107   Index cols = matrix.cols();
00108   
00109   eigen_assert(rows >= cols && "UpperBidiagonalization is only for matrices satisfying rows>=cols.");
00110   
00111   m_householder = matrix;
00112 
00113   ColVectorType temp(rows);
00114 
00115   for (Index k = 0; /* breaks at k==cols-1 below */ ; ++k)
00116   {
00117     Index remainingRows = rows - k;
00118     Index remainingCols = cols - k - 1;
00119 
00120     // construct left householder transform in-place in m_householder
00121     m_householder.col(k).tail(remainingRows)
00122                  .makeHouseholderInPlace(m_householder.coeffRef(k,k),
00123                                          m_bidiagonal.template diagonal<0>().coeffRef(k));
00124     // apply householder transform to remaining part of m_householder on the left
00125     m_householder.bottomRightCorner(remainingRows, remainingCols)
00126                  .applyHouseholderOnTheLeft(m_householder.col(k).tail(remainingRows-1),
00127                                             m_householder.coeff(k,k),
00128                                             temp.data());
00129 
00130     if(k == cols-1) break;
00131     
00132     // construct right householder transform in-place in m_householder
00133     m_householder.row(k).tail(remainingCols)
00134                  .makeHouseholderInPlace(m_householder.coeffRef(k,k+1),
00135                                          m_bidiagonal.template diagonal<1>().coeffRef(k));
00136     // apply householder transform to remaining part of m_householder on the left
00137     m_householder.bottomRightCorner(remainingRows-1, remainingCols)
00138                  .applyHouseholderOnTheRight(m_householder.row(k).tail(remainingCols-1).transpose(),
00139                                              m_householder.coeff(k,k+1),
00140                                              temp.data());
00141   }
00142   m_isInitialized = true;
00143   return *this;
00144 }
00145 
00146 #if 0
00147 
00151 template<typename Derived>
00152 const UpperBidiagonalization<typename MatrixBase<Derived>::PlainObject>
00153 MatrixBase<Derived>::bidiagonalization() const
00154 {
00155   return UpperBidiagonalization<PlainObject>(eval());
00156 }
00157 #endif
00158 
00159 } // end namespace internal
00160 
00161 } // end namespace Eigen
00162 
00163 #endif // EIGEN_BIDIAGONALIZATION_H