feat(analysis): execute solid generalized frequency pencil

src/fesa/analysis/solid_frequency_procedure.cpp

@@ -0,0 +1,351 @@
+#include "fesa/analysis/solid_frequency_procedure.h"
+
+#include <algorithm>
+#include <cmath>
+#include <limits>
+#include <memory>
+#include <new>
+#include <stdexcept>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include "fesa/analysis/analysis_model.h"
+#include "fesa/assembly/parallel_for.h"
+#include "fesa/assembly/sparse_assembler.h"
+#include "fesa/constraints/essential_constraint_policy.h"
+#include "fesa/elements/element_factory.h"
+#include "fesa/elements/hex8_element.h"
+#include "fesa/fem/dof_manager.h"
+#include "fesa/math/sparse_matrix.h"
+#include "fesa/math/vector.h"
+#include "fesa/model/domain.h"
+
+namespace fesa {
+namespace {
+/// @brief Reports a frequency failure without publishing partial state.
+Status Failure(const std::string& code, const std::string& message) {
+  return Status::Failure(
+      FailureCategory::kSolver,
+      {{Severity::kError, code, {}, "FREQUENCY", {}, message}});
+}
+/// @brief Rounds nonnegative arithmetic toward conservative bounds.
+double Up(double value) {
+  return std::nextafter(value, (std::numeric_limits<double>::infinity)());
+}
+/// @brief Rounds a lower bound downward, including underflow uncertainty.
+double Down(double value) {
+  return std::nextafter(value, -(std::numeric_limits<double>::infinity)());
+}
+/// @brief Uses robust Frobenius accumulation on original CSR storage.
+double Frobenius(const SparseMatrix& matrix) {
+  double norm = 0.0;
+  for (const double value : matrix.Values()) norm = std::hypot(norm, value);
+  return norm;
+}
+/// @brief Evaluates exact-zero denominator semantics for modal checks.
+double Ratio(double numerator, double denominator) {
+  if (denominator == 0.0)
+    return numerator == 0.0 ? 0.0 : (std::numeric_limits<double>::infinity)();
+  return numerator / denominator;
+}
+/// @brief Builds an immutable candidate; all returns precede public commit.
+Result<SolidFrequencyState> ExecuteCandidate(const Domain& domain) {
+  const auto* frequency = domain.Frequency();
+  if (!domain.IsSolid() || frequency == nullptr ||
+      frequency->requested_mode_count <= 0)
+    return Result<SolidFrequencyState>::Failure(Failure(
+        "invalid-solid-frequency-procedure", "Solid frequency step required."));
+  auto model = AnalysisModel::Create(domain);
+  if (!model.HasValue())
+    return Result<SolidFrequencyState>::Failure(model.GetStatus());
+  if (!model.Value().Step().Loads().empty() ||
+      !model.Value().ActiveSurfacePressureLoads().empty())
+    return Result<SolidFrequencyState>::Failure(
+        Failure("loaded-solid-frequency",
+                "Frequency requires unloaded reference state."));
+  std::vector<double> density(domain.Materials().Size(), 0.0);
+  for (const auto& item : domain.Densities()) {
+    if (item.material_index >= density.size() ||
+        density[item.material_index] != 0.0 || !std::isfinite(item.rho) ||
+        !(item.rho > 0.0))
+      return Result<SolidFrequencyState>::Failure(
+          Failure("invalid-material-density",
+                  "Density must be unique finite positive."));
+    density[item.material_index] = item.rho;
+  }
+  ElementFactory factory;
+  std::vector<std::unique_ptr<Element>> owned;
+  ElementView view;
+  for (const auto index : model.Value().ActiveElements()) {
+    const auto& definition = domain.Elements()[index];
+    if (definition.MaterialIndex() >= density.size() ||
+        !(density[definition.MaterialIndex()] > 0.0))
+      return Result<SolidFrequencyState>::Failure(
+          Failure("missing-frequency-density",
+                  "Every Solid material requires density."));
+    auto element = factory.Create(definition, domain);
+    if (!element.HasValue())
+      return Result<SolidFrequencyState>::Failure(element.GetStatus());
+    if (dynamic_cast<const Hex8Element*>(element.Value().get()) == nullptr)
+      return Result<SolidFrequencyState>::Failure(
+          Failure("unsupported-frequency-element", "Only C3D8 is supported."));
+    view.emplace_back(*element.Value());
+    owned.push_back(std::move(element.Value()));
+  }
+  DofManager dofs;
+  const Status mapping = dofs.Build(model.Value(), view);
+  if (!mapping.IsOk()) return Result<SolidFrequencyState>::Failure(mapping);
+  if (dofs.Profile() != DofProfile::kSolidStructural ||
+      dofs.ComponentsPerNode() != 3U)
+    return Result<SolidFrequencyState>::Failure(
+        Failure("invalid-solid-frequency-map",
+                "Frequency requires actual 3N mapping."));
+  for (std::size_t i = 0; i < dofs.PrescribedValues().Size(); ++i)
+    if (dofs.PrescribedValues()[i] != 0.0)
+      return Result<SolidFrequencyState>::Failure(Failure(
+          "unsupported-eigen-boundary", "Frequency requires homogeneous BC."));
+  const std::size_t requested =
+      static_cast<std::size_t>(frequency->requested_mode_count);
+  if (dofs.FreeDofCount() == 0U || requested > dofs.FreeDofCount())
+    return Result<SolidFrequencyState>::Failure(
+        Failure("insufficient-solid-frequency-modes",
+                "Requested modes exceed free dimension."));
+  const SerialParallelFor serial;
+  auto stiffness = SparseAssembler::Assemble(view, dofs, serial);
+  if (!stiffness.HasValue())
+    return Result<SolidFrequencyState>::Failure(stiffness.GetStatus());
+  std::vector<CooContribution> mass_coo;
+  std::vector<double> lower_by_dof(dofs.FullDofCount(), 0.0);
+  std::vector<double> magnitude_by_dof(dofs.FullDofCount(), 0.0);
+  std::vector<std::size_t> additions_by_dof(dofs.FullDofCount(), 0U);
+  constexpr double eps = std::numeric_limits<double>::epsilon();
+  for (std::size_t e = 0; e < owned.size(); ++e) {
+    const auto& hex = static_cast<const Hex8Element&>(*owned[e]);
+    auto mass = hex.ComputeMass();
+    if (!mass.HasValue())
+      return Result<SolidFrequencyState>::Failure(mass.GetStatus());
+    auto scatter = dofs.ElementScatter(hex.DofLayout());
+    if (!scatter.HasValue())
+      return Result<SolidFrequencyState>::Failure(scatter.GetStatus());
+    const auto weights = hex.ReferenceWeights();
+    const double min_weight = *std::min_element(weights.begin(), weights.end());
+    const auto element_index = model.Value().ActiveElements()[e];
+    const double rho =
+        density[domain.Elements()[element_index].MaterialIndex()];
+    const double local_lower = Down(Down(rho * min_weight) / 27.0);
+    // Relative gamma bounds require normal nonzero products. Fixed Gauss
+    // shapes exceed 1/128, so rho/16384 bounds the smallest density product.
+    if (!std::isfinite(local_lower) ||
+        local_lower < std::numeric_limits<double>::min() ||
+        rho / 16384.0 < std::numeric_limits<double>::min() ||
+        min_weight < std::numeric_limits<double>::min())
+      return Result<SolidFrequencyState>::Failure(
+          Failure("solid-mass-bound-unavailable",
+                  "Local reference mass bound is unrepresentable."));
+    // The eight fixed trilinear Gauss shapes have sigma_min(H)^2=1/27.
+    // 512 eps bounds their fixed evaluation, products and eight-term sums.
+    // Every mass term is positive; stable incidence scatter preserves the
+    // bound.
+    for (std::size_t i = 0; i < 24U; ++i) {
+      const auto row = scatter.Value()[i];
+      lower_by_dof[row] = Down(lower_by_dof[row] + local_lower);
+      double row_sum = 0.0;
+      for (std::size_t j = 0; j < 24U; ++j) {
+        const double value = mass.Value().values(i, j);
+        if (value != 0.0 &&
+            std::abs(value) < std::numeric_limits<double>::min())
+          return Result<SolidFrequencyState>::Failure(Failure(
+              "solid-mass-bound-unavailable",
+              "Subnormal local mass invalidates relative roundoff bounds."));
+        row_sum = Up(row_sum + std::abs(value));
+        mass_coo.push_back({row, scatter.Value()[j], value, e, i * 24U + j});
+      }
+      magnitude_by_dof[row] = Up(magnitude_by_dof[row] + row_sum);
+      ++additions_by_dof[row];
+    }
+  }
+  double ideal_lower = (std::numeric_limits<double>::infinity)();
+  double guard = 0.0;
+  for (const auto row : dofs.FreeDofs()) {
+    ideal_lower = (std::min)(ideal_lower, lower_by_dof[row]);
+    const double count_eps =
+        Up((512.0 + static_cast<double>(additions_by_dof[row])) * eps);
+    if (!(count_eps < 0.5))
+      return Result<SolidFrequencyState>::Failure(
+          Failure("solid-mass-bound-unavailable",
+                  "Mass reduction exceeds roundoff budget."));
+    const double gamma = Up(count_eps / (1.0 - count_eps));
+    // Convert the forward bound relative to the exact sum to an error bound
+    // relative to the available computed magnitude. Row infinity bounds also
+    // bound the symmetric perturbation's spectral norm under restriction.
+    const double computed_gamma = Up(gamma / (1.0 - gamma));
+    guard = (std::max)(guard, Up(computed_gamma * magnitude_by_dof[row]));
+  }
+  const double mass_lower = Down(ideal_lower - guard);
+  if (!std::isfinite(mass_lower) || !(mass_lower > 0.0) ||
+      !(guard < 0.01 * ideal_lower))
+    return Result<SolidFrequencyState>::Failure(
+        Failure("solid-mass-bound-unavailable",
+                "Roundoff prevents a useful SPD mass enclosure."));
+  auto mass =
+      SparseMatrix::FromCoo(dofs.FullDofCount(), dofs.FullDofCount(),
+                            std::move(mass_coo), dofs.GetSparsePattern());
+  if (!mass.HasValue())
+    return Result<SolidFrequencyState>::Failure(mass.GetStatus());
+  const EssentialConstraintPolicy partition;
+  auto k_blocks = partition.Partition(stiffness.Value(), dofs);
+  auto m_blocks = partition.Partition(mass.Value(), dofs);
+  if (!k_blocks.HasValue())
+    return Result<SolidFrequencyState>::Failure(k_blocks.GetStatus());
+  if (!m_blocks.HasValue())
+    return Result<SolidFrequencyState>::Failure(m_blocks.GetStatus());
+  const auto& k = k_blocks.Value().k_ff;
+  const auto& m = m_blocks.Value().k_ff;
+  double scale = 0.0;
+  for (std::size_t row = 0; row < k.Rows(); ++row) {
+    double kd = 0.0;
+    double md = 0.0;
+    for (std::size_t p = k.RowOffsets()[row]; p < k.RowOffsets()[row + 1U]; ++p)
+      if (k.ColumnIndices()[p] == row) kd = k.Values()[p];
+    for (std::size_t p = m.RowOffsets()[row]; p < m.RowOffsets()[row + 1U]; ++p)
+      if (m.ColumnIndices()[p] == row) md = m.Values()[p];
+    if (!(md > 0.0) || !std::isfinite(std::abs(kd) / md))
+      return Result<SolidFrequencyState>::Failure(
+          Failure("invalid-solid-eigen-scale",
+                  "Mass diagonal or classification scale is invalid."));
+    scale = (std::max)(scale, std::abs(kd) / md);
+  }
+  if (!(scale > 0.0) || !std::isfinite(scale))
+    return Result<SolidFrequencyState>::Failure(
+        Failure("invalid-solid-eigen-scale",
+                "A positive finite classification scale is required."));
+  auto backend = MklFeastEigenSolver{}.SolveGeneralized(k, m, requested,
+                                                        mass_lower, scale);
+  if (!backend.HasValue())
+    return Result<SolidFrequencyState>::Failure(backend.GetStatus());
+  SolidFrequencyState state(std::move(backend.Value()));
+  state.full_dof_count = dofs.FullDofCount();
+  state.free_dof_count = dofs.FreeDofCount();
+  state.classification_scale = scale;
+  state.mass_ideal_lower_bound = ideal_lower;
+  state.mass_roundoff_guard = guard;
+  state.mass_lower_bound = mass_lower;
+  state.constraint_mask.assign(dofs.FullDofCount(), 0U);
+  state.prescribed_displacement.assign(dofs.FullDofCount(), 0.0);
+  for (const auto row : dofs.ConstrainedDofs()) state.constraint_mask[row] = 1U;
+  const std::size_t count = state.backend.eigenvalues.size();
+  if (count < requested || state.backend.eigenvectors.Rows() != k.Rows() ||
+      state.backend.eigenvectors.Columns() != count)
+    return Result<SolidFrequencyState>::Failure(
+        Failure("insufficient-solid-frequency-pairs",
+                "Backend returned insufficient complete pairs."));
+  const double norm_k = Frobenius(k);
+  const double norm_m = Frobenius(m);
+  std::vector<Vector> phi;
+  std::vector<Vector> mass_phi;
+  phi.reserve(count);
+  mass_phi.reserve(count);
+  for (std::size_t c = 0; c < count; ++c) {
+    Vector vector(k.Rows());
+    for (std::size_t row = 0; row < k.Rows(); ++row)
+      vector[row] = state.backend.eigenvectors(row, c);
+    Vector mv = m.Multiply(vector);
+    const double raw_mass = vector.Dot(mv);
+    if (!(raw_mass > 0.0) || !std::isfinite(raw_mass))
+      return Result<SolidFrequencyState>::Failure(
+          Failure("solid-eigenpair-validation-failed",
+                  "Modal mass is not finite positive."));
+    vector.Scale(1.0 / std::sqrt(raw_mass));
+    mv = m.Multiply(vector);
+    const Vector kv = k.Multiply(vector);
+    const double lambda = state.backend.eigenvalues[c];
+    const double modal_mass = vector.Dot(mv);
+    const double energy = vector.Dot(kv);
+    Vector residual = kv;
+    residual.Axpy(-lambda, mv);
+    const double norm_phi = vector.Norm();
+    const double denominator = (norm_k + std::abs(lambda) * norm_m) * norm_phi;
+    const double eta_r = Ratio(residual.Norm(), denominator);
+    const double eta_norm = std::abs(modal_mass - 1.0);
+    const double eta_energy =
+        Ratio(std::abs(energy - lambda * modal_mass), denominator * norm_phi);
+    if (!std::isfinite(lambda) || !std::isfinite(eta_r) || eta_r > 1e-10 ||
+        !std::isfinite(eta_norm) || eta_norm > 1e-10 ||
+        !std::isfinite(eta_energy) || eta_energy > 1e-10 ||
+        !std::isfinite(energy / modal_mass) || lambda / scale < -1e-9)
+      return Result<SolidFrequencyState>::Failure(
+          Failure("solid-eigenpair-validation-failed",
+                  "Original K/M pair checks failed."));
+    state.eigenpair_checks.push_back({c + 1U, state.backend.column_ids[c],
+                                      lambda, modal_mass, energy / modal_mass,
+                                      eta_r, eta_norm, eta_energy});
+    phi.push_back(std::move(vector));
+    mass_phi.push_back(std::move(mv));
+  }
+  for (std::size_t a = 0; a < count; ++a)
+    for (std::size_t b = 0; b < count; ++b) {
+      const double dot = phi[a].Dot(mass_phi[b]);
+      state.eta_gram = std::hypot(state.eta_gram, dot - (a == b ? 1.0 : 0.0));
+      if (a != b)
+        state.max_eta_orth =
+            (std::max)(state.max_eta_orth,
+                       std::abs(dot) /
+                           std::sqrt(state.eigenpair_checks[a].modal_mass *
+                                     state.eigenpair_checks[b].modal_mass));
+    }
+  if (!std::isfinite(state.eta_gram) || state.eta_gram > 1e-10 ||
+      !std::isfinite(state.max_eta_orth) || state.max_eta_orth > 1e-10)
+    return Result<SolidFrequencyState>::Failure(
+        Failure("solid-eigenpair-validation-failed",
+                "All-pair mass Gram checks failed."));
+  state.vectors = Matrix(requested, dofs.FullDofCount());
+  std::size_t cluster = 1U;
+  for (std::size_t c = 0; c < requested; ++c) {
+    const double lambda = state.backend.eigenvalues[c];
+    if (c > 0 && (lambda - state.backend.eigenvalues[c - 1U]) / scale > 1e-8)
+      ++cluster;
+    std::optional<double> frequency_value;
+    if (lambda >= 0.0)
+      frequency_value = std::sqrt(lambda) / (2.0 * std::acos(-1.0));
+    state.modes.push_back(
+        {c + 1U, lambda,
+         static_cast<std::uint8_t>(std::abs(lambda / scale) <= 1e-9 ? 0U : 1U),
+         cluster, frequency_value});
+    // Strict > retains the first source/component entry for equal magnitudes.
+    std::size_t pivot = 0U;
+    for (std::size_t row = 1; row < phi[c].Size(); ++row)
+      if (std::abs(phi[c][row]) > std::abs(phi[c][pivot])) pivot = row;
+    const double sign = phi[c][pivot] < 0.0 ? -1.0 : 1.0;
+    for (std::size_t row = 0; row < phi[c].Size(); ++row)
+      state.vectors(c, dofs.FreeDofs()[row]) = sign * phi[c][row];
+  }
+  state.backend.eigenvectors = Matrix(0U, 0U);
+  return Result<SolidFrequencyState>::Success(std::move(state));
+}
+}  // namespace
+
+Result<SolidFrequencyState> SolidFrequencyProcedureExecutor::Execute(
+    const Domain& domain) const {
+  try {
+    auto result = ExecuteCandidate(domain);
+    if (result.HasValue()) return result;
+    auto diagnostics = result.GetStatus().Diagnostics();
+    for (auto& diagnostic : diagnostics)
+      if (diagnostic.location.file.empty())
+        diagnostic.location = domain.Frequency() != nullptr
+                                  ? domain.Frequency()->frequency_location
+                                  : SourceLocation{domain.SourcePath(), 0U};
+    return Result<SolidFrequencyState>::Failure(Status::Failure(
+        result.GetStatus().Category().value_or(FailureCategory::kSolver),
+        std::move(diagnostics)));
+  } catch (const std::bad_alloc&) {
+    return Result<SolidFrequencyState>::Failure(Failure(
+        "eigen-allocation-failed", "Solid modal candidate allocation failed."));
+  } catch (const std::length_error&) {
+    return Result<SolidFrequencyState>::Failure(
+        Failure("eigen-allocation-failed",
+                "Solid modal workspace exceeds addressable storage."));
+  }
+}
+}  // namespace fesa

