blob: d6c7861ab5bb0bba3533455a1b20c067f1ea67df [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/blas" "gonum.org/v1/gonum/blas/blas64" "gonum.org/v1/gonum/lapack" ) // Dsteqr computes the eigenvalues and optionally the eigenvectors of a symmetric // tridiagonal matrix using the implicit QL or QR method. The eigenvectors of a // full or band symmetric matrix can also be found if Dsytrd, Dsptrd, or Dsbtrd // have been used to reduce this matrix to tridiagonal form. // // d, on entry, contains the diagonal elements of the tridiagonal matrix. On exit, // d contains the eigenvalues in ascending order. d must have length n and // Dsteqr will panic otherwise. // // e, on entry, contains the off-diagonal elements of the tridiagonal matrix on // entry, and is overwritten during the call to Dsteqr. e must have length n-1 and // Dsteqr will panic otherwise. // // z, on entry, contains the n×n orthogonal matrix used in the reduction to // tridiagonal form if compz == lapack.EVOrig. On exit, if // compz == lapack.EVOrig, z contains the orthonormal eigenvectors of the // original symmetric matrix, and if compz == lapack.EVTridiag, z contains the // orthonormal eigenvectors of the symmetric tridiagonal matrix. z is not used // if compz == lapack.EVCompNone. // // work must have length at least max(1, 2*n-2) if the eigenvectors are computed, // and Dsteqr will panic otherwise. // // Dsteqr is an internal routine. It is exported for testing purposes. func (impl Implementation) Dsteqr(compz lapack.EVComp, n int, d, e, z []float64, ldz int, work []float64) (ok bool) { switch { case compz != lapack.EVCompNone && compz != lapack.EVTridiag && compz != lapack.EVOrig: panic(badEVComp) case n < 0: panic(nLT0) case ldz < 1, compz != lapack.EVCompNone && ldz < n: panic(badLdZ) } // Quick return if possible. if n == 0 { return true } switch { case len(d) < n: panic(shortD) case len(e) < n-1: panic(shortE) case compz != lapack.EVCompNone && len(z) < (n-1)*ldz+n: panic(shortZ) case compz != lapack.EVCompNone && len(work) < max(1, 2*n-2): panic(shortWork) } var icompz int if compz == lapack.EVOrig { icompz = 1 } else if compz == lapack.EVTridiag { icompz = 2 } if n == 1 { if icompz == 2 { z[0] = 1 } return true } bi := blas64.Implementation() eps := dlamchE eps2 := eps * eps safmin := dlamchS safmax := 1 / safmin ssfmax := math.Sqrt(safmax) / 3 ssfmin := math.Sqrt(safmin) / eps2 // Compute the eigenvalues and eigenvectors of the tridiagonal matrix. if icompz == 2 { impl.Dlaset(blas.All, n, n, 0, 1, z, ldz) } const maxit = 30 nmaxit := n * maxit jtot := 0 // Determine where the matrix splits and choose QL or QR iteration for each // block, according to whether top or bottom diagonal element is smaller. l1 := 0 nm1 := n - 1 type scaletype int const ( down scaletype = iota + 1 up ) var iscale scaletype for { if l1 > n-1 { // Order eigenvalues and eigenvectors. if icompz == 0 { impl.Dlasrt(lapack.SortIncreasing, n, d) } else { // TODO(btracey): Consider replacing this sort with a call to sort.Sort. for ii := 1; ii < n; ii++ { i := ii - 1 k := i p := d[i] for j := ii; j < n; j++ { if d[j] < p { k = j p = d[j] } } if k != i { d[k] = d[i] d[i] = p bi.Dswap(n, z[i:], ldz, z[k:], ldz) } } } return true } if l1 > 0 { e[l1-1] = 0 } var m int if l1 <= nm1 { for m = l1; m < nm1; m++ { test := math.Abs(e[m]) if test == 0 { break } if test <= (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 == l { continue } // Scale submatrix in rows and columns L to Lend anorm := impl.Dlanst(lapack.MaxAbs, lend-l+1, d[l:], e[l:]) switch { case anorm == 0: continue case anorm > ssfmax: iscale = down // Pretend that d and e are matrices with 1 column. impl.Dlascl(lapack.General, 0, 0, anorm, ssfmax, lend-l+1, 1, d[l:], 1) impl.Dlascl(lapack.General, 0, 0, anorm, ssfmax, lend-l, 1, e[l:], 1) case anorm < ssfmin: iscale = up impl.Dlascl(lapack.General, 0, 0, anorm, ssfmin, lend-l+1, 1, d[l:], 1) impl.Dlascl(lapack.General, 0, 0, anorm, ssfmin, lend-l, 1, e[l:], 1) } // Choose between QL and QR. if math.Abs(d[lend]) < math.Abs(d[l]) { lend = lsv l = lendsv } if lend > l { // QL Iteration. Look for small subdiagonal element. for { if l != lend { for m = l; m < lend; m++ { v := math.Abs(e[m]) if v*v <= (eps2*math.Abs(d[m]))*math.Abs(d[m+1])+safmin { 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×2, use Dlae2 to compute its eigensystem. if m == l+1 { if icompz > 0 { d[l], d[l+1], work[l], work[n-1+l] = impl.Dlaev2(d[l], e[l], d[l+1]) impl.Dlasr(blas.Right, lapack.Variable, lapack.Backward, n, 2, work[l:], work[n-1+l:], z[l:], ldz) } else { d[l], d[l+1] = impl.Dlae2(d[l], e[l], d[l+1]) } e[l] = 0 l += 2 if l > lend { break } continue } if jtot == nmaxit { break } jtot++ // Form shift g := (d[l+1] - p) / (2 * e[l]) r := impl.Dlapy2(g, 1) g = d[m] - p + e[l]/(g+math.Copysign(r, g)) s := 1.0 c := 1.0 p = 0.0 // Inner loop for i := m - 1; i >= l; i-- { f := s * e[i] b := c * e[i] c, s, r = impl.Dlartg(g, f) if i != m-1 { e[i+1] = r } g = d[i+1] - p r = (d[i]-g)*s + 2*c*b p = s * r d[i+1] = g + p g = c*r - b // If eigenvectors are desired, then save rotations. if icompz > 0 { work[i] = c work[n-1+i] = -s } } // If eigenvectors are desired, then apply saved rotations. if icompz > 0 { mm := m - l + 1 impl.Dlasr(blas.Right, lapack.Variable, lapack.Backward, n, mm, work[l:], work[n-1+l:], z[l:], ldz) } d[l] -= p e[l] = g } } else { // QR Iteration. // Look for small superdiagonal element. for { if l != lend { for m = l; m > lend; m-- { v := math.Abs(e[m-1]) if v*v <= (eps2*math.Abs(d[m])*math.Abs(d[m-1]) + safmin) { break } } } else { m = lend } 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×2, use Dlae2 to compute its eigenvalues. if m == l-1 { if icompz > 0 { d[l-1], d[l], work[m], work[n-1+m] = impl.Dlaev2(d[l-1], e[l-1], d[l]) impl.Dlasr(blas.Right, lapack.Variable, lapack.Forward, n, 2, work[m:], work[n-1+m:], z[l-1:], ldz) } else { d[l-1], d[l] = impl.Dlae2(d[l-1], e[l-1], d[l]) } e[l-1] = 0 l -= 2 if l < lend { break } continue } if jtot == nmaxit { break } jtot++ // Form shift. g := (d[l-1] - p) / (2 * e[l-1]) r := impl.Dlapy2(g, 1) g = d[m] - p + (e[l-1])/(g+math.Copysign(r, g)) s := 1.0 c := 1.0 p = 0.0 // Inner loop. for i := m; i < l; i++ { f := s * e[i] b := c * e[i] c, s, r = impl.Dlartg(g, f) if i != m { e[i-1] = r } g = d[i] - p r = (d[i+1]-g)*s + 2*c*b p = s * r d[i] = g + p g = c*r - b // If eigenvectors are desired, then save rotations. if icompz > 0 { work[i] = c work[n-1+i] = s } } // If eigenvectors are desired, then apply saved rotations. if icompz > 0 { mm := l - m + 1 impl.Dlasr(blas.Right, lapack.Variable, lapack.Forward, n, mm, work[m:], work[n-1+m:], z[m:], ldz) } d[l] -= p e[l-1] = g } } // Undo scaling if necessary. switch iscale { case down: // Pretend that d and e are matrices with 1 column. impl.Dlascl(lapack.General, 0, 0, ssfmax, anorm, lendsv-lsv+1, 1, d[lsv:], 1) impl.Dlascl(lapack.General, 0, 0, ssfmax, anorm, lendsv-lsv, 1, e[lsv:], 1) case up: impl.Dlascl(lapack.General, 0, 0, ssfmin, anorm, lendsv-lsv+1, 1, d[lsv:], 1) impl.Dlascl(lapack.General, 0, 0, ssfmin, anorm, lendsv-lsv, 1, e[lsv:], 1) } // Check for no convergence to an eigenvalue after a total of n*maxit iterations. if jtot >= nmaxit { break } } for i := 0; i < n-1; i++ { if e[i] != 0 { return false } } return true }