lapack/gonum: rewrite Dlassq and its test
MODIFY
lapack/gonum/dlassq.go
lapack/gonum/dlassq.go
@@ -6,13 +6,16 @@ package gonum
import "math"
-// Dlassq updates a sum of squares in scaled form. The input parameters scale and
-// sumsq represent the current scale and total sum of squares. These values are
-// updated with the information in the first n elements of the vector specified
-// by x and incX.
+// Dlassq updates a sum of squares represented in scaled form. Dlassq returns
+// the values scl and smsq such that
+// scl^2*smsq = X[0]^2 + ... + X[n-1]^2 + scale^2*sumsq
//
// Dlassq is an internal routine. It is exported for testing purposes.
func (impl Implementation) Dlassq(n int, x []float64, incx int, scale float64, sumsq float64) (scl, smsq 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.
+ // DOI: https://doi.org/10.1145/3061665
switch {
case n < 0:
panic(nLT0)
@@ -22,20 +25,94 @@ func (impl Implementation) Dlassq(n int, x []float64, incx int, scale float64, s
panic(shortX)
}
+ if math.IsNaN(scale) || math.IsNaN(sumsq) {
+ return scale, sumsq
+ }
+
+ if sumsq == 0 {
+ scale = 1
+ }
+ if scale == 0 {
+ scale = 1
+ sumsq = 0
+ }
+
if n == 0 {
return scale, sumsq
}
- for ix := 0; ix <= (n-1)*incx; ix += incx {
- absxi := math.Abs(x[ix])
- if absxi > 0 || math.IsNaN(absxi) {
- if scale < absxi {
- sumsq = 1 + sumsq*(scale/absxi)*(scale/absxi)
- scale = absxi
- } else {
- sumsq += (absxi / scale) * (absxi / scale)
+ // Compute the sum of squares in 3 accumulators:
+ // - abig: sum of squares scaled down to avoid overflow
+ // - asml: sum of squares scaled up to avoid underflow
+ // - amed: sum of squares that do not require scaling
+ // The thresholds and multipliers are:
+ // - values bigger than dtbig are scaled down by dsbig
+ // - values smaller than dtsml are scaled up by dssml
+ var (
+ isBig bool
+ asml, amed, abig float64
+ )
+ for i, ix := 0, 0; i < n; i++ {
+ ax := math.Abs(x[ix])
+ switch {
+ case ax > dtbig:
+ ax *= dsbig
+ abig += ax * ax
+ isBig = true
+ case ax < dtsml:
+ if !isBig {
+ ax *= dssml
+ asml += ax * ax
+ }
+ default:
+ amed += ax * ax
}
+ ix += incx
+ }
+ // Put the existing sum of squares into one of the accumulators.
+ if sumsq > 0 {
+ ax := scale * math.Sqrt(sumsq)
+ switch {
+ case ax > dtbig:
+ ax *= dsbig
+ abig += ax * ax
+ case ax < dtsml:
+ if !isBig {
+ ax *= dssml
+ asml += ax * ax
+ }
+ default:
+ amed += ax * ax
+ }
+ }
+ // Combine abig and amed or amed and asml if more than one accumulator was
+ // used.
+ switch {
+ case abig > 0:
+ // Combine abig and amed:
+ if amed > 0 || math.IsNaN(amed) {
+ abig += (amed * dsbig) * dsbig
+ }
+ scale = 1 / dsbig
+ sumsq = abig
+ case asml > 0:
+ // Combine amed and asml:
+ if amed > 0 || math.IsNaN(amed) {
+ amed = math.Sqrt(amed)
+ asml = math.Sqrt(asml) / dssml
+ ymin, ymax := asml, amed
+ if asml > amed {
+ ymin, ymax = amed, asml
+ }
+ scale = 1
+ sumsq = ymax * ymax * (1 + (ymin/ymax)*(ymin/ymax))
+ } else {
+ scale = 1 / dssml
+ sumsq = asml
}
+ default:
+ scale = 1
+ sumsq = amed
}
return scale, sumsq
}

MODIFY
lapack/testlapack/dlassq.go
lapack/testlapack/dlassq.go
@@ -19,46 +19,151 @@ type Dlassqer interface {
}
func DlassqTest(t *testing.T, impl Dlassqer) {
- const tol = 1e-14
+ const (
+ safmin = dlamchS
+ safmax = 1 / safmin
+ ulp = dlamchP
+ smlnum = safmin / ulp
+ bignum = safmax * ulp
+ )
+ values := []float64{
+ 0,
+ 2 * safmin,
+ smlnum,
+ ulp,
+ 1,
+ 1 / ulp,
+ bignum,
+ safmax,
+ math.Inf(1),
+ math.NaN(),
+ }
rnd := rand.New(rand.NewSource(1))
- for _, n := range []int{0, 1, 2, 3, 4, 5, 10} {
+ for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 20, 30, 40} {
for _, incx := range []int{1, 3} {
- name := fmt.Sprintf("n=%v,incx=%v", n, incx)
+ for cas := 0; cas < 3; cas++ {
+ for _, v0 := range values {
+ if v0 > 1 {
+ v0 *= 0.5
+ }
+ for _, v1 := range values {
+ if v1 > 1 {
+ v1 = 0.5 * v1 / math.Sqrt(float64(n+1))
+ }
+ dlassqTest(t, impl, rnd, n, incx, cas, v0, v1)
+ }
+ }
+ }
+ }
+ }
+}
- // Allocate a slice of minimum length and fill it with
- // random numbers.
- x := make([]float64, max(0, 1+(n-1)*incx))
- for i := range x {
- x[i] = rnd.Float64()
+func dlassqTest(t *testing.T, impl Dlassqer, rnd *rand.Rand, n, incx, cas int, v0, v1 float64) {
+ const (
+ rogue = 1234.5678
+ tol = 1e-14
+ )
+
+ name := fmt.Sprintf("n=%v,incx=%v,cas=%v,v0=%v,v1=%v", n, incx, cas, v0, v1)
+
+ // Generate n random values in (-1,1].
+ work := make([]float64, n)
+ for i := range work {
+ work[i] = 1 - 2*rnd.Float64()
}
- // Fill the referenced elements of x and compute the
- // expected result in a non-sophisticated way.
- scale := rnd.Float64()
- ssq := rnd.Float64()
- want := scale * scale * ssq
- for i := 0; i < n; i++ {
- xi := rnd.NormFloat64()
- x[i*incx] = xi
- want += xi * xi
+ // Compute the sum of squares by an unscaled algorithm.
+ var workssq float64
+ for _, wi := range work {
+ workssq += wi * wi
}
+ // Construct a test (n+1)-vector z as
+ // z[0] = known value
+ // z[1:n+1] = scaling factor * work[:]
+ z := make([]float64, n+1)
+ z[0] = math.Sqrt(v0)
+ for i, wi := range work {
+ z[i+1] = v1 * wi
+ }
+ // Set initial scale and ssq corresponding to z[0].
+ var scale, ssq float64
+ switch cas {
+ case 0:
+ scale = 1
+ ssq = v0
+ case 1:
+ scale = math.Sqrt(v0)
+ ssq = 1
+ case 2:
+ if v0 < 1 {
+ scale = 1.0 / 3
+ ssq = 9 * v0
+ } else {
+ scale = 3
+ ssq = v0 / 9
+ }
+ default:
+ panic("bad cas")
+ }
+
+ // Compute the expected value of the sum of squares using the (n+1)-vector z.
+ z0 := z[0]
+ var z1n float64
+ if n >= 1 {
+ z1n = v1 * math.Sqrt(workssq)
+ }
+ zmin := math.Min(z0, z1n)
+ zmax := math.Max(z0, z1n)
+ var nrmWant float64
+ switch {
+ case math.IsNaN(z0) || math.IsNaN(z1n):
+ nrmWant = math.NaN()
+ case zmin == zmax:
+ nrmWant = math.Sqrt2 * zmax
+ case zmax == 0:
+ nrmWant = 0
+ default:
+ nrmWant = zmax * math.Sqrt(1+(zmin/zmax)*(zmin/zmax))
+ }
+
+ // Allocate input slice for Dlassq and fill it with z[1:].
+ x := make([]float64, max(0, 1+(n-1)*incx))
+ for i := range x {
+ x[i] = rogue
+ }
+ for i, zi := range z[1:] {
+ x[i*incx] = zi
+ }
xCopy := make([]float64, len(x))
copy(xCopy, x)
- // Update scale and ssq so that
- // scale_out^2 * ssq_out = x[0]^2 + ... + x[n-1]^2 + scale_in^2*ssq_in
- scale, ssq = impl.Dlassq(n, x, incx, scale, ssq)
- if !floats.Equal(x, xCopy) {
+ scaleGot, ssqGot := impl.Dlassq(n, x, incx, scale, ssq)
+ nrmGot := scaleGot * math.Sqrt(ssqGot)
+
+ if !floats.Same(x, xCopy) {
t.Fatalf("%v: unexpected modification of x", name)
}
// Check the result.
- got := scale * scale * ssq
- if math.Abs(got-want) >= tol {
- t.Errorf("%v: unexpected result; got %v, want %v", name, got, want)
+ switch {
+ case math.IsNaN(nrmGot) || math.IsNaN(nrmWant):
+ if !math.IsNaN(nrmGot) {
+ t.Errorf("%v: expected NaN; got %v", name, nrmGot)
+ }
+ if !math.IsNaN(nrmWant) {
+ t.Errorf("%v: unexpected NaN; want %v", name, nrmWant)
+ }
+ case nrmGot == nrmWant:
+ case nrmWant == 0:
+ if nrmGot > tol {
+ t.Errorf("%v: unexpected result; got %v, want 0", name, nrmGot)
}
+ default:
+ diff := math.Abs(nrmGot-nrmWant) / nrmWant / math.Max(1, float64(n))
+ if math.IsNaN(diff) || diff > tol {
+ t.Errorf("%v: unexpected result; got %v, want %v", name, nrmGot, nrmWant)
}
}
}

lapack/gonum: rewrite Dlassq and its test