1 // Copyright ©2015 The Gonum Authors. All rights reserved.
2 // Use of this source code is governed by a BSD-style
3 // license that can be found in the LICENSE file.
9 // Dlasq5 computes one dqds transform in ping-pong form.
10 // i0 and n0 are zero-indexed.
12 // Dlasq5 is an internal routine. It is exported for testing purposes.
13 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) {
14 // The lapack function has inputs for ieee and eps, but Go requires ieee so
15 // these are unnecessary.
17 return i0, n0, pp, tau, sigma, dmin, dmin1, dmin2, dn, dnm1, dnm2
20 dthresh := eps * (sigma + tau)
21 if tau < dthresh*0.5 {
31 // In the reference there are code paths that actually return this value.
34 for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 {
37 tmp := z[j4+1] / z[j4-2]
39 dmin = math.Min(dmin, d)
41 emin = math.Min(z[j4], emin)
44 for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 {
47 tmp := z[j4+2] / z[j4-3]
49 dmin = math.Min(dmin, d)
51 emin = math.Min(z[j4-1], emin)
54 // Unroll the last two steps.
57 j4 = 4*((n0+1)-2) - pp - 1
59 z[j4-2] = dnm2 + z[j4p2]
60 z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2])
61 dnm1 = z[j4p2+2]*(dnm2/z[j4-2]) - tau
62 dmin = math.Min(dmin, dnm1)
67 z[j4-2] = dnm1 + z[j4p2]
68 z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2])
69 dn = z[j4p2+2]*(dnm1/z[j4-2]) - tau
70 dmin = math.Min(dmin, dn)
72 // This is the version that sets d's to zero if they are small enough.
73 j4 = 4*(i0+1) + pp - 4
77 // In the reference there are code paths that actually return this value.
80 for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 {
83 tmp := z[j4+1] / z[j4-2]
88 dmin = math.Min(dmin, d)
90 emin = math.Min(z[j4], emin)
93 for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 {
96 tmp := z[j4+2] / z[j4-3]
101 dmin = math.Min(dmin, d)
102 z[j4-1] = z[j4] * tmp
103 emin = math.Min(z[j4-1], emin)
106 // Unroll the last two steps.
109 j4 = 4*((n0+1)-2) - pp - 1
110 j4p2 := j4 + 2*pp - 1
111 z[j4-2] = dnm2 + z[j4p2]
112 z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2])
113 dnm1 = z[j4p2+2]*(dnm2/z[j4-2]) - tau
114 dmin = math.Min(dmin, dnm1)
119 z[j4-2] = dnm1 + z[j4p2]
120 z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2])
121 dn = z[j4p2+2]*(dnm1/z[j4-2]) - tau
122 dmin = math.Min(dmin, dn)
125 z[4*(n0+1)-pp-1] = emin
126 return i0, n0, pp, tau, sigma, dmin, dmin1, dmin2, dn, dnm1, dnm2