vendor/gonum.org/v1/gonum/lapack/testlapack/dsterf.go

   1 // Copyright ©2016 The Gonum Authors. All rights reserved.
   2 // Use of this source code is governed by a BSD-style
   3 // license that can be found in the LICENSE file.
   4
   5 package testlapack
   6
   7 import (
   8         "math"
   9         "sort"
  10         "testing"
  11
  12         "golang.org/x/exp/rand"
  13
  14         "gonum.org/v1/gonum/floats"
  15 )
  16
  17 type Dsterfer interface {
  18         Dgetrfer
  19         Dsterf(n int, d, e []float64) (ok bool)
  20 }
  21
  22 func DsterfTest(t *testing.T, impl Dsterfer) {
  23         // Hand coded tests.
  24         for cas, test := range []struct {
  25                 d []float64
  26                 e []float64
  27                 n int
  28
  29                 ans []float64
  30         }{
  31                 // Computed from Fortran code.
  32                 {
  33                         d:   []float64{1, 3, 4, 6},
  34                         e:   []float64{2, 4, 5},
  35                         n:   4,
  36                         ans: []float64{11.046227528488854, 4.795922173417400, -2.546379458290125, 0.704229756383872},
  37                 },
  38         } {
  39                 n := test.n
  40                 d := make([]float64, len(test.d))
  41                 copy(d, test.d)
  42                 e := make([]float64, len(test.e))
  43                 copy(e, test.e)
  44                 ok := impl.Dsterf(n, d, e)
  45                 if !ok {
  46                         t.Errorf("Case %d, Eigenvalue decomposition failed", cas)
  47                         continue
  48                 }
  49                 ans := make([]float64, len(test.ans))
  50                 copy(ans, test.ans)
  51                 sort.Float64s(ans)
  52                 if !floats.EqualApprox(ans, d, 1e-10) {
  53                         t.Errorf("eigenvalue mismatch")
  54                 }
  55         }
  56
  57         rnd := rand.New(rand.NewSource(1))
  58         // Probabilistic tests.
  59         for _, n := range []int{4, 6, 10} {
  60                 for cas := 0; cas < 10; cas++ {
  61                         d := make([]float64, n)
  62                         for i := range d {
  63                                 d[i] = rnd.NormFloat64()
  64                         }
  65                         dCopy := make([]float64, len(d))
  66                         copy(dCopy, d)
  67                         e := make([]float64, n-1)
  68                         for i := range e {
  69                                 e[i] = rnd.NormFloat64()
  70                         }
  71                         eCopy := make([]float64, len(e))
  72                         copy(eCopy, e)
  73
  74                         ok := impl.Dsterf(n, d, e)
  75                         if !ok {
  76                                 t.Errorf("Eigenvalue decomposition failed")
  77                                 continue
  78                         }
  79
  80                         // Test that the eigenvalues are sorted.
  81                         if !sort.Float64sAreSorted(d) {
  82                                 t.Errorf("Values are not sorted")
  83                         }
  84
  85                         // Construct original tridagional matrix.
  86                         lda := n
  87                         a := make([]float64, n*lda)
  88                         for i := 0; i < n; i++ {
  89                                 a[i*lda+i] = dCopy[i]
  90                                 if i != n-1 {
  91                                         a[i*lda+i+1] = eCopy[i]
  92                                         a[(i+1)*lda+i] = eCopy[i]
  93                                 }
  94                         }
  95
  96                         asub := make([]float64, len(a))
  97                         ipiv := make([]int, n)
  98
  99                         // Test that they are actually eigenvalues by computing the
 100                         // determinant of A - λI.
 101                         // TODO(btracey): Replace this test with a more numerically stable
 102                         // test.
 103                         for _, lambda := range d {
 104                                 copy(asub, a)
 105                                 for i := 0; i < n; i++ {
 106                                         asub[i*lda+i] -= lambda
 107                                 }
 108
 109                                 // Compute LU.
 110                                 ok := impl.Dgetrf(n, n, asub, lda, ipiv)
 111                                 if !ok {
 112                                         // Definitely singular.
 113                                         continue
 114                                 }
 115                                 // Compute determinant.
 116                                 var logdet float64
 117                                 for i := 0; i < n; i++ {
 118                                         v := asub[i*lda+i]
 119                                         logdet += math.Log(math.Abs(v))
 120                                 }
 121                                 if math.Exp(logdet) > 2 {
 122                                         t.Errorf("Incorrect singular value. n = %d, cas = %d, det = %v", n, cas, math.Exp(logdet))
 123                                 }
 124                         }
 125                 }
 126         }
 127 }