diff --git a/src/LinearSolvers/Krylov/CGSolvers.jl b/src/LinearSolvers/Krylov/CGSolvers.jl
new file mode 100644
index 00000000..69767ddd
--- /dev/null
+++ b/src/LinearSolvers/Krylov/CGSolvers.jl
@@ -0,0 +1,92 @@
+
+struct CGSolver <: Gridap.Algebra.LinearSolver
+  Pl       :: Gridap.Algebra.LinearSolver
+  maxiter  :: Int64
+  atol     :: Float64
+  rtol     :: Float64
+  flexible :: Bool
+  verbose  :: Bool
+end
+
+function CGSolver(Pl;maxiter=10000,atol=1e-12,rtol=1.e-6,flexible=false,verbose=false)
+  return CGSolver(Pl,maxiter,atol,rtol,flexible,verbose)
+end
+
+struct CGSymbolicSetup <: Gridap.Algebra.SymbolicSetup
+  solver
+end
+
+function Gridap.Algebra.symbolic_setup(solver::CGSolver, A::AbstractMatrix)
+  return CGSymbolicSetup(solver)
+end
+
+mutable struct CGNumericalSetup <: Gridap.Algebra.NumericalSetup
+  solver
+  A
+  Pl_ns
+  caches
+end
+
+function get_solver_caches(solver::CGSolver,A)
+  w = allocate_col_vector(A)
+  p = allocate_col_vector(A)
+  z = allocate_col_vector(A)
+  r = allocate_col_vector(A)
+  return (w,p,z,r)
+end
+
+function Gridap.Algebra.numerical_setup(ss::CGSymbolicSetup, A::AbstractMatrix)
+  solver = ss.solver
+  Pl_ns  = numerical_setup(symbolic_setup(solver.Pl,A),A)
+  caches = get_solver_caches(solver,A)
+  return CGNumericalSetup(solver,A,Pl_ns,caches)
+end
+
+function Gridap.Algebra.numerical_setup!(ns::CGNumericalSetup, A::AbstractMatrix)
+  numerical_setup!(ns.Pl_ns,A)
+  ns.A = A
+end
+
+function Gridap.Algebra.solve!(x::AbstractVector,ns::CGNumericalSetup,b::AbstractVector)
+  solver, A, Pl, caches = ns.solver, ns.A, ns.Pl_ns, ns.caches
+  maxiter, atol, rtol = solver.maxiter, solver.atol, solver.rtol
+  flexible, verbose = solver.flexible, solver.verbose
+  w,p,z,r = caches
+  verbose && println(" > Starting CG solver: ")
+
+  # Initial residual
+  mul!(w,A,x); r .= b .- w
+  fill!(p,0.0); γ = 1.0
+
+  res  = norm(r); res_0 = res
+  iter = 0; converged = false
+  while !converged && (iter < maxiter)
+    verbose && println("   > Iteration ", iter," - Residual: ", res)
+
+    if !flexible # β = (zₖ₊₁ ⋅ rₖ₊₁)/(zₖ ⋅ rₖ)
+      solve!(z, Pl, r)
+      β = γ; γ = dot(z, r); β = γ / β
+    else         # β = (zₖ₊₁ ⋅ (rₖ₊₁-rₖ))/(zₖ ⋅ rₖ)
+      β = γ; γ = dot(z, r)
+      solve!(z, Pl, r)
+      γ = dot(z, r) - γ; β = γ / β
+    end
+    p .= z .+ β .* p
+
+    # w = A⋅p
+    mul!(w,A,p)
+    α = γ / dot(p, w)
+
+    # Update solution and residual
+    x .+= α .* p
+    r .-= α .* w
+
+    res = norm(r)
+    converged = (res < atol || res < rtol*res_0)
+    iter += 1
+  end
+  verbose && println("   > Num Iter: ", iter," - Final residual: ", res)
+  verbose && println("   Exiting CG solver.")
+
+  return x
+end
diff --git a/src/LinearSolvers/Krylov/FGMRESSolvers.jl b/src/LinearSolvers/Krylov/FGMRESSolvers.jl
new file mode 100644
index 00000000..4aeae87b
--- /dev/null
+++ b/src/LinearSolvers/Krylov/FGMRESSolvers.jl
@@ -0,0 +1,127 @@
+
+# FGMRES Solver
+struct FGMRESSolver <: Gridap.Algebra.LinearSolver
+  m       :: Int
+  Pr      :: Gridap.Algebra.LinearSolver
+  Pl      :: Union{Gridap.Algebra.LinearSolver,Nothing}
+  atol    :: Float64
+  rtol    :: Float64
+  verbose :: Bool
+end
+
+function FGMRESSolver(m,Pr;Pl=nothing,atol=1e-12,rtol=1.e-6,verbose=false)
+  return FGMRESSolver(m,Pr,Pl,atol,rtol,verbose)
+end
+
+struct FGMRESSymbolicSetup <: Gridap.Algebra.SymbolicSetup
+  solver
+end
+
+function Gridap.Algebra.symbolic_setup(solver::FGMRESSolver, A::AbstractMatrix)
+  return FGMRESSymbolicSetup(solver)
+end
+
+mutable struct FGMRESNumericalSetup <: Gridap.Algebra.NumericalSetup
+  solver
+  A
+  Pr_ns
+  Pl_ns
+  caches
+end
+
+function get_solver_caches(solver::FGMRESSolver,A)
+  m = solver.m; Pl = solver.Pl
+
+  V  = [allocate_col_vector(A) for i in 1:m+1]
+  Z  = [allocate_col_vector(A) for i in 1:m]
+  zl = !isa(Pl,Nothing) ? allocate_col_vector(A) : nothing
+
+  H = zeros(m+1,m)  # Hessenberg matrix
+  g = zeros(m+1)    # Residual vector
+  c = zeros(m)      # Gibens rotation cosines
+  s = zeros(m)      # Gibens rotation sines
+  return (V,Z,zl,H,g,c,s)
+end
+
+function Gridap.Algebra.numerical_setup(ss::FGMRESSymbolicSetup, A::AbstractMatrix)
+  solver = ss.solver
+  Pr_ns  = numerical_setup(symbolic_setup(solver.Pr,A),A)
+  Pl_ns  = isa(solver.Pl,Nothing) ? nothing : numerical_setup(symbolic_setup(solver.Pl,A),A)
+  caches = get_solver_caches(solver,A)
+  return FGMRESNumericalSetup(solver,A,Pr_ns,Pl_ns,caches)
+end
+
+function Gridap.Algebra.numerical_setup!(ns::FGMRESNumericalSetup, A::AbstractMatrix)
+  numerical_setup!(ns.Pr_ns,A)
+  if !isa(ns.Pl_ns,Nothing)
+    numerical_setup!(ns.Pl_ns,A)
+  end
+  ns.A = A
+end
+
+function Gridap.Algebra.solve!(x::AbstractVector,ns::FGMRESNumericalSetup,b::AbstractVector)
+  solver, A, Pl, Pr, caches = ns.solver, ns.A, ns.Pl_ns, ns.Pr_ns, ns.caches
+  m, atol, rtol, verbose = solver.m, solver.atol, solver.rtol, solver.verbose
+  V, Z, zl, H, g, c, s = caches
+  verbose && println(" > Starting FGMRES solver: ")
+
+  # Initial residual
+  krylov_residual!(V[1],x,A,b,Pl,zl)
+
+  iter = 0
+  β    = norm(V[1]); β0 = β
+  converged = (β < atol || β < rtol*β0)
+  while !converged
+    verbose && println("   > Iteration ", iter," - Residual: ", β)
+    fill!(H,0.0)
+    
+    # Arnoldi process
+    j = 1
+    V[1] ./= β
+    fill!(g,0.0); g[1] = β
+    while ( j < m+1 && !converged )
+      verbose && println("      > Inner iteration ", j," - Residual: ", β)
+      # Arnoldi orthogonalization by Modified Gram-Schmidt
+      krylov_mul!(V[j+1],A,V[j],Pr,Pl,Z[j],zl)
+      for i in 1:j
+        H[i,j] = dot(V[j+1],V[i])
+        V[j+1] .= V[j+1] .- H[i,j] .* V[i]
+      end
+      H[j+1,j] = norm(V[j+1])
+      V[j+1] ./= H[j+1,j]
+
+      # Update QR
+      for i in 1:j-1
+        γ = c[i]*H[i,j] + s[i]*H[i+1,j]
+        H[i+1,j] = -s[i]*H[i,j] + c[i]*H[i+1,j]
+        H[i,j] = γ
+      end
+
+      # New Givens rotation, update QR and residual
+      c[j], s[j], _ = LinearAlgebra.givensAlgorithm(H[j,j],H[j+1,j])
+      H[j,j] = c[j]*H[j,j] + s[j]*H[j+1,j]; H[j+1,j] = 0.0
+      g[j+1] = -s[j]*g[j]; g[j] = c[j]*g[j]
+
+      β  = abs(g[j+1]); converged = (β < atol || β < rtol*β0)
+      j += 1
+    end
+    j = j-1
+
+    # Solve least squares problem Hy = g by backward substitution
+    for i in j:-1:1
+      g[i] = (g[i] - dot(H[i,i+1:j],g[i+1:j])) / H[i,i]
+    end
+
+    # Update solution & residual
+    for i in 1:j
+      x .+= g[i] .* Z[i]
+    end
+    krylov_residual!(V[1],x,A,b,Pl,zl)
+
+    iter += 1
+  end
+  verbose && println("   > Num Iter: ", iter," - Final residual: ", β)
+  verbose && println("   Exiting FGMRES solver.")
+
+  return x
+end
diff --git a/src/LinearSolvers/GMRESSolvers.jl b/src/LinearSolvers/Krylov/GMRESSolvers.jl
similarity index 57%
rename from src/LinearSolvers/GMRESSolvers.jl
rename to src/LinearSolvers/Krylov/GMRESSolvers.jl
index fc2de310..92234193 100644
--- a/src/LinearSolvers/GMRESSolvers.jl
+++ b/src/LinearSolvers/Krylov/GMRESSolvers.jl
@@ -1,15 +1,15 @@
-
 # GMRES Solver
 struct GMRESSolver <: Gridap.Algebra.LinearSolver
