avoid large intermediates in complex LARTG#1245
Merged
langou merged 4 commits intoJul 15, 2026
Merged
Conversation
Replace the literal 1.0 / 1.0D0 in the computation of MU = SQRT(ONE - NU**2) so the BBCSD routines consistently use the named parameter constants already defined in each source. Also reorder the parameter declarations to match the PARAMETER list and remove trailing blank lines at end of file. No functional change intended.
Rewrite the complex sine computation in CLARTG and ZLARTG to avoid forming a large intermediate value before the final multiplication by the conjugate of g.
The existing scaling keeps f2*h2 in a safe range, so this is not an overflow or underflow fix for f2*h2 itself. However, expressions such as
conjg(g) * ( f / sqrt(f2*h2) )
can still form a large intermediate value in the division before it is multiplied by conjg(g). Use the mathematically equivalent form
( f / sqrt(f2) ) * ( conjg(g) / sqrt(h2) )
so that the factors stay closer to unit scale. Apply the same rewrite to the scaled branch using fs and gs.
langou
approved these changes
Jul 15, 2026
langou
left a comment
langou
left a comment
Contributor
There was a problem hiding this comment.
This is a good idea. Thanks.
Contributor
|
Note that, in addition to the LARTG fix, there are some cosmetic changes to BBCSD too in this PR. |
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.
This PR rewrites the complex sine computation in
CLARTGandZLARTGto avoid forming unnecessarily large intermediate values.
The existing implementation uses expressions of the following form:
and similarly in the scaled branch:
The existing scaling and branch conditions already keep
f2*h2in asafe range. Therefore, this is not primarily an overflow or underflow
fix for the product
f2*h2itself.However, even when
sqrt(f2*h2)is safely representable, the divisioncan still create a large intermediate value before it is multiplied by
conjg(g). This is mathematically equivalent to multiplying normalizedfactors instead:
This PR applies that rewrite in both the unscaled and scaled branches of
CLARTGandZLARTG.The goal is to keep intermediate values closer to unit scale while
preserving the same mathematical result.
Updated files:
SRC/clartg.f90SRC/zlartg.f90