![]() |
Eigen
3.3.3
|
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 // Copyright (C) 2013-2014 Gael Guennebaud <gael.guennebaud@inria.fr> 00006 // 00007 // This Source Code Form is subject to the terms of the Mozilla 00008 // Public License v. 2.0. If a copy of the MPL was not distributed 00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00010 00011 #ifndef EIGEN_BIDIAGONALIZATION_H 00012 #define EIGEN_BIDIAGONALIZATION_H 00013 00014 namespace Eigen { 00015 00016 namespace internal { 00017 // UpperBidiagonalization will probably be replaced by a Bidiagonalization class, don't want to make it stable API. 00018 // At the same time, it's useful to keep for now as it's about the only thing that is testing the BandMatrix class. 00019 00020 template<typename _MatrixType> class UpperBidiagonalization 00021 { 00022 public: 00023 00024 typedef _MatrixType MatrixType; 00025 enum { 00026 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00027 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00028 ColsAtCompileTimeMinusOne = internal::decrement_size<ColsAtCompileTime>::ret 00029 }; 00030 typedef typename MatrixType::Scalar Scalar; 00031 typedef typename MatrixType::RealScalar RealScalar; 00032 typedef Eigen::Index Index; 00033 typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType; 00034 typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType; 00035 typedef BandMatrix<RealScalar, ColsAtCompileTime, ColsAtCompileTime, 1, 0, RowMajor> BidiagonalType; 00036 typedef Matrix<Scalar, ColsAtCompileTime, 1> DiagVectorType; 00037 typedef Matrix<Scalar, ColsAtCompileTimeMinusOne, 1> SuperDiagVectorType; 00038 typedef HouseholderSequence< 00039 const MatrixType, 00040 const typename internal::remove_all<typename Diagonal<const MatrixType,0>::ConjugateReturnType>::type 00041 > HouseholderUSequenceType; 00042 typedef HouseholderSequence< 00043 const typename internal::remove_all<typename MatrixType::ConjugateReturnType>::type, 00044 Diagonal<const MatrixType,1>, 00045 OnTheRight 00046 > HouseholderVSequenceType; 00047 00054 UpperBidiagonalization() : m_householder(), m_bidiagonal(), m_isInitialized(false) {} 00055 00056 explicit UpperBidiagonalization(const MatrixType& matrix) 00057 : m_householder(matrix.rows(), matrix.cols()), 00058 m_bidiagonal(matrix.cols(), matrix.cols()), 00059 m_isInitialized(false) 00060 { 00061 compute(matrix); 00062 } 00063 00064 UpperBidiagonalization& compute(const MatrixType& matrix); 00065 UpperBidiagonalization& computeUnblocked(const MatrixType& matrix); 00066 00067 const MatrixType& householder() const { return m_householder; } 00068 const BidiagonalType& bidiagonal() const { return m_bidiagonal; } 00069 00070 const HouseholderUSequenceType householderU() const 00071 { 00072 eigen_assert(m_isInitialized && "UpperBidiagonalization is not initialized."); 00073 return HouseholderUSequenceType(m_householder, m_householder.diagonal().conjugate()); 00074 } 00075 00076 const HouseholderVSequenceType householderV() // const here gives nasty errors and i'm lazy 00077 { 00078 eigen_assert(m_isInitialized && "UpperBidiagonalization is not initialized."); 00079 return HouseholderVSequenceType(m_householder.conjugate(), m_householder.const_derived().template diagonal<1>()) 00080 .setLength(m_householder.cols()-1) 00081 .setShift(1); 00082 } 00083 00084 protected: 00085 MatrixType m_householder; 00086 BidiagonalType m_bidiagonal; 00087 bool m_isInitialized; 00088 }; 00089 00090 // Standard upper bidiagonalization without fancy optimizations 00091 // This version should be faster for small matrix size 00092 template<typename MatrixType> 00093 void upperbidiagonalization_inplace_unblocked(MatrixType& mat, 00094 typename MatrixType::RealScalar *diagonal, 00095 typename MatrixType::RealScalar *upper_diagonal, 00096 typename MatrixType::Scalar* tempData = 0) 00097 { 00098 typedef typename MatrixType::Scalar Scalar; 00099 00100 Index rows = mat.rows(); 00101 Index cols = mat.cols(); 00102 00103 typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixType::MaxRowsAtCompileTime,1> TempType; 00104 TempType tempVector; 00105 if(tempData==0) 00106 { 00107 tempVector.resize(rows); 00108 tempData = tempVector.data(); 00109 } 00110 00111 for (Index k = 0; /* breaks at k==cols-1 below */ ; ++k) 00112 { 00113 Index remainingRows = rows - k; 00114 Index remainingCols = cols - k - 1; 00115 00116 // construct left householder transform in-place in A 00117 mat.col(k).tail(remainingRows) 00118 .makeHouseholderInPlace(mat.coeffRef(k,k), diagonal[k]); 00119 // apply householder transform to remaining part of A on the left 00120 mat.bottomRightCorner(remainingRows, remainingCols) 00121 .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), mat.coeff(k,k), tempData); 00122 00123 if(k == cols-1) break; 00124 00125 // construct right householder transform in-place in mat 00126 mat.row(k).tail(remainingCols) 00127 .makeHouseholderInPlace(mat.coeffRef(k,k+1), upper_diagonal[k]); 00128 // apply householder transform to remaining part of mat on the left 00129 mat.bottomRightCorner(remainingRows-1, remainingCols) 00130 .applyHouseholderOnTheRight(mat.row(k).tail(remainingCols-1).transpose(), mat.coeff(k,k+1), tempData); 00131 } 00132 } 00133 00151 template<typename MatrixType> 00152 void upperbidiagonalization_blocked_helper(MatrixType& A, 00153 typename MatrixType::RealScalar *diagonal, 00154 typename MatrixType::RealScalar *upper_diagonal, 00155 Index bs, 00156 Ref<Matrix<typename MatrixType::Scalar, Dynamic, Dynamic, 00157 traits<MatrixType>::Flags & RowMajorBit> > X, 00158 Ref<Matrix<typename MatrixType::Scalar, Dynamic, Dynamic, 00159 traits<MatrixType>::Flags & RowMajorBit> > Y) 00160 { 00161 typedef typename MatrixType::Scalar Scalar; 00162 enum { StorageOrder = traits<MatrixType>::Flags & RowMajorBit }; 00163 typedef InnerStride<int(StorageOrder) == int(ColMajor) ? 1 : Dynamic> ColInnerStride; 00164 typedef InnerStride<int(StorageOrder) == int(ColMajor) ? Dynamic : 1> RowInnerStride; 00165 typedef Ref<Matrix<Scalar, Dynamic, 1>, 0, ColInnerStride> SubColumnType; 00166 typedef Ref<Matrix<Scalar, 1, Dynamic>, 0, RowInnerStride> SubRowType; 00167 typedef Ref<Matrix<Scalar, Dynamic, Dynamic, StorageOrder > > SubMatType; 00168 00169 Index brows = A.rows(); 00170 Index bcols = A.cols(); 00171 00172 Scalar tau_u, tau_u_prev(0), tau_v; 00173 00174 for(Index k = 0; k < bs; ++k) 00175 { 00176 Index remainingRows = brows - k; 00177 Index remainingCols = bcols - k - 1; 00178 00179 SubMatType X_k1( X.block(k,0, remainingRows,k) ); 00180 SubMatType V_k1( A.block(k,0, remainingRows,k) ); 00181 00182 // 1 - update the k-th column of A 00183 SubColumnType v_k = A.col(k).tail(remainingRows); 00184 v_k -= V_k1 * Y.row(k).head(k).adjoint(); 00185 if(k) v_k -= X_k1 * A.col(k).head(k); 00186 00187 // 2 - construct left Householder transform in-place 00188 v_k.makeHouseholderInPlace(tau_v, diagonal[k]); 00189 00190 if(k+1<bcols) 00191 { 00192 SubMatType Y_k ( Y.block(k+1,0, remainingCols, k+1) ); 00193 SubMatType U_k1 ( A.block(0,k+1, k,remainingCols) ); 00194 00195 // this eases the application of Householder transforAions 00196 // A(k,k) will store tau_v later 00197 A(k,k) = Scalar(1); 00198 00199 // 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 ) 00200 { 00201 SubColumnType y_k( Y.col(k).tail(remainingCols) ); 00202 00203 // let's use the begining of column k of Y as a temporary vector 00204 SubColumnType tmp( Y.col(k).head(k) ); 00205 y_k.noalias() = A.block(k,k+1, remainingRows,remainingCols).adjoint() * v_k; // bottleneck 00206 tmp.noalias() = V_k1.adjoint() * v_k; 00207 y_k.noalias() -= Y_k.leftCols(k) * tmp; 00208 tmp.noalias() = X_k1.adjoint() * v_k; 00209 y_k.noalias() -= U_k1.adjoint() * tmp; 00210 y_k *= numext::conj(tau_v); 00211 } 00212 00213 // 4 - update k-th row of A (it will become u_k) 00214 SubRowType u_k( A.row(k).tail(remainingCols) ); 00215 u_k = u_k.conjugate(); 00216 { 00217 u_k -= Y_k * A.row(k).head(k+1).adjoint(); 00218 if(k) u_k -= U_k1.adjoint() * X.row(k).head(k).adjoint(); 00219 } 00220 00221 // 5 - construct right Householder transform in-place 00222 u_k.makeHouseholderInPlace(tau_u, upper_diagonal[k]); 00223 00224 // this eases the application of Householder transformations 00225 // A(k,k+1) will store tau_u later 00226 A(k,k+1) = Scalar(1); 00227 00228 // 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 ) 00229 { 00230 SubColumnType x_k ( X.col(k).tail(remainingRows-1) ); 00231 00232 // let's use the begining of column k of X as a temporary vectors 00233 // note that tmp0 and tmp1 overlaps 00234 SubColumnType tmp0 ( X.col(k).head(k) ), 00235 tmp1 ( X.col(k).head(k+1) ); 00236 00237 x_k.noalias() = A.block(k+1,k+1, remainingRows-1,remainingCols) * u_k.transpose(); // bottleneck 00238 tmp0.noalias() = U_k1 * u_k.transpose(); 00239 x_k.noalias() -= X_k1.bottomRows(remainingRows-1) * tmp0; 00240 tmp1.noalias() = Y_k.adjoint() * u_k.transpose(); 00241 x_k.noalias() -= A.block(k+1,0, remainingRows-1,k+1) * tmp1; 00242 x_k *= numext::conj(tau_u); 00243 tau_u = numext::conj(tau_u); 00244 u_k = u_k.conjugate(); 00245 } 00246 00247 if(k>0) A.coeffRef(k-1,k) = tau_u_prev; 00248 tau_u_prev = tau_u; 00249 } 00250 else 00251 A.coeffRef(k-1,k) = tau_u_prev; 00252 00253 A.coeffRef(k,k) = tau_v; 00254 } 00255 00256 if(bs<bcols) 00257 A.coeffRef(bs-1,bs) = tau_u_prev; 00258 00259 // update A22 00260 if(bcols>bs && brows>bs) 00261 { 00262 SubMatType A11( A.bottomRightCorner(brows-bs,bcols-bs) ); 00263 SubMatType A10( A.block(bs,0, brows-bs,bs) ); 00264 SubMatType A01( A.block(0,bs, bs,bcols-bs) ); 00265 Scalar tmp = A01(bs-1,0); 00266 A01(bs-1,0) = 1; 00267 A11.noalias() -= A10 * Y.topLeftCorner(bcols,bs).bottomRows(bcols-bs).adjoint(); 00268 A11.noalias() -= X.topLeftCorner(brows,bs).bottomRows(brows-bs) * A01; 00269 A01(bs-1,0) = tmp; 00270 } 00271 } 00272 00281 template<typename MatrixType, typename BidiagType> 00282 void upperbidiagonalization_inplace_blocked(MatrixType& A, BidiagType& bidiagonal, 00283 Index maxBlockSize=32, 00284 typename MatrixType::Scalar* /*tempData*/ = 0) 00285 { 00286 typedef typename MatrixType::Scalar Scalar; 00287 typedef Block<MatrixType,Dynamic,Dynamic> BlockType; 00288 00289 Index rows = A.rows(); 00290 Index cols = A.cols(); 00291 Index size = (std::min)(rows, cols); 00292 00293 // X and Y are work space 00294 enum { StorageOrder = traits<MatrixType>::Flags & RowMajorBit }; 00295 Matrix<Scalar, 00296 MatrixType::RowsAtCompileTime, 00297 Dynamic, 00298 StorageOrder, 00299 MatrixType::MaxRowsAtCompileTime> X(rows,maxBlockSize); 00300 Matrix<Scalar, 00301 MatrixType::ColsAtCompileTime, 00302 Dynamic, 00303 StorageOrder, 00304 MatrixType::MaxColsAtCompileTime> Y(cols,maxBlockSize); 00305 Index blockSize = (std::min)(maxBlockSize,size); 00306 00307 Index k = 0; 00308 for(k = 0; k < size; k += blockSize) 00309 { 00310 Index bs = (std::min)(size-k,blockSize); // actual size of the block 00311 Index brows = rows - k; // rows of the block 00312 Index bcols = cols - k; // columns of the block 00313 00314 // partition the matrix A: 00315 // 00316 // | A00 A01 A02 | 00317 // | | 00318 // A = | A10 A11 A12 | 00319 // | | 00320 // | A20 A21 A22 | 00321 // 00322 // where A11 is a bs x bs diagonal block, 00323 // and let: 00324 // | A11 A12 | 00325 // B = | | 00326 // | A21 A22 | 00327 00328 BlockType B = A.block(k,k,brows,bcols); 00329 00330 // This stage performs the bidiagonalization of A11, A21, A12, and updating of A22. 00331 // Finally, the algorithm continue on the updated A22. 00332 // 00333 // However, if B is too small, or A22 empty, then let's use an unblocked strategy 00334 if(k+bs==cols || bcols<48) // somewhat arbitrary threshold 00335 { 00336 upperbidiagonalization_inplace_unblocked(B, 00337 &(bidiagonal.template diagonal<0>().coeffRef(k)), 00338 &(bidiagonal.template diagonal<1>().coeffRef(k)), 00339 X.data() 00340 ); 00341 break; // We're done 00342 } 00343 else 00344 { 00345 upperbidiagonalization_blocked_helper<BlockType>( B, 00346 &(bidiagonal.template diagonal<0>().coeffRef(k)), 00347 &(bidiagonal.template diagonal<1>().coeffRef(k)), 00348 bs, 00349 X.topLeftCorner(brows,bs), 00350 Y.topLeftCorner(bcols,bs) 00351 ); 00352 } 00353 } 00354 } 00355 00356 template<typename _MatrixType> 00357 UpperBidiagonalization<_MatrixType>& UpperBidiagonalization<_MatrixType>::computeUnblocked(const _MatrixType& matrix) 00358 { 00359 Index rows = matrix.rows(); 00360 Index cols = matrix.cols(); 00361 EIGEN_ONLY_USED_FOR_DEBUG(cols); 00362 00363 eigen_assert(rows >= cols && "UpperBidiagonalization is only for Arices satisfying rows>=cols."); 00364 00365 m_householder = matrix; 00366 00367 ColVectorType temp(rows); 00368 00369 upperbidiagonalization_inplace_unblocked(m_householder, 00370 &(m_bidiagonal.template diagonal<0>().coeffRef(0)), 00371 &(m_bidiagonal.template diagonal<1>().coeffRef(0)), 00372 temp.data()); 00373 00374 m_isInitialized = true; 00375 return *this; 00376 } 00377 00378 template<typename _MatrixType> 00379 UpperBidiagonalization<_MatrixType>& UpperBidiagonalization<_MatrixType>::compute(const _MatrixType& matrix) 00380 { 00381 Index rows = matrix.rows(); 00382 Index cols = matrix.cols(); 00383 EIGEN_ONLY_USED_FOR_DEBUG(rows); 00384 EIGEN_ONLY_USED_FOR_DEBUG(cols); 00385 00386 eigen_assert(rows >= cols && "UpperBidiagonalization is only for Arices satisfying rows>=cols."); 00387 00388 m_householder = matrix; 00389 upperbidiagonalization_inplace_blocked(m_householder, m_bidiagonal); 00390 00391 m_isInitialized = true; 00392 return *this; 00393 } 00394 00395 #if 0 00396 00400 template<typename Derived> 00401 const UpperBidiagonalization<typename MatrixBase<Derived>::PlainObject> 00402 MatrixBase<Derived>::bidiagonalization() const 00403 { 00404 return UpperBidiagonalization<PlainObject>(eval()); 00405 } 00406 #endif 00407 00408 } // end namespace internal 00409 00410 } // end namespace Eigen 00411 00412 #endif // EIGEN_BIDIAGONALIZATION_H