github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/native/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 native
     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  		dmin1 = -z[j4]
    32  		if pp == 0 {
    33  			for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 {
    34  				j4 := j4loop - 1
    35  				z[j4-2] = d + z[j4-1]
    36  				tmp := z[j4+1] / z[j4-2]
    37  				d = d*tmp - tau
    38  				dmin = math.Min(dmin, d)
    39  				z[j4] = z[j4-1] * tmp
    40  				emin = math.Min(z[j4], emin)
    41  			}
    42  		} else {
    43  			for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 {
    44  				j4 := j4loop - 1
    45  				z[j4-3] = d + z[j4]
    46  				tmp := z[j4+2] / z[j4-3]
    47  				d = d*tmp - tau
    48  				dmin = math.Min(dmin, d)
    49  				z[j4-1] = z[j4] * tmp
    50  				emin = math.Min(z[j4-1], emin)
    51  			}
    52  		}
    53  		// Unroll the last two steps.
    54  		dnm2 = d
    55  		dmin2 = dmin
    56  		j4 = 4*((n0+1)-2) - pp - 1
    57  		j4p2 := j4 + 2*pp - 1
    58  		z[j4-2] = dnm2 + z[j4p2]
    59  		z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2])
    60  		dnm1 = z[j4p2+2]*(dnm2/z[j4-2]) - tau
    61  		dmin = math.Min(dmin, dnm1)
    62  
    63  		dmin1 = dmin
    64  		j4 += 4
    65  		j4p2 = j4 + 2*pp - 1
    66  		z[j4-2] = dnm1 + z[j4p2]
    67  		z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2])
    68  		dn = z[j4p2+2]*(dnm1/z[j4-2]) - tau
    69  		dmin = math.Min(dmin, dn)
    70  	} else {
    71  		// This is the version that sets d's to zero if they are small enough.
    72  		j4 = 4*(i0+1) + pp - 4
    73  		emin = z[j4+4]
    74  		d := z[j4] - tau
    75  		dmin = d
    76  		dmin1 = -z[j4]
    77  		if pp == 0 {
    78  			for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 {
    79  				j4 := j4loop - 1
    80  				z[j4-2] = d + z[j4-1]
    81  				tmp := z[j4+1] / z[j4-2]
    82  				d = d*tmp - tau
    83  				if d < dthresh {
    84  					d = 0
    85  				}
    86  				dmin = math.Min(dmin, d)
    87  				z[j4] = z[j4-1] * tmp
    88  				emin = math.Min(z[j4], emin)
    89  			}
    90  		} else {
    91  			for j4loop := 4 * (i0 + 1); j4loop <= 4*((n0+1)-3); j4loop += 4 {
    92  				j4 := j4loop - 1
    93  				z[j4-3] = d + z[j4]
    94  				tmp := z[j4+2] / z[j4-3]
    95  				d = d*tmp - tau
    96  				if d < dthresh {
    97  					d = 0
    98  				}
    99  				dmin = math.Min(dmin, d)
   100  				z[j4-1] = z[j4] * tmp
   101  				emin = math.Min(z[j4-1], emin)
   102  			}
   103  		}
   104  		// Unroll the last two steps.
   105  		dnm2 = d
   106  		dmin2 = dmin
   107  		j4 = 4*((n0+1)-2) - pp - 1
   108  		j4p2 := j4 + 2*pp - 1
   109  		z[j4-2] = dnm2 + z[j4p2]
   110  		z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2])
   111  		dnm1 = z[j4p2+2]*(dnm2/z[j4-2]) - tau
   112  		dmin = math.Min(dmin, dnm1)
   113  
   114  		dmin1 = dmin
   115  		j4 += 4
   116  		j4p2 = j4 + 2*pp - 1
   117  		z[j4-2] = dnm1 + z[j4p2]
   118  		z[j4] = z[j4p2+2] * (z[j4p2] / z[j4-2])
   119  		dn = z[j4p2+2]*(dnm1/z[j4-2]) - tau
   120  		dmin = math.Min(dmin, dn)
   121  	}
   122  	z[j4+2] = dn
   123  	z[4*(n0+1)-pp-1] = emin
   124  	return i0, n0, pp, tau, sigma, dmin, dmin1, dmin2, dn, dnm1, dnm2
   125  }