-  m   ::Int
-  Pl  ::Gridap.Algebra.LinearSolver
-  atol::Float64
-  rtol::Float64
-  verbose::Bool
+  m       :: Int
+  Pr      :: Union{Gridap.Algebra.LinearSolver,Nothing}
+  Pl      :: Union{Gridap.Algebra.LinearSolver,Nothing}
+  atol    :: Float64
+  rtol    :: Float64
+  verbose :: Bool
 end
 
-function GMRESSolver(m,Pl;atol=1e-12,rtol=1.e-6,verbose=false)
-  return GMRESSolver(m,Pl,atol,rtol,verbose)
+function GMRESSolver(m;Pr=nothing,Pl=nothing,atol=1e-12,rtol=1.e-6,verbose=false)
+  return GMRESSolver(m,Pr,Pl,atol,rtol,verbose)
 end
 
 struct GMRESSymbolicSetup <: Gridap.Algebra.SymbolicSetup
@@ -23,65 +23,72 @@ end
 mutable struct GMRESNumericalSetup <: Gridap.Algebra.NumericalSetup
   solver
   A
+  Pr_ns
   Pl_ns
   caches
 end
 
-function get_gmres_caches(m,A)
-  w = allocate_col_vector(A)
-  V = [allocate_col_vector(A) for i in 1:m+1]
-  Z = [allocate_col_vector(A) for i in 1:m]
+function get_solver_caches(solver::GMRESSolver,A)
+  m, Pl, Pr = solver.m, solver.Pl, solver.Pr
+
+  V  = [allocate_col_vector(A) for i in 1:m+1]
+  zr = !isa(Pr,Nothing) ? allocate_col_vector(A) : nothing
+  zl = !isa(Pr,Nothing) ? allocate_col_vector(A) : nothing
 
   H = zeros(m+1,m)  # Hessenberg matrix
   g = zeros(m+1)    # Residual vector
   c = zeros(m)      # Gibens rotation cosines
   s = zeros(m)      # Gibens rotation sines
