# Why this function is type unstable

I have a function here

``````function egvalsbyed(m, α, β, l, ϵ0)
nhos = 100
gs = 2.0
H = spzeros(Complex128, 2*nhos, 2*nhos)
for i in 1:nhos
H[i, i] += l^(-2)/m*(i-1+1/2)-gs/(4*m*l^2)
H[i+nhos, i+nhos] += l^(-2)/m*(i-1+1/2)+gs/(4*m*l^2)
if i > 1
H[i+nhos-1, i] += im*√2*α*l^(-1)*√(i-1)
H[i-1, i+nhos] += √2*β*l^(-1)*√(i-1)
end
if i < nhos
H[i+1, i+nhos] += -im*√2*α*l^(-1)*√(i)
H[i+nhos+1, i] += √2*β*l^(-1)*√(i)
end
end
return real.(eigs(H-ϵ0*speye(2*nhos), which=:SM))[1]+ϵ0
end
``````

An example is

``````julia> egvalsbyed(0.01, 0.05, 0.03, 100.0, 0.4)
6-element Array{Float64,1}:
0.404381
0.395597
0.405545
0.394328
0.38565
0.414433
``````

However, if I do @code_warntype

``````julia> @code_warntype(egvalsbyed(0.01, 0.05, 0.03, 100.0, 0.4))
Variables:
#self# <optimized out>
m::Float64
α::Float64
β::Float64
l::Float64
ϵ0::Float64
i::Int64
#temp#@_8::Int64
nhos <optimized out>
gs <optimized out>
H::SparseMatrixCSC{Complex{Float64},Int64}
#temp#@_12::Float64
#temp#@_13::Float64
#temp#@_14::Float64
#temp#@_15::Float64
#temp#@_16::Float64
#temp#@_17::Float64

