Add ELLSparsityPattern

include/NeoN/linearAlgebra/cooSparsityPattern.hpp

@@ -66,17 +66,24 @@ class CooSparsityPattern : public NeoN::SupportsCopyTo<CooSparsityPattern<IndexT
 
     [[nodiscard]] localIdx rows() const { return dimensions_.rows; };
 
+    /*@brief true count of nonzero matrix entries (no padding in COO) */
     [[nodiscard]] localIdx nnz() const { return colIdxs_.size(); };
 
+    /*@brief size of the backing storage, identical to nnz() since COO has no padding */
+    [[nodiscard]] localIdx storageSize() const { return colIdxs_.size(); };
+
     [[nodiscard]] Dimensions dimension() const { return dimensions_; };
 
+    using ViewType = SparsityView<IndexType>;
+
     /**
      * @brief Get a view representation of the matrix's data.
-     * @return MatrixView for easy access to matrix elements.
+     * @return MatrixView for easy access to matrix elements. Backed by rowOffs_,
+     * not rowIdxs_ -- entry()/findEntry() need row-range offsets.
      */
     [[nodiscard]] SparsityView<IndexType> view() const
     {
-        return SparsityView<IndexType>(colIdxs_.view(), rowIdxs_.view());
+        return SparsityView<IndexType>(colIdxs_.view(), rowOffs_.view());
     }
 
 private:

include/NeoN/linearAlgebra/csrSparsityPattern.hpp

@@ -58,9 +58,16 @@ class CsrSparsityPattern : public NeoN::SupportsCopyTo<CsrSparsityPattern<IndexT
 
     [[nodiscard]] localIdx rows() const { return dimensions_.rows; };
 
+    /*@brief true count of nonzero matrix entries (no padding in CSR) */
     [[nodiscard]] localIdx nnz() const { return colIdxs_.size(); };
 
+    /*@brief size of the backing storage, identical to nnz() since CSR has no padding */
+    [[nodiscard]] localIdx storageSize() const { return colIdxs_.size(); };
+
     [[nodiscard]] Dimensions dimension() const { return dimensions_; };
