| // 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) |
| } |
| } |
| } |
| } |