Type instability in complex schurfact

General Usage Performance
1

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

https://github.com/JuliaLang/julia/blob/v0.6.4/base/linalg/schur.jl

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
3

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]?

4

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.
https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/schur.jl

5

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