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