native,cgo,testlapack: add test for Dgeev

cgo/lapack_test.go

@@ -80,6 +80,10 @@ func TestDgecon(t *testing.T) {
 	testlapack.DgeconTest(t, impl)
 }
 
+func TestDgeev(t *testing.T) {
+	testlapack.DgeevTest(t, impl)
+}
+
 func TestDgehrd(t *testing.T) {
 	testlapack.DgehrdTest(t, impl)
 }

native/lapack_test.go

@@ -40,6 +40,10 @@ func TestDgecon(t *testing.T) {
 	testlapack.DgeconTest(t, impl)
 }
 
+func TestDgeev(t *testing.T) {
+	testlapack.DgeevTest(t, impl)
+}
+
 func TestDgehd2(t *testing.T) {
 	testlapack.Dgehd2Test(t, impl)
 }

testlapack/dgeev.go

@@ -0,0 +1,317 @@
+// 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/cmplx"
+	"math/rand"
+	"testing"
+
+	"github.com/gonum/blas/blas64"
+	"github.com/gonum/floats"
+	"github.com/gonum/lapack"
+)
+
+type Dgeever interface {
+	Dgeev(jobvl lapack.JobLeftEV, jobvr lapack.JobRightEV, n int, a []float64, lda int,
+		wr, wi []float64, vl []float64, ldvl int, vr []float64, ldvr int, work []float64, lwork int) int
+}
+
+func DgeevTest(t *testing.T, impl Dgeever) {
+	type evTestMatrix interface {
+		Matrix() blas64.General
+		Eigenvalues() []complex128
+	}
+
+	rnd := rand.New(rand.NewSource(1))
+
+	for i, test := range []evTestMatrix{
+		A123{},
+
+		NewAntisymRandom(10, rnd),
+		NewAntisymRandom(51, rnd),
+		NewAntisymRandom(100, rnd),
+
+		Circulant(2),
+		Circulant(3),
+		Circulant(4),
+		Circulant(5),
+		Circulant(10),
+		Circulant(15),
+		Circulant(30),
+		Circulant(50),
+		Circulant(101),
+		Circulant(250),
+
+		Clement(1),
+		Clement(2),
+		Clement(3),
+		Clement(4),
+		Clement(5),
+		Clement(10),
+		Clement(15),
+		Clement(30),
+
+		Creation(2),
+		Creation(3),
+		Creation(4),
+		Creation(5),
+		Creation(10),
+		Creation(15),
+		Creation(30),
+		Creation(101),
+
+		Diagonal(0),
+		Diagonal(1),
+		Diagonal(2),
+		Diagonal(3),
+		Diagonal(4),
+		Diagonal(5),
+		Diagonal(10),
+		Diagonal(15),
+		Diagonal(30),
+
+		Downshift(2),
+		Downshift(3),
+		Downshift(4),
+		Downshift(5),
+		Downshift(10),
+		Downshift(15),
+		Downshift(30),
+
+		Fibonacci(1),
+		Fibonacci(2),
+		Fibonacci(3),
+		Fibonacci(4),
+		Fibonacci(5),
+		Fibonacci(10),
+		Fibonacci(15),
+		Fibonacci(30),
+		Fibonacci(101),
+
+		Gear(1),
+		Gear(2),
+		Gear(3),
+		Gear(5),
+		Gear(10),
+		Gear(15),
+		Gear(50),
+		Gear(101),
+
+		Grcar{N: 10, K: 3},
+		Grcar{N: 10, K: 7},
+		Grcar{N: 50, K: 3},
+		Grcar{N: 50, K: 10},
+		Grcar{N: 50, K: 30},
+		Grcar{N: 150, K: 2},
+		Grcar{N: 150, K: 148},
+
+		Hanowa{N: 6, Alpha: 17},
+		Hanowa{N: 50, Alpha: -1},
+		Hanowa{N: 100, Alpha: -1},
+
+		Lesp(2),
+		Lesp(3),
+		Lesp(4),
+		Lesp(5),
+		Lesp(6),
+		Lesp(10),
+		Lesp(15),
+		Lesp(50),
+		Lesp(101),
+
+		Rutis{},
+
+		Tris{N: 74, X: 1, Y: -2, Z: 1},
+		Tris{N: 74, X: 1, Y: 2, Z: -3},
+		Tris{N: 75, X: 1, Y: 2, Z: -3},
+
+		Wilk4{},
+		Wilk12{},
+		Wilk20(0),
+		Wilk20(1e-10),
+
+		Zero(10),
+		Zero(50),
+		Zero(100),
+	} {
+		for _, jobvl := range []lapack.JobLeftEV{lapack.ComputeLeftEV, lapack.None} {
+			for _, jobvr := range []lapack.JobRightEV{lapack.ComputeRightEV, lapack.None} {
+				for _, extra := range []int{0, 11} {
+					for _, optwork := range []bool{false, true} {
+						testDgeev(t, impl, i, test.Matrix(), test.Eigenvalues(), jobvl, jobvr, extra, optwork)
+					}
+				}
+			}
+		}
+	}
+}
+
+func testDgeev(t *testing.T, impl Dgeever, itest int, a blas64.General, wantEigVal []complex128, jobvl lapack.JobLeftEV, jobvr lapack.JobRightEV, extra int, optwork bool) {
+	const tol = 1e-10
+
+	n := a.Rows
+	aCopy := cloneGeneral(a)
+
+	var vl blas64.General
+	if jobvl == lapack.ComputeLeftEV {
+		vl = nanGeneral(n, n, n)
+	}
+
+	var vr blas64.General
+	if jobvr == lapack.ComputeRightEV {
+		vr = nanGeneral(n, n, n)
+	}
+
+	wr := make([]float64, n)
+	wi := make([]float64, n)
+
+	lwork := max(1, 3*n)
+	if jobvl == lapack.ComputeLeftEV || jobvr == lapack.ComputeRightEV {
+		lwork = max(1, 4*n)
+	}
+	if optwork {
+		work := make([]float64, 1)
+		impl.Dgeev(jobvl, jobvr, n, nil, 1, nil, nil, nil, 1, nil, 1, work, -1)
+		lwork = int(work[0])
+	}
+	work := make([]float64, lwork)
+
+	firstok := impl.Dgeev(jobvl, jobvr, n, a.Data, a.Stride, wr, wi,
+		vl.Data, vl.Stride, vr.Data, vr.Stride, work, len(work))
+
+	prefix := fmt.Sprintf("Case #%v: n=%v, jobvl=%v, jobvr=%v, extra=%v, optwk=%v",
+		itest, n, jobvl, jobvr, extra, optwork)
+
+	if !generalOutsideAllNaN(vl) {
+		t.Errorf("%v: out-of-range write to VL", prefix)
+	}
+	if !generalOutsideAllNaN(vr) {
+		t.Errorf("%v: out-of-range write to VR", prefix)
+	}
+
+	if firstok > 0 {
+		t.Log("%v: all eigenvalues haven't been computed, firstok=%v", prefix, firstok)
+	}
+
+	// Check that conjugate pair eigevalues are ordered correctly.
+	for i := firstok; i < n; {
+		if wi[i] == 0 {
+			i++
+			continue
+		}
+		if wr[i] != wr[i+1] {
+			t.Errorf("%v: real parts of %vth conjugate pair not equal", prefix, i)
+		}
+		if wi[i] < 0 || wi[i+1] > 0 {
+			t.Errorf("%v: unexpected ordering of %vth conjugate pair", prefix, i)
+		}
+		i += 2
+	}
+
+	// If the eigenvectors were not computed, check the computed eigenvalues
+	// against provided known eigenvalues.
+	if firstok > 0 || (jobvl == lapack.None && jobvr == lapack.None) {
+		if wantEigVal == nil {
+			return
+		}
+		used := make([]bool, n)
+		for i := firstok; i < n; i++ {
+			gotev := complex(wr[i], wi[i])
+			idx := -1
+			for k, wantev := range wantEigVal {
+				if used[k] {
+					continue
+				}
+				if cmplx.Abs(wantev-gotev) < 1e-8 {
+					idx = k
+					used[k] = true
+					break
+				}
+			}
+			if idx == -1 {
+				t.Errorf("%v: unexpected eigenvalue %v", prefix, gotev)
+			}
+		}
+		return
+	}
+
+	for k := 0; k < n; {
+		if wi[k] == 0 {
+			if jobvl == lapack.ComputeLeftEV {
+				ev := columnOf(vl, k)
+				if !isLeftEigenvectorOf(aCopy, ev, nil, complex(wr[k], 0), tol) {
+					t.Errorf("%v: VL[:,%v] is not left real eigenvector",
+						prefix, k)
+				}
+
+				norm := floats.Norm(ev, 2)
+				if math.Abs(norm-1) >= tol {
+					t.Errorf("%v: norm of left real eigenvector %v not equal to 1: got %v",
+						prefix, k, norm)
+				}
+			}
+			if jobvr == lapack.ComputeRightEV {
+				ev := columnOf(vr, k)
+				if !isRightEigenvectorOf(aCopy, ev, nil, complex(wr[k], 0), tol) {
+					t.Errorf("%v: VR[:,%v] is not right real eigenvector",
+						prefix, k)
+				}
+
+				norm := floats.Norm(ev, 2)
+				if math.Abs(norm-1) >= tol {
+					t.Errorf("%v: norm of right real eigenvector %v not equal to 1: got %v",
+						prefix, k, norm)
+				}
+			}
+			k++
+		} else {
+			if jobvl == lapack.ComputeLeftEV {
+				evre := columnOf(vl, k)
+				evim := columnOf(vl, k+1)
+				if !isLeftEigenvectorOf(aCopy, evre, evim, complex(wr[k], wi[k]), tol) {
+					t.Errorf("%v: VL[:,%v:%v] is not left complex eigenvector",
+						prefix, k, k+1)
+				}
+				floats.Scale(-1, evim)
+				if !isLeftEigenvectorOf(aCopy, evre, evim, complex(wr[k+1], wi[k+1]), tol) {
+					t.Errorf("%v: VL[:,%v:%v] is not left complex eigenvector",
+						prefix, k, k+1)
+				}
+
+				norm := math.Hypot(floats.Norm(evre, 2), floats.Norm(evim, 2))
+				if math.Abs(norm-1) > tol {
+					t.Errorf("%v: norm of left complex eigenvector %v not equal to 1: got %v",
+						prefix, k, norm)
+				}
+			}
+			if jobvr == lapack.ComputeRightEV {
+				evre := columnOf(vr, k)
+				evim := columnOf(vr, k+1)
+				if !isRightEigenvectorOf(aCopy, evre, evim, complex(wr[k], wi[k]), tol) {
+					t.Errorf("%v: VR[:,%v:%v] is not right complex eigenvector",
+						prefix, k, k+1)
+				}
+				floats.Scale(-1, evim)
+				if !isRightEigenvectorOf(aCopy, evre, evim, complex(wr[k+1], wi[k+1]), tol) {
+					t.Errorf("%v: VR[:,%v:%v] is not right complex eigenvector",
+						prefix, k, k+1)
+				}
+
+				norm := math.Hypot(floats.Norm(evre, 2), floats.Norm(evim, 2))
+				if math.Abs(norm-1) > tol {
+					t.Errorf("%v: norm of right complex eigenvector %v not equal to 1: got %v",
+						prefix, k, norm)
+				}
+			}
+			// We don't test whether the largest component is real
+			// because checking it is flaky due to rounding errors.
+
+			k += 2
+		}
+	}
+}