From 7a6c096735f7b5086262baf0d2ebbb474201a3bc Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Wed, 21 Oct 2026 08:16:05 +0900 Subject: [PATCH] perf: reuse owned shell Jet Hessians in accumulation --- src/fesa/elements/nonlinear_mitc4_shell.cpp | 33 ++- src/fesa/math/jet_arithmetic.h | 5 +- .../nonlinear_mitc4_shell_jet_baseline.cpp | 203 ++++++++++++++++++ .../nonlinear_mitc4_shell_jet_baseline.md | 27 +++ 4 files changed, 263 insertions(+), 5 deletions(-) diff --git a/src/fesa/elements/nonlinear_mitc4_shell.cpp b/src/fesa/elements/nonlinear_mitc4_shell.cpp index 93aa2f2..4d8eb0e 100644 --- a/src/fesa/elements/nonlinear_mitc4_shell.cpp +++ b/src/fesa/elements/nonlinear_mitc4_shell.cpp @@ -386,6 +386,9 @@ class Hessian { /// @brief Identifies only the immutable zero Hessian, independently of seeds. bool IsImplicitZero() const { return !values_; } + /// @brief Allows a consumed temporary to reuse only its own materialized data. + bool HasUniqueStorage() const { return values_ && values_.use_count() == 1; } + /// @brief Gets fresh output storage without redundant zero initialization. /// @note Only fresh operator results use this pointer; all 576 entries must /// be written before any read or copy of the result. @@ -488,6 +491,30 @@ Jet operator+(const Jet& left, const Jet& right) { return result; } +/// @brief Adds to a consumed temporary without allocating another Hessian. +/// @note Shared/live aliases and unmasked exceptions retain the original path. +Jet operator+(Jet&& left, const Jet& right) { + if ((_mm_getcsr() & _MM_MASK_MASK) != _MM_MASK_MASK || &left == &right || + !left.second.HasUniqueStorage()) + return static_cast(left) + right; + left.value = left.value + right.value; + left.directional_first = left.directional_first + right.directional_first; + left.directional_second = left.directional_second + right.directional_second; + for (std::size_t index = 0U; index < kGlobalDofCount; ++index) { + left.first[index] = left.first[index] + right.first[index]; + left.mixed_first[index] = left.mixed_first[index] + right.mixed_first[index]; + left.mixed_second[index] = left.mixed_second[index] + right.mixed_second[index]; + } + left.constant_derivatives = left.constant_derivatives && + right.constant_derivatives && + std::isfinite(left.value); + const double* right_second = right.second.ReadData(); + double* left_second = left.second.Data(); + for (std::size_t index = 0U; index < kHessianSize; ++index) + left_second[index] = left_second[index] + right_second[index]; + return std::move(left); +} + /// @brief Subtracts two second-order trial-coordinate quantities. Jet operator-(const Jet& left, const Jet& right) { Jet result{}; @@ -1105,7 +1132,7 @@ Jet FrameComponent(const JetVector3& left, const JetTensor& tensor, Jet result{}; for (std::size_t row = 0U; row < 3U; ++row) { for (std::size_t column = 0U; column < 3U; ++column) { - result = result + left[row] * tensor[row][column] * right[column]; + result = std::move(result) + left[row] * tensor[row][column] * right[column]; } } return result; @@ -1183,9 +1210,9 @@ Jet PhysicalEnergy(const std::array& coordinates, shear * strain[4U]}; Jet density{}; for (std::size_t component = 0U; component < strain.size(); ++component) { - density = density + strain[component] * stress[component]; + density = std::move(density) + strain[component] * stress[component]; } - energy = energy + 0.5 * reference.reference_jacobian * density; + energy = std::move(energy) + 0.5 * reference.reference_jacobian * density; } } if (!IsFinite(energy)) { diff --git a/src/fesa/math/jet_arithmetic.h b/src/fesa/math/jet_arithmetic.h index 439e439..986e74f 100644 --- a/src/fesa/math/jet_arithmetic.h +++ b/src/fesa/math/jet_arithmetic.h @@ -3,6 +3,7 @@ #include #include +#include namespace fesa::math_internal { /// @brief Multiplies fixed Jet matrices in row/column/index order. @@ -17,7 +18,7 @@ std::array, 3> MultiplyJetMatrix( for (std::size_t column = 0; column < 3; ++column) { result[row][column] = constant(0.0); for (std::size_t index = 0; index < 3; ++index) - result[row][column] = result[row][column] + + result[row][column] = std::move(result[row][column]) + left[row][index] * right[index][column]; } return result; @@ -31,7 +32,7 @@ std::array MultiplyJetMatrixVector( std::array result{constant(0.0),constant(0.0),constant(0.0)}; for (std::size_t row = 0; row < 3; ++row) for (std::size_t column = 0; column < 3; ++column) - result[row] = result[row] + matrix[row][column] * vector[column]; + result[row] = std::move(result[row]) + matrix[row][column] * vector[column]; return result; } } // namespace fesa::math_internal diff --git a/tests/unit/elements/nonlinear_mitc4_shell_jet_baseline.cpp b/tests/unit/elements/nonlinear_mitc4_shell_jet_baseline.cpp index 656168a..a05264d 100644 --- a/tests/unit/elements/nonlinear_mitc4_shell_jet_baseline.cpp +++ b/tests/unit/elements/nonlinear_mitc4_shell_jet_baseline.cpp @@ -630,6 +630,176 @@ __declspec(noinline) double Consume(const Jet& value) { return result + value.constant_derivatives + value.second.IsImplicitZero(); } +// T3 evidence is separate from both immutable earlier snapshot formats. These +// calls also exercise the original overload when built against the T1 source. +std::string AccumulationSnapshots(bool reassociate) { + const std::array boundary{ + 0.0, -0.0, 1.0, -1.0, 0.1, -0.3, 1e16, + std::numeric_limits::denorm_min(), + std::numeric_limits::max(), + std::numeric_limits::infinity(), + std::numeric_limits::quiet_NaN(), + std::numeric_limits::signaling_NaN()}; + std::string out; + const std::array rounding{FE_TONEAREST, FE_DOWNWARD, FE_UPWARD, + FE_TOWARDZERO}; + for (const int mode : rounding) { + for (const unsigned flush : {0u, 0x8040u}) { + for (const unsigned sticky : {0u, 0x21u}) { + for (std::size_t kind = 0; kind < boundary.size(); ++kind) { + ResetFp(); + Jet left = Variable(boundary[kind], kind, -0.25); + Jet right = Variable(0.5, (kind + 1) % 24, 0.125); + if (kind % 3 == 0) left = Constant(boundary[kind]); + if (kind % 3 == 1) { + left.second.Data()[25] = boundary[kind]; + right.second.Data()[48] = -0.5; + } + left.mixed_first[3] = boundary[kind]; + left.mixed_second[4] = -0.0; + const Jet saved_left = left; + const Jet saved_right = right; + const auto before_left = Snapshot("input", left, 0, 0, 0); + const auto before_right = Snapshot("input", right, 0, 0, 0); + const auto environment = [&] { + ResetFp(); + std::fesetround(mode); + _mm_setcsr((_mm_getcsr() & ~0x8040u) | flush | sticky); + }; + const auto observe = [&](const char* name, auto expression) { + environment(); + const auto before = FpState(); + const Jet result = expression(); + const int exceptions = std::fetestexcept(FE_ALL_EXCEPT); + const unsigned status = _statusfp(), mxcsr = _mm_getcsr(); + out += std::to_string(mode) + ":" + std::to_string(flush) + ":" + + std::to_string(sticky) + ":" + std::to_string(kind) + " " + + before + " " + FpState() + " " + + Snapshot(name, result, exceptions, status, mxcsr); + }; + const JetVector3 a{left, right, Constant(-1)}; + const JetVector3 b{right, Constant(1), left}; + observe("dot", [&] { + if (reassociate) + return a[0] * b[0] + (a[1] * b[1] + a[2] * b[2]); + return Dot(a, b); + }); + for (std::size_t component = 0; component < 3; ++component) + observe("cross", [&] { return Cross(a, b)[component]; }); + const JetMatrix3 matrix{{a, b, a}}; + const JetMatrix3 other{{b, a, b}}; + // Every matrix/vector result is consumed over the whole fixture; + // boundary groups rotate the selected entry to keep evidence bounded. + observe("matrix", [&] { + return Multiply(matrix, other)[kind % 3][(kind / 3) % 3]; + }); + observe("matrix-vector", [&] { + return Multiply(matrix, b)[kind % 3]; + }); + observe("frame", [&] { return FrameComponent(a, matrix, b); }); + observe("fixed-frame", [&] { + return FrameComponent(ConstantVector(Vector3{0.5, -1, 0}), matrix, + ConstantVector(Vector3{-0.0, 0.25, 1})); + }); + observe("energy-accumulation", [&] { + Jet energy{}; + for (std::size_t i = 0; i < 3; ++i) + energy = std::move(energy) + a[i] * b[i]; + return energy; + }); + observe("shared-rvalue", [&] { + Jet value = left; + return std::move(value) + right; + }); + observe("same-object-rvalue", [&] { + Jet value = right; + return std::move(value) + value; + }); + observe("unique-rvalue", [&] { + Jet value = left * right; + return std::move(value) + right; + }); + observe("constant-product", [&] { return Constant(-0.5) * left; }); + observe("scalar-product", [&] { return -0.5 * left; }); + if (before_left != Snapshot("input", saved_left, 0, 0, 0) || + before_left != Snapshot("input", left, 0, 0, 0) || + before_right != Snapshot("input", saved_right, 0, 0, 0) || + before_right != Snapshot("input", right, 0, 0, 0)) + throw std::runtime_error("accumulation input alias modified"); + } + } + } + } + ResetFp(); + return out; +} + +std::string AccumulationAliasSnapshots(bool overwrite_input) { + std::string out; + for (int kind = 0; kind < 4; ++kind) { + Jet value = Variable(kind == 0 ? -0.0 : 0.5, 0, -0.25) * + Variable(0.75, 1, 0.5); + if (kind == 1) value = value - value; + if (kind == 2) value.second.Data()[0] = + std::numeric_limits::max(); + const auto before = Snapshot("live-input", value, 0, 0, 0); + out += Observe("unique-same-object", [&] { + if (overwrite_input) value.first[0] = 17; // Fixture-only negative control. + return std::move(value) + value; + }); + if (before != Snapshot("live-input", value, 0, 0, 0)) + throw std::runtime_error("same-object live operand modified"); + Jet survivor = value; + out += Observe("shared-move", [&] { + return std::move(value) + Constant(-1); + }); + if (before != Snapshot("live-input", survivor, 0, 0, 0)) + throw std::runtime_error("surviving shared operand modified"); + } + return out; +} + +void AccumulationBench() { + const Jet left = Variable(0.5, 0, 0.25) * Variable(0.75, 1, -0.5); + const Jet right = Variable(0.25, 2, 0.125) * Variable(-0.5, 3, 0.25); + const JetVector3 a{left, right, left}, b{right, left, right}; + const JetMatrix3 matrix{{a, b, a}}, other{{b, a, b}}; + const auto e1 = ConstantVector(Vector3{0.5, -1, 0.25}); + const auto e2 = ConstantVector(Vector3{-0.25, 0.5, 1}); + volatile double consumed = 0; + for (int kind = 0; kind < 6; ++kind) { + const auto start = std::chrono::steady_clock::now(); + for (int i = 0; i < 3000; ++i) { + if (kind == 0) consumed = consumed + Consume(Dot(a, b)); + if (kind == 1) { + const auto result = Cross(a, b); + for (const auto& value : result) consumed = consumed + Consume(value); + } + if (kind == 2) { + const auto result = Multiply(matrix, other); + for (const auto& row : result) + for (const auto& value : row) consumed = consumed + Consume(value); + } + if (kind == 3) { + const auto result = Multiply(matrix, b); + for (const auto& value : result) consumed = consumed + Consume(value); + } + if (kind == 4) + consumed = consumed + Consume(FrameComponent(e1, matrix, e2)); + if (kind == 5) { + Jet density{}; + for (std::size_t component = 0; component < 3; ++component) + density = std::move(density) + a[component] * b[component]; + consumed = consumed + Consume(0.5 * density); + } + } + std::cout << "accumulation-bench," << kind << ",3000," + << std::chrono::duration( + std::chrono::steady_clock::now() - start).count() + << ",allocation-included,setup-excluded," << consumed << '\n'; + } +} + void Bench() { volatile double consumed = 0; for (int kind = 0; kind < 8; ++kind) { @@ -671,6 +841,39 @@ int main(int argc, char** argv) { fesa::Bench(); return 0; } + if (argc == 2 && std::string(argv[1]) == "--accumulation-bench") { + fesa::AccumulationBench(); + return 0; + } + if (argc >= 3 && (std::string(argv[1]) == "--accumulation-write" || + std::string(argv[1]) == "--accumulation-compare")) { + const bool alias = argc > 3 && std::string(argv[3]) == "--unique-alias"; + auto actual = alias ? fesa::AccumulationAliasSnapshots( + argc > 4 && std::string(argv[4]) == + "--overwrite-input") + : fesa::AccumulationSnapshots( + argc > 3 && std::string(argv[3]) == + "--reassociate"); + if (std::string(argv[1]) == "--accumulation-write") { + if (std::ifstream(argv[2], std::ios::binary)) + throw std::runtime_error("accumulation snapshot already exists"); + std::ofstream output(argv[2], std::ios::binary); + output << actual; + output.close(); + if (!output) throw std::runtime_error("accumulation snapshot write failure"); + } else { + std::ifstream input(argv[2], std::ios::binary); + const std::string expected((std::istreambuf_iterator(input)), {}); + if (expected != actual) + throw std::runtime_error("accumulation snapshot mismatch"); + } + std::cout << argv[1] + << (alias ? " passed; 8 unique/same-object/shared alias cases" + : " passed; 192 boundary/environment groups, Dot/" + "Cross/matrix/vector/frame/accumulation/constants") + << ", 651 channels and FP state\n"; + return 0; + } if (argc >= 3 && (std::string(argv[1]) == "--product-write" || std::string(argv[1]) == "--product-compare")) { auto actual = fesa::ProductSnapshots( diff --git a/tests/unit/elements/nonlinear_mitc4_shell_jet_baseline.md b/tests/unit/elements/nonlinear_mitc4_shell_jet_baseline.md index 72b528c..0398788 100644 --- a/tests/unit/elements/nonlinear_mitc4_shell_jet_baseline.md +++ b/tests/unit/elements/nonlinear_mitc4_shell_jet_baseline.md @@ -65,3 +65,30 @@ T2의 무가드 단일 외적 후보는 unmasked 최초 INVALID/OVERFLOW 순서 전치 성분과 276개 상삼각 scratch를 재사용한다. Unmasked 환경과 외적 거부는 기존 두 호출, 576개 reverse 배열과 row-major 결합을 유지한다. Masked snapshot과 별도 원본/후보 SEH 관측을 함께 확인하며, 정상 경로의 임시 쓰기 감소를 전체 stack 예약 감소로 보고하지 않는다. + +T3 누적 검사는 기존 T1/T2 snapshot을 확장하지 않고 별도 파일을 쓴다. 원본 source와 +보존 archive의 출처를 확인한 실행물에서만 새 경로에 생성한다. 실제 CMake 경로는 +원본 source staging을 CMake로 configure/build한 실행물과 구분해 연결한다. + +```bat + --accumulation-write + --accumulation-compare + --accumulation-compare --reassociate + --accumulation-write --unique-alias + --accumulation-compare --unique-alias + --accumulation-compare --unique-alias --overwrite-input + --accumulation-bench +``` + +일반 그룹은 12개 경계 × 4 rounding × 2 FTZ/DAZ × 2 초기 sticky의 192개 조합으로 +Dot·Cross·3×3 matrix/vector·일반/고정 frame·energy 형태 누적·상수 Jet/scalar를 +검사한다. 별도 8개 alias 그룹은 유일한 materialized Hessian을 가진 동일 객체, +surviving copy와 공유 rvalue를 검사한다. `--reassociate`와 `--overwrite-input`는 +fixture 전용 음성 대조이며 exit 1이 기대값이다. 모든 비교는 651개 채널·provenance· +implicit Hessian·FP 상태를 포함한다. 상수 Jet와 scalar를 서로 같은 연산으로 가정하지 않는다. + +`--accumulation-bench`는 여섯 연산 종류마다 3,000회, 전체 결과 채널 checksum을 +소비한다. Fixture 준비는 제외하고 연산의 할당·해제와 checksum은 포함한다. 이 일반 +실행물의 시간과 별도 복사 source의 Hessian materialization 계측을 구분한다. +Unmasked 최초 trap은 별도 SEH/process 관측으로, 곱셈 단계와 실제 덧셈 단계의 +검증을 구분한다. 전체 모델 시간·stack 예약 감소를 뜻하지 않는다.