github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/testlapack/dsterf.go (about)

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