+++ /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
