Eigsolve and extended precision

[zmoitier]

I am trying to compute eigenvalues of matrices with extended precision but the function eigsolve return an error. For example with the following matrix

using LinearAlgebra, LinearMaps, Quadmath, KrylovKit

T = Complex{Float128}
N = 32
A = Tridiagonal(rand(T, N - 1), rand(T, N), rand(T, N - 1))

the matrix version give

julia> vals, vecs, info = eigsolve(A, 3, :LM)
ERROR: LoadError: MethodError: no method matching hschur!(::SubArray{Complex{Float128}, 2, Matrix{Complex{Float128}}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}, ::SubArray{Complex{Float128}, 2, Matrix{Complex{Float128}}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false})

and the matrix-free version give the same error

julia> A_map = LinearMap{T}(x -> A * x, N);
julia> vals, vecs, info = eigsolve(A_map, N, 3, :LM)
ERROR: MethodError: no method matching hschur!(::SubArray{Complex{Float128}, 2, Matrix{Complex{Float128}}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}, ::SubArray{Complex{Float128}, 2, Matrix{Complex{Float128}}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false})

From this discussion on julialang, I have been advise to open this issue.

---

Additional information:

• I get similar errors if I replace Float128 by Double64 from the DoubleFloats.jl package or Float64x2 from MultiFloats.jl package.
• Adding using GenericLinearAlgebra, GenericSchur still yield the error.
• I use julia Version 1.6.2 and KrylovKit v0.5.3.

[Jutho]

Yes, KrylovKit currently uses the Schur decomposition algorithm from BLAS, and thus only works with BLAS numbers. With GenericSchur.jl that could be circumvented, and I welcome PR requests to do that, but I am afraid I won't have time for this myself in the near future.

[zmoitier]

Ha ok, I note for the PR but I also do not have the time right now to immerse myself in how that works.