vendor/gonum.org/v1/gonum/lapack/gonum/dlasq6.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 // Dlasq6 computes one dqd transform in ping-pong form with protection against
  10 // overflow and underflow. z has length at least 4*(n0+1) and holds the qd array.
  11 // i0 is the zero-based first index.
  12 // n0 is the zero-based last index.
  13 //
  14 // Dlasq6 is an internal routine. It is exported for testing purposes.
  15 func (impl Implementation) Dlasq6(i0, n0 int, z []float64, pp int) (dmin, dmin1, dmin2, dn, dnm1, dnm2 float64) {
  16         if len(z) < 4*(n0+1) {
  17                 panic(badZ)
  18         }
  19         if n0-i0-1 <= 0 {
  20                 return dmin, dmin1, dmin2, dn, dnm1, dnm2
  21         }
  22         safmin := dlamchS
  23         j4 := 4*(i0+1) + pp - 4 // -4 rather than -3 for zero indexing
  24         emin := z[j4+4]
  25         d := z[j4]
  26         dmin = d
  27         if pp == 0 {
  28                 for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 {
  29                         j4 := j4loop - 1 // Translate back to zero-indexed.
  30                         z[j4-2] = d + z[j4-1]
  31                         if z[j4-2] == 0 {
  32                                 z[j4] = 0
  33                                 d = z[j4+1]
  34                                 dmin = d
  35                                 emin = 0
  36                         } else if safmin*z[j4+1] < z[j4-2] && safmin*z[j4-2] < z[j4+1] {
  37                                 tmp := z[j4+1] / z[j4-2]
  38                                 z[j4] = z[j4-1] * tmp
  39                                 d *= tmp
  40                         } else {
  41                                 z[j4] = z[j4+1] * (z[j4-1] / z[j4-2])
  42                                 d = z[j4+1] * (d / z[j4-2])
  43                         }
  44                         dmin = math.Min(dmin, d)
  45                         emin = math.Min(emin, z[j4])
  46                 }
  47         } else {
  48                 for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 {
  49                         j4 := j4loop - 1
  50                         z[j4-3] = d + z[j4]
  51                         if z[j4-3] == 0 {
  52                                 z[j4-1] = 0
  53                                 d = z[j4+2]
  54                                 dmin = d
  55                                 emin = 0
  56                         } else if safmin*z[j4+2] < z[j4-3] && safmin*z[j4-3] < z[j4+2] {
  57                                 tmp := z[j4+2] / z[j4-3]
  58                                 z[j4-1] = z[j4] * tmp
  59                                 d *= tmp
  60                         } else {
  61                                 z[j4-1] = z[j4+2] * (z[j4] / z[j4-3])
  62                                 d = z[j4+2] * (d / z[j4-3])
  63                         }
  64                         dmin = math.Min(dmin, d)
  65                         emin = math.Min(emin, z[j4-1])
  66                 }
  67         }
  68         // Unroll last two steps.
  69         dnm2 = d
  70         dmin2 = dmin
  71         j4 = 4*(n0-1) - pp - 1
  72         j4p2 := j4 + 2*pp - 1
  73         z[j4-2] = dnm2 + z[j4p2]
  74         if z[j4-2] == 0 {
  75                 z[j4] = 0
  76                 dnm1 = z[j4p2+2]
  77                 dmin = dnm1
  78                 emin = 0
  79         } else if safmin*z[j4p2+2] < z[j4-2] && safmin*z[j4-2] < z[j4p2+2] {
  80                 tmp := z[j4p2+2] / z[j4-2]
  81                 z[j4] = z[j4p2] * tmp
  82                 dnm1 = dnm2 * tmp
  83         } else {
  84                 z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2])
  85                 dnm1 = z[j4p2+2] * (dnm2 / z[j4-2])
  86         }
  87         dmin = math.Min(dmin, dnm1)
  88         dmin1 = dmin
  89         j4 += 4
  90         j4p2 = j4 + 2*pp - 1
  91         z[j4-2] = dnm1 + z[j4p2]
  92         if z[j4-2] == 0 {
  93                 z[j4] = 0
  94                 dn = z[j4p2+2]
  95                 dmin = dn
  96                 emin = 0
  97         } else if safmin*z[j4p2+2] < z[j4-2] && safmin*z[j4-2] < z[j4p2+2] {
  98                 tmp := z[j4p2+2] / z[j4-2]
  99                 z[j4] = z[j4p2] * tmp
 100                 dn = dnm1 * tmp
 101         } else {
 102                 z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2])
 103                 dn = z[j4p2+2] * (dnm1 / z[j4-2])
 104         }
 105         dmin = math.Min(dmin, dn)
 106         z[j4+2] = dn
 107         z[4*(n0+1)-pp-1] = emin
 108         return dmin, dmin1, dmin2, dn, dnm1, dnm2
 109 }