Eigne Matrix Class Library
Dependents: Eigen_test Odometry_test AttitudeEstimation_usingTicker MPU9250_Quaternion_Binary_Serial ... more
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
Generated on Tue Jul 12 2022 17:47:02 by 1.7.2