github.com/gopherd/gonum@v0.0.4/lapack/gonum/dlasq6.go (about)

     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  	switch {
    17  	case i0 < 0:
    18  		panic(i0LT0)
    19  	case n0 < 0:
    20  		panic(n0LT0)
    21  	case len(z) < 4*n0:
    22  		panic(shortZ)
    23  	case pp != 0 && pp != 1:
    24  		panic(badPp)
    25  	}
    26  
    27  	if n0-i0-1 <= 0 {
    28  		return dmin, dmin1, dmin2, dn, dnm1, dnm2
    29  	}
    30  
    31  	safmin := dlamchS
    32  	j4 := 4*(i0+1) + pp - 4 // -4 rather than -3 for zero indexing
    33  	emin := z[j4+4]
    34  	d := z[j4]
    35  	dmin = d
    36  	if pp == 0 {
    37  		for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 {
    38  			j4 := j4loop - 1 // Translate back to zero-indexed.
    39  			z[j4-2] = d + z[j4-1]
    40  			if z[j4-2] == 0 {
    41  				z[j4] = 0
    42  				d = z[j4+1]
    43  				dmin = d
    44  				emin = 0
    45  			} else if safmin*z[j4+1] < z[j4-2] && safmin*z[j4-2] < z[j4+1] {
    46  				tmp := z[j4+1] / z[j4-2]
    47  				z[j4] = z[j4-1] * tmp
    48  				d *= tmp
    49  			} else {
    50  				z[j4] = z[j4+1] * (z[j4-1] / z[j4-2])
    51  				d = z[j4+1] * (d / z[j4-2])
    52  			}
    53  			dmin = math.Min(dmin, d)
    54  			emin = math.Min(emin, z[j4])
    55  		}
    56  	} else {
    57  		for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 {
    58  			j4 := j4loop - 1
    59  			z[j4-3] = d + z[j4]
    60  			if z[j4-3] == 0 {
    61  				z[j4-1] = 0
    62  				d = z[j4+2]
    63  				dmin = d
    64  				emin = 0
    65  			} else if safmin*z[j4+2] < z[j4-3] && safmin*z[j4-3] < z[j4+2] {
    66  				tmp := z[j4+2] / z[j4-3]
    67  				z[j4-1] = z[j4] * tmp
    68  				d *= tmp
    69  			} else {
    70  				z[j4-1] = z[j4+2] * (z[j4] / z[j4-3])
    71  				d = z[j4+2] * (d / z[j4-3])
    72  			}
    73  			dmin = math.Min(dmin, d)
    74  			emin = math.Min(emin, z[j4-1])
    75  		}
    76  	}
    77  	// Unroll last two steps.
    78  	dnm2 = d
    79  	dmin2 = dmin
    80  	j4 = 4*(n0-1) - pp - 1
    81  	j4p2 := j4 + 2*pp - 1
    82  	z[j4-2] = dnm2 + z[j4p2]
    83  	if z[j4-2] == 0 {
    84  		z[j4] = 0
    85  		dnm1 = z[j4p2+2]
    86  		dmin = dnm1
    87  		emin = 0
    88  	} else if safmin*z[j4p2+2] < z[j4-2] && safmin*z[j4-2] < z[j4p2+2] {
    89  		tmp := z[j4p2+2] / z[j4-2]
    90  		z[j4] = z[j4p2] * tmp
    91  		dnm1 = dnm2 * tmp
    92  	} else {
    93  		z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2])
    94  		dnm1 = z[j4p2+2] * (dnm2 / z[j4-2])
    95  	}
    96  	dmin = math.Min(dmin, dnm1)
    97  	dmin1 = dmin
    98  	j4 += 4
    99  	j4p2 = j4 + 2*pp - 1
   100  	z[j4-2] = dnm1 + z[j4p2]
   101  	if z[j4-2] == 0 {
   102  		z[j4] = 0
   103  		dn = z[j4p2+2]
   104  		dmin = dn
   105  		emin = 0
   106  	} else if safmin*z[j4p2+2] < z[j4-2] && safmin*z[j4-2] < z[j4p2+2] {
   107  		tmp := z[j4p2+2] / z[j4-2]
   108  		z[j4] = z[j4p2] * tmp
   109  		dn = dnm1 * tmp
   110  	} else {
   111  		z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2])
   112  		dn = z[j4p2+2] * (dnm1 / z[j4-2])
   113  	}
   114  	dmin = math.Min(dmin, dn)
   115  	z[j4+2] = dn
   116  	z[4*(n0+1)-pp-1] = emin
   117  	return dmin, dmin1, dmin2, dn, dnm1, dnm2
   118  }