blob: dc1e178dfa19afd703618768cc2bdd651e2641f9 [file] [log] [blame]
// 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)
}
// Quick return if possible.
if n == 0 {
return true
}
switch {
case len(d) < n:
panic(shortD)
case len(e) < n-1:
panic(shortE)
}
if n == 1 {
return true
}
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
}