From 96772651d8d2fc6ed1f2829c206a76eeb00387ff Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Wed, 14 Oct 2026 22:37:39 +0900 Subject: [PATCH] perf(hex8): packetize ordered gradient and geometric contractions --- src/fesa/elements/hex8_high_low.cpp | 34 ++- src/fesa/elements/hex8_pair_simd.h | 24 ++ src/fesa/elements/hex8_pair_simd_avx2.cpp | 22 ++ src/fesa/elements/hex8_pair_simd_dispatch.cpp | 21 ++ tests/unit/elements/hex8_high_low_test.cpp | 225 ++++++++++++++++++ 5 files changed, 323 insertions(+), 3 deletions(-) diff --git a/src/fesa/elements/hex8_high_low.cpp b/src/fesa/elements/hex8_high_low.cpp index ccb873f..ed1f9dc 100644 --- a/src/fesa/elements/hex8_high_low.cpp +++ b/src/fesa/elements/hex8_high_low.cpp @@ -246,11 +246,38 @@ Result Hex8Element::EvaluateHighLow( gradient_stress[a][j][l] = gradients[a][j] * s[l][j]; for (std::size_t a = 0; a < 8; ++a) { std::array dots{}; + std::array weighted_geo{}; + if (use_packets) + for (std::size_t first_b = 0; first_b < 8; first_b += 4) { + hex8_internal::GeometryPacketInput input{}; + const auto raw_pair = [](HighLow pair) { + return hex8_internal::PacketPair{pair.high, pair.low}; + }; + input.weight = raw_pair(HighLow{point.weight}); + for (std::size_t j = 0; j < 3; ++j) { + input.gradient_a[j] = raw_pair(gradients[a][j]); + for (std::size_t l = 0; l < 3; ++l) + input.gradient_stress[j][l] = raw_pair(gradient_stress[a][j][l]); + for (std::size_t lane = 0; lane < 4; ++lane) + input.gradient_b[j].values[lane] = + raw_pair(gradients[first_b + lane][j]); + } + hex8_internal::GeometryPacketOutput output; + // Same selector contract as material; unsupported/linear stays below. + hex8_internal::EvaluateGeometryPacketAvx2(input, &output); + for (std::size_t lane = 0; lane < 4; ++lane) { + dots[first_b + lane] = {output.gradient_dot.values[lane].high, + output.gradient_dot.values[lane].low}; + weighted_geo[first_b + lane] = + {output.weighted_geometric.values[lane].high, + output.weighted_geometric.values[lane].low}; + } + } for (std::size_t b = 0; b < 8; ++b) { - const auto dot = Dot(gradients[a], gradients[b]); + const auto dot = use_packets ? dots[b] : Dot(gradients[a], gradients[b]); dots[b] = dot; HighLow geo; - if (nonlinear) + if (nonlinear && !use_packets) for (std::size_t j = 0; j < 3; ++j) for (std::size_t l = 0; l < 3; ++l) geo += gradient_stress[a][j][l] * gradients[b][l]; @@ -265,7 +292,8 @@ Result Hex8Element::EvaluateHighLow( material[row][column] += mat * HighLow{point.weight}; } if (i == k) - geometric[row][column] += geo * HighLow{point.weight}; + geometric[row][column] += + use_packets ? weighted_geo[b] : geo * HighLow{point.weight}; } } if (use_packets) diff --git a/src/fesa/elements/hex8_pair_simd.h b/src/fesa/elements/hex8_pair_simd.h index 6a0fbbd..4dffa07 100644 --- a/src/fesa/elements/hex8_pair_simd.h +++ b/src/fesa/elements/hex8_pair_simd.h @@ -23,6 +23,30 @@ struct MaterialPacketInput { PairPacket accumulated; PacketPair weight; }; +/// @brief One owned trial point/a and four independent b gradients, raw only. +struct GeometryPacketInput { + PacketPair gradient_a[3]; + PacketPair gradient_stress[3][3]; + PairPacket gradient_b[3]; + PacketPair weight; +}; +/// @brief Independent raw contractions; accumulator entries remain caller-owned. +struct GeometryPacketOutput { + PairPacket gradient_dot; + PairPacket weighted_geometric; +}; +using GeometryPacketEvaluator = void (*)(const GeometryPacketInput&, + GeometryPacketOutput*); +/// @brief Baseline checked zero-start Dot3/ordered j,l Geo9 and weight products. +/// @note Lane order is 0..3. Invalid pairs remain invalid for original gates. +void EvaluateGeometryPacketScalar(const GeometryPacketInput& input, + GeometryPacketOutput* output); +/// @brief Evaluates independent b lanes without horizontal reduction. +/// @pre The material selector has admitted CPU/OS/current masked FP state. +/// @note Keeps original zero-start, product association and j,l Add order. +void EvaluateGeometryPacketAvx2(const GeometryPacketInput& input, + GeometryPacketOutput* output); + using MaterialPacketEvaluator = void (*)(const MaterialPacketInput&, PairPacket*); /// @brief Evaluates each lane's original checked material graph in order 0..3. diff --git a/src/fesa/elements/hex8_pair_simd_avx2.cpp b/src/fesa/elements/hex8_pair_simd_avx2.cpp index 46e105d..3c7fd40 100644 --- a/src/fesa/elements/hex8_pair_simd_avx2.cpp +++ b/src/fesa/elements/hex8_pair_simd_avx2.cpp @@ -111,6 +111,10 @@ VectorPair Load(const PairPacket& value) { _mm256_setr_pd(value.values[0].low, value.values[1].low, value.values[2].low, value.values[3].low)}; } +/// @brief Replicates one raw first product without changing either component. +VectorPair Splat(PacketPair value) { + return {_mm256_set1_pd(value.high), _mm256_set1_pd(value.low)}; +} /// @brief Stores raw high/low fields without rounding or lane reduction. PairPacket Store(VectorPair value) { double high[4]; @@ -251,6 +255,20 @@ void EvaluateMaterialPacketAvx2(const MaterialPacketInput& input, *output = Store(Add(Load(input.accumulated), Multiply(material, Load(weight)), true)); } +void EvaluateGeometryPacketAvx2(const GeometryPacketInput& input, + GeometryPacketOutput* output) { + VectorPair dot{_mm256_setzero_pd(), _mm256_setzero_pd()}; + for (std::size_t i = 0; i < 3; ++i) + dot = Add(dot, Multiply(Splat(input.gradient_a[i]), + Load(input.gradient_b[i])), true); + VectorPair geo{_mm256_setzero_pd(), _mm256_setzero_pd()}; + for (std::size_t j = 0; j < 3; ++j) + for (std::size_t l = 0; l < 3; ++l) + geo = Add(geo, Multiply(Splat(input.gradient_stress[j][l]), + Load(input.gradient_b[l])), true); + output->gradient_dot = Store(dot); + output->weighted_geometric = Store(Multiply(geo, Splat(input.weight))); +} } // namespace fesa::hex8_internal #else namespace fesa::hex8_internal { @@ -258,5 +276,9 @@ void EvaluateMaterialPacketAvx2(const MaterialPacketInput& input, PairPacket* output) { EvaluateMaterialPacketScalar(input, output); } +void EvaluateGeometryPacketAvx2(const GeometryPacketInput& input, + GeometryPacketOutput* output) { + EvaluateGeometryPacketScalar(input, output); +} } // namespace fesa::hex8_internal #endif diff --git a/src/fesa/elements/hex8_pair_simd_dispatch.cpp b/src/fesa/elements/hex8_pair_simd_dispatch.cpp index a4838d7..97e0118 100644 --- a/src/fesa/elements/hex8_pair_simd_dispatch.cpp +++ b/src/fesa/elements/hex8_pair_simd_dispatch.cpp @@ -42,6 +42,27 @@ void ScalarPacketTwoProduct(const double left[4], const double right[4], output->values[lane] = {value.high, value.low}; } } +FESA_HEX8_BASELINE_NOINLINE +void EvaluateGeometryPacketScalar(const GeometryPacketInput& input, + GeometryPacketOutput* output) { + const auto value = [](PacketPair pair) { + return math_internal::HighLow{pair.high, pair.low}; + }; + for (std::size_t lane = 0; lane < 4; ++lane) { + math_internal::HighLow dot; + for (std::size_t i = 0; i < 3; ++i) + dot = AddPair(dot, MultiplyPair(value(input.gradient_a[i]), + value(input.gradient_b[i].values[lane]))); + math_internal::HighLow geo; + for (std::size_t j = 0; j < 3; ++j) + for (std::size_t l = 0; l < 3; ++l) + geo = AddPair(geo, MultiplyPair(value(input.gradient_stress[j][l]), + value(input.gradient_b[l].values[lane]))); + const auto weighted = MultiplyPair(geo, value(input.weight)); + output->gradient_dot.values[lane] = {dot.high, dot.low}; + output->weighted_geometric.values[lane] = {weighted.high, weighted.low}; + } +} #undef FESA_HEX8_BASELINE_NOINLINE bool SupportsMaterialPacketInstructions(std::uint32_t leaf1_ecx, diff --git a/tests/unit/elements/hex8_high_low_test.cpp b/tests/unit/elements/hex8_high_low_test.cpp index ad01adf..104979d 100644 --- a/tests/unit/elements/hex8_high_low_test.cpp +++ b/tests/unit/elements/hex8_high_low_test.cpp @@ -449,6 +449,231 @@ TEST(Hex8HighLow, PacketTailScreenRetainsExactBoundaryAndFallback) { } #endif +#if defined(_M_X64) || defined(__SSE2__) +fesa::hex8_internal::GeometryPacketOutput FrozenGeometryPacket( + const fesa::hex8_internal::GeometryPacketInput& input) { + const auto pair = [](fesa::hex8_internal::PacketPair value) { + return fesa::high_low_original::HighLow{value.high, value.low}; + }; + fesa::hex8_internal::GeometryPacketOutput output{}; + for (std::size_t lane = 0; lane < 4; ++lane) { + fesa::high_low_original::HighLow dot; + for (std::size_t i = 0; i < 3; ++i) + dot = fesa::high_low_original::Add(dot, + fesa::high_low_original::Multiply(pair(input.gradient_a[i]), + pair(input.gradient_b[i].values[lane]))); + fesa::high_low_original::HighLow geo; + for (std::size_t j = 0; j < 3; ++j) + for (std::size_t l = 0; l < 3; ++l) + geo = fesa::high_low_original::Add(geo, + fesa::high_low_original::Multiply(pair(input.gradient_stress[j][l]), + pair(input.gradient_b[l].values[lane]))); + const auto weighted = fesa::high_low_original::Multiply(geo, pair(input.weight)); + output.gradient_dot.values[lane] = {dot.high, dot.low}; + output.weighted_geometric.values[lane] = {weighted.high, weighted.low}; + } + return output; +} + +fesa::hex8_internal::GeometryPacketInput GeometryDyadicFixture() { + fesa::hex8_internal::GeometryPacketInput input{}; + input.weight = {0.125, 0}; + for (std::size_t j = 0; j < 3; ++j) { + input.gradient_a[j] = {static_cast(j + 1) * 0.25, 0}; + for (std::size_t l = 0; l < 3; ++l) { + // Deliberately asymmetric first-products and distinct b/component mapping. + input.gradient_stress[j][l] = + {static_cast(1 + 3 * j + l) * 0.125, 0}; + for (std::size_t lane = 0; lane < 4; ++lane) + input.gradient_b[l].values[lane] = + {static_cast(1 + 4 * l + lane) * 0.0625, 0}; + } + } + return input; +} + +// Original observation repeats geo*weight for three independent accumulators. +// Selected observation reuses the weighted pair but retains each checked Add. +void ExpectGeometryPacketBitsAndFlags( + const fesa::hex8_internal::GeometryPacketInput& input, bool prior_status) { + std::fenv_t original; + ASSERT_EQ(std::fegetenv(&original), 0); + const auto original_csr = _mm_getcsr(); + struct Observation { + fesa::hex8_internal::GeometryPacketOutput output; + fesa::hex8_internal::PairPacket accumulators[3]; + int flags; + unsigned int status; + }; + const auto current_pair = [](fesa::hex8_internal::PacketPair value) { + return HighLow{value.high, value.low}; + }; + const auto observe = [&](bool original_order) { + EXPECT_EQ(std::fesetenv(&original), 0); + _mm_setcsr(original_csr & ~0x3fU); + EXPECT_EQ(std::feclearexcept(FE_ALL_EXCEPT), 0); + if (prior_status) { + EXPECT_EQ(std::feraiseexcept(FE_INEXACT), 0); + _mm_setcsr(_mm_getcsr() | 0x02U); + } + Observation value{}; + if (!original_order) + fesa::hex8_internal::EvaluateGeometryPacketAvx2(input, &value.output); + for (std::size_t lane = 0; lane < 4; ++lane) { + HighLow geo; + if (original_order) { + HighLow dot; + for (std::size_t i = 0; i < 3; ++i) + dot = fesa::hex8_internal::AddPair(dot, + fesa::hex8_internal::MultiplyPair(current_pair(input.gradient_a[i]), + current_pair(input.gradient_b[i].values[lane]))); + for (std::size_t j = 0; j < 3; ++j) + for (std::size_t l = 0; l < 3; ++l) + geo = fesa::hex8_internal::AddPair(geo, + fesa::hex8_internal::MultiplyPair( + current_pair(input.gradient_stress[j][l]), + current_pair(input.gradient_b[l].values[lane]))); + value.output.gradient_dot.values[lane] = {dot.high, dot.low}; + } + const HighLow accumulators[3]{{}, + {static_cast(lane + 1), std::ldexp(1.0, -54)}, {-1, -0.0}}; + for (std::size_t component = 0; component < 3; ++component) { + const auto weighted = original_order + ? fesa::hex8_internal::MultiplyPair(geo, current_pair(input.weight)) + : current_pair(value.output.weighted_geometric.values[lane]); + if (original_order) + value.output.weighted_geometric.values[lane] = + {weighted.high, weighted.low}; + const auto accumulated = fesa::hex8_internal::AddPair( + accumulators[component], weighted); + value.accumulators[component].values[lane] = + {accumulated.high, accumulated.low}; + } + } + value.flags = std::fetestexcept(FE_ALL_EXCEPT); + value.status = _mm_getcsr() & 0x3fU; + return value; + }; + const auto baseline = observe(true); + const auto selected = observe(false); + fesa::hex8_internal::GeometryPacketOutput scalar; + fesa::hex8_internal::EvaluateGeometryPacketScalar(input, &scalar); + const auto expected = FrozenGeometryPacket(input); + EXPECT_EQ(std::fesetenv(&original), 0); + _mm_setcsr(original_csr); + EXPECT_EQ(selected.flags, baseline.flags); + EXPECT_EQ(selected.status, baseline.status); + const auto same_bits = [](const fesa::hex8_internal::PairPacket& a, + const fesa::hex8_internal::PairPacket& b) { + for (std::size_t lane = 0; lane < 4; ++lane) { + EXPECT_EQ(PairBits(a.values[lane].high), PairBits(b.values[lane].high)); + EXPECT_EQ(PairBits(a.values[lane].low), PairBits(b.values[lane].low)); + } + }; + same_bits(selected.output.gradient_dot, expected.gradient_dot); + same_bits(selected.output.weighted_geometric, expected.weighted_geometric); + same_bits(scalar.gradient_dot, expected.gradient_dot); + same_bits(scalar.weighted_geometric, expected.weighted_geometric); + for (std::size_t component = 0; component < 3; ++component) + same_bits(selected.accumulators[component], baseline.accumulators[component]); +} + +TEST(Hex8HighLow, GeometryPacketPreservesAsymmetricLanesAndUniformShapes) { + if (fesa::hex8_internal::SelectMaterialPacketEvaluator() == + fesa::hex8_internal::EvaluateMaterialPacketScalar) + GTEST_SKIP() << "Geometry AVX2 masked environment is not selected."; + RecordProperty("private_geometry_packet_backend", "avx2"); + static_assert(std::is_trivial_v); + static_assert(std::is_trivial_v); + std::size_t checks = 0; + for (int shape = 0; shape < 4; ++shape) { + auto input = GeometryDyadicFixture(); + for (std::size_t j = 0; j < 3; ++j) { + if ((shape & 1) != 0) { + input.gradient_a[j].low = std::ldexp(input.gradient_a[j].high, -55); + for (std::size_t l = 0; l < 3; ++l) + input.gradient_stress[j][l].low = + std::ldexp(input.gradient_stress[j][l].high, -55); + } + if ((shape & 2) != 0) + for (std::size_t lane = 0; lane < 4; ++lane) + input.gradient_b[j].values[lane].low = + std::ldexp(input.gradient_b[j].values[lane].high, -55); + } + ExpectGeometryPacketBitsAndFlags(input, false); + ExpectGeometryPacketBitsAndFlags(input, true); + checks += 2; + } + RecordProperty("avx2_geometry_uniform_bit_flag_checks", static_cast(checks)); +} + +TEST(Hex8HighLow, GeometryPacketRetainsSpecialLanesAndExactTailBoundary) { + if (fesa::hex8_internal::SelectMaterialPacketEvaluator() == + fesa::hex8_internal::EvaluateMaterialPacketScalar) + GTEST_SKIP() << "Geometry AVX2 masked environment is not selected."; + using fesa::hex8_internal::PacketPair; + const double tiny = std::numeric_limits::denorm_min(); + const double maximum = std::numeric_limits::max(); + const std::vector special{{0, 0}, {-0.0, 0}, {1, -0.0}, {1, 1}, + {tiny, 0}, {maximum, maximum}, + {std::numeric_limits::signaling_NaN(), 0}, + {0, std::numeric_limits::signaling_NaN()}}; + std::size_t checks = 0; + for (const auto pair : special) + for (std::size_t lane = 0; lane < 4; ++lane) { + auto input = GeometryDyadicFixture(); + input.gradient_b[1].values[lane] = pair; + ExpectGeometryPacketBitsAndFlags(input, false); + ExpectGeometryPacketBitsAndFlags(input, true); + checks += 2; + } + const auto odd = [](int exponent) { + return std::nextafter(std::ldexp(1.0, exponent), + std::numeric_limits::infinity()); + }; + for (int exponent : {-470, -471}) + for (std::size_t lane = 0; lane < 4; ++lane) { + auto input = GeometryDyadicFixture(); + input.gradient_stress[0][0] = {odd(-500), 0}; + input.gradient_b[0].values[lane] = {odd(exponent), 0}; + ExpectGeometryPacketBitsAndFlags(input, false); + checks++; + } + RecordProperty("avx2_geometry_special_tail_checks", static_cast(checks)); +} + +TEST(Hex8HighLow, GeometryPacketPreservesOrderedOverflowAndPreparedInvalid) { + if (fesa::hex8_internal::SelectMaterialPacketEvaluator() == + fesa::hex8_internal::EvaluateMaterialPacketScalar) + GTEST_SKIP() << "Geometry AVX2 masked environment is not selected."; + const double maximum = std::numeric_limits::max(); + std::size_t checks = 0; + for (bool overflow_order : {false, true}) { + fesa::hex8_internal::GeometryPacketInput input{}; + input.weight = {0.5, 0}; + for (std::size_t l = 0; l < 3; ++l) + for (std::size_t lane = 0; lane < 4; ++lane) + input.gradient_b[l].values[lane] = {1, 0}; + input.gradient_stress[0][0] = {maximum, 0}; + input.gradient_stress[0][1] = {overflow_order ? maximum : -maximum, 0}; + input.gradient_stress[0][2] = {overflow_order ? -maximum : maximum, 0}; + ExpectGeometryPacketBitsAndFlags(input, false); + ExpectGeometryPacketBitsAndFlags(input, true); + checks += 2; + } + auto input = GeometryDyadicFixture(); + const auto invalid_first = fesa::high_low_original::Multiply( + {std::ldexp(1.0, -600), 0}, {std::ldexp(1.0, -500), 0}); + input.gradient_stress[0][0] = {invalid_first.high, invalid_first.low}; + for (std::size_t lane = 0; lane < 4; ++lane) + input.gradient_b[0].values[lane] = {std::ldexp(1.0, 600), 0}; + ExpectGeometryPacketBitsAndFlags(input, false); + ExpectGeometryPacketBitsAndFlags(input, true); + checks += 2; + RecordProperty("avx2_geometry_ordered_invalid_checks", static_cast(checks)); +} +#endif + TEST(Hex8HighLow, PacketBackendSelectionRecoversAfterEnvironmentRestoration) { const auto initial = fesa::hex8_internal::SelectMaterialPacketEvaluator(); const int rounding = std::fegetround();