1// This file is part of Eigen, a lightweight C++ template library 2// for linear algebra. 3// 4// Copyright (C) 2010 Benoit Jacob <jacob.benoit.1@gmail.com> 5// Copyright (C) 2013-2014 Gael Guennebaud <gael.guennebaud@inria.fr> 6// 7// This Source Code Form is subject to the terms of the Mozilla 8// Public License v. 2.0. If a copy of the MPL was not distributed 9// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11#ifndef EIGEN_BIDIAGONALIZATION_H 12#define EIGEN_BIDIAGONALIZATION_H 13 14namespace Eigen { 15 16namespace internal { 17// UpperBidiagonalization will probably be replaced by a Bidiagonalization class, don't want to make it stable API. 18// At the same time, it's useful to keep for now as it's about the only thing that is testing the BandMatrix class. 19 20template<typename _MatrixType> class UpperBidiagonalization 21{ 22 public: 23 24 typedef _MatrixType MatrixType; 25 enum { 26 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 27 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 28 ColsAtCompileTimeMinusOne = internal::decrement_size<ColsAtCompileTime>::ret 29 }; 30 typedef typename MatrixType::Scalar Scalar; 31 typedef typename MatrixType::RealScalar RealScalar; 32 typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 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< 39 const MatrixType, 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>, 45 OnTheRight 46 > HouseholderVSequenceType; 47 48 /** 49 * \brief Default Constructor. 50 * 51 * The default constructor is useful in cases in which the user intends to 52 * perform decompositions via Bidiagonalization::compute(const MatrixType&). 53 */ 54 UpperBidiagonalization() : m_householder(), m_bidiagonal(), m_isInitialized(false) {} 55 56 explicit UpperBidiagonalization(const MatrixType& matrix) 57 : m_householder(matrix.rows(), matrix.cols()), 58 m_bidiagonal(matrix.cols(), matrix.cols()), 59 m_isInitialized(false) 60 { 61 compute(matrix); 62 } 63 64 UpperBidiagonalization& compute(const MatrixType& matrix); 65 UpperBidiagonalization& computeUnblocked(const MatrixType& matrix); 66 67 const MatrixType& householder() const { return m_householder; } 68 const BidiagonalType& bidiagonal() const { return m_bidiagonal; } 69 70 const HouseholderUSequenceType householderU() const 71 { 72 eigen_assert(m_isInitialized && "UpperBidiagonalization is not initialized."); 73 return HouseholderUSequenceType(m_householder, m_householder.diagonal().conjugate()); 74 } 75 76 const HouseholderVSequenceType householderV() // const here gives nasty errors and i'm lazy 77 { 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) 81 .setShift(1); 82 } 83 84 protected: 85 MatrixType m_householder; 86 BidiagonalType m_bidiagonal; 87 bool m_isInitialized; 88}; 89 90// Standard upper bidiagonalization without fancy optimizations 91// This version should be faster for small matrix size 92template<typename MatrixType> 93void upperbidiagonalization_inplace_unblocked(MatrixType& mat, 94 typename MatrixType::RealScalar *diagonal, 95 typename MatrixType::RealScalar *upper_diagonal, 96 typename MatrixType::Scalar* tempData = 0) 97{ 98 typedef typename MatrixType::Scalar Scalar; 99 100 Index rows = mat.rows(); 101 Index cols = mat.cols(); 102 103 typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixType::MaxRowsAtCompileTime,1> TempType; 104 TempType tempVector; 105 if(tempData==0) 106 { 107 tempVector.resize(rows); 108 tempData = tempVector.data(); 109 } 110 111 for (Index k = 0; /* breaks at k==cols-1 below */ ; ++k) 112 { 113 Index remainingRows = rows - k; 114 Index remainingCols = cols - k - 1; 115 116 // construct left householder transform in-place in A 117 mat.col(k).tail(remainingRows) 118 .makeHouseholderInPlace(mat.coeffRef(k,k), diagonal[k]); 119 // apply householder transform to remaining part of A on the left 120 mat.bottomRightCorner(remainingRows, remainingCols) 121 .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), mat.coeff(k,k), tempData); 122 123 if(k == cols-1) break; 124 125 // construct right householder transform in-place in mat 126 mat.row(k).tail(remainingCols) 127 .makeHouseholderInPlace(mat.coeffRef(k,k+1), upper_diagonal[k]); 128 // apply householder transform to remaining part of mat on the left 129 mat.bottomRightCorner(remainingRows-1, remainingCols) 130 .applyHouseholderOnTheRight(mat.row(k).tail(remainingCols-1).transpose(), mat.coeff(k,k+1), tempData); 131 } 132} 133 134/** \internal 135 * Helper routine for the block reduction to upper bidiagonal form. 136 * 137 * Let's partition the matrix A: 138 * 139 * | A00 A01 | 140 * A = | | 141 * | A10 A11 | 142 * 143 * This function reduces to bidiagonal form the left \c rows x \a blockSize vertical panel [A00/A10] 144 * and the \a blockSize x \c cols horizontal panel [A00 A01] of the matrix \a A. The bottom-right block A11 145 * is updated using matrix-matrix products: 146 * A22 -= V * Y^T - X * U^T 147 * where V and U contains the left and right Householder vectors. U and V are stored in A10, and A01 148 * respectively, and the update matrices X and Y are computed during the reduction. 149 * 150 */ 151template<typename MatrixType> 152void upperbidiagonalization_blocked_helper(MatrixType& A, 153 typename MatrixType::RealScalar *diagonal, 154 typename MatrixType::RealScalar *upper_diagonal, 155 Index bs, 156 Ref<Matrix<typename MatrixType::Scalar, Dynamic, Dynamic, 157 traits<MatrixType>::Flags & RowMajorBit> > X, 158 Ref<Matrix<typename MatrixType::Scalar, Dynamic, Dynamic, 159 traits<MatrixType>::Flags & RowMajorBit> > Y) 160{ 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; 168 169 Index brows = A.rows(); 170 Index bcols = A.cols(); 171 172 Scalar tau_u, tau_u_prev(0), tau_v; 173 174 for(Index k = 0; k < bs; ++k) 175 { 176 Index remainingRows = brows - k; 177 Index remainingCols = bcols - k - 1; 178 179 SubMatType X_k1( X.block(k,0, remainingRows,k) ); 180 SubMatType V_k1( A.block(k,0, remainingRows,k) ); 181 182 // 1 - update the k-th column of A 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); 186 187 // 2 - construct left Householder transform in-place 188 v_k.makeHouseholderInPlace(tau_v, diagonal[k]); 189 190 if(k+1<bcols) 191 { 192 SubMatType Y_k ( Y.block(k+1,0, remainingCols, k+1) ); 193 SubMatType U_k1 ( A.block(0,k+1, k,remainingCols) ); 194 195 // this eases the application of Householder transforAions 196 // A(k,k) will store tau_v later 197 A(k,k) = Scalar(1); 198 199 // 3 - Compute y_k^T = tau_v * ( A^T*v_k - Y_k-1*V_k-1^T*v_k - U_k-1*X_k-1^T*v_k ) 200 { 201 SubColumnType y_k( Y.col(k).tail(remainingCols) ); 202 203 // let's use the begining of column k of Y as a temporary vector 204 SubColumnType tmp( Y.col(k).head(k) ); 205 y_k.noalias() = A.block(k,k+1, remainingRows,remainingCols).adjoint() * v_k; // bottleneck 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); 211 } 212 213 // 4 - update k-th row of A (it will become u_k) 214 SubRowType u_k( A.row(k).tail(remainingCols) ); 215 u_k = u_k.conjugate(); 216 { 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(); 219 } 220 221 // 5 - construct right Householder transform in-place 222 u_k.makeHouseholderInPlace(tau_u, upper_diagonal[k]); 223 224 // this eases the application of Householder transformations 225 // A(k,k+1) will store tau_u later 226 A(k,k+1) = Scalar(1); 227 228 // 6 - Compute x_k = tau_u * ( A*u_k - X_k-1*U_k-1^T*u_k - V_k*Y_k^T*u_k ) 229 { 230 SubColumnType x_k ( X.col(k).tail(remainingRows-1) ); 231 232 // let's use the begining of column k of X as a temporary vectors 233 // note that tmp0 and tmp1 overlaps 234 SubColumnType tmp0 ( X.col(k).head(k) ), 235 tmp1 ( X.col(k).head(k+1) ); 236 237 x_k.noalias() = A.block(k+1,k+1, remainingRows-1,remainingCols) * u_k.transpose(); // bottleneck 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(); 245 } 246 247 if(k>0) A.coeffRef(k-1,k) = tau_u_prev; 248 tau_u_prev = tau_u; 249 } 250 else 251 A.coeffRef(k-1,k) = tau_u_prev; 252 253 A.coeffRef(k,k) = tau_v; 254 } 255 256 if(bs<bcols) 257 A.coeffRef(bs-1,bs) = tau_u_prev; 258 259 // update A22 260 if(bcols>bs && brows>bs) 261 { 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); 266 A01(bs-1,0) = 1; 267 A11.noalias() -= A10 * Y.topLeftCorner(bcols,bs).bottomRows(bcols-bs).adjoint(); 268 A11.noalias() -= X.topLeftCorner(brows,bs).bottomRows(brows-bs) * A01; 269 A01(bs-1,0) = tmp; 270 } 271} 272 273/** \internal 274 * 275 * Implementation of a block-bidiagonal reduction. 276 * It is based on the following paper: 277 * The Design of a Parallel Dense Linear Algebra Software Library: Reduction to Hessenberg, Tridiagonal, and Bidiagonal Form. 278 * by Jaeyoung Choi, Jack J. Dongarra, David W. Walker. (1995) 279 * section 3.3 280 */ 281template<typename MatrixType, typename BidiagType> 282void upperbidiagonalization_inplace_blocked(MatrixType& A, BidiagType& bidiagonal, 283 Index maxBlockSize=32, 284 typename MatrixType::Scalar* /*tempData*/ = 0) 285{ 286 typedef typename MatrixType::Scalar Scalar; 287 typedef Block<MatrixType,Dynamic,Dynamic> BlockType; 288 289 Index rows = A.rows(); 290 Index cols = A.cols(); 291 Index size = (std::min)(rows, cols); 292 293 // X and Y are work space 294 enum { StorageOrder = traits<MatrixType>::Flags & RowMajorBit }; 295 Matrix<Scalar, 296 MatrixType::RowsAtCompileTime, 297 Dynamic, 298 StorageOrder, 299 MatrixType::MaxRowsAtCompileTime> X(rows,maxBlockSize); 300 Matrix<Scalar, 301 MatrixType::ColsAtCompileTime, 302 Dynamic, 303 StorageOrder, 304 MatrixType::MaxColsAtCompileTime> Y(cols,maxBlockSize); 305 Index blockSize = (std::min)(maxBlockSize,size); 306 307 Index k = 0; 308 for(k = 0; k < size; k += blockSize) 309 { 310 Index bs = (std::min)(size-k,blockSize); // actual size of the block 311 Index brows = rows - k; // rows of the block 312 Index bcols = cols - k; // columns of the block 313 314 // partition the matrix A: 315 // 316 // | A00 A01 A02 | 317 // | | 318 // A = | A10 A11 A12 | 319 // | | 320 // | A20 A21 A22 | 321 // 322 // where A11 is a bs x bs diagonal block, 323 // and let: 324 // | A11 A12 | 325 // B = | | 326 // | A21 A22 | 327 328 BlockType B = A.block(k,k,brows,bcols); 329 330 // This stage performs the bidiagonalization of A11, A21, A12, and updating of A22. 331 // Finally, the algorithm continue on the updated A22. 332 // 333 // However, if B is too small, or A22 empty, then let's use an unblocked strategy 334 if(k+bs==cols || bcols<48) // somewhat arbitrary threshold 335 { 336 upperbidiagonalization_inplace_unblocked(B, 337 &(bidiagonal.template diagonal<0>().coeffRef(k)), 338 &(bidiagonal.template diagonal<1>().coeffRef(k)), 339 X.data() 340 ); 341 break; // We're done 342 } 343 else 344 { 345 upperbidiagonalization_blocked_helper<BlockType>( B, 346 &(bidiagonal.template diagonal<0>().coeffRef(k)), 347 &(bidiagonal.template diagonal<1>().coeffRef(k)), 348 bs, 349 X.topLeftCorner(brows,bs), 350 Y.topLeftCorner(bcols,bs) 351 ); 352 } 353 } 354} 355 356template<typename _MatrixType> 357UpperBidiagonalization<_MatrixType>& UpperBidiagonalization<_MatrixType>::computeUnblocked(const _MatrixType& matrix) 358{ 359 Index rows = matrix.rows(); 360 Index cols = matrix.cols(); 361 EIGEN_ONLY_USED_FOR_DEBUG(cols); 362 363 eigen_assert(rows >= cols && "UpperBidiagonalization is only for Arices satisfying rows>=cols."); 364 365 m_householder = matrix; 366 367 ColVectorType temp(rows); 368 369 upperbidiagonalization_inplace_unblocked(m_householder, 370 &(m_bidiagonal.template diagonal<0>().coeffRef(0)), 371 &(m_bidiagonal.template diagonal<1>().coeffRef(0)), 372 temp.data()); 373 374 m_isInitialized = true; 375 return *this; 376} 377 378template<typename _MatrixType> 379UpperBidiagonalization<_MatrixType>& UpperBidiagonalization<_MatrixType>::compute(const _MatrixType& matrix) 380{ 381 Index rows = matrix.rows(); 382 Index cols = matrix.cols(); 383 EIGEN_ONLY_USED_FOR_DEBUG(rows); 384 EIGEN_ONLY_USED_FOR_DEBUG(cols); 385 386 eigen_assert(rows >= cols && "UpperBidiagonalization is only for Arices satisfying rows>=cols."); 387 388 m_householder = matrix; 389 upperbidiagonalization_inplace_blocked(m_householder, m_bidiagonal); 390 391 m_isInitialized = true; 392 return *this; 393} 394 395#if 0 396/** \return the Householder QR decomposition of \c *this. 397 * 398 * \sa class Bidiagonalization 399 */ 400template<typename Derived> 401const UpperBidiagonalization<typename MatrixBase<Derived>::PlainObject> 402MatrixBase<Derived>::bidiagonalization() const 403{ 404 return UpperBidiagonalization<PlainObject>(eval()); 405} 406#endif 407 408} // end namespace internal 409 410} // end namespace Eigen 411 412#endif // EIGEN_BIDIAGONALIZATION_H 413