Add drscl and tests, and fix the dlamch parameters while here.

native/drscl.go

@@ -0,0 +1,41 @@
+package native
+
+import (
+	"math"
+
+	"github.com/gonum/blas/blas64"
+)
+
+// Drscl multiplies the vector x by 1/a being careful to avoid overflow or
+// underflow where possible.
+func (impl Implementation) Drscl(n int, a float64, x []float64, incX int) {
+	checkVector(n, x, incX)
+	bi := blas64.Implementation()
+	cden := a
+	cnum := 1.0
+	smlnum := dlamchS
+	bignum := 1 / smlnum
+	for {
+		cden1 := cden * smlnum
+		cnum1 := cnum / bignum
+		var mul float64
+		var done bool
+		switch {
+		case cnum != 0 && math.Abs(cden1) > math.Abs(cnum):
+			mul = smlnum
+			done = false
+			cden = cden1
+		case math.Abs(cnum1) > math.Abs(cden):
+			mul = bignum
+			done = false
+			cnum = cnum1
+		default:
+			mul = cnum / cden
+			done = true
+		}
+		bi.Dscal(n, mul, x, incX)
+		if done {
+			break
+		}
+	}
+}

native/general.go

@@ -78,18 +78,15 @@ func max(a, b int) int {
 // fixed values. Probably a way to get them as constants.
 // TODO(btracey): Is there a better way to find the smallest number such that 1+E > 1
 
-var dlamchE, dlamchS, dlamchP float64
+var (
+	// dlamchE is the machine epsilon. For IEEE this is 2^-53.
+	dlamchE = math.Float64frombits(0x3ca0000000000000)
 
-func init() {
+	// dlamchP is 2 * eps
-	onePlusEps := math.Nextafter(1, math.Inf(1))
+	dlamchP = math.Float64frombits(0x3cb0000000000000)
-	eps := (math.Nextafter(1, math.Inf(1)) - 1) * 0.5
+
-	dlamchE = eps
+	// dlamchS is the "safe min", that is, the lowest number such that 1/sfmin does
-	sfmin := math.SmallestNonzeroFloat64
+	// not overflow. The Netlib code for calculating this number is not correct --
-	small := 1 / math.MaxFloat64
+	// it overflows. Found by trial and error, it is equal to (1/math.MaxFloat64) * (1+ 6*eps)
-	if small >= sfmin {
+	dlamchS = math.Float64frombits(0x4000000000001)
-		sfmin = small * onePlusEps
+)
-	}
-	dlamchS = sfmin
-	radix := 2.0
-	dlamchP = radix * eps
-}

native/lapack_test.go

@@ -88,6 +88,10 @@ func TestDpotrf(t *testing.T) {
 	testlapack.DpotrfTest(t, impl)
 }
 
+func TestDrscl(t *testing.T) {
+	testlapack.DrsclTest(t, impl)
+}
+
 func TestIladlc(t *testing.T) {
 	testlapack.IladlcTest(t, impl)
 }

testlapack/drscl.go

@@ -0,0 +1,48 @@
+package testlapack
+
+import (
+	"math"
+	"testing"
+
+	"github.com/gonum/floats"
+)
+
+type Drscler interface {
+	Drscl(n int, a float64, x []float64, incX int)
+}
+
+func DrsclTest(t *testing.T, impl Drscler) {
+	for _, test := range []struct {
+		x []float64
+		a float64
+	}{
+		{
+			x: []float64{1, 2, 3, 4, 5},
+			a: 4,
+		},
+		{
+			x: []float64{1, 2, 3, 4, 5},
+			a: math.MaxFloat64,
+		},
+		{
+			x: []float64{1, 2, 3, 4, 5},
+			a: 1e-307,
+		},
+	} {
+		xcopy := make([]float64, len(test.x))
+		copy(xcopy, test.x)
+
+		// Cannot test the scaling directly because of floating point scaling issues
+		// (the purpose of Drscl). Instead, check that scaling and scaling back
+		// yeilds approximately x. If overflow or underflow occurs then the scaling
+		// won't match.
+		impl.Drscl(len(test.x), test.a, xcopy, 1)
+		if floats.Equal(xcopy, test.x) {
+			t.Errorf("x unchanged during call to drscl. a = %v, x = %v.", test.a, test.x)
+		}
+		impl.Drscl(len(test.x), 1/test.a, xcopy, 1)
+		if !floats.EqualApprox(xcopy, test.x, 1e-14) {
+			t.Errorf("x not equal after scaling and unscaling. a = %v, x = %v.", test.a, test.x)
+		}
+	}
+}