diff --git a/Project.toml b/Project.toml
index 474c235e..c8dc0949 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,6 +5,7 @@ version = "0.2.0"
 
 [deps]
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
+BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 GridapDistributed = "f9701e48-63b3-45aa-9a63-9bc6c271f355"
@@ -15,6 +16,8 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1"
 
 [compat]
 ArgParse = "1"
diff --git a/src/GridapSolvers.jl b/src/GridapSolvers.jl
index 988c541f..87294562 100644
--- a/src/GridapSolvers.jl
+++ b/src/GridapSolvers.jl
@@ -26,7 +26,16 @@ module GridapSolvers
   # LinearSolvers
   export JacobiLinearSolver
   export RichardsonSmoother
+  export SymGaussSeidelSmoother
   export GMGLinearSolver
+  export BlockDiagonalSmoother
+
+  export ConjugateGradientSolver
+  export IS_GMRESSolver
+  export IS_MINRESSolver
+  export IS_SSORSolver
+
+  export GMRESSolver
 
   # PatchBasedSmoothers
   export PatchDecomposition
diff --git a/src/LinearSolvers/BlockDiagonalSmoothers.jl b/src/LinearSolvers/BlockDiagonalSmoothers.jl
new file mode 100644
index 00000000..8e34511e
--- /dev/null
+++ b/src/LinearSolvers/BlockDiagonalSmoothers.jl
@@ -0,0 +1,80 @@
+struct BlockDiagonalSmoother{A,B} <: Gridap.Algebra.LinearSolver
+  blocks  :: A
+  solvers :: B
+  function BlockDiagonalSmoother(blocks ::AbstractArray{<:AbstractMatrix},
+                                 solvers::AbstractArray{<:Gridap.Algebra.LinearSolver})
+    @check length(blocks) == length(solvers)
+    A = typeof(blocks)
+    B = typeof(solvers)
+    return new{A,B}(blocks,solvers)
+  end
+end
+
+# Constructors 
+
+function BlockDiagonalSmoother(block_mat :: AbstractBlockMatrix,
+                               solvers   :: AbstractArray{<:Gridap.Algebra.LinearSolver})
+  mat_blocks = diag(blocks(block_mat))
+  return BlockDiagonalSmoother(mat_blocks,solvers)
+end
+
+function BlockDiagonalSmoother(biforms :: AbstractArray{<:Function},
+                               trials  :: AbstractArray{<:FESpace},
+                               tests   :: AbstractArray{<:FESpace},
+                               solvers :: AbstractArray{<:Gridap.Algebra.LinearSolver})
+  mat_blocks = compute_block_matrices(biforms,trials,tests)
+  return BlockDiagonalSmoother(mat_blocks,solvers)
+end
+
+function BlockDiagonalSmoother(biforms :: AbstractArray{<:Function},
+                               U       :: Union{MultiFieldFESpace,GridapDistributed.DistributedMultiFieldFESpace},
+                               V       :: Union{MultiFieldFESpace,GridapDistributed.DistributedMultiFieldFESpace},
+                               solvers :: AbstractArray{<:Gridap.Algebra.LinearSolver})
+  return BlockDiagonalSmoother(biforms,[U...],[V...],solvers)
+end
+
+function compute_block_matrices(biforms :: AbstractArray{<:Function},
+                                trials  :: AbstractArray{<:FESpace},
+                                tests   :: AbstractArray{<:FESpace})
+  @check length(biforms) == length(tests) == length(trials)
+  mat_blocks = map(assemble_matrix,biforms,tests,trials)
+  return mat_blocks
+end
+
+# Symbolic and numerical setup
+struct BlockDiagonalSmootherSS{A,B} <: Gridap.Algebra.SymbolicSetup
+  solver   :: A
+  block_ss :: B
+end
+
+function Gridap.Algebra.symbolic_setup(solver::BlockDiagonalSmoother,mat::AbstractMatrix)
+  block_ss = map(symbolic_setup,solver.solvers,solver.blocks)
+  return BlockDiagonalSmootherSS(solver,block_ss)
+end
+
+struct BlockDiagonalSmootherNS{A,B} <: Gridap.Algebra.NumericalSetup
+  solver   :: A
+  block_ns :: B
+end
+
+function Gridap.Algebra.numerical_setup(ss::BlockDiagonalSmootherSS,mat::AbstractMatrix)
+  solver   = ss.solver
+  block_ns = map(numerical_setup,ss.block_ss,solver.blocks)
+  return BlockDiagonalSmootherNS(solver,block_ns)
+end
+
+# Solve
+
+function Gridap.Algebra.solve!(x::AbstractBlockVector,ns::BlockDiagonalSmootherNS,b::AbstractBlockVector)
+  @check blocklength(x) == blocklength(b) == length(ns.block_ns)
+  for (iB,bns) in enumerate(ns.block_ns)
+    xi = x[Block(iB)]
+    bi = b[Block(iB)]
+    solve!(xi,bns,bi)
+  end
+  return x
+end
+
+function LinearAlgebra.ldiv!(x,ns::BlockDiagonalSmootherNS,b)
+  solve!(x,ns,b)
+end
\ No newline at end of file
diff --git a/src/LinearSolvers/GMGLinearSolvers.jl b/src/LinearSolvers/GMGLinearSolvers.jl
index d03ed8d3..11223024 100644
--- a/src/LinearSolvers/GMGLinearSolvers.jl
+++ b/src/LinearSolvers/GMGLinearSolvers.jl
@@ -88,7 +88,7 @@ function setup_finest_level_cache(mh::ModelHierarchy,smatrices::Vector{<:Abstrac
   parts = get_level_parts(mh,1)
   if i_am_in(parts)
     Ah = smatrices[1]
-    rh = pfill(0.0, partition(axes(Ah,2)))
+    rh = allocate_col_vector(Ah)
     cache = rh
   end
   return cache
@@ -156,14 +156,14 @@ function setup_coarsest_solver_cache(mh::ModelHierarchy,coarsest_solver::PETScLi
 end
 
 function allocate_level_work_vectors(mh::ModelHierarchy,smatrices::Vector{<:AbstractMatrix},lev::Integer)
-  dxh   = pfill(0.0,partition(axes(smatrices[lev],2)))
-  Adxh  = pfill(0.0,partition(axes(smatrices[lev],1)))
+  dxh   = allocate_col_vector(smatrices[lev])
+  Adxh  = allocate_row_vector(smatrices[lev])
 
   cparts = get_level_parts(mh,lev+1)
   if i_am_in(cparts)
     AH  = smatrices[lev+1]
-    rH  = pfill(0.0,partition(axes(AH,2)))
-    dxH = pfill(0.0,partition(axes(AH,2)))
+    rH  = allocate_col_vector(AH)
+    dxH = allocate_col_vector(AH)
   else
     rH  = nothing
     dxH = nothing
diff --git a/src/LinearSolvers/GMRESSolvers.jl b/src/LinearSolvers/GMRESSolvers.jl
new file mode 100644
index 00000000..79729f23
--- /dev/null
+++ b/src/LinearSolvers/GMRESSolvers.jl
@@ -0,0 +1,113 @@
+
+# GMRES Solver
+struct GMRESSolver <: Gridap.Algebra.LinearSolver
+  m  ::Int
+  Pl ::Gridap.Algebra.LinearSolver
+  tol::Float64
+end
+
+struct GMRESSymbolicSetup <: Gridap.Algebra.SymbolicSetup
+  solver
+end
+
+function Gridap.Algebra.symbolic_setup(solver::GMRESSolver, A::AbstractMatrix)
+  return GMRESSymbolicSetup(solver)
+end
+
+mutable struct GMRESNumericalSetup <: Gridap.Algebra.NumericalSetup
+  solver
+  A
+  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]
+
+  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)
+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)
+end
+
+function Gridap.Algebra.numerical_setup!(ns::GMRESNumericalSetup, A::AbstractMatrix)
+  numerical_setup!(ns.Pl_ns,A)
+  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
+  m, tol = solver.m, solver.tol
+  w, V, Z, H, g, c, s = caches
+  println(" > Starting GMRES solver: ")
+
+  # Initial residual
+  mul!(w,A,x); w .= b .- w
+
+  β    = norm(w)
+  iter = 0
+  while (β > tol)
+    println("   > Iteration ", iter," - Residual: ", β)
+    fill!(H,0.0)
+    
+    # Arnoldi process
+    fill!(g,0.0); g[1] = β
+    V[1] .= w ./ β
+    j = 1
+    while ( j < m+1 && β > tol )
+      println("      > Inner iteration ", j," - Residual: ", β)
+      # Arnoldi orthogonalization by Modified Gram-Schmidt
+      solve!(Z[j],Pl,V[j])
+      mul!(w,A,Z[j])
+      for i in 1:j
+        H[i,j] = dot(w,V[i])
+        w .= w .- H[i,j] .* V[i]
+      end
+      H[j+1,j] = norm(w)
+      V[j+1] = w ./ 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])
+      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
+    mul!(w,A,x); w .= b .- w
+
+    iter += 1
+  end
+  println("   Exiting GMRES solver.")
+  println("   > Num Iter: ", iter," - Final residual: ", β)
+
+  return x
+end
diff --git a/src/LinearSolvers/Helpers.jl b/src/LinearSolvers/Helpers.jl
new file mode 100644
index 00000000..013ee4f9
--- /dev/null
+++ b/src/LinearSolvers/Helpers.jl
@@ -0,0 +1,29 @@
+
+# Row/Col vector allocations for serial
+function allocate_row_vector(A::AbstractMatrix{T}) where T
+  return zeros(T,size(A,1))
+end
+
+function allocate_col_vector(A::AbstractMatrix{T}) where T
+  return zeros(T,size(A,2))
+end
+
+# Row/Col vector allocations for parallel
+function allocate_row_vector(A::PSparseMatrix)
+  T = eltype(A)
+  return pfill(zero(T),partition(axes(A,1)))
+end
+
+function allocate_col_vector(A::PSparseMatrix)
+  T = eltype(A)
+  return pfill(zero(T),partition(axes(A,2)))
+end
+
+# Row/Col vector allocations for blocks
+function allocate_row_vector(A::AbstractBlockMatrix)
+  return mortar(map(Aii->allocate_row_vector(Aii),blocks(A)[:,1]))
+end
+
+function allocate_col_vector(A::AbstractBlockMatrix)
+  return mortar(map(Aii->allocate_col_vector(Aii),blocks(A)[1,:]))
+end
diff --git a/src/LinearSolvers/IdentityLinearSolvers.jl b/src/LinearSolvers/IdentityLinearSolvers.jl
new file mode 100644
index 00000000..93a24d64
--- /dev/null
+++ b/src/LinearSolvers/IdentityLinearSolvers.jl
@@ -0,0 +1,26 @@
+
+# Identity solver, for testing purposes 
+struct IdentitySolver <: Gridap.Algebra.LinearSolver
+end
+
+struct IdentitySymbolicSetup <: Gridap.Algebra.SymbolicSetup
+  solver
+end
+
+function Gridap.Algebra.symbolic_setup(s::IdentitySolver,A::AbstractMatrix) 
+  IdentitySymbolicSetup(s)
+end
+
+struct IdentityNumericalSetup <: Gridap.Algebra.NumericalSetup
+  solver
+end
+
+function Gridap.Algebra.numerical_setup(ss::IdentitySymbolicSetup,mat::AbstractMatrix)
+  s = ss.solver
+  return IdentityNumericalSetup(s)
+end
+
+function Gridap.Algebra.solve!(x::AbstractVector,ns::IdentityNumericalSetup,y::AbstractVector)
+  copy!(x,y)
+  return x
+end
diff --git a/src/LinearSolvers/IterativeLinearSolvers.jl b/src/LinearSolvers/IterativeLinearSolvers.jl
new file mode 100644
index 00000000..49d5ff22
--- /dev/null
+++ b/src/LinearSolvers/IterativeLinearSolvers.jl
@@ -0,0 +1,180 @@
+
+abstract type IterativeLinearSolverType end
+struct CGIterativeSolverType     <: IterativeLinearSolverType end
+struct GMRESIterativeSolverType  <: IterativeLinearSolverType end
+struct MINRESIterativeSolverType <: IterativeLinearSolverType end
+struct SSORIterativeSolverType   <: IterativeLinearSolverType end
+
+# Constructors
+
+"""
+  Wrappers for [IterativeSolvers.jl](https://github.com/JuliaLinearAlgebra/IterativeSolvers.jl)
+  krylov-like iterative solvers.
+
+  Currently supported: 
+  - ConjugateGradientSolver
+  - GMRESSolver
+  - MINRESSolver
+"""
+struct IterativeLinearSolver{A} <: Gridap.Algebra.LinearSolver
+  args
+  kwargs
+
+  function IterativeLinearSolver(type::IterativeLinearSolverType,args,kwargs)
+    A = typeof(type)
+    return new{A}(args,kwargs)
+  end
+end
+
+SolverType(::IterativeLinearSolver{T}) where T = T()
+
+function IS_ConjugateGradientSolver(;kwargs...)
+  options = [:statevars,:initially_zero,:Pl,:abstol,:reltol,:maxiter,:verbose,:log]
+  @check all(map(opt -> opt ∈ options,keys(kwargs)))
+  return IterativeLinearSolver(CGIterativeSolverType(),nothing,kwargs)
+end
+
+function IS_GMRESSolver(;kwargs...)
+  options = [:initially_zero,:abstol,:reltol,:restart,:maxiter,:Pl,:Pr,:log,:verbose,:orth_meth]
+  @check all(map(opt -> opt ∈ options,keys(kwargs)))
+  return IterativeLinearSolver(GMRESIterativeSolverType(),nothing,kwargs)
+end
+
+function IS_MINRESSolver(;kwargs...)
+  options = [:initially_zero,:skew_hermitian,:abstol,:reltol,:maxiter,:log,:verbose]
+  @check all(map(opt -> opt ∈ options,keys(kwargs)))
+  return IterativeLinearSolver(MINRESIterativeSolverType(),nothing,kwargs)
+end
+
+function IS_SSORSolver(ω::Real;kwargs...)
+  options = [:maxiter]
+  @check all(map(opt -> opt ∈ options,keys(kwargs)))
+  args = Dict(:ω => ω)
+  return IterativeLinearSolver(SSORIterativeSolverType(),args,kwargs)
+end
+
+# Symbolic setup
+
+struct IterativeLinearSolverSS <: Gridap.Algebra.SymbolicSetup
+  solver
+end
+
+function Gridap.Algebra.symbolic_setup(solver::IterativeLinearSolver,A::AbstractMatrix)
+  IterativeLinearSolverSS(solver)
+end
+
+# Numerical setup
+
+struct IterativeLinearSolverNS <: Gridap.Algebra.NumericalSetup
+  solver
+  A
+  caches
+end
+
+function Gridap.Algebra.numerical_setup(ss::IterativeLinearSolverSS,A::AbstractMatrix)
+  solver_type = SolverType(ss.solver)
+  numerical_setup(solver_type,ss,A)
+end
+
+function Gridap.Algebra.numerical_setup(::IterativeLinearSolverType,
+                                        ss::IterativeLinearSolverSS,
+                                        A::AbstractMatrix)
+  IterativeLinearSolverNS(ss.solver,A,nothing)
+end
+
+function Gridap.Algebra.numerical_setup(::SSORIterativeSolverType,
+                                        ss::IterativeLinearSolverSS,
+                                        A::AbstractMatrix)
+  x = allocate_row_vector(A)
+  b = allocate_col_vector(A)
+  ω       = ss.solver.args[:ω]
+  maxiter = ss.solver.kwargs[:maxiter]
+  caches  = IterativeSolvers.ssor_iterable(x,A,b,ω;maxiter=maxiter)
+  return IterativeLinearSolverNS(ss.solver,A,caches)
+end
+
+function IterativeSolvers.ssor_iterable(x::PVector,
+                                        A::PSparseMatrix,
+                                        b::PVector,
+                                        ω::Real;
+                                        maxiter::Int = 10)
+  iterables = map(own_values(x),own_values(A),own_values(b)) do _xi,_Aii,_bi
+    xi  = Vector(_xi)
+    Aii = SparseMatrixCSC(_Aii)
+    bi  = Vector(_bi)
+    return IterativeSolvers.ssor_iterable(xi,Aii,bi,ω;maxiter=maxiter)
+  end
+  return iterables
+end
+
+# Solve
+
+function LinearAlgebra.ldiv!(x::AbstractVector,ns::IterativeLinearSolverNS,b::AbstractVector)
+  solve!(x,ns,b)
+end
+
+function Gridap.Algebra.solve!(x::AbstractVector,
+                               ns::IterativeLinearSolverNS,
+                               y::AbstractVector)
+  solver_type = SolverType(ns.solver)
+  solve!(solver_type,x,ns,y)
+end
+
+function Gridap.Algebra.solve!(::IterativeLinearSolverType,
+                               ::AbstractVector,
+                               ::IterativeLinearSolverNS,
+                               ::AbstractVector)
+  @abstractmethod
+end
+
+function Gridap.Algebra.solve!(::CGIterativeSolverType,
+                               x::AbstractVector,
+                               ns::IterativeLinearSolverNS,
+                               y::AbstractVector)
+  A, kwargs = ns.A, ns.solver.kwargs
+  return cg!(x,A,y;kwargs...)
+end
+
+function Gridap.Algebra.solve!(::GMRESIterativeSolverType,
+                               x::AbstractVector,
+                               ns::IterativeLinearSolverNS,
+                               y::AbstractVector)
+  A, kwargs = ns.A, ns.solver.kwargs
+  return gmres!(x,A,y;kwargs...)
+end
+
+function Gridap.Algebra.solve!(::MINRESIterativeSolverType,
+                               x::AbstractVector,
+                               ns::IterativeLinearSolverNS,
+                               y::AbstractVector)
+  A, kwargs = ns.A, ns.solver.kwargs
+  return minres!(x,A,y;kwargs...)
+end
+
+function Gridap.Algebra.solve!(::SSORIterativeSolverType,
+                               x::AbstractVector,
+                               ns::IterativeLinearSolverNS,
+                               y::AbstractVector)
+  iterable = ns.caches
+  iterable.x = x
+  iterable.b = y
+
+  for item = iterable end
+  return x
+end
+
+function Gridap.Algebra.solve!(::SSORIterativeSolverType,
+                               x::PVector,
+                               ns::IterativeLinearSolverNS,
+                               y::PVector)
+  iterables = ns.caches
+  map(iterables,own_values(x),own_values(y)) do iterable, xi, yi
+    iterable.x .= xi
+    iterable.b .= yi
+    for item = iterable end
+    xi .= iterable.x
+    yi .= iterable.b
+  end
+  consistent!(x) |> fetch
+  return x
+end
diff --git a/src/LinearSolvers/JacobiLinearSolvers.jl b/src/LinearSolvers/JacobiLinearSolvers.jl
index 29319cce..5f412d99 100644
--- a/src/LinearSolvers/JacobiLinearSolvers.jl
+++ b/src/LinearSolvers/JacobiLinearSolvers.jl
@@ -22,7 +22,7 @@ function Gridap.Algebra.numerical_setup!(ns::JacobiNumericalSetup, A::AbstractMa
 end
 
 function Gridap.Algebra.numerical_setup(ss::JacobiSymbolicSetup,A::PSparseMatrix)
-  inv_diag=map(own_values(A)) do a
+  inv_diag = map(own_values(A)) do a
     1.0 ./ diag(a)
   end
   return JacobiNumericalSetup(inv_diag)
diff --git a/src/LinearSolvers/LinearSolvers.jl b/src/LinearSolvers/LinearSolvers.jl
index 4fe84ae9..35009b2d 100644
--- a/src/LinearSolvers/LinearSolvers.jl
+++ b/src/LinearSolvers/LinearSolvers.jl
@@ -2,21 +2,49 @@ module LinearSolvers
 
 using Printf
 using LinearAlgebra
+using SparseArrays
+using SparseMatricesCSR
+using BlockArrays
+using IterativeSolvers
+
 using Gridap
+using Gridap.Helpers
 using Gridap.Algebra
+using Gridap.FESpaces
+using Gridap.MultiField
 using PartitionedArrays
 using GridapPETSc
 
+using GridapDistributed
+
 using GridapSolvers.MultilevelTools
 
 import LinearAlgebra: mul!, ldiv!
 
 export JacobiLinearSolver
 export RichardsonSmoother
+export SymGaussSeidelSmoother
 export GMGLinearSolver
+export BlockDiagonalSmoother
+
+# Wrappers for IterativeSolvers.jl
+export IS_ConjugateGradientSolver
+export IS_GMRESSolver
+export IS_MINRESSolver
+export IS_SSORSolver
+
+export GMRESSolver
+export SchurComplementSolver
 
+include("Helpers.jl")
+include("IdentityLinearSolvers.jl")
 include("JacobiLinearSolvers.jl")
 include("RichardsonSmoothers.jl")
+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/src/LinearSolvers/RichardsonSmoothers.jl b/src/LinearSolvers/RichardsonSmoothers.jl
index 39b45fce..86cfb176 100644
--- a/src/LinearSolvers/RichardsonSmoothers.jl
+++ b/src/LinearSolvers/RichardsonSmoothers.jl
@@ -20,28 +20,21 @@ struct RichardsonSmootherSymbolicSetup{A} <: Gridap.Algebra.SymbolicSetup
 end
 
 function Gridap.Algebra.symbolic_setup(smoother::RichardsonSmoother,mat::AbstractMatrix)
-  Mss=symbolic_setup(smoother.M,mat)
+  Mss = symbolic_setup(smoother.M,mat)
   return RichardsonSmootherSymbolicSetup(smoother,Mss)
 end
 
 mutable struct RichardsonSmootherNumericalSetup{A,B,C,D} <: Gridap.Algebra.NumericalSetup
-  smoother       :: RichardsonSmoother
-  A              :: A
-  Adx            :: B
-  dx             :: C
-  Mns            :: D
-end
-
-function Gridap.Algebra.numerical_setup(ss::RichardsonSmootherSymbolicSetup, A::AbstractMatrix{T}) where T
-  Adx = zeros(size(A,1))
-  dx  = zeros(size(A,2))
-  Mns = numerical_setup(ss.Mss,A)
-  return RichardsonSmootherNumericalSetup(ss.smoother,A,Adx,dx,Mns)
+  smoother :: RichardsonSmoother
+  A        :: A
+  Adx      :: B
+  dx       :: C
+  Mns      :: D
 end
 
-function Gridap.Algebra.numerical_setup(ss::RichardsonSmootherSymbolicSetup, A::PSparseMatrix)
-  Adx = pfill(0.0,partition(axes(A,1)))
-  dx  = pfill(0.0,partition(axes(A,2)))
+function Gridap.Algebra.numerical_setup(ss::RichardsonSmootherSymbolicSetup, A::AbstractMatrix)
+  Adx = allocate_row_vector(A)
+  dx  = allocate_col_vector(A)
   Mns = numerical_setup(ss.Mss,A)
   return RichardsonSmootherNumericalSetup(ss.smoother,A,Adx,dx,Mns)
 end
@@ -70,4 +63,3 @@ function LinearAlgebra.ldiv!(x::AbstractVector,ns::RichardsonSmootherNumericalSe
   solve!(x,ns,aux)
   return x
 end
-
diff --git a/src/LinearSolvers/SchurComplementSolvers.jl b/src/LinearSolvers/SchurComplementSolvers.jl
new file mode 100644
index 00000000..d6ae9e5f
--- /dev/null
+++ b/src/LinearSolvers/SchurComplementSolvers.jl
@@ -0,0 +1,73 @@
+
+"""
+  Schur complement solver
+
+  [A B] ^ -1    [Ip  -A^-1 B]   [A^-1     ]   [   Ip       ]
+  [C D]       = [       Iq  ] ⋅ [     S^-1] ⋅ [-C A^-1   Iq]
+
+  where S = D - C A^-1 B
+"""
+struct SchurComplementSolver{T1,T2,T3,T4} <: Gridap.Algebra.LinearSolver
+  A :: T1
+  B :: T2
+  C :: T3 
+  S :: T4
+  function SchurComplementSolver(A::Gridap.Algebra.NumericalSetup,
+                                 B::AbstractMatrix,
+                                 C::AbstractMatrix,
+                                 S::Gridap.Algebra.NumericalSetup)
+    T1 = typeof(A)
+    T2 = typeof(B)
+    T3 = typeof(C)
+    T4 = typeof(S)
+    return new{T1,T2,T3,T4}(A,B,C,S)
+  end
+end
+
+struct SchurComplementSymbolicSetup <: Gridap.Algebra.SymbolicSetup
+  solver
+end
+
+function Gridap.Algebra.symbolic_setup(s::SchurComplementSolver,A::AbstractMatrix) 
+  SchurComplementSymbolicSetup(s)
+end
+
+struct SchurComplementNumericalSetup{A,B,C} <: Gridap.Algebra.NumericalSetup
+  solver::A
+  mat   ::B
+  caches::C
+end 
+
+function get_shur_complement_caches(B::AbstractMatrix,C::AbstractMatrix)
+  du = allocate_col_vector(C)
+  bu = allocate_col_vector(C)
+  bp = allocate_col_vector(B)
+  return du,bu,bp
+end
+
+function Gridap.Algebra.numerical_setup(ss::SchurComplementSymbolicSetup,mat::AbstractMatrix)
+  s   = ss.solver
+  caches = get_shur_complement_caches(s.B,s.C)
+  return SchurComplementNumericalSetup(s,mat,caches)
+end
+
+function Gridap.Algebra.solve!(x::AbstractBlockVector,ns::SchurComplementNumericalSetup,y::AbstractBlockVector)
+  s = ns.solver
+  A,B,C,S = s.A,s.B,s.C,s.S
+  du,bu,bp = ns.caches
+
+  @check blocklength(x) == blocklength(y) == 2
+  y_u = y[Block(1)]; y_p = y[Block(2)]
+  x_u = x[Block(1)]; x_p = x[Block(2)]
+
+  # Solve Schur complement
+  solve!(x_u,A,y_u)                      # x_u = A^-1 y_u
+  copy!(bp,y_p); mul!(bp,C,du,1.0,-1.0)  # bp  = C*(A^-1 y_u) - y_p
+  solve!(x_p,S,bp)                       # x_p = S^-1 bp
+
+  mul!(bu,B,x_p)         # bu  = B*x_p
+  solve!(du,A,bu)        # du = A^-1 bu
+  x_u .-= du             # x_u  = x_u - du
+
+  return x
+end
diff --git a/src/LinearSolvers/SymGaussSeidelSmoothers.jl b/src/LinearSolvers/SymGaussSeidelSmoothers.jl
new file mode 100644
index 00000000..841eac96
--- /dev/null
+++ b/src/LinearSolvers/SymGaussSeidelSmoothers.jl
@@ -0,0 +1,187 @@
+
+# Extensions of IterativeSolvers.jl to support parallel matrices
+struct DiagonalIndices{Tv,Ti,A,B}
+  mat  :: A
+  diag :: B
+  last :: B
+  function DiagonalIndices(mat ::AbstractSparseMatrix{Tv,Ti},
+                           diag::AbstractVector{Ti},
+                           last::AbstractVector{Ti}) where {Tv,Ti}
+    A = typeof(mat)
+    B = typeof(diag)
+    @check typeof(last) == B
+    new{Tv,Ti,A,B}(mat,diag,last)
+  end
+end
+
+function DiagonalIndices(A::SparseMatrixCSR{Tv,Ti},row_range) where {Tv,Ti}
+  @notimplemented
+end
+
+function DiagonalIndices(A::SparseMatrixCSC{Tv,Ti},col_range) where {Tv,Ti}
+  n    = length(col_range)
+  diag = Vector{Ti}(undef, n)
+  last = Vector{Ti}(undef, n)
+
+  for col in col_range
+    # Diagonal index
+    r1 = Int(A.colptr[col])
+    r2 = Int(A.colptr[col + 1] - 1)
+    r1 = searchsortedfirst(A.rowval, col, r1, r2, Base.Order.Forward)
+    if r1 > r2 || A.rowval[r1] != col || iszero(A.nzval[r1])
+        throw(LinearAlgebra.SingularException(col))
+    end
+    diag[col] = r1
+
+    # Last owned index
+    r1 = Int(A.colptr[col])
+    r2 = Int(A.colptr[col + 1] - 1)
+    r1 = searchsortedfirst(A.rowval, n+1, r1, r2, Base.Order.Forward) - 1
+    last[col] = r1
+  end
+  return DiagonalIndices(A,diag,last)
+end
+
+struct LowerTriangular{Tv,Ti,A,B}
+  mat  :: A
+  diag :: DiagonalIndices{Tv,Ti,A,B}
+end
+
+struct UpperTriangular{Tv,Ti,A,B}
+  mat  :: A
+  diag :: DiagonalIndices{Tv,Ti,A,B}
+end
+
+function forward_sub!(L::LowerTriangular{Tv,Ti,<:SparseMatrixCSC},x::AbstractVector) where {Tv,Ti}
+  A, diag, last = L.mat, L.diag.diag, L.diag.last
+  n = length(diag)
+  for col = 1 : n
+    # Solve for diagonal element
+    idx = diag[col]
+    x[col] /= A.nzval[idx]
+
+    # Substitute next values involving x[col]
+    for i = idx + 1 : last[col]
+      x[A.rowval[i]] -= A.nzval[i] * x[col]
+    end
+  end
+  return x
+end
+
+function forward_sub!(L::AbstractArray{<:LowerTriangular},x::PVector)
+  map(L,own_values(x)) do L, x
+    forward_sub!(L, x)
+  end
+end
+
+function backward_sub!(U::UpperTriangular{Tv,Ti,<:SparseMatrixCSC}, x::AbstractVector) where {Tv,Ti}
+  A, diag = U.mat, U.diag.diag
+  n = length(diag)
+  for col = n : -1 : 1
+    # Solve for diagonal element
+    idx = diag[col]
+    x[col] = x[col] / A.nzval[idx]
+
+    # Substitute next values involving x[col]
+    for i = A.colptr[col] : idx - 1
+      x[A.rowval[i]] -= A.nzval[i] * x[col]
+    end
+  end
+  return x
+end
+
+function backward_sub!(U::AbstractArray{<:UpperTriangular},x::PVector)
+  map(U,own_values(x)) do U, x
+    backward_sub!(U, x)
+  end
+end
+
+# Smoother
+struct SymGaussSeidelSmoother <: Gridap.Algebra.LinearSolver
+  num_iters::Int
+end
+
+struct SymGaussSeidelSymbolicSetup <: Gridap.Algebra.SymbolicSetup
+  solver :: SymGaussSeidelSmoother
+end
+
+function Gridap.Algebra.symbolic_setup(s::SymGaussSeidelSmoother,A::AbstractMatrix)
+  SymGaussSeidelSymbolicSetup(s)
+end
+
+# Numerical setup
+
+struct SymGaussSeidelNumericalSetup{A,B,C,D} <: Gridap.Algebra.NumericalSetup
+  solver :: SymGaussSeidelSmoother
+  mat    :: A
+  L      :: B
+  U      :: C
+  caches :: D
+end
+
+function _gs_get_caches(A::AbstractMatrix)
+  dx  = allocate_col_vector(A)
+  Adx = allocate_row_vector(A)
+  return dx, Adx
+end
+
+function _gs_decompose_matrix(A::AbstractMatrix)
+  idx_range = 1:minimum(size(A))
+  D  = DiagonalIndices(A,idx_range)
+  L  = LowerTriangular(A, D)
+  U  = UpperTriangular(A, D)
+  return L,U
+end
+
+function _gs_decompose_matrix(A::PSparseMatrix{T,<:AbstractArray{MatType}}) where {T, MatType}
+  values  = partition(A)
+  indices = isa(PartitionedArrays.getany(values),SparseMatrixCSC) ? partition(axes(A,2)) : partition(axes(A,1))
+  L,U = map(values,indices) do A, indices
+    D = DiagonalIndices(A,own_to_local(indices))
+    L = LowerTriangular(A,D)
+    U = UpperTriangular(A,D)
+    return L,U
+  end |> tuple_of_arrays
+  return L,U
+end
+
+function Gridap.Algebra.numerical_setup(ss::SymGaussSeidelSymbolicSetup,A::AbstractMatrix)
+  L, U   = _gs_decompose_matrix(A)
+  caches = _gs_get_caches(A)
+  return SymGaussSeidelNumericalSetup(ss.solver,A,L,U,caches)
+end
+
+# Solve
+
+function Gridap.Algebra.solve!(x::AbstractVector, ns::SymGaussSeidelNumericalSetup, r::AbstractVector)
+  A, L, U, caches = ns.mat, ns.L, ns.U, ns.caches
+  dx, Adx = caches
+  niter   = ns.solver.num_iters
+
+  iter = 1
+  while iter <= niter
+    # Forward pass
+    copy!(dx,r)
+    forward_sub!(L, dx)
+    x  .= x .+ dx
+    mul!(Adx, A, dx)
+    r  .= r .- Adx
+
+    # Backward pass
+    copy!(dx,r)
+    backward_sub!(U, dx)
+    x  .= x .+ dx
+    mul!(Adx, A, dx)
+    r  .= r .- Adx
+
+    iter += 1
+  end
+
+  return x
+end
+
+function LinearAlgebra.ldiv!(x::AbstractVector, ns::SymGaussSeidelNumericalSetup, b::AbstractVector)
+  fill!(x,0.0)
+  aux = copy(b)
+  solve!(x,ns,aux)
+end
diff --git a/src/PatchBasedSmoothers/PatchBasedSmoothers.jl b/src/PatchBasedSmoothers/PatchBasedSmoothers.jl
index e6a03299..d4673061 100644
--- a/src/PatchBasedSmoothers/PatchBasedSmoothers.jl
+++ b/src/PatchBasedSmoothers/PatchBasedSmoothers.jl
@@ -19,6 +19,7 @@ export PatchFESpace
 export PatchBasedLinearSolver
 
 include("seq/PatchDecompositions.jl")
+include("seq/PatchTriangulations.jl")
 include("seq/PatchFESpaces.jl")
 include("seq/PatchBasedLinearSolvers.jl")
 
diff --git a/src/PatchBasedSmoothers/seq/PatchDecompositions.jl b/src/PatchBasedSmoothers/seq/PatchDecompositions.jl
index 913f0c3d..c549e76e 100644
--- a/src/PatchBasedSmoothers/seq/PatchDecompositions.jl
+++ b/src/PatchBasedSmoothers/seq/PatchDecompositions.jl
@@ -2,6 +2,8 @@ abstract type PatchBoundaryStyle end ;
 struct PatchBoundaryExclude  <: PatchBoundaryStyle end ;
 struct PatchBoundaryInclude  <: PatchBoundaryStyle end ;
 
+# TODO: Make patch_cells a Table
+
 # Question? Might a patch decomposition involve patches
 #           with roots of different topological dimension?
 #           This is not currently supported.
@@ -56,11 +58,6 @@ function PatchDecomposition(
                                    patch_cells_faces_on_boundary)
 end
 
-function Gridap.Geometry.Triangulation(a::PatchDecomposition)
-   patch_cells = Gridap.Arrays.Table(a.patch_cells)
-   return view(Triangulation(a.model),patch_cells.data)
-end
-
 function setup_patch_cells_overlapped_mesh(patch_cells)
   num_patches = length(patch_cells)
   cache = array_cache(patch_cells)
@@ -222,3 +219,194 @@ function generate_patch_boundary_faces!(patch_cells_faces_on_boundary,
     end
   end
 end
+
+# Patch cell faces: 
+#   patch_faces[pcell] = [face1, face2, ...]
+#   where face1, face2, ... are the faces of the overlapped cell `pcell` such that 
+#        - they are NOT on the boundary of the patch
+#        - they are flagged `true` in faces_mask
+function get_patch_cell_faces(PD::PatchDecomposition,Df::Integer)
+  model    = PD.model
+  topo     = get_grid_topology(model)
+  faces_mask = Fill(true,num_faces(topo,Df))
+  return get_patch_cell_faces(PD,Df,faces_mask)
+end
+
+function get_patch_cell_faces(PD::PatchDecomposition{Dc},Df::Integer,faces_mask) where Dc
+  model    = PD.model
+  topo     = get_grid_topology(model)
+
+  c2e_map  = Gridap.Geometry.get_faces(topo,Dc,Df)
+  patch_cells = Gridap.Arrays.Table(PD.patch_cells)
+  patch_cell_faces  = map(Reindex(c2e_map),patch_cells.data)
+  faces_on_boundary = PD.patch_cells_faces_on_boundary[Df+1]
+
+  patch_faces = _allocate_patch_cell_faces(patch_cell_faces,faces_on_boundary,faces_mask)
+  _generate_patch_cell_faces!(patch_faces,patch_cell_faces,faces_on_boundary,faces_mask)
+
+  return patch_faces
+end
+
+function _allocate_patch_cell_faces(patch_cell_faces,faces_on_boundary,faces_mask)
+  num_patch_cells = length(patch_cell_faces)
+
+  num_patch_faces = 0
+  patch_cells_faces_cache = array_cache(patch_cell_faces)
+  faces_on_boundary_cache = array_cache(faces_on_boundary)
+  for iC in 1:num_patch_cells
+    cell_faces  = getindex!(patch_cells_faces_cache,patch_cell_faces,iC)
+    on_boundary = getindex!(faces_on_boundary_cache,faces_on_boundary,iC)
+    for (iF,face) in enumerate(cell_faces)
+      if (!on_boundary[iF] && faces_mask[face])
+        num_patch_faces += 1
+      end
+    end
+  end
+
+  patch_faces_data = zeros(Int64,num_patch_faces)
+  patch_faces_ptrs = zeros(Int64,num_patch_cells+1)
+  return Gridap.Arrays.Table(patch_faces_data,patch_faces_ptrs)
+end
+
+function _generate_patch_cell_faces!(patch_faces,patch_cell_faces,faces_on_boundary,faces_mask)
+  num_patch_cells = length(patch_cell_faces)
+  patch_faces_data, patch_faces_ptrs = patch_faces.data, patch_faces.ptrs
+
+  pface = 1
+  patch_faces_ptrs[1] = 1
+  patch_cells_faces_cache = array_cache(patch_cell_faces)
+  faces_on_boundary_cache = array_cache(faces_on_boundary)
+  for iC in 1:num_patch_cells
+    cell_faces  = getindex!(patch_cells_faces_cache,patch_cell_faces,iC)
+    on_boundary = getindex!(faces_on_boundary_cache,faces_on_boundary,iC)
+    patch_faces_ptrs[iC+1] = patch_faces_ptrs[iC]
+    for (iF,face) in enumerate(cell_faces)
+      if (!on_boundary[iF] && faces_mask[face])
+        patch_faces_data[pface] = face
+        patch_faces_ptrs[iC+1] += 1
+        pface += 1
+      end
+    end
+  end
+
+  return patch_faces
+end
+
+# Patch faces: 
+#   patch_faces[patch] = [face1, face2, ...]
+#   where face1, face2, ... are the faces of the patch such that 
+#        - they are NOT on the boundary of the patch
+#        - they are flagged `true` in faces_mask
+function get_patch_faces(PD::PatchDecomposition{Dc},Df::Integer,faces_mask) where Dc
+  model    = PD.model
+  topo     = get_grid_topology(model)
+
+  c2e_map  = Gridap.Geometry.get_faces(topo,Dc,Df)
+  patch_cells = Gridap.Arrays.Table(PD.patch_cells)
+  patch_cell_faces  = map(Reindex(c2e_map),patch_cells.data)
+  faces_on_boundary = PD.patch_cells_faces_on_boundary[Df+1]
+
+  patch_faces = _allocate_patch_faces(patch_cells,patch_cell_faces,faces_on_boundary,faces_mask)
+  _generate_patch_faces!(patch_faces,patch_cells,patch_cell_faces,faces_on_boundary,faces_mask)
+
+  return patch_faces
+end
+
+function _allocate_patch_faces(patch_cells,patch_cell_faces,faces_on_boundary,faces_mask)
+  num_patches = length(patch_cells)
+
+  touched = Dict{Int,Bool}()
+  pcell = 1
+  num_patch_faces = 0
+  patch_cells_cache       = array_cache(patch_cells)
+  patch_cells_faces_cache = array_cache(patch_cell_faces)
+  faces_on_boundary_cache = array_cache(faces_on_boundary)
+  for patch in 1:num_patches
+    current_patch_cells = getindex!(patch_cells_cache,patch_cells,patch)
+    for iC_local in 1:length(current_patch_cells)
+      cell_faces  = getindex!(patch_cells_faces_cache,patch_cell_faces,pcell)
+      on_boundary = getindex!(faces_on_boundary_cache,faces_on_boundary,pcell)
+      for (iF,face) in enumerate(cell_faces)
+        if (!on_boundary[iF] && faces_mask[face] && !haskey(touched,face))
+          num_patch_faces += 1
+          touched[face] = true
+        end
+      end
+      pcell += 1
+    end
+    empty!(touched)
+  end
+
+  patch_faces_data = zeros(Int64,num_patch_faces)
+  patch_faces_ptrs = zeros(Int64,num_patches+1)
+  return Gridap.Arrays.Table(patch_faces_data,patch_faces_ptrs)
+end
+
+function _generate_patch_faces!(patch_faces,patch_cells,patch_cell_faces,faces_on_boundary,faces_mask)
+  num_patches = length(patch_cells)
+  patch_faces_data, patch_faces_ptrs = patch_faces.data, patch_faces.ptrs
+
+  touched = Dict{Int,Bool}()
+  pcell = 1
+  pface = 1
+  patch_faces_ptrs[1] = 1
+  patch_cells_cache       = array_cache(patch_cells)
+  patch_cells_faces_cache = array_cache(patch_cell_faces)
+  faces_on_boundary_cache = array_cache(faces_on_boundary)
+  for patch in 1:num_patches
+    current_patch_cells = getindex!(patch_cells_cache,patch_cells,patch)
+    patch_faces_ptrs[patch+1] = patch_faces_ptrs[patch]
+    for _ in 1:length(current_patch_cells)
+      cell_faces  = getindex!(patch_cells_faces_cache,patch_cell_faces,pcell)
+      on_boundary = getindex!(faces_on_boundary_cache,faces_on_boundary,pcell)
+      for (iF,face) in enumerate(cell_faces)
+        if (!on_boundary[iF] && faces_mask[face] && !haskey(touched,face))
+          patch_faces_data[pface] = face
+          patch_faces_ptrs[patch+1] += 1
+          touched[face] = true
+          pface += 1
+        end
+      end
+      pcell += 1
+    end
+    empty!(touched)
+  end
+
+  return patch_faces
+end
+
+# Face connectivity for the patches
+#    pfaces_to_pcells[pface] = [pcell1, pcell2, ...]
+# This would be the Gridap equivalent to `get_faces(patch_topology,Df,Dc)`.
+# The argument `patch_faces` allows to select only some pfaces (i.e boundary/skeleton/etc...).
+function get_pfaces_to_pcells(PD::PatchDecomposition{Dc},Df::Integer,patch_faces) where Dc
+  model    = PD.model
+  topo     = get_grid_topology(model)
+
+  faces_to_cells  = Gridap.Geometry.get_faces(topo,Df,Dc)
+  pfaces_to_cells = lazy_map(Reindex(faces_to_cells),patch_faces.data)
+  patch_cells     = Gridap.Arrays.Table(PD.patch_cells)
+  patch_cells_overlapped = PD.patch_cells_overlapped_mesh
+
+  num_patches = length(patch_cells)
+  pf2pc_ptrs  = Gridap.Adaptivity.counts_to_ptrs(map(length,pfaces_to_cells))
+  pf2pc_data  = zeros(Int64,pf2pc_ptrs[end]-1)
+
+  patch_cells_cache = array_cache(patch_cells)
+  patch_cells_overlapped_cache = array_cache(patch_cells_overlapped)
+  pfaces_to_cells_cache = array_cache(pfaces_to_cells)
+  for patch in 1:num_patches
+    cells = getindex!(patch_cells_cache,patch_cells,patch)
+    cells_overlapped = getindex!(patch_cells_overlapped_cache,patch_cells_overlapped,patch)
+    for pface in patch_faces.ptrs[patch]:patch_faces.ptrs[patch+1]-1
+      pface_to_cells = getindex!(pfaces_to_cells_cache,pfaces_to_cells,pface)
+      for (lid,cell) in enumerate(pface_to_cells)
+        lid_patch = findfirst(c->c==cell,cells)
+        pf2pc_data[pf2pc_ptrs[pface]+lid-1] = cells_overlapped[lid_patch]
+      end
+    end
+  end
+
+  pfaces_to_pcells = Gridap.Arrays.Table(pf2pc_data,pf2pc_ptrs)
+  return pfaces_to_pcells
+end
diff --git a/src/PatchBasedSmoothers/seq/PatchFESpaces.jl b/src/PatchBasedSmoothers/seq/PatchFESpaces.jl
index 9c6358f6..2605a8ab 100644
--- a/src/PatchBasedSmoothers/seq/PatchFESpaces.jl
+++ b/src/PatchBasedSmoothers/seq/PatchFESpaces.jl
@@ -77,18 +77,48 @@ end
 
 Gridap.FESpaces.get_dof_value_type(a::PatchFESpace) = Gridap.FESpaces.get_dof_value_type(a.Vh)
 Gridap.FESpaces.get_free_dof_ids(a::PatchFESpace)   = Base.OneTo(a.num_dofs)
-Gridap.FESpaces.get_cell_dof_ids(a::PatchFESpace)   = a.patch_cell_dofs_ids
-Gridap.FESpaces.get_cell_dof_ids(a::PatchFESpace,::Triangulation) = a.patch_cell_dofs_ids
 Gridap.FESpaces.get_fe_basis(a::PatchFESpace)       = get_fe_basis(a.Vh)
 Gridap.FESpaces.ConstraintStyle(::PatchFESpace)     = Gridap.FESpaces.UnConstrained()
 Gridap.FESpaces.get_vector_type(a::PatchFESpace)    = get_vector_type(a.Vh)
 Gridap.FESpaces.get_fe_dof_basis(a::PatchFESpace)   = get_fe_dof_basis(a.Vh)
 
+# get_cell_dof_ids
+
+Gridap.FESpaces.get_cell_dof_ids(a::PatchFESpace)   = a.patch_cell_dofs_ids
+Gridap.FESpaces.get_cell_dof_ids(a::PatchFESpace,::Triangulation) = @notimplemented
+
+function Gridap.FESpaces.get_cell_dof_ids(a::PatchFESpace,trian::PatchTriangulation)
+  return get_cell_dof_ids(trian.trian,a,trian)
+end
+
+function Gridap.FESpaces.get_cell_dof_ids(::Triangulation,a::PatchFESpace,trian::PatchTriangulation)
+  return a.patch_cell_dofs_ids
+end
+
+function Gridap.FESpaces.get_cell_dof_ids(::BoundaryTriangulation,a::PatchFESpace,trian::PatchTriangulation)
+  cell_dof_ids     = get_cell_dof_ids(a)
+  pfaces_to_pcells = trian.pfaces_to_pcells
+  return lazy_map(Reindex(cell_dof_ids),lazy_map(x->x[1],pfaces_to_pcells))
+end
+
+function Gridap.FESpaces.get_cell_dof_ids(::SkeletonTriangulation,a::PatchFESpace,trian::PatchTriangulation)
+  cell_dof_ids     = get_cell_dof_ids(a)
+  pfaces_to_pcells = trian.pfaces_to_pcells
+  
+  plus         = lazy_map(Reindex(cell_dof_ids),lazy_map(x->x[1],pfaces_to_pcells))
+  minus        = lazy_map(Reindex(cell_dof_ids),lazy_map(x->x[2],pfaces_to_pcells))
+  return lazy_map(Gridap.Fields.BlockMap(2,[1,2]),plus,minus)
+end
+
+# scatter dof values
+
 function Gridap.FESpaces.scatter_free_and_dirichlet_values(f::PatchFESpace,free_values,dirichlet_values)
   cell_vals = Gridap.Fields.PosNegReindex(free_values,dirichlet_values)
   return lazy_map(Broadcasting(cell_vals),f.patch_cell_dofs_ids)
 end
 
+# Construction of the patch cell dofs ids
+
 function setup_cell_reffe(model::DiscreteModel,reffe::Tuple{<:Gridap.FESpaces.ReferenceFEName,Any,Any}; kwargs...)
   basis, reffe_args,reffe_kwargs = reffe
   cell_reffe = ReferenceFE(model,basis,reffe_args...;reffe_kwargs...)
diff --git a/src/PatchBasedSmoothers/seq/PatchTriangulations.jl b/src/PatchBasedSmoothers/seq/PatchTriangulations.jl
new file mode 100644
index 00000000..b16abfd3
--- /dev/null
+++ b/src/PatchBasedSmoothers/seq/PatchTriangulations.jl
@@ -0,0 +1,99 @@
+
+struct PatchTriangulation{Dc,Dp,A,B,C,D,E} <: Gridap.Geometry.Triangulation{Dc,Dp}
+  trian            :: A
+  PD               :: B
+  patch_faces      :: C
+  pfaces_to_pcells :: D
+  mface_to_tface   :: E
+
+  function PatchTriangulation(trian::Triangulation{Dc,Dp},
+                              PD::PatchDecomposition,
+                              patch_faces,pfaces_to_pcells,mface_to_tface) where {Dc,Dp}
+    A = typeof(trian)
+    B = typeof(PD)
+    C = typeof(patch_faces)
+    D = typeof(pfaces_to_pcells)
+    E = typeof(mface_to_tface)
+    new{Dc,Dp,A,B,C,D,E}(trian,PD,patch_faces,pfaces_to_pcells,mface_to_tface)
+  end
+end
+
+# Triangulation API
+
+function Geometry.get_background_model(t::PatchTriangulation)
+  get_background_model(t.trian)
+end
+
+function Geometry.get_grid(t::PatchTriangulation)
+  get_grid(t.trian)
+end
+
+function Geometry.get_glue(t::PatchTriangulation,::Val{d}) where d
+  get_glue(t.trian,Val(d))
+end
+
+function Geometry.get_facet_normal(trian::PatchTriangulation)
+  get_facet_normal(trian.trian)
+end
+
+# Constructors 
+
+function Gridap.Geometry.Triangulation(PD::PatchDecomposition)
+  patch_cells = Gridap.Arrays.Table(PD.patch_cells)
+  trian = Triangulation(PD.model)
+  return PatchTriangulation(trian,PD,patch_cells,nothing,nothing)
+end
+
+function Gridap.Geometry.BoundaryTriangulation(PD::PatchDecomposition{Dc}) where Dc
+  Df       = Dc -1 
+  model    = PD.model
+  labeling = get_face_labeling(model)
+
+  is_boundary = get_face_mask(labeling,["boundary"],Df)
+  patch_faces = get_patch_faces(PD,Df,is_boundary)
+  pfaces_to_pcells = get_pfaces_to_pcells(PD,Df,patch_faces)
+
+  trian = BoundaryTriangulation(model)
+  glue  = get_glue(trian,Val(Df))
+  mface_to_tface = Gridap.Arrays.find_inverse_index_map(glue.tface_to_mface,num_faces(model,Df))
+  
+  return PatchTriangulation(trian,PD,patch_faces,pfaces_to_pcells,mface_to_tface)
+end
+
+function Gridap.Geometry.SkeletonTriangulation(PD::PatchDecomposition{Dc}) where Dc
+  Df       = Dc -1 
+  model    = PD.model
+  labeling = get_face_labeling(model)
+
+  is_interior = get_face_mask(labeling,["interior"],Df)
+  patch_faces = get_patch_faces(PD,Df,is_interior)
+  pfaces_to_pcells = get_pfaces_to_pcells(PD,Df,patch_faces)
+
+  trian = SkeletonTriangulation(model)
+  glue  = get_glue(trian,Val(Df))
+  mface_to_tface = Gridap.Arrays.find_inverse_index_map(glue.tface_to_mface,num_faces(model,Df))
+  
+  return PatchTriangulation(trian,PD,patch_faces,pfaces_to_pcells,mface_to_tface)
+end
+
+# Move contributions
+
+function Gridap.Geometry.move_contributions(scell_to_val::AbstractArray,strian::PatchTriangulation)
+  return move_contributions(strian.trian,scell_to_val,strian)
+end
+
+function Gridap.Geometry.move_contributions(::Triangulation,
+                                            scell_to_val::AbstractArray,
+                                            strian::PatchTriangulation)
+  patch_cells = strian.patch_faces
+  return lazy_map(Reindex(scell_to_val),patch_cells.data), strian
+end
+
+function Gridap.Geometry.move_contributions(::Union{<:BoundaryTriangulation,<:SkeletonTriangulation},
+                                            scell_to_val::AbstractArray,
+                                            strian::PatchTriangulation)
+  patch_faces      = strian.patch_faces
+  mface_to_tface   = strian.mface_to_tface
+  patch_faces_data = lazy_map(Reindex(mface_to_tface),patch_faces.data)
+  return lazy_map(Reindex(scell_to_val),patch_faces_data), strian
+end
diff --git a/test/mpi/GMGLinearSolversPoissonTests.jl b/test/mpi/GMGLinearSolversPoissonTests.jl
index 4bf6a7d3..ba90b3ef 100644
--- a/test/mpi/GMGLinearSolversPoissonTests.jl
+++ b/test/mpi/GMGLinearSolversPoissonTests.jl
@@ -37,7 +37,8 @@ function main(parts, coarse_grid_partition, num_parts_x_level, num_refs_coarse,
     smatrices, A, b = compute_hierarchy_matrices(trials,biform,liform,qdegree)
     
     # Preconditioner
-    smoothers = Fill(RichardsonSmoother(JacobiLinearSolver(),10,2.0/3.0),num_levels-1)
+    #smoothers = Fill(RichardsonSmoother(JacobiLinearSolver(),10,9.0/8.0),num_levels-1)
+    smoothers = Fill(SymGaussSeidelSmoother(5),num_levels-1)
     restrictions, prolongations = setup_transfer_operators(trials,qdegree;mode=:residual)
 
     gmg = GMGLinearSolver(mh,
diff --git a/test/mpi/RichardsonSmoothersTests.jl b/test/mpi/RichardsonSmoothersTests.jl
index d13e62e1..16074250 100644
--- a/test/mpi/RichardsonSmoothersTests.jl
+++ b/test/mpi/RichardsonSmoothersTests.jl
@@ -5,8 +5,8 @@ using MPI
 using Gridap
 using GridapDistributed
 using PartitionedArrays
-using GridapP4est
 using IterativeSolvers
+using GridapP4est
 
 using GridapSolvers
 using GridapSolvers.LinearSolvers
diff --git a/test/mpi/SymGaussSeidelSmoothersTests.jl b/test/mpi/SymGaussSeidelSmoothersTests.jl
new file mode 100644
index 00000000..86e10d1f
--- /dev/null
+++ b/test/mpi/SymGaussSeidelSmoothersTests.jl
@@ -0,0 +1,66 @@
+module RichardsonSmoothersTests
+
+using Test
+using MPI
+using Gridap
+using GridapDistributed
+using PartitionedArrays
+using IterativeSolvers
+
+using GridapSolvers
+using GridapSolvers.LinearSolvers
+
+function main(parts,num_ranks,mesh_partition)
+  domain = (0,1,0,1)
+  model  = CartesianDiscreteModel(parts,num_ranks,domain,mesh_partition)
+
+  sol(x) = x[1] + x[2]
+  f(x)   = -Δ(sol)(x)
+
+  order  = 1
+  qorder = order*2 + 1
+  reffe  = ReferenceFE(lagrangian,Float64,order)
+  Vh     = TestFESpace(model,reffe;conformity=:H1,dirichlet_tags="boundary")
+  Uh     = TrialFESpace(Vh,sol)
+  u      = interpolate(sol,Uh)
+
+  Ω      = Triangulation(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  = SymGaussSeidelSmoother(10)
+  ss = symbolic_setup(P,A)
+  ns = numerical_setup(ss,A)
+
+  x = pfill(1.0,partition(axes(A,2)))
+  x, history = IterativeSolvers.cg!(x,A,b;
+                                    verbose=i_am_main(parts),
+                                    reltol=1.0e-8,
+                                    Pl=ns,
+                                    log=true)
+
+  u  = interpolate(sol,Uh)
+  uh = FEFunction(Uh,x)
+  eh = uh - u
+  E  = sum(∫(eh*eh)*dΩ)
+  if i_am_main(parts)
+    println("L2 Error: ", E)
+  end
+  
+  @test E < 1.e-8
+end
+
+mesh_partition = (32,32)
+num_ranks = (2,2)
+parts = with_mpi() do distribute
+  distribute(LinearIndices((prod(num_ranks),)))
+end
+
+main(parts,num_ranks,mesh_partition)
+MPI.Finalize()
+
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index ae14d3e1..5a282772 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -44,7 +44,7 @@ function run_tests(testdir)
           np = 4
           extra_args = ""
         else
-          np = nprocs
+          np = 4 # nprocs
           extra_args = ""
         end
         if ! image_file_exists
@@ -63,4 +63,11 @@ end
 run_tests(joinpath(@__DIR__, "mpi"))
 
 # Sequential tests
+@time @testset "BlockDiagonalSmoothersPETScTests" begin include("seq/BlockDiagonalSmoothersPETScTests.jl") end
+@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 "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
+@time @testset "SchurComplementSolversTests" begin include("seq/SchurComplementSolversTests.jl") end
diff --git a/test/seq/BlockDiagonalSmoothersTests.jl b/test/seq/BlockDiagonalSmoothersTests.jl
new file mode 100644
index 00000000..a8a98869
--- /dev/null
+++ b/test/seq/BlockDiagonalSmoothersTests.jl
@@ -0,0 +1,149 @@
+module BlockDiagonalSmoothersTests
+
+using Test
+using Gridap
+using Gridap.MultiField
+using BlockArrays
+using LinearAlgebra
+using FillArrays
+using IterativeSolvers
+using PartitionedArrays
+
+using GridapDistributed
+using GridapSolvers
+using GridapPETSc
+
+using GridapDistributed: BlockPVector, BlockPMatrix
+
+u(x) = VectorValue(x[1],x[2])
+f(x) = VectorValue(2.0*x[2]*(1.0-x[1]*x[1]),2.0*x[1]*(1-x[2]*x[2]))
+
+p(x) = x[1] + x[2]
+g(x) = -Δ(p)(x)
+
+function set_ksp_options(ksp)
+  pc       = Ref{GridapPETSc.PETSC.PC}()
+  mumpsmat = Ref{GridapPETSc.PETSC.Mat}()
+  @check_error_code GridapPETSc.PETSC.KSPView(ksp[],C_NULL)
+  @check_error_code GridapPETSc.PETSC.KSPSetType(ksp[],GridapPETSc.PETSC.KSPPREONLY)
+  @check_error_code GridapPETSc.PETSC.KSPGetPC(ksp[],pc)
+  @check_error_code GridapPETSc.PETSC.PCSetType(pc[],GridapPETSc.PETSC.PCLU)
+  @check_error_code GridapPETSc.PETSC.PCFactorSetMatSolverType(pc[],GridapPETSc.PETSC.MATSOLVERMUMPS)
+  @check_error_code GridapPETSc.PETSC.PCFactorSetUpMatSolverType(pc[])
+  @check_error_code GridapPETSc.PETSC.PCFactorGetMatrix(pc[],mumpsmat)
+  @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[],  4, 1)
+  # percentage increase in the estimated working space
+  @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[],  14, 1000)
+  @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 28, 2)
+  @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 29, 2)
+  @check_error_code GridapPETSc.PETSC.MatMumpsSetCntl(mumpsmat[], 3, 1.0e-6)
+end
+
+function _is_same_vector(x1,x2,X1,X2)
+  res = true
+  for i in 1:length(X1)
+    x1i = restrict_to_field(X1,x1,i)
+    x2i = restrict_to_field(X2,x2,i)
+    res = res & (norm(x1i-x2i) < 1.e-5)
+  end
+  return res
+end
+
+function is_same_vector(x1::BlockVector,x2,X1,X2)
+  _is_same_vector(x1,x2,X1,X2)
+end
+
+function is_same_vector(x1::BlockPVector,x2,X1,X2)
+  _x1 = GridapDistributed.change_ghost(x1,X1.gids;make_consistent=true)
+  _x2 = GridapDistributed.change_ghost(x2,X2.gids;make_consistent=true)
+  _is_same_vector(_x1,_x2,X1,X2)
+end
+
+function main(model,use_petsc::Bool)
+  if use_petsc
+    GridapPETSc.with() do
+      solvers = Fill(PETScLinearSolver(set_ksp_options),2)
+      main(model,solvers)
+    end
+  else
+    solvers = Fill(BackslashSolver(),2)
+    main(model,solvers)
+  end
+end
+
+function main(model,solvers)
+  order  = 2
+  reffeᵤ = ReferenceFE(lagrangian,VectorValue{D,Float64},order)
+  V = TestFESpace(model,reffeᵤ,conformity=:H1,dirichlet_tags=["boundary"])
+
+  reffeₚ = ReferenceFE(lagrangian,Float64,order)
+  Q = TestFESpace(model,reffeₚ,conformity=:H1,dirichlet_tags=["boundary"])
+
+  U = TrialFESpace(V,u)
+  P = TrialFESpace(Q,p)
+
+  Y = MultiFieldFESpace([V,Q])
+  X = MultiFieldFESpace([U,P])
+
+  mfs = BlockMultiFieldStyle()
+  Yb = MultiFieldFESpace([V,Q];style=mfs)
+  Xb = MultiFieldFESpace([U,P];style=mfs)
+
+  degree = 2*(order + 1)
+  Ω  = Triangulation(model)
+  dΩ = Measure(Ω,degree)
+
+  # Global problem
+  a((u,p),(v,q)) = ∫( v⊙u + ∇(v)⊙∇(u) + q⋅p + ∇(q)⊙∇(p))dΩ
+  l((v,q)) = ∫( v⋅f + q⋅g)dΩ
+
+  op = AffineFEOperator(a,l,X,Y)
+  x_star = get_free_dof_values(solve(op))
+
+  opb = AffineFEOperator(a,l,Xb,Yb)
+  A,b = get_matrix(opb), get_vector(opb);
+
+  # Build using local weakforms
+  a1(u,v) = ∫(v⊙u + ∇(v)⊙∇(u))dΩ
+  a2(p,q) = ∫(q⋅p + ∇(q)⊙∇(p))dΩ
+  biforms = [a1,a2]
+
+  BDS   = BlockDiagonalSmoother(biforms,Xb,Yb,solvers)
+  BDSss = symbolic_setup(BDS,A)
+  BDSns = numerical_setup(BDSss,A)
+
+  x = GridapSolvers.LinearSolvers.allocate_col_vector(A)
+  x = cg!(x,A,b;verbose=true,Pl=BDSns,reltol=1.0e-12)
+  @test is_same_vector(x,x_star,Xb,X)
+
+  # Build using BlockMatrixAssemblers
+  BDS   = BlockDiagonalSmoother(A,solvers)
+  BDSss = symbolic_setup(BDS,A)
+  BDSns = numerical_setup(BDSss,A)
+
+  x = GridapSolvers.LinearSolvers.allocate_col_vector(A)
+  x = cg!(x,A,b;verbose=true,Pl=BDSns,reltol=1.0e-12)
+  @test is_same_vector(x,x_star,Xb,X)
+end
+
+num_ranks = (2,2)
+parts = with_debug() do distribute
+  distribute(LinearIndices((prod(num_ranks),)))
+end
+
+D = 2
+n = 10
+domain = Tuple(repeat([0,1],D))
+mesh_partition = (n,n)
+
+# Serial
+model = CartesianDiscreteModel(domain,mesh_partition)
+main(model,false)
+main(model,true)
+
+# Distributed, sequential
+model = CartesianDiscreteModel(parts,num_ranks,domain,mesh_partition)
+main(model,false)
+main(model,true)
+
+end
\ No newline at end of file
diff --git a/test/seq/DistributedPatchFESpacesDebuggingTests.jl b/test/seq/DistributedPatchFESpacesDebuggingTests.jl
index bea13993..a88f390b 100644
--- a/test/seq/DistributedPatchFESpacesDebuggingTests.jl
+++ b/test/seq/DistributedPatchFESpacesDebuggingTests.jl
@@ -207,7 +207,7 @@ for key in keys(res_single)
   val_m = res_multi[key]
 
   println(">>> ", key)
-  map(val_s.values) do v
+  map(partition(val_s)) do v
     println(transpose(v))
   end;
   map(own_values(val_m)) do v
diff --git a/test/seq/GMRESSolversTests.jl b/test/seq/GMRESSolversTests.jl
new file mode 100644
index 00000000..146d5e17
--- /dev/null
+++ b/test/seq/GMRESSolversTests.jl
@@ -0,0 +1,60 @@
+module GMRESSolversTests
+
+using Test
+using Gridap
+using GridapDistributed
+using PartitionedArrays
+using IterativeSolvers
+
+using GridapSolvers
+using GridapSolvers.LinearSolvers
+
+sol(x) = x[1] + x[2]
+f(x)   = -Δ(sol)(x)
+
+function main(model)
+  order  = 1
+  qorder = order*2 + 1
+  reffe  = ReferenceFE(lagrangian,Float64,order)
+  Vh     = TestFESpace(model,reffe;conformity=:H1,dirichlet_tags="boundary")
+  Uh     = TrialFESpace(Vh,sol)
+  u      = interpolate(sol,Uh)
+
+  Ω      = Triangulation(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);
+
+  Pl = JacobiLinearSolver()
+  solver = LinearSolvers.GMRESSolver(40,Pl,1.e-8)
+  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Ω)
+  return E < 1.e-8
+end
+
+# Completely serial
+mesh_partition = (10,10)
+domain = (0,1,0,1)
+model  = CartesianDiscreteModel(domain,mesh_partition)
+@test main(model)
+
+# Sequential
+num_ranks = (1,2)
+parts = with_debug() do distribute
+  distribute(LinearIndices((prod(num_ranks),)))
+end
+
+model  = CartesianDiscreteModel(parts,num_ranks,domain,mesh_partition)
+@test main(model)
+
+end
\ No newline at end of file
diff --git a/test/seq/IterativeSolversTests.jl b/test/seq/IterativeSolversTests.jl
new file mode 100644
index 00000000..6177429f
--- /dev/null
+++ b/test/seq/IterativeSolversTests.jl
@@ -0,0 +1,97 @@
+module IterativeSolversTests
+
+using Test
+using Gridap
+using IterativeSolvers
+using LinearAlgebra
+using SparseArrays
+using PartitionedArrays
+
+using GridapSolvers
+using GridapSolvers.LinearSolvers
+
+sol(x) = x[1] + x[2]
+f(x)   = -Δ(sol)(x)
+
+function l2_error(x,Uh,dΩ)
+  u  = interpolate(sol,Uh)
+  uh = FEFunction(Uh,x)
+  eh = uh - u
+  return sum(∫(eh*eh)*dΩ)
+end
+
+function main(model,is_distributed)
+  order  = 1
+  qorder = order*2 + 1
+  reffe  = ReferenceFE(lagrangian,Float64,order)
+  Vh     = TestFESpace(model,reffe;conformity=:H1,dirichlet_tags="boundary")
+  Uh     = TrialFESpace(Vh,sol)
+  u      = interpolate(sol,Uh)
+
+  Ω      = Triangulation(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);
+
+  # CG
+  solver = IS_ConjugateGradientSolver(;maxiter=100,reltol=1.e-12)
+  ss = symbolic_setup(solver,A)
+  ns = numerical_setup(ss,A)
+
+  x = LinearSolvers.allocate_col_vector(A)
+  y = copy(b)
+  solve!(x,ns,y)
+  @test l2_error(x,Uh,dΩ) < 1.e-8
+
+  # SSOR
+  solver = IS_SSORSolver(2.0/3.0;maxiter=100)
+  ss = symbolic_setup(solver,A)
+  ns = numerical_setup(ss,A)
+
+  x = LinearSolvers.allocate_row_vector(A)
+  y = copy(b)
+  solve!(x,ns,y)
+  !is_distributed && (@test l2_error(x,Uh,dΩ) < 1.e-8)
+
+    if !is_distributed
+    # GMRES
+    solver = IS_GMRESSolver(;maxiter=100,reltol=1.e-12)
+    ss = symbolic_setup(solver,A)
+    ns = numerical_setup(ss,A)
+
+    x = LinearSolvers.allocate_row_vector(A)
+    y = copy(b)
+    solve!(x,ns,y)
+    @test l2_error(x,Uh,dΩ) < 1.e-8
+
+    # MINRES
+    solver = IS_MINRESSolver(;maxiter=100,reltol=1.e-12)
+    ss = symbolic_setup(solver,A)
+    ns = numerical_setup(ss,A)
+
+    x = LinearSolvers.allocate_row_vector(A)
+    y = copy(b)
+    solve!(x,ns,y)
+    @test l2_error(x,Uh,dΩ) < 1.e-8
+  end
+end
+
+# Completely serial
+mesh_partition = (8,8)
+domain = (0,1,0,1)
+model  = CartesianDiscreteModel(domain,mesh_partition)
+main(model,false)
+
+# Sequential
+num_ranks = (1,2)
+parts = with_debug() do distribute
+  distribute(LinearIndices((prod(num_ranks),)))
+end
+
+model  = CartesianDiscreteModel(parts,num_ranks,domain,mesh_partition)
+main(model,true)
+
+end
\ No newline at end of file
diff --git a/test/seq/PatchBasedTesting.jl b/test/seq/PatchBasedTesting.jl
new file mode 100644
index 00000000..d4653a3a
--- /dev/null
+++ b/test/seq/PatchBasedTesting.jl
@@ -0,0 +1,70 @@
+
+using LinearAlgebra
+using Test
+using PartitionedArrays
+using Gridap
+using Gridap.Arrays
+using Gridap.Helpers
+using Gridap.Geometry
+using Gridap.ReferenceFEs
+using Gridap.FESpaces
+using GridapDistributed
+using FillArrays
+
+using GridapSolvers
+import GridapSolvers.PatchBasedSmoothers as PBS
+
+num_ranks = (1,2)
+parts = with_debug() do distribute
+  distribute(LinearIndices((prod(num_ranks),)))
+end
+
+domain = (0.0,1.0,0.0,1.0)
+mesh_partition = (2,4)
+model = CartesianDiscreteModel(domain,mesh_partition)
+
+order = 1; reffe = ReferenceFE(lagrangian,Float64,order;space=:P); conformity = L2Conformity();
+#order = 1; reffe = ReferenceFE(lagrangian,Float64,order); conformity = H1Conformity();
+#order = 0; reffe = ReferenceFE(raviart_thomas,Float64,order); conformity = HDivConformity();
+Vh = TestFESpace(model,reffe,conformity=conformity)
+PD = PBS.PatchDecomposition(model)
+Ph = PBS.PatchFESpace(model,reffe,conformity,PD,Vh)
+
+# ---- Assemble systems ---- #
+
+Ω  = Triangulation(model)
+dΩ = Measure(Ω,2*order+1)
+Λ  = Skeleton(model)
+dΛ = Measure(Λ,3)
+Γ  = Boundary(model)
+dΓ = Measure(Γ,3)
+
+aΩ(u,v) = ∫(v⋅u)*dΩ
+aΛ(u,v) = ∫(jump(v)⋅jump(u))*dΛ
+aΓ(u,v) = ∫(v⋅u)*dΓ
+a(u,v)  = aΩ(u,v) + aΛ(u,v) + aΓ(u,v)
+l(v)    = ∫(1*v)*dΩ
+
+assembler = SparseMatrixAssembler(Vh,Vh)
+Ah = assemble_matrix(a,assembler,Vh,Vh)
+fh = assemble_vector(l,assembler,Vh)
+
+sol_h = solve(LUSolver(),Ah,fh)
+
+Ωₚ  = Triangulation(PD)
+dΩₚ = Measure(Ωₚ,2*order+1)
+Λₚ  = SkeletonTriangulation(PD)
+dΛₚ = Measure(Λₚ,3)
+Γₚ  = BoundaryTriangulation(PD)
+dΓₚ = Measure(Γₚ,3)
+
+aΩp(u,v) = ∫(v⋅u)*dΩₚ
+aΛp(u,v) = ∫(jump(v)⋅jump(u))*dΛₚ
+aΓp(u,v) = ∫(v⋅u)*dΓₚ
+ap(u,v)  = aΩp(u,v) + aΛp(u,v) + aΓp(u,v)
+lp(v)    = ∫(1*v)*dΩₚ
+
+assembler_P = SparseMatrixAssembler(Ph,Ph)
+Ahp = assemble_matrix(ap,assembler_P,Ph,Ph)
+fhp = assemble_vector(lp,assembler_P,Ph)
+
diff --git a/test/seq/SchurComplementSolversTests.jl b/test/seq/SchurComplementSolversTests.jl
new file mode 100644
index 00000000..123ad7bf
--- /dev/null
+++ b/test/seq/SchurComplementSolversTests.jl
@@ -0,0 +1,120 @@
+module SchurComplementSolversTests
+
+using Test
+using BlockArrays
+using Gridap
+using Gridap.MultiField
+using Gridap.Algebra
+using Gridap.Geometry
+using Gridap.FESpaces
+using Gridap.ReferenceFEs
+
+using PartitionedArrays
+using GridapDistributed
+
+using GridapSolvers
+using GridapSolvers.LinearSolvers
+using GridapSolvers.MultilevelTools
+
+function l2_error(xh,sol,dΩ)
+  eh = xh - sol
+  e  = sum(∫(eh⋅eh)dΩ)
+  return e
+end
+
+function l2_error(x,sol,X,dΩ)
+  xh = FEFunction(X,x)
+  return l2_error(xh,sol,dΩ)
+end
+
+# Darcy solution
+const β_U = 50.0
+const γ = 100.0
+
+u_ref(x) = VectorValue(x[1]+x[2],-x[2])
+p_ref(x) = 2.0*x[1]-1.0
+f_ref(x) = u_ref(x) + ∇(p_ref)(x)
+
+function main(model)
+
+  labels = get_face_labeling(model)
+  add_tag_from_tags!(labels,"dirichlet",[1,2,3,4,5,6,7])
+  add_tag_from_tags!(labels,"newmann",[8,])
+
+  order  = 0
+  reffeᵤ = ReferenceFE(raviart_thomas,Float64,order)
+  V  = TestFESpace(model,reffeᵤ,conformity=:HDiv,dirichlet_tags="dirichlet")
+  U  = TrialFESpace(V,u_ref)
+
+  reffeₚ = ReferenceFE(lagrangian,Float64,order;space=:P)
+  Q = TestFESpace(model,reffeₚ,conformity=:L2)
+  P = TrialFESpace(Q,p_ref)
+
+  mfs = BlockMultiFieldStyle()
+  Y = MultiFieldFESpace([V, Q];style=mfs)
+  X = MultiFieldFESpace([U, P];style=mfs)
+
+  qdegree = 4
+  Ω   = Triangulation(model)
+  dΩ  = Measure(Ω,qdegree)
+
+  Γ_N  = BoundaryTriangulation(model;tags="newmann")
+  dΓ_N = Measure(Γ_N,qdegree)
+  n_Γ_N = get_normal_vector(Γ_N)
+
+  a(u,v) = ∫(v⊙u)dΩ + ∫(γ*(∇⋅v)*(∇⋅u))dΩ
+  b(p,v) = ∫(-(∇⋅v)*p)dΩ
+  c(u,q) = ∫(- q*(∇⋅u))dΩ
+
+  biform((u,p),(v,q)) = a(u,v) + b(p,v) + c(u,q)
+  liform((v,q)) = ∫(f_ref⋅v)dΩ - ∫((v⋅n_Γ_N)⋅p_ref)dΓ_N
+
+  op = AffineFEOperator(biform,liform,X,Y)
+  sysmat, sysvec = get_matrix(op), get_vector(op);
+
+  ############################################################################################
+  # Solve by GMRES preconditioned with inexact Schur complement
+
+  s(p,q) = ∫(γ*p*q)dΩ
+  PS = assemble_matrix(s,P,Q)
+  PS_solver = BackslashSolver()
+  PS_ns = numerical_setup(symbolic_setup(PS_solver,PS),PS)
+
+  A = sysmat[Block(1,1)]
+  A_solver = BackslashSolver()
+  A_ns = numerical_setup(symbolic_setup(A_solver,A),A)
+
+  B = sysmat[Block(1,2)]; C = sysmat[Block(2,1)]
+  psc_solver = SchurComplementSolver(A_ns,B,C,PS_ns);
+
+  gmres = GMRESSolver(20,psc_solver,1e-10)
+  gmres_ns = numerical_setup(symbolic_setup(gmres,sysmat),sysmat)
+
+  x = LinearSolvers.allocate_col_vector(sysmat)
+  solve!(x,gmres_ns,sysvec)
+
+  xh = FEFunction(X,x)
+  uh, ph = xh
+  @test l2_error(uh,u_ref,dΩ) < 1.e-4
+  @test l2_error(ph,p_ref,dΩ) < 1.e-4
+end
+
+num_ranks = (2,2)
+parts = with_debug() do distribute
+  distribute(LinearIndices((prod(num_ranks),)))
+end
+
+D = 2
+n = 60
+domain    = Tuple(repeat([0,1],D))
+mesh_partition = (n,n)
+
+# Serial
+model     = CartesianDiscreteModel(domain,mesh_partition)
+main(model)
+
+# Distributed, sequential
+model     = CartesianDiscreteModel(parts,num_ranks,domain,mesh_partition)
+main(model)
+
+end
\ No newline at end of file
diff --git a/test/seq/SymGaussSeidelSmoothersTests.jl b/test/seq/SymGaussSeidelSmoothersTests.jl
new file mode 100644
index 00000000..e29a7864
--- /dev/null
+++ b/test/seq/SymGaussSeidelSmoothersTests.jl
@@ -0,0 +1,65 @@
+module SymGaussSeidelSmoothersTests
+
+using Test
+using MPI
+using Gridap
+using GridapDistributed
+using PartitionedArrays
+using IterativeSolvers
+
+using GridapSolvers
+using GridapSolvers.LinearSolvers
+
+sol(x) = x[1] + x[2]
+f(x)   = -Δ(sol)(x)
+
+function main(model)
+  order  = 1
+  qorder = order*2 + 1
+  reffe  = ReferenceFE(lagrangian,Float64,order)
+  Vh     = TestFESpace(model,reffe;conformity=:H1,dirichlet_tags="boundary")
+  Uh     = TrialFESpace(Vh,sol)
+  u      = interpolate(sol,Uh)
+
+  Ω      = Triangulation(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  = SymGaussSeidelSmoother(10)
+  ss = symbolic_setup(P,A)
+  ns = numerical_setup(ss,A)
+
+  x = LinearSolvers.allocate_col_vector(A)
+  x, history = IterativeSolvers.cg!(x,A,b;
+                                    verbose=true,
+                                    reltol=1.0e-8,
+                                    Pl=ns,
+                                    log=true);
+
+  u  = interpolate(sol,Uh)
+  uh = FEFunction(Uh,x)
+  eh = uh - u
+  E  = sum(∫(eh*eh)*dΩ)
+  return E < 1.e-8
+end
+
+# Completely serial
+mesh_partition = (8,8)
+domain = (0,1,0,1)
+model  = CartesianDiscreteModel(domain,mesh_partition)
+@test main(model)
+
+# Sequential
+num_ranks = (1,2)
+parts = with_debug() do distribute
+  distribute(LinearIndices((prod(num_ranks),)))
+end
+
+model  = CartesianDiscreteModel(parts,num_ranks,domain,mesh_partition)
+@test main(model)
+
+end
\ No newline at end of file