Fix SparseLU and SparseQR on row-major and uncompressed sparse inputs

libeigen/eigen!2433

Closes #1964 and #798

Eigen/src/SparseLU/SparseLU.h

diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h
index 9b12f49..124becc 100644
--- a/Eigen/src/SparseLU/SparseLU.h
+++ b/Eigen/src/SparseLU/SparseLU.h
@@ -535,19 +535,21 @@
   OrderingType ord;
   ord(m_mat, m_perm_c);
 
-  // Apply the permutation to the column of the input  matrix
+  // Apply the permutation to the column of the input matrix
   if (m_perm_c.size()) {
-    m_mat.uncompress();  // NOTE: The effect of this command is only to create the InnerNonzeros pointers. FIXME : This
-                         // vector is filled but not subsequently used.
-    // Then, permute only the column pointers
+    // Switch to uncompressed mode so innerNonZeroPtr() exists and can be
+    // permuted consistently with outerIndexPtr().
+    // Downstream sparse traversals may also rely on these per-column counts
+    // while m_mat remains uncompressed.
+    m_mat.uncompress();
+    // A compressed column-major input already exposes valid column pointers.
+    // Otherwise snapshot the internal column-major structure before permuting in place.
+    const bool useInputOuterIndex = !MatrixType::IsRowMajor && mat.isCompressed();
     ei_declare_aligned_stack_constructed_variable(
-        StorageIndex, outerIndexPtr, mat.cols() + 1,
-        mat.isCompressed() ? const_cast<StorageIndex*>(mat.outerIndexPtr()) : 0);
-
-    // If the input matrix 'mat' is uncompressed, then the outer-indices do not match the ones of m_mat, and a copy is
-    // thus needed.
-    if (!mat.isCompressed())
-      IndexVector::Map(outerIndexPtr, mat.cols() + 1) = IndexVector::Map(m_mat.outerIndexPtr(), mat.cols() + 1);
+        StorageIndex, outerIndexPtr, m_mat.cols() + 1,
+        useInputOuterIndex ? const_cast<StorageIndex*>(mat.outerIndexPtr()) : 0);
+    if (!useInputOuterIndex)
+      IndexVector::Map(outerIndexPtr, m_mat.cols() + 1) = IndexVector::Map(m_mat.outerIndexPtr(), m_mat.cols() + 1);
 
     // Apply the permutation and compute the nnz per column.
     for (Index i = 0; i < mat.cols(); i++) {
@@ -618,21 +620,19 @@
   // Apply the column permutation computed in analyzepattern()
   m_mat = matrix;
   if (m_perm_c.size()) {
-    m_mat.uncompress();  // NOTE: The effect of this command is only to create the InnerNonzeros pointers.
-    // Then, permute only the column pointers
-    const StorageIndex* outerIndexPtr;
-    if (matrix.isCompressed())
-      outerIndexPtr = matrix.outerIndexPtr();
-    else {
-      StorageIndex* outerIndexPtr_t = new StorageIndex[matrix.cols() + 1];
-      for (Index i = 0; i <= matrix.cols(); i++) outerIndexPtr_t[i] = m_mat.outerIndexPtr()[i];
-      outerIndexPtr = outerIndexPtr_t;
-    }
+    // Switch to uncompressed mode so innerNonZeroPtr() exists and can be
+    // permuted consistently with outerIndexPtr().
+    m_mat.uncompress();
+    const bool useInputOuterIndex = !MatrixType::IsRowMajor && matrix.isCompressed();
+    ei_declare_aligned_stack_constructed_variable(
+        StorageIndex, outerIndexPtr, m_mat.cols() + 1,
+        useInputOuterIndex ? const_cast<StorageIndex*>(matrix.outerIndexPtr()) : 0);
+    if (!useInputOuterIndex)
+      IndexVector::Map(outerIndexPtr, m_mat.cols() + 1) = IndexVector::Map(m_mat.outerIndexPtr(), m_mat.cols() + 1);
     for (Index i = 0; i < matrix.cols(); i++) {
       m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i];
       m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i + 1] - outerIndexPtr[i];
     }
