+++ /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"
-)
-
-// Dlaqr23 performs the orthogonal similarity transformation of an n×n upper
-// Hessenberg matrix to detect and deflate fully converged eigenvalues from a
-// trailing principal submatrix using aggressive early deflation [1].
-//
-// On return, H will be overwritten by a new Hessenberg matrix that is a
-// perturbation of an orthogonal similarity transformation of H. It is hoped
-// that on output H will have many zero subdiagonal entries.
-//
-// If wantt is true, the matrix H will be fully updated so that the
-// quasi-triangular Schur factor can be computed. If wantt is false, then only
-// enough of H will be updated to preserve the eigenvalues.
-//
-// If wantz is true, the orthogonal similarity transformation will be
-// accumulated into Z[iloz:ihiz+1,ktop:kbot+1], otherwise Z is not referenced.
-//
-// ktop and kbot determine a block [ktop:kbot+1,ktop:kbot+1] along the diagonal
-// of H. It must hold that
-// 0 <= ilo <= ihi < n, if n > 0,
-// ilo == 0 and ihi == -1, if n == 0,
-// and the block must be isolated, that is, it must hold that
-// ktop == 0 or H[ktop,ktop-1] == 0,
-// kbot == n-1 or H[kbot+1,kbot] == 0,
-// otherwise Dlaqr23 will panic.
-//
-// nw is the deflation window size. It must hold that
-// 0 <= nw <= kbot-ktop+1,
-// otherwise Dlaqr23 will panic.
-//
-// iloz and ihiz specify the rows of the n×n matrix Z to which transformations
-// will be applied if wantz is true. It must hold that
-// 0 <= iloz <= ktop, and kbot <= ihiz < n,
-// otherwise Dlaqr23 will panic.
-//
-// sr and si must have length kbot+1, otherwise Dlaqr23 will panic.
-//
-// v and ldv represent an nw×nw work matrix.
-// t and ldt represent an nw×nh work matrix, and nh must be at least nw.
-// wv and ldwv represent an nv×nw work matrix.
-//
-// work must have length at least lwork and lwork must be at least max(1,2*nw),
-// otherwise Dlaqr23 will panic. Larger values of lwork may result in greater
-// efficiency. On return, work[0] will contain the optimal value of lwork.
-//
-// If lwork is -1, instead of performing Dlaqr23, the function only estimates the
-// optimal workspace size and stores it into work[0]. Neither h nor z are
-// accessed.
-//
-// recur is the non-negative recursion depth. For recur > 0, Dlaqr23 behaves
-// as DLAQR3, for recur == 0 it behaves as DLAQR2.
-//
-// On return, ns and nd will contain respectively the number of unconverged
-// (i.e., approximate) eigenvalues and converged eigenvalues that are stored in
-// sr and si.
-//
-// On return, the real and imaginary parts of approximate eigenvalues that may
-// be used for shifts will be stored respectively in sr[kbot-nd-ns+1:kbot-nd+1]
-// and si[kbot-nd-ns+1:kbot-nd+1].
-//
-// On return, the real and imaginary parts of converged eigenvalues will be
-// stored respectively in sr[kbot-nd+1:kbot+1] and si[kbot-nd+1:kbot+1].
-//
-// References:
-// [1] 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
-//
-func (impl Implementation) Dlaqr23(wantt, wantz bool, n, ktop, kbot, nw int, h []float64, ldh int, iloz, ihiz int, z []float64, ldz int, sr, si []float64, v []float64, ldv int, nh int, t []float64, ldt int, nv int, wv []float64, ldwv int, work []float64, lwork int, recur int) (ns, nd int) {
- switch {
- case ktop < 0 || max(0, n-1) < ktop:
- panic("lapack: invalid value of ktop")
- case kbot < min(ktop, n-1) || n <= kbot:
- panic("lapack: invalid value of kbot")
- case (nw < 0 || kbot-ktop+1 < nw) && lwork != -1:
- panic("lapack: invalid value of nw")
- case nh < nw:
- panic("lapack: invalid value of nh")
- case lwork < max(1, 2*nw) && lwork != -1:
- panic(badWork)
- case len(work) < lwork:
- panic(shortWork)
- case recur < 0:
- panic("lapack: recur is negative")
- }
- if wantz {
- switch {
- case iloz < 0 || ktop < iloz:
- panic("lapack: invalid value of iloz")
- case ihiz < kbot || n <= ihiz:
- panic("lapack: invalid value of ihiz")
- }
- }
- if lwork != -1 {
- // Check input slices only if not doing workspace query.
- checkMatrix(n, n, h, ldh)
- checkMatrix(nw, nw, v, ldv)
- checkMatrix(nw, nh, t, ldt)
- checkMatrix(nv, nw, wv, ldwv)
- if wantz {
- checkMatrix(n, n, z, ldz)
- }
- switch {
- case ktop > 0 && h[ktop*ldh+ktop-1] != 0:
- panic("lapack: block not isolated")
- case kbot+1 < n && h[(kbot+1)*ldh+kbot] != 0:
- panic("lapack: block not isolated")
- case len(sr) != kbot+1:
- panic("lapack: bad length of sr")
- case len(si) != kbot+1:
- panic("lapack: bad length of si")
- }
- }
-
- // Quick return for zero window size.
- if nw == 0 {
- work[0] = 1
- return 0, 0
- }
-
- // LAPACK code does not enforce the documented behavior
- // nw <= kbot-ktop+1
- // but we do (we panic above).
- jw := nw
- lwkopt := max(1, 2*nw)
- if jw > 2 {
- // Workspace query call to Dgehrd.
- impl.Dgehrd(jw, 0, jw-2, nil, 0, nil, work, -1)
- lwk1 := int(work[0])
- // Workspace query call to Dormhr.
- impl.Dormhr(blas.Right, blas.NoTrans, jw, jw, 0, jw-2, nil, 0, nil, nil, 0, work, -1)
- lwk2 := int(work[0])
- if recur > 0 {
- // Workspace query call to Dlaqr04.
- impl.Dlaqr04(true, true, jw, 0, jw-1, nil, 0, nil, nil, 0, jw-1, nil, 0, work, -1, recur-1)
- lwk3 := int(work[0])
- // Optimal workspace.
- lwkopt = max(jw+max(lwk1, lwk2), lwk3)
- } else {
- // Optimal workspace.
- lwkopt = jw + max(lwk1, lwk2)
- }
- }
- // Quick return in case of workspace query.
- if lwork == -1 {
- work[0] = float64(lwkopt)
- return 0, 0
- }
-
- // Machine constants.
- ulp := dlamchP
- smlnum := float64(n) / ulp * dlamchS
-
- // Setup deflation window.
- var s float64
- kwtop := kbot - jw + 1
- if kwtop != ktop {
- s = h[kwtop*ldh+kwtop-1]
- }
- if kwtop == kbot {
- // 1×1 deflation window.
- sr[kwtop] = h[kwtop*ldh+kwtop]
- si[kwtop] = 0
- ns = 1
- nd = 0
- if math.Abs(s) <= math.Max(smlnum, ulp*math.Abs(h[kwtop*ldh+kwtop])) {
- ns = 0
- nd = 1
- if kwtop > ktop {
- h[kwtop*ldh+kwtop-1] = 0
- }
- }
- work[0] = 1
- return ns, nd
- }
-
- // Convert to spike-triangular form. In case of a rare QR failure, this
- // routine continues to do aggressive early deflation using that part of
- // the deflation window that converged using infqr here and there to
- // keep track.
- impl.Dlacpy(blas.Upper, jw, jw, h[kwtop*ldh+kwtop:], ldh, t, ldt)
- bi := blas64.Implementation()
- bi.Dcopy(jw-1, h[(kwtop+1)*ldh+kwtop:], ldh+1, t[ldt:], ldt+1)
- impl.Dlaset(blas.All, jw, jw, 0, 1, v, ldv)
- nmin := impl.Ilaenv(12, "DLAQR3", "SV", jw, 0, jw-1, lwork)
- var infqr int
- if recur > 0 && jw > nmin {
- infqr = impl.Dlaqr04(true, true, jw, 0, jw-1, t, ldt, sr[kwtop:], si[kwtop:], 0, jw-1, v, ldv, work, lwork, recur-1)
- } else {
- infqr = impl.Dlahqr(true, true, jw, 0, jw-1, t, ldt, sr[kwtop:], si[kwtop:], 0, jw-1, v, ldv)
- }
- // Note that ilo == 0 which conveniently coincides with the success
- // value of infqr, that is, infqr as an index always points to the first
- // converged eigenvalue.
-
- // Dtrexc needs a clean margin near the diagonal.
- for j := 0; j < jw-3; j++ {
- t[(j+2)*ldt+j] = 0
- t[(j+3)*ldt+j] = 0
- }
- if jw >= 3 {
- t[(jw-1)*ldt+jw-3] = 0
- }
-
- ns = jw
- ilst := infqr
- // Deflation detection loop.
- for ilst < ns {
- bulge := false
- if ns >= 2 {
- bulge = t[(ns-1)*ldt+ns-2] != 0
- }
- if !bulge {
- // Real eigenvalue.
- abst := math.Abs(t[(ns-1)*ldt+ns-1])
- if abst == 0 {
- abst = math.Abs(s)
- }
- if math.Abs(s*v[ns-1]) <= math.Max(smlnum, ulp*abst) {
- // Deflatable.
- ns--
- } else {
- // Undeflatable, move it up out of the way.
- // Dtrexc can not fail in this case.
- _, ilst, _ = impl.Dtrexc(lapack.UpdateSchur, jw, t, ldt, v, ldv, ns-1, ilst, work)
- ilst++
- }
- continue
- }
- // Complex conjugate pair.
- abst := math.Abs(t[(ns-1)*ldt+ns-1]) + math.Sqrt(math.Abs(t[(ns-1)*ldt+ns-2]))*math.Sqrt(math.Abs(t[(ns-2)*ldt+ns-1]))
- if abst == 0 {
- abst = math.Abs(s)
- }
- if math.Max(math.Abs(s*v[ns-1]), math.Abs(s*v[ns-2])) <= math.Max(smlnum, ulp*abst) {
- // Deflatable.
- ns -= 2
- } else {
- // Undeflatable, move them up out of the way.
- // Dtrexc does the right thing with ilst in case of a
- // rare exchange failure.
- _, ilst, _ = impl.Dtrexc(lapack.UpdateSchur, jw, t, ldt, v, ldv, ns-1, ilst, work)
- ilst += 2
- }
- }
-
- // Return to Hessenberg form.
- if ns == 0 {
- s = 0
- }
- if ns < jw {
- // Sorting diagonal blocks of T improves accuracy for graded
- // matrices. Bubble sort deals well with exchange failures.
- sorted := false
- i := ns
- for !sorted {
- sorted = true
- kend := i - 1
- i = infqr
- var k int
- if i == ns-1 || t[(i+1)*ldt+i] == 0 {
- k = i + 1
- } else {
- k = i + 2
- }
- for k <= kend {
- var evi float64
- if k == i+1 {
- evi = math.Abs(t[i*ldt+i])
- } else {
- evi = math.Abs(t[i*ldt+i]) + math.Sqrt(math.Abs(t[(i+1)*ldt+i]))*math.Sqrt(math.Abs(t[i*ldt+i+1]))
- }
-
- var evk float64
- if k == kend || t[(k+1)*ldt+k] == 0 {
- evk = math.Abs(t[k*ldt+k])
- } else {
- evk = math.Abs(t[k*ldt+k]) + math.Sqrt(math.Abs(t[(k+1)*ldt+k]))*math.Sqrt(math.Abs(t[k*ldt+k+1]))
- }
-
- if evi >= evk {
- i = k
- } else {
- sorted = false
- _, ilst, ok := impl.Dtrexc(lapack.UpdateSchur, jw, t, ldt, v, ldv, i, k, work)
- if ok {
- i = ilst
- } else {
- i = k
- }
- }
- if i == kend || t[(i+1)*ldt+i] == 0 {
- k = i + 1
- } else {
- k = i + 2
- }
- }
- }
- }
-
- // Restore shift/eigenvalue array from T.
- for i := jw - 1; i >= infqr; {
- if i == infqr || t[i*ldt+i-1] == 0 {
- sr[kwtop+i] = t[i*ldt+i]
- si[kwtop+i] = 0
- i--
- continue
- }
- aa := t[(i-1)*ldt+i-1]
- bb := t[(i-1)*ldt+i]
- cc := t[i*ldt+i-1]
- dd := t[i*ldt+i]
- _, _, _, _, sr[kwtop+i-1], si[kwtop+i-1], sr[kwtop+i], si[kwtop+i], _, _ = impl.Dlanv2(aa, bb, cc, dd)
- i -= 2
- }
-
- if ns < jw || s == 0 {
- if ns > 1 && s != 0 {
- // Reflect spike back into lower triangle.
- bi.Dcopy(ns, v[:ns], 1, work[:ns], 1)
- _, tau := impl.Dlarfg(ns, work[0], work[1:ns], 1)
- work[0] = 1
- impl.Dlaset(blas.Lower, jw-2, jw-2, 0, 0, t[2*ldt:], ldt)
- impl.Dlarf(blas.Left, ns, jw, work[:ns], 1, tau, t, ldt, work[jw:])
- impl.Dlarf(blas.Right, ns, ns, work[:ns], 1, tau, t, ldt, work[jw:])
- impl.Dlarf(blas.Right, jw, ns, work[:ns], 1, tau, v, ldv, work[jw:])
- impl.Dgehrd(jw, 0, ns-1, t, ldt, work[:jw-1], work[jw:], lwork-jw)
- }
-
- // Copy updated reduced window into place.
- if kwtop > 0 {
- h[kwtop*ldh+kwtop-1] = s * v[0]
- }
- impl.Dlacpy(blas.Upper, jw, jw, t, ldt, h[kwtop*ldh+kwtop:], ldh)
- bi.Dcopy(jw-1, t[ldt:], ldt+1, h[(kwtop+1)*ldh+kwtop:], ldh+1)
-
- // Accumulate orthogonal matrix in order to update H and Z, if
- // requested.
- if ns > 1 && s != 0 {
- // work[:ns-1] contains the elementary reflectors stored
- // by a call to Dgehrd above.
- impl.Dormhr(blas.Right, blas.NoTrans, jw, ns, 0, ns-1,
- t, ldt, work[:ns-1], v, ldv, work[jw:], lwork-jw)
- }
-
- // Update vertical slab in H.
- var ltop int
- if !wantt {
- ltop = ktop
- }
- for krow := ltop; krow < kwtop; krow += nv {
- kln := min(nv, kwtop-krow)
- bi.Dgemm(blas.NoTrans, blas.NoTrans, kln, jw, jw,
- 1, h[krow*ldh+kwtop:], ldh, v, ldv,
- 0, wv, ldwv)
- impl.Dlacpy(blas.All, kln, jw, wv, ldwv, h[krow*ldh+kwtop:], ldh)
- }
-
- // Update horizontal slab in H.
- if wantt {
- for kcol := kbot + 1; kcol < n; kcol += nh {
- kln := min(nh, n-kcol)
- bi.Dgemm(blas.Trans, blas.NoTrans, jw, kln, jw,
- 1, v, ldv, h[kwtop*ldh+kcol:], ldh,
- 0, t, ldt)
- impl.Dlacpy(blas.All, jw, kln, t, ldt, h[kwtop*ldh+kcol:], ldh)
- }
- }
-
- // Update vertical slab in Z.
- if wantz {
- for krow := iloz; krow <= ihiz; krow += nv {
- kln := min(nv, ihiz-krow+1)
- bi.Dgemm(blas.NoTrans, blas.NoTrans, kln, jw, jw,
- 1, z[krow*ldz+kwtop:], ldz, v, ldv,
- 0, wv, ldwv)
- impl.Dlacpy(blas.All, kln, jw, wv, ldwv, z[krow*ldz+kwtop:], ldz)
- }
- }
- }
-
- // The number of deflations.
- nd = jw - ns
- // Shifts are converged eigenvalues that could not be deflated.
- // Subtracting infqr from the spike length takes care of the case of a
- // rare QR failure while calculating eigenvalues of the deflation
- // window.
- ns -= infqr
- work[0] = float64(lwkopt)
- return ns, nd
-}
OSDN Git Service