-  return (w,V,Z,H,g,c,s)
+  return (V,zr,zl,H,g,c,s)
 end
 
 function Gridap.Algebra.numerical_setup(ss::GMRESSymbolicSetup, A::AbstractMatrix)
   solver = ss.solver
-  Pl_ns  = numerical_setup(symbolic_setup(solver.Pl,A),A)
-  caches = get_gmres_caches(solver.m,A)
-  return GMRESNumericalSetup(solver,A,Pl_ns,caches)
+  Pr_ns  = isa(solver.Pr,Nothing) ? nothing : numerical_setup(symbolic_setup(solver.Pr,A),A)
+  Pl_ns  = isa(solver.Pl,Nothing) ? nothing : numerical_setup(symbolic_setup(solver.Pl,A),A)
+  caches = get_solver_caches(solver,A)
+  return GMRESNumericalSetup(solver,A,Pr_ns,Pl_ns,caches)
 end
 
 function Gridap.Algebra.numerical_setup!(ns::GMRESNumericalSetup, A::AbstractMatrix)
-  numerical_setup!(ns.Pl_ns,A)
+  if !isa(ns.Pr_ns,Nothing)
+    numerical_setup!(ns.Pr_ns,A)
+  end
+  if !isa(ns.Pl_ns,Nothing)
+    numerical_setup!(ns.Pl_ns,A)
+  end
   ns.A = A
 end
 
 function Gridap.Algebra.solve!(x::AbstractVector,ns::GMRESNumericalSetup,b::AbstractVector)