src/fesa/elements/hex8_element.cpp

@@ -172,6 +172,13 @@ const ElementDofLayout& Hex8Element::DofLayout() const noexcept {
 Result<ElementStiffnessContribution> Hex8Element::ComputeStiffness() const {
   return Result<ElementStiffnessContribution>::Success({layout_, stiffness_});
 }
+std::array<double, 8> Hex8Element::ReferenceWeights() const noexcept {
+  std::array<double, 8> weights{};
+  for (std::size_t g = 0; g < weights.size(); ++g)
+    weights[g] = points_[g].weight;
+  return weights;
+}
+
 Result<ElementStiffnessContribution> Hex8Element::ComputeMass() const {
   if (!density_)
     return Result<ElementStiffnessContribution>::Failure(

src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp

@@ -17,6 +17,8 @@
 #include <vector>
 
 #include "fesa/math/sparse_matrix.h"
+#include "fesa/math/vector.h"
+#include "fesa/solvers/linear/mkl_pardiso_solver.h"
 #include "solvers/eigen/mkl_feast_eigen_solver_test_seam.h"
 
 namespace fesa {
@@ -364,15 +366,57 @@ bool FillInitialBasis(std::vector<double>& vectors, const std::size_t n,
 
 }  // namespace
 
+Result<FeastEigenResult> MklFeastEigenSolver::SolveGeneralized(
+    const SparseMatrix& stiffness, const SparseMatrix& mass,
+    const std::size_t requested_modes, const double mass_lower_bound,
+    const double classification_scale) const {
+  try {
+    return SolveImpl(stiffness, &mass, requested_modes, classification_scale,
+                     mass_lower_bound);
+  } catch (const std::bad_alloc&) {
+    return Result<FeastEigenResult>::Failure(
+        Failure("eigen-allocation-failed", "generalized-workspace",
+                "Allocation of generalized sparse-pencil storage failed."));
+  }
+}
+
 Result<FeastEigenResult> MklFeastEigenSolver::Solve(
     const SparseMatrix& whitened_operator, const std::size_t requested_modes,
     const double zero_operator_scale) const {
+  return SolveImpl(whitened_operator, nullptr, requested_modes,
+                   zero_operator_scale, 1.0);
+}
+
+Result<FeastEigenResult> MklFeastEigenSolver::SolveImpl(
+    const SparseMatrix& whitened_operator, const SparseMatrix* mass,
+    const std::size_t requested_modes, const double zero_operator_scale,
+    const double mass_lower_bound) const {
+  static_assert(sizeof(MKL_INT) == 8U,
+                "FEAST requires the installed ILP64 ABI.");
   const Status csr_status = whitened_operator.Validate();
   if (!csr_status.IsOk())
     return Result<FeastEigenResult>::Failure(
         Failure("feast-invalid-csr", "public-csr",
                 "The public sparse matrix failed CSR validation."));
   const std::size_t n = whitened_operator.Rows();
+  if (mass != nullptr) {
+    if (!mass->Validate().IsOk() || mass->Rows() != n || mass->Columns() != n ||
+        !std::isfinite(mass_lower_bound) || !(mass_lower_bound > 0.0) ||
+        !std::isfinite(zero_operator_scale) || !(zero_operator_scale > 0.0)) {
+      return Result<FeastEigenResult>::Failure(Failure(
+          "feast-invalid-mass-bound", "generalized-pencil",
+          "A matching finite mass and certified positive bound are required."));
+    }
+    const Status mass_symmetry = ValidateSymmetry(*mass);
+    if (!mass_symmetry.IsOk())
+      return Result<FeastEigenResult>::Failure(mass_symmetry);
+    MklPardisoSolver spd_check;
+    const Status spd = spd_check.Factorize(*mass);
+    if (!spd.IsOk())
+      return Result<FeastEigenResult>::Failure(
+          Failure("feast-mass-not-positive-definite", "mass-spd",
+                  "Sparse SPD mass factorization failed."));
+  }
   if (n == 0U || whitened_operator.Columns() != n || requested_modes == 0U ||
       requested_modes > n) {
     return Result<FeastEigenResult>::Failure(Failure(
@@ -399,7 +443,8 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
   for (std::size_t row = 0U; row < n; ++row) {
     double sum = 0.0;
     for (std::size_t p = offsets[row]; p < offsets[row + 1U]; ++p) {
-      sum += std::abs(values[p]);
+      sum = mass != nullptr ? NextUp(sum + std::abs(values[p]))
+                            : sum + std::abs(values[p]);
     }
     if (!std::isfinite(sum)) {
       return Result<FeastEigenResult>::Failure(Failure(
@@ -416,7 +461,7 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
         "An exactly zero operator requires a finite positive physical scale."));
   }
   const double operator_scale =
-      zero_operator ? zero_operator_scale : infinity_norm;
+      mass != nullptr || zero_operator ? zero_operator_scale : infinity_norm;
   if (!std::isfinite(operator_scale) || !(operator_scale > 0.0)) {
     return Result<FeastEigenResult>::Failure(
         Failure("feast-nonfinite-operator-scale", "operator-scale",
@@ -424,7 +469,8 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
   }
   auto a = CopyUpperOneBased(whitened_operator, operator_scale);
   if (!a.HasValue()) return Result<FeastEigenResult>::Failure(a.GetStatus());
-  auto b = IdentityOneBased(n);
+  auto b =
+      mass != nullptr ? CopyUpperOneBased(*mass, 1.0) : IdentityOneBased(n);
   if (!b.HasValue()) return Result<FeastEigenResult>::Failure(b.GetStatus());
 
   double lower = -1.0;
@@ -432,7 +478,17 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
   double g_lower = 0.0;
   double g_upper = 0.0;
   double beta = 1.0;
-  if (!zero_operator) {
+  if (mass != nullptr) {
+    // ||K||2 <= ||K||inf and M >= m_lower*I bound the genuine pencil.
+    const double bound =
+        NextUp(NextUp(infinity_norm / mass_lower_bound) / operator_scale);
+    g_lower = -bound;
+    g_upper = bound;
+    beta = (std::max)(1.0, bound);
+    lower = std::nextafter(g_lower - 16.0 * Delta(g_lower, beta),
+                           -(std::numeric_limits<double>::infinity)());
+    upper = NextUp(g_upper + 16.0 * Delta(g_upper, beta));
+  } else if (!zero_operator) {
     g_lower = (std::numeric_limits<double>::infinity)();
     g_upper = -(std::numeric_limits<double>::infinity)();
     const auto& columns = whitened_operator.ColumnIndices();
@@ -467,6 +523,18 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
 
   LocalThreadGuard threads;
   std::vector<FeastBackendCall> calls;
+  if (mass != nullptr)
+    calls.push_back({0U, "mass-spd", "original-mass-factorization", 0.0,
+                     mass_lower_bound, 0U, 0, false, 0U, false, 0, false, 0.0,
+                     "accepted",
+                     "sparse-SPD; caller-certified-roundoff-bound"});
+  auto TraceFailure = [&](const std::string& code, const std::string& identity,
+                          const std::string& message) {
+    return Failure(code, identity,
+                   mass != nullptr
+                       ? message + " Backend calls=" + CallHistory(calls)
+                       : message);
+  };
   std::array<MKL_INT, 128> snapshot{};
   SetControls(snapshot, false);
   std::vector<std::int64_t> base_fpm;
@@ -482,16 +550,16 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
     count_res.assign(n, 0.0);
     count_x.assign(n * n, 0.0);
   } catch (const std::bad_alloc&) {
-    return Result<FeastEigenResult>::Failure(
-        Failure("eigen-allocation-failed", "count-workspace",
-                "Allocation of the required FEAST count workspace failed."));
+    return Result<FeastEigenResult>::Failure(TraceFailure(
+        "eigen-allocation-failed", "count-workspace",
+        "Allocation of the required FEAST count workspace failed."));
   }
 
   auto Count = [&](const double call_lower, const double call_upper,
                    const std::string& purpose, std::size_t* count) -> Status {
     if (!internal::CanCallFeast(calls.size()))
-      return Failure("eigen-search-budget-exhausted", "count-calls",
-                     "The bounded FEAST count-call budget was exhausted.");
+      return TraceFailure("eigen-search-budget-exhausted", "count-calls",
+                          "The bounded FEAST count-call budget was exhausted.");
     std::array<MKL_INT, 128> fpm{};
     SetControls(fpm, true);
     MKL_INT m0 = static_cast<MKL_INT>(n);
@@ -500,7 +568,7 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
     MKL_INT info = 0;
     double epsout = 0.0;
     if (!FillInitialBasis(count_x, n, n)) {
-      return Failure(
+      return TraceFailure(
           "eigenpair-validation-failed", "count-initial-basis",
           "The prescribed count initial basis failed QR validation.");
     }
@@ -540,11 +608,12 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
         count_available ? "accepted" : "vendor-info=" + std::to_string(info);
     calls.push_back(std::move(event));
     if (info != 0)
-      return Failure("eigen-feast-failed", "info=" + std::to_string(info),
-                     "FEAST did not certify the requested count interval.");
+      return TraceFailure(
+          "eigen-feast-failed", "info=" + std::to_string(info),
+          "FEAST did not certify the requested count interval.");
     if (m0 < 0 || static_cast<std::size_t>(m0) > n) {
-      return Failure("eigen-count-inconsistent", "count-output",
-                     "FEAST count output is outside [0,n].");
+      return TraceFailure("eigen-count-inconsistent", "count-output",
+                          "FEAST count output is outside [0,n].");
     }
     *count = static_cast<std::size_t>(m0);
     for (std::size_t i = 0U; i + 1U < calls.size(); ++i) {
@@ -555,7 +624,7 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
       if ((prior.upper == call_upper && prior.count != *count) ||
           (prior.upper < call_upper && prior.count > *count) ||
           (prior.upper > call_upper && prior.count < *count))
-        return Failure(
+        return TraceFailure(
             "eigen-count-inconsistent", "count-history",
             "Successful prefix counts violate identity or monotonicity.");
     }
@@ -599,10 +668,36 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
     return Result<FeastEigenResult>::Failure(
         !count_status.IsOk()
             ? count_status
-            : Failure(
+            : TraceFailure(
                   "eigen-count-inconsistent", "total-enclosure",
                   "The enclosing FEAST count did not equal matrix dimension."));
   }
+  if (mass != nullptr) {
+    // Count the actual pencil before removing any part of the conservative
+    // negative enclosure. A far-negative endpoint degrades contour filtering
+    // of low positive modes on a stiff mesh; an empty guarded prefix proves
+    // that moving it cannot remove a mode, including a repeated rigid cluster.
+    const double margin = Delta(0.0, beta);
+    const double candidate = -16.0 * margin;
+    const std::array<double, 3U> guards{candidate - margin, candidate,
+                                        candidate + margin};
+    if (lower < guards.front()) {
+      for (const double guard : guards) {
+        std::size_t negative_count = 0U;
+        const Status prefix_status =
+            Count(lower, guard, "lower-empty-prefix", &negative_count);
+        if (!prefix_status.IsOk())
+          return Result<FeastEigenResult>::Failure(prefix_status);
+        if (negative_count != 0U)
+          return Result<FeastEigenResult>::Failure(TraceFailure(
+              "negative-eigenvalue-spectrum", "certified-negative-prefix",
+              "The generalized pencil contains significant negative "
+              "eigenvalues."));
+      }
+      lower = candidate;
+      calls.back().decision_reason = "guarded-empty-prefix; lowest-spectrum-preserved";
+    }
+  }
   const double fixed_lower = lower;
   const double enclosure_upper = upper;
   double solve_upper = enclosure_upper;
@@ -694,7 +789,7 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
   count_status = Count(fixed_lower, solve_upper, "final-count", &count);
   if (!count_status.IsOk() || count < requested_modes ||
       count != certified_count) {
-    return Result<FeastEigenResult>::Failure(Failure(
+    return Result<FeastEigenResult>::Failure(TraceFailure(
         "eigen-count-inconsistent", "final-count",
         "The final FEAST interval count is not the certified guarded count."));
   }
@@ -709,27 +804,27 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
     if (!below_status.IsOk() || !above_status.IsOk() || below != count ||
         above != count) {
       return Result<FeastEigenResult>::Failure(
-          Failure("eigen-count-inconsistent", "final-guard",
-                  "The final interior endpoint is not guarded by equal "
-                  "successful counts."));
+          TraceFailure("eigen-count-inconsistent", "final-guard",
+                       "The final interior endpoint is not guarded by equal "
+                       "successful counts."));
     }
   }
   if (count > (std::numeric_limits<std::size_t>::max)() - 8U ||
       count > ((std::numeric_limits<std::size_t>::max)() - 1U) / 3U) {
     return Result<FeastEigenResult>::Failure(
-        Failure("eigen-allocation-failed", "solve-capacity",
-                "FEAST solve-capacity padding overflows size_t."));
+        TraceFailure("eigen-allocation-failed", "solve-capacity",
+                     "FEAST solve-capacity padding overflows size_t."));
   }
   const std::size_t padded_count =
       (std::max)(count + 8U, (3U * count + 1U) / 2U);
   std::size_t capacity = count == n ? n : (std::min)(n, padded_count);
   if (!internal::ProductFits(n, capacity))
     return Result<FeastEigenResult>::Failure(
-        Failure("feast-workspace-overflow", "solve-workspace",
-                "The FEAST solve workspace size overflows."));
+        TraceFailure("feast-workspace-overflow", "solve-workspace",
+                     "The FEAST solve workspace size overflows."));
   if (!internal::FitsDoubleStorage(n * capacity) ||
       !internal::FitsDoubleStorage(capacity)) {
-    return Result<FeastEigenResult>::Failure(Failure(
+    return Result<FeastEigenResult>::Failure(TraceFailure(
         "eigen-allocation-failed", "solve-workspace",
         "The required solve workspace exceeds addressable double storage."));
   }
@@ -742,23 +837,23 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
   std::size_t solve_capacity = capacity;
   for (std::size_t attempt = 0U; attempt < 3U; ++attempt) {
     if (!internal::CanCallFeast(calls.size()))
-      return Result<FeastEigenResult>::Failure(
-          Failure("eigen-search-budget-exhausted", "solve-call",
-                  "The total backend call budget is exhausted before solve."));
+      return Result<FeastEigenResult>::Failure(TraceFailure(
+          "eigen-search-budget-exhausted", "solve-call",
+          "The total backend call budget is exhausted before solve."));
     try {
       internal::ApplyFeastWorkspaceTestHook(false);
       eigenvalues.assign(solve_capacity, 0.0);
       residuals.assign(solve_capacity, 0.0);
       vectors.assign(n * solve_capacity, 0.0);
     } catch (const std::bad_alloc&) {
-      return Result<FeastEigenResult>::Failure(
-          Failure("eigen-allocation-failed", "solve-workspace",
-                  "Allocation of the required FEAST solve workspace failed."));
+      return Result<FeastEigenResult>::Failure(TraceFailure(
+          "eigen-allocation-failed", "solve-workspace",
+          "Allocation of the required FEAST solve workspace failed."));
     }
     if (!FillInitialBasis(vectors, n, solve_capacity)) {
-      return Result<FeastEigenResult>::Failure(
-          Failure("eigenpair-validation-failed", "solve-initial-basis",
-                  "The prescribed solve initial basis failed QR validation."));
+      return Result<FeastEigenResult>::Failure(TraceFailure(
+          "eigenpair-validation-failed", "solve-initial-basis",
+          "The prescribed solve initial basis failed QR validation."));
     }
     std::array<MKL_INT, 128> fpm{};
     SetControls(fpm, false);
@@ -794,21 +889,21 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
                                  : "vendor-info=" + std::to_string(info)});
     if (info == 0 && (!std::isfinite(epsout) || loop < 0)) {
       return Result<FeastEigenResult>::Failure(
-          Failure("invalid-eigenvalue-result", "solve-output",
-                  "A successful FEAST solve returned nonfinite epsout or "
-                  "invalid loop."));
+          TraceFailure("invalid-eigenvalue-result", "solve-output",
+                       "A successful FEAST solve returned nonfinite epsout or "
+                       "invalid loop."));
     }
     if (info == 0 && returned >= 0 &&
         static_cast<std::size_t>(returned) <= solve_capacity) {
       for (std::size_t column = 0U; column < static_cast<std::size_t>(returned);
            ++column) {
         if (!std::isfinite(eigenvalues[column]))
-          return Result<FeastEigenResult>::Failure(Failure(
+          return Result<FeastEigenResult>::Failure(TraceFailure(
               "invalid-eigenvalue-result", std::to_string(column),
               "A nonfinite eigenvalue cannot trigger an interval retry."));
         for (std::size_t row = 0U; row < n; ++row)
           if (!std::isfinite(vectors[row + column * n]))
-            return Result<FeastEigenResult>::Failure(Failure(
+            return Result<FeastEigenResult>::Failure(TraceFailure(
                 "invalid-eigenvalue-result", std::to_string(column),
                 "A nonfinite eigenvector cannot trigger an interval retry."));
       }
@@ -818,10 +913,10 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
     if (endpoint_separated) {
       for (std::size_t index = 0U; index < count; ++index) {
         if (!std::isfinite(eigenvalues[index])) {
-          return Result<FeastEigenResult>::Failure(
-              Failure("invalid-eigenvalue-result", std::to_string(index),
-                      "FEAST returned a nonfinite eigenvalue before endpoint "
-                      "validation."));
+          return Result<FeastEigenResult>::Failure(TraceFailure(
+              "invalid-eigenvalue-result", std::to_string(index),
+              "FEAST returned a nonfinite eigenvalue before endpoint "
+              "validation."));
         }
         endpoint_separated =
             endpoint_separated &&
@@ -856,7 +951,7 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
       const Status recount_status =
           Count(fixed_lower, solve_upper, "final-count", &recounted);
       if (!recount_status.IsOk() || recounted < requested_modes) {
-        return Result<FeastEigenResult>::Failure(Failure(
+        return Result<FeastEigenResult>::Failure(TraceFailure(
             "eigen-count-inconsistent", "solve-retry",
             "A larger certified solve interval could not be recounted."));
       }
@@ -871,7 +966,7 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
             fixed_lower, solve_upper + retry_margin, "final-count", &above);
         if (!below_status.IsOk() || !above_status.IsOk() || below != count ||
             above != count) {
-          return Result<FeastEigenResult>::Failure(Failure(
+          return Result<FeastEigenResult>::Failure(TraceFailure(
               "eigen-count-inconsistent", "solve-retry-guard",
               "The larger retry interval is not guarded by matching counts."));
         }
@@ -879,8 +974,8 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
       if (count > (std::numeric_limits<std::size_t>::max)() - 8U ||
           count > ((std::numeric_limits<std::size_t>::max)() - 1U) / 3U) {
         return Result<FeastEigenResult>::Failure(
-            Failure("eigen-allocation-failed", "solve-capacity",
-                    "FEAST retry capacity padding overflows size_t."));
+            TraceFailure("eigen-allocation-failed", "solve-capacity",
+                         "FEAST retry capacity padding overflows size_t."));
       }
       const std::size_t retry_padding =
           (std::max)(count + 8U, (3U * count + 1U) / 2U);
@@ -894,22 +989,22 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
   }
   if (!solve_accepted) {
     return Result<FeastEigenResult>::Failure(
-        Failure("eigen-feast-failed", "info=" + std::to_string(info),
-                "FEAST solve did not return the certified interval count "
-                "within its bounded retries."));
+        TraceFailure("eigen-feast-failed", "info=" + std::to_string(info),
+                     "FEAST solve did not return the certified interval count "
+                     "within its bounded retries."));
   }
 
   for (std::size_t index = 0U; index < count; ++index) {
     if (!std::isfinite(eigenvalues[index])) {
-      return Result<FeastEigenResult>::Failure(
-          Failure("invalid-eigenvalue-result", std::to_string(index),
-                  "FEAST returned a nonfinite eigenvalue before sorting."));
+      return Result<FeastEigenResult>::Failure(TraceFailure(
+          "invalid-eigenvalue-result", std::to_string(index),
+          "FEAST returned a nonfinite eigenvalue before sorting."));
     }
     for (std::size_t row = 0U; row < n; ++row) {
       if (!std::isfinite(vectors[row + index * n])) {
-        return Result<FeastEigenResult>::Failure(
-            Failure("invalid-eigenvalue-result", std::to_string(index),
-                    "FEAST returned a nonfinite eigenvector before sorting."));
+        return Result<FeastEigenResult>::Failure(TraceFailure(
+            "invalid-eigenvalue-result", std::to_string(index),
+            "FEAST returned a nonfinite eigenvector before sorting."));
       }
     }
   }
@@ -936,28 +1031,34 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
     const double mu = eigenvalues[input];
     if (!std::isfinite(mu)) {
       return Result<FeastEigenResult>::Failure(
-          Failure("invalid-eigenvalue-result", std::to_string(input),
-                  "FEAST returned a nonfinite eigenvalue."));
+          TraceFailure("invalid-eigenvalue-result", std::to_string(input),
+                       "FEAST returned a nonfinite eigenvalue."));
     }
     if (mu < -1.0e-9) {
-      return Result<FeastEigenResult>::Failure(Failure(
+      return Result<FeastEigenResult>::Failure(TraceFailure(
           "negative-eigenvalue-spectrum", std::to_string(input),
           "The raw normalized FEAST eigenvalue is significantly negative."));
     }
     if (mu - fixed_lower < Delta(fixed_lower, beta) ||
         solve_upper - mu < Delta(solve_upper, beta)) {
-      return Result<FeastEigenResult>::Failure(Failure(
+      return Result<FeastEigenResult>::Failure(TraceFailure(
           "eigen-endpoint-separation", std::to_string(input),
           "A returned FEAST eigenvalue is too close to a solve endpoint."));
     }
     double norm = 0.0;
     for (std::size_t row = 0U; row < n; ++row)
       norm += vectors[row + input * n] * vectors[row + input * n];
+    if (mass != nullptr) {
+      Vector raw(n);
+      for (std::size_t row = 0U; row < n; ++row)
+        raw[row] = vectors[row + input * n];
+      norm = raw.Dot(mass->Multiply(raw));
+    }
     norm = std::sqrt(norm);
     if (!(norm > 0.0) || !std::isfinite(norm))
       return Result<FeastEigenResult>::Failure(
-          Failure("eigenpair-validation-failed", std::to_string(input),
-                  "FEAST returned a non-normalizable eigenvector."));
+          TraceFailure("eigenpair-validation-failed", std::to_string(input),
+                       "FEAST returned a non-normalizable eigenvector."));
     double residual_norm = 0.0;
     for (std::size_t row = 0U; row < n; ++row) {
       double ax = 0.0;
@@ -965,12 +1066,43 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
         ax += values[p] / operator_scale *
               vectors[whitened_operator.ColumnIndices()[p] + input * n];
       }
-      const double difference = ax - mu * vectors[row + input * n];
+      double mx = vectors[row + input * n];
+      if (mass != nullptr) {
+        mx = 0.0;
+        for (std::size_t p = mass->RowOffsets()[row];
+             p < mass->RowOffsets()[row + 1U]; ++p)
+          mx +=
+              mass->Values()[p] * vectors[mass->ColumnIndices()[p] + input * n];
+      }
+      const double difference = ax - mu * mx;
       residual_norm += difference * difference;
       ordered_vectors(row, output) = vectors[row + input * n] / norm;
     }
-    const double rho_b = std::sqrt(residual_norm) / norm;
+    double residual_denominator = norm;
+    if (mass != nullptr) {
+      double norm_k = 0.0;
+      double norm_m = 0.0;
+      double euclidean = 0.0;
+      for (const double value : values)
+        norm_k = std::hypot(norm_k, value / operator_scale);
+      for (const double value : mass->Values())
+        norm_m = std::hypot(norm_m, value);
+      for (std::size_t row = 0U; row < n; ++row)
+        euclidean = std::hypot(euclidean, vectors[row + input * n]);
+      residual_denominator = (norm_k + std::abs(mu) * norm_m) * euclidean;
+    }
+    const double rho_b =
+        residual_denominator == 0.0
+            ? (residual_norm == 0.0 ? 0.0
+                                    : (std::numeric_limits<double>::infinity)())
+            : std::sqrt(residual_norm) / residual_denominator;
     if (rho_b > 1.0e-10) {
+      if (mass != nullptr)
+        return Result<FeastEigenResult>::Failure(TraceFailure(
+            "eigenpair-validation-failed", std::to_string(input),
+            "Independent original-pencil residual exceeds 1e-10: eta=" +
+                FormatDouble(rho_b) + ", mu=" + FormatDouble(mu) +
+                "; calls=" + CallHistory(calls)));
       const FeastBackendCall& last_call = calls.back();
       std::vector<double> mirrored_action(n, 0.0);
       std::vector<double> full_action(n, 0.0);
@@ -999,7 +1131,7 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
       const double rho_b_upper = std::sqrt(mirrored_residual_squared) / norm;
       const double action_difference_norm =
           std::sqrt(action_difference_squared);
-      return Result<FeastEigenResult>::Failure(Failure(
+      return Result<FeastEigenResult>::Failure(TraceFailure(
           "eigenpair-validation-failed", std::to_string(input),
           "Independent whitened residual exceeds 1e-10: rho_b=" +
               FormatDouble(rho_b) + ", mu=" + FormatDouble(mu) +
@@ -1019,16 +1151,16 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
     }
     const double lambda = mu * operator_scale;
     if (!std::isfinite(lambda) || (mu != 0.0 && lambda == 0.0)) {
-      return Result<FeastEigenResult>::Failure(
-          Failure("invalid-eigenvalue-result", std::to_string(input),
-                  "Physical eigenvalue conversion overflowed or underflowed."));
+      return Result<FeastEigenResult>::Failure(TraceFailure(
+          "invalid-eigenvalue-result", std::to_string(input),
+          "Physical eigenvalue conversion overflowed or underflowed."));
     }
     ordered_values.push_back(lambda);
     const bool backend_res_finite = std::isfinite(residuals[input]);
     if (!backend_res_finite && mu != 0.0) {
-      return Result<FeastEigenResult>::Failure(
-          Failure("invalid-eigenvalue-result", std::to_string(input),
-                  "A nonzero FEAST mode has a nonfinite backend residual."));
+      return Result<FeastEigenResult>::Failure(TraceFailure(
+          "invalid-eigenvalue-result", std::to_string(input),
+          "A nonzero FEAST mode has a nonfinite backend residual."));
     }
     ordered_backend_available.push_back(backend_res_finite);
     ordered_backend_res.push_back(backend_res_finite ? residuals[input] : 0.0);
@@ -1042,23 +1174,35 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
       double dot = 0.0;
       for (std::size_t row = 0U; row < n; ++row)
         dot += ordered_vectors(row, left) * ordered_vectors(row, right);
+      if (mass != nullptr) {
+        dot = 0.0;
+        for (std::size_t row = 0U; row < n; ++row) {
+          double action = 0.0;
+          for (std::size_t p = mass->RowOffsets()[row];
+               p < mass->RowOffsets()[row + 1U]; ++p)
+            action += mass->Values()[p] *
+                      ordered_vectors(mass->ColumnIndices()[p], right);
+          dot += ordered_vectors(row, left) * action;
+        }
+      }
       const double expected = left == right ? 1.0 : 0.0;
       const double difference = std::abs(dot - expected);
       gram_squared += difference * difference;
       gram_error = (std::max)(gram_error, difference);
       if (std::abs(dot - expected) > 1.0e-10)
-        return Result<FeastEigenResult>::Failure(Failure(
+        return Result<FeastEigenResult>::Failure(TraceFailure(
             "eigenpair-validation-failed",
             std::to_string(left) + ":" + std::to_string(right),
             "The all-c FEAST vectors fail the independent Gram check."));
     }
   }
-  gram_error = std::sqrt(gram_squared / static_cast<double>(count));
+  gram_error = std::sqrt(gram_squared /
+                         (mass != nullptr ? 1.0 : static_cast<double>(count)));
   if (gram_error > 1.0e-10) {
-    return Result<FeastEigenResult>::Failure(
-        Failure("eigenpair-validation-failed", "gram",
-                "The normalized all-c Frobenius Gram error exceeds 1e-10: " +
-                    FormatDouble(gram_error) + "."));
+    return Result<FeastEigenResult>::Failure(TraceFailure(
+        "eigenpair-validation-failed", "gram",
+        "The normalized all-c Frobenius Gram error exceeds 1e-10: " +
+            FormatDouble(gram_error) + "."));
   }
   char version[256]{};
   mkl_get_version_string(version, static_cast<int>(sizeof(version)));