+
+    using ViewType = SparsityView<IndexType>;
+
     /**
      * @brief Get a view representation of the matrix's data.
      * @return MatrixView for easy access to matrix elements.

include/NeoN/linearAlgebra/ellSparsityPattern.hpp

@@ -0,0 +1,123 @@
+// SPDX-FileCopyrightText: 2025 - 2026 NeoN authors
+//
+// SPDX-License-Identifier: MIT
+
+#pragma once
+
+#include "NeoN/core/copyTo.hpp"
+#include "NeoN/core/vector/vector.hpp"
+#include "NeoN/mesh/unstructured/unstructuredMesh.hpp"
+#include "NeoN/linearAlgebra/sparsityView.hpp"
+
+namespace NeoN::la
+{
+
+/** @class EllSparsityPattern
+ * @brief fixed-width, padded column-index representation of a matrix (ELLPACK format)
+ *
+ * Every row reserves `numStoredElementsPerRow()` column-index slots, the width of
+ * the widest row in the matrix. Rows with fewer nonzeros are padded with
+ * `EllSparsityView<IndexType>::invalidIndex()`. Slots are stored column-major,
+ * i.e. slot `s` of row `i` is stored at `i + stride() * s`, so that a
+ * one-thread-per-row kernel reads/writes memory in a coalesced fashion when
+ * iterating over slots.
+ */
+template<typename IndexType>
+class EllSparsityPattern : public NeoN::SupportsCopyTo<EllSparsityPattern<IndexType>>
+{
+
+    void validate() const;
+
+public:
+
+    using SparsityIndexType = IndexType;
+
+    /* @brief create a copy of a given EllSparsityPattern */
+    EllSparsityPattern(const EllSparsityPattern& sp);
+
+    /**
+     * @brief construct from a fully-built padded, column-major colIdx array.
+     * @param colIdx per-slot column indices, sorted ascending within each row's slots,
+     * with padding (EllSparsityView<IndexType>::invalidIndex()) trailing. Not verified.
+     * @param logicalNnz count of non-padding entries in colIdx. Caller-supplied and
+     * trusted, not recomputed from colIdx -- only checked against storage size.
+     */
+    EllSparsityPattern(
+        Vector<IndexType>&& colIdx,
+        Dimensions dim,
+        localIdx numStoredElementsPerRow,
+        localIdx stride,
+        localIdx logicalNnz
+    );
+
+    [[nodiscard]] EllSparsityPattern copyToExecutor(Executor dstExec) const override
+    {
+        return EllSparsityPattern<IndexType>(
+            colIdxs_.copyToExecutor(dstExec),
+            dimensions_,
+            numStoredElementsPerRow_,
+            stride_,
+            logicalNnz_
+        );
+    }
+
+    ~EllSparsityPattern() = default;
+
+    /*@brief getter for executor */
+    const Executor& exec() const { return colIdxs_.exec(); }
+
+    /*@brief const-only getter for colIdxs -- ELL patterns are immutable after construction
+     * so nnz() can't desync from the stored columns */
+    [[nodiscard]] const Vector<IndexType>& colIdxs() const { return colIdxs_; };
+
+    [[nodiscard]] localIdx rows() const { return dimensions_.rows; };
+
+    /**
+     * @brief size of the (padded) storage backing this pattern, i.e. the size
+     * `values_` of a `Matrix` built on top of this pattern must have.
+     * @note unlike CSR, this includes padding entries and is therefore not the
+     * count of logical/non-zero matrix entries -- see nnz().
+     */
+    [[nodiscard]] localIdx storageSize() const { return colIdxs_.size(); };
+
+    /*@brief true count of logical (non-padding) nonzero matrix entries */
+    [[nodiscard]] localIdx nnz() const { return logicalNnz_; };
+
+    /*@brief number of column-index slots stored per row, including padding */
+    [[nodiscard]] localIdx numStoredElementsPerRow() const { return numStoredElementsPerRow_; };
+
+    /*@brief stride between successive slots of the column-major-stored ELL arrays */
+    [[nodiscard]] localIdx stride() const { return stride_; };
+
+    [[nodiscard]] Dimensions dimension() const { return dimensions_; };
+
+    using ViewType = EllSparsityView<IndexType>;
+
+    /**
+     * @brief Get a view representation of the matrix's data.
+     * @return EllSparsityView for easy access to matrix elements.
+     */
+    [[nodiscard]] EllSparsityView<IndexType> view() const
+    {
+        return EllSparsityView<IndexType>(
+            colIdxs_.view(),
+            static_cast<IndexType>(numStoredElementsPerRow_),
+            static_cast<IndexType>(stride_)
+        );
+    }
+
+private:
+
+    Dimensions dimensions_;
+
+    Vector<IndexType> colIdxs_; //! padded, column-major column indices,
+                                //! size stride_ * numStoredElementsPerRow_
+
+    localIdx numStoredElementsPerRow_; //! width of the widest row, i.e. slots stored per row
+
+    localIdx stride_; //! distance between slot s and slot s+1 of the same row
+
+    localIdx logicalNnz_; //! true count of non-padding nonzero entries, <= storageSize()
+};
+
+} // namespace NeoN::la

include/NeoN/linearAlgebra/linearSystem.hpp

@@ -398,7 +398,8 @@ LinearSystem<ValueType, RHSValueType, SystemMatrixType, BoundaryMatrixType> crea
  *
  * @note templated on the full LinearSystem parameter set so it also accepts the segregated
  * vector-solve form (scalar matrix, Vec3 rhs): the scalar boundary diagonal is reversed on the
- * scalar matrix while the rhs reversal uses the field (RHS) value type. **/
+ * scalar matrix while the rhs reversal uses the field (RHS) value type.
+ * @note BoundaryMatrixType must expose per-entry row indices through rowIdxs() (currently COO). **/
 template<
     typename MatrixValueType,
     typename RHSValueType,
@@ -409,6 +410,7 @@ removeBoundaryContributions(
     const la::LinearSystem<MatrixValueType, RHSValueType, SystemMatrixType, BoundaryMatrixType>&
         lsIn
 )
