| // 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" |
| ) |
| |
| // Dtrevc3 computes some or all of the right and/or left eigenvectors of an n×n |
| // upper quasi-triangular matrix T in Schur canonical form. Matrices of this |
| // type are produced by the Schur factorization of a real general matrix A |
| // A = Q T Q^T, |
| // as computed by Dhseqr. |
| // |
| // The right eigenvector x of T corresponding to an |
| // eigenvalue λ is defined by |
| // T x = λ x, |
| // and the left eigenvector is defined by |
| // y^H T = λ y^H, |
| // where y^H is the conjugate transpose of y. |
| // |
| // The eigenvalues are read directly from the diagonal blocks of T. |
| // |
| // This routine returns the matrices X and/or Y of right and left eigenvectors |
| // of T, or the products Q*X and/or Q*Y, where Q is an input matrix. If Q is the |
| // orthogonal factor that reduces a matrix A to Schur form T, then Q*X and Q*Y |
| // are the matrices of right and left eigenvectors of A. |
| // |
| // If side == lapack.RightEV, only right eigenvectors will be computed. |
| // If side == lapack.LeftEV, only left eigenvectors will be computed. |
| // If side == lapack.RightLeftEV, both right and left eigenvectors will be computed. |
| // For other values of side, Dtrevc3 will panic. |
| // |
| // If howmny == lapack.AllEV, all right and/or left eigenvectors will be |
| // computed. |
| // If howmny == lapack.AllEVMulQ, all right and/or left eigenvectors will be |
| // computed and multiplied from left by the matrices in VR and/or VL. |
| // If howmny == lapack.SelectedEV, right and/or left eigenvectors will be |
| // computed as indicated by selected. |
| // For other values of howmny, Dtrevc3 will panic. |
| // |
| // selected specifies which eigenvectors will be computed. It must have length n |
| // if howmny == lapack.SelectedEV, and it is not referenced otherwise. |
| // If w_j is a real eigenvalue, the corresponding real eigenvector will be |
| // computed if selected[j] is true. |
| // If w_j and w_{j+1} are the real and imaginary parts of a complex eigenvalue, |
| // the corresponding complex eigenvector is computed if either selected[j] or |
| // selected[j+1] is true, and on return selected[j] will be set to true and |
| // selected[j+1] will be set to false. |
| // |
| // VL and VR are n×mm matrices. If howmny is lapack.AllEV or |
| // lapack.AllEVMulQ, mm must be at least n. If howmny == |
| // lapack.SelectedEV, mm must be large enough to store the selected |
| // eigenvectors. Each selected real eigenvector occupies one column and each |
| // selected complex eigenvector occupies two columns. If mm is not sufficiently |
| // large, Dtrevc3 will panic. |
| // |
| // On entry, if howmny == lapack.AllEVMulQ, it is assumed that VL (if side |
| // is lapack.LeftEV or lapack.RightLeftEV) contains an n×n matrix QL, |
| // and that VR (if side is lapack.LeftEV or lapack.RightLeftEV) contains |
| // an n×n matrix QR. QL and QR are typically the orthogonal matrix Q of Schur |
| // vectors returned by Dhseqr. |
| // |
| // On return, if side is lapack.LeftEV or lapack.RightLeftEV, |
| // VL will contain: |
| // if howmny == lapack.AllEV, the matrix Y of left eigenvectors of T, |
| // if howmny == lapack.AllEVMulQ, the matrix Q*Y, |
| // if howmny == lapack.SelectedEV, the left eigenvectors of T specified by |
| // selected, stored consecutively in the |
| // columns of VL, in the same order as their |
| // eigenvalues. |
| // VL is not referenced if side == lapack.RightEV. |
| // |
| // On return, if side is lapack.RightEV or lapack.RightLeftEV, |
| // VR will contain: |
| // if howmny == lapack.AllEV, the matrix X of right eigenvectors of T, |
| // if howmny == lapack.AllEVMulQ, the matrix Q*X, |
| // if howmny == lapack.SelectedEV, the left eigenvectors of T specified by |
| // selected, stored consecutively in the |
| // columns of VR, in the same order as their |
| // eigenvalues. |
| // VR is not referenced if side == lapack.LeftEV. |
| // |
| // Complex eigenvectors corresponding to a complex eigenvalue are stored in VL |
| // and VR in two consecutive columns, the first holding the real part, and the |
| // second the imaginary part. |
| // |
| // Each eigenvector will be normalized so that the element of largest magnitude |
| // has magnitude 1. Here the magnitude of a complex number (x,y) is taken to be |
| // |x| + |y|. |
| // |
| // work must have length at least lwork and lwork must be at least max(1,3*n), |
| // otherwise Dtrevc3 will panic. For optimum performance, lwork should be at |
| // least n+2*n*nb, where nb is the optimal blocksize. |
| // |
| // If lwork == -1, instead of performing Dtrevc3, the function only estimates |
| // the optimal workspace size based on n and stores it into work[0]. |
| // |
| // Dtrevc3 returns the number of columns in VL and/or VR actually used to store |
| // the eigenvectors. |
| // |
| // Dtrevc3 is an internal routine. It is exported for testing purposes. |
| 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) { |
| switch side { |
| default: |
| panic(badEVSide) |
| case lapack.RightEV, lapack.LeftEV, lapack.RightLeftEV: |
| } |
| switch howmny { |
| default: |
| panic(badHowMany) |
| case lapack.AllEV, lapack.AllEVMulQ, lapack.SelectedEV: |
| } |
| switch { |
| case n < 0: |
| panic(nLT0) |
| case len(work) < lwork: |
| panic(shortWork) |
| case lwork < max(1, 3*n) && lwork != -1: |
| panic(badWork) |
| } |
| if lwork != -1 { |
| if howmny == lapack.SelectedEV { |
| if len(selected) != n { |
| panic("lapack: bad selected length") |
| } |
| // Set m to the number of columns required to store the |
| // selected eigenvectors, and standardize the slice |
| // selected. |
| for j := 0; j < n; { |
| if j == n-1 || t[(j+1)*ldt+j] == 0 { |
| // Diagonal 1×1 block corresponding to a |
| // real eigenvalue. |
| if selected[j] { |
| m++ |
| } |
| j++ |
| } else { |
| // Diagonal 2×2 block corresponding to a |
| // complex eigenvalue. |
| if selected[j] || selected[j+1] { |
| selected[j] = true |
| selected[j+1] = false |
| m += 2 |
| } |
| j += 2 |
| } |
| } |
| } else { |
| m = n |
| } |
| if m > mm { |
| panic("lapack: insufficient number of columns") |
| } |
| checkMatrix(n, n, t, ldt) |
| if (side == lapack.RightEV || side == lapack.RightLeftEV) && m > 0 { |
| checkMatrix(n, m, vr, ldvr) |
| } |
| if (side == lapack.LeftEV || side == lapack.RightLeftEV) && m > 0 { |
| checkMatrix(n, m, vl, ldvl) |
| } |
| } |
| |
| // Quick return if possible. |
| if n == 0 { |
| work[0] = 1 |
| return m |
| } |
| |
| const ( |
| nbmin = 8 |
| nbmax = 128 |
| ) |
| nb := impl.Ilaenv(1, "DTREVC", string(side)+string(howmny), n, -1, -1, -1) |
| |
| // Quick return in case of a workspace query. |
| if lwork == -1 { |
| work[0] = float64(n + 2*n*nb) |
| return m |
| } |
| |
| // Use blocked version of back-transformation if sufficient workspace. |
| // Zero-out the workspace to avoid potential NaN propagation. |
| if howmny == lapack.AllEVMulQ && lwork >= n+2*n*nbmin { |
| nb = min((lwork-n)/(2*n), nbmax) |
| impl.Dlaset(blas.All, n, 1+2*nb, 0, 0, work[:n+2*nb*n], 1+2*nb) |
| } else { |
| nb = 1 |
| } |
| |
| // Set the constants to control overflow. |
| ulp := dlamchP |
| smlnum := float64(n) / ulp * dlamchS |
| bignum := (1 - ulp) / smlnum |
| |
| // Split work into a vector of column norms and an n×2*nb matrix b. |
| norms := work[:n] |
| ldb := 2 * nb |
| b := work[n : n+n*ldb] |
| |
| // Compute 1-norm of each column of strictly upper triangular part of T |
| // to control overflow in triangular solver. |
| norms[0] = 0 |
| for j := 1; j < n; j++ { |
| var cn float64 |
| for i := 0; i < j; i++ { |
| cn += math.Abs(t[i*ldt+j]) |
| } |
| norms[j] = cn |
| } |
| |
| bi := blas64.Implementation() |
| |
| var ( |
| x [4]float64 |
| |
| iv int // Index of column in current block. |
| is int |
| |
| // ip is used below to specify the real or complex eigenvalue: |
| // ip == 0, real eigenvalue, |
| // 1, first of conjugate complex pair (wr,wi), |
| // -1, second of conjugate complex pair (wr,wi). |
| ip int |
| iscomplex [nbmax]int // Stores ip for each column in current block. |
| ) |
| |
| if side == lapack.LeftEV { |
| goto leftev |
| } |
| |
| // Compute right eigenvectors. |
| |
| // For complex right vector, iv-1 is for real part and iv for complex |
| // part. Non-blocked version always uses iv=1, blocked version starts |
| // with iv=nb-1 and goes down to 0 or 1. |
| iv = max(2, nb) - 1 |
| ip = 0 |
| is = m - 1 |
| for ki := n - 1; ki >= 0; ki-- { |
| if ip == -1 { |
| // Previous iteration (ki+1) was second of |
| // conjugate pair, so this ki is first of |
| // conjugate pair. |
| ip = 1 |
| continue |
| } |
| |
| if ki == 0 || t[ki*ldt+ki-1] == 0 { |
| // Last column or zero on sub-diagonal, so this |
| // ki must be real eigenvalue. |
| ip = 0 |
| } else { |
| // Non-zero on sub-diagonal, so this ki is |
| // second of conjugate pair. |
| ip = -1 |
| } |
| |
| if howmny == lapack.SelectedEV { |
| if ip == 0 { |
| if !selected[ki] { |
| continue |
| } |
| } else if !selected[ki-1] { |
| continue |
| } |
| } |
| |
| // Compute the ki-th eigenvalue (wr,wi). |
| wr := t[ki*ldt+ki] |
| var wi float64 |
| if ip != 0 { |
| wi = math.Sqrt(math.Abs(t[ki*ldt+ki-1])) * math.Sqrt(math.Abs(t[(ki-1)*ldt+ki])) |
| } |
| smin := math.Max(ulp*(math.Abs(wr)+math.Abs(wi)), smlnum) |
| |
| if ip == 0 { |
| // Real right eigenvector. |
| |
| b[ki*ldb+iv] = 1 |
| // Form right-hand side. |
| for k := 0; k < ki; k++ { |
| b[k*ldb+iv] = -t[k*ldt+ki] |
| } |
| // Solve upper quasi-triangular system: |
| // [ T[0:ki,0:ki] - wr ]*X = scale*b. |
| for j := ki - 1; j >= 0; { |
| if j == 0 || t[j*ldt+j-1] == 0 { |
| // 1×1 diagonal block. |
| scale, xnorm, _ := impl.Dlaln2(false, 1, 1, smin, 1, t[j*ldt+j:], ldt, |
| 1, 1, b[j*ldb+iv:], ldb, wr, 0, x[:1], 2) |
| // Scale X[0,0] to avoid overflow when updating the |
| // right-hand side. |
| if xnorm > 1 && norms[j] > bignum/xnorm { |
| x[0] /= xnorm |
| scale /= xnorm |
| } |
| // Scale if necessary. |
| if scale != 1 { |
| bi.Dscal(ki+1, scale, b[iv:], ldb) |
| } |
| b[j*ldb+iv] = x[0] |
| // Update right-hand side. |
| bi.Daxpy(j, -x[0], t[j:], ldt, b[iv:], ldb) |
| j-- |
| } else { |
| // 2×2 diagonal block. |
| scale, xnorm, _ := impl.Dlaln2(false, 2, 1, smin, 1, t[(j-1)*ldt+j-1:], ldt, |
| 1, 1, b[(j-1)*ldb+iv:], ldb, wr, 0, x[:3], 2) |
| // Scale X[0,0] and X[1,0] to avoid overflow |
| // when updating the right-hand side. |
| if xnorm > 1 { |
| beta := math.Max(norms[j-1], norms[j]) |
| if beta > bignum/xnorm { |
| x[0] /= xnorm |
| x[2] /= xnorm |
| scale /= xnorm |
| } |
| } |
| // Scale if necessary. |
| if scale != 1 { |
| bi.Dscal(ki+1, scale, b[iv:], ldb) |
| } |
| b[(j-1)*ldb+iv] = x[0] |
| b[j*ldb+iv] = x[2] |
| // Update right-hand side. |
| bi.Daxpy(j-1, -x[0], t[j-1:], ldt, b[iv:], ldb) |
| bi.Daxpy(j-1, -x[2], t[j:], ldt, b[iv:], ldb) |
| j -= 2 |
| } |
| } |
| // Copy the vector x or Q*x to VR and normalize. |
| switch { |
| case howmny != lapack.AllEVMulQ: |
| // No back-transform: copy x to VR and normalize. |
| bi.Dcopy(ki+1, b[iv:], ldb, vr[is:], ldvr) |
| ii := bi.Idamax(ki+1, vr[is:], ldvr) |
| remax := 1 / math.Abs(vr[ii*ldvr+is]) |
| bi.Dscal(ki+1, remax, vr[is:], ldvr) |
| for k := ki + 1; k < n; k++ { |
| vr[k*ldvr+is] = 0 |
| } |
| case nb == 1: |
| // Version 1: back-transform each vector with GEMV, Q*x. |
| if ki > 0 { |
| bi.Dgemv(blas.NoTrans, n, ki, 1, vr, ldvr, b[iv:], ldb, |
| b[ki*ldb+iv], vr[ki:], ldvr) |
| } |
| ii := bi.Idamax(n, vr[ki:], ldvr) |
| remax := 1 / math.Abs(vr[ii*ldvr+ki]) |
| bi.Dscal(n, remax, vr[ki:], ldvr) |
| default: |
| // Version 2: back-transform block of vectors with GEMM. |
| // Zero out below vector. |
| for k := ki + 1; k < n; k++ { |
| b[k*ldb+iv] = 0 |
| } |
| iscomplex[iv] = ip |
| // Back-transform and normalization is done below. |
| } |
| } else { |
| // Complex right eigenvector. |
| |
| // Initial solve |
| // [ ( T[ki-1,ki-1] T[ki-1,ki] ) - (wr + i*wi) ]*X = 0. |
| // [ ( T[ki, ki-1] T[ki, ki] ) ] |
| if math.Abs(t[(ki-1)*ldt+ki]) >= math.Abs(t[ki*ldt+ki-1]) { |
| b[(ki-1)*ldb+iv-1] = 1 |
| b[ki*ldb+iv] = wi / t[(ki-1)*ldt+ki] |
| } else { |
| b[(ki-1)*ldb+iv-1] = -wi / t[ki*ldt+ki-1] |
| b[ki*ldb+iv] = 1 |
| } |
| b[ki*ldb+iv-1] = 0 |
| b[(ki-1)*ldb+iv] = 0 |
| // Form right-hand side. |
| for k := 0; k < ki-1; k++ { |
| b[k*ldb+iv-1] = -b[(ki-1)*ldb+iv-1] * t[k*ldt+ki-1] |
| b[k*ldb+iv] = -b[ki*ldb+iv] * t[k*ldt+ki] |
| } |
| // Solve upper quasi-triangular system: |
| // [ T[0:ki-1,0:ki-1] - (wr+i*wi) ]*X = scale*(b1+i*b2) |
| for j := ki - 2; j >= 0; { |
| if j == 0 || t[j*ldt+j-1] == 0 { |
| // 1×1 diagonal block. |
| |
| scale, xnorm, _ := impl.Dlaln2(false, 1, 2, smin, 1, t[j*ldt+j:], ldt, |
| 1, 1, b[j*ldb+iv-1:], ldb, wr, wi, x[:2], 2) |
| // Scale X[0,0] and X[0,1] to avoid |
| // overflow when updating the right-hand side. |
| if xnorm > 1 && norms[j] > bignum/xnorm { |
| x[0] /= xnorm |
| x[1] /= xnorm |
| scale /= xnorm |
| } |
| // Scale if necessary. |
| if scale != 1 { |
| bi.Dscal(ki+1, scale, b[iv-1:], ldb) |
| bi.Dscal(ki+1, scale, b[iv:], ldb) |
| } |
| b[j*ldb+iv-1] = x[0] |
| b[j*ldb+iv] = x[1] |
| // Update the right-hand side. |
| bi.Daxpy(j, -x[0], t[j:], ldt, b[iv-1:], ldb) |
| bi.Daxpy(j, -x[1], t[j:], ldt, b[iv:], ldb) |
| j-- |
| } else { |
| // 2×2 diagonal block. |
| |
| scale, xnorm, _ := impl.Dlaln2(false, 2, 2, smin, 1, t[(j-1)*ldt+j-1:], ldt, |
| 1, 1, b[(j-1)*ldb+iv-1:], ldb, wr, wi, x[:], 2) |
| // Scale X to avoid overflow when updating |
| // the right-hand side. |
| if xnorm > 1 { |
| beta := math.Max(norms[j-1], norms[j]) |
| if beta > bignum/xnorm { |
| rec := 1 / xnorm |
| x[0] *= rec |
| x[1] *= rec |
| x[2] *= rec |
| x[3] *= rec |
| scale *= rec |
| } |
| } |
| // Scale if necessary. |
| if scale != 1 { |
| bi.Dscal(ki+1, scale, b[iv-1:], ldb) |
| bi.Dscal(ki+1, scale, b[iv:], ldb) |
| } |
| b[(j-1)*ldb+iv-1] = x[0] |
| b[(j-1)*ldb+iv] = x[1] |
| b[j*ldb+iv-1] = x[2] |
| b[j*ldb+iv] = x[3] |
| // Update the right-hand side. |
| bi.Daxpy(j-1, -x[0], t[j-1:], ldt, b[iv-1:], ldb) |
| bi.Daxpy(j-1, -x[1], t[j-1:], ldt, b[iv:], ldb) |
| bi.Daxpy(j-1, -x[2], t[j:], ldt, b[iv-1:], ldb) |
| bi.Daxpy(j-1, -x[3], t[j:], ldt, b[iv:], ldb) |
| j -= 2 |
| } |
| } |
| |
| // Copy the vector x or Q*x to VR and normalize. |
| switch { |
| case howmny != lapack.AllEVMulQ: |
| // No back-transform: copy x to VR and normalize. |
| bi.Dcopy(ki+1, b[iv-1:], ldb, vr[is-1:], ldvr) |
| bi.Dcopy(ki+1, b[iv:], ldb, vr[is:], ldvr) |
| emax := 0.0 |
| for k := 0; k <= ki; k++ { |
| emax = math.Max(emax, math.Abs(vr[k*ldvr+is-1])+math.Abs(vr[k*ldvr+is])) |
| } |
| remax := 1 / emax |
| bi.Dscal(ki+1, remax, vr[is-1:], ldvr) |
| bi.Dscal(ki+1, remax, vr[is:], ldvr) |
| for k := ki + 1; k < n; k++ { |
| vr[k*ldvr+is-1] = 0 |
| vr[k*ldvr+is] = 0 |
| } |
| case nb == 1: |
| // Version 1: back-transform each vector with GEMV, Q*x. |
| if ki-1 > 0 { |
| bi.Dgemv(blas.NoTrans, n, ki-1, 1, vr, ldvr, b[iv-1:], ldb, |
| b[(ki-1)*ldb+iv-1], vr[ki-1:], ldvr) |
| bi.Dgemv(blas.NoTrans, n, ki-1, 1, vr, ldvr, b[iv:], ldb, |
| b[ki*ldb+iv], vr[ki:], ldvr) |
| } else { |
| bi.Dscal(n, b[(ki-1)*ldb+iv-1], vr[ki-1:], ldvr) |
| bi.Dscal(n, b[ki*ldb+iv], vr[ki:], ldvr) |
| } |
| emax := 0.0 |
| for k := 0; k < n; k++ { |
| emax = math.Max(emax, math.Abs(vr[k*ldvr+ki-1])+math.Abs(vr[k*ldvr+ki])) |
| } |
| remax := 1 / emax |
| bi.Dscal(n, remax, vr[ki-1:], ldvr) |
| bi.Dscal(n, remax, vr[ki:], ldvr) |
| default: |
| // Version 2: back-transform block of vectors with GEMM. |
| // Zero out below vector. |
| for k := ki + 1; k < n; k++ { |
| b[k*ldb+iv-1] = 0 |
| b[k*ldb+iv] = 0 |
| } |
| iscomplex[iv-1] = -ip |
| iscomplex[iv] = ip |
| iv-- |
| // Back-transform and normalization is done below. |
| } |
| } |
| if nb > 1 { |
| // Blocked version of back-transform. |
| |
| // For complex case, ki2 includes both vectors (ki-1 and ki). |
| ki2 := ki |
| if ip != 0 { |
| ki2-- |
| } |
| // Columns iv:nb of b are valid vectors. |
| // When the number of vectors stored reaches nb-1 or nb, |
| // or if this was last vector, do the Gemm. |
| if iv < 2 || ki2 == 0 { |
| bi.Dgemm(blas.NoTrans, blas.NoTrans, n, nb-iv, ki2+nb-iv, |
| 1, vr, ldvr, b[iv:], ldb, |
| 0, b[nb+iv:], ldb) |
| // Normalize vectors. |
| var remax float64 |
| for k := iv; k < nb; k++ { |
| if iscomplex[k] == 0 { |
| // Real eigenvector. |
| ii := bi.Idamax(n, b[nb+k:], ldb) |
| remax = 1 / math.Abs(b[ii*ldb+nb+k]) |
| } else if iscomplex[k] == 1 { |
| // First eigenvector of conjugate pair. |
| emax := 0.0 |
| for ii := 0; ii < n; ii++ { |
| emax = math.Max(emax, math.Abs(b[ii*ldb+nb+k])+math.Abs(b[ii*ldb+nb+k+1])) |
| } |
| remax = 1 / emax |
| // Second eigenvector of conjugate pair |
| // will reuse this value of remax. |
| } |
| bi.Dscal(n, remax, b[nb+k:], ldb) |
| } |
| impl.Dlacpy(blas.All, n, nb-iv, b[nb+iv:], ldb, vr[ki2:], ldvr) |
| iv = nb - 1 |
| } else { |
| iv-- |
| } |
| } |
| is-- |
| if ip != 0 { |
| is-- |
| } |
| } |
| |
| if side == lapack.RightEV { |
| return m |
| } |
| |
| leftev: |
| // Compute left eigenvectors. |
| |
| // For complex left vector, iv is for real part and iv+1 for complex |
| // part. Non-blocked version always uses iv=0. Blocked version starts |
| // with iv=0, goes up to nb-2 or nb-1. |
| iv = 0 |
| ip = 0 |
| is = 0 |
| for ki := 0; ki < n; ki++ { |
| if ip == 1 { |
| // Previous iteration ki-1 was first of conjugate pair, |
| // so this ki is second of conjugate pair. |
| ip = -1 |
| continue |
| } |
| |
| if ki == n-1 || t[(ki+1)*ldt+ki] == 0 { |
| // Last column or zero on sub-diagonal, so this ki must |
| // be real eigenvalue. |
| ip = 0 |
| } else { |
| // Non-zero on sub-diagonal, so this ki is first of |
| // conjugate pair. |
| ip = 1 |
| } |
| if howmny == lapack.SelectedEV && !selected[ki] { |
| continue |
| } |
| |
| // Compute the ki-th eigenvalue (wr,wi). |
| wr := t[ki*ldt+ki] |
| var wi float64 |
| if ip != 0 { |
| wi = math.Sqrt(math.Abs(t[ki*ldt+ki+1])) * math.Sqrt(math.Abs(t[(ki+1)*ldt+ki])) |
| } |
| smin := math.Max(ulp*(math.Abs(wr)+math.Abs(wi)), smlnum) |
| |
| if ip == 0 { |
| // Real left eigenvector. |
| |
| b[ki*ldb+iv] = 1 |
| // Form right-hand side. |
| for k := ki + 1; k < n; k++ { |
| b[k*ldb+iv] = -t[ki*ldt+k] |
| } |
| // Solve transposed quasi-triangular system: |
| // [ T[ki+1:n,ki+1:n] - wr ]^T * X = scale*b |
| vmax := 1.0 |
| vcrit := bignum |
| for j := ki + 1; j < n; { |
| if j == n-1 || t[(j+1)*ldt+j] == 0 { |
| // 1×1 diagonal block. |
| |
| // Scale if necessary to avoid overflow |
| // when forming the right-hand side. |
| if norms[j] > vcrit { |
| rec := 1 / vmax |
| bi.Dscal(n-ki, rec, b[ki*ldb+iv:], ldb) |
| vmax = 1 |
| vcrit = bignum |
| } |
| b[j*ldb+iv] -= bi.Ddot(j-ki-1, t[(ki+1)*ldt+j:], ldt, b[(ki+1)*ldb+iv:], ldb) |
| // Solve [ T[j,j] - wr ]^T * X = b. |
| scale, _, _ := impl.Dlaln2(false, 1, 1, smin, 1, t[j*ldt+j:], ldt, |
| 1, 1, b[j*ldb+iv:], ldb, wr, 0, x[:1], 2) |
| // Scale if necessary. |
| if scale != 1 { |
| bi.Dscal(n-ki, scale, b[ki*ldb+iv:], ldb) |
| } |
| b[j*ldb+iv] = x[0] |
| vmax = math.Max(math.Abs(b[j*ldb+iv]), vmax) |
| vcrit = bignum / vmax |
| j++ |
| } else { |
| // 2×2 diagonal block. |
| |
| // Scale if necessary to avoid overflow |
| // when forming the right-hand side. |
| beta := math.Max(norms[j], norms[j+1]) |
| if beta > vcrit { |
| bi.Dscal(n-ki+1, 1/vmax, b[ki*ldb+iv:], 1) |
| vmax = 1 |
| vcrit = bignum |
| } |
| b[j*ldb+iv] -= bi.Ddot(j-ki-1, t[(ki+1)*ldt+j:], ldt, b[(ki+1)*ldb+iv:], ldb) |
| b[(j+1)*ldb+iv] -= bi.Ddot(j-ki-1, t[(ki+1)*ldt+j+1:], ldt, b[(ki+1)*ldb+iv:], ldb) |
| // Solve |
| // [ T[j,j]-wr T[j,j+1] ]^T * X = scale*[ b1 ] |
| // [ T[j+1,j] T[j+1,j+1]-wr ] [ b2 ] |
| scale, _, _ := impl.Dlaln2(true, 2, 1, smin, 1, t[j*ldt+j:], ldt, |
| 1, 1, b[j*ldb+iv:], ldb, wr, 0, x[:3], 2) |
| // Scale if necessary. |
| if scale != 1 { |
| bi.Dscal(n-ki, scale, b[ki*ldb+iv:], ldb) |
| } |
| b[j*ldb+iv] = x[0] |
| b[(j+1)*ldb+iv] = x[2] |
| vmax = math.Max(vmax, math.Max(math.Abs(b[j*ldb+iv]), math.Abs(b[(j+1)*ldb+iv]))) |
| vcrit = bignum / vmax |
| j += 2 |
| } |
| } |
| // Copy the vector x or Q*x to VL and normalize. |
| switch { |
| case howmny != lapack.AllEVMulQ: |
| // No back-transform: copy x to VL and normalize. |
| bi.Dcopy(n-ki, b[ki*ldb+iv:], ldb, vl[ki*ldvl+is:], ldvl) |
| ii := bi.Idamax(n-ki, vl[ki*ldvl+is:], ldvl) + ki |
| remax := 1 / math.Abs(vl[ii*ldvl+is]) |
| bi.Dscal(n-ki, remax, vl[ki*ldvl+is:], ldvl) |
| for k := 0; k < ki; k++ { |
| vl[k*ldvl+is] = 0 |
| } |
| case nb == 1: |
| // Version 1: back-transform each vector with Gemv, Q*x. |
| if n-ki-1 > 0 { |
| bi.Dgemv(blas.NoTrans, n, n-ki-1, |
| 1, vl[ki+1:], ldvl, b[(ki+1)*ldb+iv:], ldb, |
| b[ki*ldb+iv], vl[ki:], ldvl) |
| } |
| ii := bi.Idamax(n, vl[ki:], ldvl) |
| remax := 1 / math.Abs(vl[ii*ldvl+ki]) |
| bi.Dscal(n, remax, vl[ki:], ldvl) |
| default: |
| // Version 2: back-transform block of vectors with Gemm |
| // zero out above vector. |
| for k := 0; k < ki; k++ { |
| b[k*ldb+iv] = 0 |
| } |
| iscomplex[iv] = ip |
| // Back-transform and normalization is done below. |
| } |
| } else { |
| // Complex left eigenvector. |
| |
| // Initial solve: |
| // [ [ T[ki,ki] T[ki,ki+1] ]^T - (wr - i* wi) ]*X = 0. |
| // [ [ T[ki+1,ki] T[ki+1,ki+1] ] ] |
| if math.Abs(t[ki*ldt+ki+1]) >= math.Abs(t[(ki+1)*ldt+ki]) { |
| b[ki*ldb+iv] = wi / t[ki*ldt+ki+1] |
| b[(ki+1)*ldb+iv+1] = 1 |
| } else { |
| b[ki*ldb+iv] = 1 |
| b[(ki+1)*ldb+iv+1] = -wi / t[(ki+1)*ldt+ki] |
| } |
| b[(ki+1)*ldb+iv] = 0 |
| b[ki*ldb+iv+1] = 0 |
| // Form right-hand side. |
| for k := ki + 2; k < n; k++ { |
| b[k*ldb+iv] = -b[ki*ldb+iv] * t[ki*ldt+k] |
| b[k*ldb+iv+1] = -b[(ki+1)*ldb+iv+1] * t[(ki+1)*ldt+k] |
| } |
| // Solve transposed quasi-triangular system: |
| // [ T[ki+2:n,ki+2:n]^T - (wr-i*wi) ]*X = b1+i*b2 |
| vmax := 1.0 |
| vcrit := bignum |
| for j := ki + 2; j < n; { |
| if j == n-1 || t[(j+1)*ldt+j] == 0 { |
| // 1×1 diagonal block. |
| |
| // Scale if necessary to avoid overflow |
| // when forming the right-hand side elements. |
| if norms[j] > vcrit { |
| rec := 1 / vmax |
| bi.Dscal(n-ki, rec, b[ki*ldb+iv:], ldb) |
| bi.Dscal(n-ki, rec, b[ki*ldb+iv+1:], ldb) |
| vmax = 1 |
| vcrit = bignum |
| } |
| b[j*ldb+iv] -= bi.Ddot(j-ki-2, t[(ki+2)*ldt+j:], ldt, b[(ki+2)*ldb+iv:], ldb) |
| b[j*ldb+iv+1] -= bi.Ddot(j-ki-2, t[(ki+2)*ldt+j:], ldt, b[(ki+2)*ldb+iv+1:], ldb) |
| // Solve [ T[j,j]-(wr-i*wi) ]*(X11+i*X12) = b1+i*b2. |
| scale, _, _ := impl.Dlaln2(false, 1, 2, smin, 1, t[j*ldt+j:], ldt, |
| 1, 1, b[j*ldb+iv:], ldb, wr, -wi, x[:2], 2) |
| // Scale if necessary. |
| if scale != 1 { |
| bi.Dscal(n-ki, scale, b[ki*ldb+iv:], ldb) |
| bi.Dscal(n-ki, scale, b[ki*ldb+iv+1:], ldb) |
| } |
| b[j*ldb+iv] = x[0] |
| b[j*ldb+iv+1] = x[1] |
| vmax = math.Max(vmax, math.Max(math.Abs(b[j*ldb+iv]), math.Abs(b[j*ldb+iv+1]))) |
| vcrit = bignum / vmax |
| j++ |
| } else { |
| // 2×2 diagonal block. |
| |
| // Scale if necessary to avoid overflow |
| // when forming the right-hand side elements. |
| if math.Max(norms[j], norms[j+1]) > vcrit { |
| rec := 1 / vmax |
| bi.Dscal(n-ki, rec, b[ki*ldb+iv:], ldb) |
| bi.Dscal(n-ki, rec, b[ki*ldb+iv+1:], ldb) |
| vmax = 1 |
| vcrit = bignum |
| } |
| b[j*ldb+iv] -= bi.Ddot(j-ki-2, t[(ki+2)*ldt+j:], ldt, b[(ki+2)*ldb+iv:], ldb) |
| b[j*ldb+iv+1] -= bi.Ddot(j-ki-2, t[(ki+2)*ldt+j:], ldt, b[(ki+2)*ldb+iv+1:], ldb) |
| b[(j+1)*ldb+iv] -= bi.Ddot(j-ki-2, t[(ki+2)*ldt+j+1:], ldt, b[(ki+2)*ldb+iv:], ldb) |
| 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) |
| // Solve 2×2 complex linear equation |
| // [ [T[j,j] T[j,j+1] ]^T - (wr-i*wi)*I ]*X = scale*b |
| // [ [T[j+1,j] T[j+1,j+1]] ] |
| scale, _, _ := impl.Dlaln2(true, 2, 2, smin, 1, t[j*ldt+j:], ldt, |
| 1, 1, b[j*ldb+iv:], ldb, wr, -wi, x[:], 2) |
| // Scale if necessary. |
| if scale != 1 { |
| bi.Dscal(n-ki, scale, b[ki*ldb+iv:], ldb) |
| bi.Dscal(n-ki, scale, b[ki*ldb+iv+1:], ldb) |
| } |
| b[j*ldb+iv] = x[0] |
| b[j*ldb+iv+1] = x[1] |
| b[(j+1)*ldb+iv] = x[2] |
| b[(j+1)*ldb+iv+1] = x[3] |
| vmax01 := math.Max(math.Abs(x[0]), math.Abs(x[1])) |
| vmax23 := math.Max(math.Abs(x[2]), math.Abs(x[3])) |
| vmax = math.Max(vmax, math.Max(vmax01, vmax23)) |
| vcrit = bignum / vmax |
| j += 2 |
| } |
| } |
| // Copy the vector x or Q*x to VL and normalize. |
| switch { |
| case howmny != lapack.AllEVMulQ: |
| // No back-transform: copy x to VL and normalize. |
| bi.Dcopy(n-ki, b[ki*ldb+iv:], ldb, vl[ki*ldvl+is:], ldvl) |
| bi.Dcopy(n-ki, b[ki*ldb+iv+1:], ldb, vl[ki*ldvl+is+1:], ldvl) |
| emax := 0.0 |
| for k := ki; k < n; k++ { |
| emax = math.Max(emax, math.Abs(vl[k*ldvl+is])+math.Abs(vl[k*ldvl+is+1])) |
| } |
| remax := 1 / emax |
| bi.Dscal(n-ki, remax, vl[ki*ldvl+is:], ldvl) |
| bi.Dscal(n-ki, remax, vl[ki*ldvl+is+1:], ldvl) |
| for k := 0; k < ki; k++ { |
| vl[k*ldvl+is] = 0 |
| vl[k*ldvl+is+1] = 0 |
| } |
| case nb == 1: |
| // Version 1: back-transform each vector with GEMV, Q*x. |
| if n-ki-2 > 0 { |
| bi.Dgemv(blas.NoTrans, n, n-ki-2, |
| 1, vl[ki+2:], ldvl, b[(ki+2)*ldb+iv:], ldb, |
| b[ki*ldb+iv], vl[ki:], ldvl) |
| bi.Dgemv(blas.NoTrans, n, n-ki-2, |
| 1, vl[ki+2:], ldvl, b[(ki+2)*ldb+iv+1:], ldb, |
| b[(ki+1)*ldb+iv+1], vl[ki+1:], ldvl) |
| } else { |
| bi.Dscal(n, b[ki*ldb+iv], vl[ki:], ldvl) |
| bi.Dscal(n, b[(ki+1)*ldb+iv+1], vl[ki+1:], ldvl) |
| } |
| emax := 0.0 |
| for k := 0; k < n; k++ { |
| emax = math.Max(emax, math.Abs(vl[k*ldvl+ki])+math.Abs(vl[k*ldvl+ki+1])) |
| } |
| remax := 1 / emax |
| bi.Dscal(n, remax, vl[ki:], ldvl) |
| bi.Dscal(n, remax, vl[ki+1:], ldvl) |
| default: |
| // Version 2: back-transform block of vectors with GEMM. |
| // Zero out above vector. |
| // Could go from ki-nv+1 to ki-1. |
| for k := 0; k < ki; k++ { |
| b[k*ldb+iv] = 0 |
| b[k*ldb+iv+1] = 0 |
| } |
| iscomplex[iv] = ip |
| iscomplex[iv+1] = -ip |
| iv++ |
| // Back-transform and normalization is done below. |
| } |
| } |
| if nb > 1 { |
| // Blocked version of back-transform. |
| // For complex case, ki2 includes both vectors ki and ki+1. |
| ki2 := ki |
| if ip != 0 { |
| ki2++ |
| } |
| // Columns [0:iv] of work are valid vectors. When the |
| // number of vectors stored reaches nb-1 or nb, or if |
| // this was last vector, do the Gemm. |
| if iv >= nb-2 || ki2 == n-1 { |
| bi.Dgemm(blas.NoTrans, blas.NoTrans, n, iv+1, n-ki2+iv, |
| 1, vl[ki2-iv:], ldvl, b[(ki2-iv)*ldb:], ldb, |
| 0, b[nb:], ldb) |
| // Normalize vectors. |
| var remax float64 |
| for k := 0; k <= iv; k++ { |
| if iscomplex[k] == 0 { |
| // Real eigenvector. |
| ii := bi.Idamax(n, b[nb+k:], ldb) |
| remax = 1 / math.Abs(b[ii*ldb+nb+k]) |
| } else if iscomplex[k] == 1 { |
| // First eigenvector of conjugate pair. |
| emax := 0.0 |
| for ii := 0; ii < n; ii++ { |
| emax = math.Max(emax, math.Abs(b[ii*ldb+nb+k])+math.Abs(b[ii*ldb+nb+k+1])) |
| } |
| remax = 1 / emax |
| // Second eigenvector of conjugate pair |
| // will reuse this value of remax. |
| } |
| bi.Dscal(n, remax, b[nb+k:], ldb) |
| } |
| impl.Dlacpy(blas.All, n, iv+1, b[nb:], ldb, vl[ki2-iv:], ldvl) |
| iv = 0 |
| } else { |
| iv++ |
| } |
| } |
| is++ |
| if ip != 0 { |
| is++ |
| } |
| } |
| |
| return m |
| } |
