diff --git a/src/fesa/elements/nonlinear_mitc4_shell.cpp b/src/fesa/elements/nonlinear_mitc4_shell.cpp index e8d41e9..93aa2f2 100644 --- a/src/fesa/elements/nonlinear_mitc4_shell.cpp +++ b/src/fesa/elements/nonlinear_mitc4_shell.cpp @@ -1,5 +1,7 @@ #include "fesa/elements/nonlinear_mitc4_shell.h" +#include + #include #include #include @@ -412,11 +414,14 @@ class Hessian { /// @brief Forms an ordinary derivative outer product without a Matrix copy. /// @note Rejected nonfinite operands propagate to the existing Jet finite gate. -void Outer(const std::array& left, +bool Outer(const std::array& left, const std::array& right, double* result) { if (!dense_blas_internal::OuterProduct(left.data(), left.size(), right.data(), - right.size(), result)) + right.size(), result)) { std::fill_n(result, kHessianSize, std::numeric_limits::quiet_NaN()); + return false; + } + return true; } /// @brief Carries exact first and second trial-coordinate variations. @@ -559,7 +564,39 @@ Jet operator*(const Jet& left, const Jet& right) { left.value * right_second[index]; return result; } - Outer(left.first, right.first, result_second); + const bool masked_exceptions = + (_mm_getcsr() & _MM_MASK_MASK) == _MM_MASK_MASK; + const bool outer_accepted = Outer(left.first, right.first, result_second); + if (masked_exceptions && outer_accepted) { + // Retain earlier upper entries for later lower entries before overwriting + // them. Reuse the MKL k=1 transpose while preserving the original row-major + // arithmetic order as well as each entry's four-term association. + std::array upper; + std::size_t saved = 0U; + for (std::size_t row = 0U; row < kGlobalDofCount; ++row) { + for (std::size_t column = 0U; column < row; ++column) { + const std::size_t index = row * kGlobalDofCount + column; + result_second[index] = left_second[index] * right.value + + result_second[index] + upper[saved++] + + left.value * right_second[index]; + } + const std::size_t diagonal = row * kGlobalDofCount + row; + const double product = result_second[diagonal]; + result_second[diagonal] = left_second[diagonal] * right.value + + product + product + left.value * right_second[diagonal]; + for (std::size_t column = row + 1U; column < kGlobalDofCount; ++column) { + const std::size_t index = row * kGlobalDofCount + column; + const std::size_t transpose = column * kGlobalDofCount + row; + const double forward = result_second[index]; + const double reverse = result_second[transpose]; + upper[column * (column - 1U) / 2U + row] = forward; + result_second[index] = left_second[index] * right.value + + forward + reverse + left.value * right_second[index]; + } + } + return result; + } + // Preserve the reverse input-validation order after a rejected first outer. std::array reverse; Outer(right.first, left.first, reverse.data()); for (std::size_t row = 0U; row < kGlobalDofCount; ++row) { diff --git a/tests/unit/elements/nonlinear_mitc4_shell_jet_baseline.cpp b/tests/unit/elements/nonlinear_mitc4_shell_jet_baseline.cpp index 746c9b8..656168a 100644 --- a/tests/unit/elements/nonlinear_mitc4_shell_jet_baseline.cpp +++ b/tests/unit/elements/nonlinear_mitc4_shell_jet_baseline.cpp @@ -320,6 +320,154 @@ std::string Verify(bool corrupt_oracle) { return out; } +// Separate supplemental evidence: never changes the original T1 snapshot. +// A changed Hessian association, a lost transpose term, or a changed MKL zero +// convention must fail the bit/classification and complete FP-state comparison. +std::string ProductSnapshots(bool omit_transpose) { + const std::array boundary{ + 0.0, -0.0, 1.0, -1.0, 0.1, -0.3, + std::numeric_limits::denorm_min(), + -std::numeric_limits::denorm_min(), + std::numeric_limits::min(), + std::nextafter(std::numeric_limits::min(), 0.0), + std::numeric_limits::max(), + -std::numeric_limits::max(), + std::numeric_limits::infinity(), + -std::numeric_limits::infinity(), + std::numeric_limits::quiet_NaN(), + std::numeric_limits::signaling_NaN()}; + std::string out; + for (const int rounding : {FE_TONEAREST, FE_DOWNWARD, FE_UPWARD, + FE_TOWARDZERO}) { + for (const unsigned flush : {0u, 0x40u, 0x8000u, 0x8040u}) { + for (const bool sticky : {false, true}) { + auto reset = [&] { + ResetFp(); + std::fesetround(rounding); + _mm_setcsr((_mm_getcsr() & ~0x8040u) | flush); + if (sticky) std::feraiseexcept(FE_ALL_EXCEPT); + return FpState(); + }; + for (std::size_t kind = 0; kind < 96; ++kind) { + ResetFp(); // Construct inputs outside the observed environment. + Jet left{}, right{}; + left.value = 0.75; + right.value = -1.25; + left.directional_first = 0.25; + right.directional_first = -0.5; + left.directional_second = -0.125; + right.directional_second = 0.375; + std::uint64_t random = 0x719b245af51d6723ULL + kind; + auto coefficient = [&] { + random ^= random << 13; + random ^= random >> 7; + random ^= random << 17; + return (static_cast(random % 2000001) - 1000000) / + 1048576.0; + }; + for (std::size_t i = 0; i < 24; ++i) { + left.first[i] = kind < 32 ? boundary[(kind + i) % 16] + : coefficient(); + right.first[i] = kind < 32 ? boundary[(kind / 2 + 3 * i) % 16] + : coefficient(); + left.mixed_first[i] = coefficient(); + right.mixed_first[i] = coefficient(); + left.mixed_second[i] = coefficient(); + right.mixed_second[i] = coefficient(); + } + // Include dense, asymmetric inputs, implicit zeros and aliased input. + if (kind % 3 != 0) { + double* l = left.second.OverwriteData(); + double* r = right.second.OverwriteData(); + for (std::size_t i = 0; i < 576; ++i) { + l[i] = coefficient(); + r[i] = coefficient(); + } + } + if (kind == 32) left = Constant(-0.0); + if (kind == 33) right = Constant(0.0); + if (kind == 34) left.value = boundary[10]; + if (kind == 35) right.value = boundary[12]; + if (kind == 36) left.second.Data()[17] = boundary[10]; + if (kind == 37) right.second.Data()[24] = boundary[10]; + if (kind == 38) left.second.Data()[17] = boundary[14]; + if (kind == 39) right = left; // Shared immutable Hessian storage. + if (kind == 40) { + left.first[0] = boundary[12]; + right.first[0] = boundary[15]; + } + if (kind == 41) { + left.first[0] = boundary[15]; + right.first[0] = boundary[12]; + } + if (kind >= 42 && kind < 64) { + for (std::size_t i = 0; i < 24; ++i) { + left.first[i] = boundary[(kind + i) % 12]; + right.first[i] = boundary[(kind / 2 + 3 * i) % 12]; + if (kind == 42) { + left.first[i] = boundary[6]; + right.first[i] = i % 2 ? -0.5 : 0.5; + } + if (kind == 43) { + left.first[i] = boundary[10]; + right.first[i] = i % 2 ? -2.0 : 2.0; + } + if (kind == 44) { + left.first[i] = i % 2 ? -0.0 : 0.0; + right.first[i] = i % 3 ? -1.0 : 1.0; + } + } + } + const auto left_before = Snapshot("input", left, 0, 0, 0); + const auto right_before = Snapshot("input", right, 0, 0, 0); + const auto before = reset(); + Jet product = left * right; + if (omit_transpose && kind >= 64) { + // Deliberately omit the reverse outer term in one finite entry. + // This mutant tests arithmetic detection, not just file integrity. + product.second.Data()[1] = + left.second.ReadData()[1] * right.value + + left.first[0] * right.first[1] + + left.value * right.second.ReadData()[1]; + } + const int exceptions = std::fetestexcept(FE_ALL_EXCEPT); + const unsigned status = _statusfp(), mxcsr = _mm_getcsr(); + const auto after = FpState(); + out += "product " + std::to_string(kind) + " initial-fp=" + before + + " final-fp=" + after + " " + + Snapshot("result", product, exceptions, status, mxcsr); + if (left_before != Snapshot("input", left, 0, 0, 0) || + right_before != Snapshot("input", right, 0, 0, 0)) + throw std::runtime_error("supplemental product modified input"); + + // Independently pin the actual MKL operand-swap/transpose contract. + // Capture FP state before inspecting values (sNaN classification can + // itself raise INVALID). Each call starts from the same environment. + Jet forward{}, reverse{}; + reset(); + Outer(left.first, right.first, forward.second.OverwriteData()); + const auto forward_fp = FpState(); + reset(); + Outer(right.first, left.first, reverse.second.OverwriteData()); + const auto reverse_fp = FpState(); + for (std::size_t row = 0; row < 24; ++row) { + for (std::size_t column = 0; column < 24; ++column) { + const double a = forward.second.ReadData()[row * 24 + column]; + const double b = reverse.second.ReadData()[column * 24 + row]; + if (!(std::isnan(a) && std::isnan(b)) && Bits(a) != Bits(b)) + throw std::runtime_error("outer transpose coefficient mismatch"); + } + } + out += "forward-outer-fp=" + forward_fp + + " reverse-outer-fp=" + reverse_fp + "\n"; + } + } + } + } + ResetFp(); + return out; +} + // Public request-mode evidence for later direction-channel separation. Existing // tangent tests check nonzero-direction derivatives but do not freeze the full // omitted-direction response inventory and invalid-direction diagnostics. @@ -523,6 +671,29 @@ int main(int argc, char** argv) { fesa::Bench(); return 0; } + if (argc >= 3 && (std::string(argv[1]) == "--product-write" || + std::string(argv[1]) == "--product-compare")) { + auto actual = fesa::ProductSnapshots( + argc > 3 && std::string(argv[3]) == "--omit-transpose"); + if (argc > 3 && std::string(argv[3]) == "--mutate") + actual[actual.find('\n') - 1] ^= 1; + if (std::string(argv[1]) == "--product-write") { + if (std::ifstream(argv[2], std::ios::binary)) + throw std::runtime_error("supplemental snapshot already exists"); + std::ofstream output(argv[2], std::ios::binary); + output << actual; + output.close(); + if (!output) throw std::runtime_error("supplemental 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("supplemental product snapshot mismatch"); + } + std::cout << argv[1] << " passed; 3072 product cases and paired MKL outer " + "checks, 651 channels, input aliases, rounding/FTZ/DAZ/sticky\n"; + return 0; + } const bool verify_only = argc == 2 && std::string(argv[1]) == "--verify"; if (argc < 3 && !verify_only) throw std::runtime_error( diff --git a/tests/unit/elements/nonlinear_mitc4_shell_jet_baseline.md b/tests/unit/elements/nonlinear_mitc4_shell_jet_baseline.md index 1b8d834..72b528c 100644 --- a/tests/unit/elements/nonlinear_mitc4_shell_jet_baseline.md +++ b/tests/unit/elements/nonlinear_mitc4_shell_jet_baseline.md @@ -38,3 +38,30 @@ snapshot은 finite 값·signed zero·Inf 비트, NaN 분류, provenance와 impli 복원한 초기 FP 환경 및 연산 뒤 FE/rounding/control word/full MXCSR를 기록한다. 공개 반환의 방향 요청 유무·실패 진단·회복량 identity도 포함한다. 수식의 정확성은 독립 oracle로 검사하고, snapshot은 기준의 계산 순서와 관측 의미 보존을 검사한다. + +T2 곱셈 경계는 기존 T1 snapshot과 별도 CLI로 보존한다. 원본 제품 source의 hash를 +확인하고 보존한 원본 archive를 연결한 실행물에서만 새 경로에 한 번 생성한다. +후보는 `--product-compare`만 실행한다. 기존 파일을 덮어쓰거나 T1 기대값을 확장하지 않는다. + +```bat + --product-write + --product-compare + --product-compare --mutate + --product-compare --omit-transpose +``` + +마지막 두 명령은 각각 출력 변조와 reverse 외적 항 누락의 fixture 전용 음성 대조로 +exit 1이 기대값이다. 이 그룹은 96개 입력을 4 rounding × 4 FTZ/DAZ × 2 초기 sticky +환경에서 실행한 3072개 곱셈이다. 전체 651 채널, 비대칭·공유·implicit Hessian과 +입력 불변성, provenance, 전체 FP 상태를 비교한다. finite 전용 ±0·subnormal·underflow· +overflow gradient와 비유한 입력 거부를 구분하며, 실제 MKL 두 외적의 operand 교환과 +transpose 성분도 비교한다. 두 호출의 FP 상태는 각각 기록하므로 입력 검사 조기 종료로 +차이가 있는 경우를 수치 동등성으로 숨기지 않는다. NaN payload는 비교하지 않는다. +이 그룹은 예외를 mask한 상태의 관측이다. Unmasked 환경의 최초 trap 순서는 별도 +SEH/process 관측이나 생성 코드 근거가 필요하며, source의 row-major 순서만으로 보장하지 않는다. + +T2의 무가드 단일 외적 후보는 unmasked 최초 INVALID/OVERFLOW 순서가 달라 비채택했다. +채택한 guarded 경로는 MXCSR의 여섯 예외가 모두 masked이고 첫 외적이 성공한 경우에만 +전치 성분과 276개 상삼각 scratch를 재사용한다. Unmasked 환경과 외적 거부는 기존 두 호출, +576개 reverse 배열과 row-major 결합을 유지한다. Masked snapshot과 별도 원본/후보 SEH +관측을 함께 확인하며, 정상 경로의 임시 쓰기 감소를 전체 stack 예약 감소로 보고하지 않는다.