JuliaMath / Polynomials.jl

Polynomial manipulations in Julia
http://juliamath.github.io/Polynomials.jl/
Other
303 stars 74 forks source link

Use LAPack's `zhseqr` / `dhseqr`? #130

Open dlfivefifty opened 7 years ago

dlfivefifty commented 7 years ago

LAPack's Hessenberg eigenvalue routines use half the memory. I'd be willing to make a pull request that moves (and modernizes) ApproxFun.jl's implementation over: https://github.com/JuliaApproximation/ApproxFun.jl/blob/master/src/LinearAlgebra/hesseneigs.jl

Note that there is an unfortunate side effect that you don't get auto-column scaling, so this may have side effects. I got around this with an adhoc column scaling.

jverzani commented 7 years ago

Should this be a special method for eigvals in base? If that isn't appropriate, I'd be interested in the PR here or in PolynomialZeros.

dlfivefifty commented 7 years ago

There isn't really a Hessenberg type in Base: the Hessenberg is in fact a Hessenberg factorization:

julia> a = rand(3,3)
3×3 Array{Float64,2}:
 0.490033  0.522445  0.854002
 0.309698  0.557504  0.495107
 0.980559  0.81369   0.961208

julia> LinAlg.Hessenberg(a)|>full
3×3 Array{Float64,2}:
 0.490033  0.522445  0.854002
 0.309698  0.557504  0.495107
 0.980559  0.81369   0.961208

julia> hessfact(a)[:H] # returns a Matrix{Float64}
3×3 Array{Float64,2}:
 0.490033   0.651441  -0.760193
 1.26273    0.539174  -0.807753
 0.0       -0.48917    0.979538
dlfivefifty commented 7 years ago

Felt motivated to create an issue to avoid your (presumably) current and my previous confusion on what Hessenberg should be.