diff --git a/Project.toml b/Project.toml
index e873ffe1..685013ea 100644
--- a/Project.toml
+++ b/Project.toml
@@ -86,7 +86,7 @@ IterativeSolvers = "0.9.3"
 JET = "0.8.28, 0.9"
 KLU = "0.6"
 KernelAbstractions = "0.9.27"
-Krylov = "0.9"
+Krylov = "0.9.8"
 KrylovKit = "0.8"
 KrylovPreconditioners = "0.3"
 LazyArrays = "1.8, 2"
diff --git a/src/iterative_wrappers.jl b/src/iterative_wrappers.jl
index f487463c..1d60075a 100644
--- a/src/iterative_wrappers.jl
+++ b/src/iterative_wrappers.jl
@@ -63,6 +63,17 @@ function KrylovJL_GMRES(args...; kwargs...)
     KrylovJL(args...; KrylovAlg = Krylov.gmres!, kwargs...)
 end
 
+"""
+```julia
+KrylovJL_BLOCK_GMRES(args...; gmres_restart = 0, window = 0, kwargs...)
+```
+
+A generic BLOCK-GMRES implementation for square non-Hermitian linear systems with multiple right-hand sides
+"""
+function KrylovJL_BLOCK_GMRES(args...; kwargs...)
+    KrylovJL(args...; KrylovAlg = Krylov.block_gmres!, kwargs...)
+end
+
 """
 ```julia
 KrylovJL_BICGSTAB(args...; kwargs...)
@@ -143,8 +154,10 @@ function get_KrylovJL_solver(KrylovAlg)
         Krylov.CrmrSolver
     elseif (KrylovAlg === Krylov.cg!)
         Krylov.CgSolver
-    elseif (KrylovAlg === Krylov.cg_lanczos!)
+    elseif (KrylovAlg === Krylov.cg_lanczos_shift!)
         Krylov.CgLanczosShiftSolver
+    elseif (KrylovAlg === Krylov.cgls_lanczos_shift!)
+        Krylov.CglsLanczosShiftSolver
     elseif (KrylovAlg === Krylov.cgls!)
         Krylov.CglsSolver
     elseif (KrylovAlg === Krylov.cg_lanczos!)
@@ -163,6 +176,8 @@ function get_KrylovJL_solver(KrylovAlg)
         Krylov.GpmrSolver
     elseif (KrylovAlg === Krylov.fom!)
         Krylov.FomSolver
+    elseif (KrylovAlg === Krylov.block_gmres!)
+        Krylov.BlockGmresSolver
     else
         error("Invalid Krylov method detected")
     end
@@ -181,7 +196,8 @@ function init_cacheval(alg::KrylovJL, A, b, u, Pl, Pr, maxiters::Int, abstol, re
                      alg.KrylovAlg === Krylov.gmres! ||
                      alg.KrylovAlg === Krylov.fgmres! ||
                      alg.KrylovAlg === Krylov.gpmr! ||
-                     alg.KrylovAlg === Krylov.fom!)
+                     alg.KrylovAlg === Krylov.fom! ||
+                     alg.KrylovAlg === Krylov.block_gmres!)
             if A isa SparseMatrixCSC
                 KS(SparseMatrixCSC(0, 0, [1], Int[], eltype(A)[]), eltype(b)[], 1)
             elseif A isa Matrix
@@ -206,7 +222,8 @@ function init_cacheval(alg::KrylovJL, A, b, u, Pl, Pr, maxiters::Int, abstol, re
                      alg.KrylovAlg === Krylov.gmres! ||
                      alg.KrylovAlg === Krylov.fgmres! ||
                      alg.KrylovAlg === Krylov.gpmr! ||
-                     alg.KrylovAlg === Krylov.fom!)
+                     alg.KrylovAlg === Krylov.fom! ||
+                     alg.KrylovAlg === Krylov.block_gmres!)
             KS(A, b, memory)
         elseif (alg.KrylovAlg === Krylov.minres! ||
                 alg.KrylovAlg === Krylov.symmlq! ||
diff --git a/src/simplegmres.jl b/src/simplegmres.jl
index c2596efc..a6bd27a9 100644
--- a/src/simplegmres.jl
+++ b/src/simplegmres.jl
@@ -88,38 +88,7 @@ function update_cacheval!(cache::LinearCache, cacheval::SimpleGMRESCache, name::
     return cacheval
 end
 
-"""
-    (c, s, ρ) = _sym_givens(a, b)
-
-Numerically stable symmetric Givens reflection.
-Given `a` and `b` reals, return `(c, s, ρ)` such that
-
-    [ c  s ] [ a ] = [ ρ ]
-    [ s -c ] [ b ] = [ 0 ].
-"""
-function _sym_givens(a::T, b::T) where {T <: AbstractFloat}
-    # This has taken from Krylov.jl
-    if b == 0
-        c = ifelse(a == 0, one(T), sign(a)) # In Julia, sign(0) = 0.
-        s = zero(T)
-        ρ = abs(a)
-    elseif a == 0
-        c = zero(T)
-        s = sign(b)
-        ρ = abs(b)
-    elseif abs(b) > abs(a)
-        t = a / b
-        s = sign(b) / sqrt(one(T) + t * t)
-        c = s * t
-        ρ = b / s  # Computationally better than ρ = a / c since |c| ≤ |s|.
-    else
-        t = b / a
-        c = sign(a) / sqrt(one(T) + t * t)
-        s = c * t
-        ρ = a / c  # Computationally better than ρ = b / s since |s| ≤ |c|
-    end
-    return (c, s, ρ)
-end
+_sym_givens(a::T, b::T) where {T <: AbstractFloat} = Krylov.sym_givens(a, b)
 
 function _sym_givens!(c, s, R, nr::Int, inner_iter::Int, bsize::Int, Hbis)
     if __is_extension_loaded(Val(:KernelAbstractions))