Body:
begin  # line 99:
H::SparseMatrixCSC{Complex{Float64},Int64} = \$(Expr(:invoke, MethodInstance for spzeros(::Type
{Complex{Float64}}, ::Type{Int64}, ::Int64, ::Int64), :(Base.SparseArrays.spzeros), Complex{Float64}
, :(Base.SparseArrays.Int), :((Base.mul_int)(2, 100)::Int64), :((Base.mul_int)(2, 100)::Int64))) # l
ine 100:
SSAValue(88) = (Base.select_value)((Base.sle_int)(1, 100)::Bool, 100, (Base.sub_int)(1, 1)::In
t64)::Int64
#temp#@_8::Int64 = 1
6:
ol goto 209
SSAValue(89) = #temp#@_8::Int64
i::Int64 = SSAValue(89)
#temp#@_8::Int64 = SSAValue(90) # line 101:
SSAValue(21) = \$(Expr(:invoke, MethodInstance for getindex(::SparseMatrixCSC{Complex{Float64},
Int64}, ::Int64, ::Int64), :(Main.getindex), :(H), :(i), :(i)))
\$(Expr(:inbounds, false))
# meta: location intfuncs.jl literal_pow 208
# meta: location math.jl ^ 701
SSAValue(16) = (Base.sitofp)(Float64, -2)::Float64
# meta: location math.jl ^ 699
SSAValue(19) = \$(Expr(:foreigncall, "llvm.pow.f64", Float64, svec(Float64, Float64), :(l), 0,
SSAValue(16), 0, :(\$(Expr(:llvmcall)))))
# meta: location math.jl nan_dom_err 300
unless (Base.and_int)((Base.ne_float)(SSAValue(19), SSAValue(19))::Bool, (Base.not_int)((Base.
ne_float)(SSAValue(20), SSAValue(20))::Bool)::Bool)::Bool goto 25
#temp#@_12::Float64 = (Base.Math.throw)(\$(QuoteNode(DomainError())))::Union{}
goto 27
25:
#temp#@_12::Float64 = SSAValue(19)
27:
# meta: pop location
# meta: pop location
# meta: pop location
# meta: pop location
\$(Expr(:inbounds, :pop))
(i, 1)::Int64)::Float64, (Base.div_float)((Base.sitofp)(Float64, 1)::Float64, (Base.sitofp)(Float64,
2)::Float64)::Float64)::Float64)::Float64, (Base.div_float)(2.0, (Base.mul_float)((Base.mul_float)(
(Base.sitofp)(Float64, 4)::Float64, m)::Float64, (Base.mul_float)(l, l)::Float64)::Float64)::Float64
)::Float64, (Core.getfield)(SSAValue(21), :re)::Float64)::Float64), :((Core.getfield)(SSAValue(21),
:im)::Float64)))
\$(Expr(:invoke, MethodInstance for setindex!(::SparseMatrixCSC{Complex{Float64},Int64}, ::Comp
lex{Float64}, ::Int64, ::Int64), :(Main.setindex!), :(H), SSAValue(2), :(i), :(i))) # line 102:
SSAValue(27) = \$(Expr(:invoke, MethodInstance for getindex(::SparseMatrixCSC{Complex{Float64},
Int64}, ::Int64, ::Int64), :(Main.getindex), :(H), SSAValue(3), SSAValue(4)))
\$(Expr(:inbounds, false))
# meta: location intfuncs.jl literal_pow 208
# meta: location math.jl ^ 701
SSAValue(22) = (Base.sitofp)(Float64, -2)::Float64
# meta: location math.jl ^ 699
SSAValue(25) = \$(Expr(:foreigncall, "llvm.pow.f64", Float64, svec(Float64, Float64), :(l), 0,
SSAValue(22), 0, :(\$(Expr(:llvmcall)))))
# meta: location math.jl nan_dom_err 300
unless (Base.and_int)((Base.ne_float)(SSAValue(25), SSAValue(25))::Bool, (Base.not_int)((Base.
ne_float)(SSAValue(26), SSAValue(26))::Bool)::Bool)::Bool goto 50
#temp#@_13::Float64 = (Base.Math.throw)(\$(QuoteNode(DomainError())))::Union{}
goto 52
50:
#temp#@_13::Float64 = SSAValue(25)
52:
# meta: pop location
# meta: pop location
# meta: pop location
# meta: pop location
\$(Expr(:inbounds, :pop))
(i, 1)::Int64)::Float64, (Base.div_float)((Base.sitofp)(Float64, 1)::Float64, (Base.sitofp)(Float64,
2)::Float64)::Float64)::Float64)::Float64, (Base.div_float)(2.0, (Base.mul_float)((Base.mul_float)(
(Base.sitofp)(Float64, 4)::Float64, m)::Float64, (Base.mul_float)(l, l)::Float64)::Float64)::Float64
)::Float64, (Core.getfield)(SSAValue(27), :re)::Float64)::Float64), :((Core.getfield)(SSAValue(27),
:im)::Float64)))
\$(Expr(:invoke, MethodInstance for setindex!(::SparseMatrixCSC{Complex{Float64},Int64}, ::Comp
lex{Float64}, ::Int64, ::Int64), :(Main.setindex!), :(H), SSAValue(5), SSAValue(3), SSAValue(4))) #
line 103:
unless (Base.slt_int)(1, i::Int64)::Bool goto 132 # line 104:
SSAValue(46) = \$(Expr(:invoke, MethodInstance for getindex(::SparseMatrixCSC{Complex{Float64},
Int64}, ::Int64, ::Int64), :(Main.getindex), :(H), SSAValue(6), :(i)))
SSAValue(33) = (Base.Math.sqrt_llvm)((Base.sitofp)(Float64, 2)::Float64)::Float64
\$(Expr(:inbounds, false))
# meta: location intfuncs.jl literal_pow 208
# meta: location math.jl ^ 701
SSAValue(28) = (Base.sitofp)(Float64, -1)::Float64
# meta: location math.jl ^ 699
SSAValue(31) = \$(Expr(:foreigncall, "llvm.pow.f64", Float64, svec(Float64, Float64), :(l), 0,
SSAValue(28), 0, :(\$(Expr(:llvmcall)))))
# meta: location math.jl nan_dom_err 300
unless (Base.and_int)((Base.ne_float)(SSAValue(31), SSAValue(31))::Bool, (Base.not_int)((Base.
ne_float)(SSAValue(32), SSAValue(32))::Bool)::Bool)::Bool goto 77
#temp#@_14::Float64 = (Base.Math.throw)(\$(QuoteNode(DomainError())))::Union{}
goto 79
77:
#temp#@_14::Float64 = SSAValue(31)
79:
# meta: pop location
# meta: pop location
# meta: pop location
# meta: pop location
\$(Expr(:inbounds, :pop))
SSAValue(91) = #temp#@_14::Float64
SSAValue(92) = (Base.Math.sqrt_llvm)((Base.sitofp)(Float64, (Base.sub_int)(i::Int64, 1)::Int64
)::Float64)::Float64
\$(Expr(:inbounds, false))
# meta: location operators.jl * 424
SSAValue(42) = (Base.select_value)((Core.getfield)(Main.im, :re)::Bool, SSAValue(33), (Base.co
pysign_float)((Base.sitofp)(Float64, 0)::Float64, SSAValue(33))::Float64)::Float64
SSAValue(43) = (Base.select_value)((Core.getfield)(Main.im, :im)::Bool, SSAValue(33), (Base.co
pysign_float)((Base.sitofp)(Float64, 0)::Float64, SSAValue(33))::Float64)::Float64
SSAValue(44) = (Base.mul_float)(α::Float64, SSAValue(42))::Float64
SSAValue(45) = (Base.mul_float)(α::Float64, SSAValue(43))::Float64
SSAValue(36) = SSAValue(91)
# meta: location operators.jl afoldl 412
SSAValue(40) = (Base.mul_float)(SSAValue(36), SSAValue(44))::Float64
SSAValue(41) = (Base.mul_float)(SSAValue(36), SSAValue(45))::Float64
SSAValue(38) = SSAValue(92)
# meta: pop location
# meta: pop location
\$(Expr(:inbounds, :pop))
SSAValue(93) = (Base.mul_float)(SSAValue(38), SSAValue(40))::Float64
SSAValue(94) = (Base.mul_float)(SSAValue(38), SSAValue(41))::Float64
64, SSAValue(94))::Float64)))
\$(Expr(:invoke, MethodInstance for setindex!(::SparseMatrixCSC{Complex{Float64},Int64}, ::Comp
lex{Float64}, ::Int64, ::Int64), :(Main.setindex!), :(H), SSAValue(7), SSAValue(6), :(i))) # line 10
5:
SSAValue(8) = (Base.sub_int)(i::Int64, 1)::Int64
SSAValue(55) = \$(Expr(:invoke, MethodInstance for getindex(::SparseMatrixCSC{Complex{Float64},
Int64}, ::Int64, ::Int64), :(Main.getindex), :(H), SSAValue(8), SSAValue(9)))
SSAValue(53) = (Base.Math.sqrt_llvm)((Base.sitofp)(Float64, 2)::Float64)::Float64
\$(Expr(:inbounds, false))
# meta: location intfuncs.jl literal_pow 208
# meta: location math.jl ^ 701
SSAValue(48) = (Base.sitofp)(Float64, -1)::Float64
# meta: location math.jl ^ 699
SSAValue(51) = \$(Expr(:foreigncall, "llvm.pow.f64", Float64, svec(Float64, Float64), :(l), 0,
SSAValue(48), 0, :(\$(Expr(:llvmcall)))))
# meta: location math.jl nan_dom_err 300
unless (Base.and_int)((Base.ne_float)(SSAValue(51), SSAValue(51))::Bool, (Base.not_int)((Base.
ne_float)(SSAValue(52), SSAValue(52))::Bool)::Bool)::Bool goto 121
#temp#@_15::Float64 = (Base.Math.throw)(\$(QuoteNode(DomainError())))::Union{}
goto 123
121:
#temp#@_15::Float64 = SSAValue(51)
123:
# meta: pop location
# meta: pop location
# meta: pop location
# meta: pop location
\$(Expr(:inbounds, :pop))
SSAValue(95) = (Base.Math.sqrt_llvm)((Base.sitofp)(Float64, (Base.sub_int)(i::Int64, 1)::Int64
)::Float64)::Float64
oat)((Base.mul_float)(SSAValue(53), β)::Float64, #temp#@_15)::Float64, SSAValue(95))::Float64, (Cor
e.getfield)(SSAValue(55), :re)::Float64)::Float64), :((Core.getfield)(SSAValue(55), :im)::Float64)))
\$(Expr(:invoke, MethodInstance for setindex!(::SparseMatrixCSC{Complex{Float64},Int64}, ::Comp
lex{Float64}, ::Int64, ::Int64), :(Main.setindex!), :(H), SSAValue(10), SSAValue(8), SSAValue(9)))
132:  # line 107:
unless (Base.slt_int)(i::Int64, 100)::Bool goto 207 # line 108:
SSAValue(75) = \$(Expr(:invoke, MethodInstance for getindex(::SparseMatrixCSC{Complex{Float64},
Int64}, ::Int64, ::Int64), :(Main.getindex), :(H), SSAValue(11), SSAValue(12)))
SSAValue(61) = (Base.Math.sqrt_llvm)((Base.sitofp)(Float64, 2)::Float64)::Float64
\$(Expr(:inbounds, false))
# meta: location intfuncs.jl literal_pow 208
# meta: location math.jl ^ 701
SSAValue(56) = (Base.sitofp)(Float64, -1)::Float64
# meta: location math.jl ^ 699
SSAValue(59) = \$(Expr(:foreigncall, "llvm.pow.f64", Float64, svec(Float64, Float64), :(l), 0,
SSAValue(56), 0, :(\$(Expr(:llvmcall)))))
# meta: location math.jl nan_dom_err 300
unless (Base.and_int)((Base.ne_float)(SSAValue(59), SSAValue(59))::Bool, (Base.not_int)((Base.
ne_float)(SSAValue(60), SSAValue(60))::Bool)::Bool)::Bool goto 151
#temp#@_16::Float64 = (Base.Math.throw)(\$(QuoteNode(DomainError())))::Union{}
goto 153
151:
#temp#@_16::Float64 = SSAValue(59)
153:
# meta: pop location
# meta: pop location
# meta: pop location
# meta: pop location
\$(Expr(:inbounds, :pop))
SSAValue(96) = (Base.neg_int)((Base.and_int)((Base.zext_int)(Int64, (Core.getfield)(Main.im, :
re)::Bool)::Int64, 1)::Int64)::Int64
SSAValue(97) = (Base.neg_int)((Base.and_int)((Base.zext_int)(Int64, (Core.getfield)(Main.im, :
im)::Bool)::Int64, 1)::Int64)::Int64
SSAValue(98) = #temp#@_16::Float64
SSAValue(99) = (Base.Math.sqrt_llvm)((Base.sitofp)(Float64, i::Int64)::Float64)::Float64
\$(Expr(:inbounds, false))
# meta: location operators.jl * 424
SSAValue(71) = (Base.mul_float)(SSAValue(61), (Base.sitofp)(Float64, SSAValue(96))::Float64)::
Float64
SSAValue(72) = (Base.mul_float)(SSAValue(61), (Base.sitofp)(Float64, SSAValue(97))::Float64)::
Float64
SSAValue(73) = (Base.mul_float)(α::Float64, SSAValue(71))::Float64
SSAValue(74) = (Base.mul_float)(α::Float64, SSAValue(72))::Float64
SSAValue(65) = SSAValue(98)
# meta: location operators.jl afoldl 412
SSAValue(69) = (Base.mul_float)(SSAValue(65), SSAValue(73))::Float64
SSAValue(70) = (Base.mul_float)(SSAValue(65), SSAValue(74))::Float64
SSAValue(67) = SSAValue(99)
# meta: pop location
# meta: pop location
\$(Expr(:inbounds, :pop))
SSAValue(100) = (Base.mul_float)(SSAValue(67), SSAValue(69))::Float64
SSAValue(101) = (Base.mul_float)(SSAValue(67), SSAValue(70))::Float64
at64, SSAValue(101))::Float64)))
\$(Expr(:invoke, MethodInstance for setindex!(::SparseMatrixCSC{Complex{Float64},Int64}, ::Comp
lex{Float64}, ::Int64, ::Int64), :(Main.setindex!), :(H), SSAValue(13), SSAValue(11), SSAValue(12)))
# line 109:
SSAValue(84) = \$(Expr(:invoke, MethodInstance for getindex(::SparseMatrixCSC{Complex{Float64},
Int64}, ::Int64, ::Int64), :(Main.getindex), :(H), SSAValue(14), :(i)))
SSAValue(82) = (Base.Math.sqrt_llvm)((Base.sitofp)(Float64, 2)::Float64)::Float64
\$(Expr(:inbounds, false))
# meta: location intfuncs.jl literal_pow 208
# meta: location math.jl ^ 701
SSAValue(77) = (Base.sitofp)(Float64, -1)::Float64
# meta: location math.jl ^ 699
SSAValue(80) = \$(Expr(:foreigncall, "llvm.pow.f64", Float64, svec(Float64, Float64), :(l), 0,
SSAValue(77), 0, :(\$(Expr(:llvmcall)))))
# meta: location math.jl nan_dom_err 300
unless (Base.and_int)((Base.ne_float)(SSAValue(80), SSAValue(80))::Bool, (Base.not_int)((Base.
ne_float)(SSAValue(81), SSAValue(81))::Bool)::Bool)::Bool goto 196
#temp#@_17::Float64 = (Base.Math.throw)(\$(QuoteNode(DomainError())))::Union{}
goto 198
196:
#temp#@_17::Float64 = SSAValue(80)
198:
# meta: pop location
# meta: pop location
# meta: pop location
# meta: pop location
\$(Expr(:inbounds, :pop))
SSAValue(102) = (Base.Math.sqrt_llvm)((Base.sitofp)(Float64, i::Int64)::Float64)::Float64
oat)((Base.mul_float)(SSAValue(82), β)::Float64, #temp#@_17)::Float64, SSAValue(102))::Float64, (Co
re.getfield)(SSAValue(84), :re)::Float64)::Float64), :((Core.getfield)(SSAValue(84), :im)::Float64))
)
\$(Expr(:invoke, MethodInstance for setindex!(::SparseMatrixCSC{Complex{Float64},Int64}, ::Comp
lex{Float64}, ::Int64, ::Int64), :(Main.setindex!), :(H), SSAValue(15), SSAValue(14), :(i)))
207:
goto 6
209:  # line 112:
SSAValue(87) = \$(Expr(:invoke, MethodInstance for vector_any(::Any, ::Vararg{Any,N} where N),
:(Base.vector_any), :(:which), :(:SM)))
SSAValue(85) = (Base.mul_int)(2, 100)::Int64
SSAValue(86) = \$(Expr(:invoke, MethodInstance for *(::Float64, ::SparseMatrixCSC{Float64,Int64
}), :(Main.*), :(ϵ0), :(\$(Expr(:invoke, MethodInstance for speye_scaled(::Type{Float64}, ::Float64,
::Int64, ::Int64), :(Base.SparseArrays.speye_scaled), Float64, :((Base.sitofp)(Float64, 1)::Float64)
, SSAValue(85), SSAValue(85))))))
return ((Main.getindex)((Base.broadcast)(Main.real, \$(Expr(:invoke, MethodInstance for (::Base
.LinAlg.#kw##eigs)(::Array{Any,1}, ::Base.LinAlg.#eigs, ::SparseMatrixCSC{Complex{Float64},Int64}),
:(\$(QuoteNode(Base.LinAlg.#eigs))), SSAValue(87), :(Main.eigs), :(\$(Expr(:invoke, MethodInstance for
map(::Base.#-, ::SparseMatrixCSC{Complex{Float64},Int64}, ::SparseMatrixCSC{Float64,Int64}), :(Base
.SparseArrays.map), :(Base.SparseArrays.-), :(H), SSAValue(86)))))))::Any, 1)::Any + ϵ0::Float64)::A
ny
end::Any
``````

Usually, I can guess the meaning of the output of `@code_warntype`, but this time I have no clues. What is happening here?

It’s from the kwarg. Should be fixed on 0.7.

Thank you. Can you point me to the specific issue or pull request?

and a bunch more referenced therein.

1 Like

Looks like `eigs` is inherently type unstable though?

On 0.7:

``````julia> @code_warntype(egvalsbyed(0.01, 0.05, 0.03, 100.0, 0.4))
Body::Any
...
│     %663 = invoke IterativeEigensolvers.:(#_eigs#9)(%657::Int64, %658::Int64, %624::Symbol, %659::Float64, %660::Int64, %661::Nothing, %644::Array{Complex{Float64},1}, %662::Bool, %627::typeof(IterativeEigensolvers._eigs), %621::SparseArrays.SparseMatrixCSC{Complex{Float64},Int64}, %628::LinearAlgebra.UniformScaling{Bool})::Tuple
...
│     %672 = %569(%671)::Any                                                                                                                                                          │
│     %673 = Base.getindex(%672, 1)::Any                                                                                                                                              │
│     %674 = Main.:+(%673, %%ϵ0)::Any                                                                                                                                                 │
└────        return %674
``````
1 Like

Really, why?

I think it is because you do not know if the resulting vectors will be Real or Complex. If you know this from your problem you can type assert in your code.

Therefore, `eig` is also type unstable.

``````@code_warntype eig(eye(2))
julia> Variables:
#self# <optimized out>
A::Array{Float64,2}
F::Any

Body:
begin
\$(Expr(:inbounds, false))
# meta: location linalg\eigen.jl #eig#40 108
F::Any = \$(Expr(:invoke, MethodInstance for (::Base.LinAlg.#kw##eigfact)(::Array{Any,1}, ::Bas
e.LinAlg.#eigfact, ::Array{Float64,2}), :(\$(QuoteNode(Base.LinAlg.#eigfact))), :(\$(Expr(:invoke, Met
hodInstance for vector_any(::Any, ::Vararg{Any,N} where N), :(Base.vector_any), :(:permute), true, :
(:scale), true))), :(Base.LinAlg.eigfact), :(A)))
# meta: pop location
\$(Expr(:inbounds, :pop))
return (Core.tuple)((Core.getfield)(F::Any, :values)::Union{Array{Complex{Float64},1}, Array{F
loat64,1}}, (Core.getfield)(F::Any, :vectors)::Union{Array{Complex{Float64},2}, Array{Float64,2}})::
Tuple{Union{Array{Complex{Float64},1}, Array{Float64,1}},Union{Array{Complex{Float64},2}, Array{Floa
t64,2}}}
end::Tuple{Union{Array{Complex{Float64},1}, Array{Float64,1}},Union{Array{Complex{Float64},2}, Arr
ay{Float64,2}}}
``````

