1 // Copyright ©2016 The Gonum Authors. All rights reserved.
2 // Use of this source code is governed by a BSD-style
3 // license that can be found in the LICENSE file.
10 "gonum.org/v1/gonum/blas"
11 "gonum.org/v1/gonum/blas/blas64"
12 "gonum.org/v1/gonum/lapack"
15 // Dsteqr computes the eigenvalues and optionally the eigenvectors of a symmetric
16 // tridiagonal matrix using the implicit QL or QR method. The eigenvectors of a
17 // full or band symmetric matrix can also be found if Dsytrd, Dsptrd, or Dsbtrd
18 // have been used to reduce this matrix to tridiagonal form.
20 // d, on entry, contains the diagonal elements of the tridiagonal matrix. On exit,
21 // d contains the eigenvalues in ascending order. d must have length n and
22 // Dsteqr will panic otherwise.
24 // e, on entry, contains the off-diagonal elements of the tridiagonal matrix on
25 // entry, and is overwritten during the call to Dsteqr. e must have length n-1 and
26 // Dsteqr will panic otherwise.
28 // z, on entry, contains the n×n orthogonal matrix used in the reduction to
29 // tridiagonal form if compz == lapack.OriginalEV. On exit, if
30 // compz == lapack.OriginalEV, z contains the orthonormal eigenvectors of the
31 // original symmetric matrix, and if compz == lapack.TridiagEV, z contains the
32 // orthonormal eigenvectors of the symmetric tridiagonal matrix. z is not used
33 // if compz == lapack.None.
35 // work must have length at least max(1, 2*n-2) if the eigenvectors are computed,
36 // and Dsteqr will panic otherwise.
38 // Dsteqr is an internal routine. It is exported for testing purposes.
39 func (impl Implementation) Dsteqr(compz lapack.EVComp, n int, d, e, z []float64, ldz int, work []float64) (ok bool) {
49 if compz != lapack.None && compz != lapack.TridiagEV && compz != lapack.OriginalEV {
52 if compz != lapack.None {
53 if len(work) < max(1, 2*n-2) {
56 checkMatrix(n, n, z, ldz)
60 if compz == lapack.OriginalEV {
62 } else if compz == lapack.TridiagEV {
76 bi := blas64.Implementation()
82 ssfmax := math.Sqrt(safmax) / 3
83 ssfmin := math.Sqrt(safmin) / eps2
85 // Compute the eigenvalues and eigenvectors of the tridiagonal matrix.
87 impl.Dlaset(blas.All, n, n, 0, 1, z, ldz)
94 // Determine where the matrix splits and choose QL or QR iteration for each
95 // block, according to whether top or bottom diagonal element is smaller.
101 down scaletype = iota + 1
108 // Order eigenvalues and eigenvectors.
110 impl.Dlasrt(lapack.SortIncreasing, n, d)
112 // TODO(btracey): Consider replacing this sort with a call to sort.Sort.
113 for ii := 1; ii < n; ii++ {
117 for j := ii; j < n; j++ {
126 bi.Dswap(n, z[i:], ldz, z[k:], ldz)
137 for m = l1; m < nm1; m++ {
138 test := math.Abs(e[m])
142 if test <= (math.Sqrt(math.Abs(d[m]))*math.Sqrt(math.Abs(d[m+1])))*eps {
157 // Scale submatrix in rows and columns L to Lend
158 anorm := impl.Dlanst(lapack.MaxAbs, lend-l+1, d[l:], e[l:])
164 // Pretend that d and e are matrices with 1 column.
165 impl.Dlascl(lapack.General, 0, 0, anorm, ssfmax, lend-l+1, 1, d[l:], 1)
166 impl.Dlascl(lapack.General, 0, 0, anorm, ssfmax, lend-l, 1, e[l:], 1)
169 impl.Dlascl(lapack.General, 0, 0, anorm, ssfmin, lend-l+1, 1, d[l:], 1)
170 impl.Dlascl(lapack.General, 0, 0, anorm, ssfmin, lend-l, 1, e[l:], 1)
173 // Choose between QL and QR.
174 if math.Abs(d[lend]) < math.Abs(d[l]) {
179 // QL Iteration. Look for small subdiagonal element.
182 for m = l; m < lend; m++ {
184 if v*v <= (eps2*math.Abs(d[m]))*math.Abs(d[m+1])+safmin {
204 // If remaining matrix is 2×2, use Dlae2 to compute its eigensystem.
207 d[l], d[l+1], work[l], work[n-1+l] = impl.Dlaev2(d[l], e[l], d[l+1])
208 impl.Dlasr(blas.Right, lapack.Variable, lapack.Backward,
209 n, 2, work[l:], work[n-1+l:], z[l:], ldz)
211 d[l], d[l+1] = impl.Dlae2(d[l], e[l], d[l+1])
227 g := (d[l+1] - p) / (2 * e[l])
228 r := impl.Dlapy2(g, 1)
229 g = d[m] - p + e[l]/(g+math.Copysign(r, g))
235 for i := m - 1; i >= l; i-- {
238 c, s, r = impl.Dlartg(g, f)
243 r = (d[i]-g)*s + 2*c*b
248 // If eigenvectors are desired, then save rotations.
254 // If eigenvectors are desired, then apply saved rotations.
257 impl.Dlasr(blas.Right, lapack.Variable, lapack.Backward,
258 n, mm, work[l:], work[n-1+l:], z[l:], ldz)
265 // Look for small superdiagonal element.
268 for m = l; m > lend; m-- {
269 v := math.Abs(e[m-1])
270 if v*v <= (eps2*math.Abs(d[m])*math.Abs(d[m-1]) + safmin) {
290 // If remaining matrix is 2×2, use Dlae2 to compute its eigenvalues.
293 d[l-1], d[l], work[m], work[n-1+m] = impl.Dlaev2(d[l-1], e[l-1], d[l])
294 impl.Dlasr(blas.Right, lapack.Variable, lapack.Forward,
295 n, 2, work[m:], work[n-1+m:], z[l-1:], ldz)
297 d[l-1], d[l] = impl.Dlae2(d[l-1], e[l-1], d[l])
312 g := (d[l-1] - p) / (2 * e[l-1])
313 r := impl.Dlapy2(g, 1)
314 g = d[m] - p + (e[l-1])/(g+math.Copysign(r, g))
320 for i := m; i < l; i++ {
323 c, s, r = impl.Dlartg(g, f)
328 r = (d[i+1]-g)*s + 2*c*b
333 // If eigenvectors are desired, then save rotations.
340 // If eigenvectors are desired, then apply saved rotations.
343 impl.Dlasr(blas.Right, lapack.Variable, lapack.Forward,
344 n, mm, work[m:], work[n-1+m:], z[m:], ldz)
351 // Undo scaling if necessary.
354 // Pretend that d and e are matrices with 1 column.
355 impl.Dlascl(lapack.General, 0, 0, ssfmax, anorm, lendsv-lsv+1, 1, d[lsv:], 1)
356 impl.Dlascl(lapack.General, 0, 0, ssfmax, anorm, lendsv-lsv, 1, e[lsv:], 1)
358 impl.Dlascl(lapack.General, 0, 0, ssfmin, anorm, lendsv-lsv+1, 1, d[lsv:], 1)
359 impl.Dlascl(lapack.General, 0, 0, ssfmin, anorm, lendsv-lsv, 1, e[lsv:], 1)
362 // Check for no convergence to an eigenvalue after a total of n*maxit iterations.
367 for i := 0; i < n-1; i++ {