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

- add _packetCoeff() to Inverse, allowing vectorization.

Browse Source 
- let Inverse take template parameter MatrixType instead
  of ExpressionType, in order to reduce executable code size
  when taking inverses of xpr's.
- introduce ei_corrected_matrix_flags : the flags template
  parameter to the Matrix class is only a suggestion. This
  is also useful in ei_eval.
This commit is contained in:
Benoit Jacob 2008-04-16 07:18:27 +00:00
parent 43e2bc14fe
commit acfd6f3bda
5 changed files with 105 additions and 87 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
Eigen/src/Core/Matrix.h
View File

@ -35,7 +35,7 @@
* specify that the number of rows is dynamic, i.e. is not fixed at compile-time. * specify that the number of rows is dynamic, i.e. is not fixed at compile-time.
* \param _Cols the number of columns at compile-time. Use the special value \a Dynamic to * \param _Cols the number of columns at compile-time. Use the special value \a Dynamic to
* specify that the number of columns is dynamic, i.e. is not fixed at compile-time. * specify that the number of columns is dynamic, i.e. is not fixed at compile-time.
* \param _Flags allows to control certain features such as storage order. See MatrixBase::Flags. * \param _SuggestedFlags allows to control certain features such as storage order. See MatrixBase::Flags.
* *
* This single class template covers all kinds of matrix and vectors that Eigen can handle. * This single class template covers all kinds of matrix and vectors that Eigen can handle.
* All matrix and vector types are just typedefs to specializations of this class template. * All matrix and vector types are just typedefs to specializations of this class template.
@ -70,8 +70,8 @@
* *
* Note that most of the API is in the base class MatrixBase. * Note that most of the API is in the base class MatrixBase.
*/ */
template<typename _Scalar, int _Rows, int _Cols, unsigned int _Flags, int _MaxRows, int _MaxCols> template<typename _Scalar, int _Rows, int _Cols, unsigned int _SuggestedFlags, int _MaxRows, int _MaxCols>
struct ei_traits<Matrix<_Scalar, _Rows, _Cols, _Flags, _MaxRows, _MaxCols> > struct ei_traits<Matrix<_Scalar, _Rows, _Cols, _SuggestedFlags, _MaxRows, _MaxCols> >
{ {
typedef _Scalar Scalar; typedef _Scalar Scalar;
enum { enum {
@ -79,11 +79,7 @@ struct ei_traits<Matrix<_Scalar, _Rows, _Cols, _Flags, _MaxRows, _MaxCols> >
ColsAtCompileTime = _Cols, ColsAtCompileTime = _Cols,
MaxRowsAtCompileTime = _MaxRows, MaxRowsAtCompileTime = _MaxRows,
MaxColsAtCompileTime = _MaxCols, MaxColsAtCompileTime = _MaxCols,
Flags = (_Flags & ~VectorizableBit) Flags = ei_corrected_matrix_flags<_Scalar, _Rows, _Cols, _SuggestedFlags>::ret,
| (
ei_is_matrix_vectorizable<Scalar, _Rows, _Cols, _Flags>::ret
? VectorizableBit : 0
),
CoeffReadCost = NumTraits<Scalar>::ReadCost CoeffReadCost = NumTraits<Scalar>::ReadCost
}; };
}; };

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

@ -511,8 +511,8 @@ template<typename Derived> class MatrixBase
* \code #include <Eigen/LU> \endcode * \code #include <Eigen/LU> \endcode
*/ */
//@{ //@{
const Inverse<Derived, true> inverse() const; const Inverse<typename ei_eval<Derived>::type, true> inverse() const;
const Inverse<Derived, false> quickInverse() const; const Inverse<typename ei_eval<Derived>::type, false> quickInverse() const;
Scalar determinant() const; Scalar determinant() const;
//@} //@}

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

