Reach-based sparse triangular solveInPlace (Gilbert–Peierls) for col-major L / Lᵀ

libeigen/eigen!2678

Eigen/SparseCore

@@ -58,6 +58,7 @@
 #include "src/SparseCore/SparseDenseProduct.h"
 #include "src/SparseCore/SparseSelfAdjointView.h"
 #include "src/SparseCore/SparseTriangularView.h"
+#include "src/SparseCore/TriangularReachSolver.h"
 #include "src/SparseCore/TriangularSolver.h"
 #include "src/SparseCore/SparsePermutation.h"
 #include "src/SparseCore/SparseFuzzy.h"

Eigen/src/SparseCore/SparseUtil.h

@@ -61,6 +61,17 @@ class SparseSymmetricPermutationProduct;
 
 namespace internal {
 
+/** \internal
+ * Compile-time test for whether a sparse expression exposes its storage directly
+ * through outerIndexPtr() / innerIndexPtr() / valuePtr() / innerNonZeroPtr()
+ * (the \ref CompressedAccessBit contract: SparseMatrix, SparseVector,
+ * Ref<Sparse>, Map<SparseMatrix>, ...). Such expressions can drive raw-pointer
+ * fast paths; anything else must be walked via InnerIterator.
+ */
+template <typename Derived>
+struct has_compressed_access : std::integral_constant<bool, (int(traits<Derived>::Flags) & CompressedAccessBit) != 0> {
+};
+
 template <typename T, int Rows, int Cols, int Flags>
 struct sparse_eval;
 

Eigen/src/SparseCore/TriangularReachSolver.h

@@ -0,0 +1,323 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+// SPDX-FileCopyrightText: The Eigen Authors
+// SPDX-License-Identifier: MPL-2.0
+
+#ifndef EIGEN_TRIANGULAR_REACH_SOLVER_H
+#define EIGEN_TRIANGULAR_REACH_SOLVER_H
+
+// IWYU pragma: private
+#include "./InternalHeaderCheck.h"
+
+namespace Eigen {
+
+namespace internal {
+
+// Gilbert-Peierls sparse-rhs triangular solve.
+//
+// Solves T x = b for a column-major sparse triangular T (lower OR upper) with sorted
+// columns, and a sparse right-hand side b. The cost is O(|reach| + flops): only the
+// columns reachable from b's pattern are touched, independent of the dimension.
+//
+// T is the triangular VIEW of a possibly-general stored matrix, so a column may hold
+// entries on BOTH sides of the diagonal; only the entries in the active triangle
+// participate. Everything is therefore direction-aware via a compile-time bool `Upper`:
+//   - the reach follows only in-triangle successors (strictly below the diagonal for
+//     lower, strictly above for upper -- one comparison per stored entry);
+//   - the numeric sweep locates the diagonal within the sorted column (first entry with
+//     index >= j for lower, last with index <= j for upper) and updates only the
+//     in-triangle off-diagonal entries.
+// For a genuinely triangular stored matrix (Eigen's SparseMatrix) the in-triangle
+// entries are the whole column and the diagonal sits at the extreme, so the boundary
+// scans are single-comparison no-ops.
+
+// ===========================================================================
+// Reach: any triangular view, raw CSC storage.
+// ===========================================================================
+
+// Computes reach_{G(T)}(pattern(b)) via a non-recursive depth-first search, emitting
+// the reached columns into xi[top..n) in topological (solve) order and returning top.
+// `mark` (length-n byte array; a 0/1 visited flag needs no more) must be all-zero on
+// entry; every reached node is flagged, and since the reached set is exactly the
+// output, the caller clears those flags again while gathering (no reset needed). `xi`
+// and `pstack` are size-n scratch; the DFS stack occupies xi[0..head] while the output
+// grows down from xi[n), and head < top holds so they never overlap. Only in-triangle
+// stored entries (index > j for lower, index < j for upper) are reach successors; the
+// diagonal and any out-of-triangle entries of the stored matrix are skipped.
+// `innerNonZeroPtr` is the per-column nonzero count: pass it for an uncompressed matrix
+// so column j ends at outerIndexPtr[j]+innerNonZeroPtr[j]; pass nullptr (compressed) to
+// end at outerIndexPtr[j+1].
+template <bool Upper, typename StorageIndex>
+Index triangular_reach(const StorageIndex* outerIndexPtr, const StorageIndex* innerIndexPtr,
+                       const StorageIndex* innerNonZeroPtr, const StorageIndex* bIdx, Index bCount, StorageIndex* xi,
+                       StorageIndex* pstack, uint8_t* mark, Index n) {
+  Index top = n;
+  for (Index r = 0; r < bCount; ++r) {
+    StorageIndex root = bIdx[r];
+    if (mark[root]) continue;
+
+    Index head = 0;
+    xi[0] = root;
+    while (head >= 0) {
+      StorageIndex j = xi[head];
+      Index colBeg = outerIndexPtr[j];
+      Index colEnd = innerNonZeroPtr ? outerIndexPtr[j] + innerNonZeroPtr[j] : outerIndexPtr[j + 1];
+      if (!mark[j]) {
+        mark[j] = 1;
+        pstack[head] = StorageIndex(colBeg);
+      }
+      bool done = true;
+      for (Index p = pstack[head]; p < colEnd; ++p) {
+        StorageIndex i = innerIndexPtr[p];
+        if (Upper ? (i >= j) : (i <= j)) continue;  // out of triangle, or the diagonal
+        if (mark[i]) continue;                      // already visited
+        pstack[head] = StorageIndex(p + 1);
+        xi[++head] = i;  // descend
+        done = false;
+        break;
+      }
+      if (done) {  // no unvisited successor: postorder j
+        xi[--top] = j;
+        --head;
+      }
+    }
+  }
+  return top;
+}
+
+// ===========================================================================
+// Numeric sweep over a precomputed reach.
+// ===========================================================================
+
+// Solves T x = b in place on the dense accumulator x (zero except where b was
+// scattered), touching only the reached columns in the order produced by the reach.
+// Columns are sorted, so the active triangle is a contiguous run: for lower it is the
+// suffix from the first entry with index >= j (the diagonal, then the sub-diagonal
+// off-diagonals); for upper it is the prefix up to the last entry with index <= j (the
+// off-diagonals, then the diagonal). Locating the boundary is O(1) for a genuinely
+// triangular column (the extreme stored entry is already the diagonal); only when a
+// stored column carries out-of-triangle entries (a general matrix seen through a
+// TriangularView) do we binary-search past them, in O(log nnz/col) rather than a
+// linear scan of the wrong-side run.
+// The stored lhs scalar (LhsScalar, read from valuePtr) and the accumulator/rhs scalar
+// (RhsScalar, held in x) are separate: a real factor applied to a complex rhs must
+// accumulate in the complex type. The arithmetic runs in RhsScalar; LhsScalar values
+// promote to it (e.g. double -> complex<double>).
+template <bool Upper, bool UnitDiag, typename StorageIndex, typename LhsScalar, typename RhsScalar>
+void triangular_solve_over_reach(const StorageIndex* outerIndexPtr, const StorageIndex* innerIndexPtr,
+                                 const LhsScalar* valuePtr, const StorageIndex* innerNonZeroPtr, const StorageIndex* xi,
+                                 Index top, Index n, RhsScalar* x) {
+  for (Index k = top; k < n; ++k) {
+    StorageIndex j = xi[k];
+    Index colBeg = outerIndexPtr[j];
+    Index colEnd = innerNonZeroPtr ? outerIndexPtr[j] + innerNonZeroPtr[j] : outerIndexPtr[j + 1];
+    RhsScalar xj;
+    Index offBeg, offEnd;
+    EIGEN_IF_CONSTEXPR (Upper) {
+      // e = one past the last in-triangle entry (index <= j)
+      Index e = colEnd;
+      if (e > colBeg && innerIndexPtr[e - 1] > j)  // wrong-side (below-diagonal) tail: skip it
+        e = std::upper_bound(innerIndexPtr + colBeg, innerIndexPtr + colEnd, j) - innerIndexPtr;
+      bool hasDiag = e > colBeg && innerIndexPtr[e - 1] == j;
+      offBeg = colBeg;
+      offEnd = hasDiag ? e - 1 : e;  // drop the diagonal slot from the update iff it is stored
+      EIGEN_IF_CONSTEXPR (!UnitDiag) {
+        eigen_assert(hasDiag && "sparse triangular solve: missing diagonal");
+        // Missing diagonal is out of contract; match the old AmbiVector path -- divide by 0
+        // for a deterministic inf/NaN, rather than reading valuePtr[e-1] out of bounds.
+        x[j] /= hasDiag ? valuePtr[e - 1] : LhsScalar(0);
+      }
+    } else {
+      // s = first in-triangle entry (index >= j)
+      Index s = colBeg;
+      if (s < colEnd && innerIndexPtr[s] < j)  // wrong-side (above-diagonal) head: skip it
+        s = std::lower_bound(innerIndexPtr + colBeg, innerIndexPtr + colEnd, j) - innerIndexPtr;
+      bool hasDiag = s < colEnd && innerIndexPtr[s] == j;
+      offBeg = hasDiag ? s + 1 : s;  // drop the diagonal slot from the update iff it is stored
+      offEnd = colEnd;
+      EIGEN_IF_CONSTEXPR (!UnitDiag) {
+        eigen_assert(hasDiag && "sparse triangular solve: missing diagonal");
+        x[j] /= hasDiag ? valuePtr[s] : LhsScalar(0);  // missing diagonal -> inf/NaN, not an OOB read
+      }
+    }
+    xj = x[j];
+    for (Index p = offBeg; p < offEnd; ++p) {
+      StorageIndex i = innerIndexPtr[p];
+      x[i] = numext::madd<RhsScalar>(-xj, RhsScalar(valuePtr[p]), x[i]);
+    }
+  }
+}
+
+// ===========================================================================
+// Borrow-a-buffer solve core.
+//
+// The caller owns the scratch and passes it in, so a reused buffer makes repeated
+// solves allocation-free. All buffers are restored on exit, so one setup suffices for
+// many solves.
+// ===========================================================================
+
+// Core: compute reach(pattern(b)) and run the numeric sweep, returning top -- WITHOUT
+// any cleanup. The rhs must ALREADY be scattered into xwork (xwork[bIdx[r]] = value) by
+// the caller; bIdx is the rhs pattern (the reach roots). Pulling the scatter out lets a
+// caller reading the rhs through an iterator scatter as it reads, dropping the separate
+// value array. On return, xi = iwork[top..n) holds the reached columns in topological
+// order, xwork holds their solution values, and mark is set on the reached set; the
+// caller consumes xwork/xi and clears them.
+// Solving T x = b for a column-major, sorted triangular T (lower or upper):
+//   - iwork: >= 2n StorageIndex, carved into xi | pstack (each length n).
+//   - mark:  >= n bytes, all-zero.
+//   - xwork: >= n Scalar, the dense accumulator, zero except b scattered on bIdx.
+// `innerNonZeroPtr` is nullptr for a compressed T, or the per-column nonzero count for
+// an uncompressed T (columns then end at outerIndexPtr[j]+innerNonZeroPtr[j]).
+template <bool Upper, bool UnitDiag, typename StorageIndex, typename LhsScalar, typename RhsScalar>
+Index reach_solve_dense(const StorageIndex* outerIndexPtr, const StorageIndex* innerIndexPtr, const LhsScalar* valuePtr,
+                        const StorageIndex* innerNonZeroPtr, Index n, const StorageIndex* bIdx, Index bCount,
+                        StorageIndex* iwork, uint8_t* mark, RhsScalar* xwork) {
+  Index top =
+      triangular_reach<Upper>(outerIndexPtr, innerIndexPtr, innerNonZeroPtr, bIdx, bCount, iwork, iwork + n, mark, n);
+  triangular_solve_over_reach<Upper, UnitDiag>(outerIndexPtr, innerIndexPtr, valuePtr, innerNonZeroPtr, iwork, top, n,
+                                               xwork);
+  return top;
+}
+
+// ---------------------------------------------------------------------------
+// Iterator-driven path: for a triangular, column-major sparse expression that does NOT
+// expose raw CSC storage (has_compressed_access is false). Columns are read through the
+// expression's evaluator InnerIterator. Because an InnerIterator can't cheaply hold DFS
+// resume state, the reach uses a mark-on-push worklist plus a final sort into
+// topological order (ascending index for lower, descending for upper -- both valid
+// topological orders for the respective solve); the log factor is empirically ~free.
+// ---------------------------------------------------------------------------
+
+// Reach via a mark-on-push worklist that opens one column at a time through
+// InnerIterator. Only in-triangle successors (index > j for lower, index < j for upper)
+// are pushed. Scratch: xi (n, shared stack+output) and mark (n bytes, all-zero in,
+// reach flagged out). Returns top; xi[top..n) is sorted into the topological order for
+// the solve direction.
+template <bool Upper, typename Eval, typename StorageIndex>
+Index triangular_reach_iter(const Eval& mat, const StorageIndex* bIdx, Index bCount, StorageIndex* xi, uint8_t* mark,
+                            Index n) {
+  Index top = n;
+  Index sp = 0;
+  for (Index r = 0; r < bCount; ++r) {
+    StorageIndex root = bIdx[r];
+    if (!mark[root]) {
+      mark[root] = 1;
+      xi[sp++] = root;
+    }
+  }
+  while (sp > 0) {
+    StorageIndex j = xi[--sp];
+    xi[--top] = j;  // collect
+    for (typename Eval::InnerIterator it(mat, j); it; ++it) {
+      StorageIndex i = StorageIndex(it.index());
+      if (Upper ? (i >= j) : (i <= j)) continue;  // out of triangle, or the diagonal
+      if (!mark[i]) {
+        mark[i] = 1;
+        xi[sp++] = i;
+      }
+    }
+  }
+  // descending for upper, ascending for lower (comparator type picked at compile time)
+  using Comp = std::conditional_t<Upper, std::greater<StorageIndex>, std::less<StorageIndex>>;
+  std::sort(xi + top, xi + n, Comp{});
+  return top;
+}
+
+// Numeric sweep over the reach, reading columns through InnerIterator. Direction-aware
+// like the pointer sweep: for lower the leading out-of-triangle entries (index < j) are
+// skipped so the diagonal is the first remaining entry; for upper the diagonal is found
+// by scan and only the entries with index < j are updated.
+template <bool Upper, bool UnitDiag, typename Eval, typename StorageIndex, typename Scalar>
+void triangular_solve_over_reach_iter(const Eval& mat, const StorageIndex* xi, Index top, Index n, Scalar* x) {
+  for (Index k = top; k < n; ++k) {
+    StorageIndex j = xi[k];
+    EIGEN_IF_CONSTEXPR (Upper) {
+      EIGEN_IF_CONSTEXPR (!UnitDiag) {
+        Scalar d(0);  // stays 0 if the diagonal is missing, so singularity surfaces as inf/NaN
+        bool hasDiag = false;
+        for (typename Eval::InnerIterator dt(mat, j); dt; ++dt)
+          if (StorageIndex(dt.index()) == j) {
+            d = dt.value();
+            hasDiag = true;
+          }
+        eigen_assert(hasDiag && "sparse triangular solve: missing diagonal");
+        x[j] /= d;  // d == 0 when the diagonal is missing -> inf/NaN, consistent with the pointer path
+      }
+      Scalar xj = x[j];
+      for (typename Eval::InnerIterator it(mat, j); it && StorageIndex(it.index()) < j; ++it)
+        x[it.index()] = numext::madd<Scalar>(-xj, it.value(), x[it.index()]);
+    } else {
+      typename Eval::InnerIterator it(mat, j);
+      while (it && StorageIndex(it.index()) < j) ++it;     // skip out-of-triangle (index < j)
+      bool hasDiag = it && StorageIndex(it.index()) == j;  // diagonal is the first in-triangle entry
+      EIGEN_IF_CONSTEXPR (!UnitDiag) {
+        eigen_assert(hasDiag && "sparse triangular solve: missing diagonal");
+        // Missing diagonal -> inf/NaN, not it.value() on an ended/wrong iterator (see pointer path).
+        x[j] /= hasDiag ? it.value() : Scalar(0);
+      }
+      if (hasDiag) ++it;  // step past the stored diagonal (the divisor above, or a unit entry)
+      Scalar xj = x[j];
+      for (; it; ++it) x[it.index()] = numext::madd<Scalar>(-xj, it.value(), x[it.index()]);
+    }
+  }
+}
+
+// Iterator core: reach + numeric, returning top (no cleanup), the iterator counterpart
+// of reach_solve_dense -- xwork must already hold the scattered rhs. Uses the same 2n /
+// n(bytes) / n workspace layout (the pstack half of iwork is left unused), so the two
+// general paths share one workspace contract.
+template <bool Upper, bool UnitDiag, typename LhsType, typename StorageIndex, typename Scalar>
+Index reach_solve_dense_iter(const LhsType& lhs, Index n, const StorageIndex* bIdx, Index bCount, StorageIndex* iwork,
+                             uint8_t* mark, Scalar* xwork) {
+  evaluator<LhsType> mat(lhs);
+  Index top = triangular_reach_iter<Upper>(mat, bIdx, bCount, iwork, mark, n);
+  triangular_solve_over_reach_iter<Upper, UnitDiag>(mat, iwork, top, n, xwork);
+  return top;
+}
+
+// Policy dispatch for the core (returns top): an expression that exposes raw storage
+// takes the pointer + DFS fast path; anything else takes the evaluator + worklist path.
+// Tag dispatch (not if-constexpr) keeps the untaken branch from being instantiated, so
+// outerIndexPtr() is never named on a type that lacks it.
+//
+// CompressedAccessBit is a compile-time capability, not a guarantee the instance is
+// compressed: an uncompressed SparseMatrix keeps per-column gaps addressed via
+// innerNonZeroPtr(), so its columns do NOT run to outerIndexPtr()[j+1]. Passing
+// innerNonZeroPtr() through keeps the raw-pointer path valid either way -- it is nullptr
+// exactly when compressed (columns end at outerIndexPtr[j+1]) and the per-column count
+// otherwise (columns end at outerIndexPtr[j]+innerNonZeroPtr[j]).
+template <bool Upper, bool UnitDiag, typename LhsType, typename StorageIndex, typename Scalar>
+Index reach_solve_dense_dispatch(std::true_type /*compressed*/, const LhsType& lhs, Index n, const StorageIndex* bIdx,
+                                 Index bCount, StorageIndex* iwork, uint8_t* mark, Scalar* xwork) {
+  return reach_solve_dense<Upper, UnitDiag>(lhs.outerIndexPtr(), lhs.innerIndexPtr(), lhs.valuePtr(),
+                                            lhs.innerNonZeroPtr(), n, bIdx, bCount, iwork, mark, xwork);
+}
+template <bool Upper, bool UnitDiag, typename LhsType, typename StorageIndex, typename Scalar>
+Index reach_solve_dense_dispatch(std::false_type /*iterator*/, const LhsType& lhs, Index n, const StorageIndex* bIdx,
+                                 Index bCount, StorageIndex* iwork, uint8_t* mark, Scalar* xwork) {
+  return reach_solve_dense_iter<Upper, UnitDiag>(lhs, n, bIdx, bCount, iwork, mark, xwork);
+}
+
+// Expression core: solve T x = b for a sparse-expression triangular T with the rhs
+// PRE-SCATTERED into xwork (bIdx is its pattern), selecting the pointer or iterator path
+// at compile time; RETURNS top with the solution left in xwork and the reach in
+// iwork[top..n) (see reach_solve_dense). This is what the sparse selector uses -- it
+// scatters the rhs as it reads it and consumes xwork directly, so no bVal/outIdx/outVal.
+template <bool Upper, bool UnitDiag, typename LhsDerived, typename StorageIndex, typename Scalar>
+Index reach_solve_dense(const SparseMatrixBase<LhsDerived>& lhs, const StorageIndex* bIdx, Index bCount,
+                        StorageIndex* iwork, uint8_t* mark, Scalar* xwork) {
+  return reach_solve_dense_dispatch<Upper, UnitDiag>(
+      std::integral_constant<bool, has_compressed_access<LhsDerived>::value>{}, lhs.derived(), lhs.rows(), bIdx, bCount,
+      iwork, mark, xwork);
+}
+
+}  // namespace internal
+
+}  // namespace Eigen
+
+#endif  // EIGEN_TRIANGULAR_REACH_SOLVER_H

Eigen/src/SparseCore/TriangularSolver.h

@@ -182,64 +182,126 @@ template <typename Lhs, typename Rhs, int Mode,
           int StorageOrder = int(Lhs::Flags) & RowMajorBit>
 struct sparse_solve_triangular_sparse_selector;
 
+// True when the rhs exposes raw CSC storage with a StorageIndex matching the lhs, so a
+// column's stored index slice can serve as the reach roots directly (no bIdx copy). A
+// SparseVector qualifies too -- it is a single compressed column, handled below via the
+// null-outerIndexPtr guard (its outerIndexPtr() is null since it has no outer array).
+template <typename Lhs, typename Rhs>
+using rhs_matching_slice = std::integral_constant<
+    bool, has_compressed_access<Rhs>::value &&
+              std::is_same<typename traits<Rhs>::StorageIndex, typename traits<Lhs>::StorageIndex>::value>;
+
+// The reach for a column arrives in one of three compile-time-known orders, but the
+// output column must store ascending inner index. reach_reorder encapsulates the
+// per-order fix-up via partial specialization: an unordered reach (Ordered == false,
+// the pointer/DFS path -- Upper irrelevant) is sorted; the iterator reach is already
+// in solve order -- ascending for lower (no-op), descending for upper (reverse).
+template <bool Ordered, bool Upper, typename StorageIndex>
+struct reach_reorder {  // Ordered == false: unordered pointer/DFS reach
+  static void run(StorageIndex* first, StorageIndex* last) { std::sort(first, last); }
+};
+template <typename StorageIndex>
+struct reach_reorder<true, false, StorageIndex> {  // iterator reach, lower: already ascending
+  static void run(StorageIndex* /*first*/, StorageIndex* /*last*/) {}
+};
+template <typename StorageIndex>
+struct reach_reorder<true, true, StorageIndex> {  // iterator reach, upper: descending
+  static void run(StorageIndex* first, StorageIndex* last) { std::reverse(first, last); }
+};
+
+// Common per-column finish: reorder the reach xi[top..n) to ascending inner index,
+// insert reading values from xwork, and clear xwork and mark for the next column.
+// The reach is a structural bound, not a numeric one: a reached coefficient can be
+// exactly zero (a zero rhs entry, or numerical cancellation), so skip exact zeros at
+// insertion. This matches the AmbiVector path, which pruned zeros, and keeps a zero rhs
+// from materializing O(|reach|) stored zeros. xwork and mark are cleared regardless.
+template <bool Ordered, bool Upper, typename Res, typename StorageIndex, typename Scalar>
+void reach_insert_column(Res& res, Index col, StorageIndex* xi, Index top, Index n, Scalar* xwork, uint8_t* mark) {
+  reach_reorder<Ordered, Upper, StorageIndex>::run(xi + top, xi + n);
+  for (Index k = top; k < n; ++k) {
+    StorageIndex j = xi[k];
+    if (!numext::is_exactly_zero(xwork[j])) res.insert(j, col) = xwork[j];
+    xwork[j] = Scalar(0);
+    mark[j] = 0;
+  }
+}
+
+// Column loop, fast path: the rhs matches (see rhs_matching_slice), so each column's
+// stored index slice is the reach root list and the value slice is scattered directly
+// -- no bIdx copy, so iwork is just 2n (xi | pstack).
+template <bool Upper, bool UnitDiag, typename Lhs, typename Rhs, typename Res, typename Scalar,
+          std::enable_if_t<rhs_matching_slice<Lhs, Rhs>::value, int> = 0>
+void reach_solve_columns(const Lhs& lhs, const Rhs& other, Res& res, uint8_t* mark, Scalar* xwork, Index n) {
+  typedef typename traits<Lhs>::StorageIndex StorageIndex;
+  Matrix<StorageIndex, Dynamic, 1> iwork(2 * n);  // xi | pstack
+  StorageIndex* xi = iwork.data();
+  for (Index col = 0; col < other.cols(); ++col) {
+    const StorageIndex* outer = other.outerIndexPtr();  // null for a SparseVector (single column)
+    const StorageIndex* nnz = other.innerNonZeroPtr();  // null when compressed
+    Index p = outer ? outer[col] : 0;
+    Index bCount = outer ? (nnz ? Index(nnz[col]) : Index(outer[col + 1]) - p) : other.nonZeros();
+    const StorageIndex* roots = other.innerIndexPtr() + p;  // the column's stored indices
+    const Scalar* vals = other.valuePtr() + p;
+    // Roots are the rhs slice directly (no bIdx copy, so iwork stays 2n). An exact-zero
+    // stored rhs entry is seeded harmlessly -- it propagates zeros and is dropped at
+    // insertion; filtering it here would cost a compacted root buffer (the 3n path).
+    for (Index r = 0; r < bCount; ++r) xwork[roots[r]] = vals[r];
+    Index top = reach_solve_dense<Upper, UnitDiag>(lhs, roots, bCount, xi, mark, xwork);
+    reach_insert_column<!has_compressed_access<Lhs>::value, Upper>(res, col, xi, top, n, xwork, mark);
+  }
+}
+
+// Column loop, fallback: read each column through the InnerIterator, copying indices
+// into the bIdx third of a 3n iwork. For a rhs without raw storage or with a
+// mismatched index type.
+template <bool Upper, bool UnitDiag, typename Lhs, typename Rhs, typename Res, typename Scalar,
+          std::enable_if_t<!rhs_matching_slice<Lhs, Rhs>::value, int> = 0>
+void reach_solve_columns(const Lhs& lhs, const Rhs& other, Res& res, uint8_t* mark, Scalar* xwork, Index n) {
+  typedef typename traits<Lhs>::StorageIndex StorageIndex;
+  Matrix<StorageIndex, Dynamic, 1> iwork(3 * n);  // xi | pstack | bIdx
+  StorageIndex* xi = iwork.data();
+  StorageIndex* bIdx = iwork.data() + 2 * n;
+  for (Index col = 0; col < other.cols(); ++col) {
+    Index bCount = 0;
+    for (typename Rhs::InnerIterator it(other, col); it; ++it) {
+      if (numext::is_exactly_zero(it.value())) continue;  // a zero root seeds nothing; xwork stays clear there
+      bIdx[bCount] = StorageIndex(it.index());
+      xwork[it.index()] = it.value();
+      ++bCount;
+    }
+    Index top = reach_solve_dense<Upper, UnitDiag>(lhs, bIdx, bCount, xi, mark, xwork);
+    reach_insert_column<!has_compressed_access<Lhs>::value, Upper>(res, col, xi, top, n, xwork, mark);
+  }
+}
+
+// Reach-based (Gilbert-Peierls) sparse triangular solve, col-major, for lower OR
+// upper. Only the columns reachable from each rhs column's pattern are touched, so
+// the cost is O(|reach| + flops) per column instead of a dense O(n)-per-column sweep
+// (which also pays a coeff(i,i) binary search per row in the upper case). It is
+// the sole col-major sparse-sparse selector, dispatching lower/upper via the UpLo
+// template argument. reach_solve_dense leaves the solution values in
+// xwork and the reached indices in iwork[top..n); reach_solve_columns (slice or
+// fallback, selected on the rhs storage) scatters each column and solves, and
+// reach_insert_column reads the values out and restores mark/xwork. Only mark and
+// xwork need zeroing -- iwork is entirely written before read.
+template <bool Upper, typename Lhs, typename Rhs, int Mode>
+void run_sparse_reach_triangular_solve(const Lhs& lhs, Rhs& other) {
+  typedef typename Rhs::Scalar Scalar;
+  Index n = lhs.rows();
+  Matrix<uint8_t, Dynamic, 1> mark = Matrix<uint8_t, Dynamic, 1>::Zero(n);
+  Matrix<Scalar, Dynamic, 1> xwork = Matrix<Scalar, Dynamic, 1>::Zero(n);
+  Rhs res(other.rows(), other.cols());
+  res.reserve(other.nonZeros());
+  reach_solve_columns<Upper, bool(Mode & UnitDiag)>(lhs, other, res, mark.data(), xwork.data(), n);
+  res.finalize();
+  other = res.markAsRValue();
+}
+
 // forward and backward substitution, col-major
 template <typename Lhs, typename Rhs, int Mode, int UpLo>
 struct sparse_solve_triangular_sparse_selector<Lhs, Rhs, Mode, UpLo, ColMajor> {
-  typedef typename Rhs::Scalar Scalar;
-  typedef typename promote_index_type<typename traits<Lhs>::StorageIndex, typename traits<Rhs>::StorageIndex>::type
-      StorageIndex;
   static void run(const Lhs& lhs, Rhs& other) {
-    const bool IsLower = (UpLo == Lower);
-    AmbiVector<Scalar, StorageIndex> tempVector(other.rows() * 2);
-    tempVector.setBounds(0, other.rows());
-
-    Rhs res(other.rows(), other.cols());
-    res.reserve(other.nonZeros());
-
-    for (Index col = 0; col < other.cols(); ++col) {
-      // FIXME: estimate the number of non-zeros per column for better allocation.
-      tempVector.init(.99 /*float(other.col(col).nonZeros())/float(other.rows())*/);
-      tempVector.setZero();
-      tempVector.restart();
-      for (typename Rhs::InnerIterator rhsIt(other, col); rhsIt; ++rhsIt) {
-        tempVector.coeffRef(rhsIt.index()) = rhsIt.value();
-      }
-
-      for (Index i = IsLower ? 0 : lhs.cols() - 1; IsLower ? i < lhs.cols() : i >= 0; i += IsLower ? 1 : -1) {
-        tempVector.restart();
-        Scalar& ci = tempVector.coeffRef(i);
-        if (!numext::is_exactly_zero(ci)) {
-          // find
-          typename Lhs::InnerIterator it(lhs, i);
-          EIGEN_IF_CONSTEXPR (!(Mode & UnitDiag)) {
-            EIGEN_IF_CONSTEXPR (IsLower) {
-              eigen_assert(it.index() == i);
-              ci /= it.value();
-            } else
-              ci /= lhs.coeff(i, i);
-          }
-          tempVector.restart();
-          EIGEN_IF_CONSTEXPR (IsLower) {
-            if (it.index() == i) ++it;
-            for (; it; ++it) {
-              tempVector.coeffRef(it.index()) = numext::madd<Scalar>(-ci, it.value(), tempVector.coeffRef(it.index()));
-            }
-          } else {
-            for (; it && it.index() < i; ++it) {
-              tempVector.coeffRef(it.index()) = numext::madd<Scalar>(-ci, it.value(), tempVector.coeffRef(it.index()));
-            }
-          }
-        }
-      }
-
-      // FIXME: compute a reference value to filter zeros.
-      for (typename AmbiVector<Scalar, StorageIndex>::Iterator it(tempVector /*,1e-12*/); it; ++it) {
-        // FIXME: use insertBack for better performance.
-        res.insert(it.index(), col) = it.value();
-      }
-    }
-    res.finalize();
-    other = res.markAsRValue();
+    run_sparse_reach_triangular_solve<UpLo == Upper, Lhs, Rhs, Mode>(lhs, other);
   }
 };
 

test/sparse_solvers.cpp

@@ -91,6 +91,144 @@ void sparse_solvers(int rows, int cols) {
     m2.template triangularView<Upper>().solveInPlace(matB);
     VERIFY_IS_APPROX(matB, refMatB);
 
+    // A triangularView is a view of the triangular PART of a possibly-general matrix,
+    // so the stored matrix need not be strictly triangular. Exercise a general lhs, a
+    // SparseVector rhs, a mismatched-StorageIndex rhs, an uncompressed lhs, an
+    // expression lhs, and a unit diagonal -- none of which the checks above cover.
+    {
+      SparseMatrix<Scalar> mg(rows, rows);
+      DenseMatrix refMatG = DenseMatrix::Zero(rows, rows);
+      initSparse<Scalar>(density, refMatG, mg, ForceNonZeroDiag);  // GENERAL (both triangles stored)
+      initSparse<Scalar>(density, refMatB, matB);
+
+      // general matrix through a lower / upper / unit-upper view, sparse rhs
+      for (int mode = 0; mode < 3; ++mode) {
+        DenseMatrix rb = refMatB;
+        SparseMatrix<Scalar> mb = matB;
+        if (mode == 0) {
+          refMatG.template triangularView<Lower>().solveInPlace(rb);
+          mg.template triangularView<Lower>().solveInPlace(mb);
+        } else if (mode == 1) {
+          refMatG.template triangularView<Upper>().solveInPlace(rb);
+          mg.template triangularView<Upper>().solveInPlace(mb);
+        } else {
+          refMatG.template triangularView<UnitUpper>().solveInPlace(rb);
+          mg.template triangularView<UnitUpper>().solveInPlace(mb);
+        }
+        VERIFY_IS_APPROX(mb.toDense(), rb);
+      }
+
+      // expression lhs (no raw storage -> iterator path)
+      {
+        DenseMatrix rb = refMatB;
+        SparseMatrix<Scalar> mb = matB;
+        refMatG.template triangularView<Upper>().solveInPlace(rb);
+        (Scalar(1) * mg).template triangularView<Upper>().solveInPlace(mb);
+        VERIFY_IS_APPROX(mb.toDense(), rb);
+      }
+
+      // uncompressed lhs (innerNonZeroPtr != null)
+      {
+        DenseMatrix rb = refMatB;
+        SparseMatrix<Scalar> mb = matB, mu = mg;
+        mu.reserve(Matrix<int, Dynamic, 1>::Constant(mu.cols(), rows));  // -> uncompressed
+        refMatG.template triangularView<Lower>().solveInPlace(rb);
+        mu.template triangularView<Lower>().solveInPlace(mb);
+        VERIFY_IS_APPROX(mb.toDense(), rb);
+      }
+
+      // mismatched-StorageIndex rhs (-> InnerIterator fallback)
+      {
+        DenseMatrix rb = refMatB;
+        SparseMatrix<Scalar, ColMajor, long> mbl = matB;
+        refMatG.template triangularView<Lower>().solveInPlace(rb);
+        mg.template triangularView<Lower>().solveInPlace(mbl);
+        VERIFY_IS_APPROX(DenseMatrix(mbl), rb);
+      }
+
+      // SparseVector rhs (sets CompressedAccessBit but its outerIndexPtr() is null)
+      {
+        DenseVector rv = DenseVector::Zero(rows);
+        SparseVector<Scalar> vb(rows);
+        for (Index i = 0; i < rows; ++i)
+          if (internal::random<int>(0, 2) == 0) {
+            Scalar s = internal::random<Scalar>();
+            vb.coeffRef(i) = s;
+            rv(i) = s;
+          }
+        DenseVector rref = refMatG.template triangularView<Lower>().solve(rv);
+        SparseVector<Scalar> vx = vb;
+        mg.template triangularView<Lower>().solveInPlace(vx);
+        VERIFY_IS_APPROX(DenseVector(vx), rref);
+      }
+
+      // explicitly-stored zero rhs entries must not expand into stored zeros: the reach
+      // is a structural bound, so a zero rhs coefficient (or one that cancels to zero)
+      // is pruned at insertion rather than materialized across the whole reach.
+      {
+        SparseMatrix<Scalar> mb(rows, matB.cols());
+        DenseMatrix rb = DenseMatrix::Zero(rows, matB.cols());
+        for (Index c = 0; c < mb.cols(); ++c)
+          for (Index i = 0; i < rows; ++i)
+            if (internal::random<int>(0, 3) == 0) {
+              Scalar s = internal::random<int>(0, 2) == 0 ? Scalar(0) : internal::random<Scalar>();  // some explicit 0
+              mb.insert(i, c) = s;
+              rb(i, c) = s;
+            }
+        mb.makeCompressed();
+        refMatG.template triangularView<Lower>().solveInPlace(rb);
+        mg.template triangularView<Lower>().solveInPlace(mb);
+        VERIFY_IS_APPROX(mb.toDense(), rb);
+        for (Index c = 0; c < mb.cols(); ++c)
+          for (typename SparseMatrix<Scalar>::InnerIterator it(mb, c); it; ++it)
+            VERIFY(!numext::is_exactly_zero(it.value()));  // no stored zeros
+      }
+
+      // A reached column with no stored diagonal (non-unit) is out of contract: it must
+      // assert in debug on every path (pointer/iterator x lower/upper), rather than the
+      // failure being silently keyed to has_compressed_access. In release these divide by
+      // zero -> inf/NaN with no out-of-bounds read (covered by the sanitizer drivers).
+      {
+        SparseMatrix<Scalar> us(3, 3);
+        us.insert(0, 1) = Scalar(1);
+        us.insert(1, 1) = Scalar(2);
+        us.insert(2, 2) = Scalar(3);
+        us.makeCompressed();  // column 0 empty -> reached from rhs(0) but no diagonal
+        SparseMatrix<Scalar> ub(3, 1);
+        ub.insert(0, 0) = Scalar(1);
+        ub.makeCompressed();
+        SparseMatrix<Scalar> up = ub, ui = ub;
+        VERIFY_RAISES_ASSERT(us.template triangularView<Upper>().solveInPlace(up));                // pointer upper
+        VERIFY_RAISES_ASSERT((Scalar(1) * us).template triangularView<Upper>().solveInPlace(ui));  // iterator upper
+
+        SparseMatrix<Scalar> ls(3, 3);
+        ls.insert(0, 0) = Scalar(2);
+        ls.insert(1, 1) = Scalar(3);
+        ls.makeCompressed();  // column 2 empty (last) -> reached from rhs(2) but no diagonal
+        SparseMatrix<Scalar> lb(3, 1);
+        lb.insert(2, 0) = Scalar(1);
+        lb.makeCompressed();
+        SparseMatrix<Scalar> lp = lb, li = lb;
+        VERIFY_RAISES_ASSERT(ls.template triangularView<Lower>().solveInPlace(lp));                // pointer lower
+        VERIFY_RAISES_ASSERT((Scalar(1) * ls).template triangularView<Lower>().solveInPlace(li));  // iterator lower
+      }
+
+      // mixed-scalar rhs: a real lhs applied to a rhs must accumulate in the rhs scalar.
+      // For real Scalar this is the ordinary path; for complex Scalar it is the
+      // real-factor / complex-data case that must both compile and be correct.
+      {
+        typedef typename NumTraits<Scalar>::Real Real;
+        SparseMatrix<Real> mr(rows, rows);
+        Matrix<Real, Dynamic, Dynamic> refMatR = Matrix<Real, Dynamic, Dynamic>::Zero(rows, rows);
+        initSparse<Real>(density, refMatR, mr, ForceNonZeroDiag);
+        DenseMatrix rb = refMatB;
+        SparseMatrix<Scalar> mb = matB;
+        refMatR.template cast<Scalar>().template triangularView<Lower>().solveInPlace(rb);
+        mr.template triangularView<Lower>().solveInPlace(mb);  // SparseMatrix<Real> lhs, <Scalar> rhs
+        VERIFY_IS_APPROX(mb.toDense(), rb);
+      }
+    }
+
     // test deprecated API
     initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
     VERIFY_IS_APPROX(refMat2.template triangularView<Lower>().solve(vec2),