perf(elements): reuse masked Jet product outer transpose

src/fesa/elements/nonlinear_mitc4_shell.cpp

@@ -1,5 +1,7 @@
 #include "fesa/elements/nonlinear_mitc4_shell.h"
 
+#include <xmmintrin.h>
+
 #include <algorithm>
 #include <array>
 #include <cmath>
@@ -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<double, kGlobalDofCount>& left,
+bool Outer(const std::array<double, kGlobalDofCount>& left,
            const std::array<double, kGlobalDofCount>& 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<double>::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<double, kGlobalDofCount * (kGlobalDofCount - 1U) / 2U> 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<double, kHessianSize> reverse;
   Outer(right.first, left.first, reverse.data());
   for (std::size_t row = 0U; row < kGlobalDofCount; ++row) {