Add dlasv2 and tests

native/dlasv2.go

@@ -0,0 +1,113 @@
+// 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"
+
+// Dlasv2 computes the singular value decomposition of a 2×2 matrix.
+//  [ csl snl] [f g] [csr -snr] = [ssmax     0]
+//  [-snl csl] [0 h] [snr  csr] = [    0 ssmin]
+// ssmax is the larger absolute singular value, and ssmin is the smaller absolute
+// singular value. [cls, snl] and [csr, snr] are the left and right singular vectors.
+func (impl Implementation) Dlasv2(f, g, h float64) (ssmin, ssmax, snr, csr, snl, csl float64) {
+	ft := f
+	fa := math.Abs(ft)
+	ht := h
+	ha := math.Abs(h)
+	// pmax points to the largest element of the matrix in terms of absolute value.
+	// 1 if F, 2 if G, 3 if H.
+	pmax := 1
+	swap := ha > fa
+	if swap {
+		pmax = 3
+		ft, ht = ht, ft
+		fa, ha = ha, fa
+	}
+	gt := g
+	ga := math.Abs(gt)
+	var clt, crt, slt, srt float64
+	if ga == 0 {
+		ssmin = ha
+		ssmax = fa
+		clt = 1
+		crt = 1
+		slt = 0
+		srt = 0
+	} else {
+		gasmall := true
+		if ga > fa {
+			pmax = 2
+			if (fa / ga) < dlamchE {
+				gasmall = false
+				ssmax = ga
+				if ha > 1 {
+					ssmin = fa / (ga / ha)
+				} else {
+					ssmin = (fa / ga) * ha
+				}
+				clt = 1
+				slt = ht / gt
+				srt = 1
+				crt = ft / gt
+			}
+		}
+		if gasmall {
+			d := fa - ha
+			l := d / fa
+			if d == fa { // deal with inf
+				l = 1
+			}
+			m := gt / ft
+			t := 2 - l
+			s := math.Hypot(t, m)
+			var r float64
+			if l == 0 {
+				r = math.Abs(m)
+			} else {
+				r = math.Hypot(l, m)
+			}
+			a := 0.5 * (s + r)
+			ssmin = ha / a
+			ssmax = fa * a
+			if m == 0 {
+				if l == 0 {
+					t = math.Copysign(2, ft) * math.Copysign(1, gt)
+				} else {
+					t = gt/math.Copysign(d, ft) + m/t
+				}
+			} else {
+				t = (m/(s+t) + m/(r+l)) * (1 + a)
+			}
+			l = math.Hypot(t, 2)
+			crt = 2 / l
+			srt = t / l
+			clt = (crt + srt*m) / a
+			slt = (ht / ft) * srt / a
+		}
+	}
+	if swap {
+		csl = srt
+		snl = crt
+		csr = slt
+		snr = clt
+	} else {
+		csl = clt
+		snl = slt
+		csr = crt
+		snr = srt
+	}
+	var tsign float64
+	switch pmax {
+	case 1:
+		tsign = math.Copysign(1, csr) * math.Copysign(1, csl) * math.Copysign(1, f)
+	case 2:
+		tsign = math.Copysign(1, snr) * math.Copysign(1, csl) * math.Copysign(1, g)
+	case 3:
+		tsign = math.Copysign(1, snr) * math.Copysign(1, snl) * math.Copysign(1, h)
+	}
+	ssmax = math.Copysign(ssmax, tsign)
+	ssmin = math.Copysign(ssmin, tsign*math.Copysign(1, f)*math.Copysign(1, h))
+	return ssmin, ssmax, snr, csr, snl, csl
+}

native/lapack_test.go

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

testlapack/dlasv2.go

@@ -0,0 +1,47 @@
+// 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/rand"
+	"testing"
+
+	"github.com/gonum/floats"
+)
+
+type Dlasv2er interface {
+	Dlasv2(f, g, h float64) (ssmin, ssmax, snr, csr, snl, csl float64)
+}
+
+func Dlasv2Test(t *testing.T, impl Dlasv2er) {
+	for i := 0; i < 100; i++ {
+		f := rand.NormFloat64()
+		g := rand.NormFloat64()
+		h := rand.NormFloat64()
+
+		ssmin, ssmax, snr, csr, snl, csl := impl.Dlasv2(f, g, h)
+
+		// tmp =
+		// [ csl snl] [f g]
+		// [-snl csl] [0 h]
+		tmp11 := csl * f
+		tmp12 := csl*g + snl*h
+		tmp21 := -snl * f
+		tmp22 := -snl*g + csl*h
+		// lhs =
+		// [tmp11 tmp12] [csr -snr]
+		// [tmp21 tmp22] [snr  csr]
+		ans11 := tmp11*csr + tmp12*snr
+		ans12 := tmp11*-snr + tmp12*csr
+		ans21 := tmp21*csr + tmp22*snr
+		ans22 := tmp21*-snr + tmp22*csr
+
+		lhs := []float64{ans11, ans12, ans21, ans22}
+		rhs := []float64{ssmax, 0, 0, ssmin}
+		if !floats.EqualApprox(rhs, lhs, 1e-12) {
+			t.Errorf("SVD mismatch. f = %v, g = %v, h = %v.\nLHS: %v\nRHS: %v", f, g, h, lhs, rhs)
+		}
+	}
+}