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
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
|
// Copyright ©2015 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/floats"
)
type Dlasq1er interface {
Dlasq1(n int, d, e, work []float64) int
Dgebrd(m, n int, a []float64, lda int, d, e, tauQ, tauP, work []float64, lwork int)
}
func Dlasq1Test(t *testing.T, impl Dlasq1er) {
const tol = 1e-14
rnd := rand.New(rand.NewSource(1))
for _, n := range []int{0, 1, 2, 3, 4, 5, 8, 10, 30, 50} {
for typ := 0; typ <= 7; typ++ {
name := fmt.Sprintf("n=%v,typ=%v", n, typ)
// Generate a diagonal matrix D with positive entries.
d := make([]float64, n)
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.
for i := 0; i < n; i++ {
if i == 0 {
d[0] = 1
} else {
d[i] = 1 - (1-dlamchE)*float64(i)/float64(n-1)
}
}
case 3, 4, 5:
// A diagonal matrix with geometrically spaced entries 1, ..., eps.
for i := 0; i < n; i++ {
if i == 0 {
d[0] = 1
} else {
d[i] = math.Pow(dlamchE, float64(i)/float64(n-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.
for i := range d {
if i == 0 {
d[i] = 1
} else {
d[i] = dlamchE
}
}
case 7:
// Diagonal matrix with random entries.
for i := range d {
d[i] = math.Abs(rnd.NormFloat64())
}
}
dWant := make([]float64, n)
copy(dWant, d)
sort.Sort(sort.Reverse(sort.Float64Slice(dWant)))
// Allocate work slice to the maximum length needed below.
work := make([]float64, max(1, 4*n))
// Generate an n×n matrix A by pre- and post-multiplying D with
// random orthogonal matrices:
// A = U*D*V.
lda := max(1, n)
a := make([]float64, n*lda)
Dlagge(n, n, 0, 0, d, a, lda, rnd, work)
// Reduce A to bidiagonal form B represented by the diagonal d and
// off-diagonal e.
tauQ := make([]float64, n)
tauP := make([]float64, n)
e := make([]float64, max(0, n-1))
impl.Dgebrd(n, n, a, lda, d, e, tauQ, tauP, work, len(work))
// Compute the singular values of B.
for i := range work {
work[i] = math.NaN()
}
info := impl.Dlasq1(n, d, e, work)
if info != 0 {
t.Fatalf("%v: Dlasq1 returned non-zero info=%v", name, info)
}
if n == 0 {
continue
}
if !sort.IsSorted(sort.Reverse(sort.Float64Slice(d))) {
t.Errorf("%v: singular values not sorted", name)
}
diff := floats.Distance(d, dWant, math.Inf(1))
if diff > tol*floats.Max(dWant) {
t.Errorf("%v: unexpected result; diff=%v", name, diff)
}
}
}
}
|
