Important changes to repositories hosted on mbed.com
Mbed hosted mercurial repositories are deprecated and are due to be permanently deleted in July 2026.
To keep a copy of this software download the repository Zip archive or clone locally using Mercurial.
It is also possible to export all your personal repositories from the account settings page.
Dependents: Eigen_test Odometry_test AttitudeEstimation_usingTicker MPU9250_Quaternion_Binary_Serial ... more
FullPivHouseholderQR.h
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr> 00005 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> 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_FULLPIVOTINGHOUSEHOLDERQR_H 00012 #define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H 00013 00014 namespace Eigen { 00015 00016 namespace internal { 00017 00018 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType; 00019 00020 template<typename MatrixType> 00021 struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> > 00022 { 00023 typedef typename MatrixType::PlainObject ReturnType; 00024 }; 00025 00026 } 00027 00028 /** \ingroup QR_Module 00029 * 00030 * \class FullPivHouseholderQR 00031 * 00032 * \brief Householder rank-revealing QR decomposition of a matrix with full pivoting 00033 * 00034 * \param MatrixType the type of the matrix of which we are computing the QR decomposition 00035 * 00036 * This class performs a rank-revealing QR decomposition of a matrix \b A into matrices \b P, \b Q and \b R 00037 * such that 00038 * \f[ 00039 * \mathbf{A} \, \mathbf{P} = \mathbf{Q} \, \mathbf{R} 00040 * \f] 00041 * by using Householder transformations. Here, \b P is a permutation matrix, \b Q a unitary matrix and \b R an 00042 * upper triangular matrix. 00043 * 00044 * This decomposition performs a very prudent full pivoting in order to be rank-revealing and achieve optimal 00045 * numerical stability. The trade-off is that it is slower than HouseholderQR and ColPivHouseholderQR. 00046 * 00047 * \sa MatrixBase::fullPivHouseholderQr() 00048 */ 00049 template<typename _MatrixType> class FullPivHouseholderQR 00050 { 00051 public: 00052 00053 typedef _MatrixType MatrixType; 00054 enum { 00055 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00056 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00057 Options = MatrixType::Options, 00058 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00059 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 00060 }; 00061 typedef typename MatrixType::Scalar Scalar; 00062 typedef typename MatrixType::RealScalar RealScalar; 00063 typedef typename MatrixType::Index Index; 00064 typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType; 00065 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; 00066 typedef Matrix<Index, 1, 00067 EIGEN_SIZE_MIN_PREFER_DYNAMIC(ColsAtCompileTime,RowsAtCompileTime), RowMajor, 1, 00068 EIGEN_SIZE_MIN_PREFER_FIXED(MaxColsAtCompileTime,MaxRowsAtCompileTime)> IntDiagSizeVectorType ; 00069 typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType ; 00070 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType; 00071 typedef typename internal::plain_col_type<MatrixType>::type ColVectorType; 00072 00073 /** \brief Default Constructor. 00074 * 00075 * The default constructor is useful in cases in which the user intends to 00076 * perform decompositions via FullPivHouseholderQR::compute(const MatrixType&). 00077 */ 00078 FullPivHouseholderQR() 00079 : m_qr(), 00080 m_hCoeffs(), 00081 m_rows_transpositions(), 00082 m_cols_transpositions(), 00083 m_cols_permutation(), 00084 m_temp(), 00085 m_isInitialized(false), 00086 m_usePrescribedThreshold(false) {} 00087 00088 /** \brief Default Constructor with memory preallocation 00089 * 00090 * Like the default constructor but with preallocation of the internal data 00091 * according to the specified problem \a size. 00092 * \sa FullPivHouseholderQR() 00093 */ 00094 FullPivHouseholderQR(Index rows, Index cols) 00095 : m_qr(rows, cols), 00096 m_hCoeffs((std::min)(rows,cols)), 00097 m_rows_transpositions((std::min)(rows,cols)), 00098 m_cols_transpositions((std::min)(rows,cols)), 00099 m_cols_permutation(cols), 00100 m_temp(cols), 00101 m_isInitialized(false), 00102 m_usePrescribedThreshold(false) {} 00103 00104 /** \brief Constructs a QR factorization from a given matrix 00105 * 00106 * This constructor computes the QR factorization of the matrix \a matrix by calling 00107 * the method compute(). It is a short cut for: 00108 * 00109 * \code 00110 * FullPivHouseholderQR<MatrixType> qr(matrix.rows(), matrix.cols()); 00111 * qr.compute(matrix); 00112 * \endcode 00113 * 00114 * \sa compute() 00115 */ 00116 FullPivHouseholderQR(const MatrixType& matrix) 00117 : m_qr(matrix.rows(), matrix.cols()), 00118 m_hCoeffs((std::min)(matrix.rows(), matrix.cols())), 00119 m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())), 00120 m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())), 00121 m_cols_permutation(matrix.cols()), 00122 m_temp(matrix.cols()), 00123 m_isInitialized(false), 00124 m_usePrescribedThreshold(false) 00125 { 00126 compute(matrix); 00127 } 00128 00129 /** This method finds a solution x to the equation Ax=b, where A is the matrix of which 00130 * \c *this is the QR decomposition. 00131 * 00132 * \param b the right-hand-side of the equation to solve. 00133 * 00134 * \returns the exact or least-square solution if the rank is greater or equal to the number of columns of A, 00135 * and an arbitrary solution otherwise. 00136 * 00137 * \note The case where b is a matrix is not yet implemented. Also, this 00138 * code is space inefficient. 00139 * 00140 * \note_about_checking_solutions 00141 * 00142 * \note_about_arbitrary_choice_of_solution 00143 * 00144 * Example: \include FullPivHouseholderQR_solve.cpp 00145 * Output: \verbinclude FullPivHouseholderQR_solve.out 00146 */ 00147 template<typename Rhs> 00148 inline const internal::solve_retval<FullPivHouseholderQR, Rhs> 00149 solve(const MatrixBase<Rhs>& b) const 00150 { 00151 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00152 return internal::solve_retval<FullPivHouseholderQR, Rhs>(*this, b.derived()); 00153 } 00154 00155 /** \returns Expression object representing the matrix Q 00156 */ 00157 MatrixQReturnType matrixQ (void) const; 00158 00159 /** \returns a reference to the matrix where the Householder QR decomposition is stored 00160 */ 00161 const MatrixType& matrixQR () const 00162 { 00163 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00164 return m_qr; 00165 } 00166 00167 FullPivHouseholderQR& compute(const MatrixType& matrix); 00168 00169 /** \returns a const reference to the column permutation matrix */ 00170 const PermutationType & colsPermutation () const 00171 { 00172 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00173 return m_cols_permutation; 00174 } 00175 00176 /** \returns a const reference to the vector of indices representing the rows transpositions */ 00177 const IntDiagSizeVectorType & rowsTranspositions () const 00178 { 00179 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00180 return m_rows_transpositions; 00181 } 00182 00183 /** \returns the absolute value of the determinant of the matrix of which 00184 * *this is the QR decomposition. It has only linear complexity 00185 * (that is, O(n) where n is the dimension of the square matrix) 00186 * as the QR decomposition has already been computed. 00187 * 00188 * \note This is only for square matrices. 00189 * 00190 * \warning a determinant can be very big or small, so for matrices 00191 * of large enough dimension, there is a risk of overflow/underflow. 00192 * One way to work around that is to use logAbsDeterminant() instead. 00193 * 00194 * \sa logAbsDeterminant(), MatrixBase::determinant() 00195 */ 00196 typename MatrixType::RealScalar absDeterminant () const; 00197 00198 /** \returns the natural log of the absolute value of the determinant of the matrix of which 00199 * *this is the QR decomposition. It has only linear complexity 00200 * (that is, O(n) where n is the dimension of the square matrix) 00201 * as the QR decomposition has already been computed. 00202 * 00203 * \note This is only for square matrices. 00204 * 00205 * \note This method is useful to work around the risk of overflow/underflow that's inherent 00206 * to determinant computation. 00207 * 00208 * \sa absDeterminant(), MatrixBase::determinant() 00209 */ 00210 typename MatrixType::RealScalar logAbsDeterminant () const; 00211 00212 /** \returns the rank of the matrix of which *this is the QR decomposition. 00213 * 00214 * \note This method has to determine which pivots should be considered nonzero. 00215 * For that, it uses the threshold value that you can control by calling 00216 * setThreshold(const RealScalar&). 00217 */ 00218 inline Index rank () const 00219 { 00220 using std::abs; 00221 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00222 RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold(); 00223 Index result = 0; 00224 for(Index i = 0; i < m_nonzero_pivots; ++i) 00225 result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold); 00226 return result; 00227 } 00228 00229 /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition. 00230 * 00231 * \note This method has to determine which pivots should be considered nonzero. 00232 * For that, it uses the threshold value that you can control by calling 00233 * setThreshold(const RealScalar&). 00234 */ 00235 inline Index dimensionOfKernel () const 00236 { 00237 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00238 return cols() - rank (); 00239 } 00240 00241 /** \returns true if the matrix of which *this is the QR decomposition represents an injective 00242 * linear map, i.e. has trivial kernel; false otherwise. 00243 * 00244 * \note This method has to determine which pivots should be considered nonzero. 00245 * For that, it uses the threshold value that you can control by calling 00246 * setThreshold(const RealScalar&). 00247 */ 00248 inline bool isInjective () const 00249 { 00250 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00251 return rank () == cols(); 00252 } 00253 00254 /** \returns true if the matrix of which *this is the QR decomposition represents a surjective 00255 * linear map; false otherwise. 00256 * 00257 * \note This method has to determine which pivots should be considered nonzero. 00258 * For that, it uses the threshold value that you can control by calling 00259 * setThreshold(const RealScalar&). 00260 */ 00261 inline bool isSurjective () const 00262 { 00263 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00264 return rank () == rows(); 00265 } 00266 00267 /** \returns true if the matrix of which *this is the QR decomposition is invertible. 00268 * 00269 * \note This method has to determine which pivots should be considered nonzero. 00270 * For that, it uses the threshold value that you can control by calling 00271 * setThreshold(const RealScalar&). 00272 */ 00273 inline bool isInvertible () const 00274 { 00275 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00276 return isInjective () && isSurjective (); 00277 } 00278 00279 /** \returns the inverse of the matrix of which *this is the QR decomposition. 00280 * 00281 * \note If this matrix is not invertible, the returned matrix has undefined coefficients. 00282 * Use isInvertible() to first determine whether this matrix is invertible. 00283 */ inline const 00284 internal::solve_retval<FullPivHouseholderQR, typename MatrixType::IdentityReturnType> 00285 inverse () const 00286 { 00287 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00288 return internal::solve_retval<FullPivHouseholderQR,typename MatrixType::IdentityReturnType> 00289 (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols())); 00290 } 00291 00292 inline Index rows() const { return m_qr.rows(); } 00293 inline Index cols() const { return m_qr.cols(); } 00294 00295 /** \returns a const reference to the vector of Householder coefficients used to represent the factor \c Q. 00296 * 00297 * For advanced uses only. 00298 */ 00299 const HCoeffsType& hCoeffs () const { return m_hCoeffs; } 00300 00301 /** Allows to prescribe a threshold to be used by certain methods, such as rank(), 00302 * who need to determine when pivots are to be considered nonzero. This is not used for the 00303 * QR decomposition itself. 00304 * 00305 * When it needs to get the threshold value, Eigen calls threshold(). By default, this 00306 * uses a formula to automatically determine a reasonable threshold. 00307 * Once you have called the present method setThreshold(const RealScalar&), 00308 * your value is used instead. 00309 * 00310 * \param threshold The new value to use as the threshold. 00311 * 00312 * A pivot will be considered nonzero if its absolute value is strictly greater than 00313 * \f$ \vert pivot \vert \leqslant threshold \times \vert maxpivot \vert \f$ 00314 * where maxpivot is the biggest pivot. 00315 * 00316 * If you want to come back to the default behavior, call setThreshold(Default_t) 00317 */ 00318 FullPivHouseholderQR& setThreshold(const RealScalar& threshold) 00319 { 00320 m_usePrescribedThreshold = true; 00321 m_prescribedThreshold = threshold; 00322 return *this; 00323 } 00324 00325 /** Allows to come back to the default behavior, letting Eigen use its default formula for 00326 * determining the threshold. 00327 * 00328 * You should pass the special object Eigen::Default as parameter here. 00329 * \code qr.setThreshold(Eigen::Default); \endcode 00330 * 00331 * See the documentation of setThreshold(const RealScalar&). 00332 */ 00333 FullPivHouseholderQR& setThreshold(Default_t) 00334 { 00335 m_usePrescribedThreshold = false; 00336 return *this; 00337 } 00338 00339 /** Returns the threshold that will be used by certain methods such as rank(). 00340 * 00341 * See the documentation of setThreshold(const RealScalar&). 00342 */ 00343 RealScalar threshold() const 00344 { 00345 eigen_assert(m_isInitialized || m_usePrescribedThreshold); 00346 return m_usePrescribedThreshold ? m_prescribedThreshold 00347 // this formula comes from experimenting (see "LU precision tuning" thread on the list) 00348 // and turns out to be identical to Higham's formula used already in LDLt. 00349 : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize()); 00350 } 00351 00352 /** \returns the number of nonzero pivots in the QR decomposition. 00353 * Here nonzero is meant in the exact sense, not in a fuzzy sense. 00354 * So that notion isn't really intrinsically interesting, but it is 00355 * still useful when implementing algorithms. 00356 * 00357 * \sa rank() 00358 */ 00359 inline Index nonzeroPivots () const 00360 { 00361 eigen_assert(m_isInitialized && "LU is not initialized."); 00362 return m_nonzero_pivots; 00363 } 00364 00365 /** \returns the absolute value of the biggest pivot, i.e. the biggest 00366 * diagonal coefficient of U. 00367 */ 00368 RealScalar maxPivot () const { return m_maxpivot; } 00369 00370 protected: 00371 00372 static void check_template_parameters() 00373 { 00374 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 00375 } 00376 00377 MatrixType m_qr; 00378 HCoeffsType m_hCoeffs; 00379 IntDiagSizeVectorType m_rows_transpositions; 00380 IntDiagSizeVectorType m_cols_transpositions; 00381 PermutationType m_cols_permutation; 00382 RowVectorType m_temp; 00383 bool m_isInitialized, m_usePrescribedThreshold; 00384 RealScalar m_prescribedThreshold, m_maxpivot; 00385 Index m_nonzero_pivots; 00386 RealScalar m_precision; 00387 Index m_det_pq; 00388 }; 00389 00390 template<typename MatrixType> 00391 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant () const 00392 { 00393 using std::abs; 00394 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00395 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); 00396 return abs(m_qr.diagonal().prod()); 00397 } 00398 00399 template<typename MatrixType> 00400 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDeterminant() const 00401 { 00402 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00403 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); 00404 return m_qr.diagonal().cwiseAbs().array().log().sum(); 00405 } 00406 00407 /** Performs the QR factorization of the given matrix \a matrix. The result of 00408 * the factorization is stored into \c *this, and a reference to \c *this 00409 * is returned. 00410 * 00411 * \sa class FullPivHouseholderQR, FullPivHouseholderQR(const MatrixType&) 00412 */ 00413 template<typename MatrixType> 00414 FullPivHouseholderQR<MatrixType> & FullPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix) 00415 { 00416 check_template_parameters(); 00417 00418 using std::abs; 00419 Index rows = matrix.rows(); 00420 Index cols = matrix.cols(); 00421 Index size = (std::min)(rows,cols); 00422 00423 m_qr = matrix; 00424 m_hCoeffs.resize(size); 00425 00426 m_temp.resize(cols); 00427 00428 m_precision = NumTraits<Scalar>::epsilon() * RealScalar(size); 00429 00430 m_rows_transpositions.resize(size); 00431 m_cols_transpositions.resize(size); 00432 Index number_of_transpositions = 0; 00433 00434 RealScalar biggest(0); 00435 00436 m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) 00437 m_maxpivot = RealScalar(0); 00438 00439 for (Index k = 0; k < size; ++k) 00440 { 00441 Index row_of_biggest_in_corner, col_of_biggest_in_corner; 00442 RealScalar biggest_in_corner; 00443 00444 biggest_in_corner = m_qr.bottomRightCorner(rows-k, cols-k) 00445 .cwiseAbs() 00446 .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner); 00447 row_of_biggest_in_corner += k; 00448 col_of_biggest_in_corner += k; 00449 if(k==0) biggest = biggest_in_corner; 00450 00451 // if the corner is negligible, then we have less than full rank, and we can finish early 00452 if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision)) 00453 { 00454 m_nonzero_pivots = k; 00455 for(Index i = k; i < size; i++) 00456 { 00457 m_rows_transpositions.coeffRef(i) = i; 00458 m_cols_transpositions.coeffRef(i) = i; 00459 m_hCoeffs.coeffRef(i) = Scalar(0); 00460 } 00461 break; 00462 } 00463 00464 m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner; 00465 m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner; 00466 if(k != row_of_biggest_in_corner) { 00467 m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k)); 00468 ++number_of_transpositions; 00469 } 00470 if(k != col_of_biggest_in_corner) { 00471 m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner)); 00472 ++number_of_transpositions; 00473 } 00474 00475 RealScalar beta; 00476 m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta); 00477 m_qr.coeffRef(k,k) = beta; 00478 00479 // remember the maximum absolute value of diagonal coefficients 00480 if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta); 00481 00482 m_qr.bottomRightCorner(rows-k, cols-k-1) 00483 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1)); 00484 } 00485 00486 m_cols_permutation.setIdentity(cols); 00487 for(Index k = 0; k < size; ++k) 00488 m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k)); 00489 00490 m_det_pq = (number_of_transpositions%2) ? -1 : 1; 00491 m_isInitialized = true; 00492 00493 return *this; 00494 } 00495 00496 namespace internal { 00497 00498 template<typename _MatrixType, typename Rhs> 00499 struct solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs> 00500 : solve_retval_base<FullPivHouseholderQR<_MatrixType>, Rhs> 00501 { 00502 EIGEN_MAKE_SOLVE_HELPERS(FullPivHouseholderQR<_MatrixType>,Rhs) 00503 00504 template<typename Dest> void evalTo(Dest& dst) const 00505 { 00506 const Index rows = dec().rows(), cols = dec().cols(); 00507 eigen_assert(rhs().rows() == rows); 00508 00509 // FIXME introduce nonzeroPivots() and use it here. and more generally, 00510 // make the same improvements in this dec as in FullPivLU. 00511 if(dec().rank()==0) 00512 { 00513 dst.setZero(); 00514 return; 00515 } 00516 00517 typename Rhs::PlainObject c(rhs()); 00518 00519 Matrix<Scalar,1,Rhs::ColsAtCompileTime> temp(rhs().cols()); 00520 for (Index k = 0; k < dec().rank(); ++k) 00521 { 00522 Index remainingSize = rows-k; 00523 c.row(k).swap(c.row(dec().rowsTranspositions().coeff(k))); 00524 c.bottomRightCorner(remainingSize, rhs().cols()) 00525 .applyHouseholderOnTheLeft(dec().matrixQR().col(k).tail(remainingSize-1), 00526 dec().hCoeffs().coeff(k), &temp.coeffRef(0)); 00527 } 00528 00529 dec().matrixQR() 00530 .topLeftCorner(dec().rank(), dec().rank()) 00531 .template triangularView<Upper>() 00532 .solveInPlace(c.topRows(dec().rank())); 00533 00534 for(Index i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i); 00535 for(Index i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero(); 00536 } 00537 }; 00538 00539 /** \ingroup QR_Module 00540 * 00541 * \brief Expression type for return value of FullPivHouseholderQR::matrixQ() 00542 * 00543 * \tparam MatrixType type of underlying dense matrix 00544 */ 00545 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType 00546 : public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> > 00547 { 00548 public: 00549 typedef typename MatrixType::Index Index; 00550 typedef typename FullPivHouseholderQR<MatrixType>::IntDiagSizeVectorType IntDiagSizeVectorType; 00551 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; 00552 typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1, 00553 MatrixType::MaxRowsAtCompileTime> WorkVectorType; 00554 00555 FullPivHouseholderQRMatrixQReturnType(const MatrixType& qr, 00556 const HCoeffsType& hCoeffs, 00557 const IntDiagSizeVectorType& rowsTranspositions) 00558 : m_qr(qr), 00559 m_hCoeffs(hCoeffs), 00560 m_rowsTranspositions(rowsTranspositions) 00561 {} 00562 00563 template <typename ResultType> 00564 void evalTo(ResultType& result) const 00565 { 00566 const Index rows = m_qr.rows(); 00567 WorkVectorType workspace(rows); 00568 evalTo(result, workspace); 00569 } 00570 00571 template <typename ResultType> 00572 void evalTo(ResultType& result, WorkVectorType& workspace) const 00573 { 00574 using numext::conj; 00575 // compute the product H'_0 H'_1 ... H'_n-1, 00576 // where H_k is the k-th Householder transformation I - h_k v_k v_k' 00577 // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...] 00578 const Index rows = m_qr.rows(); 00579 const Index cols = m_qr.cols(); 00580 const Index size = (std::min)(rows, cols); 00581 workspace.resize(rows); 00582 result.setIdentity(rows, rows); 00583 for (Index k = size-1; k >= 0; k--) 00584 { 00585 result.block(k, k, rows-k, rows-k) 00586 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k)); 00587 result.row(k).swap(result.row(m_rowsTranspositions.coeff(k))); 00588 } 00589 } 00590 00591 Index rows() const { return m_qr.rows(); } 00592 Index cols() const { return m_qr.rows(); } 00593 00594 protected: 00595 typename MatrixType::Nested m_qr; 00596 typename HCoeffsType::Nested m_hCoeffs; 00597 typename IntDiagSizeVectorType::Nested m_rowsTranspositions; 00598 }; 00599 00600 } // end namespace internal 00601 00602 template<typename MatrixType> 00603 inline typename FullPivHouseholderQR<MatrixType>::MatrixQReturnType FullPivHouseholderQR<MatrixType>::matrixQ() const 00604 { 00605 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 00606 return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions); 00607 } 00608 00609 /** \return the full-pivoting Householder QR decomposition of \c *this. 00610 * 00611 * \sa class FullPivHouseholderQR 00612 */ 00613 template<typename Derived> 00614 const FullPivHouseholderQR<typename MatrixBase<Derived>::PlainObject> 00615 MatrixBase<Derived>::fullPivHouseholderQr() const 00616 { 00617 return FullPivHouseholderQR<PlainObject>(eval()); 00618 } 00619 00620 } // end namespace Eigen 00621 00622 #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
Generated on Tue Jul 12 2022 17:46:53 by
