vendor/gonum.org/v1/gonum/lapack/gonum/dgeev.go

   1 // Copyright ©2016 The Gonum Authors. All rights reserved.
   2 // Use of this source code is governed by a BSD-style
   3 // license that can be found in the LICENSE file.
   4
   5 package gonum
   6
   7 import (
   8         "math"
   9
  10         "gonum.org/v1/gonum/blas"
  11         "gonum.org/v1/gonum/blas/blas64"
  12         "gonum.org/v1/gonum/lapack"
  13 )
  14
  15 // Dgeev computes the eigenvalues and, optionally, the left and/or right
  16 // eigenvectors for an n×n real nonsymmetric matrix A.
  17 //
  18 // The right eigenvector v_j of A corresponding to an eigenvalue λ_j
  19 // is defined by
  20 //  A v_j = λ_j v_j,
  21 // and the left eigenvector u_j corresponding to an eigenvalue λ_j is defined by
  22 //  u_j^H A = λ_j u_j^H,
  23 // where u_j^H is the conjugate transpose of u_j.
  24 //
  25 // On return, A will be overwritten and the left and right eigenvectors will be
  26 // stored, respectively, in the columns of the n×n matrices VL and VR in the
  27 // same order as their eigenvalues. If the j-th eigenvalue is real, then
  28 //  u_j = VL[:,j],
  29 //  v_j = VR[:,j],
  30 // and if it is not real, then j and j+1 form a complex conjugate pair and the
  31 // eigenvectors can be recovered as
  32 //  u_j     = VL[:,j] + i*VL[:,j+1],
  33 //  u_{j+1} = VL[:,j] - i*VL[:,j+1],
  34 //  v_j     = VR[:,j] + i*VR[:,j+1],
  35 //  v_{j+1} = VR[:,j] - i*VR[:,j+1],
  36 // where i is the imaginary unit. The computed eigenvectors are normalized to
  37 // have Euclidean norm equal to 1 and largest component real.
  38 //
  39 // Left eigenvectors will be computed only if jobvl == lapack.ComputeLeftEV,
  40 // otherwise jobvl must be lapack.None. Right eigenvectors will be computed
  41 // only if jobvr == lapack.ComputeRightEV, otherwise jobvr must be lapack.None.
  42 // For other values of jobvl and jobvr Dgeev will panic.
  43 //
  44 // wr and wi contain the real and imaginary parts, respectively, of the computed
  45 // eigenvalues. Complex conjugate pairs of eigenvalues appear consecutively with
  46 // the eigenvalue having the positive imaginary part first.
  47 // wr and wi must have length n, and Dgeev will panic otherwise.
  48 //
  49 // work must have length at least lwork and lwork must be at least max(1,4*n) if
  50 // the left or right eigenvectors are computed, and at least max(1,3*n) if no
  51 // eigenvectors are computed. For good performance, lwork must generally be
  52 // larger.  On return, optimal value of lwork will be stored in work[0].
  53 //
  54 // If lwork == -1, instead of performing Dgeev, the function only calculates the
  55 // optimal vaule of lwork and stores it into work[0].
  56 //
  57 // On return, first is the index of the first valid eigenvalue. If first == 0,
  58 // all eigenvalues and eigenvectors have been computed. If first is positive,
  59 // Dgeev failed to compute all the eigenvalues, no eigenvectors have been
  60 // computed and wr[first:] and wi[first:] contain those eigenvalues which have
  61 // converged.
  62 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) {
  63         var wantvl bool
  64         switch jobvl {
  65         default:
  66                 panic("lapack: invalid LeftEVJob")
  67         case lapack.ComputeLeftEV:
  68                 wantvl = true
  69         case lapack.None:
  70         }
  71         var wantvr bool
  72         switch jobvr {
  73         default:
  74                 panic("lapack: invalid RightEVJob")
  75         case lapack.ComputeRightEV:
  76                 wantvr = true
  77         case lapack.None:
  78         }
  79         switch {
  80         case n < 0:
  81                 panic(nLT0)
  82         case len(work) < lwork:
  83                 panic(shortWork)
  84         }
  85         var minwrk int
  86         if wantvl || wantvr {
  87                 minwrk = max(1, 4*n)
  88         } else {
  89                 minwrk = max(1, 3*n)
  90         }
  91         if lwork != -1 {
  92                 checkMatrix(n, n, a, lda)
  93                 if wantvl {
  94                         checkMatrix(n, n, vl, ldvl)
  95                 }
  96                 if wantvr {
  97                         checkMatrix(n, n, vr, ldvr)
  98                 }
  99                 switch {
 100                 case len(wr) != n:
 101                         panic("lapack: bad length of wr")
 102                 case len(wi) != n:
 103                         panic("lapack: bad length of wi")
 104                 case lwork < minwrk:
 105                         panic(badWork)
 106                 }
 107         }
 108
 109         // Quick return if possible.
 110         if n == 0 {
 111                 work[0] = 1
 112                 return 0
 113         }
 114
 115         maxwrk := 2*n + n*impl.Ilaenv(1, "DGEHRD", " ", n, 1, n, 0)
 116         if wantvl || wantvr {
 117                 maxwrk = max(maxwrk, 2*n+(n-1)*impl.Ilaenv(1, "DORGHR", " ", n, 1, n, -1))
 118                 impl.Dhseqr(lapack.EigenvaluesAndSchur, lapack.OriginalEV, n, 0, n-1,
 119                         nil, 1, nil, nil, nil, 1, work, -1)
 120                 maxwrk = max(maxwrk, max(n+1, n+int(work[0])))
 121                 side := lapack.LeftEV
 122                 if wantvr {
 123                         side = lapack.RightEV
 124                 }
 125                 impl.Dtrevc3(side, lapack.AllEVMulQ, nil, n, nil, 1, nil, 1, nil, 1,
 126                         n, work, -1)
 127                 maxwrk = max(maxwrk, n+int(work[0]))
 128                 maxwrk = max(maxwrk, 4*n)
 129         } else {
 130                 impl.Dhseqr(lapack.EigenvaluesOnly, lapack.None, n, 0, n-1,
 131                         nil, 1, nil, nil, nil, 1, work, -1)
 132                 maxwrk = max(maxwrk, max(n+1, n+int(work[0])))
 133         }
 134         maxwrk = max(maxwrk, minwrk)
 135
 136         if lwork == -1 {
 137                 work[0] = float64(maxwrk)
 138                 return 0
 139         }
 140
 141         // Get machine constants.
 142         smlnum := math.Sqrt(dlamchS) / dlamchP
 143         bignum := 1 / smlnum
 144
 145         // Scale A if max element outside range [smlnum,bignum].
 146         anrm := impl.Dlange(lapack.MaxAbs, n, n, a, lda, nil)
 147         var scalea bool
 148         var cscale float64
 149         if 0 < anrm && anrm < smlnum {
 150                 scalea = true
 151                 cscale = smlnum
 152         } else if anrm > bignum {
 153                 scalea = true
 154                 cscale = bignum
 155         }
 156         if scalea {
 157                 impl.Dlascl(lapack.General, 0, 0, anrm, cscale, n, n, a, lda)
 158         }
 159
 160         // Balance the matrix.
 161         workbal := work[:n]
 162         ilo, ihi := impl.Dgebal(lapack.PermuteScale, n, a, lda, workbal)
 163
 164         // Reduce to upper Hessenberg form.
 165         iwrk := 2 * n
 166         tau := work[n : iwrk-1]
 167         impl.Dgehrd(n, ilo, ihi, a, lda, tau, work[iwrk:], lwork-iwrk)
 168
 169         var side lapack.EVSide
 170         if wantvl {
 171                 side = lapack.LeftEV
 172                 // Copy Householder vectors to VL.
 173                 impl.Dlacpy(blas.Lower, n, n, a, lda, vl, ldvl)
 174                 // Generate orthogonal matrix in VL.
 175                 impl.Dorghr(n, ilo, ihi, vl, ldvl, tau, work[iwrk:], lwork-iwrk)
 176                 // Perform QR iteration, accumulating Schur vectors in VL.
 177                 iwrk = n
 178                 first = impl.Dhseqr(lapack.EigenvaluesAndSchur, lapack.OriginalEV, n, ilo, ihi,
 179                         a, lda, wr, wi, vl, ldvl, work[iwrk:], lwork-iwrk)
 180                 if wantvr {
 181                         // Want left and right eigenvectors.
 182                         // Copy Schur vectors to VR.
 183                         side = lapack.RightLeftEV
 184                         impl.Dlacpy(blas.All, n, n, vl, ldvl, vr, ldvr)
 185                 }
 186         } else if wantvr {
 187                 side = lapack.RightEV
 188                 // Copy Householder vectors to VR.
 189                 impl.Dlacpy(blas.Lower, n, n, a, lda, vr, ldvr)
 190                 // Generate orthogonal matrix in VR.
 191                 impl.Dorghr(n, ilo, ihi, vr, ldvr, tau, work[iwrk:], lwork-iwrk)
 192                 // Perform QR iteration, accumulating Schur vectors in VR.
 193                 iwrk = n
 194                 first = impl.Dhseqr(lapack.EigenvaluesAndSchur, lapack.OriginalEV, n, ilo, ihi,
 195                         a, lda, wr, wi, vr, ldvr, work[iwrk:], lwork-iwrk)
 196         } else {
 197                 // Compute eigenvalues only.
 198                 iwrk = n
 199                 first = impl.Dhseqr(lapack.EigenvaluesOnly, lapack.None, n, ilo, ihi,
 200                         a, lda, wr, wi, nil, 1, work[iwrk:], lwork-iwrk)
 201         }
 202
 203         if first > 0 {
 204                 if scalea {
 205                         // Undo scaling.
 206                         impl.Dlascl(lapack.General, 0, 0, cscale, anrm, n-first, 1, wr[first:], 1)
 207                         impl.Dlascl(lapack.General, 0, 0, cscale, anrm, n-first, 1, wi[first:], 1)
 208                         impl.Dlascl(lapack.General, 0, 0, cscale, anrm, ilo, 1, wr, 1)
 209                         impl.Dlascl(lapack.General, 0, 0, cscale, anrm, ilo, 1, wi, 1)
 210                 }
 211                 work[0] = float64(maxwrk)
 212                 return first
 213         }
 214
 215         if wantvl || wantvr {
 216                 // Compute left and/or right eigenvectors.
 217                 impl.Dtrevc3(side, lapack.AllEVMulQ, nil, n,
 218                         a, lda, vl, ldvl, vr, ldvr, n, work[iwrk:], lwork-iwrk)
 219         }
 220         bi := blas64.Implementation()
 221         if wantvl {
 222                 // Undo balancing of left eigenvectors.
 223                 impl.Dgebak(lapack.PermuteScale, lapack.LeftEV, n, ilo, ihi, workbal, n, vl, ldvl)
 224                 // Normalize left eigenvectors and make largest component real.
 225                 for i, wii := range wi {
 226                         if wii < 0 {
 227                                 continue
 228                         }
 229                         if wii == 0 {
 230                                 scl := 1 / bi.Dnrm2(n, vl[i:], ldvl)
 231                                 bi.Dscal(n, scl, vl[i:], ldvl)
 232                                 continue
 233                         }
 234                         scl := 1 / impl.Dlapy2(bi.Dnrm2(n, vl[i:], ldvl), bi.Dnrm2(n, vl[i+1:], ldvl))
 235                         bi.Dscal(n, scl, vl[i:], ldvl)
 236                         bi.Dscal(n, scl, vl[i+1:], ldvl)
 237                         for k := 0; k < n; k++ {
 238                                 vi := vl[k*ldvl+i]
 239                                 vi1 := vl[k*ldvl+i+1]
 240                                 work[iwrk+k] = vi*vi + vi1*vi1
 241                         }
 242                         k := bi.Idamax(n, work[iwrk:iwrk+n], 1)
 243                         cs, sn, _ := impl.Dlartg(vl[k*ldvl+i], vl[k*ldvl+i+1])
 244                         bi.Drot(n, vl[i:], ldvl, vl[i+1:], ldvl, cs, sn)
 245                         vl[k*ldvl+i+1] = 0
 246                 }
 247         }
 248         if wantvr {
 249                 // Undo balancing of right eigenvectors.
 250                 impl.Dgebak(lapack.PermuteScale, lapack.RightEV, n, ilo, ihi, workbal, n, vr, ldvr)
 251                 // Normalize right eigenvectors and make largest component real.
 252                 for i, wii := range wi {
 253                         if wii < 0 {
 254                                 continue
 255                         }
 256                         if wii == 0 {
 257                                 scl := 1 / bi.Dnrm2(n, vr[i:], ldvr)
 258                                 bi.Dscal(n, scl, vr[i:], ldvr)
 259                                 continue
 260                         }
 261                         scl := 1 / impl.Dlapy2(bi.Dnrm2(n, vr[i:], ldvr), bi.Dnrm2(n, vr[i+1:], ldvr))
 262                         bi.Dscal(n, scl, vr[i:], ldvr)
 263                         bi.Dscal(n, scl, vr[i+1:], ldvr)
 264                         for k := 0; k < n; k++ {
 265                                 vi := vr[k*ldvr+i]
 266                                 vi1 := vr[k*ldvr+i+1]
 267                                 work[iwrk+k] = vi*vi + vi1*vi1
 268                         }
 269                         k := bi.Idamax(n, work[iwrk:iwrk+n], 1)
 270                         cs, sn, _ := impl.Dlartg(vr[k*ldvr+i], vr[k*ldvr+i+1])
 271                         bi.Drot(n, vr[i:], ldvr, vr[i+1:], ldvr, cs, sn)
 272                         vr[k*ldvr+i+1] = 0
 273                 }
 274         }
 275
 276         if scalea {
 277                 // Undo scaling.
 278                 impl.Dlascl(lapack.General, 0, 0, cscale, anrm, n-first, 1, wr[first:], 1)
 279                 impl.Dlascl(lapack.General, 0, 0, cscale, anrm, n-first, 1, wi[first:], 1)
 280         }
 281
 282         work[0] = float64(maxwrk)
 283         return first
 284 }