* big rework of Inverse.h:

- remove all invertibility checking, will be redundant with LU
  - general case: adapt to matrix storage order for better perf
  - size 4 case: handle corner cases without falling back to gen case.
  - rationalize with selectors instead of compile time if
  - add C-style computeInverse()
* update inverse test.
* in snippets, default cout precision to 3 decimal places
* add some cmake module from kdelibs to support btl with cmake 2.4

Eigen/src/Core/MatrixBase.h

@@ -516,8 +516,8 @@ template<typename Derived> class MatrixBase
 
 /////////// LU module ///////////
 
-    const Inverse<typename ei_eval<Derived>::type, true> inverse() const;
+    const typename ei_eval<Derived>::type inverse() const;
-    const Inverse<typename ei_eval<Derived>::type, false> quickInverse() const;
+    void computeInverse(typename ei_eval<Derived>::type *result) const;
     Scalar determinant() const;
 
 

Eigen/src/LU/Inverse.h

@@ -25,134 +25,113 @@
 #ifndef EIGEN_INVERSE_H
 #define EIGEN_INVERSE_H
 
-/** \lu_module
+/***************************************************************************
-  *
+*** Part 1 : implementation in the general case, by Gaussian elimination ***
-  * \class Inverse
+***************************************************************************/
-  *
+
-  * \brief Inverse of a matrix
+template<typename MatrixType, int StorageOrder>
-  *
+struct ei_compute_inverse_in_general_case;
-  * \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
+template<typename MatrixType>
-  *
+struct ei_compute_inverse_in_general_case<MatrixType, RowMajor>
-  * This class represents the inverse of a matrix. It is the return
+{
-  * type of MatrixBase::inverse() and most of the time this is the only way it
+  static void run(const MatrixType& _matrix, MatrixType *result)
-  * is used.
+  {
-  *
+    typedef typename MatrixType::Scalar Scalar;
-  * \sa MatrixBase::inverse(), MatrixBase::quickInverse()
+    MatrixType matrix(_matrix);
-  */
+    MatrixType &inverse = *result;
-template<typename MatrixType, bool CheckExistence>
+    inverse.setIdentity();
-struct ei_traits<Inverse<MatrixType, CheckExistence> >
+    const int size = matrix.rows();
+    for(int k = 0; k < size-1; k++)
+    {
+      int rowOfBiggest;
+      matrix.col(k).end(size-k).cwise().abs().maxCoeff(&rowOfBiggest);
+      inverse.row(k).swap(inverse.row(k+rowOfBiggest));
+      matrix.row(k).swap(matrix.row(k+rowOfBiggest));
+
+      const Scalar d = matrix(k,k);
+      inverse.block(k+1, 0, size-k-1, size)
+        -= matrix.col(k).end(size-k-1) * (inverse.row(k) / d);
+      matrix.corner(BottomRight, size-k-1, size-k)
+        -= matrix.col(k).end(size-k-1) * (matrix.row(k).end(size-k) / d);
+    }
+
+    for(int k = 0; k < size-1; k++)
+    {
+      const Scalar d = static_cast<Scalar>(1)/matrix(k,k);
+      matrix.row(k).end(size-k) *= d;
+      inverse.row(k) *= d;
+    }
+    inverse.row(size-1) /= matrix(size-1,size-1);
+
+    for(int k = size-1; k >= 1; k--)
+    {
+      inverse.block(0,0,k,size) -= matrix.col(k).start(k) * inverse.row(k);
+    }
+  }
+};
+
+template<typename MatrixType>
+struct ei_compute_inverse_in_general_case<MatrixType, ColMajor>
+{
+  static void run(const MatrixType& _matrix, MatrixType *result)
+  {
+    typedef typename MatrixType::Scalar Scalar;
+    MatrixType matrix(_matrix);
+    MatrixType& inverse = *result;
+    inverse.setIdentity();
+    const int size = matrix.rows();
+    for(int k = 0; k < size-1; k++)
+    {
+      int colOfBiggest;
+      matrix.row(k).end(size-k).cwise().abs().maxCoeff(&colOfBiggest);
+      inverse.col(k).swap(inverse.col(k+colOfBiggest));
+      matrix.col(k).swap(matrix.col(k+colOfBiggest));
+
+      const Scalar d = matrix(k,k);
+      inverse.block(0, k+1, size, size-k-1)
+        -= (inverse.col(k) / d) * matrix.row(k).end(size-k-1);
+      matrix.corner(BottomRight, size-k, size-k-1)
+        -= (matrix.col(k).end(size-k) / d) * matrix.row(k).end(size-k-1);
+    }
+
+    for(int k = 0; k < size-1; k++)
+    {
+      const Scalar d = static_cast<Scalar>(1)/matrix(k,k);
+      matrix.col(k).end(size-k) *= d;
+      inverse.col(k) *= d;
+    }
+    inverse.col(size-1) /= matrix(size-1,size-1);
+
+    for(int k = size-1; k >= 1; k--)
+    {
+      inverse.block(0,0,size,k) -= inverse.col(k) * matrix.row(k).start(k);
+    }
+  }
+};
+
+/********************************************************************
+*** Part 2 : optimized implementations for fixed-size 2,3,4 cases ***
+********************************************************************/
+
+template<typename MatrixType>
+void ei_compute_inverse_in_size2_case(const MatrixType& matrix, MatrixType* result)
 {
   typedef typename MatrixType::Scalar Scalar;
-  enum {
+  const Scalar invdet = Scalar(1) / matrix.determinant();
-    RowsAtCompileTime = MatrixType::RowsAtCompileTime,
+  result->coeffRef(0,0) = matrix.coeff(1,1) * invdet;
-    ColsAtCompileTime = MatrixType::ColsAtCompileTime,
+  result->coeffRef(1,0) = -matrix.coeff(1,0) * invdet;
-    MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
+  result->coeffRef(0,1) = -matrix.coeff(0,1) * invdet;
-    MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
+  result->coeffRef(1,1) = matrix.coeff(0,0) * invdet;
-    Flags = MatrixType::Flags,
-    CoeffReadCost = MatrixType::CoeffReadCost
-  };
-};
-
-template<typename MatrixType, bool CheckExistence> class Inverse : ei_no_assignment_operator,
-  public MatrixBase<Inverse<MatrixType, CheckExistence> >
-{
-  public:
-
-    EIGEN_GENERIC_PUBLIC_INTERFACE(Inverse)
-
-    Inverse(const MatrixType& matrix)
-      : m_inverse(MatrixType::identity(matrix.rows(), matrix.cols()))
-    {
-      if(CheckExistence) m_exists = true;
-      ei_assert(matrix.rows() == matrix.cols());
-      _compute(matrix);
-    }
-
-    /** \returns whether or not the inverse exists.
-      *
-      * \note This method is only available if CheckExistence is set to true, which is the default value.
-      *       For instance, when using quickInverse(), this method is not available.
-      */
-    bool exists() const { assert(CheckExistence); return m_exists; }
-
-    int rows() const { return m_inverse.rows(); }
-    int cols() const { return m_inverse.cols(); }
-
-    const Scalar coeff(int row, int col) const
-    {
-      return m_inverse.coeff(row, col);
-    }
-
-    template<int LoadMode>
-    PacketScalar packet(int row, int col) const
-    {
-      return m_inverse.template packet<LoadMode>(row, col);
-    }
-
-    enum { _Size = MatrixType::RowsAtCompileTime };
-    void _compute(const MatrixType& matrix);
-    void _compute_in_general_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;
-    typename MatrixType::Eval m_inverse;
-};
-
-template<typename MatrixType, bool CheckExistence>
-void Inverse<MatrixType, CheckExistence>
-::_compute_in_general_case(const MatrixType& _matrix)
-{
-  MatrixType matrix(_matrix);
-  const RealScalar max = CheckExistence ? matrix.cwise().abs().maxCoeff()
-                                        : static_cast<RealScalar>(0);
-  const int size = matrix.rows();
-  for(int k = 0; k < size-1; k++)
-  {
-    int rowOfBiggest;
-    const RealScalar max_in_this_col
-      = matrix.col(k).end(size-k).cwise().abs().maxCoeff(&rowOfBiggest);
-    if(CheckExistence && ei_isMuchSmallerThan(max_in_this_col, max))
-    { m_exists = false; return; }
-
-    m_inverse.row(k).swap(m_inverse.row(k+rowOfBiggest));
-    matrix.row(k).swap(matrix.row(k+rowOfBiggest));
-
-    const Scalar d = matrix(k,k);
-    m_inverse.block(k+1, 0, size-k-1, size)
-      -= matrix.col(k).end(size-k-1) * (m_inverse.row(k) / d);
-    matrix.corner(BottomRight, size-k-1, size-k)
-      -= matrix.col(k).end(size-k-1) * (matrix.row(k).end(size-k) / d);
-  }
-
-  for(int k = 0; k < size-1; k++)
-  {
-    const Scalar d = static_cast<Scalar>(1)/matrix(k,k);
-    matrix.row(k).end(size-k) *= d;
-    m_inverse.row(k) *= d;
-  }
-  if(CheckExistence && ei_isMuchSmallerThan(matrix(size-1,size-1), max))
-  { m_exists = false; return; }
-  m_inverse.row(size-1) /= matrix(size-1,size-1);
-
-  for(int k = size-1; k >= 1; k--)
-  {
-    m_inverse.block(0,0,k,size) -= matrix.col(k).start(k) * m_inverse.row(k);
-  }
 }
 
