baram2584/FESADev
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

perf(hex8): prepare guarded private four-column paired tangent packet

Browse Source
This commit is contained in:
KOKO\Mimi
2026-10-14 20:12:20 +09:00
parent b474d127d8
commit 9b8e53dcfd
6 changed files with 498 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/fesa/CMakeLists.txt
+10
View File
@@ -44,6 +44,8 @@ add_library(
elements/element_factory.cpp
elements/hex8_element.cpp
elements/hex8_high_low.cpp
elements/hex8_pair_simd_dispatch.cpp
elements/hex8_pair_simd_avx2.cpp
elements/hex8_geometry.cpp
elements/euler_beam_3d.cpp
elements/finite_rotation.cpp
@@ -122,6 +124,14 @@ target_compile_options(
$<$<CXX_COMPILER_ID:MSVC>:/W4 /WX /EHsc>
)
# Keep the CPU/OS dispatch on the baseline ISA; only its private kernel uses AVX2.
if(CMAKE_CXX_COMPILER_ID STREQUAL "MSVC" AND
CMAKE_CXX_COMPILER_ARCHITECTURE_ID STREQUAL "x64")
target_compile_definitions(fesa_solver PRIVATE FESA_HEX8_PAIR_AVX2=1)
set_source_files_properties(elements/hex8_pair_simd_avx2.cpp
PROPERTIES COMPILE_OPTIONS /arch:AVX2)
endif()
add_executable(fesa_cli app/main.cpp)
set_target_properties(fesa_cli PROPERTIES OUTPUT_NAME fesa)
target_link_libraries(fesa_cli PRIVATE fesa_solver)
src/fesa/elements/hex8_high_low.cpp
+39 -5
View File
@@ -9,6 +9,7 @@
#include "hex8_geometry.h"
#include "hex8_linear_tangent_cache.h"
#include "hex8_pair_arithmetic.h"
#include "hex8_pair_simd.h"
namespace fesa {
namespace {
@@ -163,6 +164,11 @@ Result<Hex8HighLowResponse> Hex8Element::EvaluateHighLow(
return Result<Hex8HighLowResponse>::Failure(
Failure("Hex8 high+low constitutive coefficients are unsupported."));
energy = HighLow{0.5} * HighLow{kappa_} * volume * mean_trace * mean_trace;
const auto packet_evaluator = nonlinear
? hex8_internal::SelectMaterialPacketEvaluator()
: hex8_internal::EvaluateMaterialPacketScalar;
const bool use_packets =
packet_evaluator != hex8_internal::EvaluateMaterialPacketScalar;
for (std::size_t g = 0; g < 8; ++g) {
const auto& point = points_[g];
const auto& gradients = point_gradients[g];
@@ -238,9 +244,11 @@ Result<Hex8HighLowResponse> Hex8Element::EvaluateHighLow(
for (std::size_t j = 0; j < 3; ++j)
for (std::size_t l = 0; l < 3; ++l)
gradient_stress[a][j][l] = gradients[a][j] * s[l][j];
for (std::size_t a = 0; a < 8; ++a)
for (std::size_t a = 0; a < 8; ++a) {
std::array<HighLow, 8> dots{};
for (std::size_t b = 0; b < 8; ++b) {
const auto dot = Dot(gradients[a], gradients[b]);
dots[b] = dot;
HighLow geo;
if (nonlinear)
for (std::size_t j = 0; j < 3; ++j)
@@ -250,14 +258,40 @@ Result<Hex8HighLowResponse> Hex8Element::EvaluateHighLow(
for (std::size_t k = 0; k < 3; ++k) {
const auto row = 3U * a + i;
const auto column = 3U * b + k;
const auto mat = deviatoric_fg[a][i] * fg[b][k] +
shear_fg[b][i] * fg[a][k] +
shear_fft[i][k] * dot;
material[row][column] += mat * HighLow{point.weight};
if (!use_packets) {
const auto mat = deviatoric_fg[a][i] * fg[b][k] +
shear_fg[b][i] * fg[a][k] +
shear_fft[i][k] * dot;
material[row][column] += mat * HighLow{point.weight};
}
if (i == k)
geometric[row][column] += geo * HighLow{point.weight};
}
}
if (use_packets)
for (std::size_t i = 0; i < 3; ++i)
for (std::size_t first = 0; first < 24; first += 4) {
hex8_internal::MaterialPacketInput input;
input.weight = HighLow{point.weight};
const auto row = 3U * a + i;
for (std::size_t lane = 0; lane < 4; ++lane) {
const auto column = first + lane;
const auto b = column / 3U;
const auto k = column % 3U;
input.deviatoric_left[lane] = deviatoric_fg[a][i];
input.deviatoric_right[lane] = fg[b][k];
input.shear_left[lane] = shear_fg[b][i];
input.shear_right[lane] = fg[a][k];
input.tensor_left[lane] = shear_fft[i][k];
input.gradient_dot[lane] = dots[b];
input.accumulated[lane] = material[row][column];
}
hex8_internal::PairPacket output;
packet_evaluator(input, &output);
for (std::size_t lane = 0; lane < 4; ++lane)
material[row][first + lane] = output[lane];
}
}
}
if (nonlinear || !paired_linear_tangent_->ready) {
std::array<Triple, 8> bulk_fg{};
src/fesa/elements/hex8_pair_simd.h
+28
View File
@@ -0,0 +1,28 @@
#ifndef FESA_SRC_FESA_ELEMENTS_HEX8_PAIR_SIMD_H_
#define FESA_SRC_FESA_ELEMENTS_HEX8_PAIR_SIMD_H_
#include <array>
#include "../math/high_low_arithmetic.h"
namespace fesa::hex8_internal {
using PairPacket = std::array<math_internal::HighLow, 4>;
struct MaterialPacketInput {
PairPacket deviatoric_left, deviatoric_right;
PairPacket shear_left, shear_right;
PairPacket tensor_left, gradient_dot;
PairPacket accumulated;
math_internal::HighLow weight;
};
using MaterialPacketEvaluator = void (*)(const MaterialPacketInput&, PairPacket*);
void EvaluateMaterialPacketScalar(const MaterialPacketInput& input,
PairPacket* output);
void EvaluateMaterialPacketAvx2(const MaterialPacketInput& input,
PairPacket* output);
/// @brief Selects a private backend from CPU/OS and this thread's FP state.
/// @note Never mutates FP state or caches a thread's environment decision.
MaterialPacketEvaluator SelectMaterialPacketEvaluator();
} // namespace fesa::hex8_internal
#endif // FESA_SRC_FESA_ELEMENTS_HEX8_PAIR_SIMD_H_
src/fesa/elements/hex8_pair_simd_avx2.cpp
+183
View File
@@ -0,0 +1,183 @@
#include "hex8_pair_simd.h"
#include "hex8_pair_arithmetic.h"
#if defined(FESA_HEX8_PAIR_AVX2) && defined(_MSC_VER) && defined(_M_X64)
#include <immintrin.h>
namespace fesa::hex8_internal {
namespace {
struct VectorPair { __m256d high, low; };
__m256d Select(__m256d mask, __m256d selected, __m256d other) {
return _mm256_blendv_pd(other, selected, mask);
}
__m256d Finite(__m256d value) {
const auto exponent = _mm256_set1_epi64x(0x7ff0000000000000LL);
const auto nonfinite = _mm256_cmpeq_epi64(
_mm256_and_si256(_mm256_castpd_si256(value), exponent), exponent);
return _mm256_castsi256_pd(_mm256_xor_si256(
nonfinite, _mm256_set1_epi64x(-1)));
}
__m256d Zero(__m256d value) {
return _mm256_castsi256_pd(_mm256_cmpeq_epi64(
_mm256_and_si256(_mm256_castpd_si256(value),
_mm256_set1_epi64x(0x7fffffffffffffffLL)),
_mm256_setzero_si256()));
}
VectorPair InvalidOutside(VectorPair value, __m256d valid) {
return {Select(valid, value.high, _mm256_set1_pd(
std::numeric_limits<double>::quiet_NaN())),
Select(valid, value.low, _mm256_setzero_pd())};
}
VectorPair Load(const PairPacket& value) {
return {_mm256_setr_pd(value[0].high, value[1].high,
value[2].high, value[3].high),
_mm256_setr_pd(value[0].low, value[1].low,
value[2].low, value[3].low)};
}
PairPacket Store(VectorPair value) {
double high[4], low[4];
_mm256_storeu_pd(high, value.high);
_mm256_storeu_pd(low, value.low);
return {{{high[0], low[0]}, {high[1], low[1]},
{high[2], low[2]}, {high[3], low[3]}}};
}
VectorPair TwoSum(__m256d left, __m256d right) {
const auto sum = _mm256_add_pd(left, right);
const auto bridge = _mm256_sub_pd(sum, left);
const auto restored = _mm256_sub_pd(sum, bridge);
const auto error = _mm256_add_pd(_mm256_sub_pd(left, restored),
_mm256_sub_pd(right, bridge));
const auto valid = _mm256_and_pd(_mm256_and_pd(Finite(sum), Finite(bridge)),
_mm256_and_pd(Finite(restored), Finite(error)));
const auto exact_zero = _mm256_and_pd(Zero(sum), Zero(error));
return InvalidOutside({Select(exact_zero, _mm256_setzero_pd(), sum),
Select(exact_zero, _mm256_setzero_pd(), error)}, valid);
}
VectorPair ScalarAdd(VectorPair left, VectorPair right, bool pair_shortcuts) {
const auto a = Store(left), b = Store(right);
PairPacket result;
for (std::size_t lane = 0; lane < result.size(); ++lane)
result[lane] = pair_shortcuts ? AddPair(a[lane], b[lane])
: math_internal::Add(a[lane], b[lane]);
return Load(result);
}
VectorPair ScalarMultiply(const PairPacket& left, const PairPacket& right) {
PairPacket result;
for (std::size_t lane = 0; lane < result.size(); ++lane)
result[lane] = MultiplyPair(left[lane], right[lane]);
return Load(result);
}
VectorPair Add(VectorPair left, VectorPair right, bool pair_shortcuts) {
const auto a = Store(left), b = Store(right);
int left_zero = 0, right_zero = 0;
for (std::size_t lane = 0; lane < a.size(); ++lane) {
if (!math_internal::IsNormalized(a[lane]) ||
!math_internal::IsNormalized(b[lane]))
return ScalarAdd(left, right, pair_shortcuts);
if (a[lane].high == 0 && a[lane].low == 0) left_zero |= 1 << lane;
if (b[lane].high == 0 && b[lane].low == 0) right_zero |= 1 << lane;
}
if (pair_shortcuts) {
if (left_zero == 15)
return {right.high, Select(Zero(right.low), _mm256_setzero_pd(), right.low)};
if (right_zero == 15)
return {left.high, Select(Zero(left.low), _mm256_setzero_pd(), left.low)};
if (left_zero != 0 || right_zero != 0)
return ScalarAdd(left, right, true);
}
const auto head = TwoSum(left.high, right.high);
const auto tail = TwoSum(left.low, right.low);
const auto middle = TwoSum(head.low, tail.high);
const auto joined = TwoSum(head.high, middle.high);
const auto remainder = _mm256_add_pd(
_mm256_add_pd(joined.low, middle.low), tail.low);
const auto valid = _mm256_and_pd(Finite(remainder), _mm256_and_pd(
_mm256_and_pd(Finite(head.high), Finite(head.low)), _mm256_and_pd(
_mm256_and_pd(Finite(tail.high), Finite(tail.low)), _mm256_and_pd(
_mm256_and_pd(Finite(middle.high), Finite(middle.low)),
_mm256_and_pd(Finite(joined.high), Finite(joined.low))))));
// Original Add returns before its last TwoSum when a prior finite gate fails.
return InvalidOutside(TwoSum(
Select(valid, joined.high, _mm256_setzero_pd()),
Select(valid, remainder, _mm256_setzero_pd())), valid);
}
VectorPair Product(__m256d left, __m256d right) {
double a[4], b[4];
_mm256_storeu_pd(a, left);
_mm256_storeu_pd(b, right);
bool special = false;
for (std::size_t lane = 0; lane < 4; ++lane)
if (!math_internal::IsFiniteBinary64(a[lane]) ||
!math_internal::IsFiniteBinary64(b[lane]) ||
a[lane] == 0 || b[lane] == 0 ||
math_internal::LowestBitExponent(a[lane]) +
math_internal::LowestBitExponent(b[lane]) < -1074)
special = true;
if (special) {
PairPacket result;
for (std::size_t lane = 0; lane < 4; ++lane)
result[lane] = math_internal::TwoProduct(a[lane], b[lane]);
return Load(result);
}
const auto product = _mm256_mul_pd(left, right);
const auto product_valid = Finite(product);
const auto one = _mm256_set1_pd(1.0);
// Invalid products never reach the residual FMA in the original helper.
const auto error = _mm256_fmsub_pd(Select(product_valid, left, one),
Select(product_valid, right, one), Select(product_valid, product, one));
const auto valid = _mm256_and_pd(product_valid, Finite(error));
return InvalidOutside(TwoSum(
Select(valid, product, _mm256_setzero_pd()),
Select(valid, error, _mm256_setzero_pd())), valid);
}
VectorPair Multiply(VectorPair left, VectorPair right) {
const auto a = Store(left), b = Store(right);
int branch = -1;
for (std::size_t lane = 0; lane < a.size(); ++lane) {
const bool normalized = math_internal::IsNormalized(a[lane]) &&
math_internal::IsNormalized(b[lane]);
if (!normalized) return ScalarMultiply(a, b);
const int shape = (a[lane].low != 0 ? 1 : 0) |
(b[lane].low != 0 ? 2 : 0);
if (branch != -1 && shape != branch) return ScalarMultiply(a, b);
branch = shape;
}
auto product = Product(left.high, right.high);
if (branch == 3) {
product = Add(product, Product(left.high, right.low), false);
product = Add(product, Product(left.low, right.high), false);
return Add(product, Product(left.low, right.low), false);
}
const VectorPair zero{_mm256_setzero_pd(), _mm256_setzero_pd()};
if ((branch & 2) != 0)
product = Add(product, Product(left.high, right.low), true);
else if ((branch & 1) != 0)
product = Add(product, zero, true);
if ((branch & 1) != 0)
product = Add(product, Product(left.low, right.high), true);
return Add(product, zero, true);
}
} // namespace
void EvaluateMaterialPacketAvx2(const MaterialPacketInput& input,
PairPacket* output) {
const auto first = Multiply(Load(input.deviatoric_left),
Load(input.deviatoric_right));
const auto second = Multiply(Load(input.shear_left), Load(input.shear_right));
const auto third = Multiply(Load(input.tensor_left), Load(input.gradient_dot));
const auto material = Add(Add(first, second, true), third, true);
const PairPacket weight{{input.weight, input.weight, input.weight, input.weight}};
*output = Store(Add(Load(input.accumulated),
Multiply(material, Load(weight)), true));
}
} // namespace fesa::hex8_internal
#else
namespace fesa::hex8_internal {
void EvaluateMaterialPacketAvx2(const MaterialPacketInput& input,
PairPacket* output) {
EvaluateMaterialPacketScalar(input, output);
}
} // namespace fesa::hex8_internal
#endif
src/fesa/elements/hex8_pair_simd_dispatch.cpp
+53
View File
@@ -0,0 +1,53 @@
#include "hex8_pair_simd.h"
#include "hex8_pair_arithmetic.h"
#if defined(FESA_HEX8_PAIR_AVX2) && defined(_MSC_VER) && defined(_M_X64)
#include <float.h>
#include <intrin.h>
#include <xmmintrin.h>
#endif
namespace fesa::hex8_internal {
void EvaluateMaterialPacketScalar(const MaterialPacketInput& input,
PairPacket* output) {
for (std::size_t lane = 0; lane < output->size(); ++lane) {
const auto first = MultiplyPair(input.deviatoric_left[lane],
input.deviatoric_right[lane]);
const auto second = MultiplyPair(input.shear_left[lane],
input.shear_right[lane]);
const auto third = MultiplyPair(input.tensor_left[lane],
input.gradient_dot[lane]);
const auto material = AddPair(AddPair(first, second), third);
(*output)[lane] = AddPair(input.accumulated[lane],
MultiplyPair(material, input.weight));
}
}
MaterialPacketEvaluator SelectMaterialPacketEvaluator() {
#if defined(FESA_HEX8_PAIR_AVX2) && defined(_MSC_VER) && defined(_M_X64)
static const bool hardware_supported = [] {
int info[4]{};
__cpuid(info, 0);
if (info[0] < 7) return false;
__cpuidex(info, 1, 0);
constexpr int required = (1 << 12) | (1 << 26) | (1 << 27) | (1 << 28);
if ((info[2] & required) != required) return false;
if ((_xgetbv(0) & 6U) != 6U) return false;
__cpuidex(info, 7, 0);
return (info[1] & (1 << 5)) != 0;
}();
if (!hardware_supported || !math_internal::ArithmeticEnvironmentSupported())
return EvaluateMaterialPacketScalar;
if ((_mm_getcsr() & 0x1f80U) != 0x1f80U)
return EvaluateMaterialPacketScalar;
unsigned int control = 0;
if (_controlfp_s(&control, 0, 0) != 0 ||
(control & _MCW_EM) != _MCW_EM)
return EvaluateMaterialPacketScalar;
return EvaluateMaterialPacketAvx2;
#else
return EvaluateMaterialPacketScalar;
#endif
}
} // namespace fesa::hex8_internal