* 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.

Eigen/LU

@@ -1,5 +1,5 @@
-#ifndef EIGEN_LU_H
-#define EIGEN_LU_H
+#ifndef EIGEN_LU_MODULE_H
+#define EIGEN_LU_MODULE_H
 
 #include "Core"
 
@@ -16,9 +16,10 @@ namespace Eigen {
   * \endcode
   */
 
+#include "src/LU/LU.h"
 #include "src/LU/Determinant.h"
 #include "src/LU/Inverse.h"
 
 } // namespace Eigen
 
-#endif // EIGEN_LU_H
+#endif // EIGEN_LU_MODULE_H

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

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;

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());
   }

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;

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, \

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

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

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();
 }
 
 

test/CMakeLists.txt

@@ -2,7 +2,7 @@ IF(BUILD_TESTS)
 
 IF(CMAKE_COMPILER_IS_GNUCXX)
   IF(CMAKE_SYSTEM_NAME MATCHES Linux)
-    SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O1 -g1")
+    SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g2")
     SET(CMAKE_CXX_FLAGS_RELWITHDEBINFO "${CMAKE_CXX_FLAGS_RELWITHDEBINFO} -O2 -g2")
     SET(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -fno-inline-functions")
     SET(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -O0 -g2")
@@ -95,7 +95,7 @@ EI_ADD_TEST(map)
 EI_ADD_TEST(array)
 EI_ADD_TEST(triangular)
 EI_ADD_TEST(cholesky)
-# EI_ADD_TEST(determinant)
+EI_ADD_TEST(determinant)
 EI_ADD_TEST(inverse)
 EI_ADD_TEST(qr)
 EI_ADD_TEST(eigensolver)

test/determinant.cpp

@@ -25,54 +25,50 @@
 #include "main.h"
 #include <Eigen/LU>
 
-template<typename MatrixType> void nullDeterminant(const MatrixType& m)
+template<typename MatrixType> void determinant(const MatrixType& m)
 {
   /* this test covers the following files:
      Determinant.h
   */
-  int rows = m.rows();
-  int cols = m.cols();
+  int size = m.rows();
 
+  MatrixType m1(size, size), m2(size, size);
+  m1.setRandom();
+  m2.setRandom();
   typedef typename MatrixType::Scalar Scalar;
-  typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> SquareMatrixType;
-  typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
-
-  MatrixType dinv(rows, cols), dnotinv(rows, cols);
-
-  dinv.col(0).setOnes();
-  dinv.block(0,1, rows, cols-2).setRandom();
-
-  dnotinv.col(0).setOnes();
-  dnotinv.block(0,1, rows, cols-2).setRandom();
-  dnotinv.col(cols-1).setOnes();
-
-  for (int i=0 ; i<rows ; ++i)
-  {
-    dnotinv.row(i).block(0,1,1,cols-2) = ei_random<Scalar>(99.999999,100.00000001)*dnotinv.row(i).block(0,1,1,cols-2).normalized();
-    dnotinv(i,cols-1) = dnotinv.row(i).block(0,1,1,cols-2).norm2();
-    dinv(i,cols-1) = dinv.row(i).block(0,1,1,cols-2).norm2();
-  }
-
-  SquareMatrixType invertibleCovarianceMatrix = dinv.transpose() * dinv;
-  SquareMatrixType notInvertibleCovarianceMatrix = dnotinv.transpose() * dnotinv;
-
-  std::cout << notInvertibleCovarianceMatrix << "\n" << notInvertibleCovarianceMatrix.determinant() << "\n";
-
-  VERIFY_IS_MUCH_SMALLER_THAN(notInvertibleCovarianceMatrix.determinant(),
-                              notInvertibleCovarianceMatrix.cwise().abs().maxCoeff());
-
-  VERIFY(invertibleCovarianceMatrix.inverse().exists());
-
-  VERIFY(!notInvertibleCovarianceMatrix.inverse().exists());
+  Scalar x = ei_random<Scalar>();
+  VERIFY(ei_isApprox(MatrixType::Identity(size, size).determinant(), Scalar(1)));
+  VERIFY(ei_isApprox((m1*m2).determinant(), m1.determinant() * m2.determinant()));
+  if(size==1) return;
+  int i = ei_random<int>(0, size-1);
+  int j;
+  do {
+    j = ei_random<int>(0, size-1);
+  } while(j==i);
+  m2 = m1;
+  m2.row(i).swap(m2.row(j));
+  VERIFY(ei_isApprox(m2.determinant(), -m1.determinant()));
+  m2 = m1;
+  m2.col(i).swap(m2.col(j));
+  VERIFY(ei_isApprox(m2.determinant(), -m1.determinant()));
+  VERIFY(ei_isApprox(m2.determinant(), m2.transpose().determinant()));
+  VERIFY(ei_isApprox(ei_conj(m2.determinant()), m2.adjoint().determinant()));
+  m2 = m1;
+  m2.row(i) += x*m2.row(j);
+  VERIFY(ei_isApprox(m2.determinant(), m1.determinant()));
+  m2 = m1;
+  m2.row(i) *= x;
+  VERIFY(ei_isApprox(m2.determinant(), m1.determinant() * x));
 }
 
 void test_determinant()
 {
   for(int i = 0; i < g_repeat; i++) {
-    CALL_SUBTEST( nullDeterminant(Matrix<float, 30, 3>()) );
-    CALL_SUBTEST( nullDeterminant(Matrix<double, 30, 3>()) );
-    CALL_SUBTEST( nullDeterminant(Matrix<float, 20, 4>()) );
-    CALL_SUBTEST( nullDeterminant(Matrix<double, 20, 4>()) );
-//     CALL_SUBTEST( nullDeterminant(MatrixXd(20,4));
+    CALL_SUBTEST( determinant(Matrix<float, 1, 1>()) );
+    CALL_SUBTEST( determinant(Matrix<double, 2, 2>()) );
+    CALL_SUBTEST( determinant(Matrix<double, 3, 3>()) );
+    CALL_SUBTEST( determinant(Matrix<double, 4, 4>()) );
+    CALL_SUBTEST( determinant(Matrix<std::complex<double>, 10, 10>()) );
+    CALL_SUBTEST( determinant(MatrixXd(20, 20)) );
   }
 }