| // 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.OriginalEV. On exit, if |
| // compz == lapack.OriginalEV, z contains the orthonormal eigenvectors of the |
| // original symmetric matrix, and if compz == lapack.TridiagEV, z contains the |
| // orthonormal eigenvectors of the symmetric tridiagonal matrix. z is not used |
| // if compz == lapack.None. |
| // |
| // 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) { |
| if n < 0 { |
| panic(nLT0) |
| } |
| if len(d) < n { |
| panic(badD) |
| } |
| if len(e) < n-1 { |
| panic(badE) |
| } |
| if compz != lapack.None && compz != lapack.TridiagEV && compz != lapack.OriginalEV { |
| panic(badEVComp) |
| } |
| if compz != lapack.None { |
| if len(work) < max(1, 2*n-2) { |
| panic(badWork) |
| } |
| checkMatrix(n, n, z, ldz) |
| } |
| |
| var icompz int |
| if compz == lapack.OriginalEV { |
| icompz = 1 |
| } else if compz == lapack.TridiagEV { |
| icompz = 2 |
| } |
| |
| if n == 0 { |
| return true |
| } |
| 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 |
| } |
