PermutationMatrix:

* make multiplication order not be reversed
 * release-quality documentation
This commit is contained in:
Benoit Jacob 2009-11-16 15:07:33 -05:00
parent 8a1bada43d
commit e8d0dbf82e
3 changed files with 60 additions and 17 deletions

View File

@ -25,18 +25,24 @@
#ifndef EIGEN_PERMUTATIONMATRIX_H
#define EIGEN_PERMUTATIONMATRIX_H
/** \nonstableyet
* \class PermutationMatrix
/** \class PermutationMatrix
*
* \brief Permutation matrix
*
* \param SizeAtCompileTime the number of rows/cols, or Dynamic
* \param MaxSizeAtCompileTime the maximum number of rows/cols, or Dynamic. This optional parameter defaults to SizeAtCompileTime.
* \param MaxSizeAtCompileTime the maximum number of rows/cols, or Dynamic. This optional parameter defaults to SizeAtCompileTime. Most of the time, you should not have to specify it.
*
* This class represents a permutation matrix, internally stored as a vector of integers.
* The convention followed here is the same as on <a href="http://en.wikipedia.org/wiki/Permutation_matrix">Wikipedia</a>,
* namely: the matrix of permutation \a p is the matrix such that on each row \a i, the only nonzero coefficient is
* in column p(i).
* The convention followed here is that if \f$ \sigma \f$ is a permutation, the corresponding permutation matrix
* \f$ P_\sigma \f$ is such that if \f$ (e_1,\ldots,e_p) \f$ is the canonical basis, we have:
* \f[ P_\sigma(e_i) = e_{\sigma(i)}. \f]
* This convention ensures that for any two permutations \f$ \sigma, \tau \f$, we have:
* \f[ P_{\sigma\circ\tau} = P_\sigma P_\tau. \f]
*
* Permutation matrices are square and invertible.
*
* Notice that in addition to the member functions and operators listed here, there also are non-member
* operator* to multiply a PermutationMatrix with any kind of matrix expression (MatrixBase) on either side.
*
* \sa class DiagonalMatrix
*/
@ -53,6 +59,7 @@ class PermutationMatrix : public AnyMatrixBase<PermutationMatrix<SizeAtCompileTi
{
public:
#ifndef EIGEN_PARSED_BY_DOXYGEN
typedef ei_traits<PermutationMatrix> Traits;
typedef Matrix<int,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime>
DenseMatrixType;
@ -65,25 +72,37 @@ class PermutationMatrix : public AnyMatrixBase<PermutationMatrix<SizeAtCompileTi
MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
};
typedef typename Traits::Scalar Scalar;
#endif
typedef Matrix<int, RowsAtCompileTime, 1, 0, MaxRowsAtCompileTime, 1> IndicesType;
typedef Matrix<int, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
inline PermutationMatrix()
{
}
/** Copy constructor. */
template<int OtherSize, int OtherMaxSize>
inline PermutationMatrix(const PermutationMatrix<OtherSize, OtherMaxSize>& other)
: m_indices(other.indices()) {}
/** copy constructor. prevent a default copy constructor from hiding the other templated constructor */
#ifndef EIGEN_PARSED_BY_DOXYGEN
/** Standard copy constructor. Defined only to prevent a default copy constructor
* from hiding the other templated constructor */
inline PermutationMatrix(const PermutationMatrix& other) : m_indices(other.indices()) {}
#endif
/** generic constructor from expression of the indices */
/** Generic constructor from expression of the indices. The indices
* array has the meaning that the permutations sends each integer i to indices[i].
*
* \warning It is your responsibility to check that the indices array that you passes actually
* describes a permutation, \ie each value between 0 and n-1 occurs exactly once, where n is the
* array's size.
*/
template<typename Other>
explicit inline PermutationMatrix(const MatrixBase<Other>& other) : m_indices(other)
explicit inline PermutationMatrix(const MatrixBase<Other>& indices) : m_indices(indices)
{}
/** Copies the other permutation into *this */
template<int OtherSize, int OtherMaxSize>
PermutationMatrix& operator=(const PermutationMatrix<OtherSize, OtherMaxSize>& other)
{
@ -91,6 +110,7 @@ class PermutationMatrix : public AnyMatrixBase<PermutationMatrix<SizeAtCompileTi
return *this;
}
#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=.
*/
@ -99,8 +119,12 @@ class PermutationMatrix : public AnyMatrixBase<PermutationMatrix<SizeAtCompileTi
m_indices = other.m_indices();
return *this;
}
#endif
inline PermutationMatrix(int rows, int cols) : m_indices(rows)
/** Constructs an uninitialized permutation matrix of given size. Note that it is required
* that rows == cols, since permutation matrices are square. The \a cols parameter may be omitted.
*/
inline PermutationMatrix(int rows, int cols = rows) : m_indices(rows)
{
ei_assert(rows == cols);
}
@ -111,24 +135,33 @@ class PermutationMatrix : public AnyMatrixBase<PermutationMatrix<SizeAtCompileTi
/** \returns the number of columns */
inline int cols() const { return m_indices.size(); }
#ifndef EIGEN_PARSED_BY_DOXYGEN
template<typename DenseDerived>
void evalTo(MatrixBase<DenseDerived>& other) const
{
other.setZero();
for (int i=0; i<rows();++i)
other.coeffRef(i,m_indices.coeff(i)) = typename DenseDerived::Scalar(1);
other.coeffRef(m_indices.coeff(i),i) = typename DenseDerived::Scalar(1);
}
#endif
/** \returns a Matrix object initialized from this permutation matrix. Notice that it
* is inefficient to return this Matrix object by value. For efficiency, favor using
* the Matrix constructor taking AnyMatrixBase objects.
*/
DenseMatrixType toDenseMatrix() const
{
return *this;
}
/** const version of indices(). */
const IndicesType& indices() const { return m_indices; }
/** \returns a reference to the stored array representing the permutation. */
IndicesType& indices() { return m_indices; }
/**** inversion and multiplication helpers to hopefully get RVO ****/
#ifndef EIGEN_PARSED_BY_DOXYGEN
protected:
enum Inverse_t {Inverse};
PermutationMatrix(Inverse_t, const PermutationMatrix& other)
@ -143,10 +176,19 @@ class PermutationMatrix : public AnyMatrixBase<PermutationMatrix<SizeAtCompileTi
ei_assert(lhs.cols() == rhs.rows());
for (int i=0; i<rows();++i) m_indices.coeffRef(i) = lhs.m_indices.coeff(rhs.m_indices.coeff(i));
}
#endif
public:
/** \returns the inverse permutation matrix.
*
* \note \note_try_to_help_rvo
*/
inline PermutationMatrix inverse() const
{ return PermutationMatrix(Inverse, *this); }
/** \returns the product permutation matrix.
*
* \note \note_try_to_help_rvo
*/
template<int OtherSize, int OtherMaxSize>
inline PermutationMatrix operator*(const PermutationMatrix<OtherSize, OtherMaxSize>& other) const
{ return PermutationMatrix(Product, *this, other); }
@ -209,7 +251,7 @@ struct ei_permut_matrix_product_retval
Dest,
Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime,
Side==OnTheRight ? 1 : Dest::ColsAtCompileTime
>(dst, Side==OnTheRight ? m_permutation.indices().coeff(i) : i)
>(dst, Side==OnTheLeft ? m_permutation.indices().coeff(i) : i)
=
@ -217,7 +259,7 @@ struct ei_permut_matrix_product_retval
MatrixTypeNestedCleaned,
Side==OnTheLeft ? 1 : MatrixType::RowsAtCompileTime,
Side==OnTheRight ? 1 : MatrixType::ColsAtCompileTime
>(m_matrix, Side==OnTheLeft ? m_permutation.indices().coeff(i) : i);
>(m_matrix, Side==OnTheRight ? m_permutation.indices().coeff(i) : i);
}
}