@ -30,7 +30,7 @@ template<typename Lhs, typename Rhs> struct ei_product_eval_mode;
template<typename T> struct NumTraits; template<typename T> struct NumTraits;
template<typename _Scalar, int _Rows, int _Cols, template<typename _Scalar, int _Rows, int _Cols,
unsigned int _Flags = EIGEN_DEFAULT_MATRIX_FLAGS, unsigned int _SuggestedFlags = EIGEN_DEFAULT_MATRIX_FLAGS,
int _MaxRows = _Rows, int _MaxCols = _Cols> int _MaxRows = _Rows, int _MaxCols = _Cols>
class Matrix; class Matrix;

60
Eigen/src/Core/util/Meta.h
View File

@ -145,36 +145,44 @@ template<typename T> struct ei_packet_traits
enum {size=1}; enum {size=1};
}; };
template<typename Scalar, int Rows, int Cols, unsigned int Flags> template<typename Scalar, int Rows, int Cols, unsigned int SuggestedFlags>
struct ei_is_matrix_vectorizable class ei_corrected_matrix_flags
{ {
enum { ret = ei_packet_traits<Scalar>::size > 1 enum { is_vectorizable
&& Rows!=Dynamic = ei_packet_traits<Scalar>::size > 1
&& Cols!=Dynamic && Rows!=Dynamic
&& && Cols!=Dynamic
( &&
(Flags&RowMajorBit && Cols%ei_packet_traits<Scalar>::size==0) (
|| (Rows%ei_packet_traits<Scalar>::size==0) SuggestedFlags&RowMajorBit
) ? Cols%ei_packet_traits<Scalar>::size==0
}; : Rows%ei_packet_traits<Scalar>::size==0
),
_flags1 = SuggestedFlags & ~(EvalBeforeNestingBit | EvalBeforeAssigningBit)
};
public:
enum { ret = is_vectorizable
? _flags1 | VectorizableBit
: _flags1 & ~VectorizableBit
};
}; };
template<typename T> struct ei_eval template<typename T> class ei_eval
{ {
typedef typename ei_traits<T>::Scalar _Scalar; typedef typename ei_traits<T>::Scalar _Scalar;
enum { _Rows = ei_traits<T>::RowsAtCompileTime, enum { _Rows = ei_traits<T>::RowsAtCompileTime,
_Cols = ei_traits<T>::ColsAtCompileTime, _Cols = ei_traits<T>::ColsAtCompileTime,
_Flags = ei_traits<T>::Flags _Flags = ei_traits<T>::Flags
}; };
typedef Matrix<_Scalar,
_Rows, public:
_Cols, typedef Matrix<_Scalar,
(_Flags & ~(EvalBeforeNestingBit | EvalBeforeAssigningBit)) _Rows,
| _Cols,
(ei_is_matrix_vectorizable<_Scalar, _Rows, _Cols, _Flags>::ret ei_corrected_matrix_flags<_Scalar, _Rows, _Cols, _Flags>::ret,
? VectorizableBit : 0), ei_traits<T>::MaxRowsAtCompileTime,
ei_traits<T>::MaxRowsAtCompileTime, ei_traits<T>::MaxColsAtCompileTime> type;
ei_traits<T>::MaxColsAtCompileTime> type;
}; };
template<typename T> struct ei_unref { typedef T type; }; template<typename T> struct ei_unref { typedef T type; };

114
Eigen/src/LU/Inverse.h
View File

