vendor/gonum.org/v1/gonum/lapack/gonum/dlasq4.go

   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 gonum
   6
   7 import "math"
   8
   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.
  12 //
  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) {
  15         const (
  16                 cnst1 = 0.563
  17                 cnst2 = 1.01
  18                 cnst3 = 1.05
  19
  20                 cnstthird = 0.333 // TODO(btracey): Fix?
  21         )
  22         // A negative dmin forces the shift to take that absolute value
  23         // ttype records the type of shift.
  24         if dmin <= 0 {
  25                 tau = -dmin
  26                 ttype = -1
  27                 return tau, ttype, g
  28         }
  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
  31         if n0in == n0 {
  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
  39                                 var gap1 float64
  40                                 if gap2 > 0 && gap2 > b2 {
  41                                         gap1 = a2 - dn - (b2/gap2)*b2
  42                                 } else {
  43                                         gap1 = a2 - dn - (b1 + b2)
  44                                 }
  45                                 if gap1 > 0 && gap1 > b1 {
  46                                         s = math.Max(dn-(b1/gap1)*b1, 0.5*dmin)
  47                                         ttype = -2
  48                                 } else {
  49                                         s = 0
  50                                         if dn > b1 {
  51                                                 s = dn - b1
  52                                         }
  53                                         if a2 > b1+b2 {
  54                                                 s = math.Min(s, a2-(b1+b2))
  55                                         }
  56                                         s = math.Max(s, cnstthird*dmin)
  57                                         ttype = -3
  58                                 }
  59                         } else {
  60                                 ttype = -4
  61                                 s = dmin / 4
  62                                 var gam float64
  63                                 var np int
  64                                 if dmin == dn {
  65                                         gam = dn
  66                                         a2 = 0
  67                                         if z[nn-5] > z[nn-7] {
  68                                                 return tau, ttype, g
  69                                         }
  70                                         b2 = z[nn-5] / z[nn-7]
  71                                         np = nn - 9
  72                                 } else {
  73                                         np = nn - 2*pp
  74                                         gam = dn1
  75                                         if z[np-4] > z[np-2] {
  76                                                 return tau, ttype, g
  77                                         }
  78                                         a2 = z[np-4] / z[np-2]
  79                                         if z[nn-9] > z[nn-11] {
  80                                                 return tau, ttype, g
  81                                         }
  82                                         b2 = z[nn-9] / z[nn-11]
  83                                         np = nn - 13
  84                                 }
  85                                 // Approximate contribution to norm squared from i < nn-1.
  86                                 a2 += b2
  87                                 for i4loop := np + 1; i4loop >= 4*(i0+1)-1+pp; i4loop -= 4 {
  88                                         i4 := i4loop - 1
  89                                         if b2 == 0 {
  90                                                 break
  91                                         }
  92                                         b1 = b2
  93                                         if z[i4] > z[i4-2] {
  94                                                 return tau, ttype, g
  95                                         }
  96                                         b2 *= z[i4] / z[i4-2]
  97                                         a2 += b2
  98                                         if 100*math.Max(b2, b1) < a2 || cnst1 < a2 {
  99                                                 break
 100                                         }
 101                                 }
 102                                 a2 *= cnst3
 103                                 // Rayleigh quotient residual bound.
 104                                 if a2 < cnst1 {
 105                                         s = gam * (1 - math.Sqrt(a2)) / (1 + a2)
 106                                 }
 107                         }
 108                 } else if dmin == dn2 {
 109                         ttype = -5
 110                         s = dmin / 4
 111                         // Compute contribution to norm squared from i > nn-2.
 112                         np := nn - 2*pp
 113                         b1 := z[np-2]
 114                         b2 := z[np-6]
 115                         gam := dn2
 116                         if z[np-8] > b2 || z[np-4] > b1 {
 117                                 return tau, ttype, g
 118                         }
 119                         a2 := (z[np-8] / b2) * (1 + z[np-4]/b1)
 120                         // Approximate contribution to norm squared from i < nn-2.
 121                         if n0-i0 > 2 {
 122                                 b2 = z[nn-13] / z[nn-15]
 123                                 a2 += b2
 124                                 for i4loop := (nn + 1) - 17; i4loop >= 4*(i0+1)-1+pp; i4loop -= 4 {
 125                                         i4 := i4loop - 1
 126                                         if b2 == 0 {
 127                                                 break
 128                                         }
 129                                         b1 = b2
 130                                         if z[i4] > z[i4-2] {
 131                                                 return tau, ttype, g
 132                                         }
 133                                         b2 *= z[i4] / z[i4-2]
 134                                         a2 += b2
 135                                         if 100*math.Max(b2, b1) < a2 || cnst1 < a2 {
 136                                                 break
 137                                         }
 138                                 }
 139                                 a2 *= cnst3
 140                         }
 141                         if a2 < cnst1 {
 142                                 s = gam * (1 - math.Sqrt(a2)) / (1 + a2)
 143                         }
 144                 } else {
 145                         // Case 6, no information to guide us.
 146                         if ttype == -6 {
 147                                 g += cnstthird * (1 - g)
 148                         } else if ttype == -18 {
 149                                 g = cnstthird / 4
 150                         } else {
 151                                 g = 1.0 / 4
 152                         }
 153                         s = g * dmin
 154                         ttype = -6
 155                 }
 156         } else if n0in == (n0 + 1) {
 157                 // One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN.
 158                 if dmin1 == dn1 && dmin2 == dn2 {
 159                         ttype = -7
 160                         s = cnstthird * dmin1
 161                         if z[nn-5] > z[nn-7] {
 162                                 return tau, ttype, g
 163                         }
 164                         b1 := z[nn-5] / z[nn-7]
 165                         b2 := b1
 166                         if b2 != 0 {
 167                                 for i4loop := 4*(n0+1) - 9 + pp; i4loop >= 4*(i0+1)-1+pp; i4loop -= 4 {
 168                                         i4 := i4loop - 1
 169                                         a2 := b1
 170                                         if z[i4] > z[i4-2] {
 171                                                 return tau, ttype, g
 172                                         }
 173                                         b1 *= z[i4] / z[i4-2]
 174                                         b2 += b1
 175                                         if 100*math.Max(b1, a2) < b2 {
 176                                                 break
 177                                         }
 178                                 }
 179                         }
 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))
 185                         } else {
 186                                 s = math.Max(s, a2*(1-cnst2*b2))
 187                                 ttype = -8
 188                         }
 189                 } else {
 190                         s = dmin1 / 4
 191                         if dmin1 == dn1 {
 192                                 s = 0.5 * dmin1
 193                         }
 194                         ttype = -9
 195                 }
 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] {
 199                         ttype = -10
 200                         s = cnstthird * dmin2
 201                         if z[nn-5] > z[nn-7] {
 202                                 return tau, ttype, g
 203                         }
 204                         b1 := z[nn-5] / z[nn-7]
 205                         b2 := b1
 206                         if b2 != 0 {
 207                                 for i4loop := 4*(n0+1) - 9 + pp; i4loop >= 4*(i0+1)-1+pp; i4loop -= 4 {
 208                                         i4 := i4loop - 1
 209                                         if z[i4] > z[i4-2] {
 210                                                 return tau, ttype, g
 211                                         }
 212                                         b1 *= z[i4] / z[i4-2]
 213                                         b2 += b1
 214                                         if 100*b1 < b2 {
 215                                                 break
 216                                         }
 217                                 }
 218                         }
 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))
 224                         } else {
 225                                 s = math.Max(s, a2*(1-cnst2*b2))
 226                         }
 227                 } else {
 228                         s = dmin2 / 4
 229                         ttype = -11
 230                 }
 231         } else if n0in > n0+2 {
 232                 // Case 12, more than two eigenvalues deflated. No information.
 233                 s = 0
 234                 ttype = -12
 235         }
 236         tau = s
 237         return tau, ttype, g
 238 }