lapack/gonum: propagate NaN from Dcombssq
MODIFY
lapack/gonum/dcombssq.go
lapack/gonum/dcombssq.go
package gonum
+import "math"
+
// Dcombssq adds two scaled sum-of-squares quantities, V := V1 + V2,
// V_scale^2 * V_ssq := V1_scale^2 * V1_ssq + V2_scale^2 * V2_ssq
// and returns the result V.
@@ -14,7 +16,10 @@ func (Implementation) Dcombssq(scale1, ssq1, scale2, ssq2 float64) (scale, ssq f
if scale1 != 0 {
return scale1, ssq1 + (scale2/scale1)*(scale2/scale1)*ssq2
}
- // If the input is non-negative and we are here, then scale2 must inevitably be 0, too.
+ // Both scales are zero.
+ if math.IsNaN(ssq1) || math.IsNaN(ssq2) {
+ return 0, math.NaN()
+ }
return 0, 0
}
return scale2, ssq2 + (scale1/scale2)*(scale1/scale2)*ssq1

MODIFY
lapack/testlapack/dcombssq.go
lapack/testlapack/dcombssq.go
@@ -19,13 +19,28 @@ func DcombssqTest(t *testing.T, impl Dcombssqer) {
const tol = 1e-15
rnd := rand.New(rand.NewSource(1))
- for i := 0; i < 10; i++ {
- // Generate random, non-negative input.
+ for i := 0; i < 100; i++ {
+ // Generate random, non-negative input, with an occasional NaN.
+ var hasNaN bool
scale1 := rnd.Float64()
- ssq1 := rnd.Float64()
+ var ssq1 float64
+ if rnd.Float64() < 0.5 {
+ ssq1 = math.NaN()
+ hasNaN = true
+ } else {
+ ssq1 = rnd.Float64()
+ }
+
scale2 := rnd.Float64()
- ssq2 := rnd.Float64()
- switch i {
+ var ssq2 float64
+ if rnd.Float64() < 0.5 {
+ ssq2 = math.NaN()
+ hasNaN = true
+ } else {
+ ssq2 = rnd.Float64()
+ }
+
+ switch rnd.Intn(4) {
case 0:
scale1 = 0
case 1:
@@ -38,6 +53,13 @@ func DcombssqTest(t *testing.T, impl Dcombssqer) {
// scale^2 * ssq := scale1^2 * ssq1 + scale2^2 * ssq2
scale, ssq := impl.Dcombssq(scale1, ssq1, scale2, ssq2)
+ if hasNaN {
+ if !math.IsNaN(ssq) {
+ t.Errorf("Case %v: unexpected ssq; got %v, want NaN", i, ssq)
+ }
+ continue
+ }
+
// Compute the expected result in a non-sophisticated way and
// compare against the result we got.
want := scale1*scale1*ssq1 + scale2*scale2*ssq2

lapack/gonum: propagate NaN from Dcombssq