diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h index 740d6fc9a..92f726011 100644 --- a/Eigen/src/Core/Matrix.h +++ b/Eigen/src/Core/Matrix.h @@ -35,7 +35,7 @@ * 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 * 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. * 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. */ -template -struct ei_traits > +template +struct ei_traits > { typedef _Scalar Scalar; enum { @@ -79,11 +79,7 @@ struct ei_traits > ColsAtCompileTime = _Cols, MaxRowsAtCompileTime = _MaxRows, MaxColsAtCompileTime = _MaxCols, - Flags = (_Flags & ~VectorizableBit) - | ( - ei_is_matrix_vectorizable::ret - ? VectorizableBit : 0 - ), + Flags = ei_corrected_matrix_flags<_Scalar, _Rows, _Cols, _SuggestedFlags>::ret, CoeffReadCost = NumTraits::ReadCost }; }; diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 23a82051f..d0faa88bf 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -511,8 +511,8 @@ template class MatrixBase * \code #include \endcode */ //@{ - const Inverse inverse() const; - const Inverse quickInverse() const; + const Inverse::type, true> inverse() const; + const Inverse::type, false> quickInverse() const; Scalar determinant() const; //@} diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 86eaf4e54..0bdfb4ea1 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -30,7 +30,7 @@ template struct ei_product_eval_mode; template struct NumTraits; template class Matrix; diff --git a/Eigen/src/Core/util/Meta.h b/Eigen/src/Core/util/Meta.h index 8eafec2f5..ce0c4a21b 100644 --- a/Eigen/src/Core/util/Meta.h +++ b/Eigen/src/Core/util/Meta.h @@ -145,36 +145,44 @@ template struct ei_packet_traits enum {size=1}; }; -template -struct ei_is_matrix_vectorizable +template +class ei_corrected_matrix_flags { - enum { ret = ei_packet_traits::size > 1 - && Rows!=Dynamic - && Cols!=Dynamic - && - ( - (Flags&RowMajorBit && Cols%ei_packet_traits::size==0) - || (Rows%ei_packet_traits::size==0) - ) - }; + enum { is_vectorizable + = ei_packet_traits::size > 1 + && Rows!=Dynamic + && Cols!=Dynamic + && + ( + SuggestedFlags&RowMajorBit + ? Cols%ei_packet_traits::size==0 + : Rows%ei_packet_traits::size==0 + ), + _flags1 = SuggestedFlags & ~(EvalBeforeNestingBit | EvalBeforeAssigningBit) + }; + + public: + enum { ret = is_vectorizable + ? _flags1 | VectorizableBit + : _flags1 & ~VectorizableBit + }; }; -template struct ei_eval +template class ei_eval { - typedef typename ei_traits::Scalar _Scalar; - enum { _Rows = ei_traits::RowsAtCompileTime, - _Cols = ei_traits::ColsAtCompileTime, - _Flags = ei_traits::Flags - }; - typedef Matrix<_Scalar, - _Rows, - _Cols, - (_Flags & ~(EvalBeforeNestingBit | EvalBeforeAssigningBit)) - | - (ei_is_matrix_vectorizable<_Scalar, _Rows, _Cols, _Flags>::ret - ? VectorizableBit : 0), - ei_traits::MaxRowsAtCompileTime, - ei_traits::MaxColsAtCompileTime> type; + typedef typename ei_traits::Scalar _Scalar; + enum { _Rows = ei_traits::RowsAtCompileTime, + _Cols = ei_traits::ColsAtCompileTime, + _Flags = ei_traits::Flags + }; + + public: + typedef Matrix<_Scalar, + _Rows, + _Cols, + ei_corrected_matrix_flags<_Scalar, _Rows, _Cols, _Flags>::ret, + ei_traits::MaxRowsAtCompileTime, + ei_traits::MaxColsAtCompileTime> type; }; template struct ei_unref { typedef T type; }; diff --git a/Eigen/src/LU/Inverse.h b/Eigen/src/LU/Inverse.h index 84d419f40..1d4bd9bf0 100644 --- a/Eigen/src/LU/Inverse.h +++ b/Eigen/src/LU/Inverse.h @@ -29,7 +29,7 @@ * * \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 * * This class represents the inverse of a matrix. It is the return @@ -38,11 +38,10 @@ * * \sa MatrixBase::inverse(), MatrixBase::quickInverse() */ -template -struct ei_traits > +template +struct ei_traits > { - typedef typename ExpressionType::Scalar Scalar; - typedef typename ExpressionType::Eval MatrixType; + typedef typename MatrixType::Scalar Scalar; enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, @@ -53,20 +52,19 @@ struct ei_traits > }; }; -template class Inverse : ei_no_assignment_operator, - public MatrixBase > +template class Inverse : ei_no_assignment_operator, + public MatrixBase > { public: EIGEN_GENERIC_PUBLIC_INTERFACE(Inverse) - typedef typename ei_traits::MatrixType MatrixType; - Inverse(const ExpressionType& xpr) + Inverse(const MatrixType& matrix) : 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()); - _compute(xpr); + ei_assert(matrix.rows() == matrix.cols()); + _compute(matrix); } /** \returns whether or not the inverse exists. @@ -86,24 +84,29 @@ template class Inverse : ei_no_ass return m_inverse.coeff(row, col); } + PacketScalar _packetCoeff(int row, int col) const + { + return m_inverse.packetCoeff(row, col); + } + enum { _Size = MatrixType::RowsAtCompileTime }; - void _compute(const ExpressionType& xpr); - void _compute_in_general_case(const ExpressionType& xpr); - void _compute_in_size1_case(const ExpressionType& xpr); - void _compute_in_size2_case(const ExpressionType& xpr); - void _compute_in_size3_case(const ExpressionType& xpr); - void _compute_in_size4_case(const ExpressionType& xpr); + void _compute(const MatrixType& matrix); + void _compute_in_general_case(const MatrixType& matrix); + void _compute_in_size1_case(const MatrixType& matrix); + void _compute_in_size2_case(const MatrixType& matrix); + void _compute_in_size3_case(const MatrixType& matrix); + void _compute_in_size4_case(const MatrixType& matrix); protected: bool m_exists; MatrixType m_inverse; }; -template -void Inverse -::_compute_in_general_case(const ExpressionType& xpr) +template +void Inverse +::_compute_in_general_case(const MatrixType& _matrix) { - MatrixType matrix(xpr); + MatrixType matrix(_matrix); const RealScalar max = CheckExistence ? matrix.cwiseAbs().maxCoeff() : static_cast(0); const int size = matrix.rows(); @@ -141,10 +144,10 @@ void Inverse } } -template -bool ei_compute_size2_inverse(const ExpressionType& xpr, MatrixType* result) +template +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::type matrix(xpr); const Scalar det = matrix.determinant(); if(CheckExistence && ei_isMuchSmallerThan(det, matrix.cwiseAbs().maxCoeff())) @@ -157,10 +160,9 @@ bool ei_compute_size2_inverse(const ExpressionType& xpr, MatrixType* result) return true; } -template -void Inverse::_compute_in_size3_case(const ExpressionType& xpr) +template +void Inverse::_compute_in_size3_case(const MatrixType& matrix) { - const typename ei_nested::type matrix(xpr); const Scalar det_minor00 = matrix.minor(0,0).determinant(); const Scalar det_minor10 = matrix.minor(1,0).determinant(); const Scalar det_minor20 = matrix.minor(2,0).determinant(); @@ -184,25 +186,37 @@ void Inverse::_compute_in_size3_case(const Expre } } -template -void Inverse::_compute_in_size4_case(const ExpressionType& xpr) +template +void Inverse::_compute_in_size4_case(const MatrixType& matrix) { - typedef Block 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 XprBlock22; typedef typename XprBlock22::Eval Block22; - Block22 P_inverse; - if(ei_compute_size2_inverse(xpr.template block<2,2>(0,0), &P_inverse)) + if(ei_compute_size2_inverse(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 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_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; Block22 Y; - if(ei_compute_size2_inverse(X, &Y)) + if(ei_compute_size2_inverse(X, &Y)) { m_inverse.template block<2,2>(2,2) = Y; m_inverse.template block<2,2>(2,0) = - Y * R_times_P_inverse; @@ -217,16 +231,16 @@ void Inverse::_compute_in_size4_case(const Expre } else { - _compute_in_general_case(xpr); + _compute_in_general_case(matrix); } } -template -void Inverse::_compute(const ExpressionType& xpr) +template +void Inverse::_compute(const MatrixType& matrix) { if(_Size == 1) { - const Scalar x = xpr.coeff(0,0); + const Scalar x = matrix.coeff(0,0); if(CheckExistence && x == static_cast(0)) m_exists = false; else @@ -235,13 +249,13 @@ void Inverse::_compute(const ExpressionType& xpr else if(_Size == 2) { if(CheckExistence) - m_exists = ei_compute_size2_inverse(xpr, &m_inverse); + m_exists = ei_compute_size2_inverse(matrix, &m_inverse); else - ei_compute_size2_inverse(xpr, &m_inverse); + ei_compute_size2_inverse(matrix, &m_inverse); } - else if(_Size == 3) _compute_in_size3_case(xpr); - else if(_Size == 4) _compute_in_size4_case(xpr); - else _compute_in_general_case(xpr); + else if(_Size == 3) _compute_in_size3_case(matrix); + else if(_Size == 4) _compute_in_size4_case(matrix); + else _compute_in_general_case(matrix); } /** \return the matrix inverse of \c *this, if it exists. @@ -252,10 +266,10 @@ void Inverse::_compute(const ExpressionType& xpr * \sa class Inverse */ template -const Inverse +const Inverse::type, true> MatrixBase::inverse() const { - return Inverse(derived()); + return Inverse::type, true>(derived()); } /** \return the matrix inverse of \c *this, which is assumed to exist. @@ -266,10 +280,10 @@ MatrixBase::inverse() const * \sa class Inverse */ template -const Inverse +const Inverse::type, false> MatrixBase::quickInverse() const { - return Inverse(derived()); + return Inverse::type, false>(derived()); } #endif // EIGEN_INVERSE_H