fix(hex8): retain original contraction order in MKL entry products

src/fesa/elements/hex8_element.cpp

@@ -5,6 +5,7 @@
 #include <stdexcept>
 #include <utility>
 
+#include "../math/dense_blas_internal.h"
 #include "hex8_geometry.h"
 
 namespace fesa {
@@ -38,24 +39,39 @@ Matrix Outer(const std::array<Triple, 8>& values, double coefficient) {
     }
   return left.Multiply(right);
 }
-/// @brief Forms node-gradient Gram and stress congruence through dense MKL.
+/// @brief Uses MKL entry products and the original ordered Gram/stress sums.
 Matrix GradientProduct(const std::array<Triple, 8>& gradients,
                        const Tensor* stress = nullptr) {
-  Matrix left{8U, 3U};
-  Matrix transpose{3U, 8U};
+  const std::size_t terms = stress == nullptr ? 3U : 9U;
+  Vector left{64U * terms};
+  Vector right{64U * terms};
+  Vector final_factor{64U * terms};
   for (std::size_t a = 0; a < 8; ++a)
-    for (std::size_t i = 0; i < 3; ++i) {
-      left(a, i) = gradients[a][i];
-      transpose(i, a) = gradients[a][i];
-    }
-  if (stress != nullptr) {
-    Matrix stress_transpose{3U, 3U};
-    for (std::size_t i = 0; i < 3; ++i)
+    for (std::size_t b = 0; b < 8; ++b)
       for (std::size_t j = 0; j < 3; ++j)
-        stress_transpose(i, j) = (*stress)[j][i];
-    left = left.Multiply(stress_transpose);
-  }
-  return left.Multiply(transpose);
+        for (std::size_t l = 0; l < (stress == nullptr ? 1U : 3U); ++l) {
+          const auto entry = (8U * a + b) * terms +
+                             (stress == nullptr ? j : 3U * j + l);
+          left[entry] = gradients[a][j];
+          right[entry] = stress == nullptr ? gradients[b][j] : (*stress)[l][j];
+          final_factor[entry] = gradients[b][l];
+        }
+  auto products = dense_blas_internal::MultiplyDiagonal(left, right);
+  if (stress != nullptr)
+    products = dense_blas_internal::MultiplyDiagonal(products, final_factor);
+  Matrix result{8U, 8U};
+  for (std::size_t a = 0; a < 8; ++a)
+    for (std::size_t b = 0; b < 8; ++b) {
+      const auto first = (8U * a + b) * terms;
+      if (stress == nullptr) {
+        result(a, b) = (products[first] + products[first + 1U]) +
+                       products[first + 2U];
+      } else {
+        for (std::size_t term = 0; term < terms; ++term)
+          result(a, b) += products[first + term];
+      }
+    }
+  return result;
 }
 /// @brief Retains product and cancellation errors in fixed-order contractions.
 void AccumulateProduct(double left, double right, double& sum,