This is really surprising, especially because it’s a built-in function and it’s used so widely.

For what it’s worth,

``````S = Symmetric(eye(2));
@code_warntype eig(S)
H = Hermitian(eye(2));
@code_warntype eig(H)
``````

Unfortunately, while

``````julia> function egvalsbyed(m, α, β, l, ϵ0)
nhos = 100
gs = 2.0
H = spzeros(Complex128, 2*nhos, 2*nhos)
for i in 1:nhos
H[i, i] += l^(-2)/m*(i-1+1/2)-gs/(4*m*l^2)
H[i+nhos, i+nhos] += l^(-2)/m*(i-1+1/2)+gs/(4*m*l^2)
if i > 1
H[i+nhos-1, i] += im*√2*α*l^(-1)*√(i-1)
H[i-1, i+nhos] += √2*β*l^(-1)*√(i-1)
end
if i < nhos
H[i+1, i+nhos] += -im*√2*α*l^(-1)*√(i)
H[i+nhos+1, i] += √2*β*l^(-1)*√(i)
end
end
H
end
egvalsbyed (generic function with 1 method)

julia> H = egvalsbyed(0.01, 0.05, 0.03, 100.0, 0.4);

julia> ishermitian(H)
true
``````

eig for sparse Hermitian matrices is still unstable:

