// Copyright ©2016 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"

	"gonum.org/v1/gonum/lapack"
)

// Dsterf computes all eigenvalues of a symmetric tridiagonal matrix using the
// Pal-Walker-Kahan variant of the QL or QR algorithm.
//
// d contains the diagonal elements of the tridiagonal matrix on entry, and
// contains the eigenvalues in ascending order on exit. d must have length at
// least n, or Dsterf will panic.
//
// e contains the off-diagonal elements of the tridiagonal matrix on entry, and is
// overwritten during the call to Dsterf. e must have length of at least n-1 or
// Dsterf will panic.
//
// Dsterf is an internal routine. It is exported for testing purposes.
func (impl Implementation) Dsterf(n int, d, e []float64) (ok bool) {
	if n < 0 {
		panic(nLT0)
	}
	if n == 0 {
		return true
	}
	if len(d) < n {
		panic(badD)
	}
	if len(e) < n-1 {
		panic(badE)
	}

	const (
		none = 0 // The values are not scaled.
		down = 1 // The values are scaled below ssfmax threshold.
		up   = 2 // The values are scaled below ssfmin threshold.
	)

	// Determine the unit roundoff for this environment.
	eps := dlamchE
	eps2 := eps * eps
	safmin := dlamchS
	safmax := 1 / safmin
	ssfmax := math.Sqrt(safmax) / 3
	ssfmin := math.Sqrt(safmin) / eps2

	// Compute the eigenvalues of the tridiagonal matrix.
	maxit := 30
	nmaxit := n * maxit
	jtot := 0

	l1 := 0

	for {
		if l1 > n-1 {
			impl.Dlasrt(lapack.SortIncreasing, n, d)
			return true
		}
		if l1 > 0 {
			e[l1-1] = 0
		}
		var m int
		for m = l1; m < n-1; m++ {
			if math.Abs(e[m]) <= math.Sqrt(math.Abs(d[m]))*math.Sqrt(math.Abs(d[m+1]))*eps {
				e[m] = 0
				break
			}
		}

		l := l1
		lsv := l
		lend := m
		lendsv := lend
		l1 = m + 1
		if lend == 0 {
			continue
		}

		// Scale submatrix in rows and columns l to lend.
		anorm := impl.Dlanst(lapack.MaxAbs, lend-l+1, d[l:], e[l:])
		iscale := none
		if anorm == 0 {
			continue
		}
		if anorm > ssfmax {
			iscale = down
			impl.Dlascl(lapack.General, 0, 0, anorm, ssfmax, lend-l+1, 1, d[l:], n)
			impl.Dlascl(lapack.General, 0, 0, anorm, ssfmax, lend-l, 1, e[l:], n)
		} else if anorm < ssfmin {
			iscale = up
			impl.Dlascl(lapack.General, 0, 0, anorm, ssfmin, lend-l+1, 1, d[l:], n)
			impl.Dlascl(lapack.General, 0, 0, anorm, ssfmin, lend-l, 1, e[l:], n)
		}

		el := e[l:lend]
		for i, v := range el {
			el[i] *= v
		}

		// Choose between QL and QR iteration.
		if math.Abs(d[lend]) < math.Abs(d[l]) {
			lend = lsv
			l = lendsv
		}
		if lend >= l {
			// QL Iteration.
			// Look for small sub-diagonal element.
			for {
				if l != lend {
					for m = l; m < lend; m++ {
						if math.Abs(e[m]) <= eps2*(math.Abs(d[m]*d[m+1])) {
							break
						}
					}
				} else {
					m = lend
				}
				if m < lend {
					e[m] = 0
				}
				p := d[l]
				if m == l {
					// Eigenvalue found.
					l++
					if l > lend {
						break
					}
					continue
				}
				// If remaining matrix is 2 by 2, use Dlae2 to compute its eigenvalues.
				if m == l+1 {
					d[l], d[l+1] = impl.Dlae2(d[l], math.Sqrt(e[l]), d[l+1])
					e[l] = 0
					l += 2
					if l > lend {
						break
					}
					continue
				}
				if jtot == nmaxit {
					break
				}
				jtot++

				// Form shift.
				rte := math.Sqrt(e[l])
				sigma := (d[l+1] - p) / (2 * rte)
				r := impl.Dlapy2(sigma, 1)
				sigma = p - (rte / (sigma + math.Copysign(r, sigma)))

				c := 1.0
				s := 0.0
				gamma := d[m] - sigma
				p = gamma * gamma

				// Inner loop.
				for i := m - 1; i >= l; i-- {
					bb := e[i]
					r := p + bb
					if i != m-1 {
						e[i+1] = s * r
					}
					oldc := c
					c = p / r
					s = bb / r
					oldgam := gamma
					alpha := d[i]
					gamma = c*(alpha-sigma) - s*oldgam
					d[i+1] = oldgam + (alpha - gamma)
					if c != 0 {
						p = (gamma * gamma) / c
					} else {
						p = oldc * bb
					}
				}
				e[l] = s * p
				d[l] = sigma + gamma
			}
		} else {
			for {
				// QR Iteration.
				// Look for small super-diagonal element.
				for m = l; m > lend; m-- {
					if math.Abs(e[m-1]) <= eps2*math.Abs(d[m]*d[m-1]) {
						break
					}
				}
				if m > lend {
					e[m-1] = 0
				}
				p := d[l]
				if m == l {
					// Eigenvalue found.
					l--
					if l < lend {
						break
					}
					continue
				}

				// If remaining matrix is 2 by 2, use Dlae2 to compute its eigenvalues.
				if m == l-1 {
					d[l], d[l-1] = impl.Dlae2(d[l], math.Sqrt(e[l-1]), d[l-1])
					e[l-1] = 0
					l -= 2
					if l < lend {
						break
					}
					continue
				}
				if jtot == nmaxit {
					break
				}
				jtot++

				// Form shift.
				rte := math.Sqrt(e[l-1])
				sigma := (d[l-1] - p) / (2 * rte)
				r := impl.Dlapy2(sigma, 1)
				sigma = p - (rte / (sigma + math.Copysign(r, sigma)))

				c := 1.0
				s := 0.0
				gamma := d[m] - sigma
				p = gamma * gamma

				// Inner loop.
				for i := m; i < l; i++ {
					bb := e[i]
					r := p + bb
					if i != m {
						e[i-1] = s * r
					}
					oldc := c
					c = p / r
					s = bb / r
					oldgam := gamma
					alpha := d[i+1]
					gamma = c*(alpha-sigma) - s*oldgam
					d[i] = oldgam + alpha - gamma
					if c != 0 {
						p = (gamma * gamma) / c
					} else {
						p = oldc * bb
					}
				}
				e[l-1] = s * p
				d[l] = sigma + gamma
			}
		}

		// Undo scaling if necessary
		switch iscale {
		case down:
			impl.Dlascl(lapack.General, 0, 0, ssfmax, anorm, lendsv-lsv+1, 1, d[lsv:], n)
		case up:
			impl.Dlascl(lapack.General, 0, 0, ssfmin, anorm, lendsv-lsv+1, 1, d[lsv:], n)
		}

		// Check for no convergence to an eigenvalue after a total of n*maxit iterations.
		if jtot >= nmaxit {
			break
		}
	}
	for _, v := range e[:n-1] {
		if v != 0 {
			return false
		}
	}
	impl.Dlasrt(lapack.SortIncreasing, n, d)
	return true
}