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 // Dlasq4 computes an approximation to the smallest eigenvalue using values of d
10 // from the previous transform.
11 // i0, n0, and n0in are zero-indexed.
13 // Dlasq4 is an internal routine. It is exported for testing purposes.
14 func (impl Implementation) Dlasq4(i0, n0 int, z []float64, pp int, n0in int, dmin, dmin1, dmin2, dn, dn1, dn2, tau float64, ttype int, g float64) (tauOut float64, ttypeOut int, gOut float64) {
20 cnstthird = 0.333 // TODO(btracey): Fix?
22 // A negative dmin forces the shift to take that absolute value
23 // ttype records the type of shift.
29 nn := 4*(n0+1) + pp - 1 // -1 for zero indexing
30 s := math.NaN() // Poison s so that failure to take a path below is obvious
32 // No eigenvalues deflated.
33 if dmin == dn || dmin == dn1 {
34 b1 := math.Sqrt(z[nn-3]) * math.Sqrt(z[nn-5])
35 b2 := math.Sqrt(z[nn-7]) * math.Sqrt(z[nn-9])
36 a2 := z[nn-7] + z[nn-5]
37 if dmin == dn && dmin1 == dn1 {
38 gap2 := dmin2 - a2 - dmin2/4
40 if gap2 > 0 && gap2 > b2 {
41 gap1 = a2 - dn - (b2/gap2)*b2
43 gap1 = a2 - dn - (b1 + b2)
45 if gap1 > 0 && gap1 > b1 {
46 s = math.Max(dn-(b1/gap1)*b1, 0.5*dmin)
54 s = math.Min(s, a2-(b1+b2))
56 s = math.Max(s, cnstthird*dmin)
67 if z[nn-5] > z[nn-7] {
70 b2 = z[nn-5] / z[nn-7]
75 if z[np-4] > z[np-2] {
78 a2 = z[np-4] / z[np-2]
79 if z[nn-9] > z[nn-11] {
82 b2 = z[nn-9] / z[nn-11]
85 // Approximate contribution to norm squared from i < nn-1.
87 for i4loop := np + 1; i4loop >= 4*(i0+1)-1+pp; i4loop -= 4 {
98 if 100*math.Max(b2, b1) < a2 || cnst1 < a2 {
103 // Rayleigh quotient residual bound.
105 s = gam * (1 - math.Sqrt(a2)) / (1 + a2)
108 } else if dmin == dn2 {
111 // Compute contribution to norm squared from i > nn-2.
116 if z[np-8] > b2 || z[np-4] > b1 {
119 a2 := (z[np-8] / b2) * (1 + z[np-4]/b1)
120 // Approximate contribution to norm squared from i < nn-2.
122 b2 = z[nn-13] / z[nn-15]
124 for i4loop := (nn + 1) - 17; i4loop >= 4*(i0+1)-1+pp; i4loop -= 4 {
133 b2 *= z[i4] / z[i4-2]
135 if 100*math.Max(b2, b1) < a2 || cnst1 < a2 {
142 s = gam * (1 - math.Sqrt(a2)) / (1 + a2)
145 // Case 6, no information to guide us.
147 g += cnstthird * (1 - g)
148 } else if ttype == -18 {
156 } else if n0in == (n0 + 1) {
157 // One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN.
158 if dmin1 == dn1 && dmin2 == dn2 {
160 s = cnstthird * dmin1
161 if z[nn-5] > z[nn-7] {
164 b1 := z[nn-5] / z[nn-7]
167 for i4loop := 4*(n0+1) - 9 + pp; i4loop >= 4*(i0+1)-1+pp; i4loop -= 4 {
173 b1 *= z[i4] / z[i4-2]
175 if 100*math.Max(b1, a2) < b2 {
180 b2 = math.Sqrt(cnst3 * b2)
181 a2 := dmin1 / (1 + b2*b2)
182 gap2 := 0.5*dmin2 - a2
183 if gap2 > 0 && gap2 > b2*a2 {
184 s = math.Max(s, a2*(1-cnst2*a2*(b2/gap2)*b2))
186 s = math.Max(s, a2*(1-cnst2*b2))
196 } else if n0in == (n0 + 2) {
197 // Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN.
198 if dmin2 == dn2 && 2*z[nn-5] < z[nn-7] {
200 s = cnstthird * dmin2
201 if z[nn-5] > z[nn-7] {
204 b1 := z[nn-5] / z[nn-7]
207 for i4loop := 4*(n0+1) - 9 + pp; i4loop >= 4*(i0+1)-1+pp; i4loop -= 4 {
212 b1 *= z[i4] / z[i4-2]
219 b2 = math.Sqrt(cnst3 * b2)
220 a2 := dmin2 / (1 + b2*b2)
221 gap2 := z[nn-7] + z[nn-9] - math.Sqrt(z[nn-11])*math.Sqrt(z[nn-9]) - a2
222 if gap2 > 0 && gap2 > b2*a2 {
223 s = math.Max(s, a2*(1-cnst2*a2*(b2/gap2)*b2))
225 s = math.Max(s, a2*(1-cnst2*b2))
231 } else if n0in > n0+2 {
232 // Case 12, more than two eigenvalues deflated. No information.