lapack/gonum: remove Dcombssq
Remove Dcombssq and rely on Dlassq to combine scaled sums. See the
upstream change in
MODIFY
lapack/gonum/dlangb.go
lapack/gonum/dlangb.go
@@ -74,15 +74,14 @@ func (impl Implementation) Dlangb(norm lapack.MatrixNorm, m, n, kl, ku int, ab [
}
case lapack.Frobenius:
scale := 0.0
- ssq := 1.0
+ sum := 1.0
for i := 0; i < min(m, n+kl); i++ {
l := max(0, kl-i)
u := min(n+kl-i, ncol)
ilen := u - l
- rowscale, rowssq := impl.Dlassq(ilen, ab[i*ldab+l:], 1, 0, 1)
- scale, ssq = impl.Dcombssq(scale, ssq, rowscale, rowssq)
+ scale, sum = impl.Dlassq(ilen, ab[i*ldab+l:], 1, scale, sum)
}
- value = scale * math.Sqrt(ssq)
+ value = scale * math.Sqrt(sum)
}
return value
}

MODIFY
lapack/gonum/dlange.go
lapack/gonum/dlange.go
@@ -80,8 +80,7 @@ func (impl Implementation) Dlange(norm lapack.MatrixNorm, m, n int, a []float64,
scale := 0.0
sum := 1.0
for i := 0; i < m; i++ {
- rowscale, rowsum := impl.Dlassq(n, a[i*lda:], 1, 0, 1)
- scale, sum = impl.Dcombssq(scale, sum, rowscale, rowsum)
+ scale, sum = impl.Dlassq(n, a[i*lda:], 1, scale, sum)
}
return scale * math.Sqrt(sum)
}

MODIFY
lapack/gonum/dlansb.go
lapack/gonum/dlansb.go
@@ -100,35 +100,31 @@ func (impl Implementation) Dlansb(norm lapack.MatrixNorm, uplo blas.Uplo, n, kd
}
case lapack.Frobenius:
scale := 0.0
- ssq := 1.0
+ sum := 1.0
if uplo == blas.Upper {
if kd > 0 {
// Sum off-diagonals.
for i := 0; i < n-1; i++ {
ilen := min(n-i-1, kd)
- rowscale, rowssq := impl.Dlassq(ilen, ab[i*ldab+1:], 1, 0, 1)
- scale, ssq = impl.Dcombssq(scale, ssq, rowscale, rowssq)
+ scale, sum = impl.Dlassq(ilen, ab[i*ldab+1:], 1, scale, sum)
}
- ssq *= 2
+ sum *= 2
}
// Sum diagonal.
- dscale, dssq := impl.Dlassq(n, ab, ldab, 0, 1)
- scale, ssq = impl.Dcombssq(scale, ssq, dscale, dssq)
+ scale, sum = impl.Dlassq(n, ab, ldab, scale, sum)
} else {
if kd > 0 {
// Sum off-diagonals.
for i := 1; i < n; i++ {
ilen := min(i, kd)
- rowscale, rowssq := impl.Dlassq(ilen, ab[i*ldab+kd-ilen:], 1, 0, 1)
- scale, ssq = impl.Dcombssq(scale, ssq, rowscale, rowssq)
+ scale, sum = impl.Dlassq(ilen, ab[i*ldab+kd-ilen:], 1, scale, sum)
}
- ssq *= 2
+ sum *= 2
}
// Sum diagonal.
- dscale, dssq := impl.Dlassq(n, ab[kd:], ldab, 0, 1)
- scale, ssq = impl.Dcombssq(scale, ssq, dscale, dssq)
+ scale, sum = impl.Dlassq(n, ab[kd:], ldab, scale, sum)
}
- value = scale * math.Sqrt(ssq)
+ value = scale * math.Sqrt(sum)
}
return value

