Revert "Add template to specify QR permutation index type, Fix ColPivHouseholderQR Lapacke bindings"

This reverts commit be7791e097c1fc21d4f2e8713467431784f3a4fd
This commit is contained in:
Rasmus Munk Larsen 2023-01-11 18:50:52 +00:00
parent be7791e097
commit 6156797016
10 changed files with 252 additions and 309 deletions

View File

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

View File

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

View File

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

View File

@ -16,12 +16,12 @@
namespace Eigen {
namespace internal {
template<typename MatrixType_, typename StorageIndex_> struct traits<ColPivHouseholderQR<MatrixType_, StorageIndex_> >
template<typename MatrixType_> struct traits<ColPivHouseholderQR<MatrixType_> >
: traits<MatrixType_>
{
typedef MatrixXpr XprKind;
typedef SolverStorage StorageKind;
typedef StorageIndex_ StorageIndex;
typedef int StorageIndex;
enum { Flags = 0 };
};
@ -50,8 +50,8 @@ template<typename MatrixType_, typename StorageIndex_> struct traits<ColPivHouse
*
* \sa MatrixBase::colPivHouseholderQr()
*/
template<typename MatrixType_, typename StorageIndex_> class ColPivHouseholderQR
: public SolverBase<ColPivHouseholderQR<MatrixType_, StorageIndex_> >
template<typename MatrixType_> class ColPivHouseholderQR
: public SolverBase<ColPivHouseholderQR<MatrixType_> >
{
public:
@ -65,13 +65,17 @@ template<typename MatrixType_, typename StorageIndex_> class ColPivHouseholderQR
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
};
typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime, StorageIndex> PermutationType;
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,internal::remove_all_t<typename HCoeffsType::ConjugateReturnType>> HouseholderSequenceType;
typedef typename MatrixType::PlainObject PlainObject;
private:
typedef typename PermutationType::StorageIndex PermIndexType;
public:
/**
@ -100,7 +104,7 @@ template<typename MatrixType_, typename StorageIndex_> class ColPivHouseholderQR
ColPivHouseholderQR(Index rows, Index cols)
: m_qr(rows, cols),
m_hCoeffs((std::min)(rows,cols)),
m_colsPermutation(cols),
m_colsPermutation(PermIndexType(cols)),
m_colsTranspositions(cols),
m_temp(cols),
m_colNormsUpdated(cols),
@ -124,7 +128,7 @@ template<typename MatrixType_, typename StorageIndex_> class ColPivHouseholderQR
explicit ColPivHouseholderQR(const EigenBase<InputType>& matrix)
: m_qr(matrix.rows(), matrix.cols()),
m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
m_colsPermutation(matrix.cols()),
m_colsPermutation(PermIndexType(matrix.cols())),
m_colsTranspositions(matrix.cols()),
m_temp(matrix.cols()),
m_colNormsUpdated(matrix.cols()),
@ -145,7 +149,7 @@ template<typename MatrixType_, typename StorageIndex_> class ColPivHouseholderQR
explicit ColPivHouseholderQR(EigenBase<InputType>& matrix)
: m_qr(matrix.derived()),
m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
m_colsPermutation(matrix.cols()),
m_colsPermutation(PermIndexType(matrix.cols())),
m_colsTranspositions(matrix.cols()),
m_temp(matrix.cols()),
m_colNormsUpdated(matrix.cols()),
@ -437,7 +441,7 @@ template<typename MatrixType_, typename StorageIndex_> class ColPivHouseholderQR
protected:
friend class CompleteOrthogonalDecomposition<MatrixType, StorageIndex>;
friend class CompleteOrthogonalDecomposition<MatrixType>;
EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
@ -456,8 +460,8 @@ template<typename MatrixType_, typename StorageIndex_> class ColPivHouseholderQR
Index m_det_p;
};
template<typename MatrixType, typename StorageIndex>
typename MatrixType::Scalar ColPivHouseholderQR<MatrixType, StorageIndex>::determinant() const
template<typename MatrixType>
typename MatrixType::Scalar ColPivHouseholderQR<MatrixType>::determinant() 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!");
@ -466,8 +470,8 @@ typename MatrixType::Scalar ColPivHouseholderQR<MatrixType, StorageIndex>::deter
return m_qr.diagonal().prod() * detQ * Scalar(m_det_p);
}
template<typename MatrixType, typename StorageIndex>
typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType, StorageIndex>::absDeterminant() const
template<typename MatrixType>
typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const
{
using std::abs;
eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
@ -475,8 +479,8 @@ typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType, StorageIndex>::a
return abs(m_qr.diagonal().prod());
}
template<typename MatrixType, typename StorageIndex>
typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType, StorageIndex>::logAbsDeterminant() const
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!");
@ -489,20 +493,20 @@ typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType, StorageIndex>::l
*
* \sa class ColPivHouseholderQR, ColPivHouseholderQR(const MatrixType&)
*/
template<typename MatrixType, typename StorageIndex>
template<typename MatrixType>
template<typename InputType>
ColPivHouseholderQR<MatrixType, StorageIndex>& ColPivHouseholderQR<MatrixType, StorageIndex>::compute(const EigenBase<InputType>& matrix)
ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const EigenBase<InputType>& matrix)
{
m_qr = matrix.derived();
computeInPlace();
return *this;
}
template<typename MatrixType, typename StorageIndex>
void ColPivHouseholderQR<MatrixType, StorageIndex>::computeInPlace()
template<typename MatrixType>
void ColPivHouseholderQR<MatrixType>::computeInPlace()
{
eigen_assert(m_qr.cols()<=NumTraits<StorageIndex>::highest());
// the column permutation is stored as int indices, so just to be sure:
eigen_assert(m_qr.cols()<=NumTraits<int>::highest());
using std::abs;
@ -591,18 +595,18 @@ void ColPivHouseholderQR<MatrixType, StorageIndex>::computeInPlace()
}
}
m_colsPermutation.setIdentity(cols);
for(Index k = 0; k < size/*m_nonzero_pivots*/; ++k)
m_colsPermutation.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
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_p = (number_of_transpositions%2) ? -1 : 1;
m_isInitialized = true;
}
#ifndef EIGEN_PARSED_BY_DOXYGEN
template<typename MatrixType_, typename StorageIndex_>
template<typename MatrixType_>
template<typename RhsType, typename DstType>
void ColPivHouseholderQR<MatrixType_, StorageIndex_>::_solve_impl(const RhsType &rhs, DstType &dst) const
void ColPivHouseholderQR<MatrixType_>::_solve_impl(const RhsType &rhs, DstType &dst) const
{
const Index nonzero_pivots = nonzeroPivots();
@ -624,9 +628,9 @@ void ColPivHouseholderQR<MatrixType_, StorageIndex_>::_solve_impl(const RhsType
for(Index i = nonzero_pivots; i < cols(); ++i) dst.row(m_colsPermutation.indices().coeff(i)).setZero();
}
template<typename MatrixType_, typename StorageIndex_>
template<typename MatrixType_>
template<bool Conjugate, typename RhsType, typename DstType>
void ColPivHouseholderQR<MatrixType_, StorageIndex_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
void ColPivHouseholderQR<MatrixType_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
{
const Index nonzero_pivots = nonzeroPivots();
@ -652,10 +656,10 @@ void ColPivHouseholderQR<MatrixType_, StorageIndex_>::_solve_impl_transposed(con
namespace internal {
template<typename DstXprType, typename MatrixType, typename StorageIndex>
struct Assignment<DstXprType, Inverse<ColPivHouseholderQR<MatrixType, StorageIndex> >, internal::assign_op<typename DstXprType::Scalar,typename ColPivHouseholderQR<MatrixType, StorageIndex>::Scalar>, Dense2Dense>
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, StorageIndex> QrType;
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> &)
{
@ -668,8 +672,8 @@ struct Assignment<DstXprType, Inverse<ColPivHouseholderQR<MatrixType, StorageInd
/** \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 StorageIndex>
typename ColPivHouseholderQR<MatrixType, StorageIndex>::HouseholderSequenceType ColPivHouseholderQR<MatrixType, StorageIndex>
template<typename MatrixType>
typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
::householderQ() const
{
eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
@ -681,11 +685,10 @@ typename ColPivHouseholderQR<MatrixType, StorageIndex>::HouseholderSequenceType
* \sa class ColPivHouseholderQR
*/
template<typename Derived>
template<typename StorageIndex>
const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject, StorageIndex>
const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
MatrixBase<Derived>::colPivHouseholderQr() const
{
return ColPivHouseholderQR<PlainObject, StorageIndex>(eval());
return ColPivHouseholderQR<PlainObject>(eval());
}
} // end namespace Eigen

View File

@ -36,110 +36,64 @@
#include "./InternalHeaderCheck.h"
namespace Eigen {
namespace Eigen {
#if defined(EIGEN_USE_LAPACKE)
/** \internal Specialization for the data types supported by LAPACKe */
template<typename Scalar>
inline lapack_int call_geqp3(int matrix_layout, lapack_int m, lapack_int n, Scalar* a, lapack_int lda, lapack_int* jpvt, Scalar* tau);
template<>
inline lapack_int call_geqp3(int matrix_layout, lapack_int m, lapack_int n, float* a, lapack_int lda, lapack_int* jpvt, float* tau)
{ 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)
{ return LAPACKE_dgeqp3(matrix_layout, m, n, a, lda, jpvt, tau); }
template<>
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)
{ return LAPACKE_cgeqp3(matrix_layout, m, n, a, lda, jpvt, tau); }
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)
{ return LAPACKE_zgeqp3(matrix_layout, m, n, a, lda, jpvt, tau); }
#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; \
}
template <typename MatrixType>
struct ColPivHouseholderQR_LAPACKE_impl {
typedef typename MatrixType::Scalar Scalar;
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, 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)
typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
typedef PermutationMatrix<Dynamic, Dynamic, lapack_int> PermutationType;
typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
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)
static void run(MatrixType& qr, HCoeffsType& hCoeffs, PermutationType& colsPermutation,
IntRowVectorType& colsTranspositions, Index& nonzero_pivots, RealScalar& maxpivot,
bool usePrescribedThreshold, RealScalar prescribedThreshold, Index& det_p, bool& isInitialized) {
using std::abs;
isInitialized = false;
hCoeffs.resize(qr.diagonalSize());
colsTranspositions.resize(qr.cols());
nonzero_pivots = 0;
maxpivot = RealScalar(0);
colsPermutation.resize(qr.cols());
colsPermutation.indices().setZero();
} // end namespace Eigen
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 = abs(maxpivot) * threshold;
nonzero_pivots = (qr.diagonal().cwiseAbs().array() > premultiplied_threshold).count();
colsPermutation.indices().array() -= 1;
det_p = colsPermutation.determinant();
isInitialized = true;
};
};
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;
template <> inline void ColPivHouseholderQR<MatrixXf, lapack_int>::computeInPlace() {
ColPivHouseholderQR_LAPACKE_impl<MatrixType>::run(m_qr, m_hCoeffs, m_colsPermutation, m_colsTranspositions,
m_nonzero_pivots, m_maxpivot, m_usePrescribedThreshold,
m_prescribedThreshold, m_det_p, m_isInitialized); }
template <> inline void ColPivHouseholderQR<MatrixXd, lapack_int>::computeInPlace() {
ColPivHouseholderQR_LAPACKE_impl<MatrixType>::run(m_qr, m_hCoeffs, m_colsPermutation, m_colsTranspositions,
m_nonzero_pivots, m_maxpivot, m_usePrescribedThreshold,
m_prescribedThreshold, m_det_p, m_isInitialized); }
template <> inline void ColPivHouseholderQR<MatrixXcf, lapack_int>::computeInPlace() {
ColPivHouseholderQR_LAPACKE_impl<MatrixType>::run(m_qr, m_hCoeffs, m_colsPermutation, m_colsTranspositions,
m_nonzero_pivots, m_maxpivot, m_usePrescribedThreshold,
m_prescribedThreshold, m_det_p, m_isInitialized); }
template <> inline void ColPivHouseholderQR<MatrixXcd, lapack_int>::computeInPlace() {
ColPivHouseholderQR_LAPACKE_impl<MatrixType>::run(m_qr, m_hCoeffs, m_colsPermutation, m_colsTranspositions,
m_nonzero_pivots, m_maxpivot, m_usePrescribedThreshold,
m_prescribedThreshold, m_det_p, m_isInitialized); }
template <> inline void ColPivHouseholderQR<MatrixXfR, lapack_int>::computeInPlace() {
ColPivHouseholderQR_LAPACKE_impl<MatrixType>::run(m_qr, m_hCoeffs, m_colsPermutation, m_colsTranspositions,
m_nonzero_pivots, m_maxpivot, m_usePrescribedThreshold,
m_prescribedThreshold, m_det_p, m_isInitialized); }
template <> inline void ColPivHouseholderQR<MatrixXdR, lapack_int>::computeInPlace() {
ColPivHouseholderQR_LAPACKE_impl<MatrixType>::run(m_qr, m_hCoeffs, m_colsPermutation, m_colsTranspositions,
m_nonzero_pivots, m_maxpivot, m_usePrescribedThreshold,
m_prescribedThreshold, m_det_p, m_isInitialized); }
template <> inline void ColPivHouseholderQR<MatrixXcfR, lapack_int>::computeInPlace() {
ColPivHouseholderQR_LAPACKE_impl<MatrixType>::run(m_qr, m_hCoeffs, m_colsPermutation, m_colsTranspositions,
m_nonzero_pivots, m_maxpivot, m_usePrescribedThreshold,
m_prescribedThreshold, m_det_p, m_isInitialized); }
template <> inline void ColPivHouseholderQR<MatrixXcdR, lapack_int>::computeInPlace() {
ColPivHouseholderQR_LAPACKE_impl<MatrixType>::run(m_qr, m_hCoeffs, m_colsPermutation, m_colsTranspositions,
m_nonzero_pivots, m_maxpivot, m_usePrescribedThreshold,
m_prescribedThreshold, m_det_p, m_isInitialized); }
#endif
} // end namespace Eigen
#endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_LAPACKE_H
#endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_LAPACKE_H

