Generic methods for BiCGStab and GMRES

[spacedome]

Currently the implementations of BiCGStab and GMRES use methods in LinearAlgebra that do not have generic fallbacks.

BiCGStab makes an explicit call to BLAS.gemv! which is implemented only for real and complex floats. I know that the LinearAlgebra package tries to have generic fallback methods, so I am not sure if this is an issue of IterativeSolvers, but currently the code does not work unless I implement gemv! for whatever type I'm using.

GMRES also calls BLAS.gemv! but even if this is fixed, when computing the FastHessenberg the LinearAlgebra function givensAlgorithm is called, which is only defined for real and complex. Fixing the gemv! means GMRES would work with BigFloats and other real and complex types, but I am not sure how generic Givens rotations could be implimented.

Here is an example using Quaternions.jl

using IterativeSolvers
using Quaternions
import Random: AbstractRNG, SamplerType

# rand not implemented in Quaternions, needed to run, but ignore this 
Base.rand(r::AbstractRNG, ::SamplerType{Quaternion{T}}) where {T<:Real} =
    Quaternion(rand(r,T), rand(r,T), rand(r,T), rand(r,T))

q = quat(1.0, 1.0, 1.0, 1.0)
A = [q q'; q' q]
b = [quat(1.0), quat(1.0)]
bicgstabl(A, b, 2)

this gives the following error:

MethodError: no method matching gemv!(::Char, ::Quaternion{Float64}, ::SubArray{Quaternion{Float64},2,Array{Quaternion{Float64},2},Tuple{Base.Slice{Base.OneTo{Int64}},UnitRange{Int64}},true}, ::Array{Quaternion{Float64},1}, ::Quaternion{Float64}, ::SubArray{Quaternion{Float64},1,Array{Quaternion{Float64},2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true})

The same error occurs with BigFloats.

Should LinearAlgebra have a generic fallback for gemv!? Should this package have a fallback for calls to gemv! based on type? I'm not entirely sure of the best way to do the latter, but it shouldn't be difficult, I can try to make a pull request if wanted. For GMRES, I'm not aware of a way to generalize Givens rotations, so maybe it should be explicit about only accepting real and complex? Not sure how important generic implimentations are to other people using this package.

[chriscoey]

anyone have thoughts on this? I need to use BigFloat with GMRES

[mschauer]

IterativeSolvers.idrs works, which is also a Krylov subspace method for solving large nonsymmetric systems of linear equations (http://ta.twi.tudelft.nl/nw/users/gijzen/IDR.html).

[spacedome]

@chriscoey As a workaround you can impliment gemv! for BigFloat with something like this:

function BLAS.gemv!(tA, alpha, A::AbstractMatrix{BigFloat}, x, beta, y)
    y .= alpha .* (if tA == 'T' transpose(A) elseif tA == 'C' A' else A end) * x .+ beta .* y
end

After looking through more of the open issues for LinearAlgebra, the generic equivalent of gemv! is supposed to take the form of mul! or muladd! or something like that, but the name and API hasn't been finalized, see https://github.com/JuliaLang/julia/issues/23919 for details. It might make sense to wait on this issue until LinearAlgebra fixes this?

[dkarrasch]

That's right, and there is also #206 here. The "drawback" of the workaround above is that (if tA == 'T' transpose(A) elseif tA == 'C' A' else A end) * x allocates a vector at each call.