``````julia> @code_warntype eig(Hermitian(H))
Variables:
#self# <optimized out>
A::Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Int64}}
args <optimized out>
F::Any

Body:
begin
F::Any = \$(Expr(:invoke, MethodInstance for eigfact(::Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Int64}}), :(Base.LinAlg.eigfact), :(A))) # line 137:
return (Core.tuple)((Core.getfield)(F::Any, :values)::Any, (Core.getfield)(F::Any, :vectors)::Any)::Tuple{Any,Any}
end::Tuple{Any,Any}

julia> @code_warntype eigs(Hermitian(H))
Variables:
#self# <optimized out>
A::Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Int64}}
#temp#::Tuple

Body:
begin
SSAValue(1) = \$(Expr(:foreigncall, :(:jl_alloc_array_1d), Array{Any,1}, svec(Any, Int64), Array{Any,1}, 0, 0, 0))
\$(Expr(:inbounds, false))
# meta: location linalg/arnoldi.jl #eigs#99 90
# meta: location abstractarray.jl isempty 831
# meta: location abstractarray.jl _length 132
# meta: location abstractarray.jl indices 64
SSAValue(4) = (Base.arraysize)(SSAValue(1), 1)::Int64
# meta: pop location
# meta: pop location
# meta: pop location
unless ((Base.select_value)((Base.slt_int)(SSAValue(4), 0)::Bool, 0, SSAValue(4))::Int64 === 0)::Bool goto 14
#temp#::Tuple = \$(Expr(:invoke, MethodInstance for eigs(::Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Int64}}, ::UniformScaling{Int64}), :(Base.LinAlg.eigs), :(A), :(Base.LinAlg.I)))
goto 16
14:
#temp#::Tuple = \$(Expr(:invoke, MethodInstance for (::Base.LinAlg.#kw##eigs)(::Array{Any,1}, ::Base.LinAlg.#eigs, ::Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Int64}}, ::UniformScaling{Int64}), :(\$(QuoteNode(Base.LinAlg.#eigs))), :(\$(Expr(:invoke, MethodInstance for as_kwargs(::Array{Any,1}), :(Base.as_kwargs), SSAValue(1)))), :(Base.LinAlg.eigs), :(A), :(Base.LinAlg.I)))
16:
# meta: pop location
\$(Expr(:inbounds, :pop))
return #temp#::Tuple
end::Tuple