MODIFY
lapack/gonum/dlansy.go
lapack/gonum/dlansy.go
@@ -106,23 +106,20 @@ func (impl Implementation) Dlansy(norm lapack.MatrixNorm, uplo blas.Uplo, n int,
default:
// lapack.Frobenius:
scale := 0.0
- ssq := 1.0
+ sum := 1.0
// Sum off-diagonals.
if uplo == blas.Upper {
for i := 0; i < n-1; i++ {
- rowscale, rowssq := impl.Dlassq(n-i-1, a[i*lda+i+1:], 1, 0, 1)
- scale, ssq = impl.Dcombssq(scale, ssq, rowscale, rowssq)
+ scale, sum = impl.Dlassq(n-i-1, a[i*lda+i+1:], 1, scale, sum)
}
} else {
for i := 1; i < n; i++ {
- rowscale, rowssq := impl.Dlassq(i, a[i*lda:], 1, 0, 1)
- scale, ssq = impl.Dcombssq(scale, ssq, rowscale, rowssq)
+ scale, sum = impl.Dlassq(i, a[i*lda:], 1, scale, sum)
}
}
- ssq *= 2
+ sum *= 2
// Sum diagonal.
- dscale, dssq := impl.Dlassq(n, a, lda+1, 0, 1)
- scale, ssq = impl.Dcombssq(scale, ssq, dscale, dssq)
- return scale * math.Sqrt(ssq)
+ scale, sum = impl.Dlassq(n, a, lda+1, scale, sum)
+ return scale * math.Sqrt(sum)
}
}

MODIFY
lapack/gonum/dlantb.go
lapack/gonum/dlantb.go
@@ -164,50 +164,46 @@ func (impl Implementation) Dlantb(norm lapack.MatrixNorm, uplo blas.Uplo, diag b
}
}
case lapack.Frobenius:
- var scale, ssq float64
+ var scale, sum float64
switch uplo {
case blas.Upper:
if diag == blas.Unit {
scale = 1
- ssq = float64(n)
+ sum = float64(n)
if k > 0 {
for i := 0; i < n-1; i++ {
ilen := min(n-i-1, k)
- rowscale, rowssq := impl.Dlassq(ilen, a[i*lda+1:], 1, 0, 1)
- scale, ssq = impl.Dcombssq(scale, ssq, rowscale, rowssq)
+ scale, sum = impl.Dlassq(ilen, a[i*lda+1:], 1, scale, sum)
}
}
} else {
scale = 0
- ssq = 1
+ sum = 1
for i := 0; i < n; i++ {
ilen := min(n-i, k+1)
- rowscale, rowssq := impl.Dlassq(ilen, a[i*lda:], 1, 0, 1)
- scale, ssq = impl.Dcombssq(scale, ssq, rowscale, rowssq)
+ scale, sum = impl.Dlassq(ilen, a[i*lda:], 1, scale, sum)
}
}
case blas.Lower:
if diag == blas.Unit {
scale = 1
- ssq = float64(n)
+ sum = float64(n)
if k > 0 {
for i := 1; i < n; i++ {
ilen := min(i, k)
- rowscale, rowssq := impl.Dlassq(ilen, a[i*lda+k-ilen:], 1, 0, 1)
- scale, ssq = impl.Dcombssq(scale, ssq, rowscale, rowssq)
+ scale, sum = impl.Dlassq(ilen, a[i*lda+k-ilen:], 1, scale, sum)
}
}
} else {
scale = 0
- ssq = 1
+ sum = 1
for i := 0; i < n; i++ {
ilen := min(i, k) + 1
- rowscale, rowssq := impl.Dlassq(ilen, a[i*lda+k+1-ilen:], 1, 0, 1)
- scale, ssq = impl.Dcombssq(scale, ssq, rowscale, rowssq)
+ scale, sum = impl.Dlassq(ilen, a[i*lda+k+1-ilen:], 1, scale, sum)
}
}
}
- value = scale * math.Sqrt(ssq)
+ value = scale * math.Sqrt(sum)
}
return value
}

