From 681e9446705449a17d7c5ba669b91c068f9b98a0 Mon Sep 17 00:00:00 2001
From: Benoit Jacob <jacob.benoit.1@gmail.com>
Date: Sat, 9 Aug 2008 15:23:54 +0000
Subject: [PATCH] *implement LU solver (solves any rectangular system) *in
 test/CMakeLists : modify EI_ADD_TEST so that 2nd argument is additional
 compiler flags. used to add -O2 to test_product_large so it doesn't take
 forever.

---
 Eigen/src/LU/LU.h   | 71 ++++++++++++++++++++-------------------------
 test/CMakeLists.txt | 14 +++++----
 2 files changed, 41 insertions(+), 44 deletions(-)

diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h
index f8c9eedb9..6a55921fe 100644
--- a/Eigen/src/LU/LU.h
+++ b/Eigen/src/LU/LU.h
@@ -108,10 +108,12 @@ template<typename MatrixType> class LU
                  LU<MatrixType>::MaxSmallDimAtCompileTime> kernel() const;
 
     template<typename OtherDerived>
+    bool solve(
+    const MatrixBase<OtherDerived>& b,
     Matrix<typename MatrixType::Scalar,
            MatrixType::ColsAtCompileTime, OtherDerived::ColsAtCompileTime,
-           MatrixType::MaxColsAtCompileTime, OtherDerived::MaxColsAtCompileTime
-    > solve(MatrixBase<OtherDerived> *b) const;
+           MatrixType::MaxColsAtCompileTime, OtherDerived::MaxColsAtCompileTime> *result
+    ) const;
 
     /**
       * This method returns the determinant of the matrix of which
@@ -278,7 +280,6 @@ LU<MatrixType>::kernel() const
   return result;
 }
 
-#if 0
 template<typename MatrixType>
 template<typename OtherDerived>
 bool LU<MatrixType>::solve(
@@ -296,55 +297,47 @@ bool LU<MatrixType>::solve(
    * Step 4: result = Qd;
    */
 
+  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());
-  Matrix<typename MatrixType::Scalar,
-         MatrixType::ColsAtCompileTime, OtherDerived::ColsAtCompileTime,
-         MatrixType::MaxColsAtCompileTime, OtherDerived::MaxColsAtCompileTime>
-    d(m_lu.cols(), b.cols());
 
   // Step 1
-  for(int i = 0; i < dim(); i++) c.row(m_p.coeff(i)) = b.row(i);
+  for(int i = 0; i < m_lu.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());
-  l.setIdentity();
-  l.corner(Eigen::TopLeft,HEIGHT,SIZE) = lu.matrixL().corner(Eigen::TopLeft,HEIGHT,SIZE);
-  l.template marked<UnitLower>.solveInPlace(c);
+  l.setZero();
+  l.corner(Eigen::TopLeft,m_lu.rows(),smalldim)
+    = m_lu.corner(Eigen::TopLeft,m_lu.rows(),smalldim);
+  l.template marked<UnitLower>().inverseProductInPlace(c);
 
   // Step 3
-  const int bigdim = std::max(m_lu.rows(), m_lu.cols());
-  const int smalldim = std::min(m_lu.rows(), m_lu.cols());
-  Matrix<Scalar, MatrixType::BigDimAtCompileTime, MatrixType::BigDimAtCompileTime,
-         MatrixType::MaxBigDimAtCompileTime,
-         MatrixType::MaxBigDimAtCompileTime> u(bigdim, bigdim);
-  u.setZero();
-  u.corner(TopLeft, smalldim, smalldim) = m_lu.corner(TopLeft, smalldim, smalldim)
-                                              .template part<Upper>();
-  if(m_lu.cols() > m_lu.rows())
-    u.corner(BottomLeft, m_lu.cols()-m_lu.rows(), m_lu.cols()).setZero();
-  const int size = std::min(m_lu.rows(), m_lu.cols());
-  for(int i = size-1; i >= m_rank; i--)
+  if(!isSurjective())
   {
-    if(c.row(i).isMuchSmallerThan(ei_abs(m_lu.coeff(0,0))))
-    {
-      d.row(i).setConstant(Scalar(1));
-    }
-    else return false;
-  }
-  for(int i = m_rank-1; i >= 0; i--)
-  {
-    d.row(i) = c.row(i);
-    for( int j = i + 1; j <= dim() - 1; j++ )
-    {
-        rowptr += dim();
-        b[i] -= b[j] * (*rowptr);
-    }
-    b[i] /= *denomptr;
+    // is c is in the image of U ?
+    RealScalar biggest_in_c = c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff();
+    for(int col = 0; col < c.cols(); col++)
+      for(int row = m_rank; row < c.rows(); row++)
+        if(!ei_isMuchSmallerThan(c.coeff(row,col), biggest_in_c))
+          return false;
   }
+  Matrix<Scalar, Dynamic, OtherDerived::ColsAtCompileTime,
+         MatrixType::MaxRowsAtCompileTime, OtherDerived::MaxColsAtCompileTime>
+    d(c.corner(TopLeft, m_rank, c.cols()));
+  m_lu.corner(TopLeft, m_rank, m_rank)
+      .template marked<Upper>()
+      .inverseProductInPlace(d);
+
+  // Step 4
+  result->resize(m_lu.cols(), b.cols());
+  for(int i = 0; i < m_rank; i++) result->row(m_q.coeff(i)) = d.row(i);
+  for(int i = m_rank; i < m_lu.cols(); i++) result->row(m_q.coeff(i)).setZero();
+  return true;
 }
-#endif
 
 /** \return the LU decomposition of \c *this.
   *
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index f3f3a625c..aa3860575 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -38,11 +38,11 @@ MACRO(EI_ADD_TEST testname)
 
   SET(targetname test_${testname})
 
-  IF(${ARGC} EQUAL 2)
-    SET(filename ${ARGV1})
-  ELSE(${ARGC} EQUAL 2)
+#  IF(${ARGC} EQUAL 2)
+#    SET(filename ${ARGV1})
+#  ELSE(${ARGC} EQUAL 2)
     SET(filename ${testname}.cpp)
-  ENDIF(${ARGC} EQUAL 2)
+#  ENDIF(${ARGC} EQUAL 2)
   ADD_EXECUTABLE(${targetname} ${filename})
 
   IF(NOT EIGEN_NO_ASSERTION_CHECKING)
@@ -59,6 +59,10 @@ MACRO(EI_ADD_TEST testname)
 
   ENDIF(NOT EIGEN_NO_ASSERTION_CHECKING)
 
+  IF(${ARGC} GREATER 1)
+    EI_ADD_TARGET_PROPERTY(${targetname} COMPILE_FLAGS "${ARGV1}")
+  ENDIF(${ARGC} GREATER 1)
+
   EI_ADD_TARGET_PROPERTY(${targetname} COMPILE_FLAGS "-DEIGEN_TEST_FUNC=${testname}")
 
   IF(TEST_LIB)
@@ -86,7 +90,7 @@ EI_ADD_TEST(basicstuff)
 EI_ADD_TEST(linearstructure)
 EI_ADD_TEST(cwiseop)
 EI_ADD_TEST(product_small)
-EI_ADD_TEST(product_large)
+EI_ADD_TEST(product_large "-O2")
 EI_ADD_TEST(adjoint)
 EI_ADD_TEST(submatrices)
 EI_ADD_TEST(miscmatrices)