vendor/gonum.org/v1/gonum/lapack/gonum/dtrevc3.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 // Dtrevc3 computes some or all of the right and/or left eigenvectors of an n×n
  16 // upper quasi-triangular matrix T in Schur canonical form. Matrices of this
  17 // type are produced by the Schur factorization of a real general matrix A
  18 //  A = Q T Q^T,
  19 // as computed by Dhseqr.
  20 //
  21 // The right eigenvector x of T corresponding to an
  22 // eigenvalue λ is defined by
  23 //  T x = λ x,
  24 // and the left eigenvector is defined by
  25 //  y^H T = λ y^H,
  26 // where y^H is the conjugate transpose of y.
  27 //
  28 // The eigenvalues are read directly from the diagonal blocks of T.
  29 //
  30 // This routine returns the matrices X and/or Y of right and left eigenvectors
  31 // of T, or the products Q*X and/or Q*Y, where Q is an input matrix. If Q is the
  32 // orthogonal factor that reduces a matrix A to Schur form T, then Q*X and Q*Y
  33 // are the matrices of right and left eigenvectors of A.
  34 //
  35 // If side == lapack.RightEV, only right eigenvectors will be computed.
  36 // If side == lapack.LeftEV, only left eigenvectors will be computed.
  37 // If side == lapack.RightLeftEV, both right and left eigenvectors will be computed.
  38 // For other values of side, Dtrevc3 will panic.
  39 //
  40 // If howmny == lapack.AllEV, all right and/or left eigenvectors will be
  41 // computed.
  42 // If howmny == lapack.AllEVMulQ, all right and/or left eigenvectors will be
  43 // computed and multiplied from left by the matrices in VR and/or VL.
  44 // If howmny == lapack.SelectedEV, right and/or left eigenvectors will be
  45 // computed as indicated by selected.
  46 // For other values of howmny, Dtrevc3 will panic.
  47 //
  48 // selected specifies which eigenvectors will be computed. It must have length n
  49 // if howmny == lapack.SelectedEV, and it is not referenced otherwise.
  50 // If w_j is a real eigenvalue, the corresponding real eigenvector will be
  51 // computed if selected[j] is true.
  52 // If w_j and w_{j+1} are the real and imaginary parts of a complex eigenvalue,
  53 // the corresponding complex eigenvector is computed if either selected[j] or
  54 // selected[j+1] is true, and on return selected[j] will be set to true and
  55 // selected[j+1] will be set to false.
  56 //
  57 // VL and VR are n×mm matrices. If howmny is lapack.AllEV or
  58 // lapack.AllEVMulQ, mm must be at least n. If howmny ==
  59 // lapack.SelectedEV, mm must be large enough to store the selected
  60 // eigenvectors. Each selected real eigenvector occupies one column and each
  61 // selected complex eigenvector occupies two columns. If mm is not sufficiently
  62 // large, Dtrevc3 will panic.
  63 //
  64 // On entry, if howmny == lapack.AllEVMulQ, it is assumed that VL (if side
  65 // is lapack.LeftEV or lapack.RightLeftEV) contains an n×n matrix QL,
  66 // and that VR (if side is lapack.LeftEV or lapack.RightLeftEV) contains
  67 // an n×n matrix QR. QL and QR are typically the orthogonal matrix Q of Schur
  68 // vectors returned by Dhseqr.
  69 //
  70 // On return, if side is lapack.LeftEV or lapack.RightLeftEV,
  71 // VL will contain:
  72 //  if howmny == lapack.AllEV,      the matrix Y of left eigenvectors of T,
  73 //  if howmny == lapack.AllEVMulQ,  the matrix Q*Y,
  74 //  if howmny == lapack.SelectedEV, the left eigenvectors of T specified by
  75 //                                  selected, stored consecutively in the
  76 //                                  columns of VL, in the same order as their
  77 //                                  eigenvalues.
  78 // VL is not referenced if side == lapack.RightEV.
  79 //
  80 // On return, if side is lapack.RightEV or lapack.RightLeftEV,
  81 // VR will contain:
  82 //  if howmny == lapack.AllEV,      the matrix X of right eigenvectors of T,
  83 //  if howmny == lapack.AllEVMulQ,  the matrix Q*X,
  84 //  if howmny == lapack.SelectedEV, the left eigenvectors of T specified by
  85 //                                  selected, stored consecutively in the
  86 //                                  columns of VR, in the same order as their
  87 //                                  eigenvalues.
  88 // VR is not referenced if side == lapack.LeftEV.
  89 //
  90 // Complex eigenvectors corresponding to a complex eigenvalue are stored in VL
  91 // and VR in two consecutive columns, the first holding the real part, and the
  92 // second the imaginary part.
  93 //
  94 // Each eigenvector will be normalized so that the element of largest magnitude
  95 // has magnitude 1. Here the magnitude of a complex number (x,y) is taken to be
  96 // |x| + |y|.
  97 //
  98 // work must have length at least lwork and lwork must be at least max(1,3*n),
  99 // otherwise Dtrevc3 will panic. For optimum performance, lwork should be at
 100 // least n+2*n*nb, where nb is the optimal blocksize.
 101 //
 102 // If lwork == -1, instead of performing Dtrevc3, the function only estimates
 103 // the optimal workspace size based on n and stores it into work[0].
 104 //
 105 // Dtrevc3 returns the number of columns in VL and/or VR actually used to store
 106 // the eigenvectors.
 107 //
 108 // Dtrevc3 is an internal routine. It is exported for testing purposes.
 109 func (impl Implementation) Dtrevc3(side lapack.EVSide, howmny lapack.HowMany, selected []bool, n int, t []float64, ldt int, vl []float64, ldvl int, vr []float64, ldvr int, mm int, work []float64, lwork int) (m int) {
 110         switch side {
 111         default:
 112                 panic(badEVSide)
 113         case lapack.RightEV, lapack.LeftEV, lapack.RightLeftEV:
 114         }
 115         switch howmny {
 116         default:
 117                 panic(badHowMany)
 118         case lapack.AllEV, lapack.AllEVMulQ, lapack.SelectedEV:
 119         }
 120         switch {
 121         case n < 0:
 122                 panic(nLT0)
 123         case len(work) < lwork:
 124                 panic(shortWork)
 125         case lwork < max(1, 3*n) && lwork != -1:
 126                 panic(badWork)
 127         }
 128         if lwork != -1 {
 129                 if howmny == lapack.SelectedEV {
 130                         if len(selected) != n {
 131                                 panic("lapack: bad selected length")
 132                         }
 133                         // Set m to the number of columns required to store the
 134                         // selected eigenvectors, and standardize the slice
 135                         // selected.
 136                         for j := 0; j < n; {
 137                                 if j == n-1 || t[(j+1)*ldt+j] == 0 {
 138                                         // Diagonal 1×1 block corresponding to a
 139                                         // real eigenvalue.
 140                                         if selected[j] {
 141                                                 m++
 142                                         }
 143                                         j++
 144                                 } else {
 145                                         // Diagonal 2×2 block corresponding to a
 146                                         // complex eigenvalue.
 147                                         if selected[j] || selected[j+1] {
 148                                                 selected[j] = true
 149                                                 selected[j+1] = false
 150                                                 m += 2
 151                                         }
 152                                         j += 2
 153                                 }
 154                         }
 155                 } else {
 156                         m = n
 157                 }
 158                 if m > mm {
 159                         panic("lapack: insufficient number of columns")
 160                 }
 161                 checkMatrix(n, n, t, ldt)
 162                 if (side == lapack.RightEV || side == lapack.RightLeftEV) && m > 0 {
 163                         checkMatrix(n, m, vr, ldvr)
 164                 }
 165                 if (side == lapack.LeftEV || side == lapack.RightLeftEV) && m > 0 {
 166                         checkMatrix(n, m, vl, ldvl)
 167                 }
 168         }
 169
 170         // Quick return if possible.
 171         if n == 0 {
 172                 work[0] = 1
 173                 return m
 174         }
 175
 176         const (
 177                 nbmin = 8
 178                 nbmax = 128
 179         )
 180         nb := impl.Ilaenv(1, "DTREVC", string(side)+string(howmny), n, -1, -1, -1)
 181
 182         // Quick return in case of a workspace query.
 183         if lwork == -1 {
 184                 work[0] = float64(n + 2*n*nb)
 185                 return m
 186         }
 187
 188         // Use blocked version of back-transformation if sufficient workspace.
 189         // Zero-out the workspace to avoid potential NaN propagation.
 190         if howmny == lapack.AllEVMulQ && lwork >= n+2*n*nbmin {
 191                 nb = min((lwork-n)/(2*n), nbmax)
 192                 impl.Dlaset(blas.All, n, 1+2*nb, 0, 0, work[:n+2*nb*n], 1+2*nb)
 193         } else {
 194                 nb = 1
 195         }
 196
 197         // Set the constants to control overflow.
 198         ulp := dlamchP
 199         smlnum := float64(n) / ulp * dlamchS
 200         bignum := (1 - ulp) / smlnum
 201
 202         // Split work into a vector of column norms and an n×2*nb matrix b.
 203         norms := work[:n]
 204         ldb := 2 * nb
 205         b := work[n : n+n*ldb]
 206
 207         // Compute 1-norm of each column of strictly upper triangular part of T
 208         // to control overflow in triangular solver.
 209         norms[0] = 0
 210         for j := 1; j < n; j++ {
 211                 var cn float64
 212                 for i := 0; i < j; i++ {
 213                         cn += math.Abs(t[i*ldt+j])
 214                 }
 215                 norms[j] = cn
 216         }
 217
 218         bi := blas64.Implementation()
 219
 220         var (
 221                 x [4]float64
 222
 223                 iv int // Index of column in current block.
 224                 is int
 225
 226                 // ip is used below to specify the real or complex eigenvalue:
 227                 //  ip == 0, real eigenvalue,
 228                 //        1, first  of conjugate complex pair (wr,wi),
 229                 //       -1, second of conjugate complex pair (wr,wi).
 230                 ip        int
 231                 iscomplex [nbmax]int // Stores ip for each column in current block.
 232         )
 233
 234         if side == lapack.LeftEV {
 235                 goto leftev
 236         }
 237
 238         // Compute right eigenvectors.
 239
 240         // For complex right vector, iv-1 is for real part and iv for complex
 241         // part. Non-blocked version always uses iv=1, blocked version starts
 242         // with iv=nb-1 and goes down to 0 or 1.
 243         iv = max(2, nb) - 1
 244         ip = 0
 245         is = m - 1
 246         for ki := n - 1; ki >= 0; ki-- {
 247                 if ip == -1 {
 248                         // Previous iteration (ki+1) was second of
 249                         // conjugate pair, so this ki is first of
 250                         // conjugate pair.
 251                         ip = 1
 252                         continue
 253                 }
 254
 255                 if ki == 0 || t[ki*ldt+ki-1] == 0 {
 256                         // Last column or zero on sub-diagonal, so this
 257                         // ki must be real eigenvalue.
 258                         ip = 0
 259                 } else {
 260                         // Non-zero on sub-diagonal, so this ki is
 261                         // second of conjugate pair.
 262                         ip = -1
 263                 }
 264
 265                 if howmny == lapack.SelectedEV {
 266                         if ip == 0 {
 267                                 if !selected[ki] {
 268                                         continue
 269                                 }
 270                         } else if !selected[ki-1] {
 271                                 continue
 272                         }
 273                 }
 274
 275                 // Compute the ki-th eigenvalue (wr,wi).
 276                 wr := t[ki*ldt+ki]
 277                 var wi float64
 278                 if ip != 0 {
 279                         wi = math.Sqrt(math.Abs(t[ki*ldt+ki-1])) * math.Sqrt(math.Abs(t[(ki-1)*ldt+ki]))
 280                 }
 281                 smin := math.Max(ulp*(math.Abs(wr)+math.Abs(wi)), smlnum)
 282
 283                 if ip == 0 {
 284                         // Real right eigenvector.
 285
 286                         b[ki*ldb+iv] = 1
 287                         // Form right-hand side.
 288                         for k := 0; k < ki; k++ {
 289                                 b[k*ldb+iv] = -t[k*ldt+ki]
 290                         }
 291                         // Solve upper quasi-triangular system:
 292                         //  [ T[0:ki,0:ki] - wr ]*X = scale*b.
 293                         for j := ki - 1; j >= 0; {
 294                                 if j == 0 || t[j*ldt+j-1] == 0 {
 295                                         // 1×1 diagonal block.
 296                                         scale, xnorm, _ := impl.Dlaln2(false, 1, 1, smin, 1, t[j*ldt+j:], ldt,
 297                                                 1, 1, b[j*ldb+iv:], ldb, wr, 0, x[:1], 2)
 298                                         // Scale X[0,0] to avoid overflow when updating the
 299                                         // right-hand side.
 300                                         if xnorm > 1 && norms[j] > bignum/xnorm {
 301                                                 x[0] /= xnorm
 302                                                 scale /= xnorm
 303                                         }
 304                                         // Scale if necessary.
 305                                         if scale != 1 {
 306                                                 bi.Dscal(ki+1, scale, b[iv:], ldb)
 307                                         }
 308                                         b[j*ldb+iv] = x[0]
 309                                         // Update right-hand side.
 310                                         bi.Daxpy(j, -x[0], t[j:], ldt, b[iv:], ldb)
 311                                         j--
 312                                 } else {
 313                                         // 2×2 diagonal block.
 314                                         scale, xnorm, _ := impl.Dlaln2(false, 2, 1, smin, 1, t[(j-1)*ldt+j-1:], ldt,
 315                                                 1, 1, b[(j-1)*ldb+iv:], ldb, wr, 0, x[:3], 2)
 316                                         // Scale X[0,0] and X[1,0] to avoid overflow
 317                                         // when updating the right-hand side.
 318                                         if xnorm > 1 {
 319                                                 beta := math.Max(norms[j-1], norms[j])
 320                                                 if beta > bignum/xnorm {
 321                                                         x[0] /= xnorm
 322                                                         x[2] /= xnorm
 323                                                         scale /= xnorm
 324                                                 }
 325                                         }
 326                                         // Scale if necessary.
 327                                         if scale != 1 {
 328                                                 bi.Dscal(ki+1, scale, b[iv:], ldb)
 329                                         }
 330                                         b[(j-1)*ldb+iv] = x[0]
 331                                         b[j*ldb+iv] = x[2]
 332                                         // Update right-hand side.
 333                                         bi.Daxpy(j-1, -x[0], t[j-1:], ldt, b[iv:], ldb)
 334                                         bi.Daxpy(j-1, -x[2], t[j:], ldt, b[iv:], ldb)
 335                                         j -= 2
 336                                 }
 337                         }
 338                         // Copy the vector x or Q*x to VR and normalize.
 339                         switch {
 340                         case howmny != lapack.AllEVMulQ:
 341                                 // No back-transform: copy x to VR and normalize.
 342                                 bi.Dcopy(ki+1, b[iv:], ldb, vr[is:], ldvr)
 343                                 ii := bi.Idamax(ki+1, vr[is:], ldvr)
 344                                 remax := 1 / math.Abs(vr[ii*ldvr+is])
 345                                 bi.Dscal(ki+1, remax, vr[is:], ldvr)
 346                                 for k := ki + 1; k < n; k++ {
 347                                         vr[k*ldvr+is] = 0
 348                                 }
 349                         case nb == 1:
 350                                 // Version 1: back-transform each vector with GEMV, Q*x.
 351                                 if ki > 0 {
 352                                         bi.Dgemv(blas.NoTrans, n, ki, 1, vr, ldvr, b[iv:], ldb,
 353                                                 b[ki*ldb+iv], vr[ki:], ldvr)
 354                                 }
 355                                 ii := bi.Idamax(n, vr[ki:], ldvr)
 356                                 remax := 1 / math.Abs(vr[ii*ldvr+ki])
 357                                 bi.Dscal(n, remax, vr[ki:], ldvr)
 358                         default:
 359                                 // Version 2: back-transform block of vectors with GEMM.
 360                                 // Zero out below vector.
 361                                 for k := ki + 1; k < n; k++ {
 362                                         b[k*ldb+iv] = 0
 363                                 }
 364                                 iscomplex[iv] = ip
 365                                 // Back-transform and normalization is done below.
 366                         }
 367                 } else {
 368                         // Complex right eigenvector.
 369
 370                         // Initial solve
 371                         //  [ ( T[ki-1,ki-1] T[ki-1,ki] ) - (wr + i*wi) ]*X = 0.
 372                         //  [ ( T[ki,  ki-1] T[ki,  ki] )               ]
 373                         if math.Abs(t[(ki-1)*ldt+ki]) >= math.Abs(t[ki*ldt+ki-1]) {
 374                                 b[(ki-1)*ldb+iv-1] = 1
 375                                 b[ki*ldb+iv] = wi / t[(ki-1)*ldt+ki]
 376                         } else {
 377                                 b[(ki-1)*ldb+iv-1] = -wi / t[ki*ldt+ki-1]
 378                                 b[ki*ldb+iv] = 1
 379                         }
 380                         b[ki*ldb+iv-1] = 0
 381                         b[(ki-1)*ldb+iv] = 0
 382                         // Form right-hand side.
 383                         for k := 0; k < ki-1; k++ {
 384                                 b[k*ldb+iv-1] = -b[(ki-1)*ldb+iv-1] * t[k*ldt+ki-1]
 385                                 b[k*ldb+iv] = -b[ki*ldb+iv] * t[k*ldt+ki]
 386                         }
 387                         // Solve upper quasi-triangular system:
 388                         //  [ T[0:ki-1,0:ki-1] - (wr+i*wi) ]*X = scale*(b1+i*b2)
 389                         for j := ki - 2; j >= 0; {
 390                                 if j == 0 || t[j*ldt+j-1] == 0 {
 391                                         // 1×1 diagonal block.
 392
 393                                         scale, xnorm, _ := impl.Dlaln2(false, 1, 2, smin, 1, t[j*ldt+j:], ldt,
 394                                                 1, 1, b[j*ldb+iv-1:], ldb, wr, wi, x[:2], 2)
 395                                         // Scale X[0,0] and X[0,1] to avoid
 396                                         // overflow when updating the right-hand side.
 397                                         if xnorm > 1 && norms[j] > bignum/xnorm {
 398                                                 x[0] /= xnorm
 399                                                 x[1] /= xnorm
 400                                                 scale /= xnorm
 401                                         }
 402                                         // Scale if necessary.
 403                                         if scale != 1 {
 404                                                 bi.Dscal(ki+1, scale, b[iv-1:], ldb)
 405                                                 bi.Dscal(ki+1, scale, b[iv:], ldb)
 406                                         }
 407                                         b[j*ldb+iv-1] = x[0]
 408                                         b[j*ldb+iv] = x[1]
 409                                         // Update the right-hand side.
 410                                         bi.Daxpy(j, -x[0], t[j:], ldt, b[iv-1:], ldb)
 411                                         bi.Daxpy(j, -x[1], t[j:], ldt, b[iv:], ldb)
 412                                         j--
 413                                 } else {
 414                                         // 2×2 diagonal block.
 415
 416                                         scale, xnorm, _ := impl.Dlaln2(false, 2, 2, smin, 1, t[(j-1)*ldt+j-1:], ldt,
 417                                                 1, 1, b[(j-1)*ldb+iv-1:], ldb, wr, wi, x[:], 2)
 418                                         // Scale X to avoid overflow when updating
 419                                         // the right-hand side.
 420                                         if xnorm > 1 {
 421                                                 beta := math.Max(norms[j-1], norms[j])
 422                                                 if beta > bignum/xnorm {
 423                                                         rec := 1 / xnorm
 424                                                         x[0] *= rec
 425                                                         x[1] *= rec
 426                                                         x[2] *= rec
 427                                                         x[3] *= rec
 428                                                         scale *= rec
 429                                                 }
 430                                         }
 431                                         // Scale if necessary.
 432                                         if scale != 1 {
 433                                                 bi.Dscal(ki+1, scale, b[iv-1:], ldb)
 434                                                 bi.Dscal(ki+1, scale, b[iv:], ldb)
 435                                         }
 436                                         b[(j-1)*ldb+iv-1] = x[0]
 437                                         b[(j-1)*ldb+iv] = x[1]
 438                                         b[j*ldb+iv-1] = x[2]
 439                                         b[j*ldb+iv] = x[3]
 440                                         // Update the right-hand side.
 441                                         bi.Daxpy(j-1, -x[0], t[j-1:], ldt, b[iv-1:], ldb)
 442                                         bi.Daxpy(j-1, -x[1], t[j-1:], ldt, b[iv:], ldb)
 443                                         bi.Daxpy(j-1, -x[2], t[j:], ldt, b[iv-1:], ldb)
 444                                         bi.Daxpy(j-1, -x[3], t[j:], ldt, b[iv:], ldb)
 445                                         j -= 2
 446                                 }
 447                         }
 448
 449                         // Copy the vector x or Q*x to VR and normalize.
 450                         switch {
 451                         case howmny != lapack.AllEVMulQ:
 452                                 // No back-transform: copy x to VR and normalize.
 453                                 bi.Dcopy(ki+1, b[iv-1:], ldb, vr[is-1:], ldvr)
 454                                 bi.Dcopy(ki+1, b[iv:], ldb, vr[is:], ldvr)
 455                                 emax := 0.0
 456                                 for k := 0; k <= ki; k++ {
 457                                         emax = math.Max(emax, math.Abs(vr[k*ldvr+is-1])+math.Abs(vr[k*ldvr+is]))
 458                                 }
 459                                 remax := 1 / emax
 460                                 bi.Dscal(ki+1, remax, vr[is-1:], ldvr)
 461                                 bi.Dscal(ki+1, remax, vr[is:], ldvr)
 462                                 for k := ki + 1; k < n; k++ {
 463                                         vr[k*ldvr+is-1] = 0
 464                                         vr[k*ldvr+is] = 0
 465                                 }
 466                         case nb == 1:
 467                                 // Version 1: back-transform each vector with GEMV, Q*x.
 468                                 if ki-1 > 0 {
 469                                         bi.Dgemv(blas.NoTrans, n, ki-1, 1, vr, ldvr, b[iv-1:], ldb,
 470                                                 b[(ki-1)*ldb+iv-1], vr[ki-1:], ldvr)
 471                                         bi.Dgemv(blas.NoTrans, n, ki-1, 1, vr, ldvr, b[iv:], ldb,
 472                                                 b[ki*ldb+iv], vr[ki:], ldvr)
 473                                 } else {
 474                                         bi.Dscal(n, b[(ki-1)*ldb+iv-1], vr[ki-1:], ldvr)
 475                                         bi.Dscal(n, b[ki*ldb+iv], vr[ki:], ldvr)
 476                                 }
 477                                 emax := 0.0
 478                                 for k := 0; k < n; k++ {
 479                                         emax = math.Max(emax, math.Abs(vr[k*ldvr+ki-1])+math.Abs(vr[k*ldvr+ki]))
 480                                 }
 481                                 remax := 1 / emax
 482                                 bi.Dscal(n, remax, vr[ki-1:], ldvr)
 483                                 bi.Dscal(n, remax, vr[ki:], ldvr)
 484                         default:
 485                                 // Version 2: back-transform block of vectors with GEMM.
 486                                 // Zero out below vector.
 487                                 for k := ki + 1; k < n; k++ {
 488                                         b[k*ldb+iv-1] = 0
 489                                         b[k*ldb+iv] = 0
 490                                 }
 491                                 iscomplex[iv-1] = -ip
 492                                 iscomplex[iv] = ip
 493                                 iv--
 494                                 // Back-transform and normalization is done below.
 495                         }
 496                 }
 497                 if nb > 1 {
 498                         // Blocked version of back-transform.
 499
 500                         // For complex case, ki2 includes both vectors (ki-1 and ki).
 501                         ki2 := ki
 502                         if ip != 0 {
 503                                 ki2--
 504                         }
 505                         // Columns iv:nb of b are valid vectors.
 506                         // When the number of vectors stored reaches nb-1 or nb,
 507                         // or if this was last vector, do the Gemm.
 508                         if iv < 2 || ki2 == 0 {
 509                                 bi.Dgemm(blas.NoTrans, blas.NoTrans, n, nb-iv, ki2+nb-iv,
 510                                         1, vr, ldvr, b[iv:], ldb,
 511                                         0, b[nb+iv:], ldb)
 512                                 // Normalize vectors.
 513                                 var remax float64
 514                                 for k := iv; k < nb; k++ {
 515                                         if iscomplex[k] == 0 {
 516                                                 // Real eigenvector.
 517                                                 ii := bi.Idamax(n, b[nb+k:], ldb)
 518                                                 remax = 1 / math.Abs(b[ii*ldb+nb+k])
 519                                         } else if iscomplex[k] == 1 {
 520                                                 // First eigenvector of conjugate pair.
 521                                                 emax := 0.0
 522                                                 for ii := 0; ii < n; ii++ {
 523                                                         emax = math.Max(emax, math.Abs(b[ii*ldb+nb+k])+math.Abs(b[ii*ldb+nb+k+1]))
 524                                                 }
 525                                                 remax = 1 / emax
 526                                                 // Second eigenvector of conjugate pair
 527                                                 // will reuse this value of remax.
 528                                         }
 529                                         bi.Dscal(n, remax, b[nb+k:], ldb)
 530                                 }
 531                                 impl.Dlacpy(blas.All, n, nb-iv, b[nb+iv:], ldb, vr[ki2:], ldvr)
 532                                 iv = nb - 1
 533                         } else {
 534                                 iv--
 535                         }
 536                 }
 537                 is--
 538                 if ip != 0 {
 539                         is--
 540                 }
 541         }
 542
 543         if side == lapack.RightEV {
 544                 return m
 545         }
 546
 547 leftev:
 548         // Compute left eigenvectors.
 549
 550         // For complex left vector, iv is for real part and iv+1 for complex
 551         // part. Non-blocked version always uses iv=0. Blocked version starts
 552         // with iv=0, goes up to nb-2 or nb-1.
 553         iv = 0
 554         ip = 0
 555         is = 0
 556         for ki := 0; ki < n; ki++ {
 557                 if ip == 1 {
 558                         // Previous iteration ki-1 was first of conjugate pair,
 559                         // so this ki is second of conjugate pair.
 560                         ip = -1
 561                         continue
 562                 }
 563
 564                 if ki == n-1 || t[(ki+1)*ldt+ki] == 0 {
 565                         // Last column or zero on sub-diagonal, so this ki must
 566                         // be real eigenvalue.
 567                         ip = 0
 568                 } else {
 569                         // Non-zero on sub-diagonal, so this ki is first of
 570                         // conjugate pair.
 571                         ip = 1
 572                 }
 573                 if howmny == lapack.SelectedEV && !selected[ki] {
 574                         continue
 575                 }
 576
 577                 // Compute the ki-th eigenvalue (wr,wi).
 578                 wr := t[ki*ldt+ki]
 579                 var wi float64
 580                 if ip != 0 {
 581                         wi = math.Sqrt(math.Abs(t[ki*ldt+ki+1])) * math.Sqrt(math.Abs(t[(ki+1)*ldt+ki]))
 582                 }
 583                 smin := math.Max(ulp*(math.Abs(wr)+math.Abs(wi)), smlnum)
 584
 585                 if ip == 0 {
 586                         // Real left eigenvector.
 587
 588                         b[ki*ldb+iv] = 1
 589                         // Form right-hand side.
 590                         for k := ki + 1; k < n; k++ {
 591                                 b[k*ldb+iv] = -t[ki*ldt+k]
 592                         }
 593                         // Solve transposed quasi-triangular system:
 594                         //  [ T[ki+1:n,ki+1:n] - wr ]^T * X = scale*b
 595                         vmax := 1.0
 596                         vcrit := bignum
 597                         for j := ki + 1; j < n; {
 598                                 if j == n-1 || t[(j+1)*ldt+j] == 0 {
 599                                         // 1×1 diagonal block.
 600
 601                                         // Scale if necessary to avoid overflow
 602                                         // when forming the right-hand side.
 603                                         if norms[j] > vcrit {
 604                                                 rec := 1 / vmax
 605                                                 bi.Dscal(n-ki, rec, b[ki*ldb+iv:], ldb)
 606                                                 vmax = 1
 607                                                 vcrit = bignum
 608                                         }
 609                                         b[j*ldb+iv] -= bi.Ddot(j-ki-1, t[(ki+1)*ldt+j:], ldt, b[(ki+1)*ldb+iv:], ldb)
 610                                         // Solve [ T[j,j] - wr ]^T * X = b.
 611                                         scale, _, _ := impl.Dlaln2(false, 1, 1, smin, 1, t[j*ldt+j:], ldt,
 612                                                 1, 1, b[j*ldb+iv:], ldb, wr, 0, x[:1], 2)
 613                                         // Scale if necessary.
 614                                         if scale != 1 {
 615                                                 bi.Dscal(n-ki, scale, b[ki*ldb+iv:], ldb)
 616                                         }
 617                                         b[j*ldb+iv] = x[0]
 618                                         vmax = math.Max(math.Abs(b[j*ldb+iv]), vmax)
 619                                         vcrit = bignum / vmax
 620                                         j++
 621                                 } else {
 622                                         // 2×2 diagonal block.
 623
 624                                         // Scale if necessary to avoid overflow
 625                                         // when forming the right-hand side.
 626                                         beta := math.Max(norms[j], norms[j+1])
 627                                         if beta > vcrit {
 628                                                 bi.Dscal(n-ki+1, 1/vmax, b[ki*ldb+iv:], 1)
 629                                                 vmax = 1
 630                                                 vcrit = bignum
 631                                         }
 632                                         b[j*ldb+iv] -= bi.Ddot(j-ki-1, t[(ki+1)*ldt+j:], ldt, b[(ki+1)*ldb+iv:], ldb)
 633                                         b[(j+1)*ldb+iv] -= bi.Ddot(j-ki-1, t[(ki+1)*ldt+j+1:], ldt, b[(ki+1)*ldb+iv:], ldb)
 634                                         // Solve
 635                                         //  [ T[j,j]-wr  T[j,j+1]      ]^T * X = scale*[ b1 ]
 636                                         //  [ T[j+1,j]   T[j+1,j+1]-wr ]               [ b2 ]
 637                                         scale, _, _ := impl.Dlaln2(true, 2, 1, smin, 1, t[j*ldt+j:], ldt,
 638                                                 1, 1, b[j*ldb+iv:], ldb, wr, 0, x[:3], 2)
 639                                         // Scale if necessary.
 640                                         if scale != 1 {
 641                                                 bi.Dscal(n-ki, scale, b[ki*ldb+iv:], ldb)
 642                                         }
 643                                         b[j*ldb+iv] = x[0]
 644                                         b[(j+1)*ldb+iv] = x[2]
 645                                         vmax = math.Max(vmax, math.Max(math.Abs(b[j*ldb+iv]), math.Abs(b[(j+1)*ldb+iv])))
 646                                         vcrit = bignum / vmax
 647                                         j += 2
 648                                 }
 649                         }
 650                         // Copy the vector x or Q*x to VL and normalize.
 651                         switch {
 652                         case howmny != lapack.AllEVMulQ:
 653                                 // No back-transform: copy x to VL and normalize.
 654                                 bi.Dcopy(n-ki, b[ki*ldb+iv:], ldb, vl[ki*ldvl+is:], ldvl)
 655                                 ii := bi.Idamax(n-ki, vl[ki*ldvl+is:], ldvl) + ki
 656                                 remax := 1 / math.Abs(vl[ii*ldvl+is])
 657                                 bi.Dscal(n-ki, remax, vl[ki*ldvl+is:], ldvl)
 658                                 for k := 0; k < ki; k++ {
 659                                         vl[k*ldvl+is] = 0
 660                                 }
 661                         case nb == 1:
 662                                 // Version 1: back-transform each vector with Gemv, Q*x.
 663                                 if n-ki-1 > 0 {
 664                                         bi.Dgemv(blas.NoTrans, n, n-ki-1,
 665                                                 1, vl[ki+1:], ldvl, b[(ki+1)*ldb+iv:], ldb,
 666                                                 b[ki*ldb+iv], vl[ki:], ldvl)
 667                                 }
 668                                 ii := bi.Idamax(n, vl[ki:], ldvl)
 669                                 remax := 1 / math.Abs(vl[ii*ldvl+ki])
 670                                 bi.Dscal(n, remax, vl[ki:], ldvl)
 671                         default:
 672                                 // Version 2: back-transform block of vectors with Gemm
 673                                 // zero out above vector.
 674                                 for k := 0; k < ki; k++ {
 675                                         b[k*ldb+iv] = 0
 676                                 }
 677                                 iscomplex[iv] = ip
 678                                 // Back-transform and normalization is done below.
 679                         }
 680                 } else {
 681                         // Complex left eigenvector.
 682
 683                         // Initial solve:
 684                         // [ [ T[ki,ki]   T[ki,ki+1]   ]^T - (wr - i* wi) ]*X = 0.
 685                         // [ [ T[ki+1,ki] T[ki+1,ki+1] ]                  ]
 686                         if math.Abs(t[ki*ldt+ki+1]) >= math.Abs(t[(ki+1)*ldt+ki]) {
 687                                 b[ki*ldb+iv] = wi / t[ki*ldt+ki+1]
 688                                 b[(ki+1)*ldb+iv+1] = 1
 689                         } else {
 690                                 b[ki*ldb+iv] = 1
 691                                 b[(ki+1)*ldb+iv+1] = -wi / t[(ki+1)*ldt+ki]
 692                         }
 693                         b[(ki+1)*ldb+iv] = 0
 694                         b[ki*ldb+iv+1] = 0
 695                         // Form right-hand side.
 696                         for k := ki + 2; k < n; k++ {
 697                                 b[k*ldb+iv] = -b[ki*ldb+iv] * t[ki*ldt+k]
 698                                 b[k*ldb+iv+1] = -b[(ki+1)*ldb+iv+1] * t[(ki+1)*ldt+k]
 699                         }
 700                         // Solve transposed quasi-triangular system:
 701                         // [ T[ki+2:n,ki+2:n]^T - (wr-i*wi) ]*X = b1+i*b2
 702                         vmax := 1.0
 703                         vcrit := bignum
 704                         for j := ki + 2; j < n; {
 705                                 if j == n-1 || t[(j+1)*ldt+j] == 0 {
 706                                         // 1×1 diagonal block.
 707
 708                                         // Scale if necessary to avoid overflow
 709                                         // when forming the right-hand side elements.
 710                                         if norms[j] > vcrit {
 711                                                 rec := 1 / vmax
 712                                                 bi.Dscal(n-ki, rec, b[ki*ldb+iv:], ldb)
 713                                                 bi.Dscal(n-ki, rec, b[ki*ldb+iv+1:], ldb)
 714                                                 vmax = 1
 715                                                 vcrit = bignum
 716                                         }
 717                                         b[j*ldb+iv] -= bi.Ddot(j-ki-2, t[(ki+2)*ldt+j:], ldt, b[(ki+2)*ldb+iv:], ldb)
 718                                         b[j*ldb+iv+1] -= bi.Ddot(j-ki-2, t[(ki+2)*ldt+j:], ldt, b[(ki+2)*ldb+iv+1:], ldb)
 719                                         // Solve [ T[j,j]-(wr-i*wi) ]*(X11+i*X12) = b1+i*b2.
 720                                         scale, _, _ := impl.Dlaln2(false, 1, 2, smin, 1, t[j*ldt+j:], ldt,
 721                                                 1, 1, b[j*ldb+iv:], ldb, wr, -wi, x[:2], 2)
 722                                         // Scale if necessary.
 723                                         if scale != 1 {
 724                                                 bi.Dscal(n-ki, scale, b[ki*ldb+iv:], ldb)
 725                                                 bi.Dscal(n-ki, scale, b[ki*ldb+iv+1:], ldb)
 726                                         }
 727                                         b[j*ldb+iv] = x[0]
 728                                         b[j*ldb+iv+1] = x[1]
 729                                         vmax = math.Max(vmax, math.Max(math.Abs(b[j*ldb+iv]), math.Abs(b[j*ldb+iv+1])))
 730                                         vcrit = bignum / vmax
 731                                         j++
 732                                 } else {
 733                                         // 2×2 diagonal block.
 734
 735                                         // Scale if necessary to avoid overflow
 736                                         // when forming the right-hand side elements.
 737                                         if math.Max(norms[j], norms[j+1]) > vcrit {
 738                                                 rec := 1 / vmax
 739                                                 bi.Dscal(n-ki, rec, b[ki*ldb+iv:], ldb)
 740                                                 bi.Dscal(n-ki, rec, b[ki*ldb+iv+1:], ldb)
 741                                                 vmax = 1
 742                                                 vcrit = bignum
 743                                         }
 744                                         b[j*ldb+iv] -= bi.Ddot(j-ki-2, t[(ki+2)*ldt+j:], ldt, b[(ki+2)*ldb+iv:], ldb)
 745                                         b[j*ldb+iv+1] -= bi.Ddot(j-ki-2, t[(ki+2)*ldt+j:], ldt, b[(ki+2)*ldb+iv+1:], ldb)
 746                                         b[(j+1)*ldb+iv] -= bi.Ddot(j-ki-2, t[(ki+2)*ldt+j+1:], ldt, b[(ki+2)*ldb+iv:], ldb)
 747                                         b[(j+1)*ldb+iv+1] -= bi.Ddot(j-ki-2, t[(ki+2)*ldt+j+1:], ldt, b[(ki+2)*ldb+iv+1:], ldb)
 748                                         // Solve 2×2 complex linear equation
 749                                         //  [ [T[j,j]   T[j,j+1]  ]^T - (wr-i*wi)*I ]*X = scale*b
 750                                         //  [ [T[j+1,j] T[j+1,j+1]]                 ]
 751                                         scale, _, _ := impl.Dlaln2(true, 2, 2, smin, 1, t[j*ldt+j:], ldt,
 752                                                 1, 1, b[j*ldb+iv:], ldb, wr, -wi, x[:], 2)
 753                                         // Scale if necessary.
 754                                         if scale != 1 {
 755                                                 bi.Dscal(n-ki, scale, b[ki*ldb+iv:], ldb)
 756                                                 bi.Dscal(n-ki, scale, b[ki*ldb+iv+1:], ldb)
 757                                         }
 758                                         b[j*ldb+iv] = x[0]
 759                                         b[j*ldb+iv+1] = x[1]
 760                                         b[(j+1)*ldb+iv] = x[2]
 761                                         b[(j+1)*ldb+iv+1] = x[3]
 762                                         vmax01 := math.Max(math.Abs(x[0]), math.Abs(x[1]))
 763                                         vmax23 := math.Max(math.Abs(x[2]), math.Abs(x[3]))
 764                                         vmax = math.Max(vmax, math.Max(vmax01, vmax23))
 765                                         vcrit = bignum / vmax
 766                                         j += 2
 767                                 }
 768                         }
 769                         // Copy the vector x or Q*x to VL and normalize.
 770                         switch {
 771                         case howmny != lapack.AllEVMulQ:
 772                                 // No back-transform: copy x to VL and normalize.
 773                                 bi.Dcopy(n-ki, b[ki*ldb+iv:], ldb, vl[ki*ldvl+is:], ldvl)
 774                                 bi.Dcopy(n-ki, b[ki*ldb+iv+1:], ldb, vl[ki*ldvl+is+1:], ldvl)
 775                                 emax := 0.0
 776                                 for k := ki; k < n; k++ {
 777                                         emax = math.Max(emax, math.Abs(vl[k*ldvl+is])+math.Abs(vl[k*ldvl+is+1]))
 778                                 }
 779                                 remax := 1 / emax
 780                                 bi.Dscal(n-ki, remax, vl[ki*ldvl+is:], ldvl)
 781                                 bi.Dscal(n-ki, remax, vl[ki*ldvl+is+1:], ldvl)
 782                                 for k := 0; k < ki; k++ {
 783                                         vl[k*ldvl+is] = 0
 784                                         vl[k*ldvl+is+1] = 0
 785                                 }
 786                         case nb == 1:
 787                                 // Version 1: back-transform each vector with GEMV, Q*x.
 788                                 if n-ki-2 > 0 {
 789                                         bi.Dgemv(blas.NoTrans, n, n-ki-2,
 790                                                 1, vl[ki+2:], ldvl, b[(ki+2)*ldb+iv:], ldb,
 791                                                 b[ki*ldb+iv], vl[ki:], ldvl)
 792                                         bi.Dgemv(blas.NoTrans, n, n-ki-2,
 793                                                 1, vl[ki+2:], ldvl, b[(ki+2)*ldb+iv+1:], ldb,
 794                                                 b[(ki+1)*ldb+iv+1], vl[ki+1:], ldvl)
 795                                 } else {
 796                                         bi.Dscal(n, b[ki*ldb+iv], vl[ki:], ldvl)
 797                                         bi.Dscal(n, b[(ki+1)*ldb+iv+1], vl[ki+1:], ldvl)
 798                                 }
 799                                 emax := 0.0
 800                                 for k := 0; k < n; k++ {
 801                                         emax = math.Max(emax, math.Abs(vl[k*ldvl+ki])+math.Abs(vl[k*ldvl+ki+1]))
 802                                 }
 803                                 remax := 1 / emax
 804                                 bi.Dscal(n, remax, vl[ki:], ldvl)
 805                                 bi.Dscal(n, remax, vl[ki+1:], ldvl)
 806                         default:
 807                                 // Version 2: back-transform block of vectors with GEMM.
 808                                 // Zero out above vector.
 809                                 // Could go from ki-nv+1 to ki-1.
 810                                 for k := 0; k < ki; k++ {
 811                                         b[k*ldb+iv] = 0
 812                                         b[k*ldb+iv+1] = 0
 813                                 }
 814                                 iscomplex[iv] = ip
 815                                 iscomplex[iv+1] = -ip
 816                                 iv++
 817                                 // Back-transform and normalization is done below.
 818                         }
 819                 }
 820                 if nb > 1 {
 821                         // Blocked version of back-transform.
 822                         // For complex case, ki2 includes both vectors ki and ki+1.
 823                         ki2 := ki
 824                         if ip != 0 {
 825                                 ki2++
 826                         }
 827                         // Columns [0:iv] of work are valid vectors. When the
 828                         // number of vectors stored reaches nb-1 or nb, or if
 829                         // this was last vector, do the Gemm.
 830                         if iv >= nb-2 || ki2 == n-1 {
 831                                 bi.Dgemm(blas.NoTrans, blas.NoTrans, n, iv+1, n-ki2+iv,
 832                                         1, vl[ki2-iv:], ldvl, b[(ki2-iv)*ldb:], ldb,
 833                                         0, b[nb:], ldb)
 834                                 // Normalize vectors.
 835                                 var remax float64
 836                                 for k := 0; k <= iv; k++ {
 837                                         if iscomplex[k] == 0 {
 838                                                 // Real eigenvector.
 839                                                 ii := bi.Idamax(n, b[nb+k:], ldb)
 840                                                 remax = 1 / math.Abs(b[ii*ldb+nb+k])
 841                                         } else if iscomplex[k] == 1 {
 842                                                 // First eigenvector of conjugate pair.
 843                                                 emax := 0.0
 844                                                 for ii := 0; ii < n; ii++ {
 845                                                         emax = math.Max(emax, math.Abs(b[ii*ldb+nb+k])+math.Abs(b[ii*ldb+nb+k+1]))
 846                                                 }
 847                                                 remax = 1 / emax
 848                                                 // Second eigenvector of conjugate pair
 849                                                 // will reuse this value of remax.
 850                                         }
 851                                         bi.Dscal(n, remax, b[nb+k:], ldb)
 852                                 }
 853                                 impl.Dlacpy(blas.All, n, iv+1, b[nb:], ldb, vl[ki2-iv:], ldvl)
 854                                 iv = 0
 855                         } else {
 856                                 iv++
 857                         }
 858                 }
 859                 is++
 860                 if ip != 0 {
 861                         is++
 862                 }
 863         }
 864
 865         return m
 866 }