* reimplement the general case of inverse() on top of LU. Advantages:

- removes much code
  - 2.5x faster (even though LU uses complete pivoting contrary to what
inverse used to do!)
  - there _were_ numeric stability problems with the partial pivoting
approach of inverse(), with 200x200 matrices inversion failed almost
half of the times (overflow). Now these problems are solved thanks to
complete pivoting.

* remove some useless stuff in LU

Eigen/src/LU/Inverse.h

@@ -25,90 +25,8 @@
 #ifndef EIGEN_INVERSE_H
 #define EIGEN_INVERSE_H
 
-/***************************************************************************
-*** Part 1 : implementation in the general case, by Gaussian elimination ***
-***************************************************************************/
-
-template<typename MatrixType, int StorageOrder>
-struct ei_compute_inverse_in_general_case;
-
-template<typename MatrixType>
-struct ei_compute_inverse_in_general_case<MatrixType, RowMajor>
-{
-  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 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);
-      for(int row = k + 1; row < size; row++)
-        inverse.row(row) -= inverse.row(k) * (matrix.coeff(row,k)/d);
-      for(int row = k + 1; row < size; row++)
-        matrix.row(row).end(size-k) -= (matrix.row(k).end(size-k)/d) * matrix.coeff(row,k);
-    }
-
-    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--)
-      for(int row = 0; row < k; row++)
-        inverse.row(row) -= inverse.row(k) * matrix.coeff(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);
-      for(int col = k + 1; col < size; col++)
-        inverse.col(col) -= inverse.col(k) * (matrix.coeff(k,col)/d);
-      for(int col = k + 1; col < size; col++)
-        matrix.col(col).end(size-k) -= (matrix.col(k).end(size-k)/d) * matrix.coeff(k,col);
-    }
-
-    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--)
-      for(int col = 0; col < k; col++)
-        inverse.col(col) -= inverse.col(k) * matrix.coeff(k,col);
-  }
-};
-
 /********************************************************************
-*** Part 2 : optimized implementations for fixed-size 2,3,4 cases ***
+*** Part 1 : optimized implementations for fixed-size 2,3,4 cases ***
 ********************************************************************/
 
 template<typename MatrixType>
@@ -234,7 +152,7 @@ void ei_compute_inverse_in_size4_case(const MatrixType& matrix, MatrixType* resu
 }
 
 /***********************************************
-*** Part 3 : selector and MatrixBase methods ***
+*** Part 2 : selector and MatrixBase methods ***
 ***********************************************/
 
 template<typename MatrixType, int Size = MatrixType::RowsAtCompileTime>
@@ -242,8 +160,8 @@ struct ei_compute_inverse
 {
   static inline void run(const MatrixType& matrix, MatrixType* result)
   {
-    ei_compute_inverse_in_general_case<MatrixType, MatrixType::Flags&RowMajorBit>
+    LU<MatrixType> lu(matrix);
-      ::run(matrix, result);
+    lu.solve(MatrixType::Identity(matrix.rows(), matrix.cols()), result);
   }
 };
 

Eigen/src/LU/LU.h

@@ -58,16 +58,7 @@ template<typename MatrixType> class LU
 
     enum { MaxSmallDimAtCompileTime = EIGEN_ENUM_MIN(
              MatrixType::MaxColsAtCompileTime,
-             MatrixType::MaxRowsAtCompileTime),
+             MatrixType::MaxRowsAtCompileTime)
-           SmallDimAtCompileTime = EIGEN_ENUM_MIN(
-             MatrixType::ColsAtCompileTime,
-             MatrixType::RowsAtCompileTime),
-           MaxBigDimAtCompileTime = EIGEN_ENUM_MAX(
-             MatrixType::MaxColsAtCompileTime,
-             MatrixType::MaxRowsAtCompileTime),
-           BigDimAtCompileTime = EIGEN_ENUM_MAX(
-             MatrixType::ColsAtCompileTime,
-             MatrixType::RowsAtCompileTime)
     };
 
     LU(const MatrixType& matrix);
@@ -107,12 +98,10 @@ template<typename MatrixType> class LU
                  MatrixType::MaxColsAtCompileTime,
                  LU<MatrixType>::MaxSmallDimAtCompileTime> kernel() const;
 
-    template<typename OtherDerived>
+    template<typename OtherDerived, typename ResultType>
     bool solve(
     const MatrixBase<OtherDerived>& b,
-    Matrix<typename MatrixType::Scalar,
+    ResultType *result
-           MatrixType::ColsAtCompileTime, OtherDerived::ColsAtCompileTime,
-           MatrixType::MaxColsAtCompileTime, OtherDerived::MaxColsAtCompileTime> *result
     ) const;
 
     /**
@@ -281,12 +270,10 @@ LU<MatrixType>::kernel() const
 }
 
 template<typename MatrixType>
-template<typename OtherDerived>
+template<typename OtherDerived, typename ResultType>
 bool LU<MatrixType>::solve(
   const MatrixBase<OtherDerived>& b,
-  Matrix<typename MatrixType::Scalar,
+  ResultType *result
-             MatrixType::ColsAtCompileTime, OtherDerived::ColsAtCompileTime,
-             MatrixType::MaxColsAtCompileTime, OtherDerived::MaxColsAtCompileTime> *result
 ) const
 {
   /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}.
@@ -298,7 +285,6 @@ bool LU<MatrixType>::solve(
    */
 
   ei_assert(b.rows() == m_lu.rows());
-  const int bigdim = std::max(m_lu.rows(), m_lu.cols());
   const int smalldim = std::min(m_lu.rows(), m_lu.cols());
 
   typename OtherDerived::Eval c(b.rows(), b.cols());