From 40a7e6c3afe39b004093291c2920d06ee1e0b30f Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Sat, 1 Aug 2026 23:06:21 +0900 Subject: [PATCH] =?UTF-8?q?feat(deterministic-parallel-assembly):=20step?= =?UTF-8?q?=200=20=E2=80=94=20canonical-contribution-order?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- CMakeLists.txt | 1 + include/fesa/assembly/contribution.hpp | 28 ++++ phases/index.json | 2 +- src/fesa/assembly/contribution.cpp | 132 ++++++++++++++++++ tests/CMakeLists.txt | 12 ++ tests/unit/assembly/serial_assembler_test.cpp | 78 ++++++++++- 6 files changed, 251 insertions(+), 2 deletions(-) create mode 100644 include/fesa/assembly/contribution.hpp create mode 100644 src/fesa/assembly/contribution.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 68681f6..4418637 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -28,6 +28,7 @@ include(cmake/FesaDependencies.cmake) add_library(fesa_core STATIC src/fesa/analysis/linear_static_analysis.cpp src/fesa/analysis/run_solver.cpp + src/fesa/assembly/contribution.cpp src/fesa/assembly/serial_assembler.cpp src/fesa/constraints/essential_bc.cpp src/fesa/core/version.cpp diff --git a/include/fesa/assembly/contribution.hpp b/include/fesa/assembly/contribution.hpp new file mode 100644 index 0000000..e7cf83e --- /dev/null +++ b/include/fesa/assembly/contribution.hpp @@ -0,0 +1,28 @@ +#pragma once + +#include +#include +#include +#include + +#include +#include + +namespace fesa { + +struct MatrixContribution final { + std::size_t row; + std::size_t column; + ElementId element; + std::uint16_t local_order; + double value; +}; + +[[nodiscard]] std::vector canonicalize_contributions( + std::span contributions); + +[[nodiscard]] SymmetricCsr merge_contributions( + std::size_t order, + std::span canonical); + +} // namespace fesa diff --git a/phases/index.json b/phases/index.json index ea538df..e28fa22 100644 --- a/phases/index.json +++ b/phases/index.json @@ -47,4 +47,4 @@ "status": "pending" } ] -} \ No newline at end of file +} diff --git a/src/fesa/assembly/contribution.cpp b/src/fesa/assembly/contribution.cpp new file mode 100644 index 0000000..54f7bef --- /dev/null +++ b/src/fesa/assembly/contribution.cpp @@ -0,0 +1,132 @@ +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace fesa { +namespace { + +std::int32_t csr_index(const std::size_t value) { + if (value > + static_cast( + std::numeric_limits::max())) { + throw std::overflow_error{ + "Symmetric CSR exceeds the 32-bit index range."}; + } + return static_cast(value); +} + +bool contribution_less( + const MatrixContribution& left, + const MatrixContribution& right) { + return std::tuple{ + left.row, + left.column, + left.element, + left.local_order} < + std::tuple{ + right.row, + right.column, + right.element, + right.local_order}; +} + +} // namespace + +std::vector canonicalize_contributions( + const std::span contributions) { + std::vector canonical{ + contributions.begin(), contributions.end()}; + std::ranges::sort(canonical, contribution_less); + return canonical; +} + +SymmetricCsr merge_contributions( + const std::size_t order, + const std::span canonical) { + if (!std::ranges::is_sorted(canonical, contribution_less)) { + throw std::invalid_argument{ + "Matrix contributions are not in canonical order."}; + } + + using Coordinate = std::pair; + std::vector coordinates; + coordinates.reserve(order + canonical.size()); + for (std::size_t row = 0; row < order; ++row) { + coordinates.emplace_back(row, row); + } + for (const MatrixContribution& contribution : canonical) { + if (contribution.row > contribution.column || + contribution.column >= order) { + throw std::invalid_argument{ + "Matrix contribution is outside the upper triangle."}; + } + coordinates.emplace_back( + contribution.row, contribution.column); + } + + std::ranges::sort(coordinates); + coordinates.erase( + std::ranges::unique(coordinates).begin(), + coordinates.end()); + + SymmetricCsr matrix{ + order, + std::vector(order + 1, 0), + {}, + {}, + }; + matrix.column_indices.reserve(coordinates.size()); + for (const auto [row, column] : coordinates) { + ++matrix.row_offsets[row + 1]; + matrix.column_indices.push_back(csr_index(column)); + } + for (std::size_t row = 0; row < order; ++row) { + const std::size_t offset = + static_cast(matrix.row_offsets[row]) + + static_cast(matrix.row_offsets[row + 1]); + matrix.row_offsets[row + 1] = csr_index(offset); + } + matrix.values.resize(matrix.column_indices.size(), 0.0); + + std::size_t contribution_index = 0; + while (contribution_index < canonical.size()) { + const MatrixContribution& first = + canonical[contribution_index]; + double value = 0.0; + do { + value += canonical[contribution_index].value; + ++contribution_index; + } while ( + contribution_index < canonical.size() && + canonical[contribution_index].row == first.row && + canonical[contribution_index].column == first.column); + + const auto row_begin = + matrix.column_indices.begin() + matrix.row_offsets[first.row]; + const auto row_end = + matrix.column_indices.begin() + + matrix.row_offsets[first.row + 1]; + const auto entry = std::lower_bound( + row_begin, + row_end, + csr_index(first.column)); + if (entry == row_end || *entry != csr_index(first.column)) { + throw std::logic_error{ + "Numeric contribution is absent from the CSR pattern."}; + } + matrix.values[static_cast( + std::distance(matrix.column_indices.begin(), entry))] = value; + } + + return matrix; +} + +} // namespace fesa diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index e35d231..7fc1dda 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -397,6 +397,18 @@ add_test( --gtest_filter=SymmetricCsr.* ) +add_test( + NAME CanonicalContribution + COMMAND "$" + --gtest_filter=CanonicalContribution.* +) + +add_test( + NAME DeterministicMerge + COMMAND "$" + --gtest_filter=DeterministicMerge.* +) + add_executable(fesa_constraint_tests unit/constraints/essential_bc_test.cpp ) diff --git a/tests/unit/assembly/serial_assembler_test.cpp b/tests/unit/assembly/serial_assembler_test.cpp index 257fb8b..cf80eb4 100644 --- a/tests/unit/assembly/serial_assembler_test.cpp +++ b/tests/unit/assembly/serial_assembler_test.cpp @@ -1,3 +1,4 @@ +#include #include #include @@ -187,8 +188,9 @@ fesa::Domain build_rounding_domain(const bool large_element_first) { storage_order = {2, 0, 1}; } for (const std::size_t index : storage_order) { + constexpr std::array element_ids{1, 2, 0}; builder.add_beam_element({ - fesa::ElementId{static_cast(index)}, + fesa::ElementId{element_ids[index]}, fesa::EntityOrigin{ "BeamPart", "Beam-1", @@ -250,6 +252,80 @@ double csr_value( std::distance(matrix.column_indices.begin(), found))]; } +std::vector value_bits( + const fesa::SymmetricCsr& matrix) { + std::vector bits; + bits.reserve(matrix.values.size()); + for (const double value : matrix.values) { + bits.push_back(std::bit_cast(value)); + } + return bits; +} + +TEST(CanonicalContribution, OrdersByCoordinateElementAndLocalOrder) { + const std::vector contributions{ + {1, 2, fesa::ElementId{9}, 4, 9.0}, + {0, 2, fesa::ElementId{3}, 8, 3.8}, + {0, 2, fesa::ElementId{3}, 2, 3.2}, + {0, 2, fesa::ElementId{1}, 7, 1.7}, + {0, 0, fesa::ElementId{5}, 0, 5.0}, + }; + + const std::vector canonical = + fesa::canonicalize_contributions(contributions); + + ASSERT_EQ(canonical.size(), contributions.size()); + EXPECT_EQ(canonical[0].row, 0); + EXPECT_EQ(canonical[0].column, 0); + EXPECT_EQ(canonical[0].element, fesa::ElementId{5}); + EXPECT_EQ(canonical[1].element, fesa::ElementId{1}); + EXPECT_EQ(canonical[2].element, fesa::ElementId{3}); + EXPECT_EQ(canonical[2].local_order, 2); + EXPECT_EQ(canonical[3].element, fesa::ElementId{3}); + EXPECT_EQ(canonical[3].local_order, 8); + EXPECT_EQ(canonical[4].row, 1); + EXPECT_EQ(canonical[4].column, 2); +} + +TEST(DeterministicMerge, IsBitwiseStableAcrossInputPermutations) { + const std::vector contributions{ + {0, 1, fesa::ElementId{2}, 0, 1.0}, + {0, 1, fesa::ElementId{0}, 0, 1.0e16}, + {0, 1, fesa::ElementId{1}, 0, -1.0e16}, + {1, 1, fesa::ElementId{1}, 1, -0.0}, + {1, 1, fesa::ElementId{0}, 1, +0.0}, + }; + const std::vector canonical = + fesa::canonicalize_contributions(contributions); + const fesa::SymmetricCsr expected = + fesa::merge_contributions(2, canonical); + + std::array permutation{0, 1, 2, 3, 4}; + do { + std::vector shuffled; + shuffled.reserve(contributions.size()); + for (const std::size_t index : permutation) { + shuffled.push_back(contributions[index]); + } + const std::vector shuffled_canonical = + fesa::canonicalize_contributions(shuffled); + const fesa::SymmetricCsr actual = + fesa::merge_contributions(2, shuffled_canonical); + + EXPECT_EQ(actual.order, expected.order); + EXPECT_EQ(actual.row_offsets, expected.row_offsets); + EXPECT_EQ(actual.column_indices, expected.column_indices); + EXPECT_EQ(value_bits(actual), value_bits(expected)); + } while (std::ranges::next_permutation(permutation).found); + + EXPECT_EQ( + std::bit_cast(csr_value(expected, 0, 1)), + std::bit_cast(1.0)); + EXPECT_EQ( + std::bit_cast(csr_value(expected, 1, 1)), + std::bit_cast(+0.0)); +} + TEST(SparsePattern, BuildsExpectedTwoElementChainStructure) { const fesa::Domain domain = build_chain_domain(false); const fesa::DofManager dofs = fesa::DofManager::build(domain);