vendor/gonum.org/v1/gonum/lapack/gonum/dlasq5.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 // 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.
  13 func (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         if n0-i0-1 <= 0 {
  17                 return i0, n0, pp, tau, sigma, dmin, dmin1, dmin2, dn, dnm1, dnm2
  18         }
  19         eps := dlamchP
  20         dthresh := eps * (sigma + tau)
  21         if tau < dthresh*0.5 {
  22                 tau = 0
  23         }
  24         var j4 int
  25         var emin float64
  26         if tau != 0 {
  27                 j4 = 4*i0 + pp
  28                 emin = z[j4+4]
  29                 d := z[j4] - tau
  30                 dmin = d
  31                 // In the reference there are code paths that actually return this value.
  32                 // dmin1 = -z[j4]
  33                 if pp == 0 {
  34                         for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 {
  35                                 j4 := j4loop - 1
  36                                 z[j4-2] = d + z[j4-1]
  37                                 tmp := z[j4+1] / z[j4-2]
  38                                 d = d*tmp - tau
  39                                 dmin = math.Min(dmin, d)
  40                                 z[j4] = z[j4-1] * tmp
  41                                 emin = math.Min(z[j4], emin)
  42                         }
  43                 } else {
  44                         for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 {
  45                                 j4 := j4loop - 1
  46                                 z[j4-3] = d + z[j4]
  47                                 tmp := z[j4+2] / z[j4-3]
  48                                 d = d*tmp - tau
  49                                 dmin = math.Min(dmin, d)
  50                                 z[j4-1] = z[j4] * tmp
  51                                 emin = math.Min(z[j4-1], emin)
  52                         }
  53                 }
  54                 // Unroll the last two steps.
  55                 dnm2 = d
  56                 dmin2 = dmin
  57                 j4 = 4*((n0+1)-2) - pp - 1
  58                 j4p2 := j4 + 2*pp - 1
  59                 z[j4-2] = dnm2 + z[j4p2]
  60                 z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2])
  61                 dnm1 = z[j4p2+2]*(dnm2/z[j4-2]) - tau
  62                 dmin = math.Min(dmin, dnm1)
  63
  64                 dmin1 = dmin
  65                 j4 += 4
  66                 j4p2 = j4 + 2*pp - 1
  67                 z[j4-2] = dnm1 + z[j4p2]
  68                 z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2])
  69                 dn = z[j4p2+2]*(dnm1/z[j4-2]) - tau
  70                 dmin = math.Min(dmin, dn)
  71         } else {
  72                 // This is the version that sets d's to zero if they are small enough.
  73                 j4 = 4*(i0+1) + pp - 4
  74                 emin = z[j4+4]
  75                 d := z[j4] - tau
  76                 dmin = d
  77                 // In the reference there are code paths that actually return this value.
  78                 // dmin1 = -z[j4]
  79                 if pp == 0 {
  80                         for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 {
  81                                 j4 := j4loop - 1
  82                                 z[j4-2] = d + z[j4-1]
  83                                 tmp := z[j4+1] / z[j4-2]
  84                                 d = d*tmp - tau
  85                                 if d < dthresh {
  86                                         d = 0
  87                                 }
  88                                 dmin = math.Min(dmin, d)
  89                                 z[j4] = z[j4-1] * tmp
  90                                 emin = math.Min(z[j4], emin)
  91                         }
  92                 } else {
  93                         for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 {
  94                                 j4 := j4loop - 1
  95                                 z[j4-3] = d + z[j4]
  96                                 tmp := z[j4+2] / z[j4-3]
  97                                 d = d*tmp - tau
  98                                 if d < dthresh {
  99                                         d = 0
 100                                 }
 101                                 dmin = math.Min(dmin, d)
 102                                 z[j4-1] = z[j4] * tmp
 103                                 emin = math.Min(z[j4-1], emin)
 104                         }
 105                 }
 106                 // Unroll the last two steps.
 107                 dnm2 = d
 108                 dmin2 = dmin
 109                 j4 = 4*((n0+1)-2) - pp - 1
 110                 j4p2 := j4 + 2*pp - 1
 111                 z[j4-2] = dnm2 + z[j4p2]
 112                 z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2])
 113                 dnm1 = z[j4p2+2]*(dnm2/z[j4-2]) - tau
 114                 dmin = math.Min(dmin, dnm1)
 115
 116                 dmin1 = dmin
 117                 j4 += 4
 118                 j4p2 = j4 + 2*pp - 1
 119                 z[j4-2] = dnm1 + z[j4p2]
 120                 z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2])
 121                 dn = z[j4p2+2]*(dnm1/z[j4-2]) - tau
 122                 dmin = math.Min(dmin, dn)
 123         }
 124         z[j4+2] = dn
 125         z[4*(n0+1)-pp-1] = emin
 126         return i0, n0, pp, tau, sigma, dmin, dmin1, dmin2, dn, dnm1, dnm2
 127 }