MODIFY
lapack/gonum/dlantr.go
lapack/gonum/dlantr.go
@@ -221,36 +221,32 @@ func (impl Implementation) Dlantr(norm lapack.MatrixNorm, uplo blas.Uplo, diag b
}
default:
// lapack.Frobenius:
- var scale, ssq float64
+ var scale, sum float64
if diag == blas.Unit {
scale = 1
- ssq = float64(min(m, n))
+ sum = float64(min(m, n))
if uplo == blas.Upper {
for i := 0; i < min(m, n); i++ {
- rowscale, rowssq := impl.Dlassq(n-i-1, a[i*lda+i+1:], 1, 0, 1)
- scale, ssq = impl.Dcombssq(scale, ssq, rowscale, rowssq)
+ scale, sum = impl.Dlassq(n-i-1, a[i*lda+i+1:], 1, scale, sum)
}
} else {
for i := 1; i < m; i++ {
- rowscale, rowssq := impl.Dlassq(min(i, n), a[i*lda:], 1, 0, 1)
- scale, ssq = impl.Dcombssq(scale, ssq, rowscale, rowssq)
+ scale, sum = impl.Dlassq(min(i, n), a[i*lda:], 1, scale, sum)
}
}
} else {
scale = 0
- ssq = 1
+ sum = 1
if uplo == blas.Upper {
for i := 0; i < min(m, n); i++ {
- rowscale, rowssq := impl.Dlassq(n-i, a[i*lda+i:], 1, 0, 1)
- scale, ssq = impl.Dcombssq(scale, ssq, rowscale, rowssq)
+ scale, sum = impl.Dlassq(n-i, a[i*lda+i:], 1, scale, sum)
}
} else {
for i := 0; i < m; i++ {
- rowscale, rowssq := impl.Dlassq(min(i+1, n), a[i*lda:], 1, 0, 1)
- scale, ssq = impl.Dcombssq(scale, ssq, rowscale, rowssq)
+ scale, sum = impl.Dlassq(min(i+1, n), a[i*lda:], 1, scale, sum)
}
}
}
- return scale * math.Sqrt(ssq)
+ return scale * math.Sqrt(sum)
}
}

MODIFY
lapack/gonum/dlassq.go
lapack/gonum/dlassq.go
@@ -9,6 +9,7 @@ import "math"
// 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
+// The value of sumsq is assumed to be non-negative.
//
// 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) {
@@ -74,12 +75,18 @@ func (impl Implementation) Dlassq(n int, x []float64, incx int, scale float64, s
ax := scale * math.Sqrt(sumsq)
switch {
case ax > dtbig:
- ax *= dsbig
- abig += ax * ax
+ // We assume scale >= sqrt( TINY*EPS ) / dsbig, that is, if the
+ // scaled sum is big then its scaling factor should not be too
+ // small.
+ v := scale * dsbig
+ abig += (v * v) * sumsq
case ax < dtsml:
if !isBig {
- ax *= dssml
- asml += ax * ax
+ // We assume scale <= sqrt( HUGE ) / dssml, that is, if the
+ // scaled sum is small then its scaling factor should not be too
+ // big.
+ v := scale * dssml
+ asml += (v * v) * sumsq
}
default:
amed += ax * ax

MODIFY
lapack/gonum/lapack_test.go
lapack/gonum/lapack_test.go
@@ -17,11 +17,6 @@ func TestDbdsqr(t *testing.T) {
testlapack.DbdsqrTest(t, impl)
}
-func TestDcombssq(t *testing.T) {
- t.Parallel()
- testlapack.DcombssqTest(t, impl)
-}
-
func TestDhseqr(t *testing.T) {
t.Parallel()
testlapack.DhseqrTest(t, impl)

MODIFY
lapack/testlapack/dlassq.go
lapack/testlapack/dlassq.go
@@ -101,6 +101,21 @@ func dlassqTest(t *testing.T, impl Dlassqer, rnd *rand.Rand, n, incx, cas int, v
panic("bad cas")
}
+ const (
+ dtsml = 0x1p-511
+ dtbig = 0x1p486
+ dssml = 0x1p537
+ dsbig = 0x1p-538
+ )
+ if scale*math.Sqrt(ssq) > dtbig && scale < math.Sqrt(math.SmallestNonzeroFloat64/ulp)/dsbig {
+ // The scaled sum is big but the scale itself is small.
+ return
+ }
+ if scale*math.Sqrt(ssq) < dtsml && scale > math.Sqrt(math.MaxFloat64)/dssml {
+ // The scaled sum is small but the scale itself is big.
+ return
+ }
+
// Compute the expected value of the sum of squares using the (n+1)-vector z.
z0 := z[0]
var z1n float64

lapack/gonum: remove Dcombssq

Remove Dcombssq and rely on Dlassq to combine scaled sums. See the
upstream change in https://github.com/Reference-LAPACK/lapack/pull/570