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 }