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 // Dlasq6 computes one dqd transform in ping-pong form with protection against
10 // overflow and underflow. z has length at least 4*(n0+1) and holds the qd array.
11 // i0 is the zero-based first index.
12 // n0 is the zero-based last index.
14 // Dlasq6 is an internal routine. It is exported for testing purposes.
15 func (impl Implementation) Dlasq6(i0, n0 int, z []float64, pp int) (dmin, dmin1, dmin2, dn, dnm1, dnm2 float64) {
16 if len(z) < 4*(n0+1) {
20 return dmin, dmin1, dmin2, dn, dnm1, dnm2
23 j4 := 4*(i0+1) + pp - 4 // -4 rather than -3 for zero indexing
28 for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 {
29 j4 := j4loop - 1 // Translate back to zero-indexed.
36 } else if safmin*z[j4+1] < z[j4-2] && safmin*z[j4-2] < z[j4+1] {
37 tmp := z[j4+1] / z[j4-2]
41 z[j4] = z[j4+1] * (z[j4-1] / z[j4-2])
42 d = z[j4+1] * (d / z[j4-2])
44 dmin = math.Min(dmin, d)
45 emin = math.Min(emin, z[j4])
48 for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 {
56 } else if safmin*z[j4+2] < z[j4-3] && safmin*z[j4-3] < z[j4+2] {
57 tmp := z[j4+2] / z[j4-3]
61 z[j4-1] = z[j4+2] * (z[j4] / z[j4-3])
62 d = z[j4+2] * (d / z[j4-3])
64 dmin = math.Min(dmin, d)
65 emin = math.Min(emin, z[j4-1])
68 // Unroll last two steps.
71 j4 = 4*(n0-1) - pp - 1
73 z[j4-2] = dnm2 + z[j4p2]
79 } else if safmin*z[j4p2+2] < z[j4-2] && safmin*z[j4-2] < z[j4p2+2] {
80 tmp := z[j4p2+2] / z[j4-2]
84 z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2])
85 dnm1 = z[j4p2+2] * (dnm2 / z[j4-2])
87 dmin = math.Min(dmin, dnm1)
91 z[j4-2] = dnm1 + z[j4p2]
97 } else if safmin*z[j4p2+2] < z[j4-2] && safmin*z[j4-2] < z[j4p2+2] {
98 tmp := z[j4p2+2] / z[j4-2]
102 z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2])
103 dn = z[j4p2+2] * (dnm1 / z[j4-2])
105 dmin = math.Min(dmin, dn)
107 z[4*(n0+1)-pp-1] = emin
108 return dmin, dmin1, dmin2, dn, dnm1, dnm2