diff options
author | Benoit Jacob <jacob.benoit.1@gmail.com> | 2009-10-26 14:16:50 -0400 |
---|---|---|
committer | Benoit Jacob <jacob.benoit.1@gmail.com> | 2009-10-26 14:16:50 -0400 |
commit | 44cdbaba4d084b31854ed5bec58f2887f0479b81 (patch) | |
tree | 36d0a427be6ee4dcf9c61eb3a6345b703e043f40 | |
parent | 07d1bcffda9436fe8ce8c56091d68841b8c3be59 (diff) |
* make inverse() do a ReturnByValue
* add computeInverseWithCheck
* doc improvements
* update test
-rw-r--r-- | Eigen/src/Core/MatrixBase.h | 8 | ||||
-rw-r--r-- | Eigen/src/Core/util/ForwardDeclarations.h | 1 | ||||
-rw-r--r-- | Eigen/src/LU/Inverse.h | 76 | ||||
-rw-r--r-- | Eigen/src/LU/PartialLU.h | 16 | ||||
-rw-r--r-- | test/inverse.cpp | 9 |
5 files changed, 89 insertions, 21 deletions
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 24b3feb6f..a2b57f8ea 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -701,7 +701,7 @@ template<typename Derived> class MatrixBase const LU<PlainMatrixType> lu() const; const PartialLU<PlainMatrixType> partialLu() const; - const PlainMatrixType inverse() const; + const ei_inverse_impl<Derived> inverse() const; template<typename ResultType> void computeInverseAndDetWithCheck( ResultType& inverse, @@ -709,6 +709,12 @@ template<typename Derived> class MatrixBase bool& invertible, const RealScalar& absDeterminantThreshold = precision<Scalar>() ) const; + template<typename ResultType> + void computeInverseWithCheck( + ResultType& inverse, + bool& invertible, + const RealScalar& absDeterminantThreshold = precision<Scalar>() + ) const; Scalar determinant() const; /////////// Cholesky module /////////// diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 65e5ce687..01adca96d 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -116,6 +116,7 @@ template<typename MatrixType, int Direction = BothDirections> class Reverse; template<typename MatrixType> class LU; template<typename MatrixType> class PartialLU; +template<typename MatrixType> struct ei_inverse_impl; template<typename MatrixType> class HouseholderQR; template<typename MatrixType> class ColPivotingHouseholderQR; template<typename MatrixType> class FullPivotingHouseholderQR; diff --git a/Eigen/src/LU/Inverse.h b/Eigen/src/LU/Inverse.h index 71dabc663..965866da6 100644 --- a/Eigen/src/LU/Inverse.h +++ b/Eigen/src/LU/Inverse.h @@ -281,6 +281,38 @@ struct ei_compute_inverse_and_det_with_check<MatrixType, ResultType, 4> *** MatrixBase methods *** *************************/ +template<typename MatrixType> +struct ei_traits<ei_inverse_impl<MatrixType> > +{ + typedef typename MatrixType::PlainMatrixType ReturnMatrixType; +}; + +template<typename MatrixType> +struct ei_inverse_impl : public ReturnByValue<ei_inverse_impl<MatrixType> > +{ + // for 2x2, it's worth giving a chance to avoid evaluating. + // for larger sizes, evaluating has negligible cost and limits code size. + typedef typename ei_meta_if< + MatrixType::RowsAtCompileTime == 2, + typename ei_nested<MatrixType,2>::type, + typename ei_eval<MatrixType>::type + >::ret MatrixTypeNested; + typedef typename ei_cleantype<MatrixTypeNested>::type MatrixTypeNestedCleaned; + const MatrixTypeNested m_matrix; + + ei_inverse_impl(const MatrixType& matrix) + : m_matrix(matrix) + {} + + inline int rows() const { return m_matrix.rows(); } + inline int cols() const { return m_matrix.cols(); } + + template<typename Dest> inline void evalTo(Dest& dst) const + { + ei_compute_inverse<MatrixTypeNestedCleaned, Dest>::run(m_matrix, dst); + } +}; + /** \lu_module * * \returns the matrix inverse of this matrix. @@ -299,21 +331,11 @@ struct ei_compute_inverse_and_det_with_check<MatrixType, ResultType, 4> * \sa computeInverseAndDetWithCheck() */ template<typename Derived> -inline const typename MatrixBase<Derived>::PlainMatrixType MatrixBase<Derived>::inverse() const +inline const ei_inverse_impl<Derived> MatrixBase<Derived>::inverse() const { EIGEN_STATIC_ASSERT(NumTraits<Scalar>::HasFloatingPoint,NUMERIC_TYPE_MUST_BE_FLOATING_POINT) ei_assert(rows() == cols()); - typedef typename MatrixBase<Derived>::PlainMatrixType ResultType; - ResultType result(rows(), cols()); - // for 2x2, it's worth giving a chance to avoid evaluating. - // for larger sizes, evaluating has negligible cost and limits code size. - typedef typename ei_meta_if< - RowsAtCompileTime == 2, - typename ei_cleantype<typename ei_nested<Derived,2>::type>::type, - PlainMatrixType - >::ret MatrixType; - ei_compute_inverse<MatrixType, ResultType>::run(derived(), result); - return result; + return ei_inverse_impl<Derived>(derived()); } /** \lu_module @@ -329,7 +351,7 @@ inline const typename MatrixBase<Derived>::PlainMatrixType MatrixBase<Derived>:: * The matrix will be declared invertible if the absolute value of its * determinant is greater than this threshold. * - * \sa inverse() + * \sa inverse(), computeInverseWithCheck() */ template<typename Derived> template<typename ResultType> @@ -353,4 +375,32 @@ inline void MatrixBase<Derived>::computeInverseAndDetWithCheck( (derived(), absDeterminantThreshold, inverse, determinant, invertible); } +/** \lu_module + * + * Computation of matrix inverse, with invertibility check. + * + * This is only for fixed-size square matrices of size up to 4x4. + * + * \param inverse Reference to the matrix in which to store the inverse. + * \param invertible Reference to the bool variable in which to store whether the matrix is invertible. + * \param absDeterminantThreshold Optional parameter controlling the invertibility check. + * The matrix will be declared invertible if the absolute value of its + * determinant is greater than this threshold. + * + * \sa inverse(), computeInverseAndDetWithCheck() + */ +template<typename Derived> +template<typename ResultType> +inline void MatrixBase<Derived>::computeInverseWithCheck( + ResultType& inverse, + bool& invertible, + const RealScalar& absDeterminantThreshold + ) const +{ + RealScalar determinant; + // i'd love to put some static assertions there, but SFINAE means that they have no effect... + ei_assert(rows() == cols()); + computeInverseAndDetWithCheck(inverse,determinant,invertible,absDeterminantThreshold); +} + #endif // EIGEN_INVERSE_H diff --git a/Eigen/src/LU/PartialLU.h b/Eigen/src/LU/PartialLU.h index 30e633eda..e8d21e5ad 100644 --- a/Eigen/src/LU/PartialLU.h +++ b/Eigen/src/LU/PartialLU.h @@ -40,18 +40,20 @@ template<typename MatrixType, typename Rhs> struct ei_partiallu_solve_impl; * is decomposed as A = PLU where L is unit-lower-triangular, U is upper-triangular, and P * is a permutation matrix. * - * Typically, partial pivoting LU decomposition is only considered numerically stable for square invertible matrices. - * So in this class, we plainly require that and take advantage of that to do some simplifications and optimizations. - * This class will assert that the matrix is square, but it won't (actually it can't) check that the matrix is invertible: - * it is your task to check that you only use this decomposition on invertible matrices. + * Typically, partial pivoting LU decomposition is only considered numerically stable for square invertible + * matrices. Thus LAPACK's dgesv and dgesvx require the matrix to be square and invertible. The present class + * does the same. It will assert that the matrix is square, but it won't (actually it can't) check that the + * matrix is invertible: it is your task to check that you only use this decomposition on invertible matrices. * - * The guaranteed safe alternative, working for all matrices, is the full pivoting LU decomposition, provided by class LU. + * The guaranteed safe alternative, working for all matrices, is the full pivoting LU decomposition, provided + * by class LU. * * This is \b not a rank-revealing LU decomposition. Many features are intentionally absent from this class, * such as rank computation. If you need these features, use class LU. * - * This LU decomposition is suitable to invert invertible matrices. It is what MatrixBase::inverse() uses. On the other hand, - * it is \b not suitable to determine whether a given matrix is invertible. + * This LU decomposition is suitable to invert invertible matrices. It is what MatrixBase::inverse() uses + * in the general case. + * On the other hand, it is \b not suitable to determine whether a given matrix is invertible. * * The data of the LU decomposition can be directly accessed through the methods matrixLU(), permutationP(). * diff --git a/test/inverse.cpp b/test/inverse.cpp index b8170a738..269678fd4 100644 --- a/test/inverse.cpp +++ b/test/inverse.cpp @@ -68,17 +68,26 @@ template<typename MatrixType> void inverse(const MatrixType& m) //First: an invertible matrix bool invertible; RealScalar det; + + m2.setZero(); m1.computeInverseAndDetWithCheck(m2, det, invertible); VERIFY(invertible); VERIFY_IS_APPROX(identity, m1*m2); VERIFY_IS_APPROX(det, m1.determinant()); + m2.setZero(); + m1.computeInverseWithCheck(m2, invertible); + VERIFY(invertible); + VERIFY_IS_APPROX(identity, m1*m2); + //Second: a rank one matrix (not invertible, except for 1x1 matrices) VectorType v3 = VectorType::Random(rows); MatrixType m3 = v3*v3.transpose(), m4(rows,cols); m3.computeInverseAndDetWithCheck(m4, det, invertible); VERIFY( rows==1 ? invertible : !invertible ); VERIFY_IS_APPROX(det, m3.determinant()); + m3.computeInverseWithCheck(m4, invertible); + VERIFY( rows==1 ? invertible : !invertible ); #endif } |