Add Dlartg and test

native/dlartg.go

@@ -0,0 +1,80 @@
+// Copyright ©2015 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package native
+
+import "math"
+
+// Dlartg generates a plane rotation so that
+//  [cs  sn] * [f] = [r]
+//  [-sn cs]   [g] = [0]
+// This is a more accurate version of BLAS drotg, with the other differences that
+// if g = 0, then cs = 1 and sn = 0, and if f = 0 and g != 0, then cs = 0 and sn = 1.
+// If abs(f) > abs(g), cs will be positive.
+func (impl Implementation) Dlartg(f, g float64) (cs, sn, r float64) {
+	safmin := dlamchS
+	eps := dlamchE
+	safmn2 := math.Pow(float64(dlamchB), math.Log(safmin/eps))
+	safmx2 := 1 / safmn2
+	if g == 0 {
+		cs = 1
+		sn = 1
+		r = f
+		return cs, sn, r
+	}
+	if f == 0 {
+		cs = 0
+		sn = 1
+		r = g
+		return cs, sn, r
+	}
+	f1 := f
+	g1 := g
+	scale := math.Max(math.Abs(f1), math.Abs(g1))
+	if scale >= safmx2 {
+		var count int
+		for {
+			count++
+			f1 *= safmn2
+			g1 *= safmn2
+			scale = math.Max(math.Abs(f1), math.Abs(g1))
+			if scale < safmx2 {
+				break
+			}
+		}
+		r = math.Sqrt(f1*f1 + g1*g1)
+		cs = f1 / r
+		sn = g1 / r
+		for i := 0; i < count; i++ {
+			r *= safmx2
+		}
+	} else if scale <= safmn2 {
+		var count int
+		for {
+			count++
+			f1 *= safmx2
+			g1 *= safmx2
+			scale = math.Max(math.Abs(f1), math.Abs(g1))
+			if scale >= safmn2 {
+				break
+			}
+		}
+		r = math.Sqrt(f1*f1 + g1*g1)
+		cs = f1 / r
+		sn = g1 / r
+		for i := 0; i < count; i++ {
+			r *= safmn2
+		}
+	} else {
+		r = math.Sqrt(f1*f1 + g1*g1)
+		cs = f1 / r
+		sn = g1 / r
+	}
+	if math.Abs(f) > math.Abs(g) && cs < 0 {
+		cs *= -1
+		sn *= -1
+		r *= -1
+	}
+	return cs, sn, r
+}

native/general.go

@@ -90,4 +90,7 @@ var (
 
 	smlnum = dlamchS / dlamchP
 	bignum = 1 / smlnum
+
+	// dlamchB is the radix of the machine (the base of the number system).
+	dlamchB = 2
 )

native/lapack_test.go

@@ -88,6 +88,10 @@ func TestDlarft(t *testing.T) {
 	testlapack.DlarftTest(t, impl)
 }
 
+func TestDlartg(t *testing.T) {
+	testlapack.DlartgTest(t, impl)
+}
+
 func TestDorml2(t *testing.T) {
 	testlapack.Dorml2Test(t, impl)
 }

testlapack/dlartg.go

@@ -0,0 +1,32 @@
+// Copyright ©2015 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package testlapack
+
+import (
+	"math"
+	"math/rand"
+	"testing"
+)
+
+type Dlartger interface {
+	Dlartg(f, g float64) (cs, sn, r float64)
+}
+
+func DlartgTest(t *testing.T, impl Dlartger) {
+	for i := 0; i < 100; i++ {
+		f := rand.NormFloat64()
+		g := rand.NormFloat64()
+		cs, sn, r := impl.Dlartg(f, g)
+
+		rTest := cs*f + sn*g
+		zeroTest := -sn*f + cs*g
+		if math.Abs(rTest-r) > 1e-14 {
+			t.Errorf("R result mismatch. Computed %v, found %v", r, rTest)
+		}
+		if math.Abs(zeroTest) > 1e-14 {
+			t.Errorf("Zero result mismatch. Found %v", zeroTest)
+		}
+	}
+}