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
|
// 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"
"testing"
)
type Dlartger interface {
Dlartg(f, g float64) (cs, sn, r float64)
}
func DlartgTest(t *testing.T, impl Dlartger) {
const tol = 20 * ulp
values := []float64{
-safmax,
-1 / ulp,
-1,
-1.0 / 3,
-ulp,
-safmin,
0,
safmin,
ulp,
1.0 / 3,
1,
1 / ulp,
safmax,
math.Inf(-1),
math.Inf(1),
math.NaN(),
}
for _, f := range values {
for _, g := range values {
name := fmt.Sprintf("Case f=%v,g=%v", f, g)
// Generate a plane rotation so that
// [ cs sn] * [f] = [r]
// [-sn cs] [g] = [0]
// where cs*cs + sn*sn = 1.
cs, sn, r := impl.Dlartg(f, g)
switch {
case math.IsNaN(f) || math.IsNaN(g):
if !math.IsNaN(r) {
t.Errorf("%v: unexpected r=%v; want NaN", name, r)
}
case math.IsInf(f, 0) || math.IsInf(g, 0):
if !math.IsNaN(r) && !math.IsInf(r, 0) {
t.Errorf("%v: unexpected r=%v; want NaN or Inf", name, r)
}
default:
d := math.Max(math.Abs(f), math.Abs(g))
d = math.Min(math.Max(safmin, d), safmax)
fs := f / d
gs := g / d
rs := r / d
// Check that cs*f + sn*g = r.
rnorm := math.Abs(rs)
if rnorm == 0 {
rnorm = math.Max(math.Abs(fs), math.Abs(gs))
if rnorm == 0 {
rnorm = 1
}
}
resid := math.Abs(rs-(cs*fs+sn*gs)) / rnorm
if resid > tol {
t.Errorf("%v: cs*f + sn*g != r; resid=%v", name, resid)
}
// Check that -sn*f + cs*g = 0.
resid = math.Abs(-sn*fs + cs*gs)
if resid > tol {
t.Errorf("%v: -sn*f + cs*g != 0; resid=%v", name, resid)
}
// Check that cs*cs + sn*sn = 1.
resid = math.Abs(1 - (cs*cs + sn*sn))
if resid > tol {
t.Errorf("%v: cs*cs + sn*sn != 1; resid=%v", name, resid)
}
// Check that cs is non-negative.
if math.Abs(f) > math.Abs(g) && cs < 0 {
t.Errorf("%v: cs is negative; cs=%v", name, cs)
}
}
}
}
}
|
