init repo.
This commit is contained in:
tqcq
674
3party/eigen/Eigen/src/QR/ColPivHouseholderQR.h
Normal file
674
3party/eigen/Eigen/src/QR/ColPivHouseholderQR.h
Normal file
@@ -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
|
||||
97
3party/eigen/Eigen/src/QR/ColPivHouseholderQR_LAPACKE.h
Normal file
97
3party/eigen/Eigen/src/QR/ColPivHouseholderQR_LAPACKE.h
Normal file
@@ -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
|
||||
635
3party/eigen/Eigen/src/QR/CompleteOrthogonalDecomposition.h
Normal file
635
3party/eigen/Eigen/src/QR/CompleteOrthogonalDecomposition.h
Normal file
@@ -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
|
||||
713
3party/eigen/Eigen/src/QR/FullPivHouseholderQR.h
Normal file
713
3party/eigen/Eigen/src/QR/FullPivHouseholderQR.h
Normal file
@@ -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
|
||||
434
3party/eigen/Eigen/src/QR/HouseholderQR.h
Normal file
434
3party/eigen/Eigen/src/QR/HouseholderQR.h
Normal file
@@ -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
|
||||
68
3party/eigen/Eigen/src/QR/HouseholderQR_LAPACKE.h
Normal file
68
3party/eigen/Eigen/src/QR/HouseholderQR_LAPACKE.h
Normal file
@@ -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
|
||||
Reference in New Issue
Block a user