* add Transpositions to PermutationMatrix conversion

* make PartialPivLu uses the  Transpositions class
This commit is contained in:
Gael Guennebaud 2010-06-08 22:23:11 +02:00
parent 684656d41c
commit 50e43bc75a
4 changed files with 38 additions and 16 deletions

View File

@ -47,7 +47,6 @@
*
* \sa class DiagonalMatrix
*/
template<int SizeAtCompileTime, int MaxSizeAtCompileTime = SizeAtCompileTime> class PermutationMatrix;
template<typename PermutationType, typename MatrixType, int Side, bool Transposed=false> struct ei_permut_matrix_product_retval;
template<int SizeAtCompileTime, int MaxSizeAtCompileTime>
@ -78,8 +77,12 @@ class PermutationMatrix : public EigenBase<PermutationMatrix<SizeAtCompileTime,
typedef Matrix<int, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
inline PermutationMatrix()
{
}
{}
/** Constructs an uninitialized permutation matrix of given size.
*/
inline PermutationMatrix(int size) : m_indices(size)
{}
/** Copy constructor. */
template<int OtherSize, int OtherMaxSize>
@ -103,6 +106,14 @@ class PermutationMatrix : public EigenBase<PermutationMatrix<SizeAtCompileTime,
explicit inline PermutationMatrix(const MatrixBase<Other>& indices) : m_indices(indices)
{}
/** Convert the Transpositions \a tr to a permutation matrix */
template<int OtherSize, int OtherMaxSize>
explicit PermutationMatrix(const Transpositions<OtherSize,OtherMaxSize>& tr)
: m_indices(tr.size())
{
*this = tr;
}
/** Copies the other permutation into *this */
template<int OtherSize, int OtherMaxSize>
PermutationMatrix& operator=(const PermutationMatrix<OtherSize, OtherMaxSize>& other)
@ -111,6 +122,15 @@ class PermutationMatrix : public EigenBase<PermutationMatrix<SizeAtCompileTime,
return *this;
}
/** Assignment from the Transpositions \a tr */
template<int OtherSize, int OtherMaxSize>
PermutationMatrix& operator=(const Transpositions<OtherSize,OtherMaxSize>& tr)
{
setIdentity(tr.size());
for(int k=size()-1; k>=0; --k)
applyTranspositionOnTheRight(k,tr.coeff(k));
}
#ifndef EIGEN_PARSED_BY_DOXYGEN
/** This is a special case of the templated operator=. Its purpose is to
* prevent a default operator= from hiding the templated operator=.
@ -122,11 +142,6 @@ class PermutationMatrix : public EigenBase<PermutationMatrix<SizeAtCompileTime,
}
#endif
/** Constructs an uninitialized permutation matrix of given size.
*/
inline PermutationMatrix(int size) : m_indices(size)
{}
/** \returns the number of rows */
inline int rows() const { return m_indices.size(); }

View File

@ -52,7 +52,6 @@
*
* \sa class PermutationMatrix
*/
template<int SizeAtCompileTime, int MaxSizeAtCompileTime = SizeAtCompileTime> class Transpositions;
template<typename TranspositionType, typename MatrixType, int Side, bool Transposed=false> struct ei_transposition_matrix_product_retval;
template<int SizeAtCompileTime, int MaxSizeAtCompileTime>
@ -108,10 +107,18 @@ class Transpositions
/** \returns the number of transpositions */
inline Index size() const { return m_indices.size(); }
/** Direct access to the underlying index vector */
inline const Index& coeff(Index i) const { return m_indices.coeff(i); }
/** Direct access to the underlying index vector */
inline Index& coeffRef(Index i) { return m_indices.coeffRef(i); }
/** Direct access to the underlying index vector */
inline const Index& operator()(Index i) const { return m_indices(i); }
/** Direct access to the underlying index vector */
inline Index& operator()(Index i) { return m_indices(i); }
/** Direct access to the underlying index vector */
inline const Index& operator[](Index i) const { return m_indices(i); }
/** Direct access to the underlying index vector */
inline Index& operator[](Index i) { return m_indices(i); }
/** const version of indices(). */
const IndicesType& indices() const { return m_indices; }

View File

@ -77,6 +77,8 @@ template<typename _DiagonalVectorType> class DiagonalWrapper;
template<typename _Scalar, int SizeAtCompileTime, int MaxSizeAtCompileTime=SizeAtCompileTime> class DiagonalMatrix;
template<typename MatrixType, typename DiagonalType, int ProductOrder> class DiagonalProduct;
template<typename MatrixType, int Index> class Diagonal;
template<int SizeAtCompileTime, int MaxSizeAtCompileTime = SizeAtCompileTime> class PermutationMatrix;
template<int SizeAtCompileTime, int MaxSizeAtCompileTime = SizeAtCompileTime> class Transpositions;
template<int InnerStrideAtCompileTime, int OuterStrideAtCompileTime> class Stride;
template<typename MatrixType, int MapOptions=Unaligned, typename StrideType = Stride<0,0> > class Map;

View File

@ -73,8 +73,8 @@ template<typename _MatrixType> class PartialPivLU
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
typedef typename ei_traits<MatrixType>::StorageKind StorageKind;
typedef typename MatrixType::Index Index;
typedef typename ei_plain_col_type<MatrixType, Index>::type PermutationVectorType;
typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
/**
@ -186,7 +186,7 @@ template<typename _MatrixType> class PartialPivLU
protected:
MatrixType m_lu;
PermutationType m_p;
PermutationVectorType m_rowsTranspositions;
TranspositionType m_rowsTranspositions;
Index m_det_p;
bool m_isInitialized;
};
@ -389,8 +389,8 @@ struct ei_partial_lu_impl
/** \internal performs the LU decomposition with partial pivoting in-place.
*/
template<typename MatrixType, typename IntVector>
void ei_partial_lu_inplace(MatrixType& lu, IntVector& row_transpositions, typename MatrixType::Index& nb_transpositions)
template<typename MatrixType, typename TranspositionType>
void ei_partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, typename MatrixType::Index& nb_transpositions)
{
ei_assert(lu.cols() == row_transpositions.size());
ei_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1);
@ -414,9 +414,7 @@ PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(const MatrixType& ma
ei_partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions);
m_det_p = (nb_transpositions%2) ? -1 : 1;
m_p.setIdentity(size);
for(Index k = size-1; k >= 0; --k)
m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k));
m_p = m_rowsTranspositions;
m_isInitialized = true;
return *this;