# Type instability in complex schurfact

Hello everyone,

I’m trying to use a complex Schur decomposition and retrieve the generalized eigenvalues. A very much simplified example is here:

``````const A = complex(randn(10,10))
const B = complex(randn(10,10))

function schur_test(A, B)
F = schurfact(A, B)
gen_ev = F[:values]
return gen_ev
end

gen_ev = schur_test(A, B)
``````

However, the type of `F[:values]` can’t be inferred.

``````@code_warntype schur_test(A, B)

Variables:
#self# <optimized out>
A::Array{Complex{Float64},2}
B::Array{Complex{Float64},2}
F::Base.LinAlg.GeneralizedSchur{Complex{Float64},Array{Complex{Float64},2}}
gen_ev::Any

Body:
begin
\$(Expr(:inbounds, false))
# meta: location linalg/schur.jl schurfact 197
SSAValue(1) = \$(Expr(:foreigncall, :(:jl_array_copy), Ref{Array{Complex{Float64},2}}, svec(Any), :(A), 0))
SSAValue(0) = \$(Expr(:foreigncall, :(:jl_array_copy), Ref{Array{Complex{Float64},2}}, svec(Any), :(B), 0))
# meta: location linalg/schur.jl schurfact! 184
SSAValue(2) = \$(Expr(:invoke, MethodInstance for gges!(::Char, ::Char, ::Array{Complex{Float64},2}, ::Array{Complex{Float64},2}), :(\$(QuoteNode(Base.LinAlg.LAPACK.gges!))), 'V', 'V', SSAValue(1), SSAValue(0)))
# meta: pop location
# meta: pop location
\$(Expr(:inbounds, :pop))
F::Base.LinAlg.GeneralizedSchur{Complex{Float64},Array{Complex{Float64},2}} = \$(Expr(:new, Base.LinAlg.GeneralizedSchur{Complex{Float64},Array{Complex{Float64},2}},:((Core.getfield)(SSAValue(2), 1)::Array{Complex{Float64},2}), :((Core.getfield)(SSAValue(2), 2)::Array{Complex{Float64},2}), :((Core.getfield)(SSAValue(2), 3)::Array{Complex{Float64},1}), :((Core.getfield)(SSAValue(2), 4)::Array{Complex{Float64},1}), :((Core.getfield)(SSAValue(2), 5)::Array{Complex{Float64},2}), :((Core.getfield)(SSAValue(2), 6)::Array{Complex{Float64},2}))) # line 145:
gen_ev::Any = \$(Expr(:invoke, MethodInstance for getindex(::Base.LinAlg.GeneralizedSchur{Complex{Float64},Array{Complex{Float64},2}}, ::Symbol), :(Main.getindex), :(F), :(:values))) # line 146:
return gen_ev::Any
end::Any
``````

Putting type information on the function arguments or on `gen_ev = Array{Complex{Float64},1}(F[:values])` doesn’t change the fact that `gen_env` and any variables that I want to construct using it are not type stable.

Any hints are very much appreciated.

Thanks, Benjamin

System:

``````Julia Version 0.6.4
Commit 9d11f62bcb (2018-07-09 19:09 UTC)
Platform Info:
OS: macOS (x86_64-apple-darwin17.7.0)
CPU: Intel(R) Core(TM) i7-6567U CPU @ 3.30GHz
WORD_SIZE: 64
BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell MAX_THREADS=16)
LAPACK: libopenblas64_
LIBM: libopenlibm
LLVM: libLLVM-3.9.1 (ORCJIT, skylake)
``````

Basically, in the GeneralizedSchur object, `alpha` is not a concrete type:

``````struct GeneralizedSchur{Ty,M<:AbstractMatrix} <: Factorization{Ty}
S::M
T::M
alpha::Vector
beta::Vector{Ty}
Q::M
Z::M
function GeneralizedSchur{Ty,M}(S::AbstractMatrix{Ty}, T::AbstractMatrix{Ty}, alpha::Vector,
beta::Vector{Ty}, Q::AbstractMatrix{Ty}, Z::AbstractMatrix{Ty}) where {Ty,M}
new(S, T, alpha, beta, Q, Z)
end
end
``````

and getindex is defined as `F.alpha ./ F.beta`:

``````function getindex(F::GeneralizedSchur, d::Symbol)
if d == :S
return F.S
elseif d == :T
return F.T
elseif d == :alpha
return F.alpha
elseif d == :beta
return F.beta
elseif d == :values
return F.alpha./F.beta
elseif d == :Q || d == :left
return F.Q
elseif d == :Z || d == :right
return F.Z
else
throw(KeyError(d))
end
end
``````

So, you could do something like this:

``````function schur_test(A, B)
F = schurfact(A, B)
alpha::Vector{complex(promote_type(eltype(A),eltype(B)))} = F.alpha
return alpha ./ F.beta
end
``````

You could simplify the `complex(promote_type(eltype(A),eltype(B))`, but I tried to make it fairly generic.

1 Like

Thanks a lot Chris. That solved it.

Just a follow up question to understand this. Why is your solution not type stable if I exchange `F.alpha` and `F.beta` by `F[:alpha]` and `F[:beta]`?

Because in the “getindex” function I showed above, the return type is dependent on the value of the symbol `d`, not the type.
I think similar code should be type stable on Julia 0.7, because of better constant propagation.

However, in 0.7 you’ll be forced to use the `F.beta` syntax instead of `F[:beta]`, because it dropped the `getindex` notation, replacing it with getproperty overloading.

I think the type instability is intentional in the real case, because you don’t know if the eigenvalues are real or not. In the complex case however, it should be possible to make it type stable. Probably worth opening an issue.

2 Likes