diff options
Diffstat (limited to 'extern/Eigen3/Eigen/src/QR')
| -rw-r--r-- | extern/Eigen3/Eigen/src/QR/ColPivHouseholderQR.h | 107 | ||||
| -rw-r--r-- | extern/Eigen3/Eigen/src/QR/ColPivHouseholderQR_MKL.h | 7 | ||||
| -rw-r--r-- | extern/Eigen3/Eigen/src/QR/FullPivHouseholderQR.h | 92 | ||||
| -rw-r--r-- | extern/Eigen3/Eigen/src/QR/HouseholderQR.h | 47 |
4 files changed, 182 insertions, 71 deletions
diff --git a/extern/Eigen3/Eigen/src/QR/ColPivHouseholderQR.h b/extern/Eigen3/Eigen/src/QR/ColPivHouseholderQR.h index 2daa23cc354..bec85810ccc 100644 --- a/extern/Eigen3/Eigen/src/QR/ColPivHouseholderQR.h +++ b/extern/Eigen3/Eigen/src/QR/ColPivHouseholderQR.h @@ -55,7 +55,13 @@ template<typename _MatrixType> class ColPivHouseholderQR typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType; typedef typename internal::plain_row_type<MatrixType>::type RowVectorType; typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType; - typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType; + typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType; + + private: + + typedef typename PermutationType::Index PermIndexType; + + public: /** * \brief Default Constructor. @@ -81,17 +87,29 @@ template<typename _MatrixType> class ColPivHouseholderQR ColPivHouseholderQR(Index rows, Index cols) : m_qr(rows, cols), m_hCoeffs((std::min)(rows,cols)), - m_colsPermutation(cols), + m_colsPermutation(PermIndexType(cols)), m_colsTranspositions(cols), m_temp(cols), m_colSqNorms(cols), m_isInitialized(false), m_usePrescribedThreshold(false) {} + /** \brief Constructs a QR factorization from a given matrix + * + * This constructor computes the QR factorization of the matrix \a matrix by calling + * the method compute(). It is a short cut for: + * + * \code + * ColPivHouseholderQR<MatrixType> qr(matrix.rows(), matrix.cols()); + * qr.compute(matrix); + * \endcode + * + * \sa compute() + */ ColPivHouseholderQR(const MatrixType& matrix) : m_qr(matrix.rows(), matrix.cols()), m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), - m_colsPermutation(matrix.cols()), + m_colsPermutation(PermIndexType(matrix.cols())), m_colsTranspositions(matrix.cols()), m_temp(matrix.cols()), m_colSqNorms(matrix.cols()), @@ -127,6 +145,10 @@ template<typename _MatrixType> class ColPivHouseholderQR } HouseholderSequenceType householderQ(void) const; + HouseholderSequenceType matrixQ(void) const + { + return householderQ(); + } /** \returns a reference to the matrix where the Householder QR decomposition is stored */ @@ -135,9 +157,25 @@ template<typename _MatrixType> class ColPivHouseholderQR eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); return m_qr; } - + + /** \returns a reference to the matrix where the result Householder QR is stored + * \warning The strict lower part of this matrix contains internal values. + * Only the upper triangular part should be referenced. To get it, use + * \code matrixR().template triangularView<Upper>() \endcode + * For rank-deficient matrices, use + * \code + * matrixR().topLeftCorner(rank(), rank()).template triangularView<Upper>() + * \endcode + */ + const MatrixType& matrixR() const + { + eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); + return m_qr; + } + ColPivHouseholderQR& compute(const MatrixType& matrix); + /** \returns a const reference to the column permutation matrix */ const PermutationType& colsPermutation() const { eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); @@ -181,11 +219,12 @@ template<typename _MatrixType> class ColPivHouseholderQR */ inline Index rank() const { + using std::abs; eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); - RealScalar premultiplied_threshold = internal::abs(m_maxpivot) * threshold(); + RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold(); Index result = 0; for(Index i = 0; i < m_nonzero_pivots; ++i) - result += (internal::abs(m_qr.coeff(i,i)) > premultiplied_threshold); + result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold); return result; } @@ -255,6 +294,11 @@ template<typename _MatrixType> class ColPivHouseholderQR inline Index rows() const { return m_qr.rows(); } inline Index cols() const { return m_qr.cols(); } + + /** \returns a const reference to the vector of Householder coefficients used to represent the factor \c Q. + * + * For advanced uses only. + */ const HCoeffsType& hCoeffs() const { return m_hCoeffs; } /** Allows to prescribe a threshold to be used by certain methods, such as rank(), @@ -305,7 +349,7 @@ template<typename _MatrixType> class ColPivHouseholderQR return m_usePrescribedThreshold ? m_prescribedThreshold // this formula comes from experimenting (see "LU precision tuning" thread on the list) // and turns out to be identical to Higham's formula used already in LDLt. - : NumTraits<Scalar>::epsilon() * m_qr.diagonalSize(); + : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize()); } /** \returns the number of nonzero pivots in the QR decomposition. @@ -325,6 +369,18 @@ template<typename _MatrixType> class ColPivHouseholderQR * diagonal coefficient of R. */ RealScalar maxPivot() const { return m_maxpivot; } + + /** \brief Reports whether the QR factorization was succesful. + * + * \note This function always returns \c Success. It is provided for compatibility + * with other factorization routines. + * \returns \c Success + */ + ComputationInfo info() const + { + eigen_assert(m_isInitialized && "Decomposition is not initialized."); + return Success; + } protected: MatrixType m_qr; @@ -342,9 +398,10 @@ template<typename _MatrixType> class ColPivHouseholderQR template<typename MatrixType> typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const { + using std::abs; eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); - return internal::abs(m_qr.diagonal().prod()); + return abs(m_qr.diagonal().prod()); } template<typename MatrixType> @@ -355,12 +412,22 @@ typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDetermina return m_qr.diagonal().cwiseAbs().array().log().sum(); } +/** Performs the QR factorization of the given matrix \a matrix. The result of + * the factorization is stored into \c *this, and a reference to \c *this + * is returned. + * + * \sa class ColPivHouseholderQR, ColPivHouseholderQR(const MatrixType&) + */ template<typename MatrixType> ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix) { + using std::abs; Index rows = matrix.rows(); Index cols = matrix.cols(); Index size = matrix.diagonalSize(); + + // the column permutation is stored as int indices, so just to be sure: + eigen_assert(cols<=NumTraits<int>::highest()); m_qr = matrix; m_hCoeffs.resize(size); @@ -374,7 +441,7 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const for(Index k = 0; k < cols; ++k) m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm(); - RealScalar threshold_helper = m_colSqNorms.maxCoeff() * internal::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows); + RealScalar threshold_helper = m_colSqNorms.maxCoeff() * numext::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows); m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) m_maxpivot = RealScalar(0); @@ -426,7 +493,7 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const m_qr.coeffRef(k,k) = beta; // remember the maximum absolute value of diagonal coefficients - if(internal::abs(beta) > m_maxpivot) m_maxpivot = internal::abs(beta); + if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta); // apply the householder transformation m_qr.bottomRightCorner(rows-k, cols-k-1) @@ -436,9 +503,9 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2(); } - m_colsPermutation.setIdentity(cols); - for(Index k = 0; k < m_nonzero_pivots; ++k) - m_colsPermutation.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k)); + m_colsPermutation.setIdentity(PermIndexType(cols)); + for(PermIndexType k = 0; k < m_nonzero_pivots; ++k) + m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k))); m_det_pq = (number_of_transpositions%2) ? -1 : 1; m_isInitialized = true; @@ -458,8 +525,8 @@ struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs> { eigen_assert(rhs().rows() == dec().rows()); - const int cols = dec().cols(), - nonzero_pivots = dec().nonzeroPivots(); + const Index cols = dec().cols(), + nonzero_pivots = dec().nonzeroPivots(); if(nonzero_pivots == 0) { @@ -475,19 +542,11 @@ struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs> .transpose() ); - dec().matrixQR() + dec().matrixR() .topLeftCorner(nonzero_pivots, nonzero_pivots) .template triangularView<Upper>() .solveInPlace(c.topRows(nonzero_pivots)); - - typename Rhs::PlainObject d(c); - d.topRows(nonzero_pivots) - = dec().matrixQR() - .topLeftCorner(nonzero_pivots, nonzero_pivots) - .template triangularView<Upper>() - * c.topRows(nonzero_pivots); - for(Index i = 0; i < nonzero_pivots; ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i); for(Index i = nonzero_pivots; i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero(); } diff --git a/extern/Eigen3/Eigen/src/QR/ColPivHouseholderQR_MKL.h b/extern/Eigen3/Eigen/src/QR/ColPivHouseholderQR_MKL.h index 745ecf8be98..b5b19832652 100644 --- a/extern/Eigen3/Eigen/src/QR/ColPivHouseholderQR_MKL.h +++ b/extern/Eigen3/Eigen/src/QR/ColPivHouseholderQR_MKL.h @@ -41,12 +41,13 @@ namespace Eigen { /** \internal Specialization for the data types supported by MKL */ #define EIGEN_MKL_QR_COLPIV(EIGTYPE, MKLTYPE, MKLPREFIX, EIGCOLROW, MKLCOLROW) \ -template<> inline\ +template<> inline \ ColPivHouseholderQR<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic> >& \ ColPivHouseholderQR<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic> >::compute( \ const Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>& matrix) \ \ { \ + using std::abs; \ typedef Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic> MatrixType; \ typedef MatrixType::Scalar Scalar; \ typedef MatrixType::RealScalar RealScalar; \ @@ -71,10 +72,10 @@ ColPivHouseholderQR<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynami m_isInitialized = true; \ m_maxpivot=m_qr.diagonal().cwiseAbs().maxCoeff(); \ m_hCoeffs.adjointInPlace(); \ - RealScalar premultiplied_threshold = internal::abs(m_maxpivot) * threshold(); \ + RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold(); \ lapack_int *perm = m_colsPermutation.indices().data(); \ for(i=0;i<size;i++) { \ - m_nonzero_pivots += (internal::abs(m_qr.coeff(i,i)) > premultiplied_threshold);\ + m_nonzero_pivots += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);\ } \ for(i=0;i<cols;i++) perm[i]--;\ \ diff --git a/extern/Eigen3/Eigen/src/QR/FullPivHouseholderQR.h b/extern/Eigen3/Eigen/src/QR/FullPivHouseholderQR.h index 37898e77cc2..6168e7abfb4 100644 --- a/extern/Eigen3/Eigen/src/QR/FullPivHouseholderQR.h +++ b/extern/Eigen3/Eigen/src/QR/FullPivHouseholderQR.h @@ -63,9 +63,10 @@ template<typename _MatrixType> class FullPivHouseholderQR typedef typename MatrixType::Index Index; typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType; typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; - typedef Matrix<Index, 1, ColsAtCompileTime, RowMajor, 1, MaxColsAtCompileTime> IntRowVectorType; + typedef Matrix<Index, 1, + EIGEN_SIZE_MIN_PREFER_DYNAMIC(ColsAtCompileTime,RowsAtCompileTime), RowMajor, 1, + EIGEN_SIZE_MIN_PREFER_FIXED(MaxColsAtCompileTime,MaxRowsAtCompileTime)> IntDiagSizeVectorType; typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType; - typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType; typedef typename internal::plain_row_type<MatrixType>::type RowVectorType; typedef typename internal::plain_col_type<MatrixType>::type ColVectorType; @@ -93,20 +94,32 @@ template<typename _MatrixType> class FullPivHouseholderQR FullPivHouseholderQR(Index rows, Index cols) : m_qr(rows, cols), m_hCoeffs((std::min)(rows,cols)), - m_rows_transpositions(rows), - m_cols_transpositions(cols), + m_rows_transpositions((std::min)(rows,cols)), + m_cols_transpositions((std::min)(rows,cols)), m_cols_permutation(cols), - m_temp((std::min)(rows,cols)), + m_temp(cols), m_isInitialized(false), m_usePrescribedThreshold(false) {} + /** \brief Constructs a QR factorization from a given matrix + * + * This constructor computes the QR factorization of the matrix \a matrix by calling + * the method compute(). It is a short cut for: + * + * \code + * FullPivHouseholderQR<MatrixType> qr(matrix.rows(), matrix.cols()); + * qr.compute(matrix); + * \endcode + * + * \sa compute() + */ FullPivHouseholderQR(const MatrixType& matrix) : m_qr(matrix.rows(), matrix.cols()), m_hCoeffs((std::min)(matrix.rows(), matrix.cols())), - m_rows_transpositions(matrix.rows()), - m_cols_transpositions(matrix.cols()), + m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())), + m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())), m_cols_permutation(matrix.cols()), - m_temp((std::min)(matrix.rows(), matrix.cols())), + m_temp(matrix.cols()), m_isInitialized(false), m_usePrescribedThreshold(false) { @@ -114,11 +127,12 @@ template<typename _MatrixType> class FullPivHouseholderQR } /** This method finds a solution x to the equation Ax=b, where A is the matrix of which - * *this is the QR decomposition, if any exists. + * \c *this is the QR decomposition. * * \param b the right-hand-side of the equation to solve. * - * \returns a solution. + * \returns the exact or least-square solution if the rank is greater or equal to the number of columns of A, + * and an arbitrary solution otherwise. * * \note The case where b is a matrix is not yet implemented. Also, this * code is space inefficient. @@ -152,13 +166,15 @@ template<typename _MatrixType> class FullPivHouseholderQR FullPivHouseholderQR& compute(const MatrixType& matrix); + /** \returns a const reference to the column permutation matrix */ const PermutationType& colsPermutation() const { eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); return m_cols_permutation; } - const IntColVectorType& rowsTranspositions() const + /** \returns a const reference to the vector of indices representing the rows transpositions */ + const IntDiagSizeVectorType& rowsTranspositions() const { eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); return m_rows_transpositions; @@ -201,11 +217,12 @@ template<typename _MatrixType> class FullPivHouseholderQR */ inline Index rank() const { + using std::abs; eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); - RealScalar premultiplied_threshold = internal::abs(m_maxpivot) * threshold(); + RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold(); Index result = 0; for(Index i = 0; i < m_nonzero_pivots; ++i) - result += (internal::abs(m_qr.coeff(i,i)) > premultiplied_threshold); + result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold); return result; } @@ -274,6 +291,11 @@ template<typename _MatrixType> class FullPivHouseholderQR inline Index rows() const { return m_qr.rows(); } inline Index cols() const { return m_qr.cols(); } + + /** \returns a const reference to the vector of Householder coefficients used to represent the factor \c Q. + * + * For advanced uses only. + */ const HCoeffsType& hCoeffs() const { return m_hCoeffs; } /** Allows to prescribe a threshold to be used by certain methods, such as rank(), @@ -324,7 +346,7 @@ template<typename _MatrixType> class FullPivHouseholderQR return m_usePrescribedThreshold ? m_prescribedThreshold // this formula comes from experimenting (see "LU precision tuning" thread on the list) // and turns out to be identical to Higham's formula used already in LDLt. - : NumTraits<Scalar>::epsilon() * m_qr.diagonalSize(); + : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize()); } /** \returns the number of nonzero pivots in the QR decomposition. @@ -348,8 +370,8 @@ template<typename _MatrixType> class FullPivHouseholderQR protected: MatrixType m_qr; HCoeffsType m_hCoeffs; - IntColVectorType m_rows_transpositions; - IntRowVectorType m_cols_transpositions; + IntDiagSizeVectorType m_rows_transpositions; + IntDiagSizeVectorType m_cols_transpositions; PermutationType m_cols_permutation; RowVectorType m_temp; bool m_isInitialized, m_usePrescribedThreshold; @@ -362,9 +384,10 @@ template<typename _MatrixType> class FullPivHouseholderQR template<typename MatrixType> typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const { + using std::abs; eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); - return internal::abs(m_qr.diagonal().prod()); + return abs(m_qr.diagonal().prod()); } template<typename MatrixType> @@ -375,9 +398,16 @@ typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDetermin return m_qr.diagonal().cwiseAbs().array().log().sum(); } +/** Performs the QR factorization of the given matrix \a matrix. The result of + * the factorization is stored into \c *this, and a reference to \c *this + * is returned. + * + * \sa class FullPivHouseholderQR, FullPivHouseholderQR(const MatrixType&) + */ template<typename MatrixType> FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix) { + using std::abs; Index rows = matrix.rows(); Index cols = matrix.cols(); Index size = (std::min)(rows,cols); @@ -387,10 +417,10 @@ FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons m_temp.resize(cols); - m_precision = NumTraits<Scalar>::epsilon() * size; + m_precision = NumTraits<Scalar>::epsilon() * RealScalar(size); - m_rows_transpositions.resize(matrix.rows()); - m_cols_transpositions.resize(matrix.cols()); + m_rows_transpositions.resize(size); + m_cols_transpositions.resize(size); Index number_of_transpositions = 0; RealScalar biggest(0); @@ -439,7 +469,7 @@ FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons m_qr.coeffRef(k,k) = beta; // remember the maximum absolute value of diagonal coefficients - if(internal::abs(beta) > m_maxpivot) m_maxpivot = internal::abs(beta); + if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta); m_qr.bottomRightCorner(rows-k, cols-k-1) .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1)); @@ -488,17 +518,6 @@ struct solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs> dec().hCoeffs().coeff(k), &temp.coeffRef(0)); } - if(!dec().isSurjective()) - { - // is c is in the image of R ? - RealScalar biggest_in_upper_part_of_c = c.topRows( dec().rank() ).cwiseAbs().maxCoeff(); - RealScalar biggest_in_lower_part_of_c = c.bottomRows(rows-dec().rank()).cwiseAbs().maxCoeff(); - // FIXME brain dead - const RealScalar m_precision = NumTraits<Scalar>::epsilon() * (std::min)(rows,cols); - // this internal:: prefix is needed by at least gcc 3.4 and ICC - if(!internal::isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision)) - return; - } dec().matrixQR() .topLeftCorner(dec().rank(), dec().rank()) .template triangularView<Upper>() @@ -520,14 +539,14 @@ template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType { public: typedef typename MatrixType::Index Index; - typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType; + typedef typename FullPivHouseholderQR<MatrixType>::IntDiagSizeVectorType IntDiagSizeVectorType; typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1, MatrixType::MaxRowsAtCompileTime> WorkVectorType; FullPivHouseholderQRMatrixQReturnType(const MatrixType& qr, const HCoeffsType& hCoeffs, - const IntColVectorType& rowsTranspositions) + const IntDiagSizeVectorType& rowsTranspositions) : m_qr(qr), m_hCoeffs(hCoeffs), m_rowsTranspositions(rowsTranspositions) @@ -544,6 +563,7 @@ public: template <typename ResultType> void evalTo(ResultType& result, WorkVectorType& workspace) const { + using numext::conj; // compute the product H'_0 H'_1 ... H'_n-1, // where H_k is the k-th Householder transformation I - h_k v_k v_k' // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...] @@ -555,7 +575,7 @@ public: for (Index k = size-1; k >= 0; k--) { result.block(k, k, rows-k, rows-k) - .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), internal::conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k)); + .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k)); result.row(k).swap(result.row(m_rowsTranspositions.coeff(k))); } } @@ -566,7 +586,7 @@ public: protected: typename MatrixType::Nested m_qr; typename HCoeffsType::Nested m_hCoeffs; - typename IntColVectorType::Nested m_rowsTranspositions; + typename IntDiagSizeVectorType::Nested m_rowsTranspositions; }; } // end namespace internal diff --git a/extern/Eigen3/Eigen/src/QR/HouseholderQR.h b/extern/Eigen3/Eigen/src/QR/HouseholderQR.h index 5bcb32c1e18..abc61bcbbe1 100644 --- a/extern/Eigen3/Eigen/src/QR/HouseholderQR.h +++ b/extern/Eigen3/Eigen/src/QR/HouseholderQR.h @@ -57,14 +57,14 @@ template<typename _MatrixType> class HouseholderQR typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, (MatrixType::Flags&RowMajorBit) ? RowMajor : ColMajor, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType; typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; typedef typename internal::plain_row_type<MatrixType>::type RowVectorType; - typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType; + typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType; /** - * \brief Default Constructor. - * - * The default constructor is useful in cases in which the user intends to - * perform decompositions via HouseholderQR::compute(const MatrixType&). - */ + * \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via HouseholderQR::compute(const MatrixType&). + */ HouseholderQR() : m_qr(), m_hCoeffs(), m_temp(), m_isInitialized(false) {} /** \brief Default Constructor with memory preallocation @@ -79,6 +79,18 @@ template<typename _MatrixType> class HouseholderQR m_temp(cols), m_isInitialized(false) {} + /** \brief Constructs a QR factorization from a given matrix + * + * This constructor computes the QR factorization of the matrix \a matrix by calling + * the method compute(). It is a short cut for: + * + * \code + * HouseholderQR<MatrixType> qr(matrix.rows(), matrix.cols()); + * qr.compute(matrix); + * \endcode + * + * \sa compute() + */ HouseholderQR(const MatrixType& matrix) : m_qr(matrix.rows(), matrix.cols()), m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), @@ -113,6 +125,14 @@ template<typename _MatrixType> class HouseholderQR return internal::solve_retval<HouseholderQR, Rhs>(*this, b.derived()); } + /** This method returns an expression of the unitary matrix Q as a sequence of Householder transformations. + * + * The returned expression can directly be used to perform matrix products. It can also be assigned to a dense Matrix object. + * Here is an example showing how to recover the full or thin matrix Q, as well as how to perform matrix products using operator*: + * + * Example: \include HouseholderQR_householderQ.cpp + * Output: \verbinclude HouseholderQR_householderQ.out + */ HouseholderSequenceType householderQ() const { eigen_assert(m_isInitialized && "HouseholderQR is not initialized."); @@ -161,6 +181,11 @@ template<typename _MatrixType> class HouseholderQR inline Index rows() const { return m_qr.rows(); } inline Index cols() const { return m_qr.cols(); } + + /** \returns a const reference to the vector of Householder coefficients used to represent the factor \c Q. + * + * For advanced uses only. + */ const HCoeffsType& hCoeffs() const { return m_hCoeffs; } protected: @@ -173,9 +198,10 @@ template<typename _MatrixType> class HouseholderQR template<typename MatrixType> typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const { + using std::abs; eigen_assert(m_isInitialized && "HouseholderQR is not initialized."); eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); - return internal::abs(m_qr.diagonal().prod()); + return abs(m_qr.diagonal().prod()); } template<typename MatrixType> @@ -232,7 +258,6 @@ void householder_qr_inplace_blocked(MatrixQR& mat, HCoeffs& hCoeffs, { typedef typename MatrixQR::Index Index; typedef typename MatrixQR::Scalar Scalar; - typedef typename MatrixQR::RealScalar RealScalar; typedef Block<MatrixQR,Dynamic,Dynamic> BlockType; Index rows = mat.rows(); @@ -309,6 +334,12 @@ struct solve_retval<HouseholderQR<_MatrixType>, Rhs> } // end namespace internal +/** Performs the QR factorization of the given matrix \a matrix. The result of + * the factorization is stored into \c *this, and a reference to \c *this + * is returned. + * + * \sa class HouseholderQR, HouseholderQR(const MatrixType&) + */ template<typename MatrixType> HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType& matrix) { |