+    requires requires(const BoundaryMatrixType& bm) { bm.sparsity()->rowIdxs(); }
 {
     auto ls =
         la::LinearSystem<MatrixValueType, RHSValueType, SystemMatrixType, BoundaryMatrixType>(lsIn);
@@ -419,12 +421,14 @@ removeBoundaryContributions(
     auto& bRhs = lsView.boundaryRhs;
 
     const auto ma = ls.faceToMatrixAddress()->view(ls.matrix().sparsity()->rowOffs().view());
+    // Per-face owner cell -- bMatrix.sparsity.rowOffs is CSR-style range data, not a cell index.
+    const auto bRowIdxs = ls.boundaryMatrix().sparsity()->rowIdxs().view();
 
     parallelFor(
         ls.exec(),
         {0, bMatrix.values.size()},
         NEON_LAMBDA(const localIdx facei) {
-            const auto celli = bMatrix.sparsity.rowOffs[facei]; // cell index stored in rowOffs
+            const auto celli = bRowIdxs[facei];
             Kokkos::atomic_add(&matrix.values[ma.diagIdx(celli)], bMatrix.values[facei]);
             Kokkos::atomic_add(&rhs[celli], bRhs[facei]);
         },

include/NeoN/linearAlgebra/matrix.hpp

@@ -7,6 +7,7 @@
 #include "NeoN/core/vector/vector.hpp"
 #include "NeoN/linearAlgebra/cooSparsityPattern.hpp"
 #include "NeoN/linearAlgebra/csrSparsityPattern.hpp"
+#include "NeoN/linearAlgebra/ellSparsityPattern.hpp"
 #include "NeoN/linearAlgebra/faceToMatrixAddress.hpp"
 
 namespace NeoN::la
@@ -56,15 +57,15 @@ struct MatrixView
     KOKKOS_INLINE_FUNCTION
     ValueType& entry(const localIdx offset) const { return values[offset]; }
 
-    View<ValueType> values; //!< View to the values of the CSR matrix.
+    View<ValueType> values; //!< View to the non-zero values of the matrix.
     SparsityViewType sparsity;
 };
 
 /**
  * @class Matrix
- * @brief Sparse matrix class with compact storage by row (CSR) format.
+ * @brief Sparse matrix class, generic over its sparsity storage format (CSR, COO, ELL, ...).
  * @tparam ValueType The value type of the non-zero entries.
- * @tparam IndexType The index type of the rows and columns.
+ * @tparam SparsityType The sparsity pattern type backing this matrix, e.g. CsrSparsityPattern.
  */
 template<typename ValueType, typename SparsityType>
 class Matrix : public NeoN::SupportsCopyTo<Matrix<ValueType, SparsityType>>
@@ -73,8 +74,9 @@ class Matrix : public NeoN::SupportsCopyTo<Matrix<ValueType, SparsityType>>
     void validate()
     {
         NF_ASSERT(values_.exec() == sparsityPattern_->exec(), "Executors are not the same");
-        // TODO this is not necessarily true for matrix types with padding like ELL
-        NF_ASSERT(values_.size() == sparsityPattern_->nnz(), "Matrix values and columns mismatch");
+        NF_ASSERT(
+            values_.size() == sparsityPattern_->storageSize(), "Matrix values and columns mismatch"
+        );
     }
 
 public:
@@ -94,7 +96,8 @@ class Matrix : public NeoN::SupportsCopyTo<Matrix<ValueType, SparsityType>>
     }
 
     /**
-     * @brief Constructor for Matrix.
+     * @brief Constructor for Matrix from a CSR/COO-shaped (colIdxs, rowOffs) pair.
+     * Not available for ELL -- use the two-argument constructor instead.
      * @param values The non-zero values of the matrix.
      * @param colIdxs The column indices for each non-zero value.
      * @param rowOffs The starting index in values/colIdxs for each row.
@@ -105,6 +108,11 @@ class Matrix : public NeoN::SupportsCopyTo<Matrix<ValueType, SparsityType>>
         const Vector<typename SparsityType::SparsityIndexType>& rowOffs,
         Dimensions dimensions
     )
+        requires requires(
+                     Vector<typename SparsityType::SparsityIndexType> c,
+                     Vector<typename SparsityType::SparsityIndexType> r,
+                     Dimensions d
+                 ) { SparsityType(std::move(c), std::move(r), d); }
         : values_(values),
           sparsityPattern_(
               std::make_shared<const SparsityType>(Vector(colIdxs), Vector(rowOffs), dimensions)
@@ -117,7 +125,8 @@ class Matrix : public NeoN::SupportsCopyTo<Matrix<ValueType, SparsityType>>
      * @brief Constructor for Matrix with a FaceToMatrixAddress.
      *
      * The sparsity pattern and the face-to-matrix address are passed as independent objects.
-     * Only available for matrices whose index type is localIdx.
+     * Only available for localIdx-indexed, rowOffs()-capable sparsity types -- FaceToMatrixAddress
+     * offsets are CSR row-local and don't apply to ELL.
      *
      * @param values The non-zero values of the matrix.
      * @param sparsity The sparsity pattern of the matrix.
@@ -128,7 +137,8 @@ class Matrix : public NeoN::SupportsCopyTo<Matrix<ValueType, SparsityType>>
         std::shared_ptr<const SparsityType> sparsity,
         std::shared_ptr<const FaceToMatrixAddress> faceToMatrixAddress
     )
-        requires std::is_same_v<typename SparsityType::SparsityIndexType, localIdx>
+        requires(std::is_same_v<typename SparsityType::SparsityIndexType, localIdx>
+                 && requires(const SparsityType& sp) { sp.rowOffs(); })
         : values_(values), sparsityPattern_(sparsity), faceToMatrixAddress_(faceToMatrixAddress)
     {
         validate();
@@ -176,10 +186,12 @@ class Matrix : public NeoN::SupportsCopyTo<Matrix<ValueType, SparsityType>>
     }
 
     /**
-     * @brief Get a reference to row offset vector.
+     * @brief Get a reference to row offset vector. Not available for ELL -- no
+     * contiguous row range to offset into.
      * @return Vector containing the row pointers.
      */
     [[nodiscard]] const Vector<typename SparsityType::SparsityIndexType>& rowOffs() const
+        requires requires(const SparsityType& sp) { sp.rowOffs(); }
     {
         return sparsityPattern_->rowOffs();
     }
@@ -215,10 +227,9 @@ class Matrix : public NeoN::SupportsCopyTo<Matrix<ValueType, SparsityType>>
      * @brief Get a view representation of the matrix's data.
      * @return MatrixView for easy access to matrix elements.
      */
-    [[nodiscard]] MatrixView<ValueType, SparsityView<typename SparsityType::SparsityIndexType>>
-    view()
+    [[nodiscard]] MatrixView<ValueType, typename SparsityType::ViewType> view()
     {
-        return MatrixView<ValueType, SparsityView<typename SparsityType::SparsityIndexType>>(
+        return MatrixView<ValueType, typename SparsityType::ViewType>(
             values_.view(), sparsityPattern_->view()
         );
     }
@@ -227,12 +238,9 @@ class Matrix : public NeoN::SupportsCopyTo<Matrix<ValueType, SparsityType>>
      * @brief Get a const view representation of the matrix's data.
      * @return Const MatrixView for read-only access to matrix elements.
      */
-    [[nodiscard]] MatrixView<
-        const ValueType,
-        SparsityView<typename SparsityType::SparsityIndexType>>
-    view() const
+    [[nodiscard]] MatrixView<const ValueType, typename SparsityType::ViewType> view() const
     {
-        return MatrixView<const ValueType, SparsityView<typename SparsityType::SparsityIndexType>>(
+        return MatrixView<const ValueType, typename SparsityType::ViewType>(
             View<const ValueType>(values_.view()), sparsityPattern_->view()
         );
     }
@@ -245,7 +253,7 @@ class Matrix : public NeoN::SupportsCopyTo<Matrix<ValueType, SparsityType>>
 
 private:
 
-    Vector<ValueType> values_; //!< The (non-zero) values of the CSR matrix.
+    Vector<ValueType> values_; //!< The (non-zero) values of the matrix.
 
     std::shared_ptr<const SparsityType> sparsityPattern_;
 
@@ -259,6 +267,9 @@ using CSRMatrix = Matrix<ValueType, la::CsrSparsityPattern<IndexType>>;
 template<typename ValueType, typename IndexType>
 using COOMatrix = Matrix<ValueType, la::CooSparsityPattern<IndexType>>;
 
+template<typename ValueType, typename IndexType>
+using ELLMatrix = Matrix<ValueType, la::EllSparsityPattern<IndexType>>;
+
 /** @brief extract the upper triangular of the matrix
  * @note this function is meant for testing purposes, it will recompute upper offsets
  */