@ -29,7 +29,7 @@
* *
* \brief Inverse of a matrix * \brief Inverse of a matrix
* *
* \param ExpressionType the type of the matrix/expression of which we are taking the inverse * \param MatrixType the type of the matrix of which we are taking the inverse
* \param CheckExistence whether or not to check the existence of the inverse while computing it * \param CheckExistence whether or not to check the existence of the inverse while computing it
* *
* This class represents the inverse of a matrix. It is the return * This class represents the inverse of a matrix. It is the return
@ -38,11 +38,10 @@
* *
* \sa MatrixBase::inverse(), MatrixBase::quickInverse() * \sa MatrixBase::inverse(), MatrixBase::quickInverse()
*/ */
template<typename ExpressionType, bool CheckExistence> template<typename MatrixType, bool CheckExistence>
struct ei_traits<Inverse<ExpressionType, CheckExistence> > struct ei_traits<Inverse<MatrixType, CheckExistence> >
{ {
typedef typename ExpressionType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename ExpressionType::Eval MatrixType;
enum { enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime, RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime,
@ -53,20 +52,19 @@ struct ei_traits<Inverse<ExpressionType, CheckExistence> >
}; };
}; };
template<typename ExpressionType, bool CheckExistence> class Inverse : ei_no_assignment_operator, template<typename MatrixType, bool CheckExistence> class Inverse : ei_no_assignment_operator,
public MatrixBase<Inverse<ExpressionType, CheckExistence> > public MatrixBase<Inverse<MatrixType, CheckExistence> >
{ {
public: public:
EIGEN_GENERIC_PUBLIC_INTERFACE(Inverse) EIGEN_GENERIC_PUBLIC_INTERFACE(Inverse)
typedef typename ei_traits<Inverse>::MatrixType MatrixType;
Inverse(const ExpressionType& xpr) Inverse(const MatrixType& matrix)
: m_exists(true), : m_exists(true),
m_inverse(MatrixType::identity(xpr.rows(), xpr.cols())) m_inverse(MatrixType::identity(matrix.rows(), matrix.cols()))
{ {
ei_assert(xpr.rows() == xpr.cols()); ei_assert(matrix.rows() == matrix.cols());
_compute(xpr); _compute(matrix);
} }
/** \returns whether or not the inverse exists. /** \returns whether or not the inverse exists.
@ -86,24 +84,29 @@ template<typename ExpressionType, bool CheckExistence> class Inverse : ei_no_ass
return m_inverse.coeff(row, col); return m_inverse.coeff(row, col);
} }
PacketScalar _packetCoeff(int row, int col) const
{
return m_inverse.packetCoeff(row, col);
}
enum { _Size = MatrixType::RowsAtCompileTime }; enum { _Size = MatrixType::RowsAtCompileTime };
void _compute(const ExpressionType& xpr); void _compute(const MatrixType& matrix);
void _compute_in_general_case(const ExpressionType& xpr); void _compute_in_general_case(const MatrixType& matrix);
void _compute_in_size1_case(const ExpressionType& xpr); void _compute_in_size1_case(const MatrixType& matrix);
void _compute_in_size2_case(const ExpressionType& xpr); void _compute_in_size2_case(const MatrixType& matrix);
void _compute_in_size3_case(const ExpressionType& xpr); void _compute_in_size3_case(const MatrixType& matrix);
void _compute_in_size4_case(const ExpressionType& xpr); void _compute_in_size4_case(const MatrixType& matrix);
protected: protected:
bool m_exists; bool m_exists;
MatrixType m_inverse; MatrixType m_inverse;
}; };
template<typename ExpressionType, bool CheckExistence> template<typename MatrixType, bool CheckExistence>
void Inverse<ExpressionType, CheckExistence> void Inverse<MatrixType, CheckExistence>
::_compute_in_general_case(const ExpressionType& xpr) ::_compute_in_general_case(const MatrixType& _matrix)
{ {
MatrixType matrix(xpr); MatrixType matrix(_matrix);
const RealScalar max = CheckExistence ? matrix.cwiseAbs().maxCoeff() const RealScalar max = CheckExistence ? matrix.cwiseAbs().maxCoeff()
: static_cast<RealScalar>(0); : static_cast<RealScalar>(0);
const int size = matrix.rows(); const int size = matrix.rows();
@ -141,10 +144,10 @@ void Inverse<ExpressionType, CheckExistence>
} }
} }
template<typename ExpressionType, typename MatrixType, bool CheckExistence> template<typename ExpressionType, bool CheckExistence>
bool ei_compute_size2_inverse(const ExpressionType& xpr, MatrixType* result) bool ei_compute_size2_inverse(const ExpressionType& xpr, typename ExpressionType::Eval* result)
{ {
typedef typename MatrixType::Scalar Scalar; typedef typename ExpressionType::Scalar Scalar;
const typename ei_nested<ExpressionType, 1+CheckExistence>::type matrix(xpr); const typename ei_nested<ExpressionType, 1+CheckExistence>::type matrix(xpr);
const Scalar det = matrix.determinant(); const Scalar det = matrix.determinant();
if(CheckExistence && ei_isMuchSmallerThan(det, matrix.cwiseAbs().maxCoeff())) if(CheckExistence && ei_isMuchSmallerThan(det, matrix.cwiseAbs().maxCoeff()))
@ -157,10 +160,9 @@ bool ei_compute_size2_inverse(const ExpressionType& xpr, MatrixType* result)
return true; return true;
} }
template<typename ExpressionType, bool CheckExistence> template<typename MatrixType, bool CheckExistence>
void Inverse<ExpressionType, CheckExistence>::_compute_in_size3_case(const ExpressionType& xpr) void Inverse<MatrixType, CheckExistence>::_compute_in_size3_case(const MatrixType& matrix)
{ {
const typename ei_nested<ExpressionType, 2+CheckExistence>::type matrix(xpr);
const Scalar det_minor00 = matrix.minor(0,0).determinant(); const Scalar det_minor00 = matrix.minor(0,0).determinant();
const Scalar det_minor10 = matrix.minor(1,0).determinant(); const Scalar det_minor10 = matrix.minor(1,0).determinant();
const Scalar det_minor20 = matrix.minor(2,0).determinant(); const Scalar det_minor20 = matrix.minor(2,0).determinant();
@ -184,25 +186,37 @@ void Inverse<ExpressionType, CheckExistence>::_compute_in_size3_case(const Expre
} }
} }
template<typename ExpressionType, bool CheckExistence> template<typename MatrixType, bool CheckExistence>
void Inverse<ExpressionType, CheckExistence>::_compute_in_size4_case(const ExpressionType& xpr) void Inverse<MatrixType, CheckExistence>::_compute_in_size4_case(const MatrixType& matrix)
{ {
typedef Block<ExpressionType,2,2> XprBlock22; /* Let's split M into four 2x2 blocks:
* (P Q)
* (R S)
* If P is invertible, with inverse denoted by P_inverse, and if
* (S - R*P_inverse*Q) is also invertible, then the inverse of M is
* (P' Q')
* (R' S')
* where
* S' = (S - R*P_inverse*Q)^(-1)
* P' = P1 + (P1*Q) * S' *(R*P_inverse)
* Q' = -(P_inverse*Q) * S'
* R' = -S' * (R*P_inverse)
*/
typedef Block<MatrixType,2,2> XprBlock22;
typedef typename XprBlock22::Eval Block22; typedef typename XprBlock22::Eval Block22;
Block22 P_inverse; Block22 P_inverse;
if(ei_compute_size2_inverse<XprBlock22, Block22, true>(xpr.template block<2,2>(0,0), &P_inverse)) if(ei_compute_size2_inverse<XprBlock22, true>(matrix.template block<2,2>(0,0), &P_inverse))
{ {
const Block22 Q = xpr.template block<2,2>(0,2); const Block22 Q = matrix.template block<2,2>(0,2);
const Block22 P_inverse_times_Q = P_inverse * Q; const Block22 P_inverse_times_Q = P_inverse * Q;
const XprBlock22 R = xpr.template block<2,2>(2,0); const XprBlock22 R = matrix.template block<2,2>(2,0);
const Block22 R_times_P_inverse = R * P_inverse; const Block22 R_times_P_inverse = R * P_inverse;
const Block22 R_times_P_inverse_times_Q = R_times_P_inverse * Q; const Block22 R_times_P_inverse_times_Q = R_times_P_inverse * Q;
const XprBlock22 S = xpr.template block<2,2>(2,2); const XprBlock22 S = matrix.template block<2,2>(2,2);
const Block22 X = S - R_times_P_inverse_times_Q; const Block22 X = S - R_times_P_inverse_times_Q;
Block22 Y; Block22 Y;
if(ei_compute_size2_inverse<Block22, Block22, CheckExistence>(X, &Y)) if(ei_compute_size2_inverse<Block22, CheckExistence>(X, &Y))
{ {
m_inverse.template block<2,2>(2,2) = Y; m_inverse.template block<2,2>(2,2) = Y;
m_inverse.template block<2,2>(2,0) = - Y * R_times_P_inverse; m_inverse.template block<2,2>(2,0) = - Y * R_times_P_inverse;
@ -217,16 +231,16 @@ void Inverse<ExpressionType, CheckExistence>::_compute_in_size4_case(const Expre
} }
else else
{ {
_compute_in_general_case(xpr); _compute_in_general_case(matrix);
} }
} }
template<typename ExpressionType, bool CheckExistence> template<typename MatrixType, bool CheckExistence>
void Inverse<ExpressionType, CheckExistence>::_compute(const ExpressionType& xpr) void Inverse<MatrixType, CheckExistence>::_compute(const MatrixType& matrix)
{ {
if(_Size == 1) if(_Size == 1)
{ {
const Scalar x = xpr.coeff(0,0); const Scalar x = matrix.coeff(0,0);
if(CheckExistence && x == static_cast<Scalar>(0)) if(CheckExistence && x == static_cast<Scalar>(0))
m_exists = false; m_exists = false;
else else
@ -235,13 +249,13 @@ void Inverse<ExpressionType, CheckExistence>::_compute(const ExpressionType& xpr
else if(_Size == 2) else if(_Size == 2)
{ {
if(CheckExistence) if(CheckExistence)
m_exists = ei_compute_size2_inverse<ExpressionType, MatrixType, true>(xpr, &m_inverse); m_exists = ei_compute_size2_inverse<MatrixType, true>(matrix, &m_inverse);
else else
ei_compute_size2_inverse<ExpressionType, MatrixType, false>(xpr, &m_inverse); ei_compute_size2_inverse<MatrixType, false>(matrix, &m_inverse);
} }
else if(_Size == 3) _compute_in_size3_case(xpr); else if(_Size == 3) _compute_in_size3_case(matrix);
else if(_Size == 4) _compute_in_size4_case(xpr); else if(_Size == 4) _compute_in_size4_case(matrix);
else _compute_in_general_case(xpr); else _compute_in_general_case(matrix);
} }
/** \return the matrix inverse of \c *this, if it exists. /** \return the matrix inverse of \c *this, if it exists.
@ -252,10 +266,10 @@ void Inverse<ExpressionType, CheckExistence>::_compute(const ExpressionType& xpr
* \sa class Inverse * \sa class Inverse
*/ */
template<typename Derived> template<typename Derived>
const Inverse<Derived, true> const Inverse<typename ei_eval<Derived>::type, true>
MatrixBase<Derived>::inverse() const MatrixBase<Derived>::inverse() const
{ {
return Inverse<Derived, true>(derived()); return Inverse<typename ei_eval<Derived>::type, true>(derived());
} }
/** \return the matrix inverse of \c *this, which is assumed to exist. /** \return the matrix inverse of \c *this, which is assumed to exist.
@ -266,10 +280,10 @@ MatrixBase<Derived>::inverse() const
* \sa class Inverse * \sa class Inverse
*/ */
template<typename Derived> template<typename Derived>
const Inverse<Derived, false> const Inverse<typename ei_eval<Derived>::type, false>
MatrixBase<Derived>::quickInverse() const MatrixBase<Derived>::quickInverse() const
{ {
return Inverse<Derived, false>(derived()); return Inverse<typename ei_eval<Derived>::type, false>(derived());
} }
#endif // EIGEN_INVERSE_H #endif // EIGEN_INVERSE_H