-    if (!matrix.isCompressed()) delete[] outerIndexPtr;
   } else {  // FIXME This should not be needed if the empty permutation is handled transparently
     m_perm_c.resize(matrix.cols());
     for (StorageIndex i = 0; i < matrix.cols(); ++i) m_perm_c.indices()(i) = i;

Eigen/src/SparseQR/SparseQR.h

diff --git a/Eigen/src/SparseQR/SparseQR.h b/Eigen/src/SparseQR/SparseQR.h
index cd2697c..a866595 100644
--- a/Eigen/src/SparseQR/SparseQR.h
+++ b/Eigen/src/SparseQR/SparseQR.h
@@ -374,16 +374,18 @@
     m_isEtreeOk = true;
   }
 
-  m_pmat.uncompress();  // To have the innerNonZeroPtr allocated
+  // Switch to uncompressed mode so innerNonZeroPtr() exists and can be
+  // permuted consistently with outerIndexPtr().
+  m_pmat.uncompress();
 
   // Apply the fill-in reducing permutation lazily:
   {
-    // If the input is row major, copy the original column indices,
-    // otherwise directly use the input matrix
-    //
+    // A compressed column-major input already exposes valid column pointers.
+    // Otherwise snapshot the internal column-major structure before permuting in place.
     IndexVector originalOuterIndicesCpy;
-    const StorageIndex* originalOuterIndices = mat.outerIndexPtr();
-    if (MatrixType::IsRowMajor) {
+    const bool useInputOuterIndices = !MatrixType::IsRowMajor && mat.isCompressed();
+    const StorageIndex* originalOuterIndices = useInputOuterIndices ? mat.outerIndexPtr() : nullptr;
+    if (!useInputOuterIndices) {
       originalOuterIndicesCpy = IndexVector::Map(m_pmat.outerIndexPtr(), n + 1);
       originalOuterIndices = originalOuterIndicesCpy.data();
     }

test/sparselu.cpp

diff --git a/test/sparselu.cpp b/test/sparselu.cpp
index d933136..fc2de47 100644
--- a/test/sparselu.cpp
+++ b/test/sparselu.cpp
@@ -35,9 +35,74 @@
   check_sparse_square_determinant(sparselu_amd);
 }
 
+template <typename T>
+void test_sparselu_rowmajor_compressed_input() {
+  typedef SparseMatrix<T, RowMajor> RowMajorSparseMatrix;
+  typedef Matrix<T, Dynamic, 1> Vector;
+
+  Vector b(2);
+  b << T(1.1), T(3.14);
+
+  Vector expected(2);
+  expected << T(1.1 - 0.0001 * 3.14), T(3.14);
+
+  RowMajorSparseMatrix compressed(2, 2);
+  compressed.insert(0, 0) = T(1.0);
+  compressed.insert(0, 1) = T(0.0001);
+  compressed.insert(1, 1) = T(1.0);
+  compressed.makeCompressed();
+
+  RowMajorSparseMatrix uncompressed(2, 2);
+  uncompressed.insert(0, 0) = T(1.0);
+  uncompressed.insert(0, 1) = T(0.0001);
+  uncompressed.insert(1, 1) = T(1.0);
+
+  SparseLU<RowMajorSparseMatrix> compressed_solver;
+  compressed_solver.compute(compressed);
+  VERIFY_IS_EQUAL(compressed_solver.info(), Success);
+  VERIFY_IS_APPROX(compressed_solver.solve(b), expected);
+
+  SparseLU<RowMajorSparseMatrix> two_step_solver;
+  two_step_solver.analyzePattern(compressed);
+  two_step_solver.factorize(compressed);
+  VERIFY_IS_EQUAL(two_step_solver.info(), Success);
+  VERIFY_IS_APPROX(two_step_solver.solve(b), expected);
+
+  SparseLU<RowMajorSparseMatrix> uncompressed_solver;
+  uncompressed_solver.compute(uncompressed);
+  VERIFY_IS_EQUAL(uncompressed_solver.info(), Success);
+  VERIFY_IS_APPROX(uncompressed_solver.solve(b), expected);
+}
+
+template <typename T>
+void test_sparselu_colmajor_uncompressed_input() {
+  typedef SparseMatrix<T, ColMajor> ColMajorSparseMatrix;
+  typedef Matrix<T, Dynamic, 1> Vector;
+
+  Vector b(2);
+  b << T(1.1), T(3.14);
+
+  Vector expected(2);
+  expected << T(1.1 - 0.0001 * 3.14), T(3.14);
+
+  ColMajorSparseMatrix uncompressed(2, 2);
+  uncompressed.insert(0, 0) = T(1.0);
+  uncompressed.insert(0, 1) = T(0.0001);
+  uncompressed.insert(1, 1) = T(1.0);
+
+  SparseLU<ColMajorSparseMatrix> uncompressed_solver;
+  uncompressed_solver.compute(uncompressed);
+  VERIFY_IS_EQUAL(uncompressed_solver.info(), Success);
+  VERIFY_IS_APPROX(uncompressed_solver.solve(b), expected);
+}
+
 EIGEN_DECLARE_TEST(sparselu) {
   CALL_SUBTEST_1(test_sparselu_T<float>());
   CALL_SUBTEST_2(test_sparselu_T<double>());
   CALL_SUBTEST_3(test_sparselu_T<std::complex<float> >());
   CALL_SUBTEST_4(test_sparselu_T<std::complex<double> >());
+  CALL_SUBTEST_5(test_sparselu_rowmajor_compressed_input<float>());
+  CALL_SUBTEST_6(test_sparselu_rowmajor_compressed_input<double>());
+  CALL_SUBTEST_7(test_sparselu_colmajor_uncompressed_input<float>());
+  CALL_SUBTEST_8(test_sparselu_colmajor_uncompressed_input<double>());
 }

test/sparseqr.cpp

diff --git a/test/sparseqr.cpp b/test/sparseqr.cpp
index 51ed398..897f426 100644
--- a/test/sparseqr.cpp
+++ b/test/sparseqr.cpp
@@ -132,9 +132,49 @@
   dQ = solver.matrixQ();
   VERIFY_IS_APPROX(Q, dQ);
 }