View File

@ -218,7 +218,8 @@ ALIASES = "only_for_vectors=This is only for vectors (either row-
"nonstableyet=\warning This is not considered to be part of the stable public API yet. Changes may happen in future releases. See \ref Experimental \"Experimental parts of Eigen\"" \
"note_about_arbitrary_choice_of_solution=If there exists more than one solution, this method will arbitrarily choose one." \
"note_about_using_kernel_to_study_multiple_solutions=If you need a complete analysis of the space of solutions, take the one solution obtained by this method and add to it elements of the kernel, as determined by kernel()." \
"note_about_checking_solutions=This method just tries to find as good a solution as possible. If you want to check whether a solution exists or if it is accurate, just call this function to get a result and then compute the error of this result, or use MatrixBase::isApprox() directly, for instance like this: \code bool a_solution_exists = (A*result).isApprox(b, precision); \endcode This method avoids dividing by zero, so that the non-existence of a solution doesn't by itself mean that you'll get \c inf or \c nan values."
"note_about_checking_solutions=This method just tries to find as good a solution as possible. If you want to check whether a solution exists or if it is accurate, just call this function to get a result and then compute the error of this result, or use MatrixBase::isApprox() directly, for instance like this: \code bool a_solution_exists = (A*result).isApprox(b, precision); \endcode This method avoids dividing by zero, so that the non-existence of a solution doesn't by itself mean that you'll get \c inf or \c nan values." \
"note_try_to_help_rvo=This function returns the result by value. In order to make that efficient, it is implemented as just a return statement using a special constructor, hopefully allowing the compiler to perform a RVO (return value optimization)."
# Set the OPTIMIZE_OUTPUT_FOR_C tag to YES if your project consists of C
# sources only. Doxygen will then generate output that is more tailored for C.

View File

@ -66,7 +66,7 @@ template<typename MatrixType> void permutationmatrices(const MatrixType& m)
for (int i=0; i<rows; i++)
for (int j=0; j<cols; j++)
VERIFY_IS_APPROX(m_original(lv(i),j), m_permuted(i,rv(j)));
VERIFY_IS_APPROX(m_permuted(lv(i),j), m_original(i,rv(j)));
Matrix<Scalar,Rows,Rows> lm(lp);
Matrix<Scalar,Cols,Cols> rm(rp);
@ -80,7 +80,7 @@ template<typename MatrixType> void permutationmatrices(const MatrixType& m)
randomPermutationVector(lv2, rows);
LeftPermutationType lp2(lv2);
Matrix<Scalar,Rows,Rows> lm2(lp2);
VERIFY_IS_APPROX((lp*lp2).toDenseMatrix().template cast<Scalar>(), lm2*lm);
VERIFY_IS_APPROX((lp*lp2).toDenseMatrix().template cast<Scalar>(), lm*lm2);
}
void test_permutationmatrices()