-template<typename ExpressionType, bool CheckExistence>
+template<typename XprType, typename MatrixType>
-bool ei_compute_size2_inverse(const ExpressionType& xpr, typename ExpressionType::Eval* result)
+bool ei_compute_inverse_in_size2_case_with_check(const XprType& matrix, MatrixType* result)
 {
-  typedef typename ExpressionType::Scalar Scalar;
+  typedef typename MatrixType::Scalar Scalar;
-  const typename ei_nested<ExpressionType, 1+CheckExistence>::type matrix(xpr);
   const Scalar det = matrix.determinant();
-  if(CheckExistence && ei_isMuchSmallerThan(det, matrix.cwise().abs().maxCoeff()))
+  if(ei_isMuchSmallerThan(det, matrix.cwise().abs().maxCoeff())) return false;
-    return false;
+  const Scalar invdet = Scalar(1) / det;
-  const Scalar invdet = static_cast<Scalar>(1) / det;
   result->coeffRef(0,0) = matrix.coeff(1,1) * invdet;
   result->coeffRef(1,0) = -matrix.coeff(1,0) * invdet;
   result->coeffRef(0,1) = -matrix.coeff(0,1) * invdet;
@@ -160,34 +139,29 @@ bool ei_compute_size2_inverse(const ExpressionType& xpr, typename ExpressionType
   return true;
 }
 
-template<typename MatrixType, bool CheckExistence>
+template<typename MatrixType>
-void Inverse<MatrixType, CheckExistence>::_compute_in_size3_case(const MatrixType& matrix)
+void ei_compute_inverse_in_size3_case(const MatrixType& matrix, MatrixType* result)
 {
+  typedef typename MatrixType::Scalar Scalar;
   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();
-  const Scalar det = det_minor00 * matrix.coeff(0,0)
+  const Scalar invdet = Scalar(1) / ( det_minor00 * matrix.coeff(0,0)
-                   - det_minor10 * matrix.coeff(1,0)
+                                    - det_minor10 * matrix.coeff(1,0)
-                   + det_minor20 * matrix.coeff(2,0);
+                                    + det_minor20 * matrix.coeff(2,0) );
-  if(CheckExistence && ei_isMuchSmallerThan(det, matrix.cwise().abs().maxCoeff()))
+  result->coeffRef(0, 0) = det_minor00 * invdet;
-    m_exists = false;
+  result->coeffRef(0, 1) = -det_minor10 * invdet;
-  else
+  result->coeffRef(0, 2) = det_minor20 * invdet;
-  {
+  result->coeffRef(1, 0) = -matrix.minor(0,1).determinant() * invdet;
-    const Scalar invdet = static_cast<Scalar>(1) / det;
+  result->coeffRef(1, 1) = matrix.minor(1,1).determinant() * invdet;
-    m_inverse.coeffRef(0, 0) = det_minor00 * invdet;
+  result->coeffRef(1, 2) = -matrix.minor(2,1).determinant() * invdet;
-    m_inverse.coeffRef(0, 1) = -det_minor10 * invdet;
+  result->coeffRef(2, 0) = matrix.minor(0,2).determinant() * invdet;
-    m_inverse.coeffRef(0, 2) = det_minor20 * invdet;
+  result->coeffRef(2, 1) = -matrix.minor(1,2).determinant() * invdet;
-    m_inverse.coeffRef(1, 0) = -matrix.minor(0,1).determinant() * invdet;
+  result->coeffRef(2, 2) = matrix.minor(2,2).determinant() * invdet;
-    m_inverse.coeffRef(1, 1) = matrix.minor(1,1).determinant() * invdet;
-    m_inverse.coeffRef(1, 2) = -matrix.minor(2,1).determinant() * invdet;
-    m_inverse.coeffRef(2, 0) = matrix.minor(0,2).determinant() * invdet;
-    m_inverse.coeffRef(2, 1) = -matrix.minor(1,2).determinant() * invdet;
-    m_inverse.coeffRef(2, 2) = matrix.minor(2,2).determinant() * invdet;
-  }
 }
 
-template<typename MatrixType, bool CheckExistence>
+template<typename MatrixType>
-void Inverse<MatrixType, CheckExistence>::_compute_in_size4_case(const MatrixType& matrix)
+bool ei_compute_inverse_in_size4_case_helper(const MatrixType& matrix, MatrixType* result)
 {
   /* Let's split M into four 2x2 blocks:
     * (P Q)
@@ -205,8 +179,7 @@ void Inverse<MatrixType, CheckExistence>::_compute_in_size4_case(const MatrixTyp
   typedef Block<MatrixType,2,2> XprBlock22;
   typedef typename XprBlock22::Eval Block22;
   Block22 P_inverse;
-
+  if(ei_compute_inverse_in_size2_case_with_check(matrix.template block<2,2>(0,0), &P_inverse))
-  if(ei_compute_size2_inverse<XprBlock22, true>(matrix.template block<2,2>(0,0), &P_inverse))
   {
     const Block22 Q = matrix.template block<2,2>(0,2);
     const Block22 P_inverse_times_Q = P_inverse * Q;
@@ -216,78 +189,147 @@ void Inverse<MatrixType, CheckExistence>::_compute_in_size4_case(const MatrixTyp
     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<Block22, CheckExistence>(X, &Y))
+    ei_compute_inverse_in_size2_case(X, &Y);
-    {
+    result->template block<2,2>(2,2) = Y;
-      m_inverse.template block<2,2>(2,2) = Y;
+    result->template block<2,2>(2,0) = - Y * R_times_P_inverse;
-      m_inverse.template block<2,2>(2,0) = - Y * R_times_P_inverse;
+    const Block22 Z = P_inverse_times_Q * Y;
-      const Block22 Z = P_inverse_times_Q * Y;
+    result->template block<2,2>(0,2) = - Z;
-      m_inverse.template block<2,2>(0,2) = - Z;
+    result->template block<2,2>(0,0) = P_inverse + Z * R_times_P_inverse;
-      m_inverse.template block<2,2>(0,0) = P_inverse + Z * R_times_P_inverse;
+    return true;
-    }
-    else
-    {
-      m_exists = false; return;
-    }
   }
   else
   {
-    _compute_in_general_case(matrix);
+    return false;
   }
 }
 
-template<typename MatrixType, bool CheckExistence>
+template<typename MatrixType>
-void Inverse<MatrixType, CheckExistence>::_compute(const MatrixType& matrix)
+void ei_compute_inverse_in_size4_case(const MatrixType& matrix, MatrixType* result)
 {
-  if(_Size == 1)
+  if(ei_compute_inverse_in_size4_case_helper(matrix, result))
   {
-    const Scalar x = matrix.coeff(0,0);
+    // good ! The topleft 2x2 block was invertible, so the 2x2 blocks approach is successful.
-    if(CheckExistence && x == static_cast<Scalar>(0))
+    return;
-      m_exists = false;
-    else
-      m_inverse.coeffRef(0,0) = static_cast<Scalar>(1) / x;
   }
-  else if(_Size == 2)
+  else
   {
-    if(CheckExistence)
+    // rare case: the topleft 2x2 block is not invertible (but the matrix itself is assumed to be).
-      m_exists = ei_compute_size2_inverse<MatrixType, true>(matrix, &m_inverse);
+    // since this is a rare case, we don't need to optimize it. We just want to handle it with little
+    // additional code.
+    MatrixType m(matrix);
+    m.row(1).swap(m.row(2));
+    if(ei_compute_inverse_in_size4_case_helper(m, result))
+    {
+      // good, the topleft 2x2 block of m is invertible. Since m is different from matrix in that two
+      // rows were permuted, the actual inverse of matrix is derived from the inverse of m by permuting
+      // the corresponding columns.
+      result->col(1).swap(result->col(2));
+    }
     else
-      ei_compute_size2_inverse<MatrixType, false>(matrix, &m_inverse);
+    {
+      // last possible case. Since matrix is assumed to be invertible, this last case has to work.
+      m.row(1).swap(m.row(2));
+      m.row(1).swap(m.row(3));
+      ei_compute_inverse_in_size4_case_helper(m, result);
+      result->col(1).swap(result->col(3));
+    }
   }
-  else if(_Size == 3) _compute_in_size3_case(matrix);
+}
-  else if(_Size == 4) _compute_in_size4_case(matrix);
+
-  else _compute_in_general_case(matrix);
+/***********************************************
+*** Part 3 : selector and MatrixBase methods ***
+***********************************************/
+
+template<typename MatrixType, int Size = MatrixType::RowsAtCompileTime>
+struct ei_compute_inverse
+{
+  static inline void run(const MatrixType& matrix, MatrixType* result)
+  {
+    ei_compute_inverse_in_general_case<MatrixType, MatrixType::Flags&RowMajorBit>
+      ::run(matrix, result);
+  }
+};
+
+template<typename MatrixType>
+struct ei_compute_inverse<MatrixType, 1>
+{
+  static inline void run(const MatrixType& matrix, MatrixType* result)
+  {
+    typedef typename MatrixType::Scalar Scalar;
+    result->coeffRef(0,0) = Scalar(1) / matrix.coeff(0,0);
+  }
+};
+
+template<typename MatrixType>
+struct ei_compute_inverse<MatrixType, 2>
+{
+  static inline void run(const MatrixType& matrix, MatrixType* result)
+  {
+    ei_compute_inverse_in_size2_case(matrix, result);
+  }
+};
+
+template<typename MatrixType>
+struct ei_compute_inverse<MatrixType, 3>
+{
+  static inline void run(const MatrixType& matrix, MatrixType* result)
+  {
+    ei_compute_inverse_in_size3_case(matrix, result);
+  }
+};
+
+template<typename MatrixType>
+struct ei_compute_inverse<MatrixType, 4>
+{
+  static inline void run(const MatrixType& matrix, MatrixType* result)
+  {
+    ei_compute_inverse_in_size4_case(matrix, result);
+  }
+};
+
+/** \lu_module
+  *
+  * Computes the matrix inverse of this matrix.
+  *
+  * \note This matrix must be invertible, otherwise the result is undefined.
+  *
+  * \param result Pointer to the matrix in which to store the result.
+  *
+  * Example: \include MatrixBase_computeInverse.cpp
+  * Output: \verbinclude MatrixBase_computeInverse.out
+  *
+  * \sa inverse()
+  */
+template<typename Derived>
+inline void MatrixBase<Derived>::computeInverse(typename ei_eval<Derived>::type *result) const
+{
+  typedef typename ei_eval<Derived>::type MatrixType;
+  ei_assert(rows() == cols());
+  ei_assert(NumTraits<Scalar>::HasFloatingPoint);
+  ei_compute_inverse<MatrixType>::run(eval(), result);
 }
 
 /** \lu_module
   *
-  * \returns the matrix inverse of \c *this, if it exists.
+  * \returns the matrix inverse of this matrix.
+  *
+  * \note This matrix must be invertible, otherwise the result is undefined.
+  *
+  * \note This method returns a matrix by value, which can be inefficient. To avoid that overhead,
+  * use computeInverse() instead.
   *
   * Example: \include MatrixBase_inverse.cpp
   * Output: \verbinclude MatrixBase_inverse.out
   *
-  * \sa class Inverse
+  * \sa computeInverse()
   */
 template<typename Derived>
-const Inverse<typename ei_eval<Derived>::type, true>
+inline const typename ei_eval<Derived>::type MatrixBase<Derived>::inverse() const
-MatrixBase<Derived>::inverse() const
 {
-  return Inverse<typename Derived::Eval, true>(eval());
+  typedef typename ei_eval<Derived>::type MatrixType;
-}
+  MatrixType result(rows(), cols());
-
+  computeInverse(&result);
-/** \lu_module
+  return result;
-  *
-  * \returns the matrix inverse of \c *this, which is assumed to exist.
-  *
-  * Example: \include MatrixBase_quickInverse.cpp
-  * Output: \verbinclude MatrixBase_quickInverse.out
-  *
-  * \sa class Inverse
-  */
-template<typename Derived>
-const Inverse<typename ei_eval<Derived>::type, false>
-MatrixBase<Derived>::quickInverse() const
-{
-  return Inverse<typename Derived::Eval, false>(eval());
 }
 
 #endif // EIGEN_INVERSE_H

bench/btl/cmake/FindPackageHandleStandardArgs.cmake

@@ -0,0 +1,60 @@
+# FIND_PACKAGE_HANDLE_STANDARD_ARGS(NAME (DEFAULT_MSG|"Custom failure message") VAR1 ... )
+#
+# This macro is intended to be used in FindXXX.cmake modules files.
+# It handles the REQUIRED and QUIET argument to FIND_PACKAGE() and
+# it also sets the <UPPERCASED_NAME>_FOUND variable.
+# The package is found if all variables listed are TRUE.
+# Example:
+#
+#    FIND_PACKAGE_HANDLE_STANDARD_ARGS(LibXml2 DEFAULT_MSG LIBXML2_LIBRARIES LIBXML2_INCLUDE_DIR)
+#
+# LibXml2 is considered to be found, if both LIBXML2_LIBRARIES and 
+# LIBXML2_INCLUDE_DIR are valid. Then also LIBXML2_FOUND is set to TRUE.
+# If it is not found and REQUIRED was used, it fails with FATAL_ERROR, 
+# independent whether QUIET was used or not.
+#
+# If it is found, the location is reported using the VAR1 argument, so 
+# here a message "Found LibXml2: /usr/lib/libxml2.so" will be printed out.
+# If the second argument is DEFAULT_MSG, the message in the failure case will 
+# be "Could NOT find LibXml2", if you don't like this message you can specify
+# your own custom failure message there.
+
+MACRO(FIND_PACKAGE_HANDLE_STANDARD_ARGS _NAME _FAIL_MSG _VAR1 )
+
+  IF("${_FAIL_MSG}" STREQUAL "DEFAULT_MSG")
+    IF (${_NAME}_FIND_REQUIRED)
+      SET(_FAIL_MESSAGE "Could not find REQUIRED package ${_NAME}")
+    ELSE (${_NAME}_FIND_REQUIRED)
+      SET(_FAIL_MESSAGE "Could not find OPTIONAL package ${_NAME}")
+    ENDIF (${_NAME}_FIND_REQUIRED)
+  ELSE("${_FAIL_MSG}" STREQUAL "DEFAULT_MSG")
+    SET(_FAIL_MESSAGE "${_FAIL_MSG}")
+  ENDIF("${_FAIL_MSG}" STREQUAL "DEFAULT_MSG")
+
+  STRING(TOUPPER ${_NAME} _NAME_UPPER)
+
+  SET(${_NAME_UPPER}_FOUND TRUE)
+  IF(NOT ${_VAR1})
+    SET(${_NAME_UPPER}_FOUND FALSE)
+  ENDIF(NOT ${_VAR1})
+
+  FOREACH(_CURRENT_VAR ${ARGN})
+    IF(NOT ${_CURRENT_VAR})
+      SET(${_NAME_UPPER}_FOUND FALSE)
+    ENDIF(NOT ${_CURRENT_VAR})
+  ENDFOREACH(_CURRENT_VAR)
+
+  IF (${_NAME_UPPER}_FOUND)
+    IF (NOT ${_NAME}_FIND_QUIETLY)
+        MESSAGE(STATUS "Found ${_NAME}: ${${_VAR1}}")
+    ENDIF (NOT ${_NAME}_FIND_QUIETLY)
+  ELSE (${_NAME_UPPER}_FOUND)
+    IF (${_NAME}_FIND_REQUIRED)
+        MESSAGE(FATAL_ERROR "${_FAIL_MESSAGE}")
+    ELSE (${_NAME}_FIND_REQUIRED)
+      IF (NOT ${_NAME}_FIND_QUIETLY)
+        MESSAGE(STATUS "${_FAIL_MESSAGE}")
+      ENDIF (NOT ${_NAME}_FIND_QUIETLY)
+    ENDIF (${_NAME}_FIND_REQUIRED)
+  ENDIF (${_NAME_UPPER}_FOUND)
+ENDMACRO(FIND_PACKAGE_HANDLE_STANDARD_ARGS)

doc/snippets/MatrixBase_computeInverse.cpp

@@ -0,0 +1,5 @@
+Matrix3d m = Matrix3d::random();
+cout << "Here is the matrix m:" << endl << m << endl;
+Matrix3d inv;
+m.computeInverse(&inv);
+cout << "Its inverse is:" << endl << inv << endl;

doc/snippets/MatrixBase_inverse.cpp

@@ -1,7 +1,3 @@
-Matrix2d m = Matrix2d::random();
+Matrix3d m = Matrix3d::random();
 cout << "Here is the matrix m:" << endl << m << endl;
-Matrix2d::InverseType m_inv = m.inverse();
+cout << "Its inverse is:" << endl << m.inverse() << endl;
-if(m_inv.exists())
-  cout << "m is invertible, and its inverse is:" << endl << m_inv << endl;
-else
-  cout << "m is not invertible." << endl;

doc/snippets/MatrixBase_quickInverse.cpp

@@ -1,5 +0,0 @@
-Matrix4d m = Matrix4d::zero();
-m.part<Eigen::Upper>().setOnes();
-cout << "Here is the matrix m:" << endl << m << endl;
-cout << "We know for sure that it is invertible." << endl;
-cout << "Here is its inverse:" << m.quickInverse() << endl;

doc/snippets/compile_snippet.cpp.in

@@ -1,10 +1,13 @@
 #include <Eigen/Core>
 #include <Eigen/Array>
 #include <Eigen/LU>
+
 USING_PART_OF_NAMESPACE_EIGEN
 using namespace std;
+
 int main(int, char**)
 {
-${snippet_source_code}
+  cout.precision(3);
-return 0;
+  ${snippet_source_code}
+  return 0;
 }

test/inverse.cpp

@@ -38,32 +38,30 @@ template<typename MatrixType> void inverse(const MatrixType& m)
 
   MatrixType m1 = MatrixType::random(rows, cols),
              m2 = MatrixType::random(rows, cols),
-             m3(rows, cols),
              mzero = MatrixType::zero(rows, cols),
              identity = MatrixType::identity(rows, rows);
 
   m2 = m1.inverse();
   VERIFY_IS_APPROX(m1, m2.inverse() );
 
-  m3 = (m1+m2).inverse();
+  m1.computeInverse(&m2);
-  VERIFY_IS_APPROX(m3+m1, (m1+m2).inverse()+m1);
+  VERIFY_IS_APPROX(m1, m2.inverse() );
 
-  VERIFY_IS_APPROX(m1, m1.inverse().eval().inverse() );
+  VERIFY_IS_APPROX((Scalar(2)*m2).inverse(), m2.inverse()*Scalar(0.5));
 
-  VERIFY_IS_NOT_APPROX(m1, m1.inverse() );
   VERIFY_IS_APPROX(identity, m1.inverse() * m1 );
   VERIFY_IS_APPROX(identity, m1 * m1.inverse() );
 
-  // this one fails:
+  VERIFY_IS_APPROX(m1, m1.inverse().inverse() );
-  VERIFY_IS_APPROX(m1, (m1.inverse()).inverse() );
 }
 
 void test_inverse()
 {
   for(int i = 0; i < 1; i++) {
-    CALL_SUBTEST( inverse(Matrix2f()) );
+    CALL_SUBTEST( inverse(Matrix<double,1,1>()) );
+    CALL_SUBTEST( inverse(Matrix2d()) );
     CALL_SUBTEST( inverse(Matrix3f()) );
-    CALL_SUBTEST( inverse(Matrix4d()) );
+    CALL_SUBTEST( inverse(Matrix4f()) );
-//     CALL_SUBTEST( inverse(MatrixXcd(7,7)) );
+    CALL_SUBTEST( inverse(MatrixXcd(7,7)) );
   }
 }