aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Eigen/src/Core/MatrixBase.h8
-rw-r--r--Eigen/src/Core/util/ForwardDeclarations.h1
-rw-r--r--Eigen/src/LU/Inverse.h76
-rw-r--r--Eigen/src/LU/PartialLU.h16
-rw-r--r--test/inverse.cpp9
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
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-01 00:39:19 +0000