LU: Let FullPivLU accept a matrix with zero columns

libeigen/eigen!2990

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>

Eigen/src/LU/FullPivLU.h

diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h
index 3c0aad3..13dca4d 100644
--- a/Eigen/src/LU/FullPivLU.h
+++ b/Eigen/src/LU/FullPivLU.h
@@ -362,7 +362,10 @@
   eigen_assert(m_lu.rows() <= NumTraits<PermutationIndex>::highest() &&
                m_lu.cols() <= NumTraits<PermutationIndex>::highest());
 
-  m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff();
+  if (m_lu.cols() > 0)
+    m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff();
+  else
+    m_l1_norm = RealScalar(0);
 
   const Index size = m_lu.diagonalSize();
   const Index rows = m_lu.rows();

test/lu.cpp

diff --git a/test/lu.cpp b/test/lu.cpp
index a4c7512..17675d1 100644
--- a/test/lu.cpp
+++ b/test/lu.cpp
@@ -165,6 +165,33 @@
   VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10);
 }
 
+// Regression test: FullPivLU took the maximum of an empty column-sum vector when computing its l1 norm,
+// so it could not be constructed at all from a matrix with zero columns. PartialPivLU already guarded the
+// same reduction.
+template <typename MatrixType>
+void lu_empty() {
+  typedef typename MatrixType::Scalar Scalar;
+  const Index n = 5;
+
+  FullPivLU<MatrixType> lu{MatrixType(0, 0)};
+  VERIFY_IS_EQUAL(lu.rank(), Index(0));
+  VERIFY(lu.isInvertible());
+  VERIFY_IS_EQUAL(lu.determinant(), Scalar(1));
+
+  lu.compute(MatrixType(0, n));
+  VERIFY_IS_EQUAL(lu.rank(), Index(0));
+  VERIFY_IS_EQUAL(lu.dimensionOfKernel(), n);
+  VERIFY(!lu.isInjective());
+
+  lu.compute(MatrixType(n, 0));
+  VERIFY_IS_EQUAL(lu.rank(), Index(0));
+  VERIFY_IS_EQUAL(lu.dimensionOfKernel(), Index(0));
+  VERIFY(!lu.isSurjective());
+
+  PartialPivLU<MatrixType> plu{MatrixType(0, 0)};
+  VERIFY_IS_EQUAL(plu.determinant(), Scalar(1));
+}
+
 template <typename MatrixType>
 void lu_verify_assert() {
   MatrixType tmp;
@@ -302,6 +329,7 @@
     CALL_SUBTEST_4(lu_invertible<MatrixXd>());
     CALL_SUBTEST_4(lu_partial_piv<MatrixXd>(internal::random<int>(1, EIGEN_TEST_MAX_SIZE)));
     CALL_SUBTEST_4(lu_verify_assert<MatrixXd>());
+    CALL_SUBTEST_4(lu_empty<MatrixXd>());
 
     CALL_SUBTEST_5(lu_non_invertible<MatrixXcf>());
     CALL_SUBTEST_5(lu_invertible<MatrixXcf>());
@@ -311,6 +339,7 @@
     CALL_SUBTEST_6(lu_invertible<MatrixXcd>());
     CALL_SUBTEST_6(lu_partial_piv<MatrixXcd>(internal::random<int>(1, EIGEN_TEST_MAX_SIZE)));
     CALL_SUBTEST_6(lu_verify_assert<MatrixXcd>());
+    CALL_SUBTEST_6(lu_empty<MatrixXcd>());
 
     CALL_SUBTEST_7((lu_non_invertible<Matrix<float, Dynamic, 16> >()));