Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifndef EIGEN_BIDIAGONALIZATION_H
00026 #define EIGEN_BIDIAGONALIZATION_H
00027
00028 namespace Eigen {
00029
00030 namespace internal {
00031
00032
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()
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; ; ++k)
00116 {
00117 Index remainingRows = rows - k;
00118 Index remainingCols = cols - k - 1;
00119
00120
00121 m_householder.col(k).tail(remainingRows)
00122 .makeHouseholderInPlace(m_householder.coeffRef(k,k),
00123 m_bidiagonal.template diagonal<0>().coeffRef(k));
00124
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
00133 m_householder.row(k).tail(remainingCols)
00134 .makeHouseholderInPlace(m_householder.coeffRef(k,k+1),
00135 m_bidiagonal.template diagonal<1>().coeffRef(k));
00136
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 }
00160
00161 }
00162
00163 #endif // EIGEN_BIDIAGONALIZATION_H