*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

@@ -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
+           MatrixType::MaxColsAtCompileTime, OtherDerived::MaxColsAtCompileTime> *result
-    > solve(MatrixBase<OtherDerived> *b) const;
+    ) 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.setZero();
-  l.corner(Eigen::TopLeft,HEIGHT,SIZE) = lu.matrixL().corner(Eigen::TopLeft,HEIGHT,SIZE);
+  l.corner(Eigen::TopLeft,m_lu.rows(),smalldim)
-  l.template marked<UnitLower>.solveInPlace(c);
+    = 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());
+  if(!isSurjective())
-  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(c.row(i).isMuchSmallerThan(ei_abs(m_lu.coeff(0,0))))
+    // is c is in the image of U ?
-    {
+    RealScalar biggest_in_c = c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff();
-      d.row(i).setConstant(Scalar(1));
+    for(int col = 0; col < c.cols(); col++)
-    }
+      for(int row = m_rank; row < c.rows(); row++)
-    else return false;
+        if(!ei_isMuchSmallerThan(c.coeff(row,col), biggest_in_c))
-  }
+          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;
   }
+  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.
   *

test/CMakeLists.txt

@@ -38,11 +38,11 @@ MACRO(EI_ADD_TEST testname)
 
   SET(targetname test_${testname})
 
-  IF(${ARGC} EQUAL 2)
+#  IF(${ARGC} EQUAL 2)
-    SET(filename ${ARGV1})
+#    SET(filename ${ARGV1})
-  ELSE(${ARGC} EQUAL 2)
+#  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)