+++ /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/lapack"
-)
-
-// Dhseqr computes the eigenvalues of an n×n Hessenberg matrix H and,
-// optionally, the matrices T and Z from the Schur decomposition
-// H = Z T Z^T,
-// where T is an n×n upper quasi-triangular matrix (the Schur form), and Z is
-// the n×n orthogonal matrix of Schur vectors.
-//
-// Optionally Z may be postmultiplied into an input orthogonal matrix Q so that
-// this routine can give the Schur factorization of a matrix A which has been
-// reduced to the Hessenberg form H by the orthogonal matrix Q:
-// A = Q H Q^T = (QZ) T (QZ)^T.
-//
-// If job == lapack.EigenvaluesOnly, only the eigenvalues will be computed.
-// If job == lapack.EigenvaluesAndSchur, the eigenvalues and the Schur form T will
-// be computed.
-// For other values of job Dhseqr will panic.
-//
-// If compz == lapack.None, no Schur vectors will be computed and Z will not be
-// referenced.
-// If compz == lapack.HessEV, on return Z will contain the matrix of Schur
-// vectors of H.
-// If compz == lapack.OriginalEV, on entry z is assumed to contain the orthogonal
-// matrix Q that is the identity except for the submatrix
-// Q[ilo:ihi+1,ilo:ihi+1]. On return z will be updated to the product Q*Z.
-//
-// ilo and ihi determine the block of H on which Dhseqr operates. It is assumed
-// that H is already upper triangular in rows and columns [0:ilo] and [ihi+1:n],
-// although it will be only checked that the block is isolated, that is,
-// ilo == 0 or H[ilo,ilo-1] == 0,
-// ihi == n-1 or H[ihi+1,ihi] == 0,
-// and Dhseqr will panic otherwise. ilo and ihi are typically set by a previous
-// call to Dgebal, otherwise they should be set to 0 and n-1, respectively. It
-// must hold that
-// 0 <= ilo <= ihi < n, if n > 0,
-// ilo == 0 and ihi == -1, if n == 0.
-//
-// wr and wi must have length n.
-//
-// work must have length at least lwork and lwork must be at least max(1,n)
-// otherwise Dhseqr will panic. The minimum lwork delivers very good and
-// sometimes optimal performance, although lwork as large as 11*n may be
-// required. On return, work[0] will contain the optimal value of lwork.
-//
-// If lwork is -1, instead of performing Dhseqr, the function only estimates the
-// optimal workspace size and stores it into work[0]. Neither h nor z are
-// accessed.
-//
-// unconverged indicates whether Dhseqr computed all the eigenvalues.
-//
-// If unconverged == 0, all the eigenvalues have been computed and their real
-// and imaginary parts will be stored on return in wr and wi, respectively. If
-// two eigenvalues are computed as a complex conjugate pair, they are stored in
-// consecutive elements of wr and wi, say the i-th and (i+1)th, with wi[i] > 0
-// and wi[i+1] < 0.
-//
-// If unconverged == 0 and job == lapack.EigenvaluesAndSchur, on return H will
-// contain the upper quasi-triangular matrix T from the Schur decomposition (the
-// Schur form). 2×2 diagonal blocks (corresponding to complex conjugate pairs of
-// eigenvalues) will be returned in standard form, with
-// H[i,i] == H[i+1,i+1],
-// and
-// H[i+1,i]*H[i,i+1] < 0.
-// The eigenvalues will be stored in wr and wi in the same order as on the
-// diagonal of the Schur form returned in H, with
-// wr[i] = H[i,i],
-// and, if H[i:i+2,i:i+2] is a 2×2 diagonal block,
-// wi[i] = sqrt(-H[i+1,i]*H[i,i+1]),
-// wi[i+1] = -wi[i].
-//
-// If unconverged == 0 and job == lapack.EigenvaluesOnly, the contents of h
-// on return is unspecified.
-//
-// If unconverged > 0, some eigenvalues have not converged, and the blocks
-// [0:ilo] and [unconverged:n] of wr and wi will contain those eigenvalues which
-// have been successfully computed. Failures are rare.
-//
-// If unconverged > 0 and job == lapack.EigenvaluesOnly, on return the
-// remaining unconverged eigenvalues are the eigenvalues of the upper Hessenberg
-// matrix H[ilo:unconverged,ilo:unconverged].
-//
-// If unconverged > 0 and job == lapack.EigenvaluesAndSchur, then on
-// return
-// (initial H) U = U (final H), (*)
-// where U is an orthogonal matrix. The final H is upper Hessenberg and
-// H[unconverged:ihi+1,unconverged:ihi+1] is upper quasi-triangular.
-//
-// If unconverged > 0 and compz == lapack.OriginalEV, then on return
-// (final Z) = (initial Z) U,
-// where U is the orthogonal matrix in (*) regardless of the value of job.
-//
-// If unconverged > 0 and compz == lapack.HessEV, then on return
-// (final Z) = U,
-// where U is the orthogonal matrix in (*) regardless of the value of job.
-//
-// References:
-// [1] R. Byers. LAPACK 3.1 xHSEQR: Tuning and Implementation Notes on the
-// Small Bulge Multi-Shift QR Algorithm with Aggressive Early Deflation.
-// LAPACK Working Note 187 (2007)
-// URL: http://www.netlib.org/lapack/lawnspdf/lawn187.pdf
-// [2] K. Braman, R. Byers, R. Mathias. The Multishift QR Algorithm. Part I:
-// Maintaining Well-Focused Shifts and Level 3 Performance. SIAM J. Matrix
-// Anal. Appl. 23(4) (2002), pp. 929—947
-// URL: http://dx.doi.org/10.1137/S0895479801384573
-// [3] K. Braman, R. Byers, R. Mathias. The Multishift QR Algorithm. Part II:
-// Aggressive Early Deflation. SIAM J. Matrix Anal. Appl. 23(4) (2002), pp. 948—973
-// URL: http://dx.doi.org/10.1137/S0895479801384585
-//
-// Dhseqr is an internal routine. It is exported for testing purposes.
-func (impl Implementation) Dhseqr(job lapack.EVJob, compz lapack.EVComp, n, ilo, ihi int, h []float64, ldh int, wr, wi []float64, z []float64, ldz int, work []float64, lwork int) (unconverged int) {
- var wantt bool
- switch job {
- default:
- panic(badEVJob)
- case lapack.EigenvaluesOnly:
- case lapack.EigenvaluesAndSchur:
- wantt = true
- }
- var wantz bool
- switch compz {
- default:
- panic(badEVComp)
- case lapack.None:
- case lapack.HessEV, lapack.OriginalEV:
- wantz = true
- }
- switch {
- case n < 0:
- panic(nLT0)
- case ilo < 0 || max(0, n-1) < ilo:
- panic(badIlo)
- case ihi < min(ilo, n-1) || n <= ihi:
- panic(badIhi)
- case len(work) < lwork:
- panic(shortWork)
- case lwork < max(1, n) && lwork != -1:
- panic(badWork)
- }
- if lwork != -1 {
- checkMatrix(n, n, h, ldh)
- switch {
- case wantz:
- checkMatrix(n, n, z, ldz)
- case len(wr) < n:
- panic("lapack: wr has insufficient length")
- case len(wi) < n:
- panic("lapack: wi has insufficient length")
- }
- }
-
- const (
- // Matrices of order ntiny or smaller must be processed by
- // Dlahqr because of insufficient subdiagonal scratch space.
- // This is a hard limit.
- ntiny = 11
-
- // nl is the size of a local workspace to help small matrices
- // through a rare Dlahqr failure. nl > ntiny is required and
- // nl <= nmin = Ilaenv(ispec=12,...) is recommended (the default
- // value of nmin is 75). Using nl = 49 allows up to six
- // simultaneous shifts and a 16×16 deflation window.
- nl = 49
- )
-
- // Quick return if possible.
- if n == 0 {
- work[0] = 1
- return 0
- }
-
- // Quick return in case of a workspace query.
- if lwork == -1 {
- impl.Dlaqr04(wantt, wantz, n, ilo, ihi, nil, 0, nil, nil, ilo, ihi, nil, 0, work, -1, 1)
- work[0] = math.Max(float64(n), work[0])
- return 0
- }
-
- // Copy eigenvalues isolated by Dgebal.
- for i := 0; i < ilo; i++ {
- wr[i] = h[i*ldh+i]
- wi[i] = 0
- }
- for i := ihi + 1; i < n; i++ {
- wr[i] = h[i*ldh+i]
- wi[i] = 0
- }
-
- // Initialize Z to identity matrix if requested.
- if compz == lapack.HessEV {
- impl.Dlaset(blas.All, n, n, 0, 1, z, ldz)
- }
-
- // Quick return if possible.
- if ilo == ihi {
- wr[ilo] = h[ilo*ldh+ilo]
- wi[ilo] = 0
- return 0
- }
-
- // Dlahqr/Dlaqr04 crossover point.
- nmin := impl.Ilaenv(12, "DHSEQR", string(job)+string(compz), n, ilo, ihi, lwork)
- nmin = max(ntiny, nmin)
-
- if n > nmin {
- // Dlaqr0 for big matrices.
- unconverged = impl.Dlaqr04(wantt, wantz, n, ilo, ihi, h, ldh, wr[:ihi+1], wi[:ihi+1],
- ilo, ihi, z, ldz, work, lwork, 1)
- } else {
- // Dlahqr for small matrices.
- unconverged = impl.Dlahqr(wantt, wantz, n, ilo, ihi, h, ldh, wr[:ihi+1], wi[:ihi+1],
- ilo, ihi, z, ldz)
- if unconverged > 0 {
- // A rare Dlahqr failure! Dlaqr04 sometimes succeeds
- // when Dlahqr fails.
- kbot := unconverged
- if n >= nl {
- // Larger matrices have enough subdiagonal
- // scratch space to call Dlaqr04 directly.
- unconverged = impl.Dlaqr04(wantt, wantz, n, ilo, kbot, h, ldh,
- wr[:ihi+1], wi[:ihi+1], ilo, ihi, z, ldz, work, lwork, 1)
- } else {
- // Tiny matrices don't have enough subdiagonal
- // scratch space to benefit from Dlaqr04. Hence,
- // tiny matrices must be copied into a larger
- // array before calling Dlaqr04.
- var hl [nl * nl]float64
- impl.Dlacpy(blas.All, n, n, h, ldh, hl[:], nl)
- impl.Dlaset(blas.All, nl, nl-n, 0, 0, hl[n:], nl)
- var workl [nl]float64
- unconverged = impl.Dlaqr04(wantt, wantz, nl, ilo, kbot, hl[:], nl,
- wr[:ihi+1], wi[:ihi+1], ilo, ihi, z, ldz, workl[:], nl, 1)
- work[0] = workl[0]
- if wantt || unconverged > 0 {
- impl.Dlacpy(blas.All, n, n, hl[:], nl, h, ldh)
- }
- }
- }
- }
- // Zero out under the first subdiagonal, if necessary.
- if (wantt || unconverged > 0) && n > 2 {
- impl.Dlaset(blas.Lower, n-2, n-2, 0, 0, h[2*ldh:], ldh)
- }
-
- work[0] = math.Max(float64(n), work[0])
- return unconverged
-}
OSDN Git Service
