Repository navigation
RDCON - BUGFIX - Galerkin factorization failures: Cholesky falls back to LU, a failed LU stops - #301
Merged
Merged
Conversation
…rf info) gal_solve ignored info from zgbtrf/zpbtrf, so an indefinite matrix with solver='cholesky' gave a garbage Delta' with no warning (g147131 n=1, TPS settings: 2/1 = 9.7e6). Now program_stop with the info value. LU results unchanged (bit-identical Delta' on g147131). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…Galerkin matrix is not positive definite zpbtrf info was ignored, so solver='cholesky' on an indefinite matrix gave a garbage Delta' (g147131 n=1, TPS settings: 2/1 = 9.7e6). Now, when zpbtrf fails, gal_solve warns, sets solver='LU', reassembles the full band from the element matrices (gal_assemble_mat) and solves it by LU (2/1 = 2.1359; solver='LU' gives 2.1296, it does not zero the boundary columns). zgbtrf failure stops the run. solver='LU' runs and Cholesky runs that succeed are bit-identical. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This was referenced Oct 13, 2026
StuartBenjamin
marked this pull request as draft
October 13, 2026 07:24
StuartBenjamin
requested review from
matt-pharr
and removed request for
d-burg
October 13, 2026 08:08
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
StuartBenjamin
marked this pull request as ready for review
October 13, 2026 09:24
matt-pharr
approved these changes
Oct 15, 2026
matt-pharr
left a comment
matt-pharr
left a comment
Collaborator
There was a problem hiding this comment.
Good very simple change, I just added a check to switch back to LU in case someone tries to run cholesky with rpec_flag=T, which incorrectly gives Delta_coil=0 presently.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
gal_solve(rdcon/gal.f) never checkedinfofromzpbtrf(Cholesky, the default) orzgbtrf(LU).zpbtrfreads only the lower band and the real part of the diagonal (there is no Hermitian check; the resonant diagonal entries carry a small imaginary numerical error). When that matrix is not positive definite,zpbtrfstops partway andzpbtrssolves with the partial factor: Δ′ is garbage and nothing is printed.Example: g147131, n = 1, tearing-physics-suite (TPS) settings (no wall, nx 64, sing_order 8, psihigh 0.999). Develop's Cholesky gives Δ′(2/1, 3/1) = 9.7e6, 7.7e7; LU gives 2.13, −14.3.
Fix
info, switch to LU, and reassemble the full band from the element matrices (gal_assemble_mat), so LU uses the upper band and the complex diagonal, as a direct LU run does. The element matrices already carry the boundary conditions, sogal_set_boundaryis not rerun. A successful Cholesky run does no extra work.info.This follows matt-pharr's suggestion for Julia GPEC #513.
Tests
develop 1055d1b vs this branch:
After the fallback, 2/1 differs from a direct LU run by 0.3 % (3/1 by 0.07 %). Both solve the same equations, but this matrix is indefinite and badly conditioned, so round-off is amplified. With #303 merged too, this case is positive definite and no longer falls back.
Scripts: galerkin_factorization_scripts.zip (needs
g147131.02300_DIIID_KEFIT, not shipped):No dependencies; merges cleanly with #303.
🤖 Generated with Claude Code