julia> @code_warntype eig(Hermitian(Matrix(H)))
Variables:
#self# <optimized out>
A::Hermitian{Complex{Float64},Array{Complex{Float64},2}}
args <optimized out>
F::Base.LinAlg.Eigen{Complex{Float64},Float64,Array{Complex{Float64},2},Array{Float64,1}}

Body:
begin
F::Base.LinAlg.Eigen{Complex{Float64},Float64,Array{Complex{Float64},2},Array{Float64,1}} = \$(Expr(:invoke, MethodInstance for eigfact(::Hermitian{Complex{Float64},Array{Complex{Float64},2}}), :(Base.LinAlg.eigfact), :(A))) # line 137:
return (Core.tuple)((Core.getfield)(F::Base.LinAlg.Eigen{Complex{Float64},Float64,Array{Complex{Float64},2},Array{Float64,1}}, :values)::Array{Float64,1}, (Core.getfield)(F::Base.LinAlg.Eigen{Complex{Float64},Float64,Array{Complex{Float64},2},Array{Float64,1}}, :vectors)::Array{Complex{Float64},2})::Tuple{Array{Float64,1},Array{Complex{Float64},2}}
end::Tuple{Array{Float64,1},Array{Complex{Float64},2}}

julia> @code_warntype eigs(Hermitian(Matrix(H)))
Variables:
#self# <optimized out>
A::Hermitian{Complex{Float64},Array{Complex{Float64},2}}
#temp#::Tuple

