Use `block_gmres` for batched linear solves

[thorek1]

The latest Krylov.jl version has a block-gmres version which I think can simplify/speed-up some parts of the code.

[gdalle]

Thanks for the idea. I don't have much time to work on this at the moment though

[mohamed82008]

how can block-gmres help ID.jl?

[thorek1]

for example here:

    dy = ntuple(Val(N)) do k
        dₖx = partials.(x_and_dx, k)
        dₖx_vec = vec(dₖx)
        dₖc_vec = similar(y_vec)
        mul!(dₖc_vec, B_vec, dₖx_vec)
        dₖy_vec = solve(implicit.linear_solver, A_vec_presolved, -dₖc_vec)
        reshape(dₖy_vec, size(y))
    end

    y_and_dy = map(eachindex(IndexCartesian(), y)) do i
        Dual{T}(y[i], Partials(ntuple(k -> dy[k][i], Val(N))))
    end

instead of solving with gmres for each slice you could use block-gmres once. I didn't check the block-gmres implementation but I think it is worth testing whether this speeds up computation and this might also allow to solve #118

[gdalle]

This is implemented in #135 but it still requires https://github.com/JuliaSmoothOptimizers/Krylov.jl/issues/854 to be fixed