* add LU unit-test. Seems like we have very good numerical stability!

* some cleanup, and grant me a copyright line on the determinant test.

Eigen/src/LU/Inverse.h

@@ -161,7 +161,7 @@ struct ei_compute_inverse
   static inline void run(const MatrixType& matrix, MatrixType* result)
   {
     LU<MatrixType> lu(matrix);
-    lu.solve(MatrixType::Identity(matrix.rows(), matrix.cols()), result);
+    lu.computeInverse(result);
   }
 };
 

Eigen/src/LU/LU.h

@@ -38,12 +38,9 @@
   * 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
+  * the rank, invertibility, inverse, kernel, and determinant.
-  * assumed, we have a specialized variant (see MatrixBase::inverse()) achieving better performance.
   *
-  * \sa MatrixBase::lu(), MatrixBase::determinant(), MatrixBase::rank(), MatrixBase::kernelDim(),
+  * \sa MatrixBase::lu(), MatrixBase::determinant(), MatrixBase::inverse(), MatrixBase::computeInverse()
-  *     MatrixBase::kernelBasis(), MatrixBase::solve(), MatrixBase::isInvertible(),
-  *     MatrixBase::inverse(), MatrixBase::computeInverse()
   */
 template<typename MatrixType> class LU
 {
@@ -141,6 +138,18 @@ template<typename MatrixType> class LU
       return isInjective() && isSurjective();
     }
 
+    inline void computeInverse(MatrixType *result) const
+    {
+      solve(MatrixType::Identity(m_lu.rows(), m_lu.cols()), result);
+    }
+
+    inline MatrixType inverse() const
+    {
+      MatrixType result;
+      computeInverse(&result);
+      return result;
+    }
+
   protected:
     MatrixType m_lu;
     IntColVectorType m_p;
@@ -163,7 +172,7 @@ LU<MatrixType>::LU(const MatrixType& matrix)
   IntRowVectorType cols_transpositions(matrix.cols());
   int number_of_transpositions = 0;
 
