Add Dlaev2

native/dlaev2.go

@@ -0,0 +1,76 @@
+package native
+
+import "math"
+
+// Dlaev2 computes the eigendecomposition of a symmetric 2×2 matrix.
+// The matrix is given by
+//  [a b]
+//  [b c]
+// Dlaev2 returns rt1 and rt2, the eigenvalues of the matrix where |RT1| > |RT2|,
+// and [cs1, sn1] which is the unit right eigenvalue for RT1.
+//  [ cs1 sn1] [a b] [cs1 -sn1] = [rt1   0]
+//  [-sn1 cs1] [b c] [sn1  cs1]   [  0 rt2]
+func (impl Implementation) Dlaev2(a, b, c float64) (rt1, rt2, cs1, sn1 float64) {
+	sm := a + c
+	df := a - c
+	adf := math.Abs(df)
+	tb := b + b
+	ab := math.Abs(tb)
+	acmx := c
+	acmn := a
+	if math.Abs(a) > math.Abs(c) {
+		acmx = a
+		acmn = c
+	}
+	var rt float64
+	if adf > ab {
+		rt = adf * math.Sqrt(1+(ab/adf)*(ab/adf))
+	} else if adf < ab {
+		rt = ab * math.Sqrt(1+(adf/ab)*(adf/ab))
+	} else {
+		rt = ab * math.Sqrt(2)
+	}
+	var sgn1 float64
+	if sm < 0 {
+		rt1 = 0.5 * (sm - rt)
+		sgn1 = -1
+		rt2 = (acmx/rt1)*acmn - (b/rt1)*b
+	} else if sm > 0 {
+		rt1 = 0.5 * (sm + rt)
+		sgn1 = 1
+		rt2 = (acmx/rt1)*acmn - (b/rt1)*b
+	} else {
+		rt1 = 0.5 * rt
+		rt2 = -0.5 * rt
+		sgn1 = 1
+	}
+	var cs, sgn2 float64
+	if df >= 0 {
+		cs = df + rt
+		sgn2 = 1
+	} else {
+		cs = df - rt
+		sgn2 = -1
+	}
+	acs := math.Abs(cs)
+	if acs > ab {
+		ct := -tb / cs
+		sn1 = 1 / math.Sqrt(1+ct*ct)
+		cs1 = ct * sn1
+	} else {
+		if ab == 0 {
+			cs1 = 1
+			sn1 = 0
+		} else {
+			tn := -cs / tb
+			cs1 = 1 / math.Sqrt(1+tn*tn)
+			sn1 = tn * cs1
+		}
+	}
+	if sgn1 == sgn2 {
+		tn := cs1
+		cs1 = -sn1
+		sn1 = tn
+	}
+	return rt1, rt2, cs1, sn1
+}

native/lapack_test.go

@@ -48,6 +48,10 @@ func TestDgetrs(t *testing.T) {
 	testlapack.DgetrsTest(t, impl)
 }
 
+func TestDlaev2(t *testing.T) {
+	testlapack.Dlaev2Test(t, impl)
+}
+
 func TestDlange(t *testing.T) {
 	testlapack.DlangeTest(t, impl)
 }

testlapack/dlaev2.go

@@ -0,0 +1,41 @@
+package testlapack
+
+import (
+	"math"
+	"math/rand"
+	"testing"
+)
+
+type Dlaev2er interface {
+	Dlaev2(a, b, c float64) (rt1, rt2, cs1, sn1 float64)
+}
+
+func Dlaev2Test(t *testing.T, impl Dlaev2er) {
+	for trial := 0; trial < 100; trial++ {
+		a := rand.NormFloat64()
+		b := rand.NormFloat64()
+		c := rand.NormFloat64()
+
+		rt1, rt2, cs1, sn1 := impl.Dlaev2(a, b, c)
+		tmp := mul2by2([2][2]float64{{cs1, sn1}, {-sn1, cs1}}, [2][2]float64{{a, b}, {b, c}})
+		ans := mul2by2(tmp, [2][2]float64{{cs1, -sn1}, {sn1, cs1}})
+		if math.Abs(ans[0][0]-rt1) > 1e-14 {
+			t.Errorf("Largest eigenvalue mismatch. Returned %v, mul %v", rt1, ans[0][0])
+		}
+		if math.Abs(ans[1][0]) > 1e-14 || math.Abs(ans[0][1]) > 1e-14 {
+			t.Errorf("Non-zero off diagonal. ans[1][0] = %v, ans[0][1] = %v", ans[1][0], ans[0][1])
+		}
+		if math.Abs(ans[1][1]-rt2) > 1e-14 {
+			t.Errorf("Smallest eigenvalue mismatch. Returned %v, mul %v", rt2, ans[1][1])
+		}
+	}
+}
+
+func mul2by2(a, b [2][2]float64) [2][2]float64 {
+	var c [2][2]float64
+	c[0][0] = a[0][0]*b[0][0] + a[0][1]*b[1][0]
+	c[0][1] = a[0][0]*b[0][1] + a[0][1]*b[1][1]
+	c[1][0] = a[1][0]*b[0][0] + a[1][1]*b[1][0]
+	c[1][1] = a[1][0]*b[0][1] + a[1][1]*b[1][1]
+	return c
+}