hessenberg optimization

hessenberg.jl

using BenchmarkTools

import Base.BLAS: BlasInt, @blasfunc, liblapack

# subroutine dhseqr 	( 	
#       character  	JOB,
# 		character  	COMPZ,
# 		integer  	N,
# 		integer  	ILO,
# 		integer  	IHI,
# 		double precision, dimension( ldh, * )  	H,
# 		integer  	LDH,
# 		double precision, dimension( * )  	WR,
# 		double precision, dimension( * )  	WI,
# 		double precision, dimension( ldz, * )  	Z,
# 		integer  	LDZ,
# 		double precision, dimension( * )  	WORK,
# 		integer  	LWORK,
# 		integer  	INFO 
# 	) 	
function do_stuff(H::AbstractMatrix{Float64})
    job = 'E'
    compz = 'N'
    n = size(H, 1)
    ilo = 1
    ihi = n
    wr = Vector{Float64}(n)
    wi = Vector{Float64}(n)
    work = Vector{Float64}(n)
    lwork = n
    info = Ref{BlasInt}()

    ccall((@blasfunc(dhseqr_), LAPACK.liblapack), Void, (
            Ptr{UInt8}, #job 
            Ptr{UInt8}, #compz
            Ptr{BlasInt}, # N
            Ptr{BlasInt}, # ilo
            Ptr{BlasInt}, # ihi
            Ptr{Float64}, # H
            Ptr{BlasInt}, # ldh
            Ptr{Float64}, # wr
            Ptr{Float64}, # wi
            Ptr{Float64}, # z
            Ptr{BlasInt}, # ldz
            Ptr{Float64}, # work
            Ptr{BlasInt}, # lwork
            Ptr{BlasInt}, # info
        ), 
        &job, 
        &compz,
        &n,
        &ilo,
        &ihi,
        H,
        &n,
        wr,
        wi,
        H,
        &n,
        work,
        &lwork,
        info
    )
    
    return wr .+ wi .* im
end

function bench_stuff(n::Int)
    H = triu(rand(n, n), -1)

    if sort!(eigvals(H), by = abs) ≈ sort!(do_stuff(copy(H)), by = abs)
        println("we're good")
    else
        println("not approximately equal")
    end

    bench_1 = @benchmark eigvals($H)
    bench_2 = @benchmark do_stuff(myH) setup = (myH = copy($H))

    return bench_1, bench_2
end

# julia> a,b = bench_stuff(50)
# we're good
# (Trial(726.615 μs), Trial(55.815 μs))

second_attempt.jl

using BenchmarkTools

import Base.BLAS: BlasInt, @blasfunc, liblapack

# subroutine dhseqr 	( 	
#       character  	JOB,
# 		character  	COMPZ,
# 		integer  	N,
# 		integer  	ILO,
# 		integer  	IHI,
# 		double precision, dimension( ldh, * )  	H,
# 		integer  	LDH,
# 		double precision, dimension( * )  	WR,
# 		double precision, dimension( * )  	WI,
# 		double precision, dimension( ldz, * )  	Z,
# 		integer  	LDZ,
# 		double precision, dimension( * )  	WORK,
# 		integer  	LWORK,
# 		integer  	INFO 
# 	) 	
function do_stuff(H::AbstractMatrix{Float64}, Z::AbstractMatrix{Float64})
    job = 'E'
    compz = 'I'
    n = size(H, 1)
    ilo = 1
    ihi = n
    wr = Vector{Float64}(n)
    wi = Vector{Float64}(n)
    work = Vector{Float64}(n)
    lwork = n
    info = Ref{BlasInt}()

    ccall((@blasfunc(dhseqr_), liblapack), Void, (
            Ptr{UInt8}, #job 
            Ptr{UInt8}, #compz
            Ptr{BlasInt}, # N
            Ptr{BlasInt}, # ilo
            Ptr{BlasInt}, # ihi
            Ptr{Float64}, # H
            Ptr{BlasInt}, # ldh
            Ptr{Float64}, # wr
            Ptr{Float64}, # wi
            Ptr{Float64}, # z
            Ptr{BlasInt}, # ldz
            Ptr{Float64}, # work
            Ptr{BlasInt}, # lwork
            Ptr{BlasInt}, # info
        ), 
        &job, 
        &compz,
        &n,
        &ilo,
        &ihi,
        H,
        &n,
        wr,
        wi,
        Z,
        &n,
        work,
        &lwork,
        info
    )
    
    return wr .+ wi * im, H, Z
end

function bench_stuff(n::Int)
    H = triu(rand(n, n), -1)
    Z = similar(H)

    λ1 = sort!(abs.(eigvals(H)))
    λ2 = sort!(abs.(do_stuff(copy(H), copy(Z))[1]))

    @show norm(λ1 - λ2)

    bench_1 = @benchmark eigvals($H)
    bench_2 = @benchmark do_stuff(myH, myZ) setup = (myH = copy($H); myZ = copy($Z))

    return bench_1, bench_2
end

# julia> bench_stuff(50)
# norm(λ1 - λ2) = 8.771836821698865e-12
# (Trial(689.878 μs), Trial(117.148 μs))