From 63d82fd754182b7799c86812e294549ea2571ec3 Mon Sep 17 00:00:00 2001
From: btracey <tracey.brendan@gmail.com>
Date: Mon, 4 Jan 2016 11:52:31 -0700
Subject: [PATCH] Add dlae2 and test

---
 native/dlae2.go       | 47 ++++++++++++++++++++++++++++++++++++++
 native/lapack_test.go |  4 ++++
 testlapack/dlae2.go   | 53 +++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 104 insertions(+)
 create mode 100644 native/dlae2.go
 create mode 100644 testlapack/dlae2.go

diff --git a/native/dlae2.go b/native/dlae2.go
new file mode 100644
index 00000000..03ca0a12
--- /dev/null
+++ b/native/dlae2.go
@@ -0,0 +1,47 @@
+// Copyright ©2016 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"
+
+// Dlae2 computes the eigenvalues of a 2×2 symmetric matrix
+//  [a b]
+//  [b c]
+// and returns the eigenvalue with the larger absolute value as rt1 and the
+// smaller as rt2.
+func (impl Implementation) Dlae2(a, b, c float64) (rt1, rt2 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.Sqrt2
+	}
+	if sm < 0 {
+		rt1 = 0.5 * (sm - rt)
+		rt2 = (acmx/rt1)*acmn - (b/rt1)*b
+		return rt1, rt2
+	}
+	if sm > 0 {
+		rt1 = 0.5 * (sm + rt)
+		rt2 = (acmx/rt1)*acmn - (b/rt1)*b
+		return rt1, rt2
+	}
+	rt1 = 0.5 * rt
+	rt2 = -0.5 * rt
+	return rt1, rt2
+}
diff --git a/native/lapack_test.go b/native/lapack_test.go
index 60387b3c..d2f9d954 100644
--- a/native/lapack_test.go
+++ b/native/lapack_test.go
@@ -72,6 +72,10 @@ func TestDlacpy(t *testing.T) {
 	testlapack.DlacpyTest(t, impl)
 }
 
+func TestDlae2(t *testing.T) {
+	testlapack.Dlae2Test(t, impl)
+}
+
 func TestDlaev2(t *testing.T) {
 	testlapack.Dlaev2Test(t, impl)
 }
diff --git a/testlapack/dlae2.go b/testlapack/dlae2.go
new file mode 100644
index 00000000..236b3f2b
--- /dev/null
+++ b/testlapack/dlae2.go
@@ -0,0 +1,53 @@
+// Copyright ©2016 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 (
+	"fmt"
+	"math"
+	"testing"
+)
+
+type Dlae2er interface {
+	Dlae2(a, b, c float64) (rt1, rt2 float64)
+}
+
+func Dlae2Test(t *testing.T, impl Dlae2er) {
+	for _, test := range []struct {
+		a, b, c float64
+	}{
+		{-10, 5, 3},
+		{3, 5, -10},
+		{0, 3, 0},
+		{1, 3, 1},
+		{1, -3, 1},
+		{5, 0, 3},
+		{3, 0, -5},
+		{1, 3, 1.02},
+		{1.02, 3, 1},
+		{1, -3, -9},
+	} {
+		a := test.a
+		b := test.b
+		c := test.c
+		rt1, rt2 := impl.Dlae2(a, b, c)
+
+		errStr := fmt.Sprintf("a = %v, b = %v, c = %v", a, b, c)
+		// Check if rt1 and rt2 are eigenvalues by checking if det(a - λI) = 0
+		a1 := a - rt1
+		c1 := c - rt1
+		det := a1*c1 - b*b
+		if math.Abs(det) > 1e-10 {
+			t.Errorf("First eigenvalue mismatch. %s. Det = %v", errStr, det)
+		}
+
+		a2 := a - rt2
+		c2 := c - rt2
+		det = a2*c2 - b*b
+		if math.Abs(det) > 1e-10 {
+			t.Errorf("Second eigenvalue mismatch. %s. Det = %v", errStr, det)
+		}
+	}
+}