mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-03-07 18:27:40 +08:00
Code Issues Packages Projects Releases Wiki Activity

Specify Permutation Index for PartialPivLU and FullPivLU

Browse Source
This commit is contained in:
Charles Schlosser 2023-03-07 20:28:05 +00:00 committed by Rasmus Munk Larsen
parent eb4dbf6135
commit 7bf2968fed
5 changed files with 80 additions and 78 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -328,10 +328,10 @@ template<typename Derived> class MatrixBase
/////////// LU module ///////////
inline const FullPivLU<PlainObject> fullPivLu() const;
inline const PartialPivLU<PlainObject> partialPivLu() const;
template<typename PermutationIndex = DefaultPermutationIndex> inline const FullPivLU<PlainObject, PermutationIndex> fullPivLu() const;
template<typename PermutationIndex = DefaultPermutationIndex> inline const PartialPivLU<PlainObject, PermutationIndex> partialPivLu() const;
inline const PartialPivLU<PlainObject> lu() const;
template<typename PermutationIndex = DefaultPermutationIndex> inline const PartialPivLU<PlainObject, PermutationIndex> lu() const;
EIGEN_DEVICE_FUNC
inline const Inverse<Derived> inverse() const;

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

@ -270,8 +270,8 @@ typedef lapack_int DefaultPermutationIndex;
typedef int DefaultPermutationIndex;
#endif
template<typename MatrixType> class FullPivLU;
template<typename MatrixType> class PartialPivLU;
template<typename MatrixType, typename PermutationIndex = DefaultPermutationIndex> class FullPivLU;
template<typename MatrixType, typename PermutationIndex = DefaultPermutationIndex> class PartialPivLU;
namespace internal {
template<typename MatrixType> struct inverse_impl;
}

85
Eigen/src/LU/FullPivLU.h
View File

