lapack/gonum: improve accuracy of Dlartg
This incorporates modifications introduced in LAPACK 3.11.0 in
MODIFY
lapack/gonum/dlartg.go
lapack/gonum/dlartg.go
@@ -13,48 +13,61 @@ import "math"
//
// where cs*cs + sn*sn = 1.
//
-// This is a more accurate version of BLAS Drotg, with the other differences
-// that
+// This is a more accurate version of BLAS Drotg that uses scaling to avoid
+// overflow or underflow, with the other differences that
+// - cs >= 0
// - if g = 0, then cs = 1 and sn = 0
-// - if f = 0 and g != 0, then cs = 0 and sn = 1
-// - r takes the sign of f and so cs is always non-negative
+// - if f = 0 and g != 0, then cs = 0 and sn = sign(1,g)
//
// Dlartg is an internal routine. It is exported for testing purposes.
func (impl Implementation) Dlartg(f, g float64) (cs, sn, r float64) {
// Implementation based on Supplemental Material to:
- // Edward Anderson. 2017. Algorithm 978: Safe Scaling in the Level 1 BLAS.
- // ACM Trans. Math. Softw. 44, 1, Article 12 (July 2017), 28 pages.
+ //
+ // Edward Anderson
+ // Algorithm 978: Safe Scaling in the Level 1 BLAS
+ // ACM Trans. Math. Softw. 44, 1, Article 12 (2017)
// DOI: https://doi.org/10.1145/3061665
+ //
+ // For further details see:
+ //
+ // W. Pereira, A. Lotfi, J. Langou
+ // Numerical analysis of Givens rotation
+ // DOI: https://doi.org/10.48550/arXiv.2211.04010
+
+ if g == 0 {
+ return 1, 0, f
+ }
+
+ g1 := math.Abs(g)
+
+ if f == 0 {
+ return 0, math.Copysign(1, g), g1
+ }
+
const safmin = dlamchS
const safmax = 1 / safmin
+ rtmin := math.Sqrt(safmin)
+ rtmax := math.Sqrt(safmax / 2)
+
f1 := math.Abs(f)
- g1 := math.Abs(g)
- switch {
- case g == 0:
- cs = 1
- sn = 0
- r = f
- case f == 0:
- cs = 0
- sn = math.Copysign(1, g)
- r = g1
- case drtmin < f1 && f1 < drtmax && drtmin < g1 && g1 < drtmax:
+
+ if rtmin < f1 && f1 < rtmax && rtmin < g1 && g1 < rtmax {
d := math.Sqrt(f*f + g*g)
- p := 1 / d
- cs = f1 * p
- sn = g * math.Copysign(p, f)
+ cs = f1 / d
r = math.Copysign(d, f)
- default:
- maxfg := math.Max(f1, g1)
- u := math.Min(math.Max(safmin, maxfg), safmax)
- uu := 1 / u
- fs := f * uu
- gs := g * uu
- d := math.Sqrt(fs*fs + gs*gs)
- p := 1 / d
- cs = math.Abs(fs) * p
- sn = gs * math.Copysign(p, f)
- r = math.Copysign(d, f) * u
+ sn = g / r
+
+ return cs, sn, r
}
+
+ u := math.Min(math.Max(safmin, math.Max(f1, g1)), safmax)
+ fs := f / u
+ gs := g / u
+ d := math.Sqrt(fs*fs + gs*gs)
+ cs = math.Abs(fs) / d
+ r = math.Copysign(d, f)
+ sn = gs / r
+ r *= u
+
return cs, sn, r
}

MODIFY
lapack/gonum/lapack.go
lapack/gonum/lapack.go
@@ -49,12 +49,6 @@ const (
// For IEEE this is 2^{-1022}.
dlamchS = 0x1p-1022
- // (rtmin,rtmax) is a range of well-scaled numbers whose square
- // or sum of squares is also safe.
- // drtmin is sqrt(dlamchS/dlamchP)
- drtmin = 0x1p-485
- drtmax = 1 / drtmin
-
// Blue's scaling constants
//
// An n-vector x is well-scaled if

lapack/gonum: improve accuracy of Dlartg

This incorporates modifications introduced in LAPACK 3.11.0 in
https://github.com/Reference-LAPACK/lapack/pull/631