--- /dev/null
+// 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"
+ "sort"
+ "testing"
+
+ "golang.org/x/exp/rand"
+
+ "gonum.org/v1/gonum/floats"
+)
+
+type Dsterfer interface {
+ Dgetrfer
+ Dsterf(n int, d, e []float64) (ok bool)
+}
+
+func DsterfTest(t *testing.T, impl Dsterfer) {
+ // Hand coded tests.
+ for cas, test := range []struct {
+ d []float64
+ e []float64
+ n int
+
+ ans []float64
+ }{
+ // Computed from Fortran code.
+ {
+ d: []float64{1, 3, 4, 6},
+ e: []float64{2, 4, 5},
+ n: 4,
+ ans: []float64{11.046227528488854, 4.795922173417400, -2.546379458290125, 0.704229756383872},
+ },
+ } {
+ n := test.n
+ d := make([]float64, len(test.d))
+ copy(d, test.d)
+ e := make([]float64, len(test.e))
+ copy(e, test.e)
+ ok := impl.Dsterf(n, d, e)
+ if !ok {
+ t.Errorf("Case %d, Eigenvalue decomposition failed", cas)
+ continue
+ }
+ ans := make([]float64, len(test.ans))
+ copy(ans, test.ans)
+ sort.Float64s(ans)
+ if !floats.EqualApprox(ans, d, 1e-10) {
+ t.Errorf("eigenvalue mismatch")
+ }
+ }
+
+ rnd := rand.New(rand.NewSource(1))
+ // Probabilistic tests.
+ for _, n := range []int{4, 6, 10} {
+ for cas := 0; cas < 10; cas++ {
+ d := make([]float64, n)
+ for i := range d {
+ d[i] = rnd.NormFloat64()
+ }
+ dCopy := make([]float64, len(d))
+ copy(dCopy, d)
+ e := make([]float64, n-1)
+ for i := range e {
+ e[i] = rnd.NormFloat64()
+ }
+ eCopy := make([]float64, len(e))
+ copy(eCopy, e)
+
+ ok := impl.Dsterf(n, d, e)
+ if !ok {
+ t.Errorf("Eigenvalue decomposition failed")
+ continue
+ }
+
+ // Test that the eigenvalues are sorted.
+ if !sort.Float64sAreSorted(d) {
+ t.Errorf("Values are not sorted")
+ }
+
+ // Construct original tridagional matrix.
+ lda := n
+ a := make([]float64, n*lda)
+ for i := 0; i < n; i++ {
+ a[i*lda+i] = dCopy[i]
+ if i != n-1 {
+ a[i*lda+i+1] = eCopy[i]
+ a[(i+1)*lda+i] = eCopy[i]
+ }
+ }
+
+ asub := make([]float64, len(a))
+ ipiv := make([]int, n)
+
+ // Test that they are actually eigenvalues by computing the
+ // determinant of A - λI.
+ // TODO(btracey): Replace this test with a more numerically stable
+ // test.
+ for _, lambda := range d {
+ copy(asub, a)
+ for i := 0; i < n; i++ {
+ asub[i*lda+i] -= lambda
+ }
+
+ // Compute LU.
+ ok := impl.Dgetrf(n, n, asub, lda, ipiv)
+ if !ok {
+ // Definitely singular.
+ continue
+ }
+ // Compute determinant.
+ var logdet float64
+ for i := 0; i < n; i++ {
+ v := asub[i*lda+i]
+ logdet += math.Log(math.Abs(v))
+ }
+ if math.Exp(logdet) > 2 {
+ t.Errorf("Incorrect singular value. n = %d, cas = %d, det = %v", n, cas, math.Exp(logdet))
+ }
+ }
+ }
+ }
+}