-  solver, A, Pl, caches = ns.solver, ns.A, ns.Pl_ns, ns.caches
+  solver, A, Pl, Pr, caches = ns.solver, ns.A, ns.Pl_ns, ns.Pr_ns, ns.caches
   m, atol, rtol, verbose = solver.m, solver.atol, solver.rtol, solver.verbose
-  w, V, Z, H, g, c, s = caches
+  V, zr, zl, H, g, c, s = caches
   verbose && println(" > Starting GMRES solver: ")
 
   # Initial residual
-  mul!(w,A,x); w .= b .- w
-
-  β    = norm(w); β0 = β
-  converged = (β < atol || β < rtol*β0)
+  krylov_residual!(V[1],x,A,b,Pl,zl)
+  β    = norm(V[1]); β0 = β
   iter = 0
+  converged = (β < atol || β < rtol*β0)
   while !converged
     verbose && println("   > Iteration ", iter," - Residual: ", β)
     fill!(H,0.0)
     
     # Arnoldi process
-    fill!(g,0.0); g[1] = β
-    V[1] .= w ./ β
     j = 1
+    V[1] ./= β
+    fill!(g,0.0); g[1] = β
     while ( j < m+1 && !converged )
       verbose && println("      > Inner iteration ", j," - Residual: ", β)
       # Arnoldi orthogonalization by Modified Gram-Schmidt
-      solve!(Z[j],Pl,V[j])
-      mul!(w,A,Z[j])
+      krylov_mul!(V[j+1],A,V[j],Pr,Pl,zr,zl)
       for i in 1:j
-        H[i,j] = dot(w,V[i])
-        w .= w .- H[i,j] .* V[i]
+        H[i,j] = dot(V[j+1],V[i])
+        V[j+1] .= V[j+1] .- H[i,j] .* V[i]
       end
-      H[j+1,j] = norm(w)
-      V[j+1] = w ./ H[j+1,j]
+      H[j+1,j] = norm(V[j+1])
+      V[j+1] ./= H[j+1,j]
 
       # Update QR
       for i in 1:j-1
@@ -106,10 +113,19 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::GMRESNumericalSetup,b::Abst
     end
 
     # Update solution & residual
-    for i in 1:j
-      x .+= g[i] .* Z[i]
+    if isa(Pr,Nothing)
+      for i in 1:j
+        x .+= g[i] .* V[i]
+      end
+    else
+      fill!(zl,0.0)
+      for i in 1:j
+        zl .+= g[i] .* V[i]
+      end
+      solve!(zr,Pr,zl)
+      x .+= zr
     end
-    mul!(w,A,x); w .= b .- w
+    krylov_residual!(V[1],x,A,b,Pl,zl)
 
     iter += 1
   end
@@ -117,4 +133,4 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::GMRESNumericalSetup,b::Abst
   verbose && println("   Exiting GMRES solver.")
 
   return x
