perf(math): use checked MKL GER for rank-one products

tests/unit/math/dense_blas_internal_test.cpp

@@ -3,14 +3,65 @@
 #include <gtest/gtest.h>
 
 #include <array>
+#include <cfenv>
 #include <cmath>
 #include <cstddef>
+#include <cstdint>
+#include <cstring>
+#include <iostream>
 #include <limits>
 #include <stdexcept>
 
+#if defined(_M_X64) || defined(_M_IX86)
+#include <xmmintrin.h>
+#endif
+
 namespace fesa::dense_blas_internal {
 namespace {
 
+template <std::size_t N>
+void ExpectRankOneContract(const std::array<double, N>& left,
+                           const std::array<double, N>& right) {
+  const auto saved_left = left;
+  const auto saved_right = right;
+  const int rounding = std::fegetround();
+#if defined(_M_X64) || defined(_M_IX86)
+  const unsigned control = _mm_getcsr() & ~0x3fU;
+#endif
+  std::array<double, N * N> product{};
+  std::array<double, N * N> original{};
+  product.fill(std::numeric_limits<double>::quiet_NaN());
+  const ScopedLocalDenseThreads scope{};
+  ASSERT_TRUE(MultiplyTransposeRight(left.data(), N, right.data(), N, 1U,
+                                     original.data()));
+  ASSERT_TRUE(OuterProduct(left.data(), N, right.data(), N, product.data()));
+  std::size_t zero_sign_differences = 0U;
+  for (std::size_t i = 0U; i < product.size(); ++i) {
+    std::uint64_t original_bits = 0U;
+    std::uint64_t product_bits = 0U;
+    std::memcpy(&original_bits, &original[i], sizeof(original_bits));
+    std::memcpy(&product_bits, &product[i], sizeof(product_bits));
+    if (original[i] == 0.0 && product[i] == 0.0) {
+      zero_sign_differences += original_bits != product_bits ? 1U : 0U;
+      EXPECT_EQ(original[i], product[i]);
+    } else {
+      EXPECT_EQ(original_bits, product_bits)
+          << "entry=" << i << " left=" << left[i / N]
+          << " right=" << right[i % N]
+          << " original=" << original[i] << " candidate=" << product[i];
+    }
+  }
+  std::cout << "DGER boundary shape=" << N << "x" << N
+            << " zero_sign_differences=" << zero_sign_differences << '\n';
+  EXPECT_EQ(std::memcmp(saved_left.data(), left.data(), sizeof(left)), 0);
+  EXPECT_EQ(std::memcmp(saved_right.data(), right.data(), sizeof(right)), 0);
+  EXPECT_EQ(std::fegetround(), rounding);
+#if defined(_M_X64) || defined(_M_IX86)
+  // Arithmetic may raise sticky exception flags; control bits must be retained.
+  EXPECT_EQ(_mm_getcsr() & ~0x3fU, control);
+#endif
+}
+
 TEST(DenseBlasInternal,
      OuterProductUsesRectangularStorageAndPreservesRejection) {
   const std::array<double, 2U> left{2.0, -3.0};
@@ -38,6 +89,13 @@ TEST(DenseBlasInternal,
       std::numeric_limits<std::size_t>::max() / sizeof(double) + 1U;
   EXPECT_THROW(OuterProduct(nullptr, too_many_elements, nullptr, 1U, nullptr),
                std::length_error);
+  product = accepted;
+  EXPECT_THROW(OuterProduct(left.data(), too_large, right.data(), 1U,
+                            product.data()), std::length_error);
+  EXPECT_EQ(product, accepted);
+  EXPECT_THROW(OuterProduct(nullptr, 2U, right.data(), 3U, product.data()),
+               std::invalid_argument);
+  EXPECT_EQ(product, accepted);
 }
 
 TEST(DenseBlasInternal, LocalDenseScopeRestoresTlsAfterProductsAndUnwind) {
@@ -83,6 +141,60 @@ TEST(DenseBlasInternal, OuterProductLeavesFiniteInputOverflowForOwnerGate) {
   EXPECT_TRUE(std::isinf(product[0U]));
 }
 
+TEST(DenseBlasInternal, OuterProductMatchesK1ForSignedZeroAndExtremeScale) {
+  const std::array<double, 6U> left{
+      0.0, -0.0, std::numeric_limits<double>::denorm_min(),
+      std::ldexp(1.0, -500), std::numeric_limits<double>::max(),
+      -std::numeric_limits<double>::max()};
+  const std::array<double, 6U> right{
+      0.0, -0.0, 2.0, -2.0, std::ldexp(1.0, 500),
+      std::numeric_limits<double>::denorm_min()};
+  ExpectRankOneContract(left, right);
+}
+
+TEST(DenseBlasInternal, OuterProduct24PreservesNonzeroAndFailureBoundaries) {
+  std::array<double, 24U> left{};
+  std::array<double, 24U> right{};
+  for (std::size_t i = 0U; i < left.size(); ++i) {
+    left[i] = std::ldexp(static_cast<double>(i + 1U), -8);
+    right[i] = std::ldexp(-static_cast<double>(i + 2U), -10);
+  }
+  left[0U] = 0.0;
+  left[1U] = -0.0;
+  left[2U] = std::numeric_limits<double>::denorm_min();
+  left[3U] = -std::numeric_limits<double>::denorm_min();
+  left[4U] = std::ldexp(1.0, -500);
+  left[5U] = std::numeric_limits<double>::max();
+  right[0U] = -0.0;
+  right[1U] = 0.0;
+  right[2U] = std::numeric_limits<double>::denorm_min();
+  right[3U] = 2.0;
+  right[4U] = -2.0;
+  right[5U] = std::ldexp(1.0, 500);
+  ExpectRankOneContract(left, right);
+}
+
+TEST(DenseBlasInternal, OuterProductRejectsEitherNonfiniteInputBeforeWriting) {
+  std::array<double, 2U> left{2.0, -3.0};
+  std::array<double, 3U> right{4.0, -5.0, 0.25};
+  std::array<double, 6U> destination{1.0, -2.0, 3.0, -4.0, 5.0, -6.0};
+  const auto accepted = destination;
+  for (const double invalid : {std::numeric_limits<double>::quiet_NaN(),
+                               std::numeric_limits<double>::infinity(),
+                               -std::numeric_limits<double>::infinity()}) {
+    right[2U] = invalid;
+    EXPECT_FALSE(OuterProduct(left.data(), left.size(), right.data(),
+                              right.size(), destination.data()));
+    EXPECT_EQ(destination, accepted);
+    right[2U] = 0.25;
+    left[1U] = invalid;
+    EXPECT_FALSE(OuterProduct(left.data(), left.size(), right.data(),
+                              right.size(), destination.data()));
+    EXPECT_EQ(destination, accepted);
+    left[1U] = -3.0;
+  }
+}
+
 TEST(DenseBlasInternal, TransposeRightProductUsesInnerAndRectangularStrides) {
   const std::array<double, 6U> left{1, 2, 3, 4, 5, 6};
   const std::array<double, 12U> right{7,  8,  9,  10, 11, 12,