+
+void test_sparseqr_factorize_uncompressed_input() {
+  typedef SparseMatrix<double, ColMajor> MatrixType;
+  typedef VectorXd Vector;
+  typedef SparseQR<MatrixType, NaturalOrdering<int> > Solver;
+
+  MatrixType uncompressed(2, 2);
+  VectorXi reserve(2);
+  reserve << 2, 3;
+  uncompressed.reserve(reserve);
+  uncompressed.insert(0, 0) = 1.0;
+  uncompressed.insert(0, 1) = 0.5;
+  uncompressed.insert(1, 1) = 1.0;
+
+  // Poison inactive capacity so factorize() must ignore unused slots in the
+  // uncompressed input instead of treating reserved space as structural nnz.
+  uncompressed.innerIndexPtr()[1] = 1;
+  uncompressed.valuePtr()[1] = 7.0;
+
+  MatrixType compressed = uncompressed;
+  compressed.makeCompressed();
+
+  Vector b(2);
+  b << 1.0, 2.0;
+  Vector expected(2);
+  expected << 0.0, 2.0;
+
+  Solver compressed_solver;
+  compressed_solver.compute(compressed);
+  VERIFY_IS_EQUAL(compressed_solver.info(), Success);
+  VERIFY_IS_APPROX(compressed_solver.solve(b), expected);
+
+  Solver two_step_solver;
+  two_step_solver.analyzePattern(compressed);
+  two_step_solver.factorize(uncompressed);
+  VERIFY_IS_EQUAL(two_step_solver.info(), Success);
+  VERIFY_IS_APPROX(two_step_solver.solve(b), expected);
+}
+
 EIGEN_DECLARE_TEST(sparseqr) {
   for (int i = 0; i < g_repeat; ++i) {
     CALL_SUBTEST_1(test_sparseqr_scalar<double>());
     CALL_SUBTEST_2(test_sparseqr_scalar<std::complex<double> >());
   }
+  CALL_SUBTEST_3(test_sparseqr_factorize_uncompressed_input());
 }