-end
+end
\ No newline at end of file
diff --git a/src/LinearSolvers/Krylov/KrylovUtils.jl b/src/LinearSolvers/Krylov/KrylovUtils.jl
new file mode 100644
index 00000000..c8cde945
--- /dev/null
+++ b/src/LinearSolvers/Krylov/KrylovUtils.jl
@@ -0,0 +1,37 @@
+
+"""
+  Computes the Krylov matrix-vector product y = Pl⁻¹⋅A⋅Pr⁻¹⋅x
+  by solving: 
+    Pr⋅wr = x
+    wl = A⋅wr
+    Pl⋅y = wl
+"""
+function krylov_mul!(y,A,x,Pr,Pl,wr,wl)
+  solve!(wr,Pr,x)
+  mul!(wl,A,wr)
+  solve!(y,Pl,wl)
+end
+function krylov_mul!(y,A,x,Pr,Pl::Nothing,wr,wl)
+  solve!(wr,Pr,x)
+  mul!(y,A,wr)
+end
+function krylov_mul!(y,A,x,Pr::Nothing,Pl,wr,wl)
+  mul!(wl,A,x)
+  solve!(y,Pl,wl)
+end
+
+"""
+  Computes the Krylov residual r = Pl⁻¹(A⋅x - b).
+  by solving: 
+    w = A⋅x - b
+    Pl⋅r = w
+"""
+function krylov_residual!(r,x,A,b,Pl,w)
+  mul!(w,A,x)
+  w .= b .- w
+  solve!(r,Pl,w)
+end
+function krylov_residual!(r,x,A,b,Pl::Nothing,w::Nothing)
+  mul!(r,A,x)
+  r .= b .- r
+end
diff --git a/src/LinearSolvers/Krylov/MINRESSolvers.jl b/src/LinearSolvers/Krylov/MINRESSolvers.jl
new file mode 100644
index 00000000..0f010aa1
--- /dev/null
+++ b/src/LinearSolvers/Krylov/MINRESSolvers.jl
@@ -0,0 +1,118 @@
+# MINRES Solver
+struct MINRESSolver <: Gridap.Algebra.LinearSolver
+  Pr      :: Union{Gridap.Algebra.LinearSolver,Nothing}
+  Pl      :: Union{Gridap.Algebra.LinearSolver,Nothing}
+  atol    :: Float64
+  rtol    :: Float64
+  verbose :: Bool
+end
+
+function MINRESSolver(;Pr=nothing,Pl=nothing,atol=1e-12,rtol=1.e-6,verbose=false)
+  return MINRESSolver(Pr,Pl,atol,rtol,verbose)
+end
+
+struct MINRESSymbolicSetup <: Gridap.Algebra.SymbolicSetup
+  solver
+end
+
+function Gridap.Algebra.symbolic_setup(solver::MINRESSolver, A::AbstractMatrix)
+  return MINRESSymbolicSetup(solver)
+end
+
+mutable struct MINRESNumericalSetup <: Gridap.Algebra.NumericalSetup
+  solver
+  A
+  Pr_ns
+  Pl_ns
+  caches
+end
+
+function get_solver_caches(solver::MINRESSolver,A)
+  Pl, Pr = solver.Pl, solver.Pr
+
+  V  = [allocate_col_vector(A) for i in 1:3]
+  W  = [allocate_col_vector(A) for i in 1:3]
+  zr = !isa(Pr,Nothing) ? allocate_col_vector(A) : nothing
+  zl = !isa(Pl,Nothing) ? allocate_col_vector(A) : nothing
+
+  H = zeros(4) # Hessenberg matrix
+  g = zeros(2) # Residual vector
+  c = zeros(2) # Gibens rotation cosines
+  s = zeros(2) # Gibens rotation sines
+  return (V,W,zr,zl,H,g,c,s)
+end
+
+function Gridap.Algebra.numerical_setup(ss::MINRESSymbolicSetup, A::AbstractMatrix)
+  solver = ss.solver
+  Pr_ns  = isa(solver.Pr,Nothing) ? nothing : numerical_setup(symbolic_setup(solver.Pr,A),A)
+  Pl_ns  = isa(solver.Pl,Nothing) ? nothing : numerical_setup(symbolic_setup(solver.Pl,A),A)
+  caches = get_solver_caches(solver,A)
+  return MINRESNumericalSetup(solver,A,Pr_ns,Pl_ns,caches)
+end
+
+function Gridap.Algebra.numerical_setup!(ns::MINRESNumericalSetup, A::AbstractMatrix)
+  if !isa(ns.Pr_ns,Nothing)
+    numerical_setup!(ns.Pr_ns,A)
+  end
+  if !isa(ns.Pl_ns,Nothing)
+    numerical_setup!(ns.Pl_ns,A)
+  end
+  ns.A = A
+end
+
+function Gridap.Algebra.solve!(x::AbstractVector,ns::MINRESNumericalSetup,b::AbstractVector)
+  solver, A, Pl, Pr, caches = ns.solver, ns.A, ns.Pl_ns, ns.Pr_ns, ns.caches
+  atol, rtol, verbose = solver.atol, solver.rtol, solver.verbose
+  V, W, zr, zl, H, g, c, s = caches
+  verbose && println(" > Starting MINRES solver: ")
+
+  Vjm1, Vj, Vjp1 = V
+  Wjm1, Wj, Wjp1 = W
+
+  fill!(Vjm1,0.0); fill!(Vjp1,0.0); copy!(Vj,b)
+  fill!(Wjm1,0.0); fill!(Wjp1,0.0); fill!(Wj,0.0)
+  fill!(H,0.0), fill!(c,1.0); fill!(s,0.0); fill!(g,0.0)
+
+  krylov_residual!(Vj,x,A,b,Pl,zl)
+  β    = norm(Vj); β0 = β; Vj ./= β; g[1] = β
+  iter = 0
+  converged = (β < atol || β < rtol*β0)
+  while !converged
+    verbose && println("   > Iteration ", iter," - Residual: ", β)
+
+    krylov_mul!(Vjp1,A,Vj,Pr,Pl,zr,zl)
+    H[3] = dot(Vjp1,Vj)
+    Vjp1 .= Vjp1 .- H[3] .* Vj .- H[2] .* Vjm1
+    H[4] = norm(Vjp1)
+    Vjp1 ./= H[4]
+
+    # Update QR
+    H[1] = s[1]*H[2]; H[2] = c[1]*H[2]
+    γ = c[2]*H[2] + s[2]*H[3]; H[3] = -s[2]*H[2] + c[2]*H[3]; H[2] = γ
+
+    # New Givens rotation, update QR and residual
+    c[1], s[1] = c[2], s[2]
+    c[2], s[2], H[3] = LinearAlgebra.givensAlgorithm(H[3],H[4])
+    g[2] = -s[2]*g[1]; g[1] = c[2]*g[1]
+
+    # Update solution
+    Wjp1 .= Vj .- H[2] .* Wj .- H[1] .* Wjm1
+    Wjp1 ./= H[3]
+    if isa(Pr,Nothing)
+      x .+= g[1] .* Wjp1
+    else
+      solve!(zr,Pr,Wjp1)
+      x .+= g[1] .* zr
+    end
+
+    β  = abs(g[2]); converged = (β < atol || β < rtol*β0)
+    Vjm1, Vj, Vjp1 = Vj, Vjp1, Vjm1
+    Wjm1, Wj, Wjp1 = Wj, Wjp1, Wjm1
+    g[1] = g[2]; H[2] = H[4];
+    iter += 1
+  end
+  verbose && println("   > Num Iter: ", iter," - Final residual: ", β)
+  verbose && println("   Exiting MINRES solver.")
+
+  return x
+end
\ No newline at end of file
diff --git a/src/LinearSolvers/LinearSolvers.jl b/src/LinearSolvers/LinearSolvers.jl
index 14477214..0db43168 100644
--- a/src/LinearSolvers/LinearSolvers.jl
+++ b/src/LinearSolvers/LinearSolvers.jl
@@ -23,6 +23,7 @@ export RichardsonSmoother
 export SymGaussSeidelSmoother
 export GMGLinearSolver
 export BlockDiagonalSmoother
