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. 4 5package gonum 6 7import "math" 8 9// Dlasq5 computes one dqds transform in ping-pong form. 10// i0 and n0 are zero-indexed. 11// 12// Dlasq5 is an internal routine. It is exported for testing purposes. 13func (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. 16 17 switch { 18 case i0 < 0: 19 panic(i0LT0) 20 case n0 < 0: 21 panic(n0LT0) 22 case len(z) < 4*n0: 23 panic(shortZ) 24 case pp != 0 && pp != 1: 25 panic(badPp) 26 } 27 28 if n0-i0-1 <= 0 { 29 return i0, n0, pp, tau, sigma, dmin, dmin1, dmin2, dn, dnm1, dnm2 30 } 31 32 eps := dlamchP 33 dthresh := eps * (sigma + tau) 34 if tau < dthresh*0.5 { 35 tau = 0 36 } 37 var j4 int 38 var emin float64 39 if tau != 0 { 40 j4 = 4*i0 + pp 41 emin = z[j4+4] 42 d := z[j4] - tau 43 dmin = d 44 // In the reference there are code paths that actually return this value. 45 // dmin1 = -z[j4] 46 if pp == 0 { 47 for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 { 48 j4 := j4loop - 1 49 z[j4-2] = d + z[j4-1] 50 tmp := z[j4+1] / z[j4-2] 51 d = d*tmp - tau 52 dmin = math.Min(dmin, d) 53 z[j4] = z[j4-1] * tmp 54 emin = math.Min(z[j4], emin) 55 } 56 } else { 57 for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 { 58 j4 := j4loop - 1 59 z[j4-3] = d + z[j4] 60 tmp := z[j4+2] / z[j4-3] 61 d = d*tmp - tau 62 dmin = math.Min(dmin, d) 63 z[j4-1] = z[j4] * tmp 64 emin = math.Min(z[j4-1], emin) 65 } 66 } 67 // Unroll the last two steps. 68 dnm2 = d 69 dmin2 = dmin 70 j4 = 4*((n0+1)-2) - pp - 1 71 j4p2 := j4 + 2*pp - 1 72 z[j4-2] = dnm2 + z[j4p2] 73 z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2]) 74 dnm1 = z[j4p2+2]*(dnm2/z[j4-2]) - tau 75 dmin = math.Min(dmin, dnm1) 76 77 dmin1 = dmin 78 j4 += 4 79 j4p2 = j4 + 2*pp - 1 80 z[j4-2] = dnm1 + z[j4p2] 81 z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2]) 82 dn = z[j4p2+2]*(dnm1/z[j4-2]) - tau 83 dmin = math.Min(dmin, dn) 84 } else { 85 // This is the version that sets d's to zero if they are small enough. 86 j4 = 4*(i0+1) + pp - 4 87 emin = z[j4+4] 88 d := z[j4] - tau 89 dmin = d 90 // In the reference there are code paths that actually return this value. 91 // dmin1 = -z[j4] 92 if pp == 0 { 93 for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 { 94 j4 := j4loop - 1 95 z[j4-2] = d + z[j4-1] 96 tmp := z[j4+1] / z[j4-2] 97 d = d*tmp - tau 98 if d < dthresh { 99 d = 0 100 } 101 dmin = math.Min(dmin, d) 102 z[j4] = z[j4-1] * tmp 103 emin = math.Min(z[j4], emin) 104 } 105 } else { 106 for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 { 107 j4 := j4loop - 1 108 z[j4-3] = d + z[j4] 109 tmp := z[j4+2] / z[j4-3] 110 d = d*tmp - tau 111 if d < dthresh { 112 d = 0 113 } 114 dmin = math.Min(dmin, d) 115 z[j4-1] = z[j4] * tmp 116 emin = math.Min(z[j4-1], emin) 117 } 118 } 119 // Unroll the last two steps. 120 dnm2 = d 121 dmin2 = dmin 122 j4 = 4*((n0+1)-2) - pp - 1 123 j4p2 := j4 + 2*pp - 1 124 z[j4-2] = dnm2 + z[j4p2] 125 z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2]) 126 dnm1 = z[j4p2+2]*(dnm2/z[j4-2]) - tau 127 dmin = math.Min(dmin, dnm1) 128 129 dmin1 = dmin 130 j4 += 4 131 j4p2 = j4 + 2*pp - 1 132 z[j4-2] = dnm1 + z[j4p2] 133 z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2]) 134 dn = z[j4p2+2]*(dnm1/z[j4-2]) - tau 135 dmin = math.Min(dmin, dn) 136 } 137 z[j4+2] = dn 138 z[4*(n0+1)-pp-1] = emin 139 return i0, n0, pp, tau, sigma, dmin, dmin1, dmin2, dn, dnm1, dnm2 140} 141