| // 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" |
| ) |
| |
| // Dgeev computes the eigenvalues and, optionally, the left and/or right |
| // eigenvectors for an n×n real nonsymmetric matrix A. |
| // |
| // The right eigenvector v_j of A corresponding to an eigenvalue λ_j |
| // is defined by |
| // A v_j = λ_j v_j, |
| // and the left eigenvector u_j corresponding to an eigenvalue λ_j is defined by |
| // u_j^H A = λ_j u_j^H, |
| // where u_j^H is the conjugate transpose of u_j. |
| // |
| // On return, A will be overwritten and the left and right eigenvectors will be |
| // stored, respectively, in the columns of the n×n matrices VL and VR in the |
| // same order as their eigenvalues. If the j-th eigenvalue is real, then |
| // u_j = VL[:,j], |
| // v_j = VR[:,j], |
| // and if it is not real, then j and j+1 form a complex conjugate pair and the |
| // eigenvectors can be recovered as |
| // u_j = VL[:,j] + i*VL[:,j+1], |
| // u_{j+1} = VL[:,j] - i*VL[:,j+1], |
| // v_j = VR[:,j] + i*VR[:,j+1], |
| // v_{j+1} = VR[:,j] - i*VR[:,j+1], |
| // where i is the imaginary unit. The computed eigenvectors are normalized to |
| // have Euclidean norm equal to 1 and largest component real. |
| // |
| // Left eigenvectors will be computed only if jobvl == lapack.ComputeLeftEV, |
| // otherwise jobvl must be lapack.None. Right eigenvectors will be computed |
| // only if jobvr == lapack.ComputeRightEV, otherwise jobvr must be lapack.None. |
| // For other values of jobvl and jobvr Dgeev will panic. |
| // |
| // wr and wi contain the real and imaginary parts, respectively, of the computed |
| // eigenvalues. Complex conjugate pairs of eigenvalues appear consecutively with |
| // the eigenvalue having the positive imaginary part first. |
| // wr and wi must have length n, and Dgeev will panic otherwise. |
| // |
| // work must have length at least lwork and lwork must be at least max(1,4*n) if |
| // the left or right eigenvectors are computed, and at least max(1,3*n) if no |
| // eigenvectors are computed. For good performance, lwork must generally be |
| // larger. On return, optimal value of lwork will be stored in work[0]. |
| // |
| // If lwork == -1, instead of performing Dgeev, the function only calculates the |
| // optimal vaule of lwork and stores it into work[0]. |
| // |
| // On return, first is the index of the first valid eigenvalue. If first == 0, |
| // all eigenvalues and eigenvectors have been computed. If first is positive, |
| // Dgeev failed to compute all the eigenvalues, no eigenvectors have been |
| // computed and wr[first:] and wi[first:] contain those eigenvalues which have |
| // converged. |
| func (impl Implementation) Dgeev(jobvl lapack.LeftEVJob, jobvr lapack.RightEVJob, n int, a []float64, lda int, wr, wi []float64, vl []float64, ldvl int, vr []float64, ldvr int, work []float64, lwork int) (first int) { |
| var wantvl bool |
| switch jobvl { |
| default: |
| panic("lapack: invalid LeftEVJob") |
| case lapack.ComputeLeftEV: |
| wantvl = true |
| case lapack.None: |
| } |
| var wantvr bool |
| switch jobvr { |
| default: |
| panic("lapack: invalid RightEVJob") |
| case lapack.ComputeRightEV: |
| wantvr = true |
| case lapack.None: |
| } |
| switch { |
| case n < 0: |
| panic(nLT0) |
| case len(work) < lwork: |
| panic(shortWork) |
| } |
| var minwrk int |
| if wantvl || wantvr { |
| minwrk = max(1, 4*n) |
| } else { |
| minwrk = max(1, 3*n) |
| } |
| if lwork != -1 { |
| checkMatrix(n, n, a, lda) |
| if wantvl { |
| checkMatrix(n, n, vl, ldvl) |
| } |
| if wantvr { |
| checkMatrix(n, n, vr, ldvr) |
| } |
| switch { |
| case len(wr) != n: |
| panic("lapack: bad length of wr") |
| case len(wi) != n: |
| panic("lapack: bad length of wi") |
| case lwork < minwrk: |
| panic(badWork) |
| } |
| } |
| |
| // Quick return if possible. |
| if n == 0 { |
| work[0] = 1 |
| return 0 |
| } |
| |
| maxwrk := 2*n + n*impl.Ilaenv(1, "DGEHRD", " ", n, 1, n, 0) |
| if wantvl || wantvr { |
| maxwrk = max(maxwrk, 2*n+(n-1)*impl.Ilaenv(1, "DORGHR", " ", n, 1, n, -1)) |
| impl.Dhseqr(lapack.EigenvaluesAndSchur, lapack.OriginalEV, n, 0, n-1, |
| nil, 1, nil, nil, nil, 1, work, -1) |
| maxwrk = max(maxwrk, max(n+1, n+int(work[0]))) |
| side := lapack.LeftEV |
| if wantvr { |
| side = lapack.RightEV |
| } |
| impl.Dtrevc3(side, lapack.AllEVMulQ, nil, n, nil, 1, nil, 1, nil, 1, |
| n, work, -1) |
| maxwrk = max(maxwrk, n+int(work[0])) |
| maxwrk = max(maxwrk, 4*n) |
| } else { |
| impl.Dhseqr(lapack.EigenvaluesOnly, lapack.None, n, 0, n-1, |
| nil, 1, nil, nil, nil, 1, work, -1) |
| maxwrk = max(maxwrk, max(n+1, n+int(work[0]))) |
| } |
| maxwrk = max(maxwrk, minwrk) |
| |
| if lwork == -1 { |
| work[0] = float64(maxwrk) |
| return 0 |
| } |
| |
| // Get machine constants. |
| smlnum := math.Sqrt(dlamchS) / dlamchP |
| bignum := 1 / smlnum |
| |
| // Scale A if max element outside range [smlnum,bignum]. |
| anrm := impl.Dlange(lapack.MaxAbs, n, n, a, lda, nil) |
| var scalea bool |
| var cscale float64 |
| if 0 < anrm && anrm < smlnum { |
| scalea = true |
| cscale = smlnum |
| } else if anrm > bignum { |
| scalea = true |
| cscale = bignum |
| } |
| if scalea { |
| impl.Dlascl(lapack.General, 0, 0, anrm, cscale, n, n, a, lda) |
| } |
| |
| // Balance the matrix. |
| workbal := work[:n] |
| ilo, ihi := impl.Dgebal(lapack.PermuteScale, n, a, lda, workbal) |
| |
| // Reduce to upper Hessenberg form. |
| iwrk := 2 * n |
| tau := work[n : iwrk-1] |
| impl.Dgehrd(n, ilo, ihi, a, lda, tau, work[iwrk:], lwork-iwrk) |
| |
| var side lapack.EVSide |
| if wantvl { |
| side = lapack.LeftEV |
| // Copy Householder vectors to VL. |
| impl.Dlacpy(blas.Lower, n, n, a, lda, vl, ldvl) |
| // Generate orthogonal matrix in VL. |
| impl.Dorghr(n, ilo, ihi, vl, ldvl, tau, work[iwrk:], lwork-iwrk) |
| // Perform QR iteration, accumulating Schur vectors in VL. |
| iwrk = n |
| first = impl.Dhseqr(lapack.EigenvaluesAndSchur, lapack.OriginalEV, n, ilo, ihi, |
| a, lda, wr, wi, vl, ldvl, work[iwrk:], lwork-iwrk) |
| if wantvr { |
| // Want left and right eigenvectors. |
| // Copy Schur vectors to VR. |
| side = lapack.RightLeftEV |
| impl.Dlacpy(blas.All, n, n, vl, ldvl, vr, ldvr) |
| } |
| } else if wantvr { |
| side = lapack.RightEV |
| // Copy Householder vectors to VR. |
| impl.Dlacpy(blas.Lower, n, n, a, lda, vr, ldvr) |
| // Generate orthogonal matrix in VR. |
| impl.Dorghr(n, ilo, ihi, vr, ldvr, tau, work[iwrk:], lwork-iwrk) |
| // Perform QR iteration, accumulating Schur vectors in VR. |
| iwrk = n |
| first = impl.Dhseqr(lapack.EigenvaluesAndSchur, lapack.OriginalEV, n, ilo, ihi, |
| a, lda, wr, wi, vr, ldvr, work[iwrk:], lwork-iwrk) |
| } else { |
| // Compute eigenvalues only. |
| iwrk = n |
| first = impl.Dhseqr(lapack.EigenvaluesOnly, lapack.None, n, ilo, ihi, |
| a, lda, wr, wi, nil, 1, work[iwrk:], lwork-iwrk) |
| } |
| |
| if first > 0 { |
| if scalea { |
| // Undo scaling. |
| impl.Dlascl(lapack.General, 0, 0, cscale, anrm, n-first, 1, wr[first:], 1) |
| impl.Dlascl(lapack.General, 0, 0, cscale, anrm, n-first, 1, wi[first:], 1) |
| impl.Dlascl(lapack.General, 0, 0, cscale, anrm, ilo, 1, wr, 1) |
| impl.Dlascl(lapack.General, 0, 0, cscale, anrm, ilo, 1, wi, 1) |
| } |
| work[0] = float64(maxwrk) |
| return first |
| } |
| |
| if wantvl || wantvr { |
| // Compute left and/or right eigenvectors. |
| impl.Dtrevc3(side, lapack.AllEVMulQ, nil, n, |
| a, lda, vl, ldvl, vr, ldvr, n, work[iwrk:], lwork-iwrk) |
| } |
| bi := blas64.Implementation() |
| if wantvl { |
| // Undo balancing of left eigenvectors. |
| impl.Dgebak(lapack.PermuteScale, lapack.LeftEV, n, ilo, ihi, workbal, n, vl, ldvl) |
| // Normalize left eigenvectors and make largest component real. |
| for i, wii := range wi { |
| if wii < 0 { |
| continue |
| } |
| if wii == 0 { |
| scl := 1 / bi.Dnrm2(n, vl[i:], ldvl) |
| bi.Dscal(n, scl, vl[i:], ldvl) |
| continue |
| } |
| scl := 1 / impl.Dlapy2(bi.Dnrm2(n, vl[i:], ldvl), bi.Dnrm2(n, vl[i+1:], ldvl)) |
| bi.Dscal(n, scl, vl[i:], ldvl) |
| bi.Dscal(n, scl, vl[i+1:], ldvl) |
| for k := 0; k < n; k++ { |
| vi := vl[k*ldvl+i] |
| vi1 := vl[k*ldvl+i+1] |
| work[iwrk+k] = vi*vi + vi1*vi1 |
| } |
| k := bi.Idamax(n, work[iwrk:iwrk+n], 1) |
| cs, sn, _ := impl.Dlartg(vl[k*ldvl+i], vl[k*ldvl+i+1]) |
| bi.Drot(n, vl[i:], ldvl, vl[i+1:], ldvl, cs, sn) |
| vl[k*ldvl+i+1] = 0 |
| } |
| } |
| if wantvr { |
| // Undo balancing of right eigenvectors. |
| impl.Dgebak(lapack.PermuteScale, lapack.RightEV, n, ilo, ihi, workbal, n, vr, ldvr) |
| // Normalize right eigenvectors and make largest component real. |
| for i, wii := range wi { |
| if wii < 0 { |
| continue |
| } |
| if wii == 0 { |
| scl := 1 / bi.Dnrm2(n, vr[i:], ldvr) |
| bi.Dscal(n, scl, vr[i:], ldvr) |
| continue |
| } |
| scl := 1 / impl.Dlapy2(bi.Dnrm2(n, vr[i:], ldvr), bi.Dnrm2(n, vr[i+1:], ldvr)) |
| bi.Dscal(n, scl, vr[i:], ldvr) |
| bi.Dscal(n, scl, vr[i+1:], ldvr) |
| for k := 0; k < n; k++ { |
| vi := vr[k*ldvr+i] |
| vi1 := vr[k*ldvr+i+1] |
| work[iwrk+k] = vi*vi + vi1*vi1 |
| } |
| k := bi.Idamax(n, work[iwrk:iwrk+n], 1) |
| cs, sn, _ := impl.Dlartg(vr[k*ldvr+i], vr[k*ldvr+i+1]) |
| bi.Drot(n, vr[i:], ldvr, vr[i+1:], ldvr, cs, sn) |
| vr[k*ldvr+i+1] = 0 |
| } |
| } |
| |
| if scalea { |
| // Undo scaling. |
| impl.Dlascl(lapack.General, 0, 0, cscale, anrm, n-first, 1, wr[first:], 1) |
| impl.Dlascl(lapack.General, 0, 0, cscale, anrm, n-first, 1, wi[first:], 1) |
| } |
| |
| work[0] = float64(maxwrk) |
| return first |
| } |