11 #ifndef EIGEN_BIDIAGONALIZATION_H 12 #define EIGEN_BIDIAGONALIZATION_H 20 template<
typename _MatrixType>
class UpperBidiagonalization
24 typedef _MatrixType MatrixType;
26 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
27 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
28 ColsAtCompileTimeMinusOne = internal::decrement_size<ColsAtCompileTime>::ret
30 typedef typename MatrixType::Scalar Scalar;
31 typedef typename MatrixType::RealScalar RealScalar;
32 typedef Eigen::Index Index;
33 typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
34 typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType;
35 typedef BandMatrix<RealScalar, ColsAtCompileTime, ColsAtCompileTime, 1, 0, RowMajor> BidiagonalType;
36 typedef Matrix<Scalar, ColsAtCompileTime, 1> DiagVectorType;
37 typedef Matrix<Scalar, ColsAtCompileTimeMinusOne, 1> SuperDiagVectorType;
38 typedef HouseholderSequence<
40 const typename internal::remove_all<typename Diagonal<const MatrixType,0>::ConjugateReturnType>::type
41 > HouseholderUSequenceType;
42 typedef HouseholderSequence<
43 const typename internal::remove_all<typename MatrixType::ConjugateReturnType>::type,
44 Diagonal<const MatrixType,1>,
46 > HouseholderVSequenceType;
54 UpperBidiagonalization() : m_householder(), m_bidiagonal(), m_isInitialized(false) {}
56 explicit UpperBidiagonalization(
const MatrixType& matrix)
57 : m_householder(matrix.rows(), matrix.cols()),
58 m_bidiagonal(matrix.cols(), matrix.cols()),
59 m_isInitialized(false)
64 UpperBidiagonalization& compute(
const MatrixType& matrix);
65 UpperBidiagonalization& computeUnblocked(
const MatrixType& matrix);
67 const MatrixType& householder()
const {
return m_householder; }
68 const BidiagonalType& bidiagonal()
const {
return m_bidiagonal; }
70 const HouseholderUSequenceType householderU()
const 72 eigen_assert(m_isInitialized &&
"UpperBidiagonalization is not initialized.");
73 return HouseholderUSequenceType(m_householder, m_householder.diagonal().conjugate());
76 const HouseholderVSequenceType householderV()
78 eigen_assert(m_isInitialized &&
"UpperBidiagonalization is not initialized.");
79 return HouseholderVSequenceType(m_householder.conjugate(), m_householder.const_derived().template diagonal<1>())
80 .setLength(m_householder.cols()-1)
85 MatrixType m_householder;
86 BidiagonalType m_bidiagonal;
92 template<
typename MatrixType>
93 void upperbidiagonalization_inplace_unblocked(MatrixType& mat,
94 typename MatrixType::RealScalar *diagonal,
95 typename MatrixType::RealScalar *upper_diagonal,
96 typename MatrixType::Scalar* tempData = 0)
98 typedef typename MatrixType::Scalar Scalar;
100 Index rows = mat.rows();
101 Index cols = mat.cols();
103 typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixType::MaxRowsAtCompileTime,1> TempType;
107 tempVector.resize(rows);
108 tempData = tempVector.data();
111 for (Index k = 0; ; ++k)
113 Index remainingRows = rows - k;
114 Index remainingCols = cols - k - 1;
117 mat.col(k).tail(remainingRows)
118 .makeHouseholderInPlace(mat.coeffRef(k,k), diagonal[k]);
120 mat.bottomRightCorner(remainingRows, remainingCols)
121 .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), mat.coeff(k,k), tempData);
123 if(k == cols-1)
break;
126 mat.row(k).tail(remainingCols)
127 .makeHouseholderInPlace(mat.coeffRef(k,k+1), upper_diagonal[k]);
129 mat.bottomRightCorner(remainingRows-1, remainingCols)
130 .applyHouseholderOnTheRight(mat.row(k).tail(remainingCols-1).transpose(), mat.coeff(k,k+1), tempData);
151 template<
typename MatrixType>
152 void upperbidiagonalization_blocked_helper(MatrixType& A,
153 typename MatrixType::RealScalar *diagonal,
154 typename MatrixType::RealScalar *upper_diagonal,
156 Ref<Matrix<
typename MatrixType::Scalar, Dynamic, Dynamic,
158 Ref<Matrix<
typename MatrixType::Scalar, Dynamic, Dynamic,
161 typedef typename MatrixType::Scalar Scalar;
162 enum { StorageOrder = traits<MatrixType>::Flags &
RowMajorBit };
163 typedef InnerStride<int(StorageOrder) == int(ColMajor) ? 1 : Dynamic> ColInnerStride;
164 typedef InnerStride<int(StorageOrder) == int(ColMajor) ? Dynamic : 1> RowInnerStride;
165 typedef Ref<Matrix<Scalar, Dynamic, 1>, 0, ColInnerStride> SubColumnType;
166 typedef Ref<Matrix<Scalar, 1, Dynamic>, 0, RowInnerStride> SubRowType;
167 typedef Ref<Matrix<Scalar, Dynamic, Dynamic, StorageOrder > > SubMatType;
169 Index brows = A.rows();
170 Index bcols = A.cols();
172 Scalar tau_u, tau_u_prev(0), tau_v;
174 for(Index k = 0; k < bs; ++k)
176 Index remainingRows = brows - k;
177 Index remainingCols = bcols - k - 1;
179 SubMatType X_k1( X.block(k,0, remainingRows,k) );
180 SubMatType V_k1( A.block(k,0, remainingRows,k) );
183 SubColumnType v_k = A.col(k).tail(remainingRows);
184 v_k -= V_k1 * Y.row(k).head(k).adjoint();
185 if(k) v_k -= X_k1 * A.col(k).head(k);
188 v_k.makeHouseholderInPlace(tau_v, diagonal[k]);
192 SubMatType Y_k ( Y.block(k+1,0, remainingCols, k+1) );
193 SubMatType U_k1 ( A.block(0,k+1, k,remainingCols) );
201 SubColumnType y_k( Y.col(k).tail(remainingCols) );
204 SubColumnType tmp( Y.col(k).head(k) );
205 y_k.noalias() = A.block(k,k+1, remainingRows,remainingCols).adjoint() * v_k;
206 tmp.noalias() = V_k1.adjoint() * v_k;
207 y_k.noalias() -= Y_k.leftCols(k) * tmp;
208 tmp.noalias() = X_k1.adjoint() * v_k;
209 y_k.noalias() -= U_k1.adjoint() * tmp;
210 y_k *= numext::conj(tau_v);
214 SubRowType u_k( A.row(k).tail(remainingCols) );
215 u_k = u_k.conjugate();
217 u_k -= Y_k * A.row(k).head(k+1).adjoint();
218 if(k) u_k -= U_k1.adjoint() * X.row(k).head(k).adjoint();
222 u_k.makeHouseholderInPlace(tau_u, upper_diagonal[k]);
226 A(k,k+1) = Scalar(1);
230 SubColumnType x_k ( X.col(k).tail(remainingRows-1) );
234 SubColumnType tmp0 ( X.col(k).head(k) ),
235 tmp1 ( X.col(k).head(k+1) );
237 x_k.noalias() = A.block(k+1,k+1, remainingRows-1,remainingCols) * u_k.transpose();
238 tmp0.noalias() = U_k1 * u_k.transpose();
239 x_k.noalias() -= X_k1.bottomRows(remainingRows-1) * tmp0;
240 tmp1.noalias() = Y_k.adjoint() * u_k.transpose();
241 x_k.noalias() -= A.block(k+1,0, remainingRows-1,k+1) * tmp1;
242 x_k *= numext::conj(tau_u);
243 tau_u = numext::conj(tau_u);
244 u_k = u_k.conjugate();
247 if(k>0) A.coeffRef(k-1,k) = tau_u_prev;
251 A.coeffRef(k-1,k) = tau_u_prev;
253 A.coeffRef(k,k) = tau_v;
257 A.coeffRef(bs-1,bs) = tau_u_prev;
260 if(bcols>bs && brows>bs)
262 SubMatType A11( A.bottomRightCorner(brows-bs,bcols-bs) );
263 SubMatType A10( A.block(bs,0, brows-bs,bs) );
264 SubMatType A01( A.block(0,bs, bs,bcols-bs) );
265 Scalar tmp = A01(bs-1,0);
267 A11.noalias() -= A10 * Y.topLeftCorner(bcols,bs).bottomRows(bcols-bs).adjoint();
268 A11.noalias() -= X.topLeftCorner(brows,bs).bottomRows(brows-bs) * A01;
281 template<
typename MatrixType,
typename B
idiagType>
282 void upperbidiagonalization_inplace_blocked(MatrixType& A, BidiagType& bidiagonal,
283 Index maxBlockSize=32,
284 typename MatrixType::Scalar* = 0)
286 typedef typename MatrixType::Scalar Scalar;
287 typedef Block<MatrixType,Dynamic,Dynamic> BlockType;
289 Index rows = A.rows();
290 Index cols = A.cols();
291 Index size = (std::min)(rows, cols);
294 enum { StorageOrder = traits<MatrixType>::Flags & RowMajorBit };
296 MatrixType::RowsAtCompileTime,
299 MatrixType::MaxRowsAtCompileTime> X(rows,maxBlockSize);
301 MatrixType::ColsAtCompileTime,
304 MatrixType::MaxColsAtCompileTime> Y(cols,maxBlockSize);
305 Index blockSize = (std::min)(maxBlockSize,size);
308 for(k = 0; k < size; k += blockSize)
310 Index bs = (std::min)(size-k,blockSize);
311 Index brows = rows - k;
312 Index bcols = cols - k;
328 BlockType B = A.block(k,k,brows,bcols);
334 if(k+bs==cols || bcols<48)
336 upperbidiagonalization_inplace_unblocked(B,
337 &(bidiagonal.template diagonal<0>().coeffRef(k)),
338 &(bidiagonal.template diagonal<1>().coeffRef(k)),
345 upperbidiagonalization_blocked_helper<BlockType>( B,
346 &(bidiagonal.template diagonal<0>().coeffRef(k)),
347 &(bidiagonal.template diagonal<1>().coeffRef(k)),
349 X.topLeftCorner(brows,bs),
350 Y.topLeftCorner(bcols,bs)
356 template<
typename _MatrixType>
357 UpperBidiagonalization<_MatrixType>& UpperBidiagonalization<_MatrixType>::computeUnblocked(
const _MatrixType& matrix)
359 Index rows = matrix.rows();
360 Index cols = matrix.cols();
361 EIGEN_ONLY_USED_FOR_DEBUG(cols);
363 eigen_assert(rows >= cols &&
"UpperBidiagonalization is only for Arices satisfying rows>=cols.");
365 m_householder = matrix;
367 ColVectorType temp(rows);
369 upperbidiagonalization_inplace_unblocked(m_householder,
370 &(m_bidiagonal.template diagonal<0>().coeffRef(0)),
371 &(m_bidiagonal.template diagonal<1>().coeffRef(0)),
374 m_isInitialized =
true;
378 template<
typename _MatrixType>
379 UpperBidiagonalization<_MatrixType>& UpperBidiagonalization<_MatrixType>::compute(
const _MatrixType& matrix)
381 Index rows = matrix.rows();
382 Index cols = matrix.cols();
383 EIGEN_ONLY_USED_FOR_DEBUG(rows);
384 EIGEN_ONLY_USED_FOR_DEBUG(cols);
386 eigen_assert(rows >= cols &&
"UpperBidiagonalization is only for Arices satisfying rows>=cols.");
388 m_householder = matrix;
389 upperbidiagonalization_inplace_blocked(m_householder, m_bidiagonal);
391 m_isInitialized =
true;
400 template<
typename Derived>
401 const UpperBidiagonalization<typename MatrixBase<Derived>::PlainObject>
402 MatrixBase<Derived>::bidiagonalization()
const 404 return UpperBidiagonalization<PlainObject>(eval());
412 #endif // EIGEN_BIDIAGONALIZATION_H
const unsigned int RowMajorBit
Definition: Constants.h:61
Definition: Constants.h:335
Definition: Eigen_Colamd.h:54