* make inverse() do a ReturnByValue

* add computeInverseWithCheck
* doc improvements
* update test

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 ///////////

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;

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;
+  return ei_inverse_impl<Derived>(derived());
-  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;
 }
 
 /** \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

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.
+  * Typically, partial pivoting LU decomposition is only considered numerically stable for square invertible
-  * So in this class, we plainly require that and take advantage of that to do some simplifications and optimizations.
+  * matrices. Thus LAPACK's dgesv and dgesvx require the matrix to be square and invertible. The present class
-  * This class will assert that the matrix is square, but it won't (actually it can't) check that the matrix is invertible:
+  * does the same. It will assert that the matrix is square, but it won't (actually it can't) check that the
-  * it is your task to check that you only use this decomposition on invertible matrices.
+  * 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,
+  * This LU decomposition is suitable to invert invertible matrices. It is what MatrixBase::inverse() uses
-  * it is \b not suitable to determine whether a given matrix is invertible.
+  * 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().
   *

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
 }