blob: 2d546a5b71c499571096e7a246326cd2cacba4ee [file] [log] [blame]
 // 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" "sort" "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 Dsterfer interface { Dsteqrer Dlansyer Dsterf(n int, d, e []float64) (ok bool) } func DsterfTest(t *testing.T, impl Dsterfer) { const tol = 1e-14 // Tests with precomputed eigenvalues. for cas, test := range []struct { d []float64 e []float64 n int want []float64 }{ { d: []float64{1, 3, 4, 6}, e: []float64{2, 4, 5}, n: 4, // Computed from original Fortran code. want: []float64{11.046227528488854, 4.795922173417400, -2.546379458290125, 0.704229756383872}, }, } { n := test.n got := make([]float64, len(test.d)) copy(got, test.d) e := make([]float64, len(test.e)) copy(e, test.e) ok := impl.Dsterf(n, got, e) if !ok { t.Errorf("Case %d, n=%v: Dsterf failed", cas, n) continue } want := make([]float64, len(test.want)) copy(want, test.want) sort.Float64s(want) diff := floats.Distance(got, want, math.Inf(1)) if diff > tol { t.Errorf("Case %d, n=%v: unexpected result, |dGot-dWant|=%v", cas, n, diff) } } rnd := rand.New(rand.NewSource(1)) // Probabilistic tests. for _, n := range []int{0, 1, 2, 3, 4, 5, 6, 10, 50} { for typ := 0; typ <= 8; typ++ { d := make([]float64, n) var e []float64 if n > 1 { e = make([]float64, n-1) } // Generate a tridiagonal matrix A. switch typ { case 0: // The zero matrix. case 1: // The identity matrix. for i := range d { d[i] = 1 } case 2: // A diagonal matrix with evenly spaced entries // 1, ..., eps and random signs. for i := 0; i < n; i++ { if i == 0 { d[i] = 1 } else { d[i] = 1 - (1-dlamchE)*float64(i)/float64(n-1) } if rnd.Float64() < 0.5 { d[i] *= -1 } } case 3, 4, 5: // A diagonal matrix with geometrically spaced entries // 1, ..., eps and random signs. for i := 0; i < n; i++ { if i == 0 { d[i] = 1 } else { d[i] = math.Pow(dlamchE, float64(i)/float64(n-1)) } if rnd.Float64() < 0.5 { d[i] *= -1 } } switch typ { case 4: // Multiply by SQRT(overflow threshold). floats.Scale(math.Sqrt(1/dlamchS), d) case 5: // Multiply by SQRT(underflow threshold). floats.Scale(math.Sqrt(dlamchS), d) } case 6: // A diagonal matrix with "clustered" entries 1, eps, ..., eps // and random signs. for i := range d { if i == 0 { d[i] = 1 } else { d[i] = dlamchE } } for i := range d { if rnd.Float64() < 0.5 { d[i] *= -1 } } case 7: // Diagonal matrix with random entries. for i := range d { d[i] = rnd.NormFloat64() } case 8: // Random symmetric tridiagonal matrix. for i := range d { d[i] = rnd.NormFloat64() } for i := range e { e[i] = rnd.NormFloat64() } } eCopy := make([]float64, len(e)) copy(eCopy, e) name := fmt.Sprintf("n=%d,type=%d", n, typ) // Compute the eigenvalues of A using Dsterf. dGot := make([]float64, len(d)) copy(dGot, d) ok := impl.Dsterf(n, dGot, e) if !ok { t.Errorf("%v: Dsterf failed", name) continue } if n == 0 { continue } // Test that the eigenvalues are sorted. if !sort.Float64sAreSorted(dGot) { t.Errorf("%v: eigenvalues are not sorted", name) continue } // Compute the expected eigenvalues of A using Dsteqr. dWant := make([]float64, len(d)) copy(dWant, d) copy(e, eCopy) z := nanGeneral(n, n, n) ok = impl.Dsteqr(lapack.EVTridiag, n, dWant, e, z.Data, z.Stride, make([]float64, 2*n)) if !ok { t.Fatalf("%v: computing reference solution using Dsteqr failed", name) continue } if resid := residualOrthogonal(z, false); resid > tol*float64(n) { t.Errorf("%v: Z is not orthogonal; resid=%v, want<=%v", name, resid, tol*float64(n)) } // Check whether eigenvalues from Dsteqr and Dsterf (which use // different algorithms) are equal. var diff float64 for i := range dGot { diffAbs := math.Abs(dGot[i] - dWant[i]) diff = math.Max(diff, diffAbs) } if diff > tol { t.Errorf("%v: unexpected result; |dGot-dWant|=%v", name, diff) } // Construct A as a symmetric dense matrix and compute its 1-norm. copy(e, eCopy) lda := n a := make([]float64, n*lda) var anorm, tmp float64 for i := 0; i < n-1; i++ { a[i*lda+i] = d[i] a[i*lda+i+1] = e[i] tmp2 := math.Abs(e[i]) anorm = math.Max(anorm, math.Abs(d[i])+tmp+tmp2) tmp = tmp2 } a[(n-1)*lda+n-1] = d[n-1] anorm = math.Max(anorm, math.Abs(d[n-1])+tmp) // Compute A - Z D Zᵀ. The result should be the zero matrix. bi := blas64.Implementation() for i := 0; i < n; i++ { bi.Dsyr(blas.Upper, n, -dGot[i], z.Data[i:], z.Stride, a, lda) } // Compute |A - Z D Zᵀ|. wnorm := impl.Dlansy(lapack.MaxColumnSum, blas.Upper, n, a, lda, make([]float64, n)) // Compute diff := |A - Z D Zᵀ| / (|A| N). if anorm > wnorm { diff = wnorm / anorm / float64(n) } else { if anorm < 1 { diff = math.Min(wnorm, float64(n)*anorm) / anorm / float64(n) } else { diff = math.Min(wnorm/anorm, float64(n)) / float64(n) } } // Check whether diff is small. if diff > tol { t.Errorf("%v: unexpected result; |A - Z D Zᵀ|/(|A| n)=%v", name, diff) } } } }