blob: 3b3e0eb21e6b774a66184bb9058c55b66d1115e8 [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 ( "math" "testing" "golang.org/x/exp/rand" "gonum.org/v1/gonum/floats" ) type Dlaqr1er interface { Dlaqr1(n int, h []float64, ldh int, sr1, si1, sr2, si2 float64, v []float64) } func Dlaqr1Test(t *testing.T, impl Dlaqr1er) { rnd := rand.New(rand.NewSource(1)) for _, n := range []int{2, 3} { for _, ldh := range []int{n, n + 1, n + 10} { for _, cas := range []int{1, 2} { for k := 0; k < 100; k++ { v := make([]float64, n) for i := range v { v[i] = math.NaN() } h := make([]float64, n*(n-1)*ldh) for i := range h { h[i] = math.NaN() } for i := 0; i < n; i++ { for j := 0; j < n; j++ { h[i*ldh+j] = rnd.NormFloat64() } } var sr1, sr2, si1, si2 float64 if cas == 1 { sr1 = rnd.NormFloat64() sr2 = sr1 si1 = rnd.NormFloat64() si2 = -si1 } else { sr1 = rnd.NormFloat64() sr2 = rnd.NormFloat64() si1 = 0 si2 = 0 } impl.Dlaqr1(n, h, ldh, sr1, si1, sr2, si2, v) // Matrix H - s1*I. h1 := make([]complex128, n*n) for i := 0; i < n; i++ { for j := 0; j < n; j++ { h1[i*n+j] = complex(h[i*ldh+j], 0) if i == j { h1[i*n+j] -= complex(sr1, si1) } } } // First column of H - s2*I. h2 := make([]complex128, n) for i := 0; i < n; i++ { h2[i] = complex(h[i*ldh], 0) } h2[0] -= complex(sr2, si2) wantv := make([]float64, n) // Multiply (H-s1*I)*(H-s2*I) to get a tentative // wantv. for i := 0; i < n; i++ { for j := 0; j < n; j++ { wantv[i] += real(h1[i*n+j] * h2[j]) } } // Get the unknown scale. scale := v[0] / wantv[0] // Compute the actual wantv. floats.Scale(scale, wantv) // The scale must be the same for all elements. if floats.Distance(wantv, v, math.Inf(1)) > 1e-13 { t.Errorf("n = %v, ldh = %v, case = %v: Unexpected value of v: got %v, want %v", n, ldh, cas, v, wantv) } } } } } }