testlapack: add DlahqrTest

native/lapack_test.go

@@ -96,6 +96,10 @@ func TestDlaev2(t *testing.T) {
 	testlapack.Dlaev2Test(t, impl)
 }
 
+func TestDlahqr(t *testing.T) {
+	testlapack.DlahqrTest(t, impl)
+}
+
 func TestDlahr2(t *testing.T) {
 	testlapack.Dlahr2Test(t, impl)
 }

testlapack/dlahqr.go

@@ -0,0 +1,159 @@
+// 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/blas64"
+)
+
+type Dlahqrer interface {
+	Dlahqr(wantt, wantz bool, n, ilo, ihi int, h []float64, ldh int, wr, wi []float64, iloz, ihiz int, z []float64, ldz int) int
+}
+
+func DlahqrTest(t *testing.T, impl Dlahqrer) {
+	rnd := rand.New(rand.NewSource(1))
+	for _, wantt := range []bool{true, false} {
+		for _, wantz := range []bool{true, false} {
+			for _, n := range []int{1, 2, 3, 4, 5, 6, 10, 18, 31, 53} {
+				for _, extra := range []int{0, 1, 11} {
+					for cas := 0; cas < 10; cas++ {
+						ilo := rnd.Intn(n)
+						ihi := rnd.Intn(n)
+						if ilo > ihi {
+							ilo, ihi = ihi, ilo
+						}
+						iloz := rnd.Intn(ilo + 1)
+						ihiz := ihi + rnd.Intn(n-ihi)
+						testDlahqr(t, impl, wantt, wantz, n, ilo, ihi, iloz, ihiz, extra, rnd)
+					}
+				}
+			}
+		}
+	}
+	for _, wantt := range []bool{true, false} {
+		for _, wantz := range []bool{true, false} {
+			testDlahqr(t, impl, wantt, wantz, 0, 0, -1, 0, -1, 0, rnd)
+		}
+	}
+}
+
+func testDlahqr(t *testing.T, impl Dlahqrer, wantt, wantz bool, n, ilo, ihi, iloz, ihiz, extra int, rnd *rand.Rand) {
+	const tol = 1e-14
+
+	var z, zCopy blas64.General
+	if wantz {
+		z = eye(n, n+extra)
+		zCopy = cloneGeneral(z)
+	}
+	h := randomHessenberg(n, n+extra, rnd)
+	if ilo > 0 {
+		h.Data[ilo*h.Stride+ilo-1] = 0
+	}
+	wr := nanSlice(n)
+	wi := nanSlice(n)
+
+	info := impl.Dlahqr(wantt, wantz, n, ilo, ihi, h.Data, h.Stride, wr, wi, iloz, ihiz, z.Data, z.Stride)
+
+	prefix := fmt.Sprintf("Case n=%v, ilo=%v, ihi=%v, iloz=%v, ihiz=%v, wantt=%v, wantz=%v, extra=%v", n, ilo, ihi, iloz, ihiz, wantt, wantz, extra)
+
+	if !generalOutsideAllNaN(h) {
+		t.Errorf("%v: out-of-range write to H\n%v", prefix, h.Data)
+	}
+	if !generalOutsideAllNaN(z) {
+		t.Errorf("%v: out-of-range write to Z\n%v", prefix, z.Data)
+	}
+
+	if wantz {
+		// Z should contain the orthogonal matrix U.
+		if !isOrthonormal(z) {
+			t.Errorf("%v: Z is not orthogonal", prefix)
+		}
+		// Z should have been modified only in the
+		// [iloz:ihiz+1:ilo:ihi+1] block.
+		for i := 0; i < n; i++ {
+			for j := 0; j < n; j++ {
+				if iloz <= i && i <= ihiz && ilo <= j && j <= ihi {
+					continue
+				}
+				if z.Data[i*z.Stride+j] != zCopy.Data[i*zCopy.Stride+j] {
+					t.Errorf("%v: Z modified outside of [iloz:ihiz+1,ilo:ihi+1] block", prefix)
+				}
+			}
+		}
+	}
+
+	start := ilo
+	if info > 0 {
+		start = info + 1
+	}
+	// Check that wr and wi have not been modified outside [start:ihi+1].
+	for _, v := range wr[:start] {
+		if !math.IsNaN(v) {
+			t.Errorf("%v: wr modified before [ilo:ihiz+1] block", prefix)
+		}
+	}
+	for _, v := range wr[ihi+1:] {
+		if !math.IsNaN(v) {
+			t.Errorf("%v: wr modified after [ilo:ihiz+1] block", prefix)
+		}
+	}
+	for _, v := range wi[:start] {
+		if !math.IsNaN(v) {
+			t.Errorf("%v: wi modified before [ilo:ihiz+1] block", prefix)
+		}
+	}
+	for _, v := range wi[ihi+1:] {
+		if !math.IsNaN(v) {
+			t.Errorf("%v: wi modified after [ilo:ihiz+1] block", prefix)
+		}
+	}
+	var nReal int
+	for i := start; i <= ihi; {
+		if wi[i] == 0 {
+			// A real eigenvalue.
+			nReal++
+			if wantt && wr[i] != h.Data[i*h.Stride+i] {
+				t.Errorf("%v: wr[%v] != H[%v,%v]", prefix, i, i, i)
+			}
+			if i > 0 && h.Data[i*h.Stride+i-1] != 0 {
+				t.Errorf("%v: H[%v,%v] != 0", prefix, i, i-1)
+			}
+			if i < ihi && h.Data[(i+1)*h.Stride+i] != 0 {
+				t.Errorf("%v: H[%v,%v] != 0", prefix, i+1, i)
+			}
+			i++
+			continue
+		}
+		// Complex eigenvalue.
+		if wr[i] != wr[i+1] {
+			t.Errorf("%v: real part of conjugate pair not equal, i=%v", prefix, i)
+		}
+		if wi[i] < 0 {
+			t.Errorf("%v: wi[%v] not positive", prefix, i)
+		}
+		if wi[i] != -wi[i+1] {
+			t.Errorf("%v: wi[%v] != wi[%v]", prefix, i, i+1)
+		}
+		prod := math.Abs(h.Data[(i+1)*h.Stride+i] * h.Data[i*h.Stride+i+1])
+		if wantt && math.Abs(math.Sqrt(prod)-wi[0]) > tol {
+			t.Errorf("%v: unexpected value of wi[%v]: want %v, got %v", prefix, i, math.Sqrt(prod), wi[i])
+		}
+		if i > 0 && h.Data[i*h.Stride+i-1] != 0 {
+			t.Errorf("%v: H[%v,%v] != 0", prefix, i, i-1)
+		}
+		if i < ihi-1 && h.Data[(i+2)*h.Stride+i+1] != 0 {
+			t.Errorf("%v: H[%v,%v] != 0", prefix, i+2, i+1)
+		}
+		i += 2
+	}
+	if (ihi+1-start)%2 != 0 && nReal == 0 {
+		t.Errorf("%v: expected at least one real eigenvalue", prefix)
+	}
+}