github.com/gopherd/gonum@v0.0.4/lapack/testlapack/matgen_test.go (about) 1 // Copyright ©2017 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 "testing" 9 10 "math/rand" 11 12 "github.com/gopherd/gonum/blas/blas64" 13 ) 14 15 func TestDlagsy(t *testing.T) { 16 const tol = 1e-14 17 rnd := rand.New(rand.NewSource(1)) 18 for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 50} { 19 for _, lda := range []int{0, 2*n + 1} { 20 if lda == 0 { 21 lda = max(1, n) 22 } 23 // D is the identity matrix I. 24 d := make([]float64, n) 25 for i := range d { 26 d[i] = 1 27 } 28 // Allocate an n×n symmetric matrix A and fill it with NaNs. 29 a := nanSlice(n * lda) 30 work := make([]float64, 2*n) 31 // Compute A = U * D * Uᵀ where U is a random orthogonal matrix. 32 Dlagsy(n, 0, d, a, lda, rnd, work) 33 // A should be the identity matrix because 34 // A = U * D * Uᵀ = U * I * Uᵀ = U * Uᵀ = I. 35 dist := distFromIdentity(n, a, lda) 36 if dist > tol { 37 t.Errorf("Case n=%v,lda=%v: |A-I|=%v is too large", n, lda, dist) 38 } 39 } 40 } 41 } 42 43 func TestDlagge(t *testing.T) { 44 const tol = 1e-14 45 rnd := rand.New(rand.NewSource(1)) 46 for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 50} { 47 for _, lda := range []int{0, 2*n + 1} { 48 if lda == 0 { 49 lda = max(1, n) 50 } 51 d := make([]float64, n) 52 for i := range d { 53 d[i] = 1 54 } 55 a := blas64.General{ 56 Rows: n, 57 Cols: n, 58 Stride: lda, 59 Data: nanSlice(n * lda), 60 } 61 work := make([]float64, a.Rows+a.Cols) 62 63 Dlagge(a.Rows, a.Cols, 0, 0, d, a.Data, a.Stride, rnd, work) 64 65 if resid := residualOrthogonal(a, false); resid > tol { 66 t.Errorf("Case n=%v,lda=%v: unexpected result", n, lda) 67 } 68 } 69 } 70 } 71 72 func TestRandomOrthogonal(t *testing.T) { 73 const tol = 1e-14 74 rnd := rand.New(rand.NewSource(1)) 75 for n := 1; n <= 20; n++ { 76 q := randomOrthogonal(n, rnd) 77 if resid := residualOrthogonal(q, false); resid > tol { 78 t.Errorf("Case n=%v: Q not orthogonal; resid=%v, want<=%v", n, resid, tol) 79 } 80 } 81 }