From ef47c0cee7fb3b543f90bb49c224ae5667298ec9 Mon Sep 17 00:00:00 2001 From: Dheeraj Raghunathan <64003931+dheerajraghunathan@users.noreply.github.com> Date: Tue, 14 Jul 2026 10:37:44 +0200 Subject: [PATCH 1/3] Add ELLSparsityPattern --- .../NeoN/linearAlgebra/ellSparsityPattern.hpp | 105 +++++++++++++++++ include/NeoN/linearAlgebra/sparsityView.hpp | 63 +++++++++++ src/CMakeLists.txt | 1 + src/linearAlgebra/ellSparsityPattern.cpp | 43 +++++++ test/linearAlgebra/CMakeLists.txt | 1 + test/linearAlgebra/ellSparsityPattern.cpp | 106 ++++++++++++++++++ 6 files changed, 319 insertions(+) create mode 100644 include/NeoN/linearAlgebra/ellSparsityPattern.hpp create mode 100644 src/linearAlgebra/ellSparsityPattern.cpp create mode 100644 test/linearAlgebra/ellSparsityPattern.cpp diff --git a/include/NeoN/linearAlgebra/ellSparsityPattern.hpp b/include/NeoN/linearAlgebra/ellSparsityPattern.hpp new file mode 100644 index 0000000000..cbb7a2c5ee --- /dev/null +++ b/include/NeoN/linearAlgebra/ellSparsityPattern.hpp @@ -0,0 +1,105 @@ +// 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::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 +class EllSparsityPattern : public NeoN::SupportsCopyTo> +{ + + void validate() const; + +public: + + using SparsityIndexType = IndexType; + + /* @brief create a copy of a given EllSparsityPattern */ + EllSparsityPattern(const EllSparsityPattern& sp); + + EllSparsityPattern( + Vector&& colIdx, + Dimensions dim, + localIdx numStoredElementsPerRow, + localIdx stride + ); + + [[nodiscard]] EllSparsityPattern copyToExecutor(Executor dstExec) const override + { + return EllSparsityPattern( + colIdxs_.copyToExecutor(dstExec), dimensions_, numStoredElementsPerRow_, stride_ + ); + } + + ~EllSparsityPattern() = default; + + /*@brief getter for executor */ + const Executor& exec() const { return colIdxs_.exec(); } + + /*@brief getter for colIdxs */ + [[nodiscard]] const Vector& colIdxs() const { return colIdxs_; }; + + [[nodiscard]] Vector& colIdxs() { 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. + */ + [[nodiscard]] localIdx nnz() const { return colIdxs_.size(); }; + + /*@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_; }; + + /** + * @brief Get a view representation of the matrix's data. + * @return EllSparsityView for easy access to matrix elements. + */ + [[nodiscard]] EllSparsityView view() const + { + return EllSparsityView( + colIdxs_.view(), + static_cast(numStoredElementsPerRow_), + static_cast(stride_) + ); + } + +private: + + Dimensions dimensions_; + + Vector 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 +}; + +} // namespace NeoN::la diff --git a/include/NeoN/linearAlgebra/sparsityView.hpp b/include/NeoN/linearAlgebra/sparsityView.hpp index 2670678a03..5d2ab4bdbd 100644 --- a/include/NeoN/linearAlgebra/sparsityView.hpp +++ b/include/NeoN/linearAlgebra/sparsityView.hpp @@ -4,6 +4,8 @@ #pragma once +#include + #include "NeoN/core/view.hpp" namespace NeoN::la @@ -60,4 +62,65 @@ struct SparsityView View rowOffs; }; +/** + * @struct EllSparsityView + * @brief A view struct to allow easy read/write on all executors for the ELLPACK + * (fixed-width, padded) sparsity pattern. + * + * Column indices are stored column-major across the padded row slots, i.e. slot + * `s` of row `i` lives at flat offset `i + stride * s`. Rows with fewer nonzeros + * than `numStoredElementsPerRow` are padded with `invalidIndex()`. + * + * @tparam IndexType The index type of the rows and columns. + * @todo ideally this should be immutable + */ +template +struct EllSparsityView +{ + EllSparsityView( + View colIdxsView, IndexType numStoredElementsPerRowIn, IndexType strideIn + ) + : colIdxs(colIdxsView), numStoredElementsPerRow(numStoredElementsPerRowIn), + stride(strideIn) {}; + + /** + * @brief Sentinel column index marking an unused (padding) slot. + */ + KOKKOS_INLINE_FUNCTION + static constexpr IndexType invalidIndex() { return std::numeric_limits::max(); } + + /** + * @brief Flat storage offset of slot `slot` of row `i` (column-major layout). + */ + KOKKOS_INLINE_FUNCTION + IndexType linearIndex(const IndexType i, const IndexType slot) const + { + return i + stride * slot; + } + + /** + * @brief Retrieve the storage offset of the matrix element at position (i,j). + * @param i The row index. + * @param j The column index. + * @return Flat offset into colIdxs/values if it exists. + */ + KOKKOS_INLINE_FUNCTION + IndexType entry(const IndexType i, const IndexType j) const + { + for (std::remove_const_t slot = 0; slot < numStoredElementsPerRow; ++slot) + { + const IndexType idx = linearIndex(i, static_cast(slot)); + const IndexType col = colIdxs[idx]; + if (col == j) return idx; + if (col == invalidIndex() || col > j) break; + } + Kokkos::abort("Memory not allocated for ELL matrix component."); + return 0; // compiler warning suppression. + } + + View colIdxs; + IndexType numStoredElementsPerRow; + IndexType stride; +}; + } // namespace NeoN::la diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index be94af96c5..6844b10b31 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -44,6 +44,7 @@ target_sources( "linearAlgebra/faceToMatrixAddress.cpp" "linearAlgebra/cooSparsityPattern.cpp" "linearAlgebra/csrSparsityPattern.cpp" + "linearAlgebra/ellSparsityPattern.cpp" "linearAlgebra/meshIterationStrategies.cpp" "mesh/unstructured/boundaryMesh.cpp" "mesh/unstructured/unstructuredMesh.cpp" diff --git a/src/linearAlgebra/ellSparsityPattern.cpp b/src/linearAlgebra/ellSparsityPattern.cpp new file mode 100644 index 0000000000..443fed2234 --- /dev/null +++ b/src/linearAlgebra/ellSparsityPattern.cpp @@ -0,0 +1,43 @@ +// SPDX-FileCopyrightText: 2025 - 2026 NeoN authors +// +// SPDX-License-Identifier: MIT + +#include "NeoN/core/macros.hpp" +#include "NeoN/core/containerFreeFunctions.hpp" +#include "NeoN/core/array.hpp" +#include "NeoN/linearAlgebra/ellSparsityPattern.hpp" + +namespace NeoN::la +{ + +template +void EllSparsityPattern::validate() const +{ + NF_ASSERT( + colIdxs_.size() == stride_ * numStoredElementsPerRow_, + "ELL colIdxs size does not match stride * numStoredElementsPerRow" + ); + NF_ASSERT(stride_ >= dimensions_.rows, "ELL stride must be at least the number of rows"); +} + +template +EllSparsityPattern::EllSparsityPattern( + Vector&& colIdx, Dimensions dim, localIdx numStoredElementsPerRow, localIdx stride +) + : dimensions_(dim), colIdxs_(std::move(colIdx)), + numStoredElementsPerRow_(numStoredElementsPerRow), stride_(stride) +{ + validate(); +} + +template +EllSparsityPattern::EllSparsityPattern(const EllSparsityPattern& sp) + : dimensions_(sp.dimensions_), colIdxs_(sp.colIdxs_), + numStoredElementsPerRow_(sp.numStoredElementsPerRow_), stride_(sp.stride_) +{} + +#define NN_DECLARE_SPARSITY(TYPENAME) template class EllSparsityPattern; + +NN_FOR_ALL_INTEGER_TYPES(NN_DECLARE_SPARSITY); + +} // namespace NeoN::la diff --git a/test/linearAlgebra/CMakeLists.txt b/test/linearAlgebra/CMakeLists.txt index 84b081c8b5..59e4ba756a 100644 --- a/test/linearAlgebra/CMakeLists.txt +++ b/test/linearAlgebra/CMakeLists.txt @@ -7,6 +7,7 @@ neon_unit_test(matrixIterator) neon_unit_test(linearSystem) neon_unit_test(cooSparsityPattern) neon_unit_test(csrSparsityPattern) +neon_unit_test(ellSparsityPattern) neon_unit_test(faceToMatrixAddress) neon_unit_test(utilities) diff --git a/test/linearAlgebra/ellSparsityPattern.cpp b/test/linearAlgebra/ellSparsityPattern.cpp new file mode 100644 index 0000000000..d34367b8a8 --- /dev/null +++ b/test/linearAlgebra/ellSparsityPattern.cpp @@ -0,0 +1,106 @@ +// SPDX-FileCopyrightText: 2025 - 2026 NeoN authors +// +// SPDX-License-Identifier: MIT + +#define CATCH_CONFIG_RUNNER // Define this before including catch.hpp to create + // a custom main + +#include + +#include "catch2_common.hpp" + +#include "NeoN/NeoN.hpp" + +namespace NeoN +{ + +TEST_CASE("EllSparsityPattern") +{ + using EllSparsityType = NeoN::la::EllSparsityPattern; + auto [execName, exec] = GENERATE(allAvailableExecutor()); + + // clang-format off + // Matrix: + // Row/ColIdx 0 1 2 3 + // 0 x x + // 1 x x x + // 2 x x x + // 3 x x + // + // widest row has 3 entries -> numStoredElementsPerRow = 3 + // nRows = 4 -> stride = 4 (no extra padding rows) + // + // column-major, padded layout (INV marks unused slots): + // slot 0: [0, 0, 1, 2] + // slot 1: [1, 1, 2, 3] + // slot 2: [INV, 2, 3, INV] + // clang-format on + const auto INV = std::numeric_limits::max(); + const localIdx nRows = 4; + const localIdx numStoredElementsPerRow = 3; + const localIdx stride = nRows; + + auto colIdxExp = std::vector { + 0, + 0, + 1, + 2, // slot 0 + 1, + 1, + 2, + 3, // slot 1 + INV, + 2, + 3, + INV // slot 2 + }; + + Vector colIdx(exec, colIdxExp); + auto sp = std::make_shared( + std::move(colIdx), la::Dimensions {nRows, nRows}, numStoredElementsPerRow, stride + ); + + SECTION("Can store padded colIdxs and dimensions " + execName) + { + REQUIRE_THAT(sp->colIdxs(), Equals(colIdxExp, EqualInt())); + REQUIRE(sp->rows() == nRows); + REQUIRE(sp->numStoredElementsPerRow() == numStoredElementsPerRow); + REQUIRE(sp->stride() == stride); + REQUIRE(sp->nnz() == stride * numStoredElementsPerRow); + REQUIRE(sp->dimension().rows == nRows); + REQUIRE(sp->dimension().cols == nRows); + } + + SECTION("Can copy to executor " + execName) + { + auto spOnHost = sp->copyToExecutor(SerialExecutor()); + REQUIRE_THAT(spOnHost.colIdxs(), Equals(colIdxExp, EqualInt())); + REQUIRE(spOnHost.numStoredElementsPerRow() == numStoredElementsPerRow); + REQUIRE(spOnHost.stride() == stride); + } + + SECTION("Can resolve entry offsets on " + execName) + { + Vector checkOffset(exec, 4); + auto checkOffsetView = checkOffset.view(); + auto ellView = sp->view(); + parallelFor( + exec, + {0, 1}, + NEON_LAMBDA(const localIdx) { + checkOffsetView[0] = ellView.entry(0, 0); // row 0, slot 0 -> idx 0 + 4*0 = 0 + checkOffsetView[1] = ellView.entry(1, 2); // row 1, slot 2 -> idx 1 + 4*2 = 9 + checkOffsetView[2] = ellView.entry(2, 3); // row 2, slot 2 -> idx 2 + 4*2 = 10 + checkOffsetView[3] = ellView.entry(3, 3); // row 3, slot 1 -> idx 3 + 4*1 = 7 + } + ); + auto checkOffsetHost = checkOffset.copyToHost(); + auto checkOffsetHostView = checkOffsetHost.view(); + REQUIRE(checkOffsetHostView[0] == 0); + REQUIRE(checkOffsetHostView[1] == 9); + REQUIRE(checkOffsetHostView[2] == 10); + REQUIRE(checkOffsetHostView[3] == 7); + } +} + +} From cdddfa407ce9aee94ba4e970931eba1b7c401062 Mon Sep 17 00:00:00 2001 From: Dheeraj Raghunathan <64003931+dheerajraghunathan@users.noreply.github.com> Date: Thu, 16 Jul 2026 11:47:20 +0200 Subject: [PATCH 2/3] Update Matrix for ELL support --- .../NeoN/linearAlgebra/cooSparsityPattern.hpp | 11 +- .../NeoN/linearAlgebra/csrSparsityPattern.hpp | 7 + .../NeoN/linearAlgebra/ellSparsityPattern.hpp | 32 +++- include/NeoN/linearAlgebra/matrix.hpp | 47 +++-- include/NeoN/linearAlgebra/sparsityView.hpp | 58 ++++++- src/linearAlgebra/ellSparsityPattern.cpp | 12 +- src/linearAlgebra/matrix.cpp | 52 ++++-- test/linearAlgebra/cooSparsityPattern.cpp | 28 +++ test/linearAlgebra/ellSparsityPattern.cpp | 162 +++++++++++++++++- 9 files changed, 352 insertions(+), 57 deletions(-) diff --git a/include/NeoN/linearAlgebra/cooSparsityPattern.hpp b/include/NeoN/linearAlgebra/cooSparsityPattern.hpp index 0efed71b56..1a181156ad 100644 --- a/include/NeoN/linearAlgebra/cooSparsityPattern.hpp +++ b/include/NeoN/linearAlgebra/cooSparsityPattern.hpp @@ -66,17 +66,24 @@ class CooSparsityPattern : public NeoN::SupportsCopyTo; + /** * @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 view() const { - return SparsityView(colIdxs_.view(), rowIdxs_.view()); + return SparsityView(colIdxs_.view(), rowOffs_.view()); } private: diff --git a/include/NeoN/linearAlgebra/csrSparsityPattern.hpp b/include/NeoN/linearAlgebra/csrSparsityPattern.hpp index 9eb7134916..407fc29bcf 100644 --- a/include/NeoN/linearAlgebra/csrSparsityPattern.hpp +++ b/include/NeoN/linearAlgebra/csrSparsityPattern.hpp @@ -58,9 +58,16 @@ class CsrSparsityPattern : public NeoN::SupportsCopyTo; + /** * @brief Get a view representation of the matrix's data. * @return MatrixView for easy access to matrix elements. diff --git a/include/NeoN/linearAlgebra/ellSparsityPattern.hpp b/include/NeoN/linearAlgebra/ellSparsityPattern.hpp index cbb7a2c5ee..45a6b5e2f4 100644 --- a/include/NeoN/linearAlgebra/ellSparsityPattern.hpp +++ b/include/NeoN/linearAlgebra/ellSparsityPattern.hpp @@ -35,17 +35,29 @@ class EllSparsityPattern : public NeoN::SupportsCopyTo::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&& colIdx, Dimensions dim, localIdx numStoredElementsPerRow, - localIdx stride + localIdx stride, + localIdx logicalNnz ); [[nodiscard]] EllSparsityPattern copyToExecutor(Executor dstExec) const override { return EllSparsityPattern( - colIdxs_.copyToExecutor(dstExec), dimensions_, numStoredElementsPerRow_, stride_ + colIdxs_.copyToExecutor(dstExec), + dimensions_, + numStoredElementsPerRow_, + stride_, + logicalNnz_ ); } @@ -54,20 +66,22 @@ class EllSparsityPattern : public NeoN::SupportsCopyTo& colIdxs() const { return colIdxs_; }; - [[nodiscard]] Vector& colIdxs() { 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. + * count of logical/non-zero matrix entries -- see nnz(). */ - [[nodiscard]] localIdx nnz() const { return colIdxs_.size(); }; + [[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_; }; @@ -77,6 +91,8 @@ class EllSparsityPattern : public NeoN::SupportsCopyTo; + /** * @brief Get a view representation of the matrix's data. * @return EllSparsityView for easy access to matrix elements. @@ -100,6 +116,8 @@ class EllSparsityPattern : public NeoN::SupportsCopyTo values; //!< View to the values of the CSR matrix. + View 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 class Matrix : public NeoN::SupportsCopyTo> @@ -73,8 +74,9 @@ class Matrix : public NeoN::SupportsCopyTo> 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> } /** - * @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> const Vector& rowOffs, Dimensions dimensions ) + requires requires( + Vector c, + Vector r, + Dimensions d + ) { SparsityType(std::move(c), std::move(r), d); } : values_(values), sparsityPattern_( std::make_shared(Vector(colIdxs), Vector(rowOffs), dimensions) @@ -117,7 +125,8 @@ class Matrix : public NeoN::SupportsCopyTo> * @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> std::shared_ptr sparsity, std::shared_ptr faceToMatrixAddress ) - requires std::is_same_v + requires(std::is_same_v + && requires(const SparsityType& sp) { sp.rowOffs(); }) : values_(values), sparsityPattern_(sparsity), faceToMatrixAddress_(faceToMatrixAddress) { validate(); @@ -176,10 +186,12 @@ class Matrix : public NeoN::SupportsCopyTo> } /** - * @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& rowOffs() const + requires requires(const SparsityType& sp) { sp.rowOffs(); } { return sparsityPattern_->rowOffs(); } @@ -215,10 +227,9 @@ class Matrix : public NeoN::SupportsCopyTo> * @brief Get a view representation of the matrix's data. * @return MatrixView for easy access to matrix elements. */ - [[nodiscard]] MatrixView> - view() + [[nodiscard]] MatrixView view() { - return MatrixView>( + return MatrixView( values_.view(), sparsityPattern_->view() ); } @@ -227,12 +238,9 @@ class Matrix : public NeoN::SupportsCopyTo> * @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> - view() const + [[nodiscard]] MatrixView view() const { - return MatrixView>( + return MatrixView( View(values_.view()), sparsityPattern_->view() ); } @@ -245,7 +253,7 @@ class Matrix : public NeoN::SupportsCopyTo> private: - Vector values_; //!< The (non-zero) values of the CSR matrix. + Vector values_; //!< The (non-zero) values of the matrix. std::shared_ptr sparsityPattern_; @@ -259,6 +267,9 @@ using CSRMatrix = Matrix>; template using COOMatrix = Matrix>; +template +using ELLMatrix = Matrix>; + /** @brief extract the upper triangular of the matrix * @note this function is meant for testing purposes, it will recompute upper offsets */ diff --git a/include/NeoN/linearAlgebra/sparsityView.hpp b/include/NeoN/linearAlgebra/sparsityView.hpp index 5d2ab4bdbd..05967ab325 100644 --- a/include/NeoN/linearAlgebra/sparsityView.hpp +++ b/include/NeoN/linearAlgebra/sparsityView.hpp @@ -34,15 +34,21 @@ struct SparsityView SparsityView(View colIdxsView, View rowOffsView) : colIdxs(colIdxsView), rowOffs(rowOffsView) {}; + /** + * @brief Sentinel storage offset returned by findEntry() when (i,j) is not stored. + */ + KOKKOS_INLINE_FUNCTION + static constexpr IndexType invalidIndex() { return std::numeric_limits::max(); } /** - * @brief Retrieve a reference to the matrix element at position (i,j). + * @brief Retrieve the storage offset of the matrix element at position (i,j). * @param i The row index. * @param j The column index. - * @return Reference to the matrix element if it exists. + * @return Offset into colIdxs/values if it exists, invalidIndex() otherwise. + * @note assumes colIdxs is sorted ascending within each row's range; not verified. */ KOKKOS_INLINE_FUNCTION - IndexType entry(const IndexType i, const IndexType j) const + IndexType findEntry(const IndexType i, const IndexType j) const { const IndexType rowSize = rowOffs[i + 1] - rowOffs[i]; for (std::remove_const_t ic = 0; ic < rowSize; ++ic) @@ -54,8 +60,24 @@ struct SparsityView } if (colIdxs[localCol] > j) break; } - Kokkos::abort("Memory not allocated for CSR matrix component."); - return 0; // compiler warning suppression. + return invalidIndex(); + } + + /** + * @brief Retrieve a reference to the matrix element at position (i,j). + * @param i The row index. + * @param j The column index. + * @return Reference to the matrix element if it exists. + */ + KOKKOS_INLINE_FUNCTION + IndexType entry(const IndexType i, const IndexType j) const + { + const IndexType offset = findEntry(i, j); + if (offset == invalidIndex()) + { + Kokkos::abort("Memory not allocated for CSR matrix component."); + } + return offset; } View colIdxs; @@ -102,10 +124,12 @@ struct EllSparsityView * @brief Retrieve the storage offset of the matrix element at position (i,j). * @param i The row index. * @param j The column index. - * @return Flat offset into colIdxs/values if it exists. + * @return Flat offset into colIdxs/values if it exists, invalidIndex() otherwise. + * @note assumes each row's stored columns are sorted ascending with padding + * (invalidIndex()) trailing; not verified. */ KOKKOS_INLINE_FUNCTION - IndexType entry(const IndexType i, const IndexType j) const + IndexType findEntry(const IndexType i, const IndexType j) const { for (std::remove_const_t slot = 0; slot < numStoredElementsPerRow; ++slot) { @@ -114,8 +138,24 @@ struct EllSparsityView if (col == j) return idx; if (col == invalidIndex() || col > j) break; } - Kokkos::abort("Memory not allocated for ELL matrix component."); - return 0; // compiler warning suppression. + return invalidIndex(); + } + + /** + * @brief Retrieve the storage offset of the matrix element at position (i,j). + * @param i The row index. + * @param j The column index. + * @return Flat offset into colIdxs/values if it exists. + */ + KOKKOS_INLINE_FUNCTION + IndexType entry(const IndexType i, const IndexType j) const + { + const IndexType offset = findEntry(i, j); + if (offset == invalidIndex()) + { + Kokkos::abort("Memory not allocated for ELL matrix component."); + } + return offset; } View colIdxs; diff --git a/src/linearAlgebra/ellSparsityPattern.cpp b/src/linearAlgebra/ellSparsityPattern.cpp index 443fed2234..1e3eb8099a 100644 --- a/src/linearAlgebra/ellSparsityPattern.cpp +++ b/src/linearAlgebra/ellSparsityPattern.cpp @@ -18,14 +18,19 @@ void EllSparsityPattern::validate() const "ELL colIdxs size does not match stride * numStoredElementsPerRow" ); NF_ASSERT(stride_ >= dimensions_.rows, "ELL stride must be at least the number of rows"); + NF_ASSERT(logicalNnz_ <= colIdxs_.size(), "ELL logical nnz exceeds padded storage size"); } template EllSparsityPattern::EllSparsityPattern( - Vector&& colIdx, Dimensions dim, localIdx numStoredElementsPerRow, localIdx stride + Vector&& colIdx, + Dimensions dim, + localIdx numStoredElementsPerRow, + localIdx stride, + localIdx logicalNnz ) : dimensions_(dim), colIdxs_(std::move(colIdx)), - numStoredElementsPerRow_(numStoredElementsPerRow), stride_(stride) + numStoredElementsPerRow_(numStoredElementsPerRow), stride_(stride), logicalNnz_(logicalNnz) { validate(); } @@ -33,7 +38,8 @@ EllSparsityPattern::EllSparsityPattern( template EllSparsityPattern::EllSparsityPattern(const EllSparsityPattern& sp) : dimensions_(sp.dimensions_), colIdxs_(sp.colIdxs_), - numStoredElementsPerRow_(sp.numStoredElementsPerRow_), stride_(sp.stride_) + numStoredElementsPerRow_(sp.numStoredElementsPerRow_), stride_(sp.stride_), + logicalNnz_(sp.logicalNnz_) {} #define NN_DECLARE_SPARSITY(TYPENAME) template class EllSparsityPattern; diff --git a/src/linearAlgebra/matrix.cpp b/src/linearAlgebra/matrix.cpp index c17f2dbbc5..410e9eb3d1 100644 --- a/src/linearAlgebra/matrix.cpp +++ b/src/linearAlgebra/matrix.cpp @@ -15,20 +15,19 @@ Vector Matrix::diag() const { auto diag = Vector(values_.exec(), nRows()); fill(diag, zero()); - auto [diagV, rowOffsV, colIdxV, matrixV] = - views(diag, sparsityPattern_->rowOffs(), sparsityPattern_->colIdxs(), values_); + auto [diagV, matrixV] = views(diag, values_); + const auto sparsityV = sparsityPattern_->view(); + // Lenient: a missing diagonal keeps the zero fill() above instead of aborting. + // Traversal lives in *SparsityView::findEntry(), not here. parallelFor( values_.exec(), {0, nRows()}, NEON_LAMBDA(const localIdx rowi) { - for (auto i = rowOffsV[rowi]; i < rowOffsV[rowi + 1]; i++) + const auto offset = sparsityV.findEntry(rowi, rowi); + if (offset != decltype(sparsityV)::invalidIndex()) { - if (rowi == colIdxV[i]) - { - diagV[rowi] = matrixV[i]; - break; - } + diagV[rowi] = matrixV[offset]; } }, "copyDiag" @@ -45,14 +44,30 @@ Matrix Matrix::copyToExecutor( { return *this; } - return { - values_.copyToExecutor(dstExec), - std::make_shared(this->sparsityPattern_->copyToExecutor(dstExec)), - (faceToMatrixAddress_) ? std::make_shared( - this->faceToMatrixAddress_->copyToExecutor(dstExec) - ) - : nullptr - }; + auto copiedValues = values_.copyToExecutor(dstExec); + auto copiedSparsity = + std::make_shared(sparsityPattern_->copyToExecutor(dstExec)); + + if constexpr (std::is_same_v && requires(const SparsityType& sp) { + sp.rowOffs(); + }) + { + if (faceToMatrixAddress_) + { + return { + copiedValues, + copiedSparsity, + std::make_shared( + faceToMatrixAddress_->copyToExecutor(dstExec) + ) + }; + } + } + // faceToMatrixAddress_ can only be set via the constrained 3-arg constructor above, so + // it's null here whenever that constraint isn't met. Assert it so a future change can't + // silently drop it. + NF_ASSERT(!faceToMatrixAddress_, "Face address requires localIdx sparsity"); + return {copiedValues, copiedSparsity}; } @@ -349,4 +364,9 @@ void scaledInverseDiag( NN_DECLARE_MATRIX(scalar, localIdx); NN_DECLARE_MATRIX(Vec3, localIdx); +// ELL instantiated standalone, not via NN_DECLARE_MATRIX: upper() is CSR-shaped and has +// no ELL overload. +template class Matrix>; +template class Matrix>; + } diff --git a/test/linearAlgebra/cooSparsityPattern.cpp b/test/linearAlgebra/cooSparsityPattern.cpp index f10599c395..a797503ca9 100644 --- a/test/linearAlgebra/cooSparsityPattern.cpp +++ b/test/linearAlgebra/cooSparsityPattern.cpp @@ -49,6 +49,34 @@ TEST_CASE("SparsityPattern") { REQUIRE_THAT(bsp->rowIdxs(), Equals(std::vector {0, 3}, EqualInt())); } + + // Regression: view() used to expose rowIdxs_ where findEntry()/entry() need rowOffs_. + // Storage offsets 0..9, every row has a diagonal. + SECTION("view().findEntry() resolves offsets via rowOffs_, not rowIdxs_ " + execName) + { + Vector checkOffset(exec, 2); + auto checkOffsetView = checkOffset.view(); + auto cooView = sp->view(); + parallelFor( + exec, + {0, 1}, + NEON_LAMBDA(const localIdx) { + checkOffsetView[0] = cooView.findEntry(1, 1); // present -> storage offset 3 + checkOffsetView[1] = cooView.findEntry(0, 3); // absent -> invalidIndex() + } + ); + auto checkOffsetHost = checkOffset.copyToHost(); + auto checkOffsetHostView = checkOffsetHost.view(); + REQUIRE(checkOffsetHostView[0] == 3); + REQUIRE(checkOffsetHostView[1] == decltype(cooView)::invalidIndex()); + } + + SECTION("COOMatrix::diag() extracts the correct values " + execName) + { + NeoN::Vector values(exec, {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0}); + NeoN::la::COOMatrix cooMatrix(values, sp); + REQUIRE_THAT(cooMatrix.diag(), Equals(I({1.0, 4.0, 7.0, 10.0}))); + } } } diff --git a/test/linearAlgebra/ellSparsityPattern.cpp b/test/linearAlgebra/ellSparsityPattern.cpp index d34367b8a8..cebfc46e0e 100644 --- a/test/linearAlgebra/ellSparsityPattern.cpp +++ b/test/linearAlgebra/ellSparsityPattern.cpp @@ -6,6 +6,7 @@ // a custom main #include +#include #include "catch2_common.hpp" @@ -14,6 +15,17 @@ namespace NeoN { +// Freezes ELL's exclusion from the FaceToMatrixAddress constructor -- CSR row-local +// offsets don't apply to ELL's column-major storage. +static_assert( + !std::is_constructible_v< + la::ELLMatrix, + Vector, + std::shared_ptr>, + std::shared_ptr>, + "ELLMatrix must not be constructible with a FaceToMatrixAddress" +); + TEST_CASE("EllSparsityPattern") { using EllSparsityType = NeoN::la::EllSparsityPattern; @@ -39,6 +51,7 @@ TEST_CASE("EllSparsityPattern") const localIdx nRows = 4; const localIdx numStoredElementsPerRow = 3; const localIdx stride = nRows; + const localIdx logicalNnz = 10; // true nonzeros: row lengths 2+3+3+2, no padding auto colIdxExp = std::vector { 0, @@ -57,7 +70,11 @@ TEST_CASE("EllSparsityPattern") Vector colIdx(exec, colIdxExp); auto sp = std::make_shared( - std::move(colIdx), la::Dimensions {nRows, nRows}, numStoredElementsPerRow, stride + std::move(colIdx), + la::Dimensions {nRows, nRows}, + numStoredElementsPerRow, + stride, + logicalNnz ); SECTION("Can store padded colIdxs and dimensions " + execName) @@ -66,7 +83,9 @@ TEST_CASE("EllSparsityPattern") REQUIRE(sp->rows() == nRows); REQUIRE(sp->numStoredElementsPerRow() == numStoredElementsPerRow); REQUIRE(sp->stride() == stride); - REQUIRE(sp->nnz() == stride * numStoredElementsPerRow); + REQUIRE(sp->storageSize() == stride * numStoredElementsPerRow); + REQUIRE(sp->nnz() == logicalNnz); + REQUIRE(sp->nnz() != sp->storageSize()); REQUIRE(sp->dimension().rows == nRows); REQUIRE(sp->dimension().cols == nRows); } @@ -77,6 +96,8 @@ TEST_CASE("EllSparsityPattern") REQUIRE_THAT(spOnHost.colIdxs(), Equals(colIdxExp, EqualInt())); REQUIRE(spOnHost.numStoredElementsPerRow() == numStoredElementsPerRow); REQUIRE(spOnHost.stride() == stride); + REQUIRE(spOnHost.nnz() == logicalNnz); + REQUIRE(spOnHost.storageSize() == stride * numStoredElementsPerRow); } SECTION("Can resolve entry offsets on " + execName) @@ -101,6 +122,143 @@ TEST_CASE("EllSparsityPattern") REQUIRE(checkOffsetHostView[2] == 10); REQUIRE(checkOffsetHostView[3] == 7); } + + SECTION("findEntry returns invalidIndex() for missing entries " + execName) + { + Vector checkFound(exec, 2); + auto checkFoundView = checkFound.view(); + auto ellView = sp->view(); + parallelFor( + exec, + {0, 1}, + NEON_LAMBDA(const localIdx) { + // row 0 has no entry in column 3 (row 0's neighbours are {0,1}) + checkFoundView[0] = ellView.findEntry(0, 3); + // row 3 has no entry in column 0 (row 3's neighbours are {2,3}) + checkFoundView[1] = ellView.findEntry(3, 0); + } + ); + auto checkFoundHost = checkFound.copyToHost(); + auto checkFoundHostView = checkFoundHost.view(); + REQUIRE(checkFoundHostView[0] == decltype(ellView)::invalidIndex()); + REQUIRE(checkFoundHostView[1] == decltype(ellView)::invalidIndex()); + } +} + +TEST_CASE("ELLMatrix") +{ + using EllSparsityType = NeoN::la::EllSparsityPattern; + auto [execName, exec] = GENERATE(allAvailableExecutor()); + + // Same 1D 4-cell / 3-internal-face mesh as the EllSparsityPattern test above. + const localIdx nRows = 4; + const localIdx numStoredElementsPerRow = 3; + const localIdx stride = nRows; + const localIdx logicalNnz = 10; + const auto INV = std::numeric_limits::max(); + + Vector colIdx(exec, std::vector {0, 0, 1, 2, 1, 1, 2, 3, INV, 2, 3, INV}); + auto sp = std::make_shared( + std::move(colIdx), + la::Dimensions {nRows, nRows}, + numStoredElementsPerRow, + stride, + logicalNnz + ); + + // storageSize() == 12 (padded); diagonal entries sit at flat offsets 0, 5, 6, 7, + // so diag() reads 1, 6, 7, 8. + Vector values( + exec, std::vector {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0} + ); + + SECTION("Constructs when values.size() == storageSize(), not nnz() " + execName) + { + NeoN::la::ELLMatrix ellMatrix(values, sp); + REQUIRE(ellMatrix.values().size() == sp->storageSize()); + REQUIRE(ellMatrix.values().size() != sp->nnz()); + REQUIRE(ellMatrix.nNonZeros() == logicalNnz); + } + + SECTION("diag() extracts the correct values, generic across sparsity types " + execName) + { + NeoN::la::ELLMatrix ellMatrix(values, sp); + REQUIRE_THAT(ellMatrix.diag(), Equals(I({1.0, 6.0, 7.0, 8.0}))); + } + + SECTION("view().entry() reads and writes through the generic MatrixView " + execName) + { + NeoN::la::ELLMatrix ellMatrix(values, sp); + const NeoN::la::ELLMatrix ellMatrixConst(values, sp); + + Vector checkRead(exec, 4); + auto checkReadView = checkRead.view(); + auto constView = ellMatrixConst.view(); + parallelFor( + exec, + {0, 1}, + NEON_LAMBDA(const localIdx) { + checkReadView[0] = constView.entry(0, 0); // flat offset 0 -> 1.0 + checkReadView[1] = constView.entry(1, 1); // flat offset 5 -> 6.0 + checkReadView[2] = constView.entry(2, 2); // flat offset 6 -> 7.0 + checkReadView[3] = constView.entry(3, 3); // flat offset 7 -> 8.0 + } + ); + REQUIRE_THAT(checkRead, Equals(I({1.0, 6.0, 7.0, 8.0}))); + + auto mutableView = ellMatrix.view(); + parallelFor( + exec, {0, 1}, NEON_LAMBDA(const localIdx) { mutableView.entry(0, 0) = -1.0; } + ); + auto writtenHost = ellMatrix.values().copyToHost(); + REQUIRE(writtenHost.view()[0] == -1.0); + } + + SECTION("copyToExecutor() preserves values and sparsity metadata " + execName) + { + NeoN::la::ELLMatrix ellMatrix(values, sp); + auto hostMatrix = ellMatrix.copyToExecutor(SerialExecutor()); + + REQUIRE(hostMatrix.values().size() == sp->storageSize()); + REQUIRE(hostMatrix.nNonZeros() == logicalNnz); + REQUIRE(hostMatrix.sparsity()->stride() == stride); + REQUIRE(hostMatrix.sparsity()->numStoredElementsPerRow() == numStoredElementsPerRow); + REQUIRE_THAT(hostMatrix.diag(), Equals(I({1.0, 6.0, 7.0, 8.0}))); + + Vector checkValue(SerialExecutor(), 1); + auto checkValueView = checkValue.view(); + auto hostView = hostMatrix.view(); + parallelFor( + SerialExecutor(), + {0, 1}, + NEON_LAMBDA(const localIdx) { checkValueView[0] = hostView.entry(1, 2); } + ); + // (row 1, col 2) sits at flat offset 1 + stride*2 = 9 -> values[9] == 10.0 + REQUIRE(checkValueView[0] == 10.0); + } + + // Same sparse example as the CSR "Can extract diagonal" test in matrix.cpp: + // [ 1 . . ] + // [ . 5 6 ] + // [ . 8 . ] + // row 2 has no diagonal -- diag() should leave it at zero. + SECTION("diag() leaves a missing diagonal at zero " + execName) + { + const localIdx smallNRows = 3; + const localIdx smallWidth = 2; + const localIdx smallStride = smallNRows; + Vector smallColIdx(exec, std::vector {0, 1, 1, INV, 2, INV}); + auto smallSp = std::make_shared( + std::move(smallColIdx), + la::Dimensions {smallNRows, smallNRows}, + smallWidth, + smallStride, + 4 + ); + Vector smallValues(exec, std::vector {1.0, 5.0, 8.0, 0.0, 6.0, 0.0}); + NeoN::la::ELLMatrix smallMatrix(smallValues, smallSp); + REQUIRE_THAT(smallMatrix.diag(), Equals(I({1.0, 5.0, 0.0}))); + } } } From 9239482c6ce145fba60312c1710a98ab46e31750 Mon Sep 17 00:00:00 2001 From: Dheeraj Raghunathan <64003931+dheerajraghunathan@users.noreply.github.com> Date: Thu, 16 Jul 2026 15:09:25 +0200 Subject: [PATCH 3/3] fix boundary-cell mapping --- include/NeoN/linearAlgebra/linearSystem.hpp | 8 ++++-- test/linearAlgebra/linearSystem.cpp | 30 +++++++++++++++++++++ 2 files changed, 36 insertions(+), 2 deletions(-) diff --git a/include/NeoN/linearAlgebra/linearSystem.hpp b/include/NeoN/linearAlgebra/linearSystem.hpp index c775ec5212..3d9495032f 100644 --- a/include/NeoN/linearAlgebra/linearSystem.hpp +++ b/include/NeoN/linearAlgebra/linearSystem.hpp @@ -398,7 +398,8 @@ LinearSystem 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& lsIn ) + requires requires(const BoundaryMatrixType& bm) { bm.sparsity()->rowIdxs(); } { auto ls = la::LinearSystem(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]); }, diff --git a/test/linearAlgebra/linearSystem.cpp b/test/linearAlgebra/linearSystem.cpp index 83b216834f..2f22c5738e 100644 --- a/test/linearAlgebra/linearSystem.cpp +++ b/test/linearAlgebra/linearSystem.cpp @@ -115,6 +115,36 @@ TEMPLATE_TEST_CASE("LinearSystem", "[template]", NeoN::scalar) REQUIRE(linearSystem.getMeshIterator()->name() == "CellBased"); } + // Regression: each boundary face must land on its own owner cell, not another one -- + // CooSparsityPattern::view() briefly exposed row-range data here instead. + SECTION("removeBoundaryContributions applies each correction to its own owner cell " + execName) + { + auto nCells = 4; + auto mesh = create1DUniformMesh(exec, nCells); + auto ls = NeoN::la::createEmptyLinearSystem(mesh); + + REQUIRE(ls.boundaryMatrix().values().size() == 2); + Vector boundaryValues(exec, {10.0, 20.0}); + ls.boundaryMatrix().values() = boundaryValues; + + auto bRowIdxsHost = ls.boundaryMatrix().sparsity()->rowIdxs().copyToHost(); + auto owner0 = bRowIdxsHost.view()[0]; + auto owner1 = bRowIdxsHost.view()[1]; + REQUIRE(owner0 != owner1); // sanity: the two boundary faces belong to different cells + + auto lsNoBnd = NeoN::la::removeBoundaryContributions(ls); + auto diagHost = lsNoBnd.matrix().diag().copyToHost(); + auto diagView = diagHost.view(); + + for (localIdx i = 0; i < nCells; ++i) + { + scalar expected = 0.0; + if (i == owner0) expected += 10.0; + if (i == owner1) expected += 20.0; + REQUIRE(diagView[i] == expected); + } + } + SECTION("view read/write " + execName) {