GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-04-22 01:29:35 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix QR, again

Browse Source
This commit is contained in:
Charles Schlosser 2023-01-13 03:23:17 +00:00 committed by Rasmus Munk Larsen
parent 4d05765345
commit 68082b8226
10 changed files with 341 additions and 295 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
Eigen/src/Cholesky/LLT_LAPACKE.h
View File

@ -70,6 +70,7 @@ namespace lapacke_helpers {
template<typename Scalar, UpLoType Mode> template<typename Scalar, UpLoType Mode>
struct lapacke_llt { struct lapacke_llt {
EIGEN_STATIC_ASSERT(((Mode == Lower) || (Mode == Upper)),MODE_MUST_BE_UPPER_OR_LOWER)
template<typename MatrixType> template<typename MatrixType>
static Index blocked(MatrixType& m) static Index blocked(MatrixType& m)
{ {
@ -80,10 +81,11 @@ namespace lapacke_helpers {
/* Set up parameters for ?potrf */ /* Set up parameters for ?potrf */
lapack_int size = to_lapack(m.rows()); lapack_int size = to_lapack(m.rows());
lapack_int matrix_order = lapack_storage_of(m); lapack_int matrix_order = lapack_storage_of(m);
constexpr char uplo = Mode == Upper ? 'U' : 'L';
Scalar* a = &(m.coeffRef(0,0)); Scalar* a = &(m.coeffRef(0,0));
lapack_int lda = to_lapack(m.outerStride()); lapack_int lda = to_lapack(m.outerStride());
lapack_int info = potrf(matrix_order, translate_mode<Mode>, size, to_lapack(a), lda ); lapack_int info = potrf(matrix_order, uplo, size, to_lapack(a), lda );
info = (info==0) ? -1 : info>0 ? info-1 : size; info = (info==0) ? -1 : info>0 ? info-1 : size;
return info; return info;
} }

6
Eigen/src/Core/MatrixBase.h
View File

@ -362,9 +362,9 @@ template<typename Derived> class MatrixBase
/////////// QR module /////////// /////////// QR module ///////////
inline const HouseholderQR<PlainObject> householderQr() const; inline const HouseholderQR<PlainObject> householderQr() const;
inline const ColPivHouseholderQR<PlainObject> colPivHouseholderQr() const; template<typename PermutationIndex = DefaultPermutationIndex> inline const ColPivHouseholderQR<PlainObject, PermutationIndex> colPivHouseholderQr() const;
inline const FullPivHouseholderQR<PlainObject> fullPivHouseholderQr() const; template<typename PermutationIndex = DefaultPermutationIndex> inline const FullPivHouseholderQR<PlainObject, PermutationIndex> fullPivHouseholderQr() const;
inline const CompleteOrthogonalDecomposition<PlainObject> completeOrthogonalDecomposition() const; template<typename PermutationIndex = DefaultPermutationIndex> inline const CompleteOrthogonalDecomposition<PlainObject, PermutationIndex> completeOrthogonalDecomposition() const;
/////////// Eigenvalues module /////////// /////////// Eigenvalues module ///////////

13
Eigen/src/Core/util/ForwardDeclarations.h
View File

@ -252,15 +252,22 @@ template<typename ExpressionType, int Direction> class VectorwiseOp;
template<typename MatrixType,int RowFactor,int ColFactor> class Replicate; template<typename MatrixType,int RowFactor,int ColFactor> class Replicate;
template<typename MatrixType, int Direction = BothDirections> class Reverse; template<typename MatrixType, int Direction = BothDirections> class Reverse;
#if defined(EIGEN_USE_LAPACKE)
// Lapacke interface requires StorageIndex to be lapack_int
typedef lapack_int DefaultPermutationIndex;
#else
typedef int DefaultPermutationIndex;
#endif
template<typename MatrixType> class FullPivLU; template<typename MatrixType> class FullPivLU;
template<typename MatrixType> class PartialPivLU; template<typename MatrixType> class PartialPivLU;
namespace internal { namespace internal {
template<typename MatrixType> struct inverse_impl; template<typename MatrixType> struct inverse_impl;
} }
template<typename MatrixType> class HouseholderQR; template<typename MatrixType> class HouseholderQR;
template<typename MatrixType> class ColPivHouseholderQR; template<typename MatrixType, typename PermutationIndex = DefaultPermutationIndex> class ColPivHouseholderQR;
template<typename MatrixType> class FullPivHouseholderQR; template<typename MatrixType, typename PermutationIndex = DefaultPermutationIndex> class FullPivHouseholderQR;
template<typename MatrixType> class CompleteOrthogonalDecomposition; template<typename MatrixType, typename PermutationIndex = DefaultPermutationIndex> class CompleteOrthogonalDecomposition;
template<typename MatrixType> class SVDBase; template<typename MatrixType> class SVDBase;
template<typename MatrixType, int Options = 0> class JacobiSVD; template<typename MatrixType, int Options = 0> class JacobiSVD;
template<typename MatrixType, int Options = 0> class BDCSVD; template<typename MatrixType, int Options = 0> class BDCSVD;

128
Eigen/src/QR/ColPivHouseholderQR.h
View File

@ -16,12 +16,12 @@
namespace Eigen { namespace Eigen {
namespace internal { namespace internal {
template<typename MatrixType_> struct traits<ColPivHouseholderQR<MatrixType_> > template<typename MatrixType_, typename PermutationIndex_> struct traits<ColPivHouseholderQR<MatrixType_, PermutationIndex_> >
: traits<MatrixType_> : traits<MatrixType_>
{ {
typedef MatrixXpr XprKind; typedef MatrixXpr XprKind;
typedef SolverStorage StorageKind; typedef SolverStorage StorageKind;
typedef int StorageIndex; typedef PermutationIndex_ PermutationIndex;
enum { Flags = 0 }; enum { Flags = 0 };
}; };
@ -50,31 +50,41 @@ template<typename MatrixType_> struct traits<ColPivHouseholderQR<MatrixType_> >
* *
* \sa MatrixBase::colPivHouseholderQr() * \sa MatrixBase::colPivHouseholderQr()
*/ */
template<typename MatrixType_> class ColPivHouseholderQR template<typename MatrixType_, typename PermutationIndex_> class ColPivHouseholderQR
: public SolverBase<ColPivHouseholderQR<MatrixType_> > : public SolverBase<ColPivHouseholderQR<MatrixType_, PermutationIndex_> >
{ {
public: public:
typedef MatrixType_ MatrixType; typedef MatrixType_ MatrixType;
typedef SolverBase<ColPivHouseholderQR> Base; typedef SolverBase<ColPivHouseholderQR> Base;
friend class SolverBase<ColPivHouseholderQR>; friend class SolverBase<ColPivHouseholderQR>;
typedef PermutationIndex_ PermutationIndex;
EIGEN_GENERIC_PUBLIC_INTERFACE(ColPivHouseholderQR) EIGEN_GENERIC_PUBLIC_INTERFACE(ColPivHouseholderQR)
enum { enum {
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
}; };
typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType; typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime, PermutationIndex> PermutationType;
typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType; typedef typename internal::plain_row_type<MatrixType, PermutationIndex>::type IntRowVectorType;
typedef typename internal::plain_row_type<MatrixType>::type RowVectorType; typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType; typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType;
typedef HouseholderSequence<MatrixType,internal::remove_all_t<typename HCoeffsType::ConjugateReturnType>> HouseholderSequenceType; typedef HouseholderSequence<MatrixType,internal::remove_all_t<typename HCoeffsType::ConjugateReturnType>> HouseholderSequenceType;
typedef typename MatrixType::PlainObject PlainObject; typedef typename MatrixType::PlainObject PlainObject;
private: private:
void init(Index rows, Index cols) {
typedef typename PermutationType::StorageIndex PermIndexType; Index diag = numext::mini(rows, cols);
m_hCoeffs = HCoeffsType(diag);
m_colsPermutation = PermutationType(cols);
m_colsTranspositions = IntRowVectorType(cols);
m_temp = RealRowVectorType(cols);
m_colNormsUpdated = RealRowVectorType(cols);
m_colNormsDirect = RealRowVectorType(cols);
m_isInitialized = false;
m_usePrescribedThreshold = false;
}
public: public:
@ -101,16 +111,7 @@ template<typename MatrixType_> class ColPivHouseholderQR
* according to the specified problem \a size. * according to the specified problem \a size.
* \sa ColPivHouseholderQR() * \sa ColPivHouseholderQR()
*/ */
ColPivHouseholderQR(Index rows, Index cols) ColPivHouseholderQR(Index rows, Index cols) : m_qr(rows, cols) { init(rows, 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 /** \brief Constructs a QR factorization from a given matrix
* *
@ -124,18 +125,9 @@ template<typename MatrixType_> class ColPivHouseholderQR
* *
* \sa compute() * \sa compute()
*/ */
template<typename InputType> template <typename InputType>
explicit ColPivHouseholderQR(const EigenBase<InputType>& matrix) explicit ColPivHouseholderQR(const EigenBase<InputType>& matrix) : m_qr(matrix.rows(), matrix.cols()) {
: m_qr(matrix.rows(), matrix.cols()), init(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()); compute(matrix.derived());
} }
@ -145,18 +137,9 @@ template<typename MatrixType_> class ColPivHouseholderQR
* *
* \sa ColPivHouseholderQR(const EigenBase&) * \sa ColPivHouseholderQR(const EigenBase&)
*/ */
template<typename InputType> template <typename InputType>
explicit ColPivHouseholderQR(EigenBase<InputType>& matrix) explicit ColPivHouseholderQR(EigenBase<InputType>& matrix) : m_qr(matrix.derived()) {
: m_qr(matrix.derived()), init(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)
{
computeInPlace(); computeInPlace();
} }
@ -441,7 +424,7 @@ template<typename MatrixType_> class ColPivHouseholderQR
protected: protected:
friend class CompleteOrthogonalDecomposition<MatrixType>; friend class CompleteOrthogonalDecomposition<MatrixType, PermutationIndex>;
EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar) EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
@ -460,8 +443,8 @@ template<typename MatrixType_> class ColPivHouseholderQR
Index m_det_p; Index m_det_p;
}; };
template<typename MatrixType> template<typename MatrixType, typename PermutationIndex>
typename MatrixType::Scalar ColPivHouseholderQR<MatrixType>::determinant() const typename MatrixType::Scalar ColPivHouseholderQR<MatrixType, PermutationIndex>::determinant() const
{ {
eigen_assert(m_isInitialized && "HouseholderQR is not initialized."); 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!"); eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
@ -470,8 +453,8 @@ typename MatrixType::Scalar ColPivHouseholderQR<MatrixType>::determinant() const
return m_qr.diagonal().prod() * detQ * Scalar(m_det_p); return m_qr.diagonal().prod() * detQ * Scalar(m_det_p);
} }
template<typename MatrixType> template<typename MatrixType, typename PermutationIndex>
typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType, PermutationIndex>::absDeterminant() const
{ {
using std::abs; using std::abs;
eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
@ -479,8 +462,8 @@ typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant(
return abs(m_qr.diagonal().prod()); return abs(m_qr.diagonal().prod());
} }
template<typename MatrixType> template<typename MatrixType, typename PermutationIndex>
typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType, PermutationIndex>::logAbsDeterminant() const
{ {
eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 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!"); eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
@ -493,20 +476,20 @@ typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDetermina
* *
* \sa class ColPivHouseholderQR, ColPivHouseholderQR(const MatrixType&) * \sa class ColPivHouseholderQR, ColPivHouseholderQR(const MatrixType&)
*/ */
template<typename MatrixType> template<typename MatrixType, typename PermutationIndex>
template<typename InputType> template<typename InputType>
ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const EigenBase<InputType>& matrix) ColPivHouseholderQR<MatrixType, PermutationIndex>& ColPivHouseholderQR<MatrixType, PermutationIndex>::compute(const EigenBase<InputType>& matrix)
{ {
m_qr = matrix.derived(); m_qr = matrix.derived();
computeInPlace(); computeInPlace();
return *this; return *this;
} }
template<typename MatrixType> template<typename MatrixType, typename PermutationIndex>
void ColPivHouseholderQR<MatrixType>::computeInPlace() void ColPivHouseholderQR<MatrixType, PermutationIndex>::computeInPlace()
{ {
// the column permutation is stored as int indices, so just to be sure:
eigen_assert(m_qr.cols()<=NumTraits<int>::highest()); eigen_assert(m_qr.cols()<=NumTraits<PermutationIndex>::highest());
using std::abs; using std::abs;
@ -549,7 +532,7 @@ void ColPivHouseholderQR<MatrixType>::computeInPlace()
m_nonzero_pivots = k; m_nonzero_pivots = k;
// apply the transposition to the columns // apply the transposition to the columns
m_colsTranspositions.coeffRef(k) = biggest_col_index; m_colsTranspositions.coeffRef(k) = static_cast<PermutationIndex>(biggest_col_index);
if(k != biggest_col_index) { if(k != biggest_col_index) {
m_qr.col(k).swap(m_qr.col(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_colNormsUpdated.coeffRef(k), m_colNormsUpdated.coeffRef(biggest_col_index));
@ -595,18 +578,18 @@ void ColPivHouseholderQR<MatrixType>::computeInPlace()
} }
} }
m_colsPermutation.setIdentity(PermIndexType(cols)); m_colsPermutation.setIdentity(cols);
for(PermIndexType k = 0; k < size/*m_nonzero_pivots*/; ++k) for(Index k = 0; k < size/*m_nonzero_pivots*/; ++k)
m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k))); m_colsPermutation.applyTranspositionOnTheRight(k, static_cast<Index>(m_colsTranspositions.coeff(k)));
m_det_p = (number_of_transpositions%2) ? -1 : 1; m_det_p = (number_of_transpositions%2) ? -1 : 1;
m_isInitialized = true; m_isInitialized = true;
} }
#ifndef EIGEN_PARSED_BY_DOXYGEN #ifndef EIGEN_PARSED_BY_DOXYGEN
template<typename MatrixType_> template<typename MatrixType_, typename PermutationIndex_>
template<typename RhsType, typename DstType> template<typename RhsType, typename DstType>
void ColPivHouseholderQR<MatrixType_>::_solve_impl(const RhsType &rhs, DstType &dst) const void ColPivHouseholderQR<MatrixType_, PermutationIndex_>::_solve_impl(const RhsType &rhs, DstType &dst) const
{ {
const Index nonzero_pivots = nonzeroPivots(); const Index nonzero_pivots = nonzeroPivots();
@ -628,9 +611,9 @@ void ColPivHouseholderQR<MatrixType_>::_solve_impl(const RhsType &rhs, DstType &
for(Index i = nonzero_pivots; i < cols(); ++i) dst.row(m_colsPermutation.indices().coeff(i)).setZero(); for(Index i = nonzero_pivots; i < cols(); ++i) dst.row(m_colsPermutation.indices().coeff(i)).setZero();
} }
template<typename MatrixType_> template<typename MatrixType_, typename PermutationIndex_>
template<bool Conjugate, typename RhsType, typename DstType> template<bool Conjugate, typename RhsType, typename DstType>
void ColPivHouseholderQR<MatrixType_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const void ColPivHouseholderQR<MatrixType_, PermutationIndex_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
{ {
const Index nonzero_pivots = nonzeroPivots(); const Index nonzero_pivots = nonzeroPivots();
@ -656,10 +639,10 @@ void ColPivHouseholderQR<MatrixType_>::_solve_impl_transposed(const RhsType &rhs
namespace internal { namespace internal {
template<typename DstXprType, typename MatrixType> template<typename DstXprType, typename MatrixType, typename PermutationIndex>
struct Assignment<DstXprType, Inverse<ColPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename ColPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense> struct Assignment<DstXprType, Inverse<ColPivHouseholderQR<MatrixType, PermutationIndex> >, internal::assign_op<typename DstXprType::Scalar,typename ColPivHouseholderQR<MatrixType, PermutationIndex>::Scalar>, Dense2Dense>
{ {
typedef ColPivHouseholderQR<MatrixType> QrType; typedef ColPivHouseholderQR<MatrixType, PermutationIndex> QrType;
typedef Inverse<QrType> SrcXprType; typedef Inverse<QrType> SrcXprType;
static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &) static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &)
{ {
@ -672,8 +655,8 @@ struct Assignment<DstXprType, Inverse<ColPivHouseholderQR<MatrixType> >, interna
/** \returns the matrix Q as a sequence of householder transformations. /** \returns the matrix Q as a sequence of householder transformations.
* You can extract the meaningful part only by using: * You can extract the meaningful part only by using:
* \code qr.householderQ().setLength(qr.nonzeroPivots()) \endcode*/ * \code qr.householderQ().setLength(qr.nonzeroPivots()) \endcode*/
template<typename MatrixType> template<typename MatrixType, typename PermutationIndex>
typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType> typename ColPivHouseholderQR<MatrixType, PermutationIndex>::HouseholderSequenceType ColPivHouseholderQR<MatrixType, PermutationIndex>
::householderQ() const ::householderQ() const
{ {
eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
@ -685,10 +668,11 @@ typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHousehol
* \sa class ColPivHouseholderQR * \sa class ColPivHouseholderQR
*/ */
template<typename Derived> template<typename Derived>
const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject> template<typename PermutationIndexType>
const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject, PermutationIndexType>
MatrixBase<Derived>::colPivHouseholderQr() const MatrixBase<Derived>::colPivHouseholderQr() const
{ {
return ColPivHouseholderQR<PlainObject>(eval()); return ColPivHouseholderQR<PlainObject, PermutationIndexType>(eval());
} }
} // end namespace Eigen } // end namespace Eigen

161
Eigen/src/QR/ColPivHouseholderQR_LAPACKE.h
View File

@ -38,62 +38,115 @@
namespace Eigen { namespace Eigen {
/** \internal Specialization for the data types supported by LAPACKe */ #if defined(EIGEN_USE_LAPACKE)
#define EIGEN_LAPACKE_QR_COLPIV(EIGTYPE, LAPACKE_TYPE, LAPACKE_PREFIX, EIGCOLROW, LAPACKE_COLROW) \ template<typename Scalar>
template<> template<typename InputType> inline \ inline lapack_int call_geqp3(int matrix_layout, lapack_int m, lapack_int n, Scalar* a, lapack_int lda, lapack_int* jpvt, Scalar* tau);
ColPivHouseholderQR<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic> >& \ template<>
ColPivHouseholderQR<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic> >::compute( \ inline lapack_int call_geqp3(int matrix_layout, lapack_int m, lapack_int n, float* a, lapack_int lda, lapack_int* jpvt, float* tau)
const EigenBase<InputType>& matrix) \ { return LAPACKE_sgeqp3(matrix_layout, m, n, a, lda, jpvt, tau); }
\ template<>
{ \ inline lapack_int call_geqp3(int matrix_layout, lapack_int m, lapack_int n, double* a, lapack_int lda, lapack_int* jpvt, double* tau)
using std::abs; \ { return LAPACKE_dgeqp3(matrix_layout, m, n, a, lda, jpvt, tau); }
typedef Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic> MatrixType; \ template<>
typedef MatrixType::RealScalar RealScalar; \ inline lapack_int call_geqp3(int matrix_layout, lapack_int m, lapack_int n, lapack_complex_float* a, lapack_int lda, lapack_int* jpvt, lapack_complex_float* tau)
Index rows = matrix.rows();\ { return LAPACKE_cgeqp3(matrix_layout, m, n, a, lda, jpvt, tau); }
Index cols = matrix.cols();\ template<>
\ inline lapack_int call_geqp3(int matrix_layout, lapack_int m, lapack_int n, lapack_complex_double* a, lapack_int lda, lapack_int* jpvt, lapack_complex_double* tau)
m_qr = matrix;\ { return LAPACKE_zgeqp3(matrix_layout, m, n, a, lda, jpvt, tau); }
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) template <typename MatrixType>
EIGEN_LAPACKE_QR_COLPIV(float, float, s, ColMajor, LAPACK_COL_MAJOR) struct ColPivHouseholderQR_LAPACKE_impl {
EIGEN_LAPACKE_QR_COLPIV(dcomplex, lapack_complex_double, z, ColMajor, LAPACK_COL_MAJOR) typedef typename MatrixType::Scalar Scalar;
EIGEN_LAPACKE_QR_COLPIV(scomplex, lapack_complex_float, c, ColMajor, LAPACK_COL_MAJOR) typedef typename MatrixType::RealScalar RealScalar;
typedef typename internal::lapacke_helpers::translate_type_imp<Scalar>::type LapackeType;
static constexpr int LapackeStorage = MatrixType::IsRowMajor ? (LAPACK_ROW_MAJOR) : (LAPACK_COL_MAJOR);
EIGEN_LAPACKE_QR_COLPIV(double, double, d, RowMajor, LAPACK_ROW_MAJOR) typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
EIGEN_LAPACKE_QR_COLPIV(float, float, s, RowMajor, LAPACK_ROW_MAJOR) typedef PermutationMatrix<Dynamic, Dynamic, lapack_int> PermutationType;
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 static void run(MatrixType& qr, HCoeffsType& hCoeffs, PermutationType& colsPermutation, Index& nonzero_pivots,
RealScalar& maxpivot, bool usePrescribedThreshold, RealScalar prescribedThreshold, Index& det_p,
bool& isInitialized) {
#endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_LAPACKE_H isInitialized = false;
hCoeffs.resize(qr.diagonalSize());
nonzero_pivots = 0;
maxpivot = RealScalar(0);
colsPermutation.resize(qr.cols());
colsPermutation.indices().setZero();
lapack_int rows = internal::lapacke_helpers::to_lapack(qr.rows());
lapack_int cols = internal::lapacke_helpers::to_lapack(qr.cols());
LapackeType* qr_data = (LapackeType*)(qr.data());
lapack_int lda = internal::lapacke_helpers::to_lapack(qr.outerStride());
lapack_int* perm_data = colsPermutation.indices().data();
LapackeType* hCoeffs_data = (LapackeType*)(hCoeffs.data());
lapack_int info = call_geqp3(LapackeStorage, rows, cols, qr_data, lda, perm_data, hCoeffs_data);
if (info != 0) return;
maxpivot = qr.diagonal().cwiseAbs().maxCoeff();
hCoeffs.adjointInPlace();
RealScalar defaultThreshold = NumTraits<RealScalar>::epsilon() * RealScalar(qr.diagonalSize());
RealScalar threshold = usePrescribedThreshold ? prescribedThreshold : defaultThreshold;
RealScalar premultiplied_threshold = maxpivot * threshold;
nonzero_pivots = (qr.diagonal().cwiseAbs().array() > premultiplied_threshold).count();
colsPermutation.indices().array() -= 1;
det_p = colsPermutation.determinant();
isInitialized = true;
};
static void init(Index rows, Index cols, HCoeffsType& hCoeffs, PermutationType& colsPermutation,
bool& usePrescribedThreshold, bool& isInitialized) {
Index diag = numext::mini(rows, cols);
hCoeffs = HCoeffsType(diag);
colsPermutation = PermutationType(cols);
usePrescribedThreshold = false;
isInitialized = false;
}
};
#define COLPIVQR_LAPACKE_COMPUTEINPLACE(EIGTYPE) \
template <> inline void ColPivHouseholderQR<EIGTYPE, lapack_int>::computeInPlace() { \
ColPivHouseholderQR_LAPACKE_impl<MatrixType>::run(m_qr, m_hCoeffs, m_colsPermutation, m_nonzero_pivots, \
m_maxpivot, m_usePrescribedThreshold, m_prescribedThreshold, \
m_det_p, m_isInitialized); } \
#define COLPIVQR_LAPACKE_INIT(EIGTYPE) \
template <> inline void ColPivHouseholderQR<EIGTYPE, lapack_int>::init(Index rows, Index cols) { \
ColPivHouseholderQR_LAPACKE_impl<MatrixType>::init(rows, cols, m_hCoeffs, m_colsPermutation, m_isInitialized, \
m_usePrescribedThreshold); } \
#define COLPIVQR_LAPACKE(EIGTYPE) \
COLPIVQR_LAPACKE_COMPUTEINPLACE(EIGTYPE) \
COLPIVQR_LAPACKE_INIT(EIGTYPE) \
typedef Matrix<float, Dynamic, Dynamic, RowMajor> MatrixXfR;
typedef Matrix<double, Dynamic, Dynamic, RowMajor> MatrixXdR;
typedef Matrix<scomplex, Dynamic, Dynamic, RowMajor> MatrixXcfR;
typedef Matrix<dcomplex, Dynamic, Dynamic, RowMajor> MatrixXcdR;
COLPIVQR_LAPACKE(MatrixXf)
COLPIVQR_LAPACKE(MatrixXd)
COLPIVQR_LAPACKE(MatrixXcf)
COLPIVQR_LAPACKE(MatrixXcd)
COLPIVQR_LAPACKE(MatrixXfR)
COLPIVQR_LAPACKE(MatrixXdR)
COLPIVQR_LAPACKE(MatrixXcfR)
COLPIVQR_LAPACKE(MatrixXcdR)
COLPIVQR_LAPACKE(Ref<MatrixXf>)
COLPIVQR_LAPACKE(Ref<MatrixXd>)
COLPIVQR_LAPACKE(Ref<MatrixXcf>)
COLPIVQR_LAPACKE(Ref<MatrixXcd>)
COLPIVQR_LAPACKE(Ref<MatrixXfR>)
COLPIVQR_LAPACKE(Ref<MatrixXdR>)
COLPIVQR_LAPACKE(Ref<MatrixXcfR>)
COLPIVQR_LAPACKE(Ref<MatrixXcdR>)
#endif
} // end namespace Eigen
#endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_LAPACKE_H

73
Eigen/src/QR/CompleteOrthogonalDecomposition.h
View File

@ -15,12 +15,12 @@
namespace Eigen { namespace Eigen {
namespace internal { namespace internal {
template <typename MatrixType_> template <typename MatrixType_, typename PermutationIndex_>
struct traits<CompleteOrthogonalDecomposition<MatrixType_> > struct traits<CompleteOrthogonalDecomposition<MatrixType_, PermutationIndex_> >
: traits<MatrixType_> { : traits<MatrixType_> {
typedef MatrixXpr XprKind; typedef MatrixXpr XprKind;
typedef SolverStorage StorageKind; typedef SolverStorage StorageKind;
typedef int StorageIndex; typedef PermutationIndex_ PermutationIndex;
enum { Flags = 0 }; enum { Flags = 0 };
}; };
@ -49,8 +49,8 @@ struct traits<CompleteOrthogonalDecomposition<MatrixType_> >
* *
* \sa MatrixBase::completeOrthogonalDecomposition() * \sa MatrixBase::completeOrthogonalDecomposition()
*/ */
template <typename MatrixType_> class CompleteOrthogonalDecomposition template <typename MatrixType_, typename PermutationIndex_> class CompleteOrthogonalDecomposition
: public SolverBase<CompleteOrthogonalDecomposition<MatrixType_> > : public SolverBase<CompleteOrthogonalDecomposition<MatrixType_, PermutationIndex_> >
{ {
public: public:
typedef MatrixType_ MatrixType; typedef MatrixType_ MatrixType;
@ -58,14 +58,14 @@ template <typename MatrixType_> class CompleteOrthogonalDecomposition
template<typename Derived> template<typename Derived>
friend struct internal::solve_assertion; friend struct internal::solve_assertion;
typedef PermutationIndex_ PermutationIndex;
EIGEN_GENERIC_PUBLIC_INTERFACE(CompleteOrthogonalDecomposition) EIGEN_GENERIC_PUBLIC_INTERFACE(CompleteOrthogonalDecomposition)
enum { enum {
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
}; };
typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime, PermutationIndex>
PermutationType; PermutationType;
typedef typename internal::plain_row_type<MatrixType, Index>::type typedef typename internal::plain_row_type<MatrixType, Index>::type
IntRowVectorType; IntRowVectorType;
@ -78,9 +78,6 @@ template <typename MatrixType_> class CompleteOrthogonalDecomposition
HouseholderSequenceType; HouseholderSequenceType;
typedef typename MatrixType::PlainObject PlainObject; typedef typename MatrixType::PlainObject PlainObject;
private:
typedef typename PermutationType::Index PermIndexType;
public: public:
/** /**
* \brief Default Constructor. * \brief Default Constructor.
@ -417,26 +414,26 @@ template <typename MatrixType_> class CompleteOrthogonalDecomposition
template <typename Rhs> template <typename Rhs>
void applyZAdjointOnTheLeftInPlace(Rhs& rhs) const; void applyZAdjointOnTheLeftInPlace(Rhs& rhs) const;
ColPivHouseholderQR<MatrixType> m_cpqr; ColPivHouseholderQR<MatrixType, PermutationIndex> m_cpqr;
HCoeffsType m_zCoeffs; HCoeffsType m_zCoeffs;
RowVectorType m_temp; RowVectorType m_temp;
}; };
template <typename MatrixType> template <typename MatrixType, typename PermutationIndex>
typename MatrixType::Scalar typename MatrixType::Scalar
CompleteOrthogonalDecomposition<MatrixType>::determinant() const { CompleteOrthogonalDecomposition<MatrixType, PermutationIndex>::determinant() const {
return m_cpqr.determinant(); return m_cpqr.determinant();
} }
template <typename MatrixType> template <typename MatrixType, typename PermutationIndex>
typename MatrixType::RealScalar typename MatrixType::RealScalar
CompleteOrthogonalDecomposition<MatrixType>::absDeterminant() const { CompleteOrthogonalDecomposition<MatrixType, PermutationIndex>::absDeterminant() const {
return m_cpqr.absDeterminant(); return m_cpqr.absDeterminant();
} }
template <typename MatrixType> template <typename MatrixType, typename PermutationIndex>
typename MatrixType::RealScalar typename MatrixType::RealScalar
CompleteOrthogonalDecomposition<MatrixType>::logAbsDeterminant() const { CompleteOrthogonalDecomposition<MatrixType, PermutationIndex>::logAbsDeterminant() const {
return m_cpqr.logAbsDeterminant(); return m_cpqr.logAbsDeterminant();
} }
@ -447,11 +444,10 @@ CompleteOrthogonalDecomposition<MatrixType>::logAbsDeterminant() const {
* \sa class CompleteOrthogonalDecomposition, * \sa class CompleteOrthogonalDecomposition,
* CompleteOrthogonalDecomposition(const MatrixType&) * CompleteOrthogonalDecomposition(const MatrixType&)
*/ */
template <typename MatrixType> template <typename MatrixType, typename PermutationIndex>
void CompleteOrthogonalDecomposition<MatrixType>::computeInPlace() void CompleteOrthogonalDecomposition<MatrixType, PermutationIndex>::computeInPlace()
{ {
// the column permutation is stored as int indices, so just to be sure: eigen_assert(m_cpqr.cols() <= NumTraits<PermutationIndex>::highest());
eigen_assert(m_cpqr.cols() <= NumTraits<int>::highest());
const Index rank = m_cpqr.rank(); const Index rank = m_cpqr.rank();
const Index cols = m_cpqr.cols(); const Index cols = m_cpqr.cols();
@ -503,9 +499,9 @@ void CompleteOrthogonalDecomposition<MatrixType>::computeInPlace()
} }
} }
template <typename MatrixType> template <typename MatrixType, typename PermutationIndex>
template <bool Conjugate, typename Rhs> template <bool Conjugate, typename Rhs>
void CompleteOrthogonalDecomposition<MatrixType>::applyZOnTheLeftInPlace( void CompleteOrthogonalDecomposition<MatrixType, PermutationIndex>::applyZOnTheLeftInPlace(
Rhs& rhs) const { Rhs& rhs) const {
const Index cols = this->cols(); const Index cols = this->cols();
const Index nrhs = rhs.cols(); const Index nrhs = rhs.cols();
@ -525,9 +521,9 @@ void CompleteOrthogonalDecomposition<MatrixType>::applyZOnTheLeftInPlace(
} }
} }
template <typename MatrixType> template <typename MatrixType, typename PermutationIndex>
template <typename Rhs> template <typename Rhs>
void CompleteOrthogonalDecomposition<MatrixType>::applyZAdjointOnTheLeftInPlace( void CompleteOrthogonalDecomposition<MatrixType, PermutationIndex>::applyZAdjointOnTheLeftInPlace(
Rhs& rhs) const { Rhs& rhs) const {
const Index cols = this->cols(); const Index cols = this->cols();
const Index nrhs = rhs.cols(); const Index nrhs = rhs.cols();
@ -548,9 +544,9 @@ void CompleteOrthogonalDecomposition<MatrixType>::applyZAdjointOnTheLeftInPlace(
} }
#ifndef EIGEN_PARSED_BY_DOXYGEN #ifndef EIGEN_PARSED_BY_DOXYGEN
template <typename MatrixType_> template <typename MatrixType_, typename PermutationIndex_>
template <typename RhsType, typename DstType> template <typename RhsType, typename DstType>
void CompleteOrthogonalDecomposition<MatrixType_>::_solve_impl( void CompleteOrthogonalDecomposition<MatrixType_, PermutationIndex_>::_solve_impl(
const RhsType& rhs, DstType& dst) const { const RhsType& rhs, DstType& dst) const {
const Index rank = this->rank(); const Index rank = this->rank();
if (rank == 0) { if (rank == 0) {
@ -580,9 +576,9 @@ void CompleteOrthogonalDecomposition<MatrixType_>::_solve_impl(
dst = colsPermutation() * dst; dst = colsPermutation() * dst;
} }
template<typename MatrixType_> template<typename MatrixType_, typename PermutationIndex_>
template<bool Conjugate, typename RhsType, typename DstType> template<bool Conjugate, typename RhsType, typename DstType>
void CompleteOrthogonalDecomposition<MatrixType_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const void CompleteOrthogonalDecomposition<MatrixType_, PermutationIndex_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
{ {
const Index rank = this->rank(); const Index rank = this->rank();
@ -611,17 +607,17 @@ void CompleteOrthogonalDecomposition<MatrixType_>::_solve_impl_transposed(const
namespace internal { namespace internal {
template<typename MatrixType> template<typename MatrixType, typename PermutationIndex>
struct traits<Inverse<CompleteOrthogonalDecomposition<MatrixType> > > struct traits<Inverse<CompleteOrthogonalDecomposition<MatrixType, PermutationIndex> > >
: traits<typename Transpose<typename MatrixType::PlainObject>::PlainObject> : traits<typename Transpose<typename MatrixType::PlainObject>::PlainObject>
{ {
enum { Flags = 0 }; enum { Flags = 0 };
}; };
template<typename DstXprType, typename MatrixType> template<typename DstXprType, typename MatrixType, typename PermutationIndex>
struct Assignment<DstXprType, Inverse<CompleteOrthogonalDecomposition<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename CompleteOrthogonalDecomposition<MatrixType>::Scalar>, Dense2Dense> struct Assignment<DstXprType, Inverse<CompleteOrthogonalDecomposition<MatrixType, PermutationIndex> >, internal::assign_op<typename DstXprType::Scalar,typename CompleteOrthogonalDecomposition<MatrixType, PermutationIndex>::Scalar>, Dense2Dense>
{ {
typedef CompleteOrthogonalDecomposition<MatrixType> CodType; typedef CompleteOrthogonalDecomposition<MatrixType, PermutationIndex> CodType;
typedef Inverse<CodType> SrcXprType; typedef Inverse<CodType> SrcXprType;
static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename CodType::Scalar> &) static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename CodType::Scalar> &)
{ {
@ -633,9 +629,9 @@ struct Assignment<DstXprType, Inverse<CompleteOrthogonalDecomposition<MatrixType
} // end namespace internal } // end namespace internal
/** \returns the matrix Q as a sequence of householder transformations */ /** \returns the matrix Q as a sequence of householder transformations */
template <typename MatrixType> template <typename MatrixType, typename PermutationIndex>
typename CompleteOrthogonalDecomposition<MatrixType>::HouseholderSequenceType typename CompleteOrthogonalDecomposition<MatrixType, PermutationIndex>::HouseholderSequenceType
CompleteOrthogonalDecomposition<MatrixType>::householderQ() const { CompleteOrthogonalDecomposition<MatrixType, PermutationIndex>::householderQ() const {
return m_cpqr.householderQ(); return m_cpqr.householderQ();
} }
@ -644,7 +640,8 @@ CompleteOrthogonalDecomposition<MatrixType>::householderQ() const {
* \sa class CompleteOrthogonalDecomposition * \sa class CompleteOrthogonalDecomposition
*/ */
template <typename Derived> template <typename Derived>
const CompleteOrthogonalDecomposition<typename MatrixBase<Derived>::PlainObject> template <typename PermutationIndex>
const CompleteOrthogonalDecomposition<typename MatrixBase<Derived>::PlainObject, PermutationIndex>
MatrixBase<Derived>::completeOrthogonalDecomposition() const { MatrixBase<Derived>::completeOrthogonalDecomposition() const {
return CompleteOrthogonalDecomposition<PlainObject>(eval()); return CompleteOrthogonalDecomposition<PlainObject>(eval());
} }

81
Eigen/src/QR/FullPivHouseholderQR.h
View File

@ -17,19 +17,19 @@ namespace Eigen {
namespace internal { namespace internal {
template<typename MatrixType_> struct traits<FullPivHouseholderQR<MatrixType_> > template<typename MatrixType_, typename PermutationIndex_> struct traits<FullPivHouseholderQR<MatrixType_, PermutationIndex_> >
: traits<MatrixType_> : traits<MatrixType_>
{ {
typedef MatrixXpr XprKind; typedef MatrixXpr XprKind;
typedef SolverStorage StorageKind; typedef SolverStorage StorageKind;
typedef int StorageIndex; typedef PermutationIndex_ PermutationIndex;
enum { Flags = 0 }; enum { Flags = 0 };
}; };
template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType; template<typename MatrixType, typename PermutationIndex> struct FullPivHouseholderQRMatrixQReturnType;
template<typename MatrixType> template<typename MatrixType, typename PermutationIndex>
struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> > struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType, PermutationIndex> >
{ {
typedef typename MatrixType::PlainObject ReturnType; typedef typename MatrixType::PlainObject ReturnType;
}; };
@ -59,26 +59,27 @@ struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
* *
* \sa MatrixBase::fullPivHouseholderQr() * \sa MatrixBase::fullPivHouseholderQr()
*/ */
template<typename MatrixType_> class FullPivHouseholderQR template<typename MatrixType_, typename PermutationIndex_> class FullPivHouseholderQR
: public SolverBase<FullPivHouseholderQR<MatrixType_> > : public SolverBase<FullPivHouseholderQR<MatrixType_, PermutationIndex_> >
{ {
public: public:
typedef MatrixType_ MatrixType; typedef MatrixType_ MatrixType;
typedef SolverBase<FullPivHouseholderQR> Base; typedef SolverBase<FullPivHouseholderQR> Base;
friend class SolverBase<FullPivHouseholderQR>; friend class SolverBase<FullPivHouseholderQR>;
typedef PermutationIndex_ PermutationIndex;
EIGEN_GENERIC_PUBLIC_INTERFACE(FullPivHouseholderQR) EIGEN_GENERIC_PUBLIC_INTERFACE(FullPivHouseholderQR)
enum { enum {
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
}; };
typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType; typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType, PermutationIndex> MatrixQReturnType;
typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
typedef Matrix<StorageIndex, 1, typedef Matrix<PermutationIndex, 1,
internal::min_size_prefer_dynamic(ColsAtCompileTime,RowsAtCompileTime), RowMajor, 1, internal::min_size_prefer_dynamic(ColsAtCompileTime,RowsAtCompileTime), RowMajor, 1,
internal::min_size_prefer_fixed(MaxColsAtCompileTime, MaxRowsAtCompileTime)> IntDiagSizeVectorType; internal::min_size_prefer_fixed(MaxColsAtCompileTime, MaxRowsAtCompileTime)> IntDiagSizeVectorType;
typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType; typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime, PermutationIndex> PermutationType;
typedef typename internal::plain_row_type<MatrixType>::type RowVectorType; typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
typedef typename internal::plain_col_type<MatrixType>::type ColVectorType; typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
typedef typename MatrixType::PlainObject PlainObject; typedef typename MatrixType::PlainObject PlainObject;
@ -437,8 +438,8 @@ template<typename MatrixType_> class FullPivHouseholderQR
Index m_det_p; Index m_det_p;
}; };
template<typename MatrixType> template<typename MatrixType, typename PermutationIndex>
typename MatrixType::Scalar FullPivHouseholderQR<MatrixType>::determinant() const typename MatrixType::Scalar FullPivHouseholderQR<MatrixType, PermutationIndex>::determinant() const
{ {
eigen_assert(m_isInitialized && "HouseholderQR is not initialized."); 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!"); eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
@ -447,8 +448,8 @@ typename MatrixType::Scalar FullPivHouseholderQR<MatrixType>::determinant() cons
return m_qr.diagonal().prod() * detQ * Scalar(m_det_p); return m_qr.diagonal().prod() * detQ * Scalar(m_det_p);
} }
template<typename MatrixType> template<typename MatrixType, typename PermutationIndex>
typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType, PermutationIndex>::absDeterminant() const
{ {
using std::abs; using std::abs;
eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
@ -456,8 +457,8 @@ typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant
return abs(m_qr.diagonal().prod()); return abs(m_qr.diagonal().prod());
} }
template<typename MatrixType> template<typename MatrixType, typename PermutationIndex>
typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDeterminant() const typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType, PermutationIndex>::logAbsDeterminant() const
{ {
eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); 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!"); eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
@ -470,18 +471,19 @@ typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDetermin
* *
* \sa class FullPivHouseholderQR, FullPivHouseholderQR(const MatrixType&) * \sa class FullPivHouseholderQR, FullPivHouseholderQR(const MatrixType&)
*/ */
template<typename MatrixType> template<typename MatrixType, typename PermutationIndex>
template<typename InputType> template<typename InputType>
FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(const EigenBase<InputType>& matrix) FullPivHouseholderQR<MatrixType, PermutationIndex>& FullPivHouseholderQR<MatrixType, PermutationIndex>::compute(const EigenBase<InputType>& matrix)
{ {
m_qr = matrix.derived(); m_qr = matrix.derived();
computeInPlace(); computeInPlace();
return *this; return *this;
} }
template<typename MatrixType> template<typename MatrixType, typename PermutationIndex>
void FullPivHouseholderQR<MatrixType>::computeInPlace() void FullPivHouseholderQR<MatrixType, PermutationIndex>::computeInPlace()
{ {
eigen_assert(m_qr.cols() <= NumTraits<PermutationIndex>::highest());
using std::abs; using std::abs;
Index rows = m_qr.rows(); Index rows = m_qr.rows();
Index cols = m_qr.cols(); Index cols = m_qr.cols();
@ -523,15 +525,15 @@ void FullPivHouseholderQR<MatrixType>::computeInPlace()
m_nonzero_pivots = k; m_nonzero_pivots = k;
for(Index i = k; i < size; i++) for(Index i = k; i < size; i++)
{ {
m_rows_transpositions.coeffRef(i) = internal::convert_index<StorageIndex>(i); m_rows_transpositions.coeffRef(i) = internal::convert_index<PermutationIndex>(i);
m_cols_transpositions.coeffRef(i) = internal::convert_index<StorageIndex>(i); m_cols_transpositions.coeffRef(i) = internal::convert_index<PermutationIndex>(i);
m_hCoeffs.coeffRef(i) = Scalar(0); m_hCoeffs.coeffRef(i) = Scalar(0);
} }
break; break;
} }
m_rows_transpositions.coeffRef(k) = internal::convert_index<StorageIndex>(row_of_biggest_in_corner); m_rows_transpositions.coeffRef(k) = internal::convert_index<PermutationIndex>(row_of_biggest_in_corner);
m_cols_transpositions.coeffRef(k) = internal::convert_index<StorageIndex>(col_of_biggest_in_corner); m_cols_transpositions.coeffRef(k) = internal::convert_index<PermutationIndex>(col_of_biggest_in_corner);
if(k != row_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)); m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
++number_of_transpositions; ++number_of_transpositions;
@ -561,9 +563,9 @@ void FullPivHouseholderQR<MatrixType>::computeInPlace()
} }
#ifndef EIGEN_PARSED_BY_DOXYGEN #ifndef EIGEN_PARSED_BY_DOXYGEN
template<typename MatrixType_> template<typename MatrixType_, typename PermutationIndex_>
template<typename RhsType, typename DstType> template<typename RhsType, typename DstType>
void FullPivHouseholderQR<MatrixType_>::_solve_impl(const RhsType &rhs, DstType &dst) const void FullPivHouseholderQR<MatrixType_, PermutationIndex_>::_solve_impl(const RhsType &rhs, DstType &dst) const
{ {
const Index l_rank = rank(); const Index l_rank = rank();
@ -595,9 +597,9 @@ void FullPivHouseholderQR<MatrixType_>::_solve_impl(const RhsType &rhs, DstType
for(Index i = l_rank; i < cols(); ++i) dst.row(m_cols_permutation.indices().coeff(i)).setZero(); for(Index i = l_rank; i < cols(); ++i) dst.row(m_cols_permutation.indices().coeff(i)).setZero();
} }
template<typename MatrixType_> template<typename MatrixType_, typename PermutationIndex_>
template<bool Conjugate, typename RhsType, typename DstType> template<bool Conjugate, typename RhsType, typename DstType>
void FullPivHouseholderQR<MatrixType_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const void FullPivHouseholderQR<MatrixType_, PermutationIndex_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
{ {
const Index l_rank = rank(); const Index l_rank = rank();
@ -634,10 +636,10 @@ void FullPivHouseholderQR<MatrixType_>::_solve_impl_transposed(const RhsType &rh
namespace internal { namespace internal {
template<typename DstXprType, typename MatrixType> template<typename DstXprType, typename MatrixType, typename PermutationIndex>
struct Assignment<DstXprType, Inverse<FullPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense> struct Assignment<DstXprType, Inverse<FullPivHouseholderQR<MatrixType, PermutationIndex> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivHouseholderQR<MatrixType, PermutationIndex>::Scalar>, Dense2Dense>
{ {
typedef FullPivHouseholderQR<MatrixType> QrType; typedef FullPivHouseholderQR<MatrixType, PermutationIndex> QrType;
typedef Inverse<QrType> SrcXprType; typedef Inverse<QrType> SrcXprType;
static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &) static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &)
{ {
@ -651,11 +653,11 @@ struct Assignment<DstXprType, Inverse<FullPivHouseholderQR<MatrixType> >, intern
* *
* \tparam MatrixType type of underlying dense matrix * \tparam MatrixType type of underlying dense matrix
*/ */
template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType template<typename MatrixType, typename PermutationIndex> struct FullPivHouseholderQRMatrixQReturnType
: public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> > : public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType, PermutationIndex> >
{ {
public: public:
typedef typename FullPivHouseholderQR<MatrixType>::IntDiagSizeVectorType IntDiagSizeVectorType; typedef typename FullPivHouseholderQR<MatrixType, PermutationIndex>::IntDiagSizeVectorType IntDiagSizeVectorType;
typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1, typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1,
MatrixType::MaxRowsAtCompileTime> WorkVectorType; MatrixType::MaxRowsAtCompileTime> WorkVectorType;
@ -712,8 +714,8 @@ protected:
} // end namespace internal } // end namespace internal
template<typename MatrixType> template<typename MatrixType, typename PermutationIndex>
inline typename FullPivHouseholderQR<MatrixType>::MatrixQReturnType FullPivHouseholderQR<MatrixType>::matrixQ() const inline typename FullPivHouseholderQR<MatrixType, PermutationIndex>::MatrixQReturnType FullPivHouseholderQR<MatrixType, PermutationIndex>::matrixQ() const
{ {
eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions); return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions);
@ -724,10 +726,11 @@ inline typename FullPivHouseholderQR<MatrixType>::MatrixQReturnType FullPivHouse
* \sa class FullPivHouseholderQR * \sa class FullPivHouseholderQR
*/ */
template<typename Derived> template<typename Derived>
const FullPivHouseholderQR<typename MatrixBase<Derived>::PlainObject> template<typename PermutationIndex>
const FullPivHouseholderQR<typename MatrixBase<Derived>::PlainObject, PermutationIndex>
MatrixBase<Derived>::fullPivHouseholderQr() const MatrixBase<Derived>::fullPivHouseholderQr() const
{ {
return FullPivHouseholderQR<PlainObject>(eval()); return FullPivHouseholderQR<PlainObject, PermutationIndex>(eval());
} }
} // end namespace Eigen } // end namespace Eigen

6
Eigen/src/misc/lapacke_helpers.h
View File

@ -76,12 +76,6 @@ EIGEN_ALWAYS_INLINE EIGEN_CONSTEXPR lapack_int lapack_storage_of(const EigenBase
return Derived::IsRowMajor ? LAPACK_ROW_MAJOR : LAPACK_COL_MAJOR; return Derived::IsRowMajor ? LAPACK_ROW_MAJOR : LAPACK_COL_MAJOR;
} }
/// translate UpLo type to the corresponding letter code
template<UpLoType mode> char translate_mode;
template<> constexpr char translate_mode<Lower> = 'L';
template<> constexpr char translate_mode<Upper> = 'U';
// --------------------------------------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------------------------------------
// Automatic generation of low-level wrappers // Automatic generation of low-level wrappers
// --------------------------------------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------------------------------------

105
test/qr_colpivoting.cpp
View File

@ -13,9 +13,8 @@
#include <Eigen/SVD> #include <Eigen/SVD>
#include "solverbase.h" #include "solverbase.h"
template <typename MatrixType> template <typename MatrixType, typename PermutationIndex>
void cod() { void cod() {
STATIC_CHECK(( internal::is_same<typename CompleteOrthogonalDecomposition<MatrixType>::StorageIndex,int>::value ));
Index rows = internal::random<Index>(2, EIGEN_TEST_MAX_SIZE); Index rows = internal::random<Index>(2, EIGEN_TEST_MAX_SIZE);
Index cols = internal::random<Index>(2, EIGEN_TEST_MAX_SIZE); Index cols = internal::random<Index>(2, EIGEN_TEST_MAX_SIZE);
@ -28,7 +27,7 @@ void cod() {
MatrixQType; MatrixQType;
MatrixType matrix; MatrixType matrix;
createRandomPIMatrixOfRank(rank, rows, cols, matrix); createRandomPIMatrixOfRank(rank, rows, cols, matrix);
CompleteOrthogonalDecomposition<MatrixType> cod(matrix); CompleteOrthogonalDecomposition<MatrixType, PermutationIndex> cod(matrix);
VERIFY(rank == cod.rank()); VERIFY(rank == cod.rank());
VERIFY(cols - cod.rank() == cod.dimensionOfKernel()); VERIFY(cols - cod.rank() == cod.dimensionOfKernel());
VERIFY(!cod.isInjective()); VERIFY(!cod.isInjective());
@ -63,14 +62,14 @@ void cod() {
VERIFY_IS_APPROX(cod_solution, pinv * rhs); VERIFY_IS_APPROX(cod_solution, pinv * rhs);
} }
template <typename MatrixType, int Cols2> template <typename MatrixType, typename PermutationIndex, int Cols2>
void cod_fixedsize() { void cod_fixedsize() {
enum { enum {
Rows = MatrixType::RowsAtCompileTime, Rows = MatrixType::RowsAtCompileTime,
Cols = MatrixType::ColsAtCompileTime Cols = MatrixType::ColsAtCompileTime
}; };
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef CompleteOrthogonalDecomposition<Matrix<Scalar, Rows, Cols> > COD; typedef CompleteOrthogonalDecomposition<Matrix<Scalar, Rows, Cols>, PermutationIndex> COD;
int rank = internal::random<int>(1, (std::min)(int(Rows), int(Cols)) - 1); int rank = internal::random<int>(1, (std::min)(int(Rows), int(Cols)) - 1);
Matrix<Scalar, Rows, Cols> matrix; Matrix<Scalar, Rows, Cols> matrix;
createRandomPIMatrixOfRank(rank, Rows, Cols, matrix); createRandomPIMatrixOfRank(rank, Rows, Cols, matrix);
@ -96,12 +95,10 @@ void cod_fixedsize() {
VERIFY_IS_APPROX(cod_solution, pinv * rhs); VERIFY_IS_APPROX(cod_solution, pinv * rhs);
} }
template<typename MatrixType> void qr() template<typename MatrixType, typename PermutationIndex> void qr()
{ {
using std::sqrt; using std::sqrt;
STATIC_CHECK(( internal::is_same<typename ColPivHouseholderQR<MatrixType>::StorageIndex,int>::value ));
Index rows = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE), cols = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE), cols2 = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE); Index rows = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE), cols = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE), cols2 = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE);
Index rank = internal::random<Index>(1, (std::min)(rows, cols)-1); Index rank = internal::random<Index>(1, (std::min)(rows, cols)-1);
@ -110,7 +107,7 @@ template<typename MatrixType> void qr()
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType; typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType;
MatrixType m1; MatrixType m1;
createRandomPIMatrixOfRank(rank,rows,cols,m1); createRandomPIMatrixOfRank(rank,rows,cols,m1);
ColPivHouseholderQR<MatrixType> qr(m1); ColPivHouseholderQR<MatrixType, PermutationIndex> qr(m1);
VERIFY_IS_EQUAL(rank, qr.rank()); VERIFY_IS_EQUAL(rank, qr.rank());
VERIFY_IS_EQUAL(cols - qr.rank(), qr.dimensionOfKernel()); VERIFY_IS_EQUAL(cols - qr.rank(), qr.dimensionOfKernel());
VERIFY(!qr.isInjective()); VERIFY(!qr.isInjective());
@ -158,7 +155,7 @@ template<typename MatrixType> void qr()
} }
} }
template<typename MatrixType, int Cols2> void qr_fixedsize() template<typename MatrixType, typename PermutationIndex, int Cols2> void qr_fixedsize()
{ {
using std::sqrt; using std::sqrt;
using std::abs; using std::abs;
@ -168,7 +165,7 @@ template<typename MatrixType, int Cols2> void qr_fixedsize()
int rank = internal::random<int>(1, (std::min)(int(Rows), int(Cols))-1); int rank = internal::random<int>(1, (std::min)(int(Rows), int(Cols))-1);
Matrix<Scalar,Rows,Cols> m1; Matrix<Scalar,Rows,Cols> m1;
createRandomPIMatrixOfRank(rank,Rows,Cols,m1); createRandomPIMatrixOfRank(rank,Rows,Cols,m1);
ColPivHouseholderQR<Matrix<Scalar,Rows,Cols> > qr(m1); ColPivHouseholderQR<Matrix<Scalar,Rows,Cols>, PermutationIndex> qr(m1);
VERIFY_IS_EQUAL(rank, qr.rank()); VERIFY_IS_EQUAL(rank, qr.rank());
VERIFY_IS_EQUAL(Cols - qr.rank(), qr.dimensionOfKernel()); VERIFY_IS_EQUAL(Cols - qr.rank(), qr.dimensionOfKernel());
VERIFY_IS_EQUAL(qr.isInjective(), (rank == Rows)); VERIFY_IS_EQUAL(qr.isInjective(), (rank == Rows));
@ -207,7 +204,7 @@ template<typename MatrixType, int Cols2> void qr_fixedsize()
// for rank-revealing QR. See // for rank-revealing QR. See
// http://www.netlib.org/lapack/lawnspdf/lawn176.pdf // http://www.netlib.org/lapack/lawnspdf/lawn176.pdf
// page 3 for more detail. // page 3 for more detail.
template<typename MatrixType> void qr_kahan_matrix() template<typename MatrixType, typename PermutationIndex> void qr_kahan_matrix()
{ {
using std::sqrt; using std::sqrt;
using std::abs; using std::abs;
@ -227,7 +224,7 @@ template<typename MatrixType> void qr_kahan_matrix()
pow_s_i *= s; pow_s_i *= s;
} }
m1 = (m1 + m1.transpose()).eval(); m1 = (m1 + m1.transpose()).eval();
ColPivHouseholderQR<MatrixType> qr(m1); ColPivHouseholderQR<MatrixType, PermutationIndex> qr(m1);
MatrixType r = qr.matrixQR().template triangularView<Upper>(); MatrixType r = qr.matrixQR().template triangularView<Upper>();
RealScalar threshold = RealScalar threshold =
@ -247,7 +244,7 @@ template<typename MatrixType> void qr_kahan_matrix()
} }
} }
template<typename MatrixType> void qr_invertible() template<typename MatrixType, typename PermutationIndex> void qr_invertible()
{ {
using std::log; using std::log;
using std::abs; using std::abs;
@ -266,7 +263,7 @@ template<typename MatrixType> void qr_invertible()
m1 += a * a.adjoint(); m1 += a * a.adjoint();
} }
ColPivHouseholderQR<MatrixType> qr(m1); ColPivHouseholderQR<MatrixType, PermutationIndex> qr(m1);
check_solverbase<MatrixType, MatrixType>(m1, qr, size, size, size); check_solverbase<MatrixType, MatrixType>(m1, qr, size, size, size);
@ -283,11 +280,11 @@ template<typename MatrixType> void qr_invertible()
VERIFY_IS_APPROX(log(absdet), qr.logAbsDeterminant()); VERIFY_IS_APPROX(log(absdet), qr.logAbsDeterminant());
} }
template<typename MatrixType> void qr_verify_assert() template<typename MatrixType, typename PermutationIndex> void qr_verify_assert()
{ {
MatrixType tmp; MatrixType tmp;
ColPivHouseholderQR<MatrixType> qr; ColPivHouseholderQR<MatrixType, PermutationIndex> qr;
VERIFY_RAISES_ASSERT(qr.matrixQR()) VERIFY_RAISES_ASSERT(qr.matrixQR())
VERIFY_RAISES_ASSERT(qr.solve(tmp)) VERIFY_RAISES_ASSERT(qr.solve(tmp))
VERIFY_RAISES_ASSERT(qr.transpose().solve(tmp)) VERIFY_RAISES_ASSERT(qr.transpose().solve(tmp))
@ -303,11 +300,11 @@ template<typename MatrixType> void qr_verify_assert()
VERIFY_RAISES_ASSERT(qr.logAbsDeterminant()) VERIFY_RAISES_ASSERT(qr.logAbsDeterminant())
} }
template<typename MatrixType> void cod_verify_assert() template<typename MatrixType, typename PermutationIndex> void cod_verify_assert()
{ {
MatrixType tmp; MatrixType tmp;
CompleteOrthogonalDecomposition<MatrixType> cod; CompleteOrthogonalDecomposition<MatrixType, PermutationIndex> cod;
VERIFY_RAISES_ASSERT(cod.matrixQTZ()) VERIFY_RAISES_ASSERT(cod.matrixQTZ())
VERIFY_RAISES_ASSERT(cod.solve(tmp)) VERIFY_RAISES_ASSERT(cod.solve(tmp))
VERIFY_RAISES_ASSERT(cod.transpose().solve(tmp)) VERIFY_RAISES_ASSERT(cod.transpose().solve(tmp))
@ -323,50 +320,58 @@ template<typename MatrixType> void cod_verify_assert()
VERIFY_RAISES_ASSERT(cod.logAbsDeterminant()) VERIFY_RAISES_ASSERT(cod.logAbsDeterminant())
} }
EIGEN_DECLARE_TEST(qr_colpivoting) EIGEN_DECLARE_TEST(qr_colpivoting)
{ {
#if defined(EIGEN_USE_LAPACKE)
typedef lapack_int PermutationIndex;
#else
typedef int PermutationIndex;
#endif
for(int i = 0; i < g_repeat; i++) { for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1( qr<MatrixXf>() ); CALL_SUBTEST_1( (qr<MatrixXf, PermutationIndex>)() );
CALL_SUBTEST_2( qr<MatrixXd>() ); CALL_SUBTEST_2( (qr<MatrixXf, PermutationIndex>)() );
CALL_SUBTEST_3( qr<MatrixXcd>() ); CALL_SUBTEST_3( (qr<MatrixXf, PermutationIndex>)() );
CALL_SUBTEST_4(( qr_fixedsize<Matrix<float,3,5>, 4 >() )); CALL_SUBTEST_4(( (qr_fixedsize<Matrix<float,3,5>, PermutationIndex, 4>)() ));
CALL_SUBTEST_5(( qr_fixedsize<Matrix<double,6,2>, 3 >() )); CALL_SUBTEST_5(( (qr_fixedsize<Matrix<double,6,2>, PermutationIndex, 3>)() ));
CALL_SUBTEST_5(( qr_fixedsize<Matrix<double,1,1>, 1 >() )); CALL_SUBTEST_5(( (qr_fixedsize<Matrix<double,1,1>, PermutationIndex, 1>)() ));
} }
for(int i = 0; i < g_repeat; i++) { for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1( cod<MatrixXf>() ); CALL_SUBTEST_1( (cod<MatrixXf, PermutationIndex>)() );
CALL_SUBTEST_2( cod<MatrixXd>() ); CALL_SUBTEST_2( (cod<MatrixXf, PermutationIndex>)() );
CALL_SUBTEST_3( cod<MatrixXcd>() ); CALL_SUBTEST_3( (cod<MatrixXf, PermutationIndex>)() );
CALL_SUBTEST_4(( cod_fixedsize<Matrix<float,3,5>, 4 >() )); CALL_SUBTEST_4(( (cod_fixedsize<Matrix<float,3,5>, PermutationIndex, 4>)() ));
CALL_SUBTEST_5(( cod_fixedsize<Matrix<double,6,2>, 3 >() )); CALL_SUBTEST_5(( (cod_fixedsize<Matrix<double,6,2>, PermutationIndex, 3>)() ));
CALL_SUBTEST_5(( cod_fixedsize<Matrix<double,1,1>, 1 >() )); CALL_SUBTEST_5(( (cod_fixedsize<Matrix<double,1,1>, PermutationIndex, 1>)() ));
} }
for(int i = 0; i < g_repeat; i++) { for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1( qr_invertible<MatrixXf>() ); CALL_SUBTEST_1( (qr_invertible<MatrixXf, PermutationIndex>)() );
CALL_SUBTEST_2( qr_invertible<MatrixXd>() ); CALL_SUBTEST_2( (qr_invertible<MatrixXd, PermutationIndex>)() );
CALL_SUBTEST_6( qr_invertible<MatrixXcf>() ); CALL_SUBTEST_6( (qr_invertible<MatrixXcf, PermutationIndex>)() );
CALL_SUBTEST_3( qr_invertible<MatrixXcd>() ); CALL_SUBTEST_3( (qr_invertible<MatrixXcd, PermutationIndex>)() );
} }
CALL_SUBTEST_7(qr_verify_assert<Matrix3f>()); CALL_SUBTEST_7( (qr_verify_assert<Matrix3f, PermutationIndex>) ());
CALL_SUBTEST_8(qr_verify_assert<Matrix3d>()); CALL_SUBTEST_8( (qr_verify_assert<Matrix3d, PermutationIndex>)());
CALL_SUBTEST_1(qr_verify_assert<MatrixXf>()); CALL_SUBTEST_1( (qr_verify_assert<MatrixXf, PermutationIndex>)());
CALL_SUBTEST_2(qr_verify_assert<MatrixXd>()); CALL_SUBTEST_2( (qr_verify_assert<MatrixXd, PermutationIndex>)());
CALL_SUBTEST_6(qr_verify_assert<MatrixXcf>()); CALL_SUBTEST_6( (qr_verify_assert<MatrixXcf, PermutationIndex>)());
CALL_SUBTEST_3(qr_verify_assert<MatrixXcd>()); CALL_SUBTEST_3( (qr_verify_assert<MatrixXcd, PermutationIndex>)());
CALL_SUBTEST_7(cod_verify_assert<Matrix3f>()); CALL_SUBTEST_7( (cod_verify_assert<Matrix3f, PermutationIndex>)());
CALL_SUBTEST_8(cod_verify_assert<Matrix3d>()); CALL_SUBTEST_8( (cod_verify_assert<Matrix3d, PermutationIndex>)());
CALL_SUBTEST_1(cod_verify_assert<MatrixXf>()); CALL_SUBTEST_1( (cod_verify_assert<MatrixXf, PermutationIndex>)());
CALL_SUBTEST_2(cod_verify_assert<MatrixXd>()); CALL_SUBTEST_2( (cod_verify_assert<MatrixXd, PermutationIndex>)());
CALL_SUBTEST_6(cod_verify_assert<MatrixXcf>()); CALL_SUBTEST_6( (cod_verify_assert<MatrixXcf, PermutationIndex>)());
CALL_SUBTEST_3(cod_verify_assert<MatrixXcd>()); CALL_SUBTEST_3( (cod_verify_assert<MatrixXcd, PermutationIndex>)());
// Test problem size constructors // Test problem size constructors
CALL_SUBTEST_9(ColPivHouseholderQR<MatrixXf>(10, 20)); CALL_SUBTEST_9((ColPivHouseholderQR<MatrixXf, PermutationIndex>(10, 20)));
CALL_SUBTEST_1( qr_kahan_matrix<MatrixXf>() ); CALL_SUBTEST_1( (qr_kahan_matrix<MatrixXf, PermutationIndex>)() );
CALL_SUBTEST_2( qr_kahan_matrix<MatrixXd>() ); CALL_SUBTEST_2( (qr_kahan_matrix<MatrixXd, PermutationIndex>)() );
} }

57
test/qr_fullpivoting.cpp
View File

@ -12,9 +12,8 @@
#include <Eigen/QR> #include <Eigen/QR>
#include "solverbase.h" #include "solverbase.h"
template<typename MatrixType> void qr() template<typename MatrixType, typename PermutationIndex> void qr()
{ {
STATIC_CHECK(( internal::is_same<typename FullPivHouseholderQR<MatrixType>::StorageIndex,int>::value ));
static const int Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime; static const int Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime;
Index max_size = EIGEN_TEST_MAX_SIZE; Index max_size = EIGEN_TEST_MAX_SIZE;
@ -28,7 +27,7 @@ template<typename MatrixType> void qr()
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType; typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType;
MatrixType m1; MatrixType m1;
createRandomPIMatrixOfRank(rank,rows,cols,m1); createRandomPIMatrixOfRank(rank,rows,cols,m1);
FullPivHouseholderQR<MatrixType> qr(m1); FullPivHouseholderQR<MatrixType, PermutationIndex> qr(m1);
VERIFY_IS_EQUAL(rank, qr.rank()); VERIFY_IS_EQUAL(rank, qr.rank());
VERIFY_IS_EQUAL(cols - qr.rank(), qr.dimensionOfKernel()); VERIFY_IS_EQUAL(cols - qr.rank(), qr.dimensionOfKernel());
VERIFY(!qr.isInjective()); VERIFY(!qr.isInjective());
@ -67,7 +66,7 @@ template<typename MatrixType> void qr()
} }
} }
template<typename MatrixType> void qr_invertible() template<typename MatrixType, typename PermutationIndex> void qr_invertible()
{ {
using std::log; using std::log;
using std::abs; using std::abs;
@ -88,7 +87,7 @@ template<typename MatrixType> void qr_invertible()
m1 += a * a.adjoint(); m1 += a * a.adjoint();
} }
FullPivHouseholderQR<MatrixType> qr(m1); FullPivHouseholderQR<MatrixType, PermutationIndex> qr(m1);
VERIFY(qr.isInjective()); VERIFY(qr.isInjective());
VERIFY(qr.isInvertible()); VERIFY(qr.isInvertible());
VERIFY(qr.isSurjective()); VERIFY(qr.isSurjective());
@ -108,11 +107,11 @@ template<typename MatrixType> void qr_invertible()
VERIFY_IS_APPROX(log(absdet), qr.logAbsDeterminant()); VERIFY_IS_APPROX(log(absdet), qr.logAbsDeterminant());
} }
template<typename MatrixType> void qr_verify_assert() template<typename MatrixType, typename PermutationIndex> void qr_verify_assert()
{ {
MatrixType tmp; MatrixType tmp;
FullPivHouseholderQR<MatrixType> qr; FullPivHouseholderQR<MatrixType, PermutationIndex> qr;
VERIFY_RAISES_ASSERT(qr.matrixQR()) VERIFY_RAISES_ASSERT(qr.matrixQR())
VERIFY_RAISES_ASSERT(qr.solve(tmp)) VERIFY_RAISES_ASSERT(qr.solve(tmp))
VERIFY_RAISES_ASSERT(qr.transpose().solve(tmp)) VERIFY_RAISES_ASSERT(qr.transpose().solve(tmp))
@ -130,33 +129,35 @@ template<typename MatrixType> void qr_verify_assert()
EIGEN_DECLARE_TEST(qr_fullpivoting) EIGEN_DECLARE_TEST(qr_fullpivoting)
{ {
typedef int PermutationIndex;
for(int i = 0; i < 1; i++) { for(int i = 0; i < 1; i++) {
CALL_SUBTEST_5( qr<Matrix3f>() ); CALL_SUBTEST_5( (qr<Matrix3f, PermutationIndex>()) );
CALL_SUBTEST_6( qr<Matrix3d>() ); CALL_SUBTEST_6( (qr<Matrix3d, PermutationIndex>()) );
CALL_SUBTEST_8( qr<Matrix2f>() ); CALL_SUBTEST_8( (qr<Matrix2f, PermutationIndex>()) );
CALL_SUBTEST_1( qr<MatrixXf>() ); CALL_SUBTEST_1( (qr<MatrixXf, PermutationIndex>()) );
CALL_SUBTEST_2( qr<MatrixXd>() ); CALL_SUBTEST_2( (qr<MatrixXd, PermutationIndex>()) );
CALL_SUBTEST_3( qr<MatrixXcd>() ); CALL_SUBTEST_3( (qr<MatrixXcd, PermutationIndex>()) );
} }
for(int i = 0; i < g_repeat; i++) { for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1( qr_invertible<MatrixXf>() ); CALL_SUBTEST_1( (qr_invertible<MatrixXf, PermutationIndex>()) );
CALL_SUBTEST_2( qr_invertible<MatrixXd>() ); CALL_SUBTEST_2( (qr_invertible<MatrixXd, PermutationIndex>()) );
CALL_SUBTEST_4( qr_invertible<MatrixXcf>() ); CALL_SUBTEST_4( (qr_invertible<MatrixXcf, PermutationIndex>()) );
CALL_SUBTEST_3( qr_invertible<MatrixXcd>() ); CALL_SUBTEST_3( (qr_invertible<MatrixXcd, PermutationIndex>()) );
} }
CALL_SUBTEST_5(qr_verify_assert<Matrix3f>()); CALL_SUBTEST_5( (qr_verify_assert<Matrix3f, PermutationIndex>()) );
CALL_SUBTEST_6(qr_verify_assert<Matrix3d>()); CALL_SUBTEST_6( (qr_verify_assert<Matrix3d, PermutationIndex>()) );
CALL_SUBTEST_1(qr_verify_assert<MatrixXf>()); CALL_SUBTEST_1( (qr_verify_assert<MatrixXf, PermutationIndex>()) );
CALL_SUBTEST_2(qr_verify_assert<MatrixXd>()); CALL_SUBTEST_2( (qr_verify_assert<MatrixXd, PermutationIndex>()) );
CALL_SUBTEST_4(qr_verify_assert<MatrixXcf>()); CALL_SUBTEST_4( (qr_verify_assert<MatrixXcf, PermutationIndex>()) );
CALL_SUBTEST_3(qr_verify_assert<MatrixXcd>()); CALL_SUBTEST_3( (qr_verify_assert<MatrixXcd, PermutationIndex>()) );
// Test problem size constructors // Test problem size constructors
CALL_SUBTEST_7(FullPivHouseholderQR<MatrixXf>(10, 20)); CALL_SUBTEST_7( (FullPivHouseholderQR<MatrixXf,PermutationIndex>(10, 20)));
CALL_SUBTEST_7((FullPivHouseholderQR<Matrix<float,10,20> >(10,20))); CALL_SUBTEST_7( (FullPivHouseholderQR<Matrix<float, 10, 20>, PermutationIndex>(10, 20)));
CALL_SUBTEST_7((FullPivHouseholderQR<Matrix<float,10,20> >(Matrix<float,10,20>::Random()))); CALL_SUBTEST_7( (FullPivHouseholderQR<Matrix<float, 10, 20>, PermutationIndex>(Matrix<float,10,20>::Random())));
CALL_SUBTEST_7((FullPivHouseholderQR<Matrix<float,20,10> >(20,10))); CALL_SUBTEST_7( (FullPivHouseholderQR<Matrix<float, 20, 10>, PermutationIndex>(20, 10)));
CALL_SUBTEST_7((FullPivHouseholderQR<Matrix<float,20,10> >(Matrix<float,20,10>::Random()))); CALL_SUBTEST_7( (FullPivHouseholderQR<Matrix<float, 20, 10>, PermutationIndex>(Matrix<float,20,10>::Random())));
} }