@ -15,12 +15,12 @@
namespace Eigen {
namespace internal {
template<typename MatrixType_> struct traits<FullPivLU<MatrixType_> >
template<typename MatrixType_, typename PermutationIndex_> struct traits<FullPivLU<MatrixType_, PermutationIndex_> >
: traits<MatrixType_>
{
typedef MatrixXpr XprKind;
typedef SolverStorage StorageKind;
typedef int StorageIndex;
typedef PermutationIndex_ StorageIndex;
enum { Flags = 0 };
};
@ -59,8 +59,8 @@ template<typename MatrixType_> struct traits<FullPivLU<MatrixType_> >
*
* \sa MatrixBase::fullPivLu(), MatrixBase::determinant(), MatrixBase::inverse()
*/
template<typename MatrixType_> class FullPivLU
: public SolverBase<FullPivLU<MatrixType_> >
template<typename MatrixType_, typename PermutationIndex_> class FullPivLU
: public SolverBase<FullPivLU<MatrixType_, PermutationIndex_> >
{
public:
typedef MatrixType_ MatrixType;
@ -72,10 +72,11 @@ template<typename MatrixType_> class FullPivLU
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
};
typedef typename internal::plain_row_type<MatrixType, StorageIndex>::type IntRowVectorType;
typedef typename internal::plain_col_type<MatrixType, StorageIndex>::type IntColVectorType;
typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationQType;
typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationPType;
using PermutationIndex = PermutationIndex_;
typedef typename internal::plain_row_type<MatrixType, PermutationIndex>::type IntRowVectorType;
typedef typename internal::plain_col_type<MatrixType, PermutationIndex>::type IntColVectorType;
typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime, PermutationIndex> PermutationQType;
typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime, PermutationIndex> PermutationPType;
typedef typename MatrixType::PlainObject PlainObject;
/**
@ -437,14 +438,14 @@ template<typename MatrixType_> class FullPivLU
bool m_isInitialized, m_usePrescribedThreshold;
};
template<typename MatrixType>
FullPivLU<MatrixType>::FullPivLU()
template<typename MatrixType, typename PermutationIndex>
FullPivLU<MatrixType, PermutationIndex>::FullPivLU()
: m_isInitialized(false), m_usePrescribedThreshold(false)
{
}
template<typename MatrixType>
FullPivLU<MatrixType>::FullPivLU(Index rows, Index cols)
template<typename MatrixType, typename PermutationIndex>
FullPivLU<MatrixType, PermutationIndex>::FullPivLU(Index rows, Index cols)
: m_lu(rows, cols),
m_p(rows),
m_q(cols),
@ -455,9 +456,9 @@ FullPivLU<MatrixType>::FullPivLU(Index rows, Index cols)
{
}
template<typename MatrixType>
template<typename MatrixType, typename PermutationIndex>
template<typename InputType>
FullPivLU<MatrixType>::FullPivLU(const EigenBase<InputType>& matrix)
FullPivLU<MatrixType, PermutationIndex>::FullPivLU(const EigenBase<InputType>& matrix)
: m_lu(matrix.rows(), matrix.cols()),
m_p(matrix.rows()),
m_q(matrix.cols()),
@ -469,9 +470,9 @@ FullPivLU<MatrixType>::FullPivLU(const EigenBase<InputType>& matrix)
compute(matrix.derived());
}
template<typename MatrixType>
template<typename MatrixType, typename PermutationIndex>
template<typename InputType>
FullPivLU<MatrixType>::FullPivLU(EigenBase<InputType>& matrix)
FullPivLU<MatrixType, PermutationIndex>::FullPivLU(EigenBase<InputType>& matrix)
: m_lu(matrix.derived()),
m_p(matrix.rows()),
m_q(matrix.cols()),
@ -483,11 +484,10 @@ FullPivLU<MatrixType>::FullPivLU(EigenBase<InputType>& matrix)
computeInPlace();
}
template<typename MatrixType>
void FullPivLU<MatrixType>::computeInPlace()
template<typename MatrixType, typename PermutationIndex>
void FullPivLU<MatrixType, PermutationIndex>::computeInPlace()
{
// the permutations are stored as int indices, so just to be sure:
eigen_assert(m_lu.rows()<=NumTraits<int>::highest() && m_lu.cols()<=NumTraits<int>::highest());
eigen_assert(m_lu.rows()<=NumTraits<PermutationIndex>::highest() && m_lu.cols()<=NumTraits<PermutationIndex>::highest());
m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff();
@ -574,8 +574,8 @@ void FullPivLU<MatrixType>::computeInPlace()
m_isInitialized = true;
}
template<typename MatrixType>
typename internal::traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() const
template<typename MatrixType, typename PermutationIndex>
typename internal::traits<MatrixType>::Scalar FullPivLU<MatrixType, PermutationIndex>::determinant() const
{
eigen_assert(m_isInitialized && "LU is not initialized.");
eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!");
@ -585,8 +585,8 @@ typename internal::traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant
/** \returns the matrix represented by the decomposition,
* i.e., it returns the product: \f$ P^{-1} L U Q^{-1} \f$.
* This function is provided for debug purposes. */
template<typename MatrixType>
MatrixType FullPivLU<MatrixType>::reconstructedMatrix() const
template<typename MatrixType, typename PermutationIndex>
MatrixType FullPivLU<MatrixType, PermutationIndex>::reconstructedMatrix() const
{
eigen_assert(m_isInitialized && "LU is not initialized.");
const Index smalldim = (std::min)(m_lu.rows(), m_lu.cols());
@ -610,11 +610,12 @@ MatrixType FullPivLU<MatrixType>::reconstructedMatrix() const
/********* Implementation of kernel() **************************************************/
namespace internal {
template<typename MatrixType_>
struct kernel_retval<FullPivLU<MatrixType_> >
: kernel_retval_base<FullPivLU<MatrixType_> >
template<typename MatrixType_, typename PermutationIndex_>
struct kernel_retval<FullPivLU<MatrixType_, PermutationIndex_> >
: kernel_retval_base<FullPivLU<MatrixType_, PermutationIndex_> >
{
EIGEN_MAKE_KERNEL_HELPERS(FullPivLU<MatrixType_>)
using DecompositionType = FullPivLU<MatrixType_, PermutationIndex_>;
EIGEN_MAKE_KERNEL_HELPERS(DecompositionType)
enum { MaxSmallDimAtCompileTime = min_size_prefer_fixed(
MatrixType::MaxColsAtCompileTime,
@ -696,11 +697,12 @@ struct kernel_retval<FullPivLU<MatrixType_> >
/***** Implementation of image() *****************************************************/
template<typename MatrixType_>
struct image_retval<FullPivLU<MatrixType_> >
: image_retval_base<FullPivLU<MatrixType_> >
template<typename MatrixType_, typename PermutationIndex_>
struct image_retval<FullPivLU<MatrixType_, PermutationIndex_> >
: image_retval_base<FullPivLU<MatrixType_, PermutationIndex_> >
{
EIGEN_MAKE_IMAGE_HELPERS(FullPivLU<MatrixType_>)
using DecompositionType = FullPivLU<MatrixType_, PermutationIndex_>;
EIGEN_MAKE_IMAGE_HELPERS(DecompositionType)
enum { MaxSmallDimAtCompileTime = min_size_prefer_fixed(
MatrixType::MaxColsAtCompileTime,
@ -737,9 +739,9 @@ struct image_retval<FullPivLU<MatrixType_> >
} // end namespace internal
#ifndef EIGEN_PARSED_BY_DOXYGEN
template<typename MatrixType_>
template<typename MatrixType_, typename PermutationIndex_>
template<typename RhsType, typename DstType>
void FullPivLU<MatrixType_>::_solve_impl(const RhsType &rhs, DstType &dst) const
void FullPivLU<MatrixType_, PermutationIndex_>::_solve_impl(const RhsType &rhs, DstType &dst) const
{
/* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}.
* So we proceed as follows:
@ -784,9 +786,9 @@ void FullPivLU<MatrixType_>::_solve_impl(const RhsType &rhs, DstType &dst) const
dst.row(permutationQ().indices().coeff(i)).setZero();
}
template<typename MatrixType_>
template<typename MatrixType_, typename PermutationIndex_>
template<bool Conjugate, typename RhsType, typename DstType>
void FullPivLU<MatrixType_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
void FullPivLU<MatrixType_, PermutationIndex_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
{
/* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1},
* and since permutations are real and unitary, we can write this
@ -842,10 +844,10 @@ namespace internal {
/***** Implementation of inverse() *****************************************************/
template<typename DstXprType, typename MatrixType>
struct Assignment<DstXprType, Inverse<FullPivLU<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivLU<MatrixType>::Scalar>, Dense2Dense>
template<typename DstXprType, typename MatrixType, typename PermutationIndex>
struct Assignment<DstXprType, Inverse<FullPivLU<MatrixType, PermutationIndex> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivLU<MatrixType, PermutationIndex>::Scalar>, Dense2Dense>
{
typedef FullPivLU<MatrixType> LuType;
typedef FullPivLU<MatrixType, PermutationIndex> LuType;
typedef Inverse<LuType> SrcXprType;
static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename MatrixType::Scalar> &)
{
@ -863,10 +865,11 @@ struct Assignment<DstXprType, Inverse<FullPivLU<MatrixType> >, internal::assign_
* \sa class FullPivLU
*/
template<typename Derived>
inline const FullPivLU<typename MatrixBase<Derived>::PlainObject>
template<typename PermutationIndex>
inline const FullPivLU<typename MatrixBase<Derived>::PlainObject, PermutationIndex>
MatrixBase<Derived>::fullPivLu() const
{
return FullPivLU<PlainObject>(eval());
return FullPivLU<PlainObject, PermutationIndex>(eval());
}
} // end namespace Eigen

60
Eigen/src/LU/PartialPivLU.h
View File

@ -16,12 +16,12 @@
namespace Eigen {
namespace internal {
template<typename MatrixType_> struct traits<PartialPivLU<MatrixType_> >
template<typename MatrixType_, typename PermutationIndex_> struct traits<PartialPivLU<MatrixType_, PermutationIndex_> >
: traits<MatrixType_>
{
typedef MatrixXpr XprKind;
typedef SolverStorage StorageKind;
typedef int StorageIndex;
typedef PermutationIndex_ StorageIndex;
typedef traits<MatrixType_> BaseTraits;
enum {
Flags = BaseTraits::Flags & RowMajorBit,
@ -75,8 +75,8 @@ struct enable_if_ref<Ref<T>,Derived> {
*
* \sa MatrixBase::partialPivLu(), MatrixBase::determinant(), MatrixBase::inverse(), MatrixBase::computeInverse(), class FullPivLU
*/
template<typename MatrixType_> class PartialPivLU
: public SolverBase<PartialPivLU<MatrixType_> >
template<typename MatrixType_, typename PermutationIndex_> class PartialPivLU
: public SolverBase<PartialPivLU<MatrixType_, PermutationIndex_> >
{
public:
@ -89,8 +89,9 @@ template<typename MatrixType_> class PartialPivLU
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
};
typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
using PermutationIndex = PermutationIndex_;
typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime, PermutationIndex> PermutationType;
typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime, PermutationIndex> TranspositionType;
typedef typename MatrixType::PlainObject PlainObject;
/**
@ -279,8 +280,8 @@ template<typename MatrixType_> class PartialPivLU
bool m_isInitialized;
};
template<typename MatrixType>
PartialPivLU<MatrixType>::PartialPivLU()
template<typename MatrixType, typename PermutationIndex>
PartialPivLU<MatrixType, PermutationIndex>::PartialPivLU()
: m_lu(),
m_p(),
m_rowsTranspositions(),
@ -290,8 +291,8 @@ PartialPivLU<MatrixType>::PartialPivLU()
{
}
template<typename MatrixType>
PartialPivLU<MatrixType>::PartialPivLU(Index size)
template<typename MatrixType, typename PermutationIndex>
PartialPivLU<MatrixType, PermutationIndex>::PartialPivLU(Index size)
: m_lu(size, size),
m_p(size),
m_rowsTranspositions(size),
@ -301,9 +302,9 @@ PartialPivLU<MatrixType>::PartialPivLU(Index size)
{
}
template<typename MatrixType>
template<typename MatrixType, typename PermutationIndex>
template<typename InputType>
PartialPivLU<MatrixType>::PartialPivLU(const EigenBase<InputType>& matrix)
PartialPivLU<MatrixType, PermutationIndex>::PartialPivLU(const EigenBase<InputType>& matrix)
: m_lu(matrix.rows(),matrix.cols()),
m_p(matrix.rows()),
m_rowsTranspositions(matrix.rows()),
@ -314,9 +315,9 @@ PartialPivLU<MatrixType>::PartialPivLU(const EigenBase<InputType>& matrix)
compute(matrix.derived());
}
template<typename MatrixType>
template<typename MatrixType, typename PermutationIndex>
template<typename InputType>
PartialPivLU<MatrixType>::PartialPivLU(EigenBase<InputType>& matrix)
PartialPivLU<MatrixType, PermutationIndex>::PartialPivLU(EigenBase<InputType>& matrix)
: m_lu(matrix.derived()),
m_p(matrix.rows()),
m_rowsTranspositions(matrix.rows()),
@ -520,11 +521,10 @@ void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, t
} // end namespace internal
template<typename MatrixType>
void PartialPivLU<MatrixType>::compute()
template<typename MatrixType, typename PermutationIndex>
void PartialPivLU<MatrixType, PermutationIndex>::compute()
{
// the row permutation is stored as int indices, so just to be sure:
eigen_assert(m_lu.rows()<NumTraits<int>::highest());
eigen_assert(m_lu.rows()<NumTraits<PermutationIndex>::highest());
if(m_lu.cols()>0)
m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff();
@ -545,8 +545,8 @@ void PartialPivLU<MatrixType>::compute()
m_isInitialized = true;
}
template<typename MatrixType>
typename PartialPivLU<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() const
template<typename MatrixType, typename PermutationIndex>
typename PartialPivLU<MatrixType, PermutationIndex>::Scalar PartialPivLU<MatrixType, PermutationIndex>::determinant() const
{
eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
return Scalar(m_det_p) * m_lu.diagonal().prod();
@ -555,8 +555,8 @@ typename PartialPivLU<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant(
/** \returns the matrix represented by the decomposition,
* i.e., it returns the product: P^{-1} L U.
* This function is provided for debug purpose. */
template<typename MatrixType>
MatrixType PartialPivLU<MatrixType>::reconstructedMatrix() const
template<typename MatrixType, typename PermutationIndex>
MatrixType PartialPivLU<MatrixType, PermutationIndex>::reconstructedMatrix() const
{
eigen_assert(m_isInitialized && "LU is not initialized.");
// LU
@ -574,10 +574,10 @@ MatrixType PartialPivLU<MatrixType>::reconstructedMatrix() const
namespace internal {
/***** Implementation of inverse() *****************************************************/
template<typename DstXprType, typename MatrixType>
struct Assignment<DstXprType, Inverse<PartialPivLU<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename PartialPivLU<MatrixType>::Scalar>, Dense2Dense>
template<typename DstXprType, typename MatrixType, typename PermutationIndex>
struct Assignment<DstXprType, Inverse<PartialPivLU<MatrixType, PermutationIndex> >, internal::assign_op<typename DstXprType::Scalar,typename PartialPivLU<MatrixType, PermutationIndex>::Scalar>, Dense2Dense>
{
typedef PartialPivLU<MatrixType> LuType;
typedef PartialPivLU<MatrixType, PermutationIndex> LuType;
typedef Inverse<LuType> SrcXprType;
static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename LuType::Scalar> &)
{
@ -595,10 +595,11 @@ struct Assignment<DstXprType, Inverse<PartialPivLU<MatrixType> >, internal::assi
* \sa class PartialPivLU
*/
template<typename Derived>
inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
template<typename PermutationIndex>
inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject, PermutationIndex>
MatrixBase<Derived>::partialPivLu() const
{
return PartialPivLU<PlainObject>(eval());
return PartialPivLU<PlainObject, PermutationIndex>(eval());
}
/** \lu_module
@ -610,10 +611,11 @@ MatrixBase<Derived>::partialPivLu() const
* \sa class PartialPivLU
*/
template<typename Derived>
inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
template<typename PermutationIndex>
inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject, PermutationIndex>
MatrixBase<Derived>::lu() const
{
return PartialPivLU<PlainObject>(eval());
return PartialPivLU<PlainObject, PermutationIndex>(eval());
}
} // end namespace Eigen

3
test/lu.cpp
View File

@ -19,7 +19,6 @@ typename MatrixType::RealScalar matrix_l1_norm(const MatrixType& m) {
template<typename MatrixType> void lu_non_invertible()
{
STATIC_CHECK(( internal::is_same<typename FullPivLU<MatrixType>::StorageIndex,int>::value ));
typedef typename MatrixType::RealScalar RealScalar;
/* this test covers the following files:
@ -163,8 +162,6 @@ template<typename MatrixType> void lu_partial_piv(Index size = MatrixType::ColsA
m1.setRandom();
PartialPivLU<MatrixType> plu(m1);
STATIC_CHECK(( internal::is_same<typename PartialPivLU<MatrixType>::StorageIndex,int>::value ));
VERIFY_IS_APPROX(m1, plu.reconstructedMatrix());
check_solverbase<MatrixType, MatrixType>(m1, plu, size, size, size);