View File

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

View File

@ -17,19 +17,19 @@ namespace Eigen {
namespace internal {
template<typename MatrixType_, typename StorageIndex_> struct traits<FullPivHouseholderQR<MatrixType_, StorageIndex_> >
template<typename MatrixType_> struct traits<FullPivHouseholderQR<MatrixType_> >
: traits<MatrixType_>
{
typedef MatrixXpr XprKind;
typedef SolverStorage StorageKind;
typedef StorageIndex_ StorageIndex;
typedef int StorageIndex;
enum { Flags = 0 };
};
template<typename MatrixType, typename StorageIndex> struct FullPivHouseholderQRMatrixQReturnType;
template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType;
template<typename MatrixType, typename StorageIndex>
struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType, StorageIndex> >
template<typename MatrixType>
struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
{
typedef typename MatrixType::PlainObject ReturnType;
};
@ -59,8 +59,8 @@ struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType, StorageIndex> >
*
* \sa MatrixBase::fullPivHouseholderQr()
*/
template<typename MatrixType_, typename StorageIndex_> class FullPivHouseholderQR
: public SolverBase<FullPivHouseholderQR<MatrixType_, StorageIndex_> >
template<typename MatrixType_> class FullPivHouseholderQR
: public SolverBase<FullPivHouseholderQR<MatrixType_> >
{
public:
@ -73,12 +73,12 @@ template<typename MatrixType_, typename StorageIndex_> class FullPivHouseholderQ
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
};
typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType, StorageIndex_> MatrixQReturnType;
typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType;
typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
typedef Matrix<StorageIndex, 1,
internal::min_size_prefer_dynamic(ColsAtCompileTime,RowsAtCompileTime), RowMajor, 1,
internal::min_size_prefer_fixed(MaxColsAtCompileTime, MaxRowsAtCompileTime)> IntDiagSizeVectorType;
typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime, StorageIndex> PermutationType;
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;
@ -437,8 +437,8 @@ template<typename MatrixType_, typename StorageIndex_> class FullPivHouseholderQ
Index m_det_p;
};
template<typename MatrixType, typename StorageIndex>
typename MatrixType::Scalar FullPivHouseholderQR<MatrixType, StorageIndex>::determinant() const
template<typename MatrixType>
typename MatrixType::Scalar FullPivHouseholderQR<MatrixType>::determinant() 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!");
@ -447,8 +447,8 @@ typename MatrixType::Scalar FullPivHouseholderQR<MatrixType, StorageIndex>::dete
return m_qr.diagonal().prod() * detQ * Scalar(m_det_p);
}
template<typename MatrixType, typename StorageIndex>
typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType, StorageIndex>::absDeterminant() const
template<typename MatrixType>
typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const
{
using std::abs;
eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
@ -456,8 +456,8 @@ typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType, StorageIndex>::
return abs(m_qr.diagonal().prod());
}
template<typename MatrixType, typename StorageIndex>
typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType, StorageIndex>::logAbsDeterminant() const
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!");
@ -470,19 +470,18 @@ typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType, StorageIndex>::
*
* \sa class FullPivHouseholderQR, FullPivHouseholderQR(const MatrixType&)
*/
template<typename MatrixType, typename StorageIndex>
template<typename MatrixType>
template<typename InputType>
FullPivHouseholderQR<MatrixType, StorageIndex>& FullPivHouseholderQR<MatrixType, StorageIndex>::compute(const EigenBase<InputType>& matrix)
FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(const EigenBase<InputType>& matrix)
{
m_qr = matrix.derived();
computeInPlace();
return *this;
}
template<typename MatrixType, typename StorageIndex>
void FullPivHouseholderQR<MatrixType, StorageIndex>::computeInPlace()
template<typename MatrixType>
void FullPivHouseholderQR<MatrixType>::computeInPlace()
{
eigen_assert(m_qr.cols() <= NumTraits<StorageIndex>::highest());
using std::abs;
Index rows = m_qr.rows();
Index cols = m_qr.cols();
@ -562,9 +561,9 @@ void FullPivHouseholderQR<MatrixType, StorageIndex>::computeInPlace()
}
#ifndef EIGEN_PARSED_BY_DOXYGEN
template<typename MatrixType_, typename StorageIndex_>
template<typename MatrixType_>
template<typename RhsType, typename DstType>
void FullPivHouseholderQR<MatrixType_, StorageIndex_>::_solve_impl(const RhsType &rhs, DstType &dst) const
void FullPivHouseholderQR<MatrixType_>::_solve_impl(const RhsType &rhs, DstType &dst) const
{
const Index l_rank = rank();
@ -596,9 +595,9 @@ void FullPivHouseholderQR<MatrixType_, StorageIndex_>::_solve_impl(const RhsType
for(Index i = l_rank; i < cols(); ++i) dst.row(m_cols_permutation.indices().coeff(i)).setZero();
}
template<typename MatrixType_, typename StorageIndex_>
template<typename MatrixType_>
template<bool Conjugate, typename RhsType, typename DstType>
void FullPivHouseholderQR<MatrixType_, StorageIndex_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
void FullPivHouseholderQR<MatrixType_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
{
const Index l_rank = rank();
@ -635,10 +634,10 @@ void FullPivHouseholderQR<MatrixType_, StorageIndex_>::_solve_impl_transposed(co
namespace internal {
template<typename DstXprType, typename MatrixType, typename StorageIndex>
struct Assignment<DstXprType, Inverse<FullPivHouseholderQR<MatrixType, StorageIndex> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivHouseholderQR<MatrixType, StorageIndex>::Scalar>, Dense2Dense>
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, StorageIndex> QrType;
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> &)
{
@ -652,11 +651,11 @@ struct Assignment<DstXprType, Inverse<FullPivHouseholderQR<MatrixType, StorageIn
*
* \tparam MatrixType type of underlying dense matrix
*/
template<typename MatrixType, typename StorageIndex> struct FullPivHouseholderQRMatrixQReturnType
: public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType, StorageIndex> >
template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType
: public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
{
public:
typedef typename FullPivHouseholderQR<MatrixType, StorageIndex>::IntDiagSizeVectorType IntDiagSizeVectorType;
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;
@ -713,8 +712,8 @@ protected:
} // end namespace internal
template<typename MatrixType, typename StorageIndex>
inline typename FullPivHouseholderQR<MatrixType, StorageIndex>::MatrixQReturnType FullPivHouseholderQR<MatrixType, StorageIndex>::matrixQ() const
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);
@ -725,11 +724,10 @@ inline typename FullPivHouseholderQR<MatrixType, StorageIndex>::MatrixQReturnTyp
* \sa class FullPivHouseholderQR
*/
template<typename Derived>
template<typename StorageIndex>
const FullPivHouseholderQR<typename MatrixBase<Derived>::PlainObject, StorageIndex>
const FullPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
MatrixBase<Derived>::fullPivHouseholderQr() const
{
return FullPivHouseholderQR<PlainObject, StorageIndex>(eval());
return FullPivHouseholderQR<PlainObject>(eval());
}
} // end namespace Eigen

View File

@ -77,9 +77,9 @@ EIGEN_ALWAYS_INLINE EIGEN_CONSTEXPR lapack_int lapack_storage_of(const EigenBase
}
/// 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';
template<UpLoType mode> char translate_mode;
template<> constexpr char translate_mode<Lower> = 'L';
template<> constexpr char translate_mode<Upper> = 'U';
// ---------------------------------------------------------------------------------------------------------------------

View File

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

View File

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