* LU decomposition, supporting all rectangular matrices, with full

pivoting for better numerical stability. For now the only application is
determinant.
* New determinant unit-test.
* Disable most of Swap.h for now as it makes LU fail (mysterious).
Anyway Swap needs a big overhaul as proposed on IRC.
* Remnants of old class Inverse removed.
* Some warnings fixed.

8 files changed, 197 insertions, 7 deletions

diff --git a/Eigen/src/Core/DummyPacketMath.h b/Eigen/src/Core/DummyPacketMath.h
index 4abd6e997..313001f2f 100644
--- a/Eigen/src/Core/DummyPacketMath.h
+++ b/Eigen/src/Core/DummyPacketMath.h
@@ -139,7 +139,7 @@ template<int Offset,typename PacketType>
 struct ei_palign_impl
 {
   // by default data are aligned, so there is nothing to be done :)
-  inline static void run(PacketType& first, const PacketType& second) {}
+  inline static void run(PacketType&, const PacketType&) {}
 };
 
 /** \internal update \a first using the concatenation of the \a Offset last elements
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index f15e560cb..c2b12a0a3 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -526,6 +526,7 @@ template<typename Derived> class MatrixBase
 
 /////////// LU module ///////////
 
+    const LU<EvalType> lu() const;
     const EvalType inverse() const;
     void computeInverse(EvalType *result) const;
     Scalar determinant() const;
diff --git a/Eigen/src/Core/Swap.h b/Eigen/src/Core/Swap.h
index 35590b56f..3b864789e 100644
--- a/Eigen/src/Core/Swap.h
+++ b/Eigen/src/Core/Swap.h
@@ -42,7 +42,9 @@ template<typename OtherDerived>
 void MatrixBase<Derived>::swap(const MatrixBase<OtherDerived>& other)
 {
   MatrixBase<OtherDerived> *_other = const_cast<MatrixBase<OtherDerived>*>(&other);
-  if(SizeAtCompileTime == Dynamic)
+
+  // disable that path: it makes LU decomposition fail ! I can't see the bug though.
+  if(false /*SizeAtCompileTime == Dynamic*/)
   {
     ei_swap_selector<Derived,OtherDerived>::run(derived(),other.const_cast_derived());
   }
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h
index 76a2f60cf..067ccd0b0 100644
--- a/Eigen/src/Core/util/ForwardDeclarations.h
+++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -95,7 +95,7 @@ void ei_cache_friendly_product(
   bool _rhsRowMajor, const Scalar* _rhs, int _rhsStride,
   bool resRowMajor, Scalar* res, int resStride);
 
-template<typename ExpressionType, bool CheckExistence = true> class Inverse;
+template<typename MatrixType> class LU;
 template<typename MatrixType> class QR;
 template<typename MatrixType> class Cholesky;
 template<typename MatrixType> class CholeskyWithoutSquareRoot;
diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h
index eaa5df27f..40db9a0a7 100644
--- a/Eigen/src/Core/util/Macros.h
+++ b/Eigen/src/Core/util/Macros.h
@@ -135,7 +135,6 @@ typedef typename Eigen::NumTraits<Scalar>::Real RealScalar; \
 typedef typename Base::PacketScalar PacketScalar; \
 typedef typename Eigen::ei_nested<Derived>::type Nested; \
 typedef typename Eigen::ei_eval<Derived>::type Eval; \
-typedef typename Eigen::Inverse<Eval> InverseType; \
 enum { RowsAtCompileTime = Eigen::ei_traits<Derived>::RowsAtCompileTime, \
        ColsAtCompileTime = Eigen::ei_traits<Derived>::ColsAtCompileTime, \
        MaxRowsAtCompileTime = Eigen::ei_traits<Derived>::MaxRowsAtCompileTime, \
diff --git a/Eigen/src/LU/Determinant.h b/Eigen/src/LU/Determinant.h
index fdb1673d8..96e400564 100644
--- a/Eigen/src/LU/Determinant.h
+++ b/Eigen/src/LU/Determinant.h
@@ -63,7 +63,7 @@ const typename Derived::Scalar ei_bruteforce_det(const MatrixBase<Derived>& m)
            - ei_bruteforce_det4_helper(m,1,3,0,2)
            + ei_bruteforce_det4_helper(m,2,3,0,1);
     default:
-      assert(false);
+      ei_internal_assert(false);
   }
 }
 
