diff options
Diffstat (limited to 'src/3rdparty/eigen/Eigen/src/QR')
6 files changed, 2621 insertions, 0 deletions
diff --git a/src/3rdparty/eigen/Eigen/src/QR/ColPivHouseholderQR.h b/src/3rdparty/eigen/Eigen/src/QR/ColPivHouseholderQR.h new file mode 100644 index 000000000..9b677e9bf --- /dev/null +++ b/src/3rdparty/eigen/Eigen/src/QR/ColPivHouseholderQR.h @@ -0,0 +1,674 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_COLPIVOTINGHOUSEHOLDERQR_H +#define EIGEN_COLPIVOTINGHOUSEHOLDERQR_H + +namespace Eigen { + +namespace internal { +template<typename _MatrixType> struct traits<ColPivHouseholderQR<_MatrixType> > + : traits<_MatrixType> +{ + typedef MatrixXpr XprKind; + typedef SolverStorage StorageKind; + typedef int StorageIndex; + enum { Flags = 0 }; +}; + +} // end namespace internal + +/** \ingroup QR_Module + * + * \class ColPivHouseholderQR + * + * \brief Householder rank-revealing QR decomposition of a matrix with column-pivoting + * + * \tparam _MatrixType the type of the matrix of which we are computing the QR decomposition + * + * This class performs a rank-revealing QR decomposition of a matrix \b A into matrices \b P, \b Q and \b R + * such that + * \f[ + * \mathbf{A} \, \mathbf{P} = \mathbf{Q} \, \mathbf{R} + * \f] + * by using Householder transformations. Here, \b P is a permutation matrix, \b Q a unitary matrix and \b R an + * upper triangular matrix. + * + * This decomposition performs column pivoting in order to be rank-revealing and improve + * numerical stability. It is slower than HouseholderQR, and faster than FullPivHouseholderQR. + * + * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. + * + * \sa MatrixBase::colPivHouseholderQr() + */ +template<typename _MatrixType> class ColPivHouseholderQR + : public SolverBase<ColPivHouseholderQR<_MatrixType> > +{ + public: + + typedef _MatrixType MatrixType; + typedef SolverBase<ColPivHouseholderQR> Base; + friend class SolverBase<ColPivHouseholderQR>; + + EIGEN_GENERIC_PUBLIC_INTERFACE(ColPivHouseholderQR) + enum { + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; + typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; + typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType; + 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 HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType; + typedef typename MatrixType::PlainObject PlainObject; + + private: + + typedef typename PermutationType::StorageIndex PermIndexType; + + public: + + /** + * \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via ColPivHouseholderQR::compute(const MatrixType&). + */ + ColPivHouseholderQR() + : m_qr(), + m_hCoeffs(), + m_colsPermutation(), + m_colsTranspositions(), + m_temp(), + m_colNormsUpdated(), + m_colNormsDirect(), + m_isInitialized(false), + m_usePrescribedThreshold(false) {} + + /** \brief Default Constructor with memory preallocation + * + * Like the default constructor but with preallocation of the internal data + * according to the specified problem \a size. + * \sa ColPivHouseholderQR() + */ + ColPivHouseholderQR(Index rows, Index cols) + : m_qr(rows, cols), + m_hCoeffs((std::min)(rows,cols)), + m_colsPermutation(PermIndexType(cols)), + m_colsTranspositions(cols), + m_temp(cols), + m_colNormsUpdated(cols), + m_colNormsDirect(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() + */ + template<typename InputType> + explicit ColPivHouseholderQR(const EigenBase<InputType>& matrix) + : m_qr(matrix.rows(), matrix.cols()), + m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), + m_colsPermutation(PermIndexType(matrix.cols())), + m_colsTranspositions(matrix.cols()), + m_temp(matrix.cols()), + m_colNormsUpdated(matrix.cols()), + m_colNormsDirect(matrix.cols()), + m_isInitialized(false), + m_usePrescribedThreshold(false) + { + compute(matrix.derived()); + } + + /** \brief Constructs a QR factorization from a given matrix + * + * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref. + * + * \sa ColPivHouseholderQR(const EigenBase&) + */ + template<typename InputType> + explicit ColPivHouseholderQR(EigenBase<InputType>& matrix) + : m_qr(matrix.derived()), + m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), + m_colsPermutation(PermIndexType(matrix.cols())), + m_colsTranspositions(matrix.cols()), + m_temp(matrix.cols()), + m_colNormsUpdated(matrix.cols()), + m_colNormsDirect(matrix.cols()), + m_isInitialized(false), + m_usePrescribedThreshold(false) + { + computeInPlace(); + } + + #ifdef EIGEN_PARSED_BY_DOXYGEN + /** 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. + * + * \param b the right-hand-side of the equation to solve. + * + * \returns a solution. + * + * \note_about_checking_solutions + * + * \note_about_arbitrary_choice_of_solution + * + * Example: \include ColPivHouseholderQR_solve.cpp + * Output: \verbinclude ColPivHouseholderQR_solve.out + */ + template<typename Rhs> + inline const Solve<ColPivHouseholderQR, Rhs> + solve(const MatrixBase<Rhs>& b) const; + #endif + + HouseholderSequenceType householderQ() const; + HouseholderSequenceType matrixQ() const + { + return householderQ(); + } + + /** \returns a reference to the matrix where the Householder QR decomposition is stored + */ + const MatrixType& matrixQR() const + { + 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; + } + + template<typename InputType> + ColPivHouseholderQR& compute(const EigenBase<InputType>& matrix); + + /** \returns a const reference to the column permutation matrix */ + const PermutationType& colsPermutation() const + { + eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); + return m_colsPermutation; + } + + /** \returns the absolute value of the determinant of the matrix of which + * *this is the QR decomposition. It has only linear complexity + * (that is, O(n) where n is the dimension of the square matrix) + * as the QR decomposition has already been computed. + * + * \note This is only for square matrices. + * + * \warning a determinant can be very big or small, so for matrices + * of large enough dimension, there is a risk of overflow/underflow. + * One way to work around that is to use logAbsDeterminant() instead. + * + * \sa logAbsDeterminant(), MatrixBase::determinant() + */ + typename MatrixType::RealScalar absDeterminant() const; + + /** \returns the natural log of the absolute value of the determinant of the matrix of which + * *this is the QR decomposition. It has only linear complexity + * (that is, O(n) where n is the dimension of the square matrix) + * as the QR decomposition has already been computed. + * + * \note This is only for square matrices. + * + * \note This method is useful to work around the risk of overflow/underflow that's inherent + * to determinant computation. + * + * \sa absDeterminant(), MatrixBase::determinant() + */ + typename MatrixType::RealScalar logAbsDeterminant() const; + + /** \returns the rank of the matrix of which *this is the QR decomposition. + * + * \note This method has to determine which pivots should be considered nonzero. + * For that, it uses the threshold value that you can control by calling + * setThreshold(const RealScalar&). + */ + inline Index rank() const + { + using std::abs; + eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); + RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold(); + Index result = 0; + for(Index i = 0; i < m_nonzero_pivots; ++i) + result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold); + return result; + } + + /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition. + * + * \note This method has to determine which pivots should be considered nonzero. + * For that, it uses the threshold value that you can control by calling + * setThreshold(const RealScalar&). + */ + inline Index dimensionOfKernel() const + { + eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); + return cols() - rank(); + } + + /** \returns true if the matrix of which *this is the QR decomposition represents an injective + * linear map, i.e. has trivial kernel; false otherwise. + * + * \note This method has to determine which pivots should be considered nonzero. + * For that, it uses the threshold value that you can control by calling + * setThreshold(const RealScalar&). + */ + inline bool isInjective() const + { + eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); + return rank() == cols(); + } + + /** \returns true if the matrix of which *this is the QR decomposition represents a surjective + * linear map; false otherwise. + * + * \note This method has to determine which pivots should be considered nonzero. + * For that, it uses the threshold value that you can control by calling + * setThreshold(const RealScalar&). + */ + inline bool isSurjective() const + { + eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); + return rank() == rows(); + } + + /** \returns true if the matrix of which *this is the QR decomposition is invertible. + * + * \note This method has to determine which pivots should be considered nonzero. + * For that, it uses the threshold value that you can control by calling + * setThreshold(const RealScalar&). + */ + inline bool isInvertible() const + { + eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); + return isInjective() && isSurjective(); + } + + /** \returns the inverse of the matrix of which *this is the QR decomposition. + * + * \note If this matrix is not invertible, the returned matrix has undefined coefficients. + * Use isInvertible() to first determine whether this matrix is invertible. + */ + inline const Inverse<ColPivHouseholderQR> inverse() const + { + eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); + return Inverse<ColPivHouseholderQR>(*this); + } + + 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(), + * who need to determine when pivots are to be considered nonzero. This is not used for the + * QR decomposition itself. + * + * When it needs to get the threshold value, Eigen calls threshold(). By default, this + * uses a formula to automatically determine a reasonable threshold. + * Once you have called the present method setThreshold(const RealScalar&), + * your value is used instead. + * + * \param threshold The new value to use as the threshold. + * + * A pivot will be considered nonzero if its absolute value is strictly greater than + * \f$ \vert pivot \vert \leqslant threshold \times \vert maxpivot \vert \f$ + * where maxpivot is the biggest pivot. + * + * If you want to come back to the default behavior, call setThreshold(Default_t) + */ + ColPivHouseholderQR& setThreshold(const RealScalar& threshold) + { + m_usePrescribedThreshold = true; + m_prescribedThreshold = threshold; + return *this; + } + + /** Allows to come back to the default behavior, letting Eigen use its default formula for + * determining the threshold. + * + * You should pass the special object Eigen::Default as parameter here. + * \code qr.setThreshold(Eigen::Default); \endcode + * + * See the documentation of setThreshold(const RealScalar&). + */ + ColPivHouseholderQR& setThreshold(Default_t) + { + m_usePrescribedThreshold = false; + return *this; + } + + /** Returns the threshold that will be used by certain methods such as rank(). + * + * See the documentation of setThreshold(const RealScalar&). + */ + RealScalar threshold() const + { + eigen_assert(m_isInitialized || m_usePrescribedThreshold); + 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() * RealScalar(m_qr.diagonalSize()); + } + + /** \returns the number of nonzero pivots in the QR decomposition. + * Here nonzero is meant in the exact sense, not in a fuzzy sense. + * So that notion isn't really intrinsically interesting, but it is + * still useful when implementing algorithms. + * + * \sa rank() + */ + inline Index nonzeroPivots() const + { + eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); + return m_nonzero_pivots; + } + + /** \returns the absolute value of the biggest pivot, i.e. the biggest + * diagonal coefficient of R. + */ + RealScalar maxPivot() const { return m_maxpivot; } + + /** \brief Reports whether the QR factorization was successful. + * + * \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; + } + + #ifndef EIGEN_PARSED_BY_DOXYGEN + template<typename RhsType, typename DstType> + void _solve_impl(const RhsType &rhs, DstType &dst) const; + + template<bool Conjugate, typename RhsType, typename DstType> + void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const; + #endif + + protected: + + friend class CompleteOrthogonalDecomposition<MatrixType>; + + static void check_template_parameters() + { + EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); + } + + void computeInPlace(); + + MatrixType m_qr; + HCoeffsType m_hCoeffs; + PermutationType m_colsPermutation; + IntRowVectorType m_colsTranspositions; + RowVectorType m_temp; + RealRowVectorType m_colNormsUpdated; + RealRowVectorType m_colNormsDirect; + bool m_isInitialized, m_usePrescribedThreshold; + RealScalar m_prescribedThreshold, m_maxpivot; + Index m_nonzero_pivots; + Index m_det_pq; +}; + +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 abs(m_qr.diagonal().prod()); +} + +template<typename MatrixType> +typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const +{ + 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 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> +template<typename InputType> +ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const EigenBase<InputType>& matrix) +{ + m_qr = matrix.derived(); + computeInPlace(); + return *this; +} + +template<typename MatrixType> +void ColPivHouseholderQR<MatrixType>::computeInPlace() +{ + check_template_parameters(); + + // the column permutation is stored as int indices, so just to be sure: + eigen_assert(m_qr.cols()<=NumTraits<int>::highest()); + + using std::abs; + + Index rows = m_qr.rows(); + Index cols = m_qr.cols(); + Index size = m_qr.diagonalSize(); + + m_hCoeffs.resize(size); + + m_temp.resize(cols); + + m_colsTranspositions.resize(m_qr.cols()); + Index number_of_transpositions = 0; + + m_colNormsUpdated.resize(cols); + m_colNormsDirect.resize(cols); + for (Index k = 0; k < cols; ++k) { + // colNormsDirect(k) caches the most recent directly computed norm of + // column k. + m_colNormsDirect.coeffRef(k) = m_qr.col(k).norm(); + m_colNormsUpdated.coeffRef(k) = m_colNormsDirect.coeffRef(k); + } + + RealScalar threshold_helper = numext::abs2<RealScalar>(m_colNormsUpdated.maxCoeff() * NumTraits<RealScalar>::epsilon()) / RealScalar(rows); + RealScalar norm_downdate_threshold = numext::sqrt(NumTraits<RealScalar>::epsilon()); + + m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) + m_maxpivot = RealScalar(0); + + for(Index k = 0; k < size; ++k) + { + // first, we look up in our table m_colNormsUpdated which column has the biggest norm + Index biggest_col_index; + RealScalar biggest_col_sq_norm = numext::abs2(m_colNormsUpdated.tail(cols-k).maxCoeff(&biggest_col_index)); + biggest_col_index += k; + + // Track the number of meaningful pivots but do not stop the decomposition to make + // sure that the initial matrix is properly reproduced. See bug 941. + if(m_nonzero_pivots==size && biggest_col_sq_norm < threshold_helper * RealScalar(rows-k)) + m_nonzero_pivots = k; + + // apply the transposition to the columns + m_colsTranspositions.coeffRef(k) = biggest_col_index; + if(k != biggest_col_index) { + m_qr.col(k).swap(m_qr.col(biggest_col_index)); + std::swap(m_colNormsUpdated.coeffRef(k), m_colNormsUpdated.coeffRef(biggest_col_index)); + std::swap(m_colNormsDirect.coeffRef(k), m_colNormsDirect.coeffRef(biggest_col_index)); + ++number_of_transpositions; + } + + // generate the householder vector, store it below the diagonal + RealScalar beta; + m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta); + + // apply the householder transformation to the diagonal coefficient + m_qr.coeffRef(k,k) = beta; + + // remember the maximum absolute value of diagonal coefficients + if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta); + + // apply the householder transformation + 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)); + + // update our table of norms of the columns + for (Index j = k + 1; j < cols; ++j) { + // The following implements the stable norm downgrade step discussed in + // http://www.netlib.org/lapack/lawnspdf/lawn176.pdf + // and used in LAPACK routines xGEQPF and xGEQP3. + // See lines 278-297 in http://www.netlib.org/lapack/explore-html/dc/df4/sgeqpf_8f_source.html + if (m_colNormsUpdated.coeffRef(j) != RealScalar(0)) { + RealScalar temp = abs(m_qr.coeffRef(k, j)) / m_colNormsUpdated.coeffRef(j); + temp = (RealScalar(1) + temp) * (RealScalar(1) - temp); + temp = temp < RealScalar(0) ? RealScalar(0) : temp; + RealScalar temp2 = temp * numext::abs2<RealScalar>(m_colNormsUpdated.coeffRef(j) / + m_colNormsDirect.coeffRef(j)); + if (temp2 <= norm_downdate_threshold) { + // The updated norm has become too inaccurate so re-compute the column + // norm directly. + m_colNormsDirect.coeffRef(j) = m_qr.col(j).tail(rows - k - 1).norm(); + m_colNormsUpdated.coeffRef(j) = m_colNormsDirect.coeffRef(j); + } else { + m_colNormsUpdated.coeffRef(j) *= numext::sqrt(temp); + } + } + } + } + + m_colsPermutation.setIdentity(PermIndexType(cols)); + for(PermIndexType k = 0; k < size/*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; +} + +#ifndef EIGEN_PARSED_BY_DOXYGEN +template<typename _MatrixType> +template<typename RhsType, typename DstType> +void ColPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const +{ + const Index nonzero_pivots = nonzeroPivots(); + + if(nonzero_pivots == 0) + { + dst.setZero(); + return; + } + + typename RhsType::PlainObject c(rhs); + + c.applyOnTheLeft(householderQ().setLength(nonzero_pivots).adjoint() ); + + m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots) + .template triangularView<Upper>() + .solveInPlace(c.topRows(nonzero_pivots)); + + for(Index i = 0; i < nonzero_pivots; ++i) dst.row(m_colsPermutation.indices().coeff(i)) = c.row(i); + for(Index i = nonzero_pivots; i < cols(); ++i) dst.row(m_colsPermutation.indices().coeff(i)).setZero(); +} + +template<typename _MatrixType> +template<bool Conjugate, typename RhsType, typename DstType> +void ColPivHouseholderQR<_MatrixType>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const +{ + const Index nonzero_pivots = nonzeroPivots(); + + if(nonzero_pivots == 0) + { + dst.setZero(); + return; + } + + typename RhsType::PlainObject c(m_colsPermutation.transpose()*rhs); + + m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots) + .template triangularView<Upper>() + .transpose().template conjugateIf<Conjugate>() + .solveInPlace(c.topRows(nonzero_pivots)); + + dst.topRows(nonzero_pivots) = c.topRows(nonzero_pivots); + dst.bottomRows(rows()-nonzero_pivots).setZero(); + + dst.applyOnTheLeft(householderQ().setLength(nonzero_pivots).template conjugateIf<!Conjugate>() ); +} +#endif + +namespace internal { + +template<typename DstXprType, typename MatrixType> +struct Assignment<DstXprType, Inverse<ColPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename ColPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense> +{ + typedef ColPivHouseholderQR<MatrixType> QrType; + typedef Inverse<QrType> SrcXprType; + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &) + { + dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols())); + } +}; + +} // end namespace internal + +/** \returns the matrix Q as a sequence of householder transformations. + * You can extract the meaningful part only by using: + * \code qr.householderQ().setLength(qr.nonzeroPivots()) \endcode*/ +template<typename MatrixType> +typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType> + ::householderQ() const +{ + eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); + return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()); +} + +/** \return the column-pivoting Householder QR decomposition of \c *this. + * + * \sa class ColPivHouseholderQR + */ +template<typename Derived> +const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject> +MatrixBase<Derived>::colPivHouseholderQr() const +{ + return ColPivHouseholderQR<PlainObject>(eval()); +} + +} // end namespace Eigen + +#endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H diff --git a/src/3rdparty/eigen/Eigen/src/QR/ColPivHouseholderQR_LAPACKE.h b/src/3rdparty/eigen/Eigen/src/QR/ColPivHouseholderQR_LAPACKE.h new file mode 100644 index 000000000..4e9651f83 --- /dev/null +++ b/src/3rdparty/eigen/Eigen/src/QR/ColPivHouseholderQR_LAPACKE.h @@ -0,0 +1,97 @@ +/* + Copyright (c) 2011, Intel Corporation. All rights reserved. + + Redistribution and use in source and binary forms, with or without modification, + are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + * Neither the name of Intel Corporation nor the names of its contributors may + be used to endorse or promote products derived from this software without + specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR + ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON + ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + ******************************************************************************** + * Content : Eigen bindings to LAPACKe + * Householder QR decomposition of a matrix with column pivoting based on + * LAPACKE_?geqp3 function. + ******************************************************************************** +*/ + +#ifndef EIGEN_COLPIVOTINGHOUSEHOLDERQR_LAPACKE_H +#define EIGEN_COLPIVOTINGHOUSEHOLDERQR_LAPACKE_H + +namespace Eigen { + +/** \internal Specialization for the data types supported by LAPACKe */ + +#define EIGEN_LAPACKE_QR_COLPIV(EIGTYPE, LAPACKE_TYPE, LAPACKE_PREFIX, EIGCOLROW, LAPACKE_COLROW) \ +template<> template<typename InputType> inline \ +ColPivHouseholderQR<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic> >& \ +ColPivHouseholderQR<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic> >::compute( \ + const EigenBase<InputType>& matrix) \ +\ +{ \ + using std::abs; \ + typedef Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic> MatrixType; \ + typedef MatrixType::RealScalar RealScalar; \ + Index rows = matrix.rows();\ + Index cols = matrix.cols();\ +\ + m_qr = matrix;\ + Index size = m_qr.diagonalSize();\ + m_hCoeffs.resize(size);\ +\ + m_colsTranspositions.resize(cols);\ + /*Index number_of_transpositions = 0;*/ \ +\ + m_nonzero_pivots = 0; \ + m_maxpivot = RealScalar(0);\ + m_colsPermutation.resize(cols); \ + m_colsPermutation.indices().setZero(); \ +\ + lapack_int lda = internal::convert_index<lapack_int,Index>(m_qr.outerStride()); \ + lapack_int matrix_order = LAPACKE_COLROW; \ + LAPACKE_##LAPACKE_PREFIX##geqp3( matrix_order, internal::convert_index<lapack_int,Index>(rows), internal::convert_index<lapack_int,Index>(cols), \ + (LAPACKE_TYPE*)m_qr.data(), lda, (lapack_int*)m_colsPermutation.indices().data(), (LAPACKE_TYPE*)m_hCoeffs.data()); \ + m_isInitialized = true; \ + m_maxpivot=m_qr.diagonal().cwiseAbs().maxCoeff(); \ + m_hCoeffs.adjointInPlace(); \ + RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold(); \ + lapack_int *perm = m_colsPermutation.indices().data(); \ + for(Index i=0;i<size;i++) { \ + m_nonzero_pivots += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);\ + } \ + for(Index i=0;i<cols;i++) perm[i]--;\ +\ + /*m_det_pq = (number_of_transpositions%2) ? -1 : 1; // TODO: It's not needed now; fix upon availability in Eigen */ \ +\ + return *this; \ +} + +EIGEN_LAPACKE_QR_COLPIV(double, double, d, ColMajor, LAPACK_COL_MAJOR) +EIGEN_LAPACKE_QR_COLPIV(float, float, s, ColMajor, LAPACK_COL_MAJOR) +EIGEN_LAPACKE_QR_COLPIV(dcomplex, lapack_complex_double, z, ColMajor, LAPACK_COL_MAJOR) +EIGEN_LAPACKE_QR_COLPIV(scomplex, lapack_complex_float, c, ColMajor, LAPACK_COL_MAJOR) + +EIGEN_LAPACKE_QR_COLPIV(double, double, d, RowMajor, LAPACK_ROW_MAJOR) +EIGEN_LAPACKE_QR_COLPIV(float, float, s, RowMajor, LAPACK_ROW_MAJOR) +EIGEN_LAPACKE_QR_COLPIV(dcomplex, lapack_complex_double, z, RowMajor, LAPACK_ROW_MAJOR) +EIGEN_LAPACKE_QR_COLPIV(scomplex, lapack_complex_float, c, RowMajor, LAPACK_ROW_MAJOR) + +} // end namespace Eigen + +#endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_LAPACKE_H diff --git a/src/3rdparty/eigen/Eigen/src/QR/CompleteOrthogonalDecomposition.h b/src/3rdparty/eigen/Eigen/src/QR/CompleteOrthogonalDecomposition.h new file mode 100644 index 000000000..486d3373a --- /dev/null +++ b/src/3rdparty/eigen/Eigen/src/QR/CompleteOrthogonalDecomposition.h @@ -0,0 +1,635 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2016 Rasmus Munk Larsen <rmlarsen@google.com> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H +#define EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H + +namespace Eigen { + +namespace internal { +template <typename _MatrixType> +struct traits<CompleteOrthogonalDecomposition<_MatrixType> > + : traits<_MatrixType> { + typedef MatrixXpr XprKind; + typedef SolverStorage StorageKind; + typedef int StorageIndex; + enum { Flags = 0 }; +}; + +} // end namespace internal + +/** \ingroup QR_Module + * + * \class CompleteOrthogonalDecomposition + * + * \brief Complete orthogonal decomposition (COD) of a matrix. + * + * \param MatrixType the type of the matrix of which we are computing the COD. + * + * This class performs a rank-revealing complete orthogonal decomposition of a + * matrix \b A into matrices \b P, \b Q, \b T, and \b Z such that + * \f[ + * \mathbf{A} \, \mathbf{P} = \mathbf{Q} \, + * \begin{bmatrix} \mathbf{T} & \mathbf{0} \\ + * \mathbf{0} & \mathbf{0} \end{bmatrix} \, \mathbf{Z} + * \f] + * by using Householder transformations. Here, \b P is a permutation matrix, + * \b Q and \b Z are unitary matrices and \b T an upper triangular matrix of + * size rank-by-rank. \b A may be rank deficient. + * + * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. + * + * \sa MatrixBase::completeOrthogonalDecomposition() + */ +template <typename _MatrixType> class CompleteOrthogonalDecomposition + : public SolverBase<CompleteOrthogonalDecomposition<_MatrixType> > +{ + public: + typedef _MatrixType MatrixType; + typedef SolverBase<CompleteOrthogonalDecomposition> Base; + + template<typename Derived> + friend struct internal::solve_assertion; + + EIGEN_GENERIC_PUBLIC_INTERFACE(CompleteOrthogonalDecomposition) + enum { + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; + typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; + typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> + PermutationType; + 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 HouseholderSequence< + MatrixType, typename internal::remove_all< + typename HCoeffsType::ConjugateReturnType>::type> + HouseholderSequenceType; + typedef typename MatrixType::PlainObject PlainObject; + + private: + typedef typename PermutationType::Index PermIndexType; + + public: + /** + * \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via + * \c CompleteOrthogonalDecomposition::compute(const* MatrixType&). + */ + CompleteOrthogonalDecomposition() : m_cpqr(), m_zCoeffs(), m_temp() {} + + /** \brief Default Constructor with memory preallocation + * + * Like the default constructor but with preallocation of the internal data + * according to the specified problem \a size. + * \sa CompleteOrthogonalDecomposition() + */ + CompleteOrthogonalDecomposition(Index rows, Index cols) + : m_cpqr(rows, cols), m_zCoeffs((std::min)(rows, cols)), m_temp(cols) {} + + /** \brief Constructs a complete orthogonal decomposition from a given + * matrix. + * + * This constructor computes the complete orthogonal decomposition of the + * matrix \a matrix by calling the method compute(). The default + * threshold for rank determination will be used. It is a short cut for: + * + * \code + * CompleteOrthogonalDecomposition<MatrixType> cod(matrix.rows(), + * matrix.cols()); + * cod.setThreshold(Default); + * cod.compute(matrix); + * \endcode + * + * \sa compute() + */ + template <typename InputType> + explicit CompleteOrthogonalDecomposition(const EigenBase<InputType>& matrix) + : m_cpqr(matrix.rows(), matrix.cols()), + m_zCoeffs((std::min)(matrix.rows(), matrix.cols())), + m_temp(matrix.cols()) + { + compute(matrix.derived()); + } + + /** \brief Constructs a complete orthogonal decomposition from a given matrix + * + * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref. + * + * \sa CompleteOrthogonalDecomposition(const EigenBase&) + */ + template<typename InputType> + explicit CompleteOrthogonalDecomposition(EigenBase<InputType>& matrix) + : m_cpqr(matrix.derived()), + m_zCoeffs((std::min)(matrix.rows(), matrix.cols())), + m_temp(matrix.cols()) + { + computeInPlace(); + } + + #ifdef EIGEN_PARSED_BY_DOXYGEN + /** This method computes the minimum-norm solution X to a least squares + * problem \f[\mathrm{minimize} \|A X - B\|, \f] where \b A is the matrix of + * which \c *this is the complete orthogonal decomposition. + * + * \param b the right-hand sides of the problem to solve. + * + * \returns a solution. + * + */ + template <typename Rhs> + inline const Solve<CompleteOrthogonalDecomposition, Rhs> solve( + const MatrixBase<Rhs>& b) const; + #endif + + HouseholderSequenceType householderQ(void) const; + HouseholderSequenceType matrixQ(void) const { return m_cpqr.householderQ(); } + + /** \returns the matrix \b Z. + */ + MatrixType matrixZ() const { + MatrixType Z = MatrixType::Identity(m_cpqr.cols(), m_cpqr.cols()); + applyZOnTheLeftInPlace<false>(Z); + return Z; + } + + /** \returns a reference to the matrix where the complete orthogonal + * decomposition is stored + */ + const MatrixType& matrixQTZ() const { return m_cpqr.matrixQR(); } + + /** \returns a reference to the matrix where the complete orthogonal + * decomposition is stored. + * \warning The strict lower part and \code cols() - rank() \endcode right + * columns of this matrix contains internal values. + * Only the upper triangular part should be referenced. To get it, use + * \code matrixT().template triangularView<Upper>() \endcode + * For rank-deficient matrices, use + * \code + * matrixR().topLeftCorner(rank(), rank()).template triangularView<Upper>() + * \endcode + */ + const MatrixType& matrixT() const { return m_cpqr.matrixQR(); } + + template <typename InputType> + CompleteOrthogonalDecomposition& compute(const EigenBase<InputType>& matrix) { + // Compute the column pivoted QR factorization A P = Q R. + m_cpqr.compute(matrix); + computeInPlace(); + return *this; + } + + /** \returns a const reference to the column permutation matrix */ + const PermutationType& colsPermutation() const { + return m_cpqr.colsPermutation(); + } + + /** \returns the absolute value of the determinant of the matrix of which + * *this is the complete orthogonal decomposition. It has only linear + * complexity (that is, O(n) where n is the dimension of the square matrix) + * as the complete orthogonal decomposition has already been computed. + * + * \note This is only for square matrices. + * + * \warning a determinant can be very big or small, so for matrices + * of large enough dimension, there is a risk of overflow/underflow. + * One way to work around that is to use logAbsDeterminant() instead. + * + * \sa logAbsDeterminant(), MatrixBase::determinant() + */ + typename MatrixType::RealScalar absDeterminant() const; + + /** \returns the natural log of the absolute value of the determinant of the + * matrix of which *this is the complete orthogonal decomposition. It has + * only linear complexity (that is, O(n) where n is the dimension of the + * square matrix) as the complete orthogonal decomposition has already been + * computed. + * + * \note This is only for square matrices. + * + * \note This method is useful to work around the risk of overflow/underflow + * that's inherent to determinant computation. + * + * \sa absDeterminant(), MatrixBase::determinant() + */ + typename MatrixType::RealScalar logAbsDeterminant() const; + + /** \returns the rank of the matrix of which *this is the complete orthogonal + * decomposition. + * + * \note This method has to determine which pivots should be considered + * nonzero. For that, it uses the threshold value that you can control by + * calling setThreshold(const RealScalar&). + */ + inline Index rank() const { return m_cpqr.rank(); } + + /** \returns the dimension of the kernel of the matrix of which *this is the + * complete orthogonal decomposition. + * + * \note This method has to determine which pivots should be considered + * nonzero. For that, it uses the threshold value that you can control by + * calling setThreshold(const RealScalar&). + */ + inline Index dimensionOfKernel() const { return m_cpqr.dimensionOfKernel(); } + + /** \returns true if the matrix of which *this is the decomposition represents + * an injective linear map, i.e. has trivial kernel; false otherwise. + * + * \note This method has to determine which pivots should be considered + * nonzero. For that, it uses the threshold value that you can control by + * calling setThreshold(const RealScalar&). + */ + inline bool isInjective() const { return m_cpqr.isInjective(); } + + /** \returns true if the matrix of which *this is the decomposition represents + * a surjective linear map; false otherwise. + * + * \note This method has to determine which pivots should be considered + * nonzero. For that, it uses the threshold value that you can control by + * calling setThreshold(const RealScalar&). + */ + inline bool isSurjective() const { return m_cpqr.isSurjective(); } + + /** \returns true if the matrix of which *this is the complete orthogonal + * decomposition is invertible. + * + * \note This method has to determine which pivots should be considered + * nonzero. For that, it uses the threshold value that you can control by + * calling setThreshold(const RealScalar&). + */ + inline bool isInvertible() const { return m_cpqr.isInvertible(); } + + /** \returns the pseudo-inverse of the matrix of which *this is the complete + * orthogonal decomposition. + * \warning: Do not compute \c this->pseudoInverse()*rhs to solve a linear systems. + * It is more efficient and numerically stable to call \c this->solve(rhs). + */ + inline const Inverse<CompleteOrthogonalDecomposition> pseudoInverse() const + { + eigen_assert(m_cpqr.m_isInitialized && "CompleteOrthogonalDecomposition is not initialized."); + return Inverse<CompleteOrthogonalDecomposition>(*this); + } + + inline Index rows() const { return m_cpqr.rows(); } + inline Index cols() const { return m_cpqr.cols(); } + + /** \returns a const reference to the vector of Householder coefficients used + * to represent the factor \c Q. + * + * For advanced uses only. + */ + inline const HCoeffsType& hCoeffs() const { return m_cpqr.hCoeffs(); } + + /** \returns a const reference to the vector of Householder coefficients + * used to represent the factor \c Z. + * + * For advanced uses only. + */ + const HCoeffsType& zCoeffs() const { return m_zCoeffs; } + + /** Allows to prescribe a threshold to be used by certain methods, such as + * rank(), who need to determine when pivots are to be considered nonzero. + * Most be called before calling compute(). + * + * When it needs to get the threshold value, Eigen calls threshold(). By + * default, this uses a formula to automatically determine a reasonable + * threshold. Once you have called the present method + * setThreshold(const RealScalar&), your value is used instead. + * + * \param threshold The new value to use as the threshold. + * + * A pivot will be considered nonzero if its absolute value is strictly + * greater than + * \f$ \vert pivot \vert \leqslant threshold \times \vert maxpivot \vert \f$ + * where maxpivot is the biggest pivot. + * + * If you want to come back to the default behavior, call + * setThreshold(Default_t) + */ + CompleteOrthogonalDecomposition& setThreshold(const RealScalar& threshold) { + m_cpqr.setThreshold(threshold); + return *this; + } + + /** Allows to come back to the default behavior, letting Eigen use its default + * formula for determining the threshold. + * + * You should pass the special object Eigen::Default as parameter here. + * \code qr.setThreshold(Eigen::Default); \endcode + * + * See the documentation of setThreshold(const RealScalar&). + */ + CompleteOrthogonalDecomposition& setThreshold(Default_t) { + m_cpqr.setThreshold(Default); + return *this; + } + + /** Returns the threshold that will be used by certain methods such as rank(). + * + * See the documentation of setThreshold(const RealScalar&). + */ + RealScalar threshold() const { return m_cpqr.threshold(); } + + /** \returns the number of nonzero pivots in the complete orthogonal + * decomposition. Here nonzero is meant in the exact sense, not in a + * fuzzy sense. So that notion isn't really intrinsically interesting, + * but it is still useful when implementing algorithms. + * + * \sa rank() + */ + inline Index nonzeroPivots() const { return m_cpqr.nonzeroPivots(); } + + /** \returns the absolute value of the biggest pivot, i.e. the biggest + * diagonal coefficient of R. + */ + inline RealScalar maxPivot() const { return m_cpqr.maxPivot(); } + + /** \brief Reports whether the complete orthogonal decomposition was + * successful. + * + * \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_cpqr.m_isInitialized && "Decomposition is not initialized."); + return Success; + } + +#ifndef EIGEN_PARSED_BY_DOXYGEN + template <typename RhsType, typename DstType> + void _solve_impl(const RhsType& rhs, DstType& dst) const; + + template<bool Conjugate, typename RhsType, typename DstType> + void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const; +#endif + + protected: + static void check_template_parameters() { + EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); + } + + template<bool Transpose_, typename Rhs> + void _check_solve_assertion(const Rhs& b) const { + EIGEN_ONLY_USED_FOR_DEBUG(b); + eigen_assert(m_cpqr.m_isInitialized && "CompleteOrthogonalDecomposition is not initialized."); + eigen_assert((Transpose_?derived().cols():derived().rows())==b.rows() && "CompleteOrthogonalDecomposition::solve(): invalid number of rows of the right hand side matrix b"); + } + + void computeInPlace(); + + /** Overwrites \b rhs with \f$ \mathbf{Z} * \mathbf{rhs} \f$ or + * \f$ \mathbf{\overline Z} * \mathbf{rhs} \f$ if \c Conjugate + * is set to \c true. + */ + template <bool Conjugate, typename Rhs> + void applyZOnTheLeftInPlace(Rhs& rhs) const; + + /** Overwrites \b rhs with \f$ \mathbf{Z}^* * \mathbf{rhs} \f$. + */ + template <typename Rhs> + void applyZAdjointOnTheLeftInPlace(Rhs& rhs) const; + + ColPivHouseholderQR<MatrixType> m_cpqr; + HCoeffsType m_zCoeffs; + RowVectorType m_temp; +}; + +template <typename MatrixType> +typename MatrixType::RealScalar +CompleteOrthogonalDecomposition<MatrixType>::absDeterminant() const { + return m_cpqr.absDeterminant(); +} + +template <typename MatrixType> +typename MatrixType::RealScalar +CompleteOrthogonalDecomposition<MatrixType>::logAbsDeterminant() const { + return m_cpqr.logAbsDeterminant(); +} + +/** Performs the complete orthogonal decomposition 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 CompleteOrthogonalDecomposition, + * CompleteOrthogonalDecomposition(const MatrixType&) + */ +template <typename MatrixType> +void CompleteOrthogonalDecomposition<MatrixType>::computeInPlace() +{ + check_template_parameters(); + + // the column permutation is stored as int indices, so just to be sure: + eigen_assert(m_cpqr.cols() <= NumTraits<int>::highest()); + + const Index rank = m_cpqr.rank(); + const Index cols = m_cpqr.cols(); + const Index rows = m_cpqr.rows(); + m_zCoeffs.resize((std::min)(rows, cols)); + m_temp.resize(cols); + + if (rank < cols) { + // We have reduced the (permuted) matrix to the form + // [R11 R12] + // [ 0 R22] + // where R11 is r-by-r (r = rank) upper triangular, R12 is + // r-by-(n-r), and R22 is empty or the norm of R22 is negligible. + // We now compute the complete orthogonal decomposition by applying + // Householder transformations from the right to the upper trapezoidal + // matrix X = [R11 R12] to zero out R12 and obtain the factorization + // [R11 R12] = [T11 0] * Z, where T11 is r-by-r upper triangular and + // Z = Z(0) * Z(1) ... Z(r-1) is an n-by-n orthogonal matrix. + // We store the data representing Z in R12 and m_zCoeffs. + for (Index k = rank - 1; k >= 0; --k) { + if (k != rank - 1) { + // Given the API for Householder reflectors, it is more convenient if + // we swap the leading parts of columns k and r-1 (zero-based) to form + // the matrix X_k = [X(0:k, k), X(0:k, r:n)] + m_cpqr.m_qr.col(k).head(k + 1).swap( + m_cpqr.m_qr.col(rank - 1).head(k + 1)); + } + // Construct Householder reflector Z(k) to zero out the last row of X_k, + // i.e. choose Z(k) such that + // [X(k, k), X(k, r:n)] * Z(k) = [beta, 0, .., 0]. + RealScalar beta; + m_cpqr.m_qr.row(k) + .tail(cols - rank + 1) + .makeHouseholderInPlace(m_zCoeffs(k), beta); + m_cpqr.m_qr(k, rank - 1) = beta; + if (k > 0) { + // Apply Z(k) to the first k rows of X_k + m_cpqr.m_qr.topRightCorner(k, cols - rank + 1) + .applyHouseholderOnTheRight( + m_cpqr.m_qr.row(k).tail(cols - rank).adjoint(), m_zCoeffs(k), + &m_temp(0)); + } + if (k != rank - 1) { + // Swap X(0:k,k) back to its proper location. + m_cpqr.m_qr.col(k).head(k + 1).swap( + m_cpqr.m_qr.col(rank - 1).head(k + 1)); + } + } + } +} + +template <typename MatrixType> +template <bool Conjugate, typename Rhs> +void CompleteOrthogonalDecomposition<MatrixType>::applyZOnTheLeftInPlace( + Rhs& rhs) const { + const Index cols = this->cols(); + const Index nrhs = rhs.cols(); + const Index rank = this->rank(); + Matrix<typename Rhs::Scalar, Dynamic, 1> temp((std::max)(cols, nrhs)); + for (Index k = rank-1; k >= 0; --k) { + if (k != rank - 1) { + rhs.row(k).swap(rhs.row(rank - 1)); + } + rhs.middleRows(rank - 1, cols - rank + 1) + .applyHouseholderOnTheLeft( + matrixQTZ().row(k).tail(cols - rank).transpose().template conjugateIf<!Conjugate>(), zCoeffs().template conjugateIf<Conjugate>()(k), + &temp(0)); + if (k != rank - 1) { + rhs.row(k).swap(rhs.row(rank - 1)); + } + } +} + +template <typename MatrixType> +template <typename Rhs> +void CompleteOrthogonalDecomposition<MatrixType>::applyZAdjointOnTheLeftInPlace( + Rhs& rhs) const { + const Index cols = this->cols(); + const Index nrhs = rhs.cols(); + const Index rank = this->rank(); + Matrix<typename Rhs::Scalar, Dynamic, 1> temp((std::max)(cols, nrhs)); + for (Index k = 0; k < rank; ++k) { + if (k != rank - 1) { + rhs.row(k).swap(rhs.row(rank - 1)); + } + rhs.middleRows(rank - 1, cols - rank + 1) + .applyHouseholderOnTheLeft( + matrixQTZ().row(k).tail(cols - rank).adjoint(), zCoeffs()(k), + &temp(0)); + if (k != rank - 1) { + rhs.row(k).swap(rhs.row(rank - 1)); + } + } +} + +#ifndef EIGEN_PARSED_BY_DOXYGEN +template <typename _MatrixType> +template <typename RhsType, typename DstType> +void CompleteOrthogonalDecomposition<_MatrixType>::_solve_impl( + const RhsType& rhs, DstType& dst) const { + const Index rank = this->rank(); + if (rank == 0) { + dst.setZero(); + return; + } + + // Compute c = Q^* * rhs + typename RhsType::PlainObject c(rhs); + c.applyOnTheLeft(matrixQ().setLength(rank).adjoint()); + + // Solve T z = c(1:rank, :) + dst.topRows(rank) = matrixT() + .topLeftCorner(rank, rank) + .template triangularView<Upper>() + .solve(c.topRows(rank)); + + const Index cols = this->cols(); + if (rank < cols) { + // Compute y = Z^* * [ z ] + // [ 0 ] + dst.bottomRows(cols - rank).setZero(); + applyZAdjointOnTheLeftInPlace(dst); + } + + // Undo permutation to get x = P^{-1} * y. + dst = colsPermutation() * dst; +} + +template<typename _MatrixType> +template<bool Conjugate, typename RhsType, typename DstType> +void CompleteOrthogonalDecomposition<_MatrixType>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const +{ + const Index rank = this->rank(); + + if (rank == 0) { + dst.setZero(); + return; + } + + typename RhsType::PlainObject c(colsPermutation().transpose()*rhs); + + if (rank < cols()) { + applyZOnTheLeftInPlace<!Conjugate>(c); + } + + matrixT().topLeftCorner(rank, rank) + .template triangularView<Upper>() + .transpose().template conjugateIf<Conjugate>() + .solveInPlace(c.topRows(rank)); + + dst.topRows(rank) = c.topRows(rank); + dst.bottomRows(rows()-rank).setZero(); + + dst.applyOnTheLeft(householderQ().setLength(rank).template conjugateIf<!Conjugate>() ); +} +#endif + +namespace internal { + +template<typename MatrixType> +struct traits<Inverse<CompleteOrthogonalDecomposition<MatrixType> > > + : traits<typename Transpose<typename MatrixType::PlainObject>::PlainObject> +{ + enum { Flags = 0 }; +}; + +template<typename DstXprType, typename MatrixType> +struct Assignment<DstXprType, Inverse<CompleteOrthogonalDecomposition<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename CompleteOrthogonalDecomposition<MatrixType>::Scalar>, Dense2Dense> +{ + typedef CompleteOrthogonalDecomposition<MatrixType> CodType; + typedef Inverse<CodType> SrcXprType; + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename CodType::Scalar> &) + { + typedef Matrix<typename CodType::Scalar, CodType::RowsAtCompileTime, CodType::RowsAtCompileTime, 0, CodType::MaxRowsAtCompileTime, CodType::MaxRowsAtCompileTime> IdentityMatrixType; + dst = src.nestedExpression().solve(IdentityMatrixType::Identity(src.cols(), src.cols())); + } +}; + +} // end namespace internal + +/** \returns the matrix Q as a sequence of householder transformations */ +template <typename MatrixType> +typename CompleteOrthogonalDecomposition<MatrixType>::HouseholderSequenceType +CompleteOrthogonalDecomposition<MatrixType>::householderQ() const { + return m_cpqr.householderQ(); +} + +/** \return the complete orthogonal decomposition of \c *this. + * + * \sa class CompleteOrthogonalDecomposition + */ +template <typename Derived> +const CompleteOrthogonalDecomposition<typename MatrixBase<Derived>::PlainObject> +MatrixBase<Derived>::completeOrthogonalDecomposition() const { + return CompleteOrthogonalDecomposition<PlainObject>(eval()); +} + +} // end namespace Eigen + +#endif // EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H diff --git a/src/3rdparty/eigen/Eigen/src/QR/FullPivHouseholderQR.h b/src/3rdparty/eigen/Eigen/src/QR/FullPivHouseholderQR.h new file mode 100644 index 000000000..d0664a1d8 --- /dev/null +++ b/src/3rdparty/eigen/Eigen/src/QR/FullPivHouseholderQR.h @@ -0,0 +1,713 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H +#define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H + +namespace Eigen { + +namespace internal { + +template<typename _MatrixType> struct traits<FullPivHouseholderQR<_MatrixType> > + : traits<_MatrixType> +{ + typedef MatrixXpr XprKind; + typedef SolverStorage StorageKind; + typedef int StorageIndex; + enum { Flags = 0 }; +}; + +template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType; + +template<typename MatrixType> +struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> > +{ + typedef typename MatrixType::PlainObject ReturnType; +}; + +} // end namespace internal + +/** \ingroup QR_Module + * + * \class FullPivHouseholderQR + * + * \brief Householder rank-revealing QR decomposition of a matrix with full pivoting + * + * \tparam _MatrixType the type of the matrix of which we are computing the QR decomposition + * + * This class performs a rank-revealing QR decomposition of a matrix \b A into matrices \b P, \b P', \b Q and \b R + * such that + * \f[ + * \mathbf{P} \, \mathbf{A} \, \mathbf{P}' = \mathbf{Q} \, \mathbf{R} + * \f] + * by using Householder transformations. Here, \b P and \b P' are permutation matrices, \b Q a unitary matrix + * and \b R an upper triangular matrix. + * + * This decomposition performs a very prudent full pivoting in order to be rank-revealing and achieve optimal + * numerical stability. The trade-off is that it is slower than HouseholderQR and ColPivHouseholderQR. + * + * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. + * + * \sa MatrixBase::fullPivHouseholderQr() + */ +template<typename _MatrixType> class FullPivHouseholderQR + : public SolverBase<FullPivHouseholderQR<_MatrixType> > +{ + public: + + typedef _MatrixType MatrixType; + typedef SolverBase<FullPivHouseholderQR> Base; + friend class SolverBase<FullPivHouseholderQR>; + + EIGEN_GENERIC_PUBLIC_INTERFACE(FullPivHouseholderQR) + enum { + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; + typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType; + typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; + typedef Matrix<StorageIndex, 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_row_type<MatrixType>::type RowVectorType; + typedef typename internal::plain_col_type<MatrixType>::type ColVectorType; + typedef typename MatrixType::PlainObject PlainObject; + + /** \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via FullPivHouseholderQR::compute(const MatrixType&). + */ + FullPivHouseholderQR() + : m_qr(), + m_hCoeffs(), + m_rows_transpositions(), + m_cols_transpositions(), + m_cols_permutation(), + m_temp(), + m_isInitialized(false), + m_usePrescribedThreshold(false) {} + + /** \brief Default Constructor with memory preallocation + * + * Like the default constructor but with preallocation of the internal data + * according to the specified problem \a size. + * \sa FullPivHouseholderQR() + */ + FullPivHouseholderQR(Index rows, Index cols) + : m_qr(rows, cols), + m_hCoeffs((std::min)(rows,cols)), + m_rows_transpositions((std::min)(rows,cols)), + m_cols_transpositions((std::min)(rows,cols)), + m_cols_permutation(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() + */ + template<typename InputType> + explicit FullPivHouseholderQR(const EigenBase<InputType>& matrix) + : m_qr(matrix.rows(), matrix.cols()), + m_hCoeffs((std::min)(matrix.rows(), 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(matrix.cols()), + m_isInitialized(false), + m_usePrescribedThreshold(false) + { + compute(matrix.derived()); + } + + /** \brief Constructs a QR factorization from a given matrix + * + * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref. + * + * \sa FullPivHouseholderQR(const EigenBase&) + */ + template<typename InputType> + explicit FullPivHouseholderQR(EigenBase<InputType>& matrix) + : m_qr(matrix.derived()), + m_hCoeffs((std::min)(matrix.rows(), 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(matrix.cols()), + m_isInitialized(false), + m_usePrescribedThreshold(false) + { + computeInPlace(); + } + + #ifdef EIGEN_PARSED_BY_DOXYGEN + /** This method finds a solution x to the equation Ax=b, where A is the matrix of which + * \c *this is the QR decomposition. + * + * \param b the right-hand-side of the equation to solve. + * + * \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_about_checking_solutions + * + * \note_about_arbitrary_choice_of_solution + * + * Example: \include FullPivHouseholderQR_solve.cpp + * Output: \verbinclude FullPivHouseholderQR_solve.out + */ + template<typename Rhs> + inline const Solve<FullPivHouseholderQR, Rhs> + solve(const MatrixBase<Rhs>& b) const; + #endif + + /** \returns Expression object representing the matrix Q + */ + MatrixQReturnType matrixQ(void) const; + + /** \returns a reference to the matrix where the Householder QR decomposition is stored + */ + const MatrixType& matrixQR() const + { + eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); + return m_qr; + } + + template<typename InputType> + FullPivHouseholderQR& compute(const EigenBase<InputType>& 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; + } + + /** \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; + } + + /** \returns the absolute value of the determinant of the matrix of which + * *this is the QR decomposition. It has only linear complexity + * (that is, O(n) where n is the dimension of the square matrix) + * as the QR decomposition has already been computed. + * + * \note This is only for square matrices. + * + * \warning a determinant can be very big or small, so for matrices + * of large enough dimension, there is a risk of overflow/underflow. + * One way to work around that is to use logAbsDeterminant() instead. + * + * \sa logAbsDeterminant(), MatrixBase::determinant() + */ + typename MatrixType::RealScalar absDeterminant() const; + + /** \returns the natural log of the absolute value of the determinant of the matrix of which + * *this is the QR decomposition. It has only linear complexity + * (that is, O(n) where n is the dimension of the square matrix) + * as the QR decomposition has already been computed. + * + * \note This is only for square matrices. + * + * \note This method is useful to work around the risk of overflow/underflow that's inherent + * to determinant computation. + * + * \sa absDeterminant(), MatrixBase::determinant() + */ + typename MatrixType::RealScalar logAbsDeterminant() const; + + /** \returns the rank of the matrix of which *this is the QR decomposition. + * + * \note This method has to determine which pivots should be considered nonzero. + * For that, it uses the threshold value that you can control by calling + * setThreshold(const RealScalar&). + */ + inline Index rank() const + { + using std::abs; + eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); + RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold(); + Index result = 0; + for(Index i = 0; i < m_nonzero_pivots; ++i) + result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold); + return result; + } + + /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition. + * + * \note This method has to determine which pivots should be considered nonzero. + * For that, it uses the threshold value that you can control by calling + * setThreshold(const RealScalar&). + */ + inline Index dimensionOfKernel() const + { + eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); + return cols() - rank(); + } + + /** \returns true if the matrix of which *this is the QR decomposition represents an injective + * linear map, i.e. has trivial kernel; false otherwise. + * + * \note This method has to determine which pivots should be considered nonzero. + * For that, it uses the threshold value that you can control by calling + * setThreshold(const RealScalar&). + */ + inline bool isInjective() const + { + eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); + return rank() == cols(); + } + + /** \returns true if the matrix of which *this is the QR decomposition represents a surjective + * linear map; false otherwise. + * + * \note This method has to determine which pivots should be considered nonzero. + * For that, it uses the threshold value that you can control by calling + * setThreshold(const RealScalar&). + */ + inline bool isSurjective() const + { + eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); + return rank() == rows(); + } + + /** \returns true if the matrix of which *this is the QR decomposition is invertible. + * + * \note This method has to determine which pivots should be considered nonzero. + * For that, it uses the threshold value that you can control by calling + * setThreshold(const RealScalar&). + */ + inline bool isInvertible() const + { + eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); + return isInjective() && isSurjective(); + } + + /** \returns the inverse of the matrix of which *this is the QR decomposition. + * + * \note If this matrix is not invertible, the returned matrix has undefined coefficients. + * Use isInvertible() to first determine whether this matrix is invertible. + */ + inline const Inverse<FullPivHouseholderQR> inverse() const + { + eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); + return Inverse<FullPivHouseholderQR>(*this); + } + + 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(), + * who need to determine when pivots are to be considered nonzero. This is not used for the + * QR decomposition itself. + * + * When it needs to get the threshold value, Eigen calls threshold(). By default, this + * uses a formula to automatically determine a reasonable threshold. + * Once you have called the present method setThreshold(const RealScalar&), + * your value is used instead. + * + * \param threshold The new value to use as the threshold. + * + * A pivot will be considered nonzero if its absolute value is strictly greater than + * \f$ \vert pivot \vert \leqslant threshold \times \vert maxpivot \vert \f$ + * where maxpivot is the biggest pivot. + * + * If you want to come back to the default behavior, call setThreshold(Default_t) + */ + FullPivHouseholderQR& setThreshold(const RealScalar& threshold) + { + m_usePrescribedThreshold = true; + m_prescribedThreshold = threshold; + return *this; + } + + /** Allows to come back to the default behavior, letting Eigen use its default formula for + * determining the threshold. + * + * You should pass the special object Eigen::Default as parameter here. + * \code qr.setThreshold(Eigen::Default); \endcode + * + * See the documentation of setThreshold(const RealScalar&). + */ + FullPivHouseholderQR& setThreshold(Default_t) + { + m_usePrescribedThreshold = false; + return *this; + } + + /** Returns the threshold that will be used by certain methods such as rank(). + * + * See the documentation of setThreshold(const RealScalar&). + */ + RealScalar threshold() const + { + eigen_assert(m_isInitialized || m_usePrescribedThreshold); + 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() * RealScalar(m_qr.diagonalSize()); + } + + /** \returns the number of nonzero pivots in the QR decomposition. + * Here nonzero is meant in the exact sense, not in a fuzzy sense. + * So that notion isn't really intrinsically interesting, but it is + * still useful when implementing algorithms. + * + * \sa rank() + */ + inline Index nonzeroPivots() const + { + eigen_assert(m_isInitialized && "LU is not initialized."); + return m_nonzero_pivots; + } + + /** \returns the absolute value of the biggest pivot, i.e. the biggest + * diagonal coefficient of U. + */ + RealScalar maxPivot() const { return m_maxpivot; } + + #ifndef EIGEN_PARSED_BY_DOXYGEN + template<typename RhsType, typename DstType> + void _solve_impl(const RhsType &rhs, DstType &dst) const; + + template<bool Conjugate, typename RhsType, typename DstType> + void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const; + #endif + + protected: + + static void check_template_parameters() + { + EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); + } + + void computeInPlace(); + + MatrixType m_qr; + HCoeffsType m_hCoeffs; + IntDiagSizeVectorType m_rows_transpositions; + IntDiagSizeVectorType m_cols_transpositions; + PermutationType m_cols_permutation; + RowVectorType m_temp; + bool m_isInitialized, m_usePrescribedThreshold; + RealScalar m_prescribedThreshold, m_maxpivot; + Index m_nonzero_pivots; + RealScalar m_precision; + Index m_det_pq; +}; + +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 abs(m_qr.diagonal().prod()); +} + +template<typename MatrixType> +typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDeterminant() const +{ + 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 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> +template<typename InputType> +FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(const EigenBase<InputType>& matrix) +{ + m_qr = matrix.derived(); + computeInPlace(); + return *this; +} + +template<typename MatrixType> +void FullPivHouseholderQR<MatrixType>::computeInPlace() +{ + check_template_parameters(); + + using std::abs; + Index rows = m_qr.rows(); + Index cols = m_qr.cols(); + Index size = (std::min)(rows,cols); + + + m_hCoeffs.resize(size); + + m_temp.resize(cols); + + m_precision = NumTraits<Scalar>::epsilon() * RealScalar(size); + + m_rows_transpositions.resize(size); + m_cols_transpositions.resize(size); + Index number_of_transpositions = 0; + + RealScalar biggest(0); + + m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) + m_maxpivot = RealScalar(0); + + for (Index k = 0; k < size; ++k) + { + Index row_of_biggest_in_corner, col_of_biggest_in_corner; + typedef internal::scalar_score_coeff_op<Scalar> Scoring; + typedef typename Scoring::result_type Score; + + Score score = m_qr.bottomRightCorner(rows-k, cols-k) + .unaryExpr(Scoring()) + .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner); + row_of_biggest_in_corner += k; + col_of_biggest_in_corner += k; + RealScalar biggest_in_corner = internal::abs_knowing_score<Scalar>()(m_qr(row_of_biggest_in_corner, col_of_biggest_in_corner), score); + if(k==0) biggest = biggest_in_corner; + + // if the corner is negligible, then we have less than full rank, and we can finish early + if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision)) + { + m_nonzero_pivots = k; + for(Index i = k; i < size; i++) + { + m_rows_transpositions.coeffRef(i) = internal::convert_index<StorageIndex>(i); + m_cols_transpositions.coeffRef(i) = internal::convert_index<StorageIndex>(i); + m_hCoeffs.coeffRef(i) = Scalar(0); + } + break; + } + + m_rows_transpositions.coeffRef(k) = internal::convert_index<StorageIndex>(row_of_biggest_in_corner); + m_cols_transpositions.coeffRef(k) = internal::convert_index<StorageIndex>(col_of_biggest_in_corner); + if(k != row_of_biggest_in_corner) { + m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k)); + ++number_of_transpositions; + } + if(k != col_of_biggest_in_corner) { + m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner)); + ++number_of_transpositions; + } + + RealScalar beta; + m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta); + m_qr.coeffRef(k,k) = beta; + + // remember the maximum absolute value of diagonal coefficients + 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)); + } + + m_cols_permutation.setIdentity(cols); + for(Index k = 0; k < size; ++k) + m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k)); + + m_det_pq = (number_of_transpositions%2) ? -1 : 1; + m_isInitialized = true; +} + +#ifndef EIGEN_PARSED_BY_DOXYGEN +template<typename _MatrixType> +template<typename RhsType, typename DstType> +void FullPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const +{ + const Index l_rank = rank(); + + // FIXME introduce nonzeroPivots() and use it here. and more generally, + // make the same improvements in this dec as in FullPivLU. + if(l_rank==0) + { + dst.setZero(); + return; + } + + typename RhsType::PlainObject c(rhs); + + Matrix<typename RhsType::Scalar,1,RhsType::ColsAtCompileTime> temp(rhs.cols()); + for (Index k = 0; k < l_rank; ++k) + { + Index remainingSize = rows()-k; + c.row(k).swap(c.row(m_rows_transpositions.coeff(k))); + c.bottomRightCorner(remainingSize, rhs.cols()) + .applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingSize-1), + m_hCoeffs.coeff(k), &temp.coeffRef(0)); + } + + m_qr.topLeftCorner(l_rank, l_rank) + .template triangularView<Upper>() + .solveInPlace(c.topRows(l_rank)); + + for(Index i = 0; i < l_rank; ++i) dst.row(m_cols_permutation.indices().coeff(i)) = c.row(i); + for(Index i = l_rank; i < cols(); ++i) dst.row(m_cols_permutation.indices().coeff(i)).setZero(); +} + +template<typename _MatrixType> +template<bool Conjugate, typename RhsType, typename DstType> +void FullPivHouseholderQR<_MatrixType>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const +{ + const Index l_rank = rank(); + + if(l_rank == 0) + { + dst.setZero(); + return; + } + + typename RhsType::PlainObject c(m_cols_permutation.transpose()*rhs); + + m_qr.topLeftCorner(l_rank, l_rank) + .template triangularView<Upper>() + .transpose().template conjugateIf<Conjugate>() + .solveInPlace(c.topRows(l_rank)); + + dst.topRows(l_rank) = c.topRows(l_rank); + dst.bottomRows(rows()-l_rank).setZero(); + + Matrix<Scalar, 1, DstType::ColsAtCompileTime> temp(dst.cols()); + const Index size = (std::min)(rows(), cols()); + for (Index k = size-1; k >= 0; --k) + { + Index remainingSize = rows()-k; + + dst.bottomRightCorner(remainingSize, dst.cols()) + .applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingSize-1).template conjugateIf<!Conjugate>(), + m_hCoeffs.template conjugateIf<Conjugate>().coeff(k), &temp.coeffRef(0)); + + dst.row(k).swap(dst.row(m_rows_transpositions.coeff(k))); + } +} +#endif + +namespace internal { + +template<typename DstXprType, typename MatrixType> +struct Assignment<DstXprType, Inverse<FullPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense> +{ + typedef FullPivHouseholderQR<MatrixType> QrType; + typedef Inverse<QrType> SrcXprType; + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &) + { + dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols())); + } +}; + +/** \ingroup QR_Module + * + * \brief Expression type for return value of FullPivHouseholderQR::matrixQ() + * + * \tparam MatrixType type of underlying dense matrix + */ +template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType + : public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> > +{ +public: + 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 IntDiagSizeVectorType& rowsTranspositions) + : m_qr(qr), + m_hCoeffs(hCoeffs), + m_rowsTranspositions(rowsTranspositions) + {} + + template <typename ResultType> + void evalTo(ResultType& result) const + { + const Index rows = m_qr.rows(); + WorkVectorType workspace(rows); + evalTo(result, workspace); + } + + 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), ...] + const Index rows = m_qr.rows(); + const Index cols = m_qr.cols(); + const Index size = (std::min)(rows, cols); + workspace.resize(rows); + result.setIdentity(rows, rows); + 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), conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k)); + result.row(k).swap(result.row(m_rowsTranspositions.coeff(k))); + } + } + + Index rows() const { return m_qr.rows(); } + Index cols() const { return m_qr.rows(); } + +protected: + typename MatrixType::Nested m_qr; + typename HCoeffsType::Nested m_hCoeffs; + typename IntDiagSizeVectorType::Nested m_rowsTranspositions; +}; + +// template<typename MatrixType> +// struct evaluator<FullPivHouseholderQRMatrixQReturnType<MatrixType> > +// : public evaluator<ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> > > +// {}; + +} // end namespace internal + +template<typename MatrixType> +inline typename FullPivHouseholderQR<MatrixType>::MatrixQReturnType FullPivHouseholderQR<MatrixType>::matrixQ() const +{ + eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); + return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions); +} + +/** \return the full-pivoting Householder QR decomposition of \c *this. + * + * \sa class FullPivHouseholderQR + */ +template<typename Derived> +const FullPivHouseholderQR<typename MatrixBase<Derived>::PlainObject> +MatrixBase<Derived>::fullPivHouseholderQr() const +{ + return FullPivHouseholderQR<PlainObject>(eval()); +} + +} // end namespace Eigen + +#endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H diff --git a/src/3rdparty/eigen/Eigen/src/QR/HouseholderQR.h b/src/3rdparty/eigen/Eigen/src/QR/HouseholderQR.h new file mode 100644 index 000000000..801739fbd --- /dev/null +++ b/src/3rdparty/eigen/Eigen/src/QR/HouseholderQR.h @@ -0,0 +1,434 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// Copyright (C) 2010 Vincent Lejeune +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_QR_H +#define EIGEN_QR_H + +namespace Eigen { + +namespace internal { +template<typename _MatrixType> struct traits<HouseholderQR<_MatrixType> > + : traits<_MatrixType> +{ + typedef MatrixXpr XprKind; + typedef SolverStorage StorageKind; + typedef int StorageIndex; + enum { Flags = 0 }; +}; + +} // end namespace internal + +/** \ingroup QR_Module + * + * + * \class HouseholderQR + * + * \brief Householder QR decomposition of a matrix + * + * \tparam _MatrixType the type of the matrix of which we are computing the QR decomposition + * + * This class performs a QR decomposition of a matrix \b A into matrices \b Q and \b R + * such that + * \f[ + * \mathbf{A} = \mathbf{Q} \, \mathbf{R} + * \f] + * by using Householder transformations. Here, \b Q a unitary matrix and \b R an upper triangular matrix. + * The result is stored in a compact way compatible with LAPACK. + * + * Note that no pivoting is performed. This is \b not a rank-revealing decomposition. + * If you want that feature, use FullPivHouseholderQR or ColPivHouseholderQR instead. + * + * This Householder QR decomposition is faster, but less numerically stable and less feature-full than + * FullPivHouseholderQR or ColPivHouseholderQR. + * + * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. + * + * \sa MatrixBase::householderQr() + */ +template<typename _MatrixType> class HouseholderQR + : public SolverBase<HouseholderQR<_MatrixType> > +{ + public: + + typedef _MatrixType MatrixType; + typedef SolverBase<HouseholderQR> Base; + friend class SolverBase<HouseholderQR>; + + EIGEN_GENERIC_PUBLIC_INTERFACE(HouseholderQR) + enum { + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; + 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 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&). + */ + HouseholderQR() : m_qr(), m_hCoeffs(), m_temp(), m_isInitialized(false) {} + + /** \brief Default Constructor with memory preallocation + * + * Like the default constructor but with preallocation of the internal data + * according to the specified problem \a size. + * \sa HouseholderQR() + */ + HouseholderQR(Index rows, Index cols) + : m_qr(rows, cols), + m_hCoeffs((std::min)(rows,cols)), + 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() + */ + template<typename InputType> + explicit HouseholderQR(const EigenBase<InputType>& matrix) + : m_qr(matrix.rows(), matrix.cols()), + m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), + m_temp(matrix.cols()), + m_isInitialized(false) + { + compute(matrix.derived()); + } + + + /** \brief Constructs a QR factorization from a given matrix + * + * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when + * \c MatrixType is a Eigen::Ref. + * + * \sa HouseholderQR(const EigenBase&) + */ + template<typename InputType> + explicit HouseholderQR(EigenBase<InputType>& matrix) + : m_qr(matrix.derived()), + m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), + m_temp(matrix.cols()), + m_isInitialized(false) + { + computeInPlace(); + } + + #ifdef EIGEN_PARSED_BY_DOXYGEN + /** 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. + * + * \param b the right-hand-side of the equation to solve. + * + * \returns a solution. + * + * \note_about_checking_solutions + * + * \note_about_arbitrary_choice_of_solution + * + * Example: \include HouseholderQR_solve.cpp + * Output: \verbinclude HouseholderQR_solve.out + */ + template<typename Rhs> + inline const Solve<HouseholderQR, Rhs> + solve(const MatrixBase<Rhs>& b) const; + #endif + + /** 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."); + return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()); + } + + /** \returns a reference to the matrix where the Householder QR decomposition is stored + * in a LAPACK-compatible way. + */ + const MatrixType& matrixQR() const + { + eigen_assert(m_isInitialized && "HouseholderQR is not initialized."); + return m_qr; + } + + template<typename InputType> + HouseholderQR& compute(const EigenBase<InputType>& matrix) { + m_qr = matrix.derived(); + computeInPlace(); + return *this; + } + + /** \returns the absolute value of the determinant of the matrix of which + * *this is the QR decomposition. It has only linear complexity + * (that is, O(n) where n is the dimension of the square matrix) + * as the QR decomposition has already been computed. + * + * \note This is only for square matrices. + * + * \warning a determinant can be very big or small, so for matrices + * of large enough dimension, there is a risk of overflow/underflow. + * One way to work around that is to use logAbsDeterminant() instead. + * + * \sa logAbsDeterminant(), MatrixBase::determinant() + */ + typename MatrixType::RealScalar absDeterminant() const; + + /** \returns the natural log of the absolute value of the determinant of the matrix of which + * *this is the QR decomposition. It has only linear complexity + * (that is, O(n) where n is the dimension of the square matrix) + * as the QR decomposition has already been computed. + * + * \note This is only for square matrices. + * + * \note This method is useful to work around the risk of overflow/underflow that's inherent + * to determinant computation. + * + * \sa absDeterminant(), MatrixBase::determinant() + */ + typename MatrixType::RealScalar logAbsDeterminant() const; + + 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; } + + #ifndef EIGEN_PARSED_BY_DOXYGEN + template<typename RhsType, typename DstType> + void _solve_impl(const RhsType &rhs, DstType &dst) const; + + template<bool Conjugate, typename RhsType, typename DstType> + void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const; + #endif + + protected: + + static void check_template_parameters() + { + EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); + } + + void computeInPlace(); + + MatrixType m_qr; + HCoeffsType m_hCoeffs; + RowVectorType m_temp; + bool m_isInitialized; +}; + +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 abs(m_qr.diagonal().prod()); +} + +template<typename MatrixType> +typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const +{ + 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 m_qr.diagonal().cwiseAbs().array().log().sum(); +} + +namespace internal { + +/** \internal */ +template<typename MatrixQR, typename HCoeffs> +void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0) +{ + typedef typename MatrixQR::Scalar Scalar; + typedef typename MatrixQR::RealScalar RealScalar; + Index rows = mat.rows(); + Index cols = mat.cols(); + Index size = (std::min)(rows,cols); + + eigen_assert(hCoeffs.size() == size); + + typedef Matrix<Scalar,MatrixQR::ColsAtCompileTime,1> TempType; + TempType tempVector; + if(tempData==0) + { + tempVector.resize(cols); + tempData = tempVector.data(); + } + + for(Index k = 0; k < size; ++k) + { + Index remainingRows = rows - k; + Index remainingCols = cols - k - 1; + + RealScalar beta; + mat.col(k).tail(remainingRows).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta); + mat.coeffRef(k,k) = beta; + + // apply H to remaining part of m_qr from the left + mat.bottomRightCorner(remainingRows, remainingCols) + .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), hCoeffs.coeffRef(k), tempData+k+1); + } +} + +/** \internal */ +template<typename MatrixQR, typename HCoeffs, + typename MatrixQRScalar = typename MatrixQR::Scalar, + bool InnerStrideIsOne = (MatrixQR::InnerStrideAtCompileTime == 1 && HCoeffs::InnerStrideAtCompileTime == 1)> +struct householder_qr_inplace_blocked +{ + // This is specialized for LAPACK-supported Scalar types in HouseholderQR_LAPACKE.h + static void run(MatrixQR& mat, HCoeffs& hCoeffs, Index maxBlockSize=32, + typename MatrixQR::Scalar* tempData = 0) + { + typedef typename MatrixQR::Scalar Scalar; + typedef Block<MatrixQR,Dynamic,Dynamic> BlockType; + + Index rows = mat.rows(); + Index cols = mat.cols(); + Index size = (std::min)(rows, cols); + + typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType; + TempType tempVector; + if(tempData==0) + { + tempVector.resize(cols); + tempData = tempVector.data(); + } + + Index blockSize = (std::min)(maxBlockSize,size); + + Index k = 0; + for (k = 0; k < size; k += blockSize) + { + Index bs = (std::min)(size-k,blockSize); // actual size of the block + Index tcols = cols - k - bs; // trailing columns + Index brows = rows-k; // rows of the block + + // partition the matrix: + // A00 | A01 | A02 + // mat = A10 | A11 | A12 + // A20 | A21 | A22 + // and performs the qr dec of [A11^T A12^T]^T + // and update [A21^T A22^T]^T using level 3 operations. + // Finally, the algorithm continue on A22 + + BlockType A11_21 = mat.block(k,k,brows,bs); + Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs); + + householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData); + + if(tcols) + { + BlockType A21_22 = mat.block(k,k+bs,brows,tcols); + apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment, false); // false == backward + } + } + } +}; + +} // end namespace internal + +#ifndef EIGEN_PARSED_BY_DOXYGEN +template<typename _MatrixType> +template<typename RhsType, typename DstType> +void HouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const +{ + const Index rank = (std::min)(rows(), cols()); + + typename RhsType::PlainObject c(rhs); + + c.applyOnTheLeft(householderQ().setLength(rank).adjoint() ); + + m_qr.topLeftCorner(rank, rank) + .template triangularView<Upper>() + .solveInPlace(c.topRows(rank)); + + dst.topRows(rank) = c.topRows(rank); + dst.bottomRows(cols()-rank).setZero(); +} + +template<typename _MatrixType> +template<bool Conjugate, typename RhsType, typename DstType> +void HouseholderQR<_MatrixType>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const +{ + const Index rank = (std::min)(rows(), cols()); + + typename RhsType::PlainObject c(rhs); + + m_qr.topLeftCorner(rank, rank) + .template triangularView<Upper>() + .transpose().template conjugateIf<Conjugate>() + .solveInPlace(c.topRows(rank)); + + dst.topRows(rank) = c.topRows(rank); + dst.bottomRows(rows()-rank).setZero(); + + dst.applyOnTheLeft(householderQ().setLength(rank).template conjugateIf<!Conjugate>() ); +} +#endif + +/** 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> +void HouseholderQR<MatrixType>::computeInPlace() +{ + check_template_parameters(); + + Index rows = m_qr.rows(); + Index cols = m_qr.cols(); + Index size = (std::min)(rows,cols); + + m_hCoeffs.resize(size); + + m_temp.resize(cols); + + internal::householder_qr_inplace_blocked<MatrixType, HCoeffsType>::run(m_qr, m_hCoeffs, 48, m_temp.data()); + + m_isInitialized = true; +} + +/** \return the Householder QR decomposition of \c *this. + * + * \sa class HouseholderQR + */ +template<typename Derived> +const HouseholderQR<typename MatrixBase<Derived>::PlainObject> +MatrixBase<Derived>::householderQr() const +{ + return HouseholderQR<PlainObject>(eval()); +} + +} // end namespace Eigen + +#endif // EIGEN_QR_H diff --git a/src/3rdparty/eigen/Eigen/src/QR/HouseholderQR_LAPACKE.h b/src/3rdparty/eigen/Eigen/src/QR/HouseholderQR_LAPACKE.h new file mode 100644 index 000000000..1dc7d5363 --- /dev/null +++ b/src/3rdparty/eigen/Eigen/src/QR/HouseholderQR_LAPACKE.h @@ -0,0 +1,68 @@ +/* + Copyright (c) 2011, Intel Corporation. All rights reserved. + + Redistribution and use in source and binary forms, with or without modification, + are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + * Neither the name of Intel Corporation nor the names of its contributors may + be used to endorse or promote products derived from this software without + specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR + ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON + ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + ******************************************************************************** + * Content : Eigen bindings to LAPACKe + * Householder QR decomposition of a matrix w/o pivoting based on + * LAPACKE_?geqrf function. + ******************************************************************************** +*/ + +#ifndef EIGEN_QR_LAPACKE_H +#define EIGEN_QR_LAPACKE_H + +namespace Eigen { + +namespace internal { + +/** \internal Specialization for the data types supported by LAPACKe */ + +#define EIGEN_LAPACKE_QR_NOPIV(EIGTYPE, LAPACKE_TYPE, LAPACKE_PREFIX) \ +template<typename MatrixQR, typename HCoeffs> \ +struct householder_qr_inplace_blocked<MatrixQR, HCoeffs, EIGTYPE, true> \ +{ \ + static void run(MatrixQR& mat, HCoeffs& hCoeffs, Index = 32, \ + typename MatrixQR::Scalar* = 0) \ + { \ + lapack_int m = (lapack_int) mat.rows(); \ + lapack_int n = (lapack_int) mat.cols(); \ + lapack_int lda = (lapack_int) mat.outerStride(); \ + lapack_int matrix_order = (MatrixQR::IsRowMajor) ? LAPACK_ROW_MAJOR : LAPACK_COL_MAJOR; \ + LAPACKE_##LAPACKE_PREFIX##geqrf( matrix_order, m, n, (LAPACKE_TYPE*)mat.data(), lda, (LAPACKE_TYPE*)hCoeffs.data()); \ + hCoeffs.adjointInPlace(); \ + } \ +}; + +EIGEN_LAPACKE_QR_NOPIV(double, double, d) +EIGEN_LAPACKE_QR_NOPIV(float, float, s) +EIGEN_LAPACKE_QR_NOPIV(dcomplex, lapack_complex_double, z) +EIGEN_LAPACKE_QR_NOPIV(scomplex, lapack_complex_float, c) + +} // end namespace internal + +} // end namespace Eigen + +#endif // EIGEN_QR_LAPACKE_H |