dgeev.go

// 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"
	"strconv"
	"testing"

	"golang.org/x/exp/rand"

	"gonum.org/v1/gonum/blas"
	"gonum.org/v1/gonum/blas/blas64"
	"gonum.org/v1/gonum/floats"
	"gonum.org/v1/gonum/lapack"
)

type Dgeever interface {
	Dgeev(jobvl lapack.LeftEVJob, jobvr lapack.RightEVJob, n int, a []float64, lda int,
		wr, wi []float64, vl []float64, ldvl int, vr []float64, ldvr int, work []float64, lwork int) int
}

type dgeevTest struct {
	a      blas64.General
	evWant []complex128 // If nil, the eigenvalues are not known.
	valTol float64      // Tolerance for eigenvalue checks.
	vecTol float64      // Tolerance for eigenvector checks.
}

func DgeevTest(t *testing.T, impl Dgeever) {
	rnd := rand.New(rand.NewSource(1))

	for i, test := range []dgeevTest{
		{
			a:      A123{}.Matrix(),
			evWant: A123{}.Eigenvalues(),
		},

		dgeevTestForAntisymRandom(10, rnd),
		dgeevTestForAntisymRandom(11, rnd),
		dgeevTestForAntisymRandom(50, rnd),
		dgeevTestForAntisymRandom(51, rnd),
		dgeevTestForAntisymRandom(100, rnd),
		dgeevTestForAntisymRandom(101, rnd),

		{
			a:      Circulant(2).Matrix(),
			evWant: Circulant(2).Eigenvalues(),
		},
		{
			a:      Circulant(3).Matrix(),
			evWant: Circulant(3).Eigenvalues(),
		},
		{
			a:      Circulant(4).Matrix(),
			evWant: Circulant(4).Eigenvalues(),
		},
		{
			a:      Circulant(5).Matrix(),
			evWant: Circulant(5).Eigenvalues(),
		},
		{
			a:      Circulant(10).Matrix(),
			evWant: Circulant(10).Eigenvalues(),
		},
		{
			a:      Circulant(15).Matrix(),
			evWant: Circulant(15).Eigenvalues(),
			valTol: 1e-12,
		},
		{
			a:      Circulant(30).Matrix(),
			evWant: Circulant(30).Eigenvalues(),
			valTol: 1e-11,
			vecTol: 1e-12,
		},
		{
			a:      Circulant(50).Matrix(),
			evWant: Circulant(50).Eigenvalues(),
			valTol: 1e-11,
			vecTol: 1e-12,
		},
		{
			a:      Circulant(101).Matrix(),
			evWant: Circulant(101).Eigenvalues(),
			valTol: 1e-10,
			vecTol: 1e-11,
		},
		{
			a:      Circulant(150).Matrix(),
			evWant: Circulant(150).Eigenvalues(),
			valTol: 1e-9,
			vecTol: 1e-10,
		},

		{
			a:      Clement(2).Matrix(),
			evWant: Clement(2).Eigenvalues(),
		},
		{
			a:      Clement(3).Matrix(),
			evWant: Clement(3).Eigenvalues(),
		},
		{
			a:      Clement(4).Matrix(),
			evWant: Clement(4).Eigenvalues(),
		},
		{
			a:      Clement(5).Matrix(),
			evWant: Clement(5).Eigenvalues(),
		},
		{
			a:      Clement(10).Matrix(),
			evWant: Clement(10).Eigenvalues(),
		},
		{
			a:      Clement(15).Matrix(),
			evWant: Clement(15).Eigenvalues(),
		},
		{
			a:      Clement(30).Matrix(),
			evWant: Clement(30).Eigenvalues(),
			valTol: 1e-11,
		},
		{
			a:      Clement(50).Matrix(),
			evWant: Clement(50).Eigenvalues(),
			valTol: 1e-7,
			vecTol: 1e-11,
		},

		{
			a:      Creation(2).Matrix(),
			evWant: Creation(2).Eigenvalues(),
		},
		{
			a:      Creation(3).Matrix(),
			evWant: Creation(3).Eigenvalues(),
		},
		{
			a:      Creation(4).Matrix(),
			evWant: Creation(4).Eigenvalues(),
		},
		{
			a:      Creation(5).Matrix(),
			evWant: Creation(5).Eigenvalues(),
		},
		{
			a:      Creation(10).Matrix(),
			evWant: Creation(10).Eigenvalues(),
		},
		{
			a:      Creation(15).Matrix(),
			evWant: Creation(15).Eigenvalues(),
		},
		{
			a:      Creation(30).Matrix(),
			evWant: Creation(30).Eigenvalues(),
		},
		{
			a:      Creation(50).Matrix(),
			evWant: Creation(50).Eigenvalues(),
		},
		{
			a:      Creation(101).Matrix(),
			evWant: Creation(101).Eigenvalues(),
		},
		{
			a:      Creation(150).Matrix(),
			evWant: Creation(150).Eigenvalues(),
		},

		{
			a:      Diagonal(0).Matrix(),
			evWant: Diagonal(0).Eigenvalues(),
		},
		{
			a:      Diagonal(10).Matrix(),
			evWant: Diagonal(10).Eigenvalues(),
		},
		{
			a:      Diagonal(50).Matrix(),
			evWant: Diagonal(50).Eigenvalues(),
		},
		{
			a:      Diagonal(151).Matrix(),
			evWant: Diagonal(151).Eigenvalues(),
		},

		{
			a:      Downshift(2).Matrix(),
			evWant: Downshift(2).Eigenvalues(),
		},
		{
			a:      Downshift(3).Matrix(),
			evWant: Downshift(3).Eigenvalues(),
		},
		{
			a:      Downshift(4).Matrix(),
			evWant: Downshift(4).Eigenvalues(),
		},
		{
			a:      Downshift(5).Matrix(),
			evWant: Downshift(5).Eigenvalues(),
		},
		{
			a:      Downshift(10).Matrix(),
			evWant: Downshift(10).Eigenvalues(),
		},
		{
			a:      Downshift(15).Matrix(),
			evWant: Downshift(15).Eigenvalues(),
		},
		{
			a:      Downshift(30).Matrix(),
			evWant: Downshift(30).Eigenvalues(),
		},
		{
			a:      Downshift(50).Matrix(),
			evWant: Downshift(50).Eigenvalues(),
		},
		{
			a:      Downshift(101).Matrix(),
			evWant: Downshift(101).Eigenvalues(),
		},
		{
			a:      Downshift(150).Matrix(),
			evWant: Downshift(150).Eigenvalues(),
		},

		{
			a:      Fibonacci(2).Matrix(),
			evWant: Fibonacci(2).Eigenvalues(),
		},
		{
			a:      Fibonacci(3).Matrix(),
			evWant: Fibonacci(3).Eigenvalues(),
		},
		{
			a:      Fibonacci(4).Matrix(),
			evWant: Fibonacci(4).Eigenvalues(),
		},
		{
			a:      Fibonacci(5).Matrix(),
			evWant: Fibonacci(5).Eigenvalues(),
		},
		{
			a:      Fibonacci(10).Matrix(),
			evWant: Fibonacci(10).Eigenvalues(),
		},
		{
			a:      Fibonacci(15).Matrix(),
			evWant: Fibonacci(15).Eigenvalues(),
		},
		{
			a:      Fibonacci(30).Matrix(),
			evWant: Fibonacci(30).Eigenvalues(),
		},
		{
			a:      Fibonacci(50).Matrix(),
			evWant: Fibonacci(50).Eigenvalues(),
		},
		{
			a:      Fibonacci(101).Matrix(),
			evWant: Fibonacci(101).Eigenvalues(),
		},
		{
			a:      Fibonacci(150).Matrix(),
			evWant: Fibonacci(150).Eigenvalues(),
		},

		{
			a:      Gear(2).Matrix(),
			evWant: Gear(2).Eigenvalues(),
		},
		{
			a:      Gear(3).Matrix(),
			evWant: Gear(3).Eigenvalues(),
		},
		{
			a:      Gear(4).Matrix(),
			evWant: Gear(4).Eigenvalues(),
			valTol: 1e-7,
		},
		{
			a:      Gear(5).Matrix(),
			evWant: Gear(5).Eigenvalues(),
		},
		{
			a:      Gear(10).Matrix(),
			evWant: Gear(10).Eigenvalues(),
			valTol: 1e-8,
		},
		{
			a:      Gear(15).Matrix(),
			evWant: Gear(15).Eigenvalues(),
		},
		{
			a:      Gear(30).Matrix(),
			evWant: Gear(30).Eigenvalues(),
			valTol: 1e-8,
		},
		{
			a:      Gear(50).Matrix(),
			evWant: Gear(50).Eigenvalues(),
			valTol: 1e-8,
		},
		{
			a:      Gear(101).Matrix(),
			evWant: Gear(101).Eigenvalues(),
		},
		{
			a:      Gear(150).Matrix(),
			evWant: Gear(150).Eigenvalues(),
			valTol: 1e-8,
		},

		{
			a:      Grcar{N: 10, K: 3}.Matrix(),
			evWant: Grcar{N: 10, K: 3}.Eigenvalues(),
		},
		{
			a:      Grcar{N: 10, K: 7}.Matrix(),
			evWant: Grcar{N: 10, K: 7}.Eigenvalues(),
		},
		{
			a:      Grcar{N: 11, K: 7}.Matrix(),
			evWant: Grcar{N: 11, K: 7}.Eigenvalues(),
		},
		{
			a:      Grcar{N: 50, K: 3}.Matrix(),
			evWant: Grcar{N: 50, K: 3}.Eigenvalues(),
		},
		{
			a:      Grcar{N: 51, K: 3}.Matrix(),
			evWant: Grcar{N: 51, K: 3}.Eigenvalues(),
		},
		{
			a:      Grcar{N: 50, K: 10}.Matrix(),
			evWant: Grcar{N: 50, K: 10}.Eigenvalues(),
		},
		{
			a:      Grcar{N: 51, K: 10}.Matrix(),
			evWant: Grcar{N: 51, K: 10}.Eigenvalues(),
		},
		{
			a:      Grcar{N: 50, K: 30}.Matrix(),
			evWant: Grcar{N: 50, K: 30}.Eigenvalues(),
		},
		{
			a:      Grcar{N: 150, K: 2}.Matrix(),
			evWant: Grcar{N: 150, K: 2}.Eigenvalues(),
		},
		{
			a:      Grcar{N: 150, K: 148}.Matrix(),
			evWant: Grcar{N: 150, K: 148}.Eigenvalues(),
		},

		{
			a:      Hanowa{N: 6, Alpha: 17}.Matrix(),
			evWant: Hanowa{N: 6, Alpha: 17}.Eigenvalues(),
		},
		{
			a:      Hanowa{N: 50, Alpha: -1}.Matrix(),
			evWant: Hanowa{N: 50, Alpha: -1}.Eigenvalues(),
		},
		{
			a:      Hanowa{N: 100, Alpha: -1}.Matrix(),
			evWant: Hanowa{N: 100, Alpha: -1}.Eigenvalues(),
		},

		{
			a:      Lesp(2).Matrix(),
			evWant: Lesp(2).Eigenvalues(),
		},
		{
			a:      Lesp(3).Matrix(),
			evWant: Lesp(3).Eigenvalues(),
		},
		{
			a:      Lesp(4).Matrix(),
			evWant: Lesp(4).Eigenvalues(),
		},
		{
			a:      Lesp(5).Matrix(),
			evWant: Lesp(5).Eigenvalues(),
		},
		{
			a:      Lesp(10).Matrix(),
			evWant: Lesp(10).Eigenvalues(),
		},
		{
			a:      Lesp(15).Matrix(),
			evWant: Lesp(15).Eigenvalues(),
		},
		{
			a:      Lesp(30).Matrix(),
			evWant: Lesp(30).Eigenvalues(),
		},
		{
			a:      Lesp(50).Matrix(),
			evWant: Lesp(50).Eigenvalues(),
			valTol: 1e-12,
			vecTol: 1e-12,
		},
		{
			a:      Lesp(101).Matrix(),
			evWant: Lesp(101).Eigenvalues(),
			valTol: 1e-12,
			vecTol: 1e-12,
		},
		{
			a:      Lesp(150).Matrix(),
			evWant: Lesp(150).Eigenvalues(),
			valTol: 1e-12,
			vecTol: 1e-12,
		},

		{
			a:      Rutis{}.Matrix(),
			evWant: Rutis{}.Eigenvalues(),
		},

		{
			a:      Tris{N: 74, X: 1, Y: -2, Z: 1}.Matrix(),
			evWant: Tris{N: 74, X: 1, Y: -2, Z: 1}.Eigenvalues(),
		},
		{
			a:      Tris{N: 74, X: 1, Y: 2, Z: -3}.Matrix(),
			evWant: Tris{N: 74, X: 1, Y: 2, Z: -3}.Eigenvalues(),
		},
		{
			a:      Tris{N: 75, X: 1, Y: 2, Z: -3}.Matrix(),
			evWant: Tris{N: 75, X: 1, Y: 2, Z: -3}.Eigenvalues(),
		},

		{
			a:      Wilk4{}.Matrix(),
			evWant: Wilk4{}.Eigenvalues(),
		},
		{
			a:      Wilk12{}.Matrix(),
			evWant: Wilk12{}.Eigenvalues(),
			valTol: 1e-7,
		},
		{
			a:      Wilk20(0).Matrix(),
			evWant: Wilk20(0).Eigenvalues(),
		},
		{
			a:      Wilk20(1e-10).Matrix(),
			evWant: Wilk20(1e-10).Eigenvalues(),
			valTol: 1e-12,
			vecTol: 1e-12,
		},

		{
			a:      Zero(1).Matrix(),
			evWant: Zero(1).Eigenvalues(),
		},
		{
			a:      Zero(10).Matrix(),
			evWant: Zero(10).Eigenvalues(),
		},
		{
			a:      Zero(50).Matrix(),
			evWant: Zero(50).Eigenvalues(),
		},
		{
			a:      Zero(100).Matrix(),
			evWant: Zero(100).Eigenvalues(),
		},
	} {
		for _, jobvl := range []lapack.LeftEVJob{lapack.LeftEVCompute, lapack.LeftEVNone} {
			for _, jobvr := range []lapack.RightEVJob{lapack.RightEVCompute, lapack.RightEVNone} {
				for _, extra := range []int{0, 11} {
					for _, wl := range []worklen{minimumWork, mediumWork, optimumWork} {
						testDgeev(t, impl, strconv.Itoa(i), test, jobvl, jobvr, extra, wl)
					}
				}
			}
		}
	}

	for _, n := range []int{2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 20, 50, 51, 100, 101} {
		for _, jobvl := range []lapack.LeftEVJob{lapack.LeftEVCompute, lapack.LeftEVNone} {
			for _, jobvr := range []lapack.RightEVJob{lapack.RightEVCompute, lapack.RightEVNone} {
				for cas := 0; cas < 10; cas++ {
					// Create a block diagonal matrix with
					// random eigenvalues of random multiplicity.
					ev := make([]complex128, n)
					tmat := zeros(n, n, n)
					for i := 0; i < n; {
						re := rnd.NormFloat64()
						if i == n-1 || rnd.Float64() < 0.5 {
							// Real eigenvalue.
							nb := rnd.Intn(min(4, n-i)) + 1
							for k := 0; k < nb; k++ {
								tmat.Data[i*tmat.Stride+i] = re
								ev[i] = complex(re, 0)
								i++
							}
							continue
						}
						// Complex eigenvalue.
						im := rnd.NormFloat64()
						nb := rnd.Intn(min(4, (n-i)/2)) + 1
						for k := 0; k < nb; k++ {
							// 2×2 block for the complex eigenvalue.
							tmat.Data[i*tmat.Stride+i] = re
							tmat.Data[(i+1)*tmat.Stride+i+1] = re
							tmat.Data[(i+1)*tmat.Stride+i] = -im
							tmat.Data[i*tmat.Stride+i+1] = im
							ev[i] = complex(re, im)
							ev[i+1] = complex(re, -im)
							i += 2
						}
					}

					// Compute A = Q T Qᵀ where Q is an
					// orthogonal matrix.
					q := randomOrthogonal(n, rnd)
					tq := zeros(n, n, n)
					blas64.Gemm(blas.NoTrans, blas.Trans, 1, tmat, q, 0, tq)
					a := zeros(n, n, n)
					blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, q, tq, 0, a)

					test := dgeevTest{
						a:      a,
						evWant: ev,
						valTol: 1e-12,
						vecTol: 1e-7,
					}
					testDgeev(t, impl, "random", test, jobvl, jobvr, 0, optimumWork)
				}
			}
		}
	}
}

