diff --git a/native/dlartg.go b/native/dlartg.go new file mode 100644 index 00000000..ac8eaeca --- /dev/null +++ b/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 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 +} diff --git a/native/general.go b/native/general.go index a8b6d935..bf6a6d98 100644 --- a/native/general.go +++ b/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 ) diff --git a/native/lapack_test.go b/native/lapack_test.go index 54c4e546..fdd34dc7 100644 --- a/native/lapack_test.go +++ b/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) } diff --git a/testlapack/dlartg.go b/testlapack/dlartg.go new file mode 100644 index 00000000..a29abd42 --- /dev/null +++ b/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) + } + } +}