From 44cdbaba4d084b31854ed5bec58f2887f0479b81 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Mon, 26 Oct 2009 14:16:50 -0400 Subject: [PATCH] * make inverse() do a ReturnByValue * add computeInverseWithCheck * doc improvements * update test --- Eigen/src/Core/MatrixBase.h | 8 ++- Eigen/src/Core/util/ForwardDeclarations.h | 1 + Eigen/src/LU/Inverse.h | 76 +++++++++++++++++++---- Eigen/src/LU/PartialLU.h | 16 ++--- 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 class MatrixBase const LU lu() const; const PartialLU partialLu() const; - const PlainMatrixType inverse() const; + const ei_inverse_impl inverse() const; template void computeInverseAndDetWithCheck( ResultType& inverse, @@ -709,6 +709,12 @@ template class MatrixBase bool& invertible, const RealScalar& absDeterminantThreshold = precision() ) const; + template + void computeInverseWithCheck( + ResultType& inverse, + bool& invertible, + const RealScalar& absDeterminantThreshold = precision() + ) 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 class Reverse; template class LU; template class PartialLU; +template struct ei_inverse_impl; template class HouseholderQR; template class ColPivotingHouseholderQR; template 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 *** MatrixBase methods *** *************************/ +template +struct ei_traits > +{ + typedef typename MatrixType::PlainMatrixType ReturnMatrixType; +}; + +template +struct ei_inverse_impl : public ReturnByValue > +{ + // 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::type, + typename ei_eval::type + >::ret MatrixTypeNested; + typedef typename ei_cleantype::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 inline void evalTo(Dest& dst) const + { + ei_compute_inverse::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 * \sa computeInverseAndDetWithCheck() */ template -inline const typename MatrixBase::PlainMatrixType MatrixBase::inverse() const +inline const ei_inverse_impl MatrixBase::inverse() const { EIGEN_STATIC_ASSERT(NumTraits::HasFloatingPoint,NUMERIC_TYPE_MUST_BE_FLOATING_POINT) ei_assert(rows() == cols()); - typedef typename MatrixBase::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::type>::type, - PlainMatrixType - >::ret MatrixType; - ei_compute_inverse::run(derived(), result); - return result; + return ei_inverse_impl(derived()); } /** \lu_module @@ -329,7 +351,7 @@ inline const typename MatrixBase::PlainMatrixType MatrixBase:: * The matrix will be declared invertible if the absolute value of its * determinant is greater than this threshold. * - * \sa inverse() + * \sa inverse(), computeInverseWithCheck() */ template template @@ -353,4 +375,32 @@ inline void MatrixBase::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 +template +inline void MatrixBase::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 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 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 }