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 internal {
00029 // UpperBidiagonalization will probably be replaced by a Bidiagonalization class, don't want to make it stable API.
00030 // At the same time, it's useful to keep for now as it's about the only thing that is testing the BandMatrix class.
00031 
00032 template<typename _MatrixType> class UpperBidiagonalization
00033 {
00034   public:
00035 
00036     typedef _MatrixType MatrixType;
00037     enum {
00038       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00039       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00040       ColsAtCompileTimeMinusOne = internal::decrement_size<ColsAtCompileTime>::ret
00041     };
00042     typedef typename MatrixType::Scalar Scalar;
00043     typedef typename MatrixType::RealScalar RealScalar;
00044     typedef typename MatrixType::Index Index;
00045     typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
00046     typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType;
00047     typedef BandMatrix<RealScalar, ColsAtCompileTime, ColsAtCompileTime, 1, 0> BidiagonalType;
00048     typedef Matrix<Scalar, ColsAtCompileTime, 1> DiagVectorType;
00049     typedef Matrix<Scalar, ColsAtCompileTimeMinusOne, 1> SuperDiagVectorType;
00050     typedef HouseholderSequence<
00051               const MatrixType,
00052               CwiseUnaryOp<internal::scalar_conjugate_op<Scalar>, const Diagonal<const MatrixType,0> >
00053             > HouseholderUSequenceType;
00054     typedef HouseholderSequence<
00055               const MatrixType,
00056               Diagonal<const MatrixType,1>,
00057               OnTheRight
00058             > HouseholderVSequenceType;
00059     
00066     UpperBidiagonalization() : m_householder(), m_bidiagonal(), m_isInitialized(false) {}
00067 
00068     UpperBidiagonalization(const MatrixType& matrix)
00069       : m_householder(matrix.rows(), matrix.cols()),
00070         m_bidiagonal(matrix.cols(), matrix.cols()),
00071         m_isInitialized(false)
00072     {
00073       compute(matrix);
00074     }
00075     
00076     UpperBidiagonalization& compute(const MatrixType& matrix);
00077     
00078     const MatrixType& householder() const { return m_householder; }
00079     const BidiagonalType& bidiagonal() const { return m_bidiagonal; }
00080     
00081     const HouseholderUSequenceType householderU() const
00082     {
00083       eigen_assert(m_isInitialized && "UpperBidiagonalization is not initialized.");
00084       return HouseholderUSequenceType(m_householder, m_householder.diagonal().conjugate());
00085     }
00086 
00087     const HouseholderVSequenceType householderV() // const here gives nasty errors and i'm lazy
00088     {
00089       eigen_assert(m_isInitialized && "UpperBidiagonalization is not initialized.");
00090       return HouseholderVSequenceType(m_householder, m_householder.const_derived().template diagonal<1>())
00091              .setLength(m_householder.cols()-1)
00092              .setShift(1);
00093     }
00094     
00095   protected:
00096     MatrixType m_householder;
00097     BidiagonalType m_bidiagonal;
00098     bool m_isInitialized;
00099 };
00100 
00101 template<typename _MatrixType>
00102 UpperBidiagonalization<_MatrixType>& UpperBidiagonalization<_MatrixType>::compute(const _MatrixType& matrix)
00103 {
00104   Index rows = matrix.rows();
00105   Index cols = matrix.cols();
00106   
00107   eigen_assert(rows >= cols && "UpperBidiagonalization is only for matrices satisfying rows>=cols.");
00108   
00109   m_householder = matrix;
00110 
00111   ColVectorType temp(rows);
00112 
00113   for (Index k = 0; /* breaks at k==cols-1 below */ ; ++k)
00114   {
00115     Index remainingRows = rows - k;
00116     Index remainingCols = cols - k - 1;
00117 
00118     // construct left householder transform in-place in m_householder
00119     m_householder.col(k).tail(remainingRows)
00120                  .makeHouseholderInPlace(m_householder.coeffRef(k,k),
00121                                          m_bidiagonal.template diagonal<0>().coeffRef(k));
00122     // apply householder transform to remaining part of m_householder on the left
00123     m_householder.bottomRightCorner(remainingRows, remainingCols)
00124                  .applyHouseholderOnTheLeft(m_householder.col(k).tail(remainingRows-1),
00125                                             m_householder.coeff(k,k),
00126                                             temp.data());
00127 
00128     if(k == cols-1) break;
00129     
00130     // construct right householder transform in-place in m_householder
00131     m_householder.row(k).tail(remainingCols)
00132                  .makeHouseholderInPlace(m_householder.coeffRef(k,k+1),
00133                                          m_bidiagonal.template diagonal<1>().coeffRef(k));
00134     // apply householder transform to remaining part of m_householder on the left
00135     m_householder.bottomRightCorner(remainingRows-1, remainingCols)
00136                  .applyHouseholderOnTheRight(m_householder.row(k).tail(remainingCols-1).transpose(),
00137                                              m_householder.coeff(k,k+1),
00138                                              temp.data());
00139   }
00140   m_isInitialized = true;
00141   return *this;
00142 }
00143 
00144 #if 0
00145 
00149 template<typename Derived>
00150 const UpperBidiagonalization<typename MatrixBase<Derived>::PlainObject>
00151 MatrixBase<Derived>::bidiagonalization() const
00152 {
00153   return UpperBidiagonalization<PlainObject>(eval());
00154 }
00155 #endif
00156 
00157 } // end namespace internal
00158 
00159 #endif // EIGEN_BIDIAGONALIZATION_H


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:34:02