+export SchurComplementSolver
 
 # Wrappers for IterativeSolvers.jl
 export IS_ConjugateGradientSolver
@@ -30,8 +31,17 @@ export IS_GMRESSolver
 export IS_MINRESSolver
 export IS_SSORSolver
 
+# Krylov solvers 
+export CGSolver
 export GMRESSolver
-export SchurComplementSolver
+export FGMRESSolver
+export MINRESSolver
+
+include("Krylov/KrylovUtils.jl")
+include("Krylov/CGSolvers.jl")
+include("Krylov/GMRESSolvers.jl")
+include("Krylov/FGMRESSolvers.jl")
+include("Krylov/MINRESSolvers.jl")
 
 include("IdentityLinearSolvers.jl")
 include("JacobiLinearSolvers.jl")
@@ -40,7 +50,6 @@ include("SymGaussSeidelSmoothers.jl")
 include("GMGLinearSolvers.jl")
 include("BlockDiagonalSmoothers.jl")
 include("IterativeLinearSolvers.jl")
-include("GMRESSolvers.jl")
 include("SchurComplementSolvers.jl")
 
 end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index fb4835f7..4fd2e586 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -65,7 +65,7 @@ run_tests(joinpath(@__DIR__, "mpi"))
 # Sequential tests
 @time @testset "BlockDiagonalSmoothersTests" begin include("seq/BlockDiagonalSmoothersTests.jl") end
 @time @testset "DistributedPatchFESpacesTests" begin include("seq/DistributedPatchFESpacesTests.jl") end
