From 39ee5d4efa73ea6afc692115167db6ac5f0fde1d Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Fri, 10 Jun 2016 14:13:28 +0900
Subject: [PATCH] native:dlasy2

---
 native/dlasy2.go      |   8 ++--
 native/lapack_test.go |   4 ++
 testlapack/dlasy2.go  | 104 ++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 112 insertions(+), 4 deletions(-)
 create mode 100644 testlapack/dlasy2.go

diff --git a/native/dlasy2.go b/native/dlasy2.go
index 603a180e..cbdae6f7 100644
--- a/native/dlasy2.go
+++ b/native/dlasy2.go
@@ -21,10 +21,10 @@ import (
 // isgn must be 1 or -1, and n1 and n2 must be 0, 1, or 2, but these conditions
 // are not checked.
 //
-// Dlasy2 returns a scale factor that is chosen less than or equal to 1 to
-// prevent the solution overflowing, the infinity norm of the solution, and an
-// indicator of success. If ok is false, TL and TR have too close eigenvalues,
-// so TL or TR is perturbed to get a non-singular equation.
+// Dlasy2 returns three values, a scale factor that is chosen less than or equal
+// to 1 to prevent the solution overflowing, the infinity norm of the solution,
+// and an indicator of success. If ok is false, TL and TR have eigenvalues that
+// are too close, so TL or TR is perturbed to get a non-singular equation.
 //
 // Dlasy2 is an internal routine. It is exported for testing purposes.
 func (impl Implementation) Dlasy2(tranl, tranr bool, isgn, n1, n2 int, tl []float64, ldtl int, tr []float64, ldtr int, b []float64, ldb int, x []float64, ldx int) (scale, xnorm float64, ok bool) {
diff --git a/native/lapack_test.go b/native/lapack_test.go
index d9506cf8..7525eccb 100644
--- a/native/lapack_test.go
+++ b/native/lapack_test.go
@@ -108,6 +108,10 @@ func TestDlas2(t *testing.T) {
 	testlapack.Dlas2Test(t, impl)
 }
 
+func TestDlasy2(t *testing.T) {
+	testlapack.Dlasy2Test(t, impl)
+}
+
 func TestDlanst(t *testing.T) {
 	testlapack.DlanstTest(t, impl)
 }
diff --git a/testlapack/dlasy2.go b/testlapack/dlasy2.go
new file mode 100644
index 00000000..44ea65cb
--- /dev/null
+++ b/testlapack/dlasy2.go
@@ -0,0 +1,104 @@
+// 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"
+	"math/rand"
+	"testing"
+
+	"github.com/gonum/blas"
+	"github.com/gonum/blas/blas64"
+)
+
+type Dlasy2er interface {
+	Dlasy2(tranl, tranr bool, isgn, n1, n2 int, tl []float64, ldtl int, tr []float64, ldtr int, b []float64, ldb int, x []float64, ldx int) (scale, xnorm float64, ok bool)
+}
+
+func Dlasy2Test(t *testing.T, impl Dlasy2er) {
+	rnd := rand.New(rand.NewSource(1))
+	for _, tranl := range []bool{true, false} {
+		for _, tranr := range []bool{true, false} {
+			for _, isgn := range []int{1, -1} {
+				for _, n1 := range []int{0, 1, 2} {
+					for _, n2 := range []int{0, 1, 2} {
+						for _, extra := range []int{0, 1, 2, 13} {
+							for cas := 0; cas < 1000; cas++ {
+								testDlasy2(t, impl, tranl, tranr, isgn, n1, n2, extra, rnd)
+							}
+						}
+					}
+				}
+			}
+		}
+	}
+}
+
+func testDlasy2(t *testing.T, impl Dlasy2er, tranl, tranr bool, isgn, n1, n2, extra int, rnd *rand.Rand) {
+	const tol = 1e-11
+
+	tl := randomGeneral(n1, n1, n1+extra, rnd)
+	tr := randomGeneral(n2, n2, n2+extra, rnd)
+	b := randomGeneral(n1, n2, n2+extra, rnd)
+	x := randomGeneral(n1, n2, n2+extra, rnd)
+
+	scale, xnorm, ok := impl.Dlasy2(tranl, tranr, isgn, n1, n2, tl.Data, tl.Stride, tr.Data, tr.Stride, b.Data, b.Stride, x.Data, x.Stride)
+	if scale > 1 {
+		t.Errorf("invalid value of scale, want <= 1, got %v", scale)
+	}
+	if n1 == 0 || n2 == 0 {
+		return
+	}
+
+	prefix := fmt.Sprintf("Case n1=%v, n2=%v, isgn=%v", n1, n2, isgn)
+
+	// Check any invalid modifications of x.
+	if !generalOutsideAllNaN(x) {
+		t.Errorf("%v: out-of-range write to x\n%v", prefix, x.Data)
+	}
+
+	var xnormWant float64
+	for i := 0; i < n1; i++ {
+		var rowsum float64
+		for j := 0; j < n2; j++ {
+			rowsum += math.Abs(x.Data[i*x.Stride+j])
+		}
+		if rowsum > xnormWant {
+			xnormWant = rowsum
+		}
+	}
+	if xnormWant != xnorm {
+		t.Errorf("%v: unexpected xnorm: want %v, got %v", prefix, xnormWant, xnorm)
+	}
+
+	// Multiply b by scale to get the wanted right-hand side.
+	for i := 0; i < n1; i++ {
+		for j := 0; j < n2; j++ {
+			b.Data[i*b.Stride+j] *= scale
+		}
+	}
+	// Compute the wanted left-hand side.
+	lhsWant := randomGeneral(n1, n2, n2, rnd)
+	if tranl {
+		blas64.Gemm(blas.Trans, blas.NoTrans, 1, tl, x, 0, lhsWant)
+	} else {
+		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, tl, x, 0, lhsWant)
+	}
+	if tranr {
+		blas64.Gemm(blas.NoTrans, blas.Trans, float64(isgn), x, tr, 1, lhsWant)
+	} else {
+		blas64.Gemm(blas.NoTrans, blas.NoTrans, float64(isgn), x, tr, 1, lhsWant)
+	}
+	// Compare them.
+	for i := 0; i < n1; i++ {
+		for j := 0; j < n2; j++ {
+			diff := lhsWant.Data[i*lhsWant.Stride+j] - b.Data[i*b.Stride+j]
+			if math.Abs(diff) > tol && ok {
+				t.Errorf("%v: unexpected result, diff[%v,%v]=%v", prefix, i, j, diff)
+			}
+		}
+	}
+}