vendor/gonum.org/v1/gonum/lapack/gonum/dlasq3.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 // 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 }