sylvester.jl

function fast_solve_sylvester(A, B, C) # 6 m^3
    n = size(A, 1)
    R,U = schur(complex(A)) # m^3
    S,V = schur(complex(B)) # m^3
    M   = U'*C*V 
    
    Y = zeros(n, n)
    m = zeros(n)
    
    for i = 1:n # 2 m^3
        local m = M[:, i] + Y[:, 1:i-1] * S[1:i-1, i] # m^2
        Y[:, i] = (R - S[i, i] .* I) \ m # m^2
    end
    
    U * Y * V' # 2 m^3
end