Body:
begin
SSAValue(1) = \$(Expr(:foreigncall, :(:jl_alloc_array_1d), Array{Any,1}, svec(Any, Int64), Array{Any,1}, 0, 0, 0))
\$(Expr(:inbounds, false))
# meta: location linalg/arnoldi.jl #eigs#99 90
# meta: location abstractarray.jl isempty 831
# meta: location abstractarray.jl _length 132
# meta: location abstractarray.jl indices 64
SSAValue(4) = (Base.arraysize)(SSAValue(1), 1)::Int64
# meta: pop location
# meta: pop location
# meta: pop location
unless ((Base.select_value)((Base.slt_int)(SSAValue(4), 0)::Bool, 0, SSAValue(4))::Int64 === 0)::Bool goto 14
#temp#::Tuple = \$(Expr(:invoke, MethodInstance for eigs(::Hermitian{Complex{Float64},Array{Complex{Float64},2}}, ::UniformScaling{Int64}), :(Base.LinAlg.eigs), :(A), :(Base.LinAlg.I)))
goto 16
14:
#temp#::Tuple = \$(Expr(:invoke, MethodInstance for (::Base.LinAlg.#kw##eigs)(::Array{Any,1}, ::Base.LinAlg.#eigs, ::Hermitian{Complex{Float64},Array{Complex{Float64},2}}, ::UniformScaling{Int64}), :(\$(QuoteNode(Base.LinAlg.#eigs))), :(\$(Expr(:invoke, MethodInstance for as_kwargs(::Array{Any,1}), :(Base.as_kwargs), SSAValue(1)))), :(Base.LinAlg.eigs), :(A), :(Base.LinAlg.I)))
16:
# meta: pop location
\$(Expr(:inbounds, :pop))
return #temp#::Tuple
end::Tuple
``````

