github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/native/dlasq3.go (about) 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 5 package native 6 7 import "math" 8 9 // Dlasq3 checks for deflation, computes a shift (tau) and calls dqds. 10 // In case of failure it changes shifts, and tries again until output 11 // is positive. 12 // 13 // Dlasq3 is an internal routine. It is exported for testing purposes. 14 func (impl Implementation) Dlasq3(i0, n0 int, z []float64, pp int, dmin, sigma, desig, qmax float64, nFail, iter, nDiv int, ttype int, dmin1, dmin2, dn, dn1, dn2, g, tau float64) ( 15 i0Out, n0Out, ppOut int, dminOut, sigmaOut, desigOut, qmaxOut float64, nFailOut, iterOut, nDivOut, ttypeOut int, dmin1Out, dmin2Out, dnOut, dn1Out, dn2Out, gOut, tauOut float64) { 16 const cbias = 1.5 17 18 n0in := n0 19 eps := dlamchP 20 tol := eps * 100 21 tol2 := tol * tol 22 var nn int 23 var t float64 24 for { 25 if n0 < i0 { 26 return i0, n0, pp, dmin, sigma, desig, qmax, nFail, iter, nDiv, ttype, dmin1, dmin2, dn, dn1, dn2, g, tau 27 } 28 if n0 == i0 { 29 z[4*(n0+1)-4] = z[4*(n0+1)+pp-4] + sigma 30 n0-- 31 continue 32 } 33 nn = 4*(n0+1) + pp - 1 34 if n0 != i0+1 { 35 // Check whether e[n0-1] is negligible, 1 eigenvalue. 36 if z[nn-5] > tol2*(sigma+z[nn-3]) && z[nn-2*pp-4] > tol2*z[nn-7] { 37 // Check whether e[n0-2] is negligible, 2 eigenvalues. 38 if z[nn-9] > tol2*sigma && z[nn-2*pp-8] > tol2*z[nn-11] { 39 break 40 } 41 } else { 42 z[4*(n0+1)-4] = z[4*(n0+1)+pp-4] + sigma 43 n0-- 44 continue 45 } 46 } 47 if z[nn-3] > z[nn-7] { 48 z[nn-3], z[nn-7] = z[nn-7], z[nn-3] 49 } 50 t = 0.5 * (z[nn-7] - z[nn-3] + z[nn-5]) 51 if z[nn-5] > z[nn-3]*tol2 && t != 0 { 52 s := z[nn-3] * (z[nn-5] / t) 53 if s <= t { 54 s = z[nn-3] * (z[nn-5] / (t * (1 + math.Sqrt(1+s/t)))) 55 } else { 56 s = z[nn-3] * (z[nn-5] / (t + math.Sqrt(t)*math.Sqrt(t+s))) 57 } 58 t = z[nn-7] + (s + z[nn-5]) 59 z[nn-3] *= z[nn-7] / t 60 z[nn-7] = t 61 } 62 z[4*(n0+1)-8] = z[nn-7] + sigma 63 z[4*(n0+1)-4] = z[nn-3] + sigma 64 n0 -= 2 65 } 66 if pp == 2 { 67 pp = 0 68 } 69 70 // Reverse the qd-array, if warranted. 71 if dmin <= 0 || n0 < n0in { 72 if cbias*z[4*(i0+1)+pp-4] < z[4*(n0+1)+pp-4] { 73 ipn4Out := 4 * (i0 + n0 + 2) 74 for j4loop := 4 * (i0 + 1); j4loop <= 2*((i0+1)+(n0+1)-1); j4loop += 4 { 75 ipn4 := ipn4Out - 1 76 j4 := j4loop - 1 77 78 z[j4-3], z[ipn4-j4-4] = z[ipn4-j4-4], z[j4-3] 79 z[j4-2], z[ipn4-j4-3] = z[ipn4-j4-3], z[j4-2] 80 z[j4-1], z[ipn4-j4-6] = z[ipn4-j4-6], z[j4-1] 81 z[j4], z[ipn4-j4-5] = z[ipn4-j4-5], z[j4] 82 } 83 if n0-i0 <= 4 { 84 z[4*(n0+1)+pp-2] = z[4*(i0+1)+pp-2] 85 z[4*(n0+1)-pp-1] = z[4*(i0+1)-pp-1] 86 } 87 dmin2 = math.Min(dmin2, z[4*(i0+1)-pp-2]) 88 z[4*(n0+1)+pp-2] = math.Min(math.Min(z[4*(n0+1)+pp-2], z[4*(i0+1)+pp-2]), z[4*(i0+1)+pp+2]) 89 z[4*(n0+1)-pp-1] = math.Min(math.Min(z[4*(n0+1)-pp-1], z[4*(i0+1)-pp-1]), z[4*(i0+1)-pp+3]) 90 qmax = math.Max(math.Max(qmax, z[4*(i0+1)+pp-4]), z[4*(i0+1)+pp]) 91 dmin = math.Copysign(0, -1) // Fortran code has -zero, but -0 in go is 0 92 } 93 } 94 95 // Choose a shift. 96 tau, ttype, g = impl.Dlasq4(i0, n0, z, pp, n0in, dmin, dmin1, dmin2, dn, dn1, dn2, tau, ttype, g) 97 98 // Call dqds until dmin > 0. 99 loop: 100 for { 101 i0, n0, pp, tau, sigma, dmin, dmin1, dmin2, dn, dn1, dn2 = impl.Dlasq5(i0, n0, z, pp, tau, sigma) 102 103 nDiv += n0 - i0 + 2 104 iter++ 105 switch { 106 case dmin >= 0 && dmin1 >= 0: 107 // Success. 108 goto done 109 110 case dmin < 0 && dmin1 > 0 && z[4*n0-pp-1] < tol*(sigma+dn1) && math.Abs(dn) < tol*sigma: 111 // Convergence hidden by negative dn. 112 z[4*n0-pp+1] = 0 113 dmin = 0 114 goto done 115 116 case dmin < 0: 117 // Tau too big. Select new Tau and try again. 118 nFail++ 119 if ttype < -22 { 120 // Failed twice. Play it safe. 121 tau = 0 122 } else if dmin1 > 0 { 123 // Late failure. Gives excellent shift. 124 tau = (tau + dmin) * (1 - 2*eps) 125 ttype -= 11 126 } else { 127 // Early failure. Divide by 4. 128 tau = tau / 4 129 ttype -= 12 130 } 131 132 case math.IsNaN(dmin): 133 if tau == 0 { 134 break loop 135 } 136 tau = 0 137 138 default: 139 // Possible underflow. Play it safe. 140 break loop 141 } 142 } 143 144 // Risk of underflow. 145 dmin, dmin1, dmin2, dn, dn1, dn2 = impl.Dlasq6(i0, n0, z, pp) 146 nDiv += n0 - i0 + 2 147 iter++ 148 tau = 0 149 150 done: 151 if tau < sigma { 152 desig += tau 153 t = sigma + desig 154 desig -= t - sigma 155 } else { 156 t = sigma + tau 157 desig += sigma - (t - tau) 158 } 159 sigma = t 160 return i0, n0, pp, dmin, sigma, desig, qmax, nFail, iter, nDiv, ttype, dmin1, dmin2, dn, dn1, dn2, g, tau 161 }