func testDgeev(t *testing.T, impl Dgeever, tc string, test dgeevTest, jobvl lapack.LeftEVJob, jobvr lapack.RightEVJob, extra int, wl worklen) {
	const defaultTol = 1e-12
	valTol := test.valTol
	if valTol == 0 {
		valTol = defaultTol
	}
	vecTol := test.vecTol
	if vecTol == 0 {
		vecTol = defaultTol
	}

	a := cloneGeneral(test.a)
	n := a.Rows

	var vl blas64.General
	if jobvl == lapack.LeftEVCompute {
		vl = nanGeneral(n, n, n)
	} else {
		vl.Stride = 1
	}

	var vr blas64.General
	if jobvr == lapack.RightEVCompute {
		vr = nanGeneral(n, n, n)
	} else {
		vr.Stride = 1
	}

	wr := make([]float64, n)
	wi := make([]float64, n)

	var lwork int
	switch wl {
	case minimumWork:
		if jobvl == lapack.LeftEVCompute || jobvr == lapack.RightEVCompute {
			lwork = max(1, 4*n)
		} else {
			lwork = max(1, 3*n)
		}
	case mediumWork:
		work := make([]float64, 1)
		impl.Dgeev(jobvl, jobvr, n, a.Data, a.Stride, wr, wi, vl.Data, vl.Stride, vr.Data, vr.Stride, work, -1)
		if jobvl == lapack.LeftEVCompute || jobvr == lapack.RightEVCompute {
			lwork = (int(work[0]) + 4*n) / 2
		} else {
			lwork = (int(work[0]) + 3*n) / 2
		}
		lwork = max(1, lwork)
	case optimumWork:
		work := make([]float64, 1)
		impl.Dgeev(jobvl, jobvr, n, a.Data, a.Stride, wr, wi, vl.Data, vl.Stride, vr.Data, vr.Stride, work, -1)
		lwork = int(work[0])
	}
	work := make([]float64, lwork)

	first := 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, work=%v",
		tc, n, jobvl, jobvr, extra, wl)

	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 first > 0 {
		t.Logf("%v: all eigenvalues haven't been computed, first=%v", prefix, first)
	}

	// Check that conjugate pair eigevalues are ordered correctly.
	for i := first; 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
	}

	// Check the computed eigenvalues against provided known eigenvalues.
	if test.evWant != nil {
		used := make([]bool, n)
		for i := first; i < n; i++ {
			evGot := complex(wr[i], wi[i])
			idx := -1
			for k, evWant := range test.evWant {
				if !used[k] && cmplx.Abs(evWant-evGot) < valTol {
					idx = k
					used[k] = true
					break
				}
			}
			if idx == -1 {
				t.Errorf("%v: unexpected eigenvalue %v", prefix, evGot)
			}
		}
	}

	if first > 0 || (jobvl == lapack.LeftEVNone && jobvr == lapack.RightEVNone) {
		// No eigenvectors have been computed.
		return
	}

	// Check that the columns of VL and VR are eigenvectors that correspond
	// to the computed eigenvalues.
	for k := 0; k < n; {
		if wi[k] == 0 {
			if jobvl == lapack.LeftEVCompute {
				ev := columnOf(vl, k)
				if !isLeftEigenvectorOf(test.a, ev, nil, complex(wr[k], 0), vecTol) {
					t.Errorf("%v: VL[:,%v] is not left real eigenvector",
						prefix, k)
				}

				norm := floats.Norm(ev, 2)
				if math.Abs(norm-1) >= defaultTol {
					t.Errorf("%v: norm of left real eigenvector %v not equal to 1: got %v",
						prefix, k, norm)
				}
			}
			if jobvr == lapack.RightEVCompute {
				ev := columnOf(vr, k)
				if !isRightEigenvectorOf(test.a, ev, nil, complex(wr[k], 0), vecTol) {
					t.Errorf("%v: VR[:,%v] is not right real eigenvector",
						prefix, k)
				}

				norm := floats.Norm(ev, 2)
				if math.Abs(norm-1) >= defaultTol {
					t.Errorf("%v: norm of right real eigenvector %v not equal to 1: got %v",
						prefix, k, norm)
				}
			}
			k++
		} else {
			if jobvl == lapack.LeftEVCompute {
				evre := columnOf(vl, k)
				evim := columnOf(vl, k+1)
				if !isLeftEigenvectorOf(test.a, evre, evim, complex(wr[k], wi[k]), vecTol) {
					t.Errorf("%v: VL[:,%v:%v] is not left complex eigenvector",
						prefix, k, k+1)
				}
				floats.Scale(-1, evim)
				if !isLeftEigenvectorOf(test.a, evre, evim, complex(wr[k+1], wi[k+1]), vecTol) {
					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) > defaultTol {
					t.Errorf("%v: norm of left complex eigenvector %v not equal to 1: got %v",
						prefix, k, norm)
				}
			}
			if jobvr == lapack.RightEVCompute {
				evre := columnOf(vr, k)
				evim := columnOf(vr, k+1)
				if !isRightEigenvectorOf(test.a, evre, evim, complex(wr[k], wi[k]), vecTol) {
					t.Errorf("%v: VR[:,%v:%v] is not right complex eigenvector",
						prefix, k, k+1)
				}
				floats.Scale(-1, evim)
				if !isRightEigenvectorOf(test.a, evre, evim, complex(wr[k+1], wi[k+1]), vecTol) {
					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) > defaultTol {
					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
		}
	}
}

func dgeevTestForAntisymRandom(n int, rnd *rand.Rand) dgeevTest {
	a := NewAntisymRandom(n, rnd)
	return dgeevTest{
		a:      a.Matrix(),
		evWant: a.Eigenvalues(),
	}
}