--- /dev/null
+// 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
+}
OSDN Git Service
