github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/lapack/gonum/dlasq3.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  // Dlasq3 checks for deflation, computes a shift (tau) and calls dqds.
    10  // In case of failure it changes shifts, and tries again until output
    11  // is positive.
    12  //
    13  // Dlasq3 is an internal routine. It is exported for testing purposes.
    14  func (impl Implementation) Dlasq3(i0, n0 int, z []float64, pp int, dmin, sigma, desig, qmax float64, nFail, iter, nDiv int, ttype int, dmin1, dmin2, dn, dn1, dn2, g, tau float64) (
    15  	i0Out, n0Out, ppOut int, dminOut, sigmaOut, desigOut, qmaxOut float64, nFailOut, iterOut, nDivOut, ttypeOut int, dmin1Out, dmin2Out, dnOut, dn1Out, dn2Out, gOut, tauOut 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 && pp != 2:
    24  		panic(badPp)
    25  	}
    26  
    27  	const cbias = 1.5
    28  
    29  	n0in := n0
    30  	eps := dlamchP
    31  	tol := eps * 100
    32  	tol2 := tol * tol
    33  	var nn int
    34  	var t float64
    35  	for {
    36  		if n0 < i0 {
    37  			return i0, n0, pp, dmin, sigma, desig, qmax, nFail, iter, nDiv, ttype, dmin1, dmin2, dn, dn1, dn2, g, tau
    38  		}
    39  		if n0 == i0 {
    40  			z[4*(n0+1)-4] = z[4*(n0+1)+pp-4] + sigma
    41  			n0--
    42  			continue
    43  		}
    44  		nn = 4*(n0+1) + pp - 1
    45  		if n0 != i0+1 {
    46  			// Check whether e[n0-1] is negligible, 1 eigenvalue.
    47  			if z[nn-5] > tol2*(sigma+z[nn-3]) && z[nn-2*pp-4] > tol2*z[nn-7] {
    48  				// Check whether e[n0-2] is negligible, 2 eigenvalues.
    49  				if z[nn-9] > tol2*sigma && z[nn-2*pp-8] > tol2*z[nn-11] {
    50  					break
    51  				}
    52  			} else {
    53  				z[4*(n0+1)-4] = z[4*(n0+1)+pp-4] + sigma
    54  				n0--
    55  				continue
    56  			}
    57  		}
    58  		if z[nn-3] > z[nn-7] {
    59  			z[nn-3], z[nn-7] = z[nn-7], z[nn-3]
    60  		}
    61  		t = 0.5 * (z[nn-7] - z[nn-3] + z[nn-5])
    62  		if z[nn-5] > z[nn-3]*tol2 && t != 0 {
    63  			s := z[nn-3] * (z[nn-5] / t)
    64  			if s <= t {
    65  				s = z[nn-3] * (z[nn-5] / (t * (1 + math.Sqrt(1+s/t))))
    66  			} else {
    67  				s = z[nn-3] * (z[nn-5] / (t + math.Sqrt(t)*math.Sqrt(t+s)))
    68  			}
    69  			t = z[nn-7] + (s + z[nn-5])
    70  			z[nn-3] *= z[nn-7] / t
    71  			z[nn-7] = t
    72  		}
    73  		z[4*(n0+1)-8] = z[nn-7] + sigma
    74  		z[4*(n0+1)-4] = z[nn-3] + sigma
    75  		n0 -= 2
    76  	}
    77  	if pp == 2 {
    78  		pp = 0
    79  	}
    80  
    81  	// Reverse the qd-array, if warranted.
    82  	if dmin <= 0 || n0 < n0in {
    83  		if cbias*z[4*(i0+1)+pp-4] < z[4*(n0+1)+pp-4] {
    84  			ipn4Out := 4 * (i0 + n0 + 2)
    85  			for j4loop := 4 * (i0 + 1); j4loop <= 2*((i0+1)+(n0+1)-1); j4loop += 4 {
    86  				ipn4 := ipn4Out - 1
    87  				j4 := j4loop - 1
    88  
    89  				z[j4-3], z[ipn4-j4-4] = z[ipn4-j4-4], z[j4-3]
    90  				z[j4-2], z[ipn4-j4-3] = z[ipn4-j4-3], z[j4-2]
    91  				z[j4-1], z[ipn4-j4-6] = z[ipn4-j4-6], z[j4-1]
    92  				z[j4], z[ipn4-j4-5] = z[ipn4-j4-5], z[j4]
    93  			}
    94  			if n0-i0 <= 4 {
    95  				z[4*(n0+1)+pp-2] = z[4*(i0+1)+pp-2]
    96  				z[4*(n0+1)-pp-1] = z[4*(i0+1)-pp-1]
    97  			}
    98  			dmin2 = math.Min(dmin2, z[4*(i0+1)-pp-2])
    99  			z[4*(n0+1)+pp-2] = math.Min(math.Min(z[4*(n0+1)+pp-2], z[4*(i0+1)+pp-2]), z[4*(i0+1)+pp+2])
   100  			z[4*(n0+1)-pp-1] = math.Min(math.Min(z[4*(n0+1)-pp-1], z[4*(i0+1)-pp-1]), z[4*(i0+1)-pp+3])
   101  			qmax = math.Max(math.Max(qmax, z[4*(i0+1)+pp-4]), z[4*(i0+1)+pp])
   102  			dmin = math.Copysign(0, -1) // Fortran code has -zero, but -0 in go is 0
   103  		}
   104  	}
   105  
   106  	// Choose a shift.
   107  	tau, ttype, g = impl.Dlasq4(i0, n0, z, pp, n0in, dmin, dmin1, dmin2, dn, dn1, dn2, tau, ttype, g)
   108  
   109  	// Call dqds until dmin > 0.
   110  loop:
   111  	for {
   112  		i0, n0, pp, tau, sigma, dmin, dmin1, dmin2, dn, dn1, dn2 = impl.Dlasq5(i0, n0, z, pp, tau, sigma)
   113  
   114  		nDiv += n0 - i0 + 2
   115  		iter++
   116  		switch {
   117  		case dmin >= 0 && dmin1 >= 0:
   118  			// Success.
   119  			goto done
   120  
   121  		case dmin < 0 && dmin1 > 0 && z[4*n0-pp-1] < tol*(sigma+dn1) && math.Abs(dn) < tol*sigma:
   122  			// Convergence hidden by negative dn.
   123  			z[4*n0-pp+1] = 0
   124  			dmin = 0
   125  			goto done
   126  
   127  		case dmin < 0:
   128  			// Tau too big. Select new Tau and try again.
   129  			nFail++
   130  			if ttype < -22 {
   131  				// Failed twice. Play it safe.
   132  				tau = 0
   133  			} else if dmin1 > 0 {
   134  				// Late failure. Gives excellent shift.
   135  				tau = (tau + dmin) * (1 - 2*eps)
   136  				ttype -= 11
   137  			} else {
   138  				// Early failure. Divide by 4.
   139  				tau = tau / 4
   140  				ttype -= 12
   141  			}
   142  
   143  		case math.IsNaN(dmin):
   144  			if tau == 0 {
   145  				break loop
   146  			}
   147  			tau = 0
   148  
   149  		default:
   150  			// Possible underflow. Play it safe.
   151  			break loop
   152  		}
   153  	}
   154  
   155  	// Risk of underflow.
   156  	dmin, dmin1, dmin2, dn, dn1, dn2 = impl.Dlasq6(i0, n0, z, pp)
   157  	nDiv += n0 - i0 + 2
   158  	iter++
   159  	tau = 0
   160  
   161  done:
   162  	if tau < sigma {
   163  		desig += tau
   164  		t = sigma + desig
   165  		desig -= t - sigma
   166  	} else {
   167  		t = sigma + tau
   168  		desig += sigma - (t - tau)
   169  	}
   170  	sigma = t
   171  	return i0, n0, pp, dmin, sigma, desig, qmax, nFail, iter, nDiv, ttype, dmin1, dmin2, dn, dn1, dn2, g, tau
   172  }