Why this function is type unstable

Chong_Wang · May 29, 2018, 1:57am

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:
      unless (Base.not_int)((#temp#@_8::Int64 === (Base.add_int)(SSAValue(88), 1)::Int64)::Bool)::Bo
ol goto 209
      SSAValue(89) = #temp#@_8::Int64
      SSAValue(90) = (Base.add_int)(#temp#@_8::Int64, 1)::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)))))
      SSAValue(20) = (Base.add_float)(l::Float64, SSAValue(16))::Float64
      # 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))
      SSAValue(2) = $(Expr(:new, Complex{Float64}, :((Base.add_float)((Base.sub_float)((Base.mul_flo
at)((Base.div_float)(#temp#@_12, m)::Float64, (Base.add_float)((Base.sitofp)(Float64, (Base.sub_int)
(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(3) = (Base.add_int)(i::Int64, 100)::Int64
      SSAValue(4) = (Base.add_int)(i::Int64, 100)::Int64
      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)))))
      SSAValue(26) = (Base.add_float)(l::Float64, SSAValue(22))::Float64
      # 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))
      SSAValue(5) = $(Expr(:new, Complex{Float64}, :((Base.add_float)((Base.add_float)((Base.mul_flo
at)((Base.div_float)(#temp#@_13, m)::Float64, (Base.add_float)((Base.sitofp)(Float64, (Base.sub_int)
(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(6) = (Base.sub_int)((Base.add_int)(i::Int64, 100)::Int64, 1)::Int64
      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)))))
      SSAValue(32) = (Base.add_float)(l::Float64, SSAValue(28))::Float64
      # 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
      SSAValue(7) = $(Expr(:new, Complex{Float64}, :((Base.add_float)((Core.getfield)(SSAValue(46),
:re)::Float64, SSAValue(93))::Float64), :((Base.add_float)((Core.getfield)(SSAValue(46), :im)::Float
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(9) = (Base.add_int)(i::Int64, 100)::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)))))
      SSAValue(52) = (Base.add_float)(l::Float64, SSAValue(48))::Float64
      # 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
      SSAValue(10) = $(Expr(:new, Complex{Float64}, :((Base.add_float)((Base.mul_float)((Base.mul_fl
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(11) = (Base.add_int)(i::Int64, 1)::Int64
      SSAValue(12) = (Base.add_int)(i::Int64, 100)::Int64
      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)))))
      SSAValue(60) = (Base.add_float)(l::Float64, SSAValue(56))::Float64
      # 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
      SSAValue(13) = $(Expr(:new, Complex{Float64}, :((Base.add_float)((Core.getfield)(SSAValue(75),
 :re)::Float64, SSAValue(100))::Float64), :((Base.add_float)((Core.getfield)(SSAValue(75), :im)::Flo
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(14) = (Base.add_int)((Base.add_int)(i::Int64, 100)::Int64, 1)::Int64
      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)))))
      SSAValue(81) = (Base.add_float)(l::Float64, SSAValue(77))::Float64
      # 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
      SSAValue(15) = $(Expr(:new, Complex{Float64}, :((Base.add_float)((Base.mul_float)((Base.mul_fl
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?

yuyichao · May 29, 2018, 2:26am

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

Chong_Wang · May 29, 2018, 3:40am

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

mauro3 · May 29, 2018, 7:37am

https://github.com/JuliaLang/julia/pull/24795
and a bunch more referenced therein.

kristoffer.carlsson · May 29, 2018, 8:31am

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
...
   │     %670 = Base.Broadcast.combine_styles(%663)::Union{Base.Broadcast.Style{Tuple}, Base.Broadcast.Unknown}                                                                          │╻            broadcasted
   │     %671 = Base.Broadcast.broadcasted(%670, %669, %663)::Union{Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple},Nothing,typeof(real),_1} where _1, Base.Broadcast.Broadcasted{Base.Broadcast.Unknown,Nothing,typeof(real),_1} where _1}
   │     %672 = %569(%671)::Any                                                                                                                                                          │
   │     %673 = Base.getindex(%672, 1)::Any                                                                                                                                              │
   │     %674 = Main.:+(%673, %%ϵ0)::Any                                                                                                                                                 │
   └────        return %674

Chong_Wang · May 29, 2018, 12:10pm

Really, why?

kristoffer.carlsson · May 29, 2018, 12:18pm

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.

Chong_Wang · May 30, 2018, 3:26am

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.

Elrod · May 30, 2018, 3:38am

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.

Tamas_Papp · May 30, 2018, 5:54am

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.

Chong_Wang · May 30, 2018, 6:15am

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 .

Elrod · May 30, 2018, 6:18am

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).

Tamas_Papp · May 30, 2018, 6:36am

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.