25 #ifndef EIGEN_BIDIAGONALIZATION_H
26 #define EIGEN_BIDIAGONALIZATION_H
34 template<
typename _MatrixType>
class UpperBidiagonalization
38 typedef _MatrixType MatrixType;
40 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
41 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
42 ColsAtCompileTimeMinusOne = internal::decrement_size<ColsAtCompileTime>::ret
44 typedef typename MatrixType::Scalar Scalar;
45 typedef typename MatrixType::RealScalar RealScalar;
46 typedef typename MatrixType::Index Index;
47 typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
48 typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType;
49 typedef BandMatrix<RealScalar, ColsAtCompileTime, ColsAtCompileTime, 1, 0> BidiagonalType;
50 typedef Matrix<Scalar, ColsAtCompileTime, 1> DiagVectorType;
51 typedef Matrix<Scalar, ColsAtCompileTimeMinusOne, 1> SuperDiagVectorType;
52 typedef HouseholderSequence<
54 CwiseUnaryOp<internal::scalar_conjugate_op<Scalar>,
const Diagonal<const MatrixType,0> >
55 > HouseholderUSequenceType;
56 typedef HouseholderSequence<
58 Diagonal<const MatrixType,1>,
60 > HouseholderVSequenceType;
68 UpperBidiagonalization() : m_householder(), m_bidiagonal(), m_isInitialized(false) {}
70 UpperBidiagonalization(
const MatrixType& matrix)
71 : m_householder(matrix.rows(), matrix.cols()),
72 m_bidiagonal(matrix.cols(), matrix.cols()),
73 m_isInitialized(false)
78 UpperBidiagonalization& compute(
const MatrixType& matrix);
80 const MatrixType& householder()
const {
return m_householder; }
81 const BidiagonalType& bidiagonal()
const {
return m_bidiagonal; }
83 const HouseholderUSequenceType householderU()
const
85 eigen_assert(m_isInitialized &&
"UpperBidiagonalization is not initialized.");
86 return HouseholderUSequenceType(m_householder, m_householder.diagonal().conjugate());
89 const HouseholderVSequenceType householderV()
91 eigen_assert(m_isInitialized &&
"UpperBidiagonalization is not initialized.");
92 return HouseholderVSequenceType(m_householder, m_householder.const_derived().template diagonal<1>())
93 .setLength(m_householder.cols()-1)
98 MatrixType m_householder;
99 BidiagonalType m_bidiagonal;
100 bool m_isInitialized;
103 template<
typename _MatrixType>
104 UpperBidiagonalization<_MatrixType>& UpperBidiagonalization<_MatrixType>::compute(
const _MatrixType& matrix)
106 Index rows = matrix.rows();
107 Index cols = matrix.cols();
109 eigen_assert(rows >= cols &&
"UpperBidiagonalization is only for matrices satisfying rows>=cols.");
111 m_householder = matrix;
113 ColVectorType temp(rows);
115 for (Index k = 0; ; ++k)
117 Index remainingRows = rows - k;
118 Index remainingCols = cols - k - 1;
121 m_householder.col(k).tail(remainingRows)
122 .makeHouseholderInPlace(m_householder.coeffRef(k,k),
123 m_bidiagonal.template diagonal<0>().coeffRef(k));
125 m_householder.bottomRightCorner(remainingRows, remainingCols)
126 .applyHouseholderOnTheLeft(m_householder.col(k).tail(remainingRows-1),
127 m_householder.coeff(k,k),
130 if(k == cols-1)
break;
133 m_householder.row(k).tail(remainingCols)
134 .makeHouseholderInPlace(m_householder.coeffRef(k,k+1),
135 m_bidiagonal.template diagonal<1>().coeffRef(k));
137 m_householder.bottomRightCorner(remainingRows-1, remainingCols)
138 .applyHouseholderOnTheRight(m_householder.row(k).tail(remainingCols-1).transpose(),
139 m_householder.coeff(k,k+1),
142 m_isInitialized =
true;
151 template<
typename Derived>
152 const UpperBidiagonalization<typename MatrixBase<Derived>::PlainObject>
153 MatrixBase<Derived>::bidiagonalization()
const
155 return UpperBidiagonalization<PlainObject>(eval());
163 #endif // EIGEN_BIDIAGONALIZATION_H