baram2584/FESADev
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

perf(shell): skip proven zero derivative Hessian work

Browse Source
This commit is contained in:
KOKO\Mimi
2026-10-14 20:14:18 +09:00
parent 8bfe97dbf8
commit c397a147b0
2 changed files with 91 additions and 22 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/fesa/elements/nonlinear_mitc4_shell.cpp
+75 -22
View File
@@ -417,6 +417,9 @@ class Hessian {
const double* ReadData() const { return Storage().data(); }
/// @brief Identifies only the immutable zero Hessian, independently of seeds.
bool IsImplicitZero() const { return !values_; }
/// @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.
@@ -461,6 +464,9 @@ struct Jet {
double directional_second{};
std::array<double, kGlobalDofCount> mixed_first{};
std::array<double, kGlobalDofCount> mixed_second{};
// Provenance for numerically zero derivatives; a null Hessian alone is not
// constant provenance because Variable retains a nonzero gradient seed.
bool constant_derivatives{false};
};
using JetVector3 = std::array<Jet, 3>;
@@ -471,6 +477,7 @@ using JetTensor = std::array<std::array<Jet, 3>, 3>;
Jet Constant(const double value) {
Jet result{};
result.value = value;
result.constant_derivatives = true;
return result;
}
@@ -478,6 +485,7 @@ Jet Constant(const double value) {
Jet Variable(const double value, const std::size_t index,
const double direction = 0.0) {
Jet result = Constant(value);
result.constant_derivatives = false;
result.first[index] = 1.0;
result.directional_first = direction;
return result;
@@ -497,6 +505,11 @@ Jet operator+(const Jet& left, const Jet& right) {
result.mixed_second[index] =
left.mixed_second[index] + right.mixed_second[index];
}
result.constant_derivatives = left.constant_derivatives &&
right.constant_derivatives &&
std::isfinite(result.value);
if (left.second.IsImplicitZero() && right.second.IsImplicitZero())
return result;
const double* left_second = left.second.ReadData();
const double* right_second = right.second.ReadData();
double* result_second = result.second.OverwriteData();
@@ -520,6 +533,11 @@ Jet operator-(const Jet& left, const Jet& right) {
result.mixed_second[index] =
left.mixed_second[index] - right.mixed_second[index];
}
result.constant_derivatives = left.constant_derivatives &&
right.constant_derivatives &&
std::isfinite(result.value);
if (left.second.IsImplicitZero() && right.second.IsImplicitZero())
return result;
const double* left_second = left.second.ReadData();
const double* right_second = right.second.ReadData();
double* result_second = result.second.OverwriteData();
@@ -554,9 +572,29 @@ Jet operator*(const Jet& left, const Jet& right) {
2.0 * left.directional_first * right.mixed_first[index] +
left.value * right.mixed_second[index];
}
const bool finite_values = std::isfinite(left.value) &&
std::isfinite(right.value) &&
std::isfinite(result.value);
const auto finite_gradient = [](const auto& gradient) {
return std::all_of(gradient.begin(), gradient.end(),
[](const double value) { return std::isfinite(value); });
};
const bool zero_outer = finite_values &&
((left.constant_derivatives && finite_gradient(right.first)) ||
(right.constant_derivatives && finite_gradient(left.first)));
result.constant_derivatives = finite_values && left.constant_derivatives &&
right.constant_derivatives;
if (zero_outer && left.second.IsImplicitZero() &&
right.second.IsImplicitZero()) return result;
const double* left_second = left.second.ReadData();
const double* right_second = right.second.ReadData();
double* result_second = result.second.OverwriteData();
if (zero_outer) {
for (std::size_t index = 0U; index < kHessianSize; ++index)
result_second[index] = left_second[index] * right.value + 0.0 + 0.0 +
left.value * right_second[index];
return result;
}
Outer(left.first, right.first, result_second);
std::array<double, kHessianSize> reverse;
Outer(right.first, left.first, reverse.data());
@@ -582,6 +620,10 @@ Jet operator*(const double scalar, const Jet& value) {
result.mixed_first[index] = scalar * value.mixed_first[index];
result.mixed_second[index] = scalar * value.mixed_second[index];
}
result.constant_derivatives = value.constant_derivatives &&
std::isfinite(scalar) &&
std::isfinite(result.value);
if (std::isfinite(scalar) && value.second.IsImplicitZero()) return result;
const double* value_second = value.second.ReadData();
double* result_second = result.second.OverwriteData();
for (std::size_t index = 0U; index < kHessianSize; ++index) {
@@ -614,6 +656,10 @@ Jet Compose(const Jet& input, const double value, const double first,
input.first[index] * input.directional_second) +
first * input.mixed_second[index];
}
result.constant_derivatives = input.constant_derivatives &&
std::isfinite(input.value) && std::isfinite(value) &&
std::isfinite(first) && std::isfinite(second) && std::isfinite(third);
if (result.constant_derivatives) return result;
std::array<double, kGlobalDofCount> scaled_first{};
for (std::size_t row = 0U; row < kGlobalDofCount; ++row)
scaled_first[row] = second * input.first[row];
@@ -673,32 +719,38 @@ Jet Atan2(const Jet& y, const Jet& x) {
const double derivative_yx =
(y.value * y.value - x.value * x.value) / radius_fourth;
const double derivative_xx = 2.0 * x.value * y.value / radius_fourth;
std::array<double, kGlobalDofCount> scaled_y{};
std::array<double, kGlobalDofCount> scaled_x{};
for (std::size_t i = 0U; i < kGlobalDofCount; ++i) {
scaled_y[i] = derivative_yy * y.first[i];
scaled_x[i] = derivative_xx * x.first[i];
}
std::array<std::array<double, kHessianSize>, 4U> products;
Outer(scaled_y, y.first, products[0U].data());
Outer(y.first, x.first, products[1U].data());
Outer(x.first, y.first, products[2U].data());
Outer(scaled_x, x.first, products[3U].data());
const bool zero_hessian = x.constant_derivatives && y.constant_derivatives &&
std::isfinite(result.value) && std::isfinite(derivative_y) &&
std::isfinite(derivative_x) && std::isfinite(derivative_yy) &&
std::isfinite(derivative_yx) && std::isfinite(derivative_xx);
for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
result.first[index] =
derivative_y * y.first[index] + derivative_x * x.first[index];
}
const double* y_second = y.second.ReadData();
const double* x_second = x.second.ReadData();
double* result_second = result.second.OverwriteData();
for (std::size_t row = 0U; row < kGlobalDofCount; ++row) {
for (std::size_t column = 0U; column < kGlobalDofCount; ++column) {
const std::size_t index = row * kGlobalDofCount + column;
result_second[index] =
derivative_y * y_second[index] + derivative_x * x_second[index] +
products[0U][index] +
derivative_yx * (products[1U][index] + products[2U][index]) +
products[3U][index];
if (!zero_hessian) {
std::array<double, kGlobalDofCount> scaled_y{};
std::array<double, kGlobalDofCount> scaled_x{};
for (std::size_t i = 0U; i < kGlobalDofCount; ++i) {
scaled_y[i] = derivative_yy * y.first[i];
scaled_x[i] = derivative_xx * x.first[i];
}
std::array<std::array<double, kHessianSize>, 4U> products;
Outer(scaled_y, y.first, products[0U].data());
Outer(y.first, x.first, products[1U].data());
Outer(x.first, y.first, products[2U].data());
Outer(scaled_x, x.first, products[3U].data());
const double* y_second = y.second.ReadData();
const double* x_second = x.second.ReadData();
double* result_second = result.second.OverwriteData();
for (std::size_t row = 0U; row < kGlobalDofCount; ++row) {
for (std::size_t column = 0U; column < kGlobalDofCount; ++column) {
const std::size_t index = row * kGlobalDofCount + column;
result_second[index] =
derivative_y * y_second[index] + derivative_x * x_second[index] +
products[0U][index] +
derivative_yx * (products[1U][index] + products[2U][index]) +
products[3U][index];
}
}
}
// Use a bounded ratio for the directional third derivative only. Retain
@@ -717,6 +769,7 @@ Jet Atan2(const Jet& y, const Jet& x) {
result.mixed_first[i] = sign * angle.mixed_first[i];
result.mixed_second[i] = sign * angle.mixed_second[i];
}
result.constant_derivatives = zero_hessian && angle.constant_derivatives;
return result;
}