native: add test for Dlaqr1

native/lapack_test.go

@@ -108,6 +108,10 @@ func TestDlantr(t *testing.T) {
 	testlapack.DlantrTest(t, impl)
 }
 
+func TestDlaqr1(t *testing.T) {
+	testlapack.Dlaqr1Test(t, impl)
+}
+
 func TestDlarfb(t *testing.T) {
 	testlapack.DlarfbTest(t, impl)
 }

testlapack/dlaqr1.go

@@ -0,0 +1,195 @@
+// 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 (
+	"math"
+	"math/rand"
+	"testing"
+)
+
+type Dlaqr1er interface {
+	Dlaqr1(n int, h []float64, ldh int, sr1, si1, sr2, si2 float64, v []float64)
+}
+
+func Dlaqr1Test(t *testing.T, impl Dlaqr1er) {
+	rnd := rand.New(rand.NewSource(1))
+
+	n := 2
+	for k := 0; k < 100; k++ {
+		for _, ldh := range []int{2, 3, 7} {
+			v := make([]float64, n)
+			h := make([]float64, n*(n-1)*ldh)
+
+			// Case sr1 = sr2, si1 = -si2.
+
+			for i := range v {
+				v[i] = math.NaN()
+			}
+			for i := range h {
+				h[i] = math.NaN()
+			}
+			for i := 0; i < 2; i++ {
+				for j := 0; j < 2; j++ {
+					h[i*ldh+j] = rnd.NormFloat64()
+				}
+			}
+
+			sr1 := rnd.NormFloat64()
+			si := rnd.NormFloat64()
+			impl.Dlaqr1(n, h, ldh, sr1, si, sr1, -si, v)
+
+			// Matrix H - s1*I.
+			h1a := complex(h[0], 0) - complex(sr1, si)
+			h1b := complex(h[1], 0)
+			h1c := complex(h[ldh], 0)
+			h1d := complex(h[ldh+1], 0) - complex(sr1, si)
+			// First column of H - s2*I.
+			h2a := complex(h[0], 0) - complex(sr1, -si)
+			h2c := complex(h[ldh], 0)
+
+			// Multiply (H-s1*I)*(H-s2*I) to get a tentative v[0].
+			wantv0 := h1a*h2a + h1b*h2c
+			// Get the unknown scale.
+			scale := v[0] / real(wantv0)
+
+			// Compute v[1].
+			wantv1 := scale * real(h1c*h2a+h1d*h2c)
+			// The scale must be the same for both elements.
+			if math.Abs(wantv1-v[1]) > 1e-12 {
+				t.Errorf("n = %v, ldh = %v, si = -si: Unexpected value of v[1]: got %v, want %v", n, ldh, v[1], wantv1)
+			}
+
+			// Case sr1 != sr2, si1 = si2 = 0.
+
+			for i := range v {
+				v[i] = math.NaN()
+			}
+			for i := range h {
+				h[i] = math.NaN()
+			}
+			for i := 0; i < 2; i++ {
+				for j := 0; j < 2; j++ {
+					h[i*ldh+j] = rnd.NormFloat64()
+				}
+			}
+
+			sr2 := rnd.NormFloat64()
+			impl.Dlaqr1(n, h, ldh, sr1, 0, sr2, 0, v)
+
+			h1a = complex(h[0], 0) - complex(sr1, 0)
+			h1b = complex(h[1], 0)
+			h1c = complex(h[ldh], 0)
+			h1d = complex(h[ldh+1], 0) - complex(sr1, 0)
+			h2a = complex(h[0], 0) - complex(sr2, 0)
+			h2c = complex(h[ldh], 0)
+
+			wantv0 = h1a*h2a + h1b*h2c
+			scale = v[0] / real(wantv0)
+
+			wantv1 = scale * real(h1c*h2a+h1d*h2c)
+			if math.Abs(wantv1-v[1]) > 1e-12 {
+				t.Errorf("n = %v, ldh = %v, si = 0: Unexpected value of v[1]: got %v, want %v", n, ldh, v[1], wantv1)
+			}
+		}
+	}
+
+	n = 3
+	for k := 0; k < 100; k++ {
+		for _, ldh := range []int{3, 4, 7} {
+			v := make([]float64, n)
+			h := make([]float64, n*(n-1)*ldh)
+
+			// Case sr1 = sr2, si1 = -si2.
+
+			for i := range v {
+				v[i] = math.NaN()
+			}
+			for i := range h {
+				h[i] = math.NaN()
+			}
+			for i := 0; i < 3; i++ {
+				for j := 0; j < 3; j++ {
+					h[i*ldh+j] = rnd.NormFloat64()
+				}
+			}
+
+			sr1 := rnd.NormFloat64()
+			si := rnd.NormFloat64()
+			impl.Dlaqr1(n, h, ldh, sr1, si, sr1, -si, v)
+
+			h1a := complex(h[0], 0) - complex(sr1, si)
+			h1b := complex(h[1], 0)
+			h1c := complex(h[2], 0)
+			h1d := complex(h[ldh], 0)
+			h1e := complex(h[ldh+1], 0) - complex(sr1, si)
+			h1f := complex(h[ldh+2], 0)
+			h1g := complex(h[2*ldh], 0)
+			h1h := complex(h[2*ldh+1], 0)
+			h1i := complex(h[2*ldh+2], 0) - complex(sr1, si)
+
+			h2a := complex(h[0], 0) - complex(sr1, -si)
+			h2d := complex(h[ldh], 0)
+			h2g := complex(h[2*ldh], 0)
+
+			wantv0 := h1a*h2a + h1b*h2d + h1c*h2g
+			scale := v[0] / real(wantv0)
+
+			wantv1 := scale * real(h1d*h2a+h1e*h2d+h1f*h2g)
+			if math.Abs(wantv1-v[1]) > 1e-12 {
+				t.Errorf("n = %v, ldh = %v, si = -si: Unexpected value of v[1]: got %v, want %v", n, ldh, v[1], wantv1)
+			}
+
+			wantv2 := scale * real(h1g*h2a+h1h*h2d+h1i*h2g)
+			if math.Abs(wantv2-v[2]) > 1e-12 {
+				t.Errorf("n = %v, ldh = %v, si = -si: Unexpected value of v[2]: got %v, want %v", n, ldh, v[2], wantv2)
+			}
+
+			// Case sr1 != sr2, si1 = si2 = 0.
+
+			for i := range h {
+				h[i] = math.NaN()
+			}
+			for i := range v {
+				v[i] = math.NaN()
+			}
+			for i := 0; i < 3; i++ {
+				for j := 0; j < 3; j++ {
+					h[i*ldh+j] = rnd.NormFloat64()
+				}
+			}
+
+			sr2 := rnd.NormFloat64()
+			impl.Dlaqr1(n, h, ldh, sr1, 0, sr2, 0, v)
+
+			h1a = complex(h[0], 0) - complex(sr1, 0)
+			h1b = complex(h[1], 0)
+			h1c = complex(h[2], 0)
+			h1d = complex(h[ldh], 0)
+			h1e = complex(h[ldh+1], 0) - complex(sr1, 0)
+			h1f = complex(h[ldh+2], 0)
+			h1g = complex(h[2*ldh], 0)
+			h1h = complex(h[2*ldh+1], 0)
+			h1i = complex(h[2*ldh+2], 0) - complex(sr1, 0)
+
+			h2a = complex(h[0], 0) - complex(sr2, 0)
+			h2d = complex(h[ldh], 0)
+			h2g = complex(h[2*ldh], 0)
+
+			wantv0 = h1a*h2a + h1b*h2d + h1c*h2g
+			scale = v[0] / real(wantv0)
+
+			wantv1 = scale * real(h1d*h2a+h1e*h2d+h1f*h2g)
+			if math.Abs(wantv1-v[1]) > 1e-12 {
+				t.Errorf("n = %v, ldh = %v, si = 0: Unexpected value of v[1]: got %v, want %v", n, ldh, v[1], wantv1)
+			}
+
+			wantv2 = scale * real(h1g*h2a+h1h*h2d+h1i*h2g)
+			if math.Abs(wantv2-v[2]) > 1e-12 {
+				t.Errorf("n = %v, ldh = %v, si = 0: Unexpected value of v[2]: got %v, want %v", n, ldh, v[2], wantv2)
+			}
+		}
+	}
+}