-@time @testset "GMRESSolversTests" begin include("seq/GMRESSolversTests.jl") end
+@time @testset "KrylovSolversTests" begin include("seq/KrylovSolversTests.jl") end
 @time @testset "IterativeSolversTests" begin include("seq/IterativeSolversTests.jl") end
 @time @testset "PatchLinearSolverTests" begin include("seq/PatchLinearSolverTests.jl") end
 @time @testset "SymGaussSeidelSmoothersTests" begin include("seq/SymGaussSeidelSmoothersTests.jl") end
diff --git a/test/seq/GMRESSolversTests.jl b/test/seq/KrylovSolversTests.jl
similarity index 65%
rename from test/seq/GMRESSolversTests.jl
rename to test/seq/KrylovSolversTests.jl
index 146d5e17..d5452065 100644
--- a/test/seq/GMRESSolversTests.jl
+++ b/test/seq/KrylovSolversTests.jl
@@ -12,6 +12,20 @@ using GridapSolvers.LinearSolvers
 sol(x) = x[1] + x[2]
 f(x)   = -Δ(sol)(x)
 
+function test_solver(solver,op,Uh,dΩ)
+  A, b = get_matrix(op), get_vector(op);
+  ns = numerical_setup(symbolic_setup(solver,A),A)
+
+  x = LinearSolvers.allocate_col_vector(A)
+  solve!(x,ns,b)
+
+  u  = interpolate(sol,Uh)
+  uh = FEFunction(Uh,x)
+  eh = uh - u
+  E  = sum(∫(eh*eh)*dΩ)
+  @test E < 1.e-8
+end
+
 function main(model)
   order  = 1
   qorder = order*2 + 1
@@ -24,29 +38,31 @@ function main(model)
   dΩ     = Measure(Ω,qorder)
   a(u,v) = ∫(∇(v)⋅∇(u))*dΩ
   l(v)   = ∫(v⋅f)*dΩ
-
   op = AffineFEOperator(a,l,Uh,Vh)
-  A, b = get_matrix(op), get_vector(op);
+  
+  P = JacobiLinearSolver()
 
-  Pl = JacobiLinearSolver()
-  solver = LinearSolvers.GMRESSolver(40,Pl,1.e-8)
-  ns = numerical_setup(symbolic_setup(solver,A),A)
+  gmres = LinearSolvers.GMRESSolver(40;Pr=P,Pl=P,rtol=1.e-8,verbose=true)
+  test_solver(gmres,op,Uh,dΩ)
 
-  x = LinearSolvers.allocate_col_vector(A)
-  solve!(x,ns,b)
+  fgmres = LinearSolvers.FGMRESSolver(40,P;rtol=1.e-8,verbose=true)
+  test_solver(fgmres,op,Uh,dΩ)
 
-  u  = interpolate(sol,Uh)
-  uh = FEFunction(Uh,x)
-  eh = uh - u
-  E  = sum(∫(eh*eh)*dΩ)
-  return E < 1.e-8
+  pcg = LinearSolvers.CGSolver(P;rtol=1.e-8,verbose=true)
+  test_solver(pcg,op,Uh,dΩ)
+
+  fpcg = LinearSolvers.CGSolver(P;flexible=true,rtol=1.e-8,verbose=true)
+  test_solver(fpcg,op,Uh,dΩ)
+
+  minres = LinearSolvers.MINRESSolver(;Pl=P,Pr=P,rtol=1.e-8,verbose=true)
+  test_solver(minres,op,Uh,dΩ)
 end
 
 # Completely serial
 mesh_partition = (10,10)
 domain = (0,1,0,1)
 model  = CartesianDiscreteModel(domain,mesh_partition)
-@test main(model)
+main(model)
 
 # Sequential
 num_ranks = (1,2)
@@ -55,6 +71,6 @@ parts = with_debug() do distribute
 end
 
 model  = CartesianDiscreteModel(parts,num_ranks,domain,mesh_partition)
-@test main(model)
+main(model)
 
 end
\ No newline at end of file