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
128
129
130
131
132
133
134
135
136
137
138
139
140
|
// 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 gonum
import "math"
// Dlasq5 computes one dqds transform in ping-pong form.
// i0 and n0 are zero-indexed.
//
// Dlasq5 is an internal routine. It is exported for testing purposes.
func (impl Implementation) Dlasq5(i0, n0 int, z []float64, pp int, tau, sigma float64) (i0Out, n0Out, ppOut int, tauOut, sigmaOut, dmin, dmin1, dmin2, dn, dnm1, dnm2 float64) {
// The lapack function has inputs for ieee and eps, but Go requires ieee so
// these are unnecessary.
switch {
case i0 < 0:
panic(i0LT0)
case n0 < 0:
panic(n0LT0)
case len(z) < 4*n0:
panic(shortZ)
case pp != 0 && pp != 1:
panic(badPp)
}
if n0-i0-1 <= 0 {
return i0, n0, pp, tau, sigma, dmin, dmin1, dmin2, dn, dnm1, dnm2
}
eps := dlamchP
dthresh := eps * (sigma + tau)
if tau < dthresh*0.5 {
tau = 0
}
var j4 int
var emin float64
if tau != 0 {
j4 = 4*i0 + pp
emin = z[j4+4]
d := z[j4] - tau
dmin = d
// In the reference there are code paths that actually return this value.
// dmin1 = -z[j4]
if pp == 0 {
for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 {
j4 := j4loop - 1
z[j4-2] = d + z[j4-1]
tmp := z[j4+1] / z[j4-2]
d = d*tmp - tau
dmin = math.Min(dmin, d)
z[j4] = z[j4-1] * tmp
emin = math.Min(z[j4], emin)
}
} else {
for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 {
j4 := j4loop - 1
z[j4-3] = d + z[j4]
tmp := z[j4+2] / z[j4-3]
d = d*tmp - tau
dmin = math.Min(dmin, d)
z[j4-1] = z[j4] * tmp
emin = math.Min(z[j4-1], emin)
}
}
// Unroll the last two steps.
dnm2 = d
dmin2 = dmin
j4 = 4*((n0+1)-2) - pp - 1
j4p2 := j4 + 2*pp - 1
z[j4-2] = dnm2 + z[j4p2]
z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2])
dnm1 = z[j4p2+2]*(dnm2/z[j4-2]) - tau
dmin = math.Min(dmin, dnm1)
dmin1 = dmin
j4 += 4
j4p2 = j4 + 2*pp - 1
z[j4-2] = dnm1 + z[j4p2]
z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2])
dn = z[j4p2+2]*(dnm1/z[j4-2]) - tau
dmin = math.Min(dmin, dn)
} else {
// This is the version that sets d's to zero if they are small enough.
j4 = 4*(i0+1) + pp - 4
emin = z[j4+4]
d := z[j4] - tau
dmin = d
// In the reference there are code paths that actually return this value.
// dmin1 = -z[j4]
if pp == 0 {
for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 {
j4 := j4loop - 1
z[j4-2] = d + z[j4-1]
tmp := z[j4+1] / z[j4-2]
d = d*tmp - tau
if d < dthresh {
d = 0
}
dmin = math.Min(dmin, d)
z[j4] = z[j4-1] * tmp
emin = math.Min(z[j4], emin)
}
} else {
for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 {
j4 := j4loop - 1
z[j4-3] = d + z[j4]
tmp := z[j4+2] / z[j4-3]
d = d*tmp - tau
if d < dthresh {
d = 0
}
dmin = math.Min(dmin, d)
z[j4-1] = z[j4] * tmp
emin = math.Min(z[j4-1], emin)
}
}
// Unroll the last two steps.
dnm2 = d
dmin2 = dmin
j4 = 4*((n0+1)-2) - pp - 1
j4p2 := j4 + 2*pp - 1
z[j4-2] = dnm2 + z[j4p2]
z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2])
dnm1 = z[j4p2+2]*(dnm2/z[j4-2]) - tau
dmin = math.Min(dmin, dnm1)
dmin1 = dmin
j4 += 4
j4p2 = j4 + 2*pp - 1
z[j4-2] = dnm1 + z[j4p2]
z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2])
dn = z[j4p2+2]*(dnm1/z[j4-2]) - tau
dmin = math.Min(dmin, dn)
}
z[j4+2] = dn
z[4*(n0+1)-pp-1] = emin
return i0, n0, pp, tau, sigma, dmin, dmin1, dmin2, dn, dnm1, dnm2
}
|