-  RealScalar biggest;
+  RealScalar biggest = RealScalar(0);
   for(int k = 0; k < size; k++)
   {
     int row_of_biggest_in_corner, col_of_biggest_in_corner;
@@ -224,7 +233,7 @@ void LU<MatrixType>::computeKernel(Matrix<typename MatrixType::Scalar,
                                    > *result) const
 {
   ei_assert(!isInvertible());
-  const int dimker = dimensionOfKernel(), rows = m_lu.rows(), cols = m_lu.cols();
+  const int dimker = dimensionOfKernel(), cols = m_lu.cols();
   result->resize(cols, dimker);
 
   /* Let us use the following lemma:
@@ -284,21 +293,22 @@ bool LU<MatrixType>::solve(
    * Step 4: result = Qd;
    */
 
-  ei_assert(b.rows() == m_lu.rows());
+  const int rows = m_lu.rows();
-  const int smalldim = std::min(m_lu.rows(), m_lu.cols());
+  ei_assert(b.rows() == rows);
+  const int smalldim = std::min(rows, m_lu.cols());
 
   typename OtherDerived::Eval c(b.rows(), b.cols());
 
   // Step 1
-  for(int i = 0; i < m_lu.rows(); i++) c.row(m_p.coeff(i)) = b.row(i);
+  for(int i = 0; i < rows; i++) c.row(m_p.coeff(i)) = b.row(i);
 
   // Step 2
   Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime,
          MatrixType::MaxRowsAtCompileTime,
-         MatrixType::MaxRowsAtCompileTime> l(m_lu.rows(), m_lu.rows());
+         MatrixType::MaxRowsAtCompileTime> l(rows, rows);
   l.setZero();
-  l.corner(Eigen::TopLeft,m_lu.rows(),smalldim)
+  l.corner(Eigen::TopLeft,rows,smalldim)
-    = m_lu.corner(Eigen::TopLeft,m_lu.rows(),smalldim);
+    = m_lu.corner(Eigen::TopLeft,rows,smalldim);
   l.template marked<UnitLower>().inverseProductInPlace(c);
 
   // Step 3
@@ -330,7 +340,7 @@ bool LU<MatrixType>::solve(
   * \sa class LU
   */
 template<typename Derived>
-const LU<typename MatrixBase<Derived>::EvalType>
+inline const LU<typename MatrixBase<Derived>::EvalType>
 MatrixBase<Derived>::lu() const
 {
   return eval();

test/CMakeLists.txt

@@ -99,6 +99,7 @@ EI_ADD_TEST(map)
 EI_ADD_TEST(array)
 EI_ADD_TEST(triangular)
 EI_ADD_TEST(cholesky)
+EI_ADD_TEST(lu "-O2")
 EI_ADD_TEST(determinant)
 EI_ADD_TEST(inverse)
 EI_ADD_TEST(qr)

test/determinant.cpp

@@ -1,6 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra. Eigen itself is part of the KDE project.
 //
+// Copyright (C) 2008 Benoit Jacob <jacob@math.jussieu.fr>
 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
 //
 // Eigen is free software; you can redistribute it and/or

test/lu.cpp

@@ -0,0 +1,119 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 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/>.
+
+#include "main.h"
+#include <Eigen/LU>
+
+template<typename Derived>
+void doSomeRankPreservingOperations(Eigen::MatrixBase<Derived>& m)
+{
+  for(int a = 0; a < 3*(m.rows()+m.cols()); a++)
+  {
+    double d = Eigen::ei_random<double>(-1,1);
+    int i = Eigen::ei_random<int>(0,m.rows()-1); // i is a random row number
+    int j;
+    do {
+      j = Eigen::ei_random<int>(0,m.rows()-1);
+    } while (i==j); // j is another one (must be different)
+    m.row(i) += d * m.row(j);
+
+    i = Eigen::ei_random<int>(0,m.cols()-1); // i is a random column number
+    do {
+      j = Eigen::ei_random<int>(0,m.cols()-1);
+    } while (i==j); // j is another one (must be different)
+    m.col(i) += d * m.col(j);
+  }
+}
+
+template<typename MatrixType> void lu_non_invertible()
+{
+  /* this test covers the following files:
+     LU.h
+  */
+  int rows = ei_random<int>(10,200), cols = ei_random<int>(10,200), cols2 = ei_random<int>(10,200);
+  int rank = ei_random<int>(1, std::min(rows, cols)-1);
+
+  MatrixType m1(rows, cols), m2(cols, cols2), m3(rows, cols2), k(1,1);
+  m1.setRandom();
+  if(rows <= cols)
+    for(int i = rank; i < rows; i++) m1.row(i).setZero();
+  else
+    for(int i = rank; i < cols; i++) m1.col(i).setZero();
+  doSomeRankPreservingOperations(m1);
+
+  LU<MatrixType> lu(m1);
+  VERIFY(cols - rank == lu.dimensionOfKernel());
+  VERIFY(rank == lu.rank());
+  VERIFY(!lu.isInjective());
+  VERIFY(!lu.isInvertible());
+  VERIFY(lu.isSurjective() == (lu.rank() == rows));
+  VERIFY((m1 * lu.kernel()).isMuchSmallerThan(m1));
+  lu.computeKernel(&k);
+  VERIFY((m1 * k).isMuchSmallerThan(m1));
+  m2.setRandom();
+  m3 = m1*m2;
+  m2.setRandom();
+  lu.solve(m3, &m2);
+  VERIFY_IS_APPROX(m3, m1*m2);
+  m3.setRandom();
+  VERIFY(!lu.solve(m3, &m2));
+}
+
+template<typename MatrixType> void lu_invertible()
+{
+  /* this test covers the following files:
+     LU.h
+  */
+  int size = ei_random<int>(10,200);
+
+  MatrixType m1(size, size), m2(size, size), m3(size, size);
+  m1.setRandom();
+
+  LU<MatrixType> lu(m1);
+  VERIFY(0 == lu.dimensionOfKernel());
+  VERIFY(size == lu.rank());
+  VERIFY(lu.isInjective());
+  VERIFY(lu.isSurjective());
+  VERIFY(lu.isInvertible());
+  m3.setRandom();
+  lu.solve(m3, &m2);
+  VERIFY_IS_APPROX(m3, m1*m2);
+  VERIFY_IS_APPROX(m2, lu.inverse()*m3);
+  m3.setRandom();
+  VERIFY(lu.solve(m3, &m2));
+}
+
+void test_lu()
+{
+  for(int i = 0; i < g_repeat; i++) {
+    CALL_SUBTEST( lu_non_invertible<MatrixXf>() );
+    CALL_SUBTEST( lu_non_invertible<MatrixXd>() );
+    CALL_SUBTEST( lu_non_invertible<MatrixXcf>() );
+    CALL_SUBTEST( lu_non_invertible<MatrixXcd>() );
+    CALL_SUBTEST( lu_invertible<MatrixXf>() );
+    CALL_SUBTEST( lu_invertible<MatrixXd>() );
+    CALL_SUBTEST( lu_invertible<MatrixXcf>() );
+    CALL_SUBTEST( lu_invertible<MatrixXcd>() );
+  }
+}