From b4846de641e04bc5c4c29fd950b19831e5d57728 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 28 Oct 2026 16:28:19 +0100 Subject: [PATCH] ModelScaleGPU: a Levenberg-Marquardt fit whose J^T J has underflowed to 0 ends at its best so far, as on the CPU Since fdafc009b J^T J and J^T r are summed in float. A fit that walks to a nonsensical scale (k_overall < 0, B of thousands of A^2, as fits of an unrelated model can in the solvent grid) makes |F| k_aniso small enough that every square underflows to exactly 0 in float, while J^T r keeps a residue near 1e-36. JordanSolve then threw "Trying to reverse singular matrix. Column 0 is zeroed." out of ValidateAgainstModel (ModelValidation_CCModelUnrelatedModelConfirmsNothing). In double the squares stayed near 1e-72 and the solve went through. LevMarFit in ModelScaling.cpp already stops at alpha[0] == 0 for the same case in double; the GPU run now does the same check before the solve. Only fits that would have thrown are affected; no kernel changes, deterministic as before. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- image_analysis/structure_refinement/ModelScaleGPU.cu | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/image_analysis/structure_refinement/ModelScaleGPU.cu b/image_analysis/structure_refinement/ModelScaleGPU.cu index de11ff941..35ebf187f 100644 --- a/image_analysis/structure_refinement/ModelScaleGPU.cu +++ b/image_analysis/structure_refinement/ModelScaleGPU.cu @@ -373,6 +373,14 @@ struct ModelScaleGPU::LevMarRun { done = true; return; } + // As LevMarFit (ModelScaling.cpp): a fit that has walked to where k_overall exp(-B s^2) |F| underflows at + // every reflection has no step left - the squares in alpha are all 0 while beta keeps a residue, and the + // solve throws on that. It ends here, at its best so far. J^T J is summed in float, so its squares + // underflow once |F| k_aniso is below about 1e-22 everywhere, still far from any scale a fit can end at. + if (alpha[0] == 0) { + done = true; + return; + } const size_t na = beta.size(); std::vector temp_alpha = alpha; for (size_t j = 0; j < na; j++)