@@ -85,7 +85,7 @@ typename ei_traits<Derived>::Scalar MatrixBase<Derived>::determinant() const
       return derived().diagonal().redux(ei_scalar_product_op<Scalar>());
   }
   else if(rows() <= 4) return ei_bruteforce_det(derived());
-  else assert(false); // unimplemented for now
+  else return lu().determinant();
 }
 
 #endif // EIGEN_DETERMINANT_H
diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h
new file mode 100644
index 000000000..aa701f2ce
--- /dev/null
+++ b/Eigen/src/LU/LU.h
@@ -0,0 +1,188 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2006-2008 Benoit Jacob <jacob@math.jussieu.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_LU_H
+#define EIGEN_LU_H
+
+/** \ingroup LU_Module
+  *
+  * \class LU
+  *
+  * \brief LU decomposition of a matrix with complete pivoting, and associated features
+  *
+  * \param MatrixType the type of the matrix of which we are computing the LU decomposition
+  *
+  * This class performs a LU decomposition of any matrix, with complete pivoting: the matrix A
+  * is decomposed as A = PLUQ where L is unit-lower-triangular, U is upper-triangular, and P and Q
+  * are permutation matrices.
+  *
+  * This decomposition provides the generic approach to solving systems of linear equations, computing
+  * the rank, invertibility, inverse, and determinant. However for the case when invertibility is
+  * assumed, we have a specialized variant (see MatrixBase::inverse()) achieving better performance.
+  *
+  * \sa MatrixBase::lu(), MatrixBase::determinant(), MatrixBase::rank(), MatrixBase::kernelDim(),
+  *     MatrixBase::kernelBasis(), MatrixBase::solve(), MatrixBase::isInvertible(),
+  *     MatrixBase::inverse(), MatrixBase::computeInverse()
+  */
+template<typename MatrixType> class LU
+{
+  public:
+
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
+    typedef Matrix<int, MatrixType::ColsAtCompileTime, 1> IntRowVectorType;
+    typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
+
+    LU(const MatrixType& matrix);
+
+    const MatrixType& matrixLU() const
+    {
+      return m_lu;
+    }
+
+    const Part<MatrixType, UnitLower> matrixL() const
+    {
+      return m_lu;
+    }
+
+    const Part<MatrixType, Upper> matrixU() const
+    {
+      return m_lu;
+    }
+
+    const IntColVectorType& permutationP() const
+    {
+      return m_p;
+    }
+
+    const IntRowVectorType& permutationQ() const
+    {
+      return m_q;
+    }
+
+    template<typename OtherDerived>
+    typename ProductReturnType<Transpose<MatrixType>, OtherDerived>::Type::Eval
+    solve(const MatrixBase<MatrixType> &b) const;
+
+    /**
+      * This method returns the determinant of the matrix of which
+      * *this is the LU decomposition. It has only linear complexity
+      * (that is, O(n) where n is the dimension of the square matrix)
+      * as the LU decomposition has already been computed.
+      *
+      * Warning: a determinant can be very big or small, so for matrices
+      * of large dimension (like a 50-by-50 matrix) there can be a risk of
+      * overflow/underflow.
+      */
+    typename ei_traits<MatrixType>::Scalar determinant() const;
+
+  protected:
+    MatrixType m_lu;
+    IntColVectorType m_p;
+    IntRowVectorType m_q;
+    int m_det_pq;
+    Scalar m_biggest_eigenvalue_of_u;
+    int m_dimension_of_kernel;
+};
+
+template<typename MatrixType>
+LU<MatrixType>::LU(const MatrixType& matrix)
+  : m_lu(matrix),
+    m_p(matrix.rows()),
+    m_q(matrix.cols())
+{
+  const int size = matrix.diagonal().size();
+  const int rows = matrix.rows();
+  const int cols = matrix.cols();
+
+  IntColVectorType rows_transpositions(matrix.rows());
+  IntRowVectorType cols_transpositions(matrix.cols());
+  int number_of_transpositions = 0;
+
+  for(int k = 0; k < size; k++)
+  {
+    int row_of_biggest, col_of_biggest;
+    const Scalar biggest = m_lu.corner(Eigen::BottomRight, rows-k, cols-k)
+                               .cwise().abs()
+                               .maxCoeff(&row_of_biggest, &col_of_biggest);
+    row_of_biggest += k;
+    col_of_biggest += k;
+    rows_transpositions.coeffRef(k) = row_of_biggest;
+    cols_transpositions.coeffRef(k) = col_of_biggest;
+    if(k != row_of_biggest) {
+      m_lu.row(k).swap(m_lu.row(row_of_biggest));
+      number_of_transpositions++;
+    }
+    if(k != col_of_biggest) {
+      m_lu.col(k).swap(m_lu.col(col_of_biggest));
+      number_of_transpositions++;
+    }
+    const Scalar lu_k_k = m_lu.coeff(k,k);
+    if(ei_isMuchSmallerThan(lu_k_k, biggest)) continue;
+    if(k<rows-1)
+      m_lu.col(k).end(rows-k-1) /= lu_k_k;
+    if(k<size-1)
+      m_lu.corner(BottomRight, rows-k-1, cols-k-1)
+          -= m_lu.col(k).end(rows-k-1) * m_lu.row(k).end(cols-k-1);
+  }
+
+  for(int k = 0; k < matrix.rows(); k++) m_p.coeffRef(k) = k;
+  for(int k = size-1; k >= 0; k--)
+    std::swap(m_p.coeffRef(k), m_p.coeffRef(rows_transpositions.coeff(k)));
+
+  for(int k = 0; k < matrix.cols(); k++) m_q.coeffRef(k) = k;
+  for(int k = 0; k < size; k++)
+    std::swap(m_q.coeffRef(k), m_q.coeffRef(cols_transpositions.coeff(k)));
+
+  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
+
+  int index_of_biggest;
+  m_lu.diagonal().cwise().abs().maxCoeff(&index_of_biggest);
+  m_biggest_eigenvalue_of_u = m_lu.diagonal().coeff(index_of_biggest);
+
+  m_dimension_of_kernel = 0;
+  for(int k = 0; k < size; k++)
+    m_dimension_of_kernel += ei_isMuchSmallerThan(m_lu.diagonal().coeff(k), m_biggest_eigenvalue_of_u);
+}
+
+template<typename MatrixType>
+typename ei_traits<MatrixType>::Scalar LU<MatrixType>::determinant() const
+{
+  Scalar res = m_det_pq;
+  for(int k = 0; k < m_lu.diagonal().size(); k++) res *= m_lu.diagonal().coeff(k);
+  return res;
+}
+
+/** \return the LU decomposition of \c *this.
+  *
+  * \sa class LU
+  */
+template<typename Derived>
+const LU<typename MatrixBase<Derived>::EvalType>
+MatrixBase<Derived>::lu() const
+{
+  return eval();
+}
+
+#endif // EIGEN_LU_H
diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h
index 8ca787c2a..a6cd31f12 100644
--- a/Eigen/src/QR/QR.h
+++ b/Eigen/src/QR/QR.h
@@ -171,7 +171,7 @@ template<typename Derived>
 const QR<typename MatrixBase<Derived>::EvalType>
 MatrixBase<Derived>::qr() const
 {
-  return QR<typename ei_eval<Derived>::type>(derived());
+  return eval();
 }