eig for dense Hermitians is the exception.

It should not be, many functions are not type stable, eg `\`. This is natural, as the type system doesn’t (and possibly can’t, in a sensible design) contain all the information necessary for the branches in the algorithm.

Type stability is an important concept to get performance in some contexts, not an end in itself. Especially since small type unions are handled well in `v0.7`.

1 Like

I understand the design is difficult but type instability can grow rather badly in a complex system (Anything calling this function becomes type unstable). There are methods to avoid this problem, at least for `eig`. For example, always return complex eigenvalues and complex eigenvectors even for Hermitian matrices. It’s up to the programmer to decide whether to convert the eigenvalues to real arrays.

As for the comment “small type unions are handled well in v0.7”, can you point me to the issue or pull request? I always want to look at the situation myself .

Unfortunately, we don’t have a small union yet. On Julia 0.7 (`H` is the matrix from egvalsbyed):

``````julia> @code_warntype eigen(H)
Body::Union{}
1012 1 ─ %1 = new(Core.ErrorException, "Use IterativeEigensolvers.eigs() instead of eigen() for sparse matrices.")::ErrorException                        │╻╷ error
│        Base.throw(%1)                                                                                                                              ││
└──      unreachable                                                                                                                                 ││
2 ─      unreachable

julia> using IterativeEigensolvers

julia> @code_warntype IterativeEigensolvers.eigs(H)
Body::Tuple
47 1 ─ %1  = IterativeEigensolvers._eigs::typeof(IterativeEigensolvers._eigs)                                                                   │╻╷╷      #eigs#1
│   %2  = IterativeEigensolvers.I::UniformScaling{Bool}                                                                                      ││┃│       eigs
│   %3  = Base.ifelse(true, 20, 13)::Int64                                                                                                   │││┃│╷      #eigs#2
│   %4  = IterativeEigensolvers.nothing::Nothing                                                                                             ││││┃        _eigs
│   %5  = \$(Expr(:foreigncall, :(:jl_alloc_array_1d), Array{Complex{Float64},1}, svec(Any, Int64), :(:ccall), 2, Array{Complex{Float64},1}, 0, 0))::Array{Complex{Float64},1}
│   %6  = new(Complex{Float64}, 0.0, 0.0)::Complex{Float64}                                                                                  ││││││╻╷╷      zero
│   %7  = invoke Base.fill!(%5::Array{Complex{Float64},1}, %6::Complex{Float64})::Array{Complex{Float64},1}                                  ││││││
│   %8  = π (6, Int64)                                                                                                                       │││││
│   %9  = π (%3, Int64)                                                                                                                      │││││
│   %10 = π (:(:LM), Symbol)                                                                                                                 │││││
│   %11 = π (0.0, Float64)                                                                                                                   │││││
│   %12 = π (300, Int64)                                                                                                                     │││││
│   %13 = π (true, Bool)                                                                                                                     │││││
│   %14 = invoke IterativeEigensolvers.:(#_eigs#9)(%8::Int64, %9::Int64, %10::Symbol, %11::Float64, %12::Int64, %4::Nothing, %7::Array{Complex{Float64},1}, %13::Bool, %1::typeof(IterativeEigensolvers._eigs), %%A::SparseArrays.SparseMatrixCSC{Complex{Float64},Int64}, %2::UniformScaling{Bool})::Tuple
└──       return %14
``````

However, I wouldn’t be surprised if this changes soon.

Also, regarding “except for Hermitian matrices”, we do already have that eigs for a real `Symmetric` matrix (eg, matrix of type `Symmetric{Float64,Array{Float64,2}}`) is type stable (always real).

Using function barriers, I found that I can contain it when necessary, leaving type unstable code mostly in the high-level logic where it has no significant performance impact.

Eg here, but the whole topic is informative.

2 Likes