1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
|
// 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)
}
}
}
}
}
}
|
