+++ /dev/null
-// Copyright ©2015 The Gonum Authors. All rights reserved.
-// Use of this source code is governed by a BSD-style
-// license that can be found in the LICENSE file.
-
-package gonum
-
-import "math"
-
-// Dlasq3 checks for deflation, computes a shift (tau) and calls dqds.
-// In case of failure it changes shifts, and tries again until output
-// is positive.
-//
-// Dlasq3 is an internal routine. It is exported for testing purposes.
-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) (
- i0Out, n0Out, ppOut int, dminOut, sigmaOut, desigOut, qmaxOut float64, nFailOut, iterOut, nDivOut, ttypeOut int, dmin1Out, dmin2Out, dnOut, dn1Out, dn2Out, gOut, tauOut float64) {
- const cbias = 1.5
-
- n0in := n0
- eps := dlamchP
- tol := eps * 100
- tol2 := tol * tol
- var nn int
- var t float64
- for {
- if n0 < i0 {
- return i0, n0, pp, dmin, sigma, desig, qmax, nFail, iter, nDiv, ttype, dmin1, dmin2, dn, dn1, dn2, g, tau
- }
- if n0 == i0 {
- z[4*(n0+1)-4] = z[4*(n0+1)+pp-4] + sigma
- n0--
- continue
- }
- nn = 4*(n0+1) + pp - 1
- if n0 != i0+1 {
- // Check whether e[n0-1] is negligible, 1 eigenvalue.
- if z[nn-5] > tol2*(sigma+z[nn-3]) && z[nn-2*pp-4] > tol2*z[nn-7] {
- // Check whether e[n0-2] is negligible, 2 eigenvalues.
- if z[nn-9] > tol2*sigma && z[nn-2*pp-8] > tol2*z[nn-11] {
- break
- }
- } else {
- z[4*(n0+1)-4] = z[4*(n0+1)+pp-4] + sigma
- n0--
- continue
- }
- }
- if z[nn-3] > z[nn-7] {
- z[nn-3], z[nn-7] = z[nn-7], z[nn-3]
- }
- t = 0.5 * (z[nn-7] - z[nn-3] + z[nn-5])
- if z[nn-5] > z[nn-3]*tol2 && t != 0 {
- s := z[nn-3] * (z[nn-5] / t)
- if s <= t {
- s = z[nn-3] * (z[nn-5] / (t * (1 + math.Sqrt(1+s/t))))
- } else {
- s = z[nn-3] * (z[nn-5] / (t + math.Sqrt(t)*math.Sqrt(t+s)))
- }
- t = z[nn-7] + (s + z[nn-5])
- z[nn-3] *= z[nn-7] / t
- z[nn-7] = t
- }
- z[4*(n0+1)-8] = z[nn-7] + sigma
- z[4*(n0+1)-4] = z[nn-3] + sigma
- n0 -= 2
- }
- if pp == 2 {
- pp = 0
- }
-
- // Reverse the qd-array, if warranted.
- if dmin <= 0 || n0 < n0in {
- if cbias*z[4*(i0+1)+pp-4] < z[4*(n0+1)+pp-4] {
- ipn4Out := 4 * (i0 + n0 + 2)
- for j4loop := 4 * (i0 + 1); j4loop <= 2*((i0+1)+(n0+1)-1); j4loop += 4 {
- ipn4 := ipn4Out - 1
- j4 := j4loop - 1
-
- z[j4-3], z[ipn4-j4-4] = z[ipn4-j4-4], z[j4-3]
- z[j4-2], z[ipn4-j4-3] = z[ipn4-j4-3], z[j4-2]
- z[j4-1], z[ipn4-j4-6] = z[ipn4-j4-6], z[j4-1]
- z[j4], z[ipn4-j4-5] = z[ipn4-j4-5], z[j4]
- }
- if n0-i0 <= 4 {
- z[4*(n0+1)+pp-2] = z[4*(i0+1)+pp-2]
- z[4*(n0+1)-pp-1] = z[4*(i0+1)-pp-1]
- }
- dmin2 = math.Min(dmin2, z[4*(i0+1)-pp-2])
- 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])
- 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])
- qmax = math.Max(math.Max(qmax, z[4*(i0+1)+pp-4]), z[4*(i0+1)+pp])
- dmin = math.Copysign(0, -1) // Fortran code has -zero, but -0 in go is 0
- }
- }
-
- // Choose a shift.
- tau, ttype, g = impl.Dlasq4(i0, n0, z, pp, n0in, dmin, dmin1, dmin2, dn, dn1, dn2, tau, ttype, g)
-
- // Call dqds until dmin > 0.
-loop:
- for {
- i0, n0, pp, tau, sigma, dmin, dmin1, dmin2, dn, dn1, dn2 = impl.Dlasq5(i0, n0, z, pp, tau, sigma)
-
- nDiv += n0 - i0 + 2
- iter++
- switch {
- case dmin >= 0 && dmin1 >= 0:
- // Success.
- goto done
-
- case dmin < 0 && dmin1 > 0 && z[4*n0-pp-1] < tol*(sigma+dn1) && math.Abs(dn) < tol*sigma:
- // Convergence hidden by negative dn.
- z[4*n0-pp+1] = 0
- dmin = 0
- goto done
-
- case dmin < 0:
- // Tau too big. Select new Tau and try again.
- nFail++
- if ttype < -22 {
- // Failed twice. Play it safe.
- tau = 0
- } else if dmin1 > 0 {
- // Late failure. Gives excellent shift.
- tau = (tau + dmin) * (1 - 2*eps)
- ttype -= 11
- } else {
- // Early failure. Divide by 4.
- tau = tau / 4
- ttype -= 12
- }
-
- case math.IsNaN(dmin):
- if tau == 0 {
- break loop
- }
- tau = 0
-
- default:
- // Possible underflow. Play it safe.
- break loop
- }
- }
-
- // Risk of underflow.
- dmin, dmin1, dmin2, dn, dn1, dn2 = impl.Dlasq6(i0, n0, z, pp)
- nDiv += n0 - i0 + 2
- iter++
- tau = 0
-
-done:
- if tau < sigma {
- desig += tau
- t = sigma + desig
- desig -= t - sigma
- } else {
- t = sigma + tau
- desig += sigma - (t - tau)
- }
- sigma = t
- return i0, n0, pp, dmin, sigma, desig, qmax, nFail, iter, nDiv, ttype, dmin1, dmin2, dn, dn1, dn2, g, tau
-}
OSDN Git Service
