diff --git a/src/BlockSolvers/BlockDiagonalSolvers.jl b/src/BlockSolvers/BlockDiagonalSolvers.jl
new file mode 100644
index 00000000..30787edb
--- /dev/null
+++ b/src/BlockSolvers/BlockDiagonalSolvers.jl
@@ -0,0 +1,116 @@
+
+struct BlockDiagonalSolver{N,A,B} <: Gridap.Algebra.LinearSolver
+  blocks  :: A
+  solvers :: B
+  function BlockDiagonalSolver(
+    blocks  :: AbstractVector{<:SolverBlock},
+    solvers :: AbstractVector{<:Gridap.Algebra.LinearSolver}
+    )
+    N = length(solvers)
+    @check length(blocks) == N
+
+    A = typeof(blocks)
+    B = typeof(solvers)
+    return new{N,A,B}(blocks,solvers)
+  end
+end
+
+# Constructors
+
+function BlockDiagonalSolver(solvers::AbstractVector{<:Gridap.Algebra.LinearSolver}; 
+                             is_nonlinear::Vector{Bool}=fill(false,length(solvers)))
+  blocks = map(nl -> nl ? NonlinearSystemBlock() : LinearSystemBlock(),is_nonlinear)
+  return BlockDiagonalSolver(blocks,solvers)
+end
+
+function BlockDiagonalSolver(funcs   :: AbstractArray{<:Function},
+                             trials  :: AbstractArray{<:FESpace},
+                             tests   :: AbstractArray{<:FESpace},
+                             solvers :: AbstractArray{<:Gridap.Algebra.LinearSolver};
+                             is_nonlinear::Vector{Bool}=fill(false,length(solvers)))
+  blocks = map(funcs,trials,tests,is_nonlinear) do f,trial,test,nl
+    nl ? TriformBlock(f,trial,test) : BiformBlock(f,trial,test)
+  end
+  return BlockDiagonalSolver(blocks,solvers)
+end
+
+function BlockDiagonalSolver(mats::AbstractVector{<:AbstractMatrix},
+                             solvers::AbstractVector{<:Gridap.Algebra.LinearSolver})
+  blocks = map(MatrixBlock,mats)
+  return BlockDiagonalSolver(blocks,solvers)
+end
+
+# Symbolic setup
+
+struct BlockDiagonalSolverSS{A,B,C} <: Gridap.Algebra.SymbolicSetup
+  solver       :: A
+  block_ss     :: B
+  block_caches :: C
+end
+
+function Gridap.Algebra.symbolic_setup(solver::BlockDiagonalSolver,mat::AbstractBlockMatrix)
+  mat_blocks   = diag(blocks(mat))
+  block_caches = map(instantiate_block_cache,solver.blocks,mat_blocks)
+  block_ss     = map(symbolic_setup,solver.solvers,block_caches)
+  return BlockDiagonalSolverSS(solver,block_ss,block_caches)
+end
+
+function Gridap.Algebra.symbolic_setup(solver::BlockDiagonalSolver,mat::AbstractBlockMatrix,x::AbstractBlockVector)
+  mat_blocks   = diag(blocks(mat))
+  vec_blocks   = blocks(x)
+  block_caches = map(instantiate_block_cache,solver.blocks,mat_blocks,vec_blocks)
+  block_ss     = map(symbolic_setup,solver.solvers,block_caches,vec_blocks)
+  return BlockDiagonalSolverSS(solver,block_ss,block_caches)
+end
+
+# Numerical setup
+
+struct BlockDiagonalSolverNS{A,B,C} <: Gridap.Algebra.NumericalSetup
+  solver       :: A
+  block_ns     :: B
+  block_caches :: C
+end
+
+function Gridap.Algebra.numerical_setup(ss::BlockDiagonalSolverSS,mat::AbstractBlockMatrix)
+  solver     = ss.solver
+  block_ns   = map(numerical_setup,ss.block_ss,ss.block_caches)
+  return BlockDiagonalSolverNS(solver,block_ns,ss.block_caches)
+end
+
+function Gridap.Algebra.numerical_setup(ss::BlockDiagonalSolverSS,mat::AbstractBlockMatrix,x::AbstractBlockVector)
+  solver     = ss.solver
+  vec_blocks = blocks(x)
+  block_ns   = map(numerical_setup,ss.block_ss,ss.block_caches,vec_blocks)
+  return BlockDiagonalSolverNS(solver,block_ns,ss.block_caches)
+end
+
+function Gridap.Algebra.numerical_setup!(ns::BlockDiagonalSolverNS,mat::AbstractBlockMatrix)
+  solver       = ns.solver
+  mat_blocks   = diag(blocks(mat))
+  block_caches = map(update_block_cache!,ns.block_caches,solver.blocks,mat_blocks)
+  map(numerical_setup!,ns.block_ns,block_caches)
+  return ns
+end
+
+function Gridap.Algebra.numerical_setup!(ns::BlockDiagonalSolverNS,mat::AbstractBlockMatrix,x::AbstractBlockVector)
+  solver       = ns.solver
+  mat_blocks   = diag(blocks(mat))
+  vec_blocks   = blocks(x)
+  block_caches = map(update_block_cache!,ns.block_caches,solver.blocks,mat_blocks,vec_blocks)
+  map(numerical_setup!,ns.block_ns,block_caches,vec_blocks)
+  return ns
+end
+
+function Gridap.Algebra.solve!(x::AbstractBlockVector,ns::BlockDiagonalSolverNS,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::BlockDiagonalSolverNS,b)
+  solve!(x,ns,b)
+end
diff --git a/src/BlockSolvers/BlockFEOperators.jl b/src/BlockSolvers/BlockFEOperators.jl
new file mode 100644
index 00000000..9904586e
--- /dev/null
+++ b/src/BlockSolvers/BlockFEOperators.jl
@@ -0,0 +1,142 @@
+
+struct BlockFEOperator{NB,SB,P} <: FEOperator
+  global_op    :: FEOperator
+  block_ops    :: Matrix{<:Union{<:FEOperator,Missing,Nothing}}
+  is_nonlinear :: Matrix{Bool}
+end
+
+const BlockFESpaceTypes{NB,SB,P} = Union{<:MultiFieldFESpace{<:BlockMultiFieldStyle{NB,SB,P}},<:GridapDistributed.DistributedMultiFieldFESpace{<:BlockMultiFieldStyle{NB,SB,P}}}
+
+function BlockFEOperator(
+  res::Matrix{<:Union{<:Function,Missing,Nothing}},
+  jac::Matrix{<:Union{<:Function,Missing,Nothing}},
+  trial::BlockFESpaceTypes,
+  test::BlockFESpaceTypes;
+  kwargs...
+)
+  assem = SparseMatrixAssembler(test,trial)
+  return BlockFEOperator(res,jac,trial,test,assem)
+end
+
+function BlockFEOperator(
+  res::Matrix{<:Union{<:Function,Missing,Nothing}},
+  jac::Matrix{<:Union{<:Function,Missing,Nothing}},
+  trial::BlockFESpaceTypes{NB,SB,P},
+  test::BlockFESpaceTypes{NB,SB,P},
+  assem::MultiField.BlockSparseMatrixAssembler{NB,NV,SB,P};
+  is_nonlinear::Matrix{Bool}=fill(true,(NB,NB))
+) where {NB,NV,SB,P}
+  @check size(res,1) == size(jac,1) == NB
+  @check size(res,2) == size(jac,2) == NB
+
+  global_res = residual_from_blocks(NB,SB,P,res)
+  global_jac = jacobian_from_blocks(NB,SB,P,jac)
+  global_op  = FEOperator(global_res,global_jac,trial,test,assem)
+
+  trial_blocks = blocks(trial)
+  test_blocks  = blocks(test)
+  assem_blocks = blocks(assem)
+  block_ops = map(CartesianIndices(res)) do I
+    if !ismissing(res[I]) && !isnothing(res[I])
+      FEOperator(res[I],jac[I],test_blocks[I[1]],trial_blocks[I[2]],assem_blocks[I])
+    end
+  end
+  return BlockFEOperator{NB,SB,P}(global_op,block_ops,is_nonlinear)
+end
+
+# BlockArrays API
+
+BlockArrays.blocks(op::BlockFEOperator) = op.block_ops
+
+# FEOperator API
+
+FESpaces.get_test(op::BlockFEOperator) = get_test(op.global_op)
+FESpaces.get_trial(op::BlockFEOperator) = get_trial(op.global_op)
+Algebra.allocate_residual(op::BlockFEOperator,u) = allocate_residual(op.global_op,u)
+Algebra.residual(op::BlockFEOperator,u) = residual(op.global_op,u)
+Algebra.allocate_jacobian(op::BlockFEOperator,u) = allocate_jacobian(op.global_op,u)
+Algebra.jacobian(op::BlockFEOperator,u) = jacobian(op.global_op,u)
+Algebra.residual!(b::AbstractVector,op::BlockFEOperator,u) = residual!(b,op.global_op,u)
+
+function Algebra.jacobian!(A::AbstractBlockMatrix,op::BlockFEOperator{NB},u) where NB
+  map(blocks(A),blocks(op),op.is_nonlinear) do A,op,nl
+    if nl
+      residual!(A,op,u)
+    end
+  end
+  return A
+end
+
+# Private methods
+
+function residual_from_blocks(NB,SB,P,block_residuals)
+  function res(u,v)
+    block_ranges = MultiField.get_block_ranges(NB,SB,P)
+    block_u = map(r -> (length(r) == 1) ? u[r[1]] : Tuple(u[r]), block_ranges)
+    block_v = map(r -> (length(r) == 1) ? v[r[1]] : Tuple(v[r]), block_ranges)
+    block_contrs = map(CartesianIndices(block_residuals)) do I
+      if !ismissing(block_residuals[I]) && !isnothing(block_residuals[I])
+        block_residuals[I](block_u[I[2]],block_v[I[1]])
+      end
+    end
+    return add_block_contribs(block_contrs)
+  end
+  return res
+end
+
+function jacobian_from_blocks(NB,SB,P,block_jacobians)
+  function jac(u,du,dv)
+    block_ranges = MultiField.get_block_ranges(NB,SB,P)
+    block_u  = map(r -> (length(r) == 1) ?  u[r[1]] : Tuple(u[r]) , block_ranges)
+    block_du = map(r -> (length(r) == 1) ? du[r[1]] : Tuple(du[r]), block_ranges)
+    block_dv = map(r -> (length(r) == 1) ? dv[r[1]] : Tuple(dv[r]), block_ranges)
+    block_contrs = map(CartesianIndices(block_jacobians)) do I
+      if !ismissing(block_jacobians[I]) && !isnothing(block_jacobians[I])
+        block_jacobians[I](block_u[I[2]],block_du[I[2]],block_dv[I[1]])
+      end
+    end
+    return add_block_contribs(block_contrs)
+  end
+  return jac
+end
+
+function add_block_contribs(contrs)
+  c = contrs[1]
+  for ci in contrs[2:end]
+    if !ismissing(ci) && !isnothing(ci)
+      c = c + ci
+    end
+  end
+  return c
+end
+
+function BlockArrays.blocks(a::MultiField.BlockSparseMatrixAssembler)
+  return a.block_assemblers
+end
+
+function BlockArrays.blocks(f::MultiFieldFESpace{<:BlockMultiFieldStyle{NB,SB,P}}) where {NB,SB,P}
+  block_ranges = MultiField.get_block_ranges(NB,SB,P)
+  block_spaces = map(block_ranges) do range
+    (length(range) == 1) ? f[range[1]] : MultiFieldFESpace(f[range])
+  end
+  return block_spaces
+end
+
+function BlockArrays.blocks(f::GridapDistributed.DistributedMultiFieldFESpace{<:BlockMultiFieldStyle{NB,SB,P}}) where {NB,SB,P}
+  block_gids   = blocks(get_free_dof_ids(f))
+
+  block_ranges = MultiField.get_block_ranges(NB,SB,P)
+  block_spaces = map(block_ranges,block_gids) do range, gids
+    if (length(range) == 1) 
+      space = f[range[1]]
+    else
+      global_sf_spaces = f[range]
+      local_sf_spaces  = GridapDistributed.to_parray_of_arrays(map(local_views,global_sf_spaces))
+      local_mf_spaces  = map(MultiFieldFESpace,local_sf_spaces)
+      vector_type = GridapDistributed._find_vector_type(local_mf_spaces,gids)
+      space = MultiFieldFESpace(global_sf_spaces,local_mf_spaces,gids,vector_type)
+    end
+    space
+  end
+  return block_spaces
+end
diff --git a/src/BlockSolvers/BlockSolverInterfaces.jl b/src/BlockSolvers/BlockSolverInterfaces.jl
new file mode 100644
index 00000000..0708b076
--- /dev/null
+++ b/src/BlockSolvers/BlockSolverInterfaces.jl
@@ -0,0 +1,95 @@
+
+abstract type SolverBlock end
+
+abstract type LinearSolverBlock <: SolverBlock end
+abstract type NonlinearSolverBlock <: SolverBlock end
+
+is_nonlinear(::LinearSolverBlock) = false
+is_nonlinear(::NonlinearSolverBlock) = true
+
+function instantiate_block_cache(block::LinearSolverBlock,mat::AbstractMatrix)
+  @abstractmethod
+end
+function instantiate_block_cache(block::NonlinearSolverBlock,mat::AbstractMatrix,x::AbstractVector)
+  @abstractmethod
+end
+function instantiate_block_cache(block::LinearSolverBlock,mat::AbstractMatrix,x::AbstractVector)
+  instantiate_block_cache(block,mat)
+end
+
+function update_block_cache!(cache,block::LinearSolverBlock,mat::AbstractMatrix)
+  return cache
+end
+function update_block_cache!(cache,block::NonlinearSolverBlock,mat::AbstractMatrix,x::AbstractVector)
+  @abstractmethod
+end
+function update_block_cache!(cache,block::LinearSolverBlock,mat::AbstractMatrix,x::AbstractVector)
+  update_block_cache!(cache,block,mat)
+end
+
+# MatrixBlock
+
+struct MatrixBlock{A} <: LinearSolverBlock
+  mat :: A
+  function MatrixBlock(mat::AbstractMatrix)
+    A = typeof(mat)
+    return new{A}(mat)
+  end
+end
+
+instantiate_block_cache(block::MatrixBlock,::AbstractMatrix) = block.mat
+
+# SystemBlocks
+
+struct LinearSystemBlock <: LinearSolverBlock end
+struct NonlinearSystemBlock <: NonlinearSolverBlock end
+
+instantiate_block_cache(block::LinearSystemBlock,mat::AbstractMatrix) = mat
+instantiate_block_cache(block::NonlinearSystemBlock,mat::AbstractMatrix,::AbstractVector) = mat
+update_block_cache!(cache,block::NonlinearSystemBlock,mat::AbstractMatrix,::AbstractVector) = mat
+
+# BiformBlock/TriformBlock
+struct BiformBlock <: LinearSolverBlock
+  f     :: Function
+  trial :: FESpace
+  test  :: FESpace
+  assem :: Assembler
+  function BiformBlock(f::Function,
+                       trial::FESpace,
+                       test::FESpace,
+                       assem=SparseMatrixAssembler(trial,test))
+    return new(f,trial,test,assem)
+  end
+end
+
+struct TriformBlock <: NonlinearSolverBlock
+  f     :: Function
+  trial :: FESpace
+  test  :: FESpace
+  assem :: Assembler
+  function TriformBlock(f::Function,
+                        trial::FESpace,
+                        test::FESpace,
+                        assem=SparseMatrixAssembler(trial,test))
+    return new(f,trial,test,assem)
+  end
+end
+
+function instantiate_block_cache(block::BiformBlock,mat::AbstractMatrix)
+  return assemble_matrix(block.f,block.assem,block.trial,block.test)
+end
+
+function instantiate_block_cache(block::TriformBlock,mat::AbstractMatrix,x::AbstractVector)
+  uh = FEFunction(block.trial,x)
+  f(u,v) = block.f(uh,u,v)
+  return assemble_matrix(f,block.assem,block.trial,block.test)
+end
+
+function update_block_cache!(cache,block::TriformBlock,mat::AbstractMatrix,x::AbstractVector)
+  uh = FEFunction(block.trial,x)
+  f(u,v) = block.f(uh,u,v)
+  assemble_matrix!(mat,f,block.assem,block.trial,block.test)
+end
+
+# CompositeBlock
+# How do we deal with different sparsity patterns? Not trivial...
diff --git a/src/BlockSolvers/BlockSolvers.jl b/src/BlockSolvers/BlockSolvers.jl
new file mode 100644
index 00000000..921f01b7
--- /dev/null
+++ b/src/BlockSolvers/BlockSolvers.jl
@@ -0,0 +1,28 @@
+module BlockSolvers
+  using LinearAlgebra
+  using SparseArrays
+  using SparseMatricesCSR
+  using BlockArrays
+  using IterativeSolvers
+
+  using Gridap
+  using Gridap.Helpers, Gridap.Algebra, Gridap.CellData, Gridap.Arrays, Gridap.FESpaces, Gridap.MultiField
+  using PartitionedArrays
+  using GridapDistributed
+
+  using GridapSolvers.MultilevelTools
+  using GridapSolvers.SolverInterfaces
+
+  include("BlockFEOperators.jl")
+
+  include("BlockSolverInterfaces.jl")
+  include("BlockDiagonalSolvers.jl")
+  include("BlockTriangularSolvers.jl")
+
+  export BlockFEOperator
+
+  export MatrixBlock, LinearSystemBlock, NonlinearSystemBlock, BiformBlock, TriformBlock
+
+  export BlockDiagonalSolver
+  export BlockTriangularSolver
+end
diff --git a/src/BlockSolvers/BlockTriangularSolvers.jl b/src/BlockSolvers/BlockTriangularSolvers.jl
new file mode 100644
index 00000000..0fa9b622
--- /dev/null
+++ b/src/BlockSolvers/BlockTriangularSolvers.jl
@@ -0,0 +1,169 @@
+struct BlockTriangularSolver{T,N,A,B,C} <: Gridap.Algebra.LinearSolver
+  blocks  :: A
+  solvers :: B
+  coeffs  :: C
+  function BlockTriangularSolver(
+    blocks  :: AbstractMatrix{<:SolverBlock},
+    solvers :: AbstractVector{<:Gridap.Algebra.LinearSolver},
+    coeffs = fill(1.0,size(blocks)),
+    half   = :upper
+    )
+    N = length(solvers)
+    @check size(blocks,1) == size(blocks,2) == N
+    @check size(coeffs,1) == size(coeffs,2) == N
+    @check half ∈ (:upper,:lower)
+
+    A = typeof(blocks)
+    B = typeof(solvers)
+    C = typeof(coeffs)
+    return new{Val{half},N,A,B,C}(blocks,solvers,coeffs)
+  end
+end
+
+function BlockTriangularSolver(solvers::AbstractVector{<:Gridap.Algebra.LinearSolver}; 
+                               is_nonlinear::Matrix{Bool}=fill(false,(length(solvers),length(solvers))),
+                               coeffs=fill(1.0,size(is_nonlinear)),
+                               half=:upper)
+  blocks = map(nl -> nl ? NonlinearSystemBlock() : LinearSystemBlock(),is_nonlinear)
+  return BlockTriangularSolver(blocks,solvers,coeffs,half)
+end
+
+# Symbolic setup
+
+struct BlockTriangularSolverSS{A,B,C} <: Gridap.Algebra.SymbolicSetup
+  solver       :: A
+  block_ss     :: B
+  block_caches :: C
+end
+
+function Gridap.Algebra.symbolic_setup(solver::BlockTriangularSolver,mat::AbstractBlockMatrix)
+  mat_blocks   = blocks(mat)
+  block_caches = map(instantiate_block_cache,solver.blocks,mat_blocks)
+  block_ss     = map(symbolic_setup,solver.solvers,diag(block_caches))
+  return BlockTriangularSolverSS(solver,block_ss,block_caches)
+end
+
+function Gridap.Algebra.symbolic_setup(solver::BlockTriangularSolver{T,N},mat::AbstractBlockMatrix,x::AbstractBlockVector) where {T,N}
+  mat_blocks   = blocks(mat)
+  vec_blocks   = blocks(x)
+  block_caches = map(CartesianIndices(solver.blocks)) do I
+    instantiate_block_cache(solver.blocks[I],mat_blocks[I],vec_blocks[I[2]])
+  end
+  block_ss     = map(symbolic_setup,solver.solvers,diag(block_caches),vec_blocks)
+  return BlockTriangularSolverSS(solver,block_ss,block_caches)
+end
+
+# Numerical setup
+
+struct BlockTriangularSolverNS{T,A,B,C,D} <: Gridap.Algebra.NumericalSetup
+  solver       :: A
+  block_ns     :: B
+  block_caches :: C
+  work_caches  :: D
+  function BlockTriangularSolverNS(
+    solver::BlockTriangularSolver{T},
+    block_ns,block_caches,work_caches
+  ) where T
+    A = typeof(solver) 
+    B = typeof(block_ns)
+    C = typeof(block_caches)
+    D = typeof(work_caches)
+    return new{T,A,B,C,D}(solver,block_ns,block_caches,work_caches)
+  end
+end
+
+function Gridap.Algebra.numerical_setup(ss::BlockTriangularSolverSS,mat::AbstractBlockMatrix)
+  solver      = ss.solver
+  block_ns    = map(numerical_setup,ss.block_ss,diag(ss.block_caches))
+  work_caches = allocate_in_range(mat)
+  return BlockTriangularSolverNS(solver,block_ns,ss.block_caches,work_caches)
+end
+
+function Gridap.Algebra.numerical_setup(ss::BlockTriangularSolverSS,mat::AbstractBlockMatrix,x::AbstractBlockVector)
+  solver      = ss.solver
+  vec_blocks  = blocks(x)
+  block_ns    = map(numerical_setup,ss.block_ss,diag(ss.block_caches),vec_blocks)
+  work_caches = allocate_in_range(mat)
+  return BlockTriangularSolverNS(solver,block_ns,ss.block_caches,work_caches)
+end
+
+function Gridap.Algebra.numerical_setup!(ns::BlockTriangularSolverNS,mat::AbstractBlockMatrix)
+  solver       = ns.solver
+  mat_blocks   = blocks(mat)
+  block_caches = map(update_block_cache!,ns.block_caches,solver.blocks,mat_blocks)
+  map(diag(solver.blocks),ns.block_ns,diag(block_caches)) do bi, nsi, ci
+    if is_nonlinear(bi)
+      numerical_setup!(nsi,ci)
+    end
+  end
+  return ns
+end
+
+function Gridap.Algebra.numerical_setup!(ns::BlockTriangularSolverNS,mat::AbstractBlockMatrix,x::AbstractBlockVector)
+  solver       = ns.solver
+  mat_blocks   = blocks(mat)
+  vec_blocks   = blocks(x)
+  block_caches = map(CartesianIndices(solver.blocks)) do I
+    update_block_cache!(ns.block_caches[I],mat_blocks[I],vec_blocks[I[2]])
+  end
+  map(diag(solver.blocks),ns.block_ns,diag(block_caches),vec_blocks) do bi, nsi, ci, xi
+    if is_nonlinear(bi)
+      numerical_setup!(nsi,ci,xi)
+    end
+  end
+  return ns
+end
+
+function Gridap.Algebra.solve!(x::AbstractBlockVector,ns::BlockTriangularSolverNS{Val{:lower}},b::AbstractBlockVector)
+  @check blocklength(x) == blocklength(b) == length(ns.block_ns)
+  NB = length(ns.block_ns)
+  c, w = ns.solver.coeffs, ns.work_caches
+  mats = ns.block_caches
+  for iB in 1:NB
+    # Add lower off-diagonal contributions
+    wi  = w[Block(iB)]
+    copy!(wi,b[Block(iB)])
+    for jB in 1:iB-1
+      cij = c[iB,jB]
+      if abs(cij) > eps(cij)
+        xj = x[Block(jB)]
+        mul!(wi,mats[iB,jB],xj,-cij,1.0)
+      end
+    end
+
+    # Solve diagonal block
+    nsi = ns.block_ns[iB]
+    xi  = x[Block(iB)]
+    solve!(xi,nsi,wi)
+  end
+  return x
+end
+
+function Gridap.Algebra.solve!(x::AbstractBlockVector,ns::BlockTriangularSolverNS{Val{:upper}},b::AbstractBlockVector)
+  @check blocklength(x) == blocklength(b) == length(ns.block_ns)
+  NB = length(ns.block_ns)
+  c, w = ns.solver.coeffs, ns.work_caches
+  mats = ns.block_caches
+  for iB in NB:-1:1
+    # Add upper off-diagonal contributions
+    wi  = w[Block(iB)]
+    copy!(wi,b[Block(iB)])
+    for jB in iB+1:NB
+      cij = c[iB,jB]
+      if abs(cij) > eps(cij)
+        xj = x[Block(jB)]
+        mul!(wi,mats[iB,jB],xj,-cij,1.0)
+      end
+    end
+
+    # Solve diagonal block
+    nsi = ns.block_ns[iB]
+    xi = x[Block(iB)]
+    solve!(xi,nsi,wi)
+  end
+  return x
+end
+
+function LinearAlgebra.ldiv!(x,ns::BlockTriangularSolverNS,b)
+  solve!(x,ns,b)
+end
diff --git a/src/GridapSolvers.jl b/src/GridapSolvers.jl
index 12870f53..516819f2 100644
--- a/src/GridapSolvers.jl
+++ b/src/GridapSolvers.jl
@@ -2,13 +2,17 @@ module GridapSolvers
 
   include("SolverInterfaces/SolverInterfaces.jl")
   include("MultilevelTools/MultilevelTools.jl")
+  include("BlockSolvers/BlockSolvers.jl")
   include("LinearSolvers/LinearSolvers.jl")
   include("PatchBasedSmoothers/PatchBasedSmoothers.jl")
+  include("NonlinearSolvers/NonlinearSolvers.jl")
 
   using GridapSolvers.SolverInterfaces
   using GridapSolvers.MultilevelTools
+  using GridapSolvers.BlockSolvers
   using GridapSolvers.LinearSolvers
   using GridapSolvers.PatchBasedSmoothers
+  using GridapSolvers.NonlinearSolvers
 
   # MultilevelTools
   export get_parts, generate_level_parts, generate_subparts
@@ -25,6 +29,9 @@ module GridapSolvers
   export RestrictionOperator, ProlongationOperator
   export setup_transfer_operators
 
+  # BlockSolvers
+  export BlockDiagonalSolver
+
   # LinearSolvers
   export JacobiLinearSolver
   export RichardsonSmoother
@@ -37,7 +44,10 @@ module GridapSolvers
   export IS_MINRESSolver
   export IS_SSORSolver
 
+  export CGSolver
+  export MINRESSolver
   export GMRESSolver
+  export FGMRESSolver
 
   # PatchBasedSmoothers
   export PatchDecomposition
diff --git a/src/LinearSolvers/BlockDiagonalSmoothers.jl b/src/LinearSolvers/BlockDiagonalSmoothers.jl
deleted file mode 100644
index 8e34511e..00000000
--- a/src/LinearSolvers/BlockDiagonalSmoothers.jl
+++ /dev/null
@@ -1,80 +0,0 @@
-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 4b8651d7..8ca782f1 100644
--- a/src/LinearSolvers/GMGLinearSolvers.jl
+++ b/src/LinearSolvers/GMGLinearSolvers.jl
@@ -1,40 +1,36 @@
-struct GMGLinearSolver{A,B,C,D,E,F,G,H} <: Gridap.Algebra.LinearSolver
-  mh              :: ModelHierarchy
-  smatrices       :: A
-  interp          :: B
-  restrict        :: C
-  pre_smoothers   :: D
-  post_smoothers  :: E
-  coarsest_solver :: F
-  maxiter         :: G
-  rtol            :: H
-  verbose         :: Bool
+struct GMGLinearSolver{A,B,C,D,E,F,G} <: Gridap.Algebra.LinearSolver
+  mh              :: A
+  smatrices       :: B
+  interp          :: C
+  restrict        :: D
+  pre_smoothers   :: E
+  post_smoothers  :: F
+  coarsest_solver :: G
   mode            :: Symbol
-end
-
-function GMGLinearSolver(mh,smatrices,interp,restrict;
-      pre_smoothers   = Fill(RichardsonSmoother(JacobiLinearSolver(),10),num_levels(mh)-1),
-      post_smoothers  = pre_smoothers,
-      coarsest_solver = Gridap.Algebra.BackslashSolver(),
-      maxiter         = 100,
-      rtol            = 1.0e-06,
-      verbose::Bool   = false,
-      mode            = :preconditioner)
-
-  Gridap.Helpers.@check mode ∈ [:preconditioner,:solver]
-  Gridap.Helpers.@check isa(maxiter,Integer)
-  Gridap.Helpers.@check isa(rtol,Real)
-
-  A=typeof(smatrices)
-  B=typeof(interp)
-  C=typeof(restrict)
-  D=typeof(pre_smoothers)
-  E=typeof(post_smoothers)
-  F=typeof(coarsest_solver)
-  G=typeof(maxiter)
-  H=typeof(rtol)
-  return GMGLinearSolver{A,B,C,D,E,F,G,H}(mh,smatrices,interp,restrict,pre_smoothers,post_smoothers,
-                                          coarsest_solver,maxiter,rtol,verbose,mode)
+  log             :: ConvergenceLog{Float64}
+end
+
+function GMGLinearSolver(
+  mh,smatrices,interp,restrict;
+  pre_smoothers   = Fill(RichardsonSmoother(JacobiLinearSolver(),10),num_levels(mh)-1),
+  post_smoothers  = pre_smoothers,
+  coarsest_solver = Gridap.Algebra.LUSolver(),
+  mode::Symbol    = :preconditioner,
+  maxiter = 100, atol = 1.0e-14, rtol = 1.0e-08, verbose = false,
+)
+  @check mode ∈ [:preconditioner,:solver]
+  tols = SolverTolerances{Float64}(;maxiter=maxiter,atol=atol,rtol=rtol)
+  log  = ConvergenceLog("GMG",tols;verbose=verbose)
+  
+  A = typeof(mh)
+  B = typeof(smatrices)
+  C = typeof(interp)
+  D = typeof(restrict)
+  E = typeof(pre_smoothers)
+  F = typeof(post_smoothers)
+  G = typeof(coarsest_solver)
+  return GMGLinearSolver{A,B,C,D,E,F,G}(mh,smatrices,interp,restrict,pre_smoothers,post_smoothers,
+                                        coarsest_solver,mode,log)
 end
 
 struct GMGSymbolicSetup <: Gridap.Algebra.SymbolicSetup
@@ -52,39 +48,35 @@ struct GMGNumericalSetup{A,B,C,D,E} <: Gridap.Algebra.NumericalSetup
   post_smoothers_caches  :: C
   coarsest_solver_cache  :: D
   work_vectors           :: E
+end
 
-  function GMGNumericalSetup(ss::GMGSymbolicSetup)
-    mh              = ss.solver.mh
-    pre_smoothers   = ss.solver.pre_smoothers
-    post_smoothers  = ss.solver.post_smoothers
-    smatrices       = ss.solver.smatrices
-    coarsest_solver = ss.solver.coarsest_solver
-
-    finest_level_cache = setup_finest_level_cache(mh,smatrices)
-    work_vectors = allocate_work_vectors(mh,smatrices)
-    pre_smoothers_caches = setup_smoothers_caches(mh,pre_smoothers,smatrices)
-    if (!(pre_smoothers === post_smoothers))
-      post_smoothers_caches = setup_smoothers_caches(mh,post_smoothers,smatrices)
-    else
-      post_smoothers_caches = pre_smoothers_caches
-    end
-    #coarsest_solver_cache = setup_coarsest_solver_cache(mh,coarsest_solver,smatrices)
-    coarsest_solver_cache = coarse_solver_caches(mh,coarsest_solver,smatrices)
-
-    A = typeof(finest_level_cache)
-    B = typeof(pre_smoothers_caches)
-    C = typeof(post_smoothers_caches)
-    D = typeof(coarsest_solver_cache)
-    E = typeof(work_vectors)
-    return new{A,B,C,D,E}(ss.solver,finest_level_cache,pre_smoothers_caches,post_smoothers_caches,coarsest_solver_cache,work_vectors)
+function Gridap.Algebra.numerical_setup(ss::GMGSymbolicSetup,mat::AbstractMatrix)
+  mh              = ss.solver.mh
+  pre_smoothers   = ss.solver.pre_smoothers
+  post_smoothers  = ss.solver.post_smoothers
+  smatrices       = ss.solver.smatrices
+  coarsest_solver = ss.solver.coarsest_solver
+
+  smatrices[1] = mat
+  finest_level_cache = gmg_finest_level_cache(mh,smatrices)
+  work_vectors = gmg_work_vectors(mh,smatrices)
+  pre_smoothers_caches = gmg_smoothers_caches(mh,pre_smoothers,smatrices)
+  if !(pre_smoothers === post_smoothers)
+    post_smoothers_caches = gmg_smoothers_caches(mh,post_smoothers,smatrices)
+  else
+    post_smoothers_caches = pre_smoothers_caches
   end
+  coarsest_solver_cache = gmg_coarse_solver_caches(mh,coarsest_solver,smatrices,work_vectors)
+
+  return GMGNumericalSetup(ss.solver,finest_level_cache,pre_smoothers_caches,post_smoothers_caches,coarsest_solver_cache,work_vectors)
 end
 
-function Gridap.Algebra.numerical_setup(ss::GMGSymbolicSetup,mat::AbstractMatrix)
-  return GMGNumericalSetup(ss)
+function Gridap.Algebra.numerical_setup!(ss::GMGNumericalSetup,mat::AbstractMatrix)
+  # TODO: This does not modify all matrices... How should we deal with this?
+  ns.solver.smatrices[1] = mat
 end
 
-function setup_finest_level_cache(mh::ModelHierarchy,smatrices::Vector{<:AbstractMatrix})
+function gmg_finest_level_cache(mh::ModelHierarchy,smatrices::Vector{<:AbstractMatrix})
   cache = nothing
   parts = get_level_parts(mh,1)
   if i_am_in(parts)
@@ -95,7 +87,7 @@ function setup_finest_level_cache(mh::ModelHierarchy,smatrices::Vector{<:Abstrac
   return cache
 end
 
-function setup_smoothers_caches(mh::ModelHierarchy,smoothers::AbstractVector{<:LinearSolver},smatrices::Vector{<:AbstractMatrix})
+function gmg_smoothers_caches(mh::ModelHierarchy,smoothers::AbstractVector{<:LinearSolver},smatrices::Vector{<:AbstractMatrix})
   Gridap.Helpers.@check length(smoothers) == num_levels(mh)-1
   nlevs = num_levels(mh)
   # Last (i.e., coarsest) level does not need pre-/post-smoothing
@@ -110,64 +102,34 @@ function setup_smoothers_caches(mh::ModelHierarchy,smoothers::AbstractVector{<:L
   return caches
 end
 
-function coarse_solver_caches(mh,s,mats)
+function gmg_coarse_solver_caches(mh,s,mats,work_vectors)
   cache = nothing
   nlevs = num_levels(mh)
   parts = get_level_parts(mh,nlevs)
   if i_am_in(parts)
     mat = mats[nlevs]
+    _, _, xH, rH = work_vectors[nlevs-1]
     cache = numerical_setup(symbolic_setup(s, mat), mat)
-  end
-  return cache
-end
-
-function setup_coarsest_solver_cache(mh::ModelHierarchy,coarsest_solver::LinearSolver,smatrices::Vector{<:AbstractMatrix})
-  cache = nothing
-  nlevs = num_levels(mh)
-  parts = get_level_parts(mh,nlevs)
-  if i_am_in(parts)
-    mat = smatrices[nlevs]
-    if (num_parts(parts) == 1) # Serial
-      cache = map(own_values(mat)) do Ah
-        ss  = symbolic_setup(coarsest_solver, Ah)
-        numerical_setup(ss, Ah)
-      end
-      cache = PartitionedArrays.getany(cache)
-    else # Parallel
-      ss = symbolic_setup(coarsest_solver, mat)
-      cache = numerical_setup(ss, mat)
+    if isa(s,PETScLinearSolver)
+      cache = CachedPETScNS(cache, xH, rH)
     end
   end
   return cache
 end
 
-function setup_coarsest_solver_cache(mh::ModelHierarchy,coarsest_solver::PETScLinearSolver,smatrices::Vector{<:AbstractMatrix})
-  cache = nothing
+function gmg_work_vectors(mh::ModelHierarchy,smatrices::Vector{<:AbstractMatrix})
   nlevs = num_levels(mh)
-  parts = get_level_parts(mh,nlevs)
-  if i_am_in(parts)
-    mat   = smatrices[nlevs]
-    if (num_parts(parts) == 1) # Serial
-      cache = map(own_values(mat)) do Ah
-        rh  = convert(PETScVector,fill(0.0,size(A,2)))
-        xh  = convert(PETScVector,fill(0.0,size(A,2)))
-        ss  = symbolic_setup(coarsest_solver, Ah)
-        ns  = numerical_setup(ss, Ah)
-        return ns, xh, rh
-      end
-      cache = PartitionedArrays.getany(cache)
-    else # Parallel
-      rh = convert(PETScVector,pfill(0.0,partition(axes(mat,2))))
-      xh = convert(PETScVector,pfill(0.0,partition(axes(mat,2))))
-      ss = symbolic_setup(coarsest_solver, mat)
-      ns = numerical_setup(ss, mat)
-      cache = ns, xh, rh
+  work_vectors = Vector{Any}(undef,nlevs-1)
+  for i = 1:nlevs-1
+    parts = get_level_parts(mh,i)
+    if i_am_in(parts)
+      work_vectors[i] = gmg_work_vectors(mh,smatrices,i)
     end
   end
-  return cache
+  return work_vectors
 end
 
-function allocate_level_work_vectors(mh::ModelHierarchy,smatrices::Vector{<:AbstractMatrix},lev::Integer)
+function gmg_work_vectors(mh::ModelHierarchy,smatrices::Vector{<:AbstractMatrix},lev::Integer)
   dxh   = allocate_in_domain(smatrices[lev])
   Adxh  = allocate_in_range(smatrices[lev])
 
@@ -183,52 +145,12 @@ function allocate_level_work_vectors(mh::ModelHierarchy,smatrices::Vector{<:Abst
   return dxh, Adxh, dxH, rH
 end
 
-function allocate_work_vectors(mh::ModelHierarchy,smatrices::Vector{<:AbstractMatrix})
-  nlevs = num_levels(mh)
-  work_vectors = Vector{Any}(undef,nlevs-1)
-  for i = 1:nlevs-1
-    parts = get_level_parts(mh,i)
-    if i_am_in(parts)
-      work_vectors[i] = allocate_level_work_vectors(mh,smatrices,i)
-    end
-  end
-  return work_vectors
-end
-
-function solve_coarsest_level!(parts::AbstractArray,::LinearSolver,xh::PVector,rh::PVector,caches)
-  if (num_parts(parts) == 1)
-    map(own_values(xh),own_values(rh)) do xh, rh
-      solve!(xh,caches,rh)
-    end
-  else
-    solve!(xh,caches,rh)
-  end
-end
-
-function solve_coarsest_level!(parts::AbstractArray,::PETScLinearSolver,xh::PVector,rh::PVector,caches)
-  solver_ns, xh_petsc, rh_petsc = caches
-  if (num_parts(parts) == 1)
-    map(own_values(xh),own_values(rh)) do xh, rh
-      copy!(rh_petsc,rh)
-      solve!(xh_petsc,solver_ns,rh_petsc)
-      copy!(xh,xh_petsc)
-    end
-  else
-    copy!(rh_petsc,rh)
-    solve!(xh_petsc,solver_ns,rh_petsc)
-    copy!(xh,xh_petsc)
-  end
-end
-
-function apply_GMG_level!(lev::Integer,xh::Union{PVector,Nothing},rh::Union{PVector,Nothing},ns::GMGNumericalSetup;verbose=false)
+function apply_GMG_level!(lev::Integer,xh::Union{PVector,Nothing},rh::Union{PVector,Nothing},ns::GMGNumericalSetup)
   mh = ns.solver.mh
   parts = get_level_parts(mh,lev)
   if i_am_in(parts)
     if (lev == num_levels(mh)) 
       ## Coarsest level
-      #coarsest_solver = ns.solver.coarsest_solver
-      #coarsest_solver_cache = ns.coarsest_solver_cache
-      #solve_coarsest_level!(parts,coarsest_solver,xh,rh,coarsest_solver_cache)
       solve!(xh, ns.coarsest_solver_cache, rh)
     else 
       ## General case
@@ -244,7 +166,7 @@ function apply_GMG_level!(lev::Integer,xh::Union{PVector,Nothing},rh::Union{PVec
 
       # Apply next_level
       !isa(dxH,Nothing) && fill!(dxH,0.0)
-      apply_GMG_level!(lev+1,dxH,rH,ns;verbose=verbose)
+      apply_GMG_level!(lev+1,dxH,rH,ns)
 
       # Interpolate dxH in finer space
       mul!(dxh,interp,dxH)
@@ -261,14 +183,9 @@ function apply_GMG_level!(lev::Integer,xh::Union{PVector,Nothing},rh::Union{PVec
 end
 
 function Gridap.Algebra.solve!(x::AbstractVector,ns::GMGNumericalSetup,b::AbstractVector)
+  mode = ns.solver.mode
+  log  = ns.solver.log
 
-  mh      = ns.solver.mh
-  maxiter = ns.solver.maxiter
-  rtol    = ns.solver.rtol
-  verbose = ns.solver.verbose
-  mode    = ns.solver.mode
-
-  # TODO: When running in preconditioner mode, do we really need to compute the norm? It's a global com....
   rh = ns.finest_level_cache
   if (mode == :preconditioner)
     fill!(x,0.0)
@@ -279,29 +196,16 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::GMGNumericalSetup,b::Abstra
     rh .= b .- rh
   end
 
-  nrm_r0 = norm(rh)
-  nrm_r  = nrm_r0
-  current_iter = 0
-  rel_res = nrm_r / nrm_r0
-  parts = get_level_parts(mh,1)
-
-  if i_am_main(parts) && verbose
-    @printf "%6s  %12s" "Iter" "Rel res\n"
-    @printf "%6i  %12.4e\n" current_iter rel_res
+  res  = norm(rh)
+  done = init!(log,res)
+  while !done
+    apply_GMG_level!(1,x,rh,ns)
+    res  = norm(rh)
+    done = update!(log,res)
   end
 
-  while (current_iter < maxiter) && (rel_res > rtol)
-    apply_GMG_level!(1,x,rh,ns;verbose=verbose)
-
-    nrm_r   = norm(rh)
-    rel_res = nrm_r / nrm_r0
-    current_iter += 1
-    if i_am_main(parts) && verbose
-      @printf "%6i  %12.4e\n" current_iter rel_res
-    end
-  end
-  converged = (rel_res < rtol)
-  return current_iter, converged
+  finalize!(log,res)
+  return x
 end
 
 function LinearAlgebra.ldiv!(x::AbstractVector,ns::GMGNumericalSetup,b::AbstractVector)
diff --git a/src/LinearSolvers/Krylov/CGSolvers.jl b/src/LinearSolvers/Krylov/CGSolvers.jl
index ed35246f..081da65d 100644
--- a/src/LinearSolvers/Krylov/CGSolvers.jl
+++ b/src/LinearSolvers/Krylov/CGSolvers.jl
@@ -46,6 +46,13 @@ end
 function Gridap.Algebra.numerical_setup!(ns::CGNumericalSetup, A::AbstractMatrix)
   numerical_setup!(ns.Pl_ns,A)
   ns.A = A
+  return ns
+end
+
+function Gridap.Algebra.numerical_setup!(ns::CGNumericalSetup, A::AbstractMatrix, x::AbstractVector)
+  numerical_setup!(ns.Pl_ns,A,x)
+  ns.A = A
+  return ns
 end
 
 function Gridap.Algebra.solve!(x::AbstractVector,ns::CGNumericalSetup,b::AbstractVector)
@@ -55,7 +62,9 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::CGNumericalSetup,b::Abstrac
 
   # Initial residual
   mul!(w,A,x); r .= b .- w
-  fill!(p,0.0); γ = 1.0
+  fill!(p,zero(eltype(p)))
+  fill!(z,zero(eltype(z)))
+  γ = one(eltype(p))
 
   res  = norm(r)
   done = init!(log,res)
diff --git a/src/LinearSolvers/Krylov/FGMRESSolvers.jl b/src/LinearSolvers/Krylov/FGMRESSolvers.jl
index f9cc9f7a..412f8c24 100644
--- a/src/LinearSolvers/Krylov/FGMRESSolvers.jl
+++ b/src/LinearSolvers/Krylov/FGMRESSolvers.jl
@@ -2,18 +2,25 @@
 # FGMRES Solver
 struct FGMRESSolver <: Gridap.Algebra.LinearSolver
   m         :: Int
+  restart   :: Bool
+  m_add     :: Int
   Pr        :: Gridap.Algebra.LinearSolver
   Pl        :: Union{Gridap.Algebra.LinearSolver,Nothing}
-  outer_log :: ConvergenceLog{Float64}
-  inner_log :: ConvergenceLog{Float64}
+  log :: ConvergenceLog{Float64}
 end
 
-function FGMRESSolver(m,Pr;Pl=nothing,maxiter=100,atol=1e-12,rtol=1.e-6,verbose=false,name="FGMRES")
-  outer_tols = SolverTolerances{Float64}(maxiter=maxiter,atol=atol,rtol=rtol)
-  outer_log  = ConvergenceLog(name,outer_tols,verbose=verbose)
-  inner_tols = SolverTolerances{Float64}(maxiter=m,atol=atol,rtol=rtol)
-  inner_log  = ConvergenceLog("$(name)_inner",inner_tols,verbose=verbose,nested=true)
-  return FGMRESSolver(m,Pr,Pl,outer_log,inner_log)
+function FGMRESSolver(m,Pr;Pl=nothing,restart=false,m_add=1,maxiter=100,atol=1e-12,rtol=1.e-6,verbose=false,name="FGMRES")
+  tols = SolverTolerances{Float64}(maxiter=maxiter,atol=atol,rtol=rtol)
+  log  = ConvergenceLog(name,tols,verbose=verbose)
+  return FGMRESSolver(m,restart,m_add,Pr,Pl,log)
+end
+
+function restart(s::FGMRESSolver,k::Int)
+  if s.restart && (k > s.m)
+    print_message(s.log,"Restarting Krylov basis.")
+    return true
+  end
+  return false
 end
 
 AbstractTrees.children(s::FGMRESSolver) = [s.Pr,s.Pl]
@@ -48,10 +55,42 @@ function get_solver_caches(solver::FGMRESSolver,A)
   return (V,Z,zl,H,g,c,s)
 end
 
+function krylov_cache_length(ns::FGMRESNumericalSetup)
+  V, _, _, _, _, _, _ = ns.caches
+  return length(V) - 1
+end
+
+function expand_krylov_caches!(ns::FGMRESNumericalSetup)
+  V, Z, zl, H, g, c, s = ns.caches
+
+  m = krylov_cache_length(ns)
+  m_add = ns.solver.m_add
+  m_new = m + m_add
+
+  for _ in 1:m_add
+    push!(V,allocate_in_domain(ns.A))
+    push!(Z,allocate_in_domain(ns.A))
+  end
+  H_new = zeros(eltype(H),m_new+1,m_new); H_new[1:m+1,1:m] .= H
+  g_new = zeros(eltype(g),m_new+1); g_new[1:m+1] .= g
+  c_new = zeros(eltype(c),m_new); c_new[1:m] .= c
+  s_new = zeros(eltype(s),m_new); s_new[1:m] .= s
+  ns.caches = (V,Z,zl,H_new,g_new,c_new,s_new)
+  return H_new,g_new,c_new,s_new
+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)
+  Pl_ns  = !isnothing(solver.Pl) ? numerical_setup(symbolic_setup(solver.Pl,A),A) : nothing
+  caches = get_solver_caches(solver,A)
+  return FGMRESNumericalSetup(solver,A,Pr_ns,Pl_ns,caches)
+end
+
+function Gridap.Algebra.numerical_setup(ss::FGMRESSymbolicSetup, A::AbstractMatrix, x::AbstractVector)
+  solver = ss.solver
+  Pr_ns  = numerical_setup(symbolic_setup(solver.Pr,A,x),A,x)
+  Pl_ns  = !isnothing(solver.Pl) ? numerical_setup(symbolic_setup(solver.Pl,A,x),A,x) : nothing
   caches = get_solver_caches(solver,A)
   return FGMRESNumericalSetup(solver,A,Pr_ns,Pl_ns,caches)
 end
@@ -62,12 +101,26 @@ function Gridap.Algebra.numerical_setup!(ns::FGMRESNumericalSetup, A::AbstractMa
     numerical_setup!(ns.Pl_ns,A)
   end
   ns.A = A
+  return ns
+end
+
+function Gridap.Algebra.numerical_setup!(ns::FGMRESNumericalSetup, A::AbstractMatrix, x::AbstractVector)
+  numerical_setup!(ns.Pr_ns,A,x)
+  if !isa(ns.Pl_ns,Nothing)
+    numerical_setup!(ns.Pl_ns,A,x)
+  end
+  ns.A = A
+  return ns
 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
-  log, ilog = solver.outer_log, solver.inner_log
   V, Z, zl, H, g, c, s = caches
+  m   = krylov_cache_length(ns)
+  log = solver.log
+
+  fill!(V[1],zero(eltype(V[1])))
+  fill!(zl,zero(eltype(zl)))
 
   # Initial residual
   krylov_residual!(V[1],x,A,b,Pl,zl)
@@ -79,9 +132,16 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::FGMRESNumericalSetup,b::Abs
     V[1] ./= β
     fill!(H,0.0)
     fill!(g,0.0); g[1] = β
-    idone = init!(ilog,β)
-    while !idone
+    while !done && !restart(solver,j)
+      # Expand Krylov basis if needed
+      if j > m  
+        H, g, c, s = expand_krylov_caches!(ns)
+        m = krylov_cache_length(ns)
+      end
+
       # Arnoldi orthogonalization by Modified Gram-Schmidt
+      fill!(V[j+1],zero(eltype(V[j+1])))
+      fill!(Z[j],zero(eltype(Z[j])))
       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])
@@ -104,7 +164,7 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::FGMRESNumericalSetup,b::Abs
 
       β  = abs(g[j+1])
       j += 1
-      idone = update!(ilog,β)
+      done = update!(log,β)
     end
     j = j-1
 
@@ -118,7 +178,6 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::FGMRESNumericalSetup,b::Abs
       x .+= g[i] .* Z[i]
     end
     krylov_residual!(V[1],x,A,b,Pl,zl)
-    done = update!(log,β)
   end
 
   finalize!(log,β)
diff --git a/src/LinearSolvers/Krylov/GMRESSolvers.jl b/src/LinearSolvers/Krylov/GMRESSolvers.jl
index 9a816b56..b7afa0bf 100644
--- a/src/LinearSolvers/Krylov/GMRESSolvers.jl
+++ b/src/LinearSolvers/Krylov/GMRESSolvers.jl
@@ -1,18 +1,25 @@
 # GMRES Solver
 struct GMRESSolver <: Gridap.Algebra.LinearSolver
-  m         :: Int
-  Pr        :: Union{Gridap.Algebra.LinearSolver,Nothing}
-  Pl        :: Union{Gridap.Algebra.LinearSolver,Nothing}
-  outer_log :: ConvergenceLog{Float64}
-  inner_log :: ConvergenceLog{Float64}
+  m       :: Int
+  restart :: Bool
+  m_add   :: Int
+  Pr      :: Union{Gridap.Algebra.LinearSolver,Nothing}
+  Pl      :: Union{Gridap.Algebra.LinearSolver,Nothing}
+  log     :: ConvergenceLog{Float64}
 end
 
-function GMRESSolver(m;Pr=nothing,Pl=nothing,maxiter=100,atol=1e-12,rtol=1.e-6,verbose=false,name="GMRES")
-  outer_tols = SolverTolerances{Float64}(maxiter=maxiter,atol=atol,rtol=rtol)
-  outer_log  = ConvergenceLog(name,outer_tols,verbose=verbose)
-  inner_tols = SolverTolerances{Float64}(maxiter=m,atol=atol,rtol=rtol)
-  inner_log  = ConvergenceLog("$(name)_inner",inner_tols,verbose=verbose,nested=true)
-  return GMRESSolver(m,Pr,Pl,outer_log,inner_log)
+function GMRESSolver(m;Pr=nothing,Pl=nothing,restart=false,m_add=1,maxiter=100,atol=1e-12,rtol=1.e-6,verbose=false,name="GMRES")
+  tols = SolverTolerances{Float64}(maxiter=maxiter,atol=atol,rtol=rtol)
+  log  = ConvergenceLog(name,tols,verbose=verbose)
+  return GMRESSolver(m,restart,m_add,Pr,Pl,log)
+end
+
+function restart(s::GMRESSolver,k::Int)
+  if s.restart && (k > s.m)
+    print_message(s.log,"Restarting Krylov basis.")
+    return true
+  end
+  return false
 end
 
 AbstractTrees.children(s::GMRESSolver) = [s.Pr,s.Pl]
@@ -37,7 +44,7 @@ function get_solver_caches(solver::GMRESSolver,A)
   m, Pl, Pr = solver.m, solver.Pl, solver.Pr
 
   V  = [allocate_in_domain(A) for i in 1:m+1]
-  zr = !isa(Pr,Nothing) ? allocate_in_domain(A) : nothing
+  zr = !isnothing(Pr) ? allocate_in_domain(A) : nothing
   zl = allocate_in_domain(A)
 
   H = zeros(m+1,m)  # Hessenberg matrix
@@ -47,10 +54,33 @@ function get_solver_caches(solver::GMRESSolver,A)
   return (V,zr,zl,H,g,c,s)
 end
 
+function krylov_cache_length(ns::GMRESNumericalSetup)
+  V, _, _, _, _, _, _ = ns.caches
+  return length(V) - 1
+end
+
+function expand_krylov_caches!(ns::GMRESNumericalSetup)
+  V, zr, zl, H, g, c, s = ns.caches
+
+  m = krylov_cache_length(ns)
+  m_add = ns.solver.m_add
+  m_new = m + m_add
+
+  for _ in 1:m_add
+    push!(V,allocate_in_domain(ns.A))
+  end
+  H_new = zeros(eltype(H),m_new+1,m_new); H_new[1:m+1,1:m] .= H
+  g_new = zeros(eltype(g),m_new+1); g_new[1:m+1] .= g
+  c_new = zeros(eltype(c),m_new); c_new[1:m] .= c
+  s_new = zeros(eltype(s),m_new); s_new[1:m] .= s
+  ns.caches = (V,zr,zl,H_new,g_new,c_new,s_new)
+  return H_new,g_new,c_new,s_new
+end
+
 function Gridap.Algebra.numerical_setup(ss::GMRESSymbolicSetup, 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)
+  Pr_ns  = !isnothing(solver.Pr) ? numerical_setup(symbolic_setup(solver.Pr,A),A) : nothing
+  Pl_ns  = !isnothing(solver.Pl) ? numerical_setup(symbolic_setup(solver.Pl,A),A) : nothing
   caches = get_solver_caches(solver,A)
   return GMRESNumericalSetup(solver,A,Pr_ns,Pl_ns,caches)
 end
@@ -63,12 +93,29 @@ function Gridap.Algebra.numerical_setup!(ns::GMRESNumericalSetup, A::AbstractMat
     numerical_setup!(ns.Pl_ns,A)
   end
   ns.A = A
+  return ns
+end
+
+function Gridap.Algebra.numerical_setup!(ns::GMRESNumericalSetup, A::AbstractMatrix, x::AbstractVector)
+  if !isa(ns.Pr_ns,Nothing)
+    numerical_setup!(ns.Pr_ns,A,x)
+  end
+  if !isa(ns.Pl_ns,Nothing)
+    numerical_setup!(ns.Pl_ns,A,x)
+  end
+  ns.A = A
+  return ns
 end
 
 function Gridap.Algebra.solve!(x::AbstractVector,ns::GMRESNumericalSetup,b::AbstractVector)
   solver, A, Pl, Pr, caches = ns.solver, ns.A, ns.Pl_ns, ns.Pr_ns, ns.caches
-  log, ilog = solver.outer_log, solver.inner_log
   V, zr, zl, H, g, c, s = caches
+  m   = krylov_cache_length(ns)
+  log = solver.log
+
+  fill!(V[1],zero(eltype(V[1])))
+  !isnothing(zr) && fill!(zr,zero(eltype(zr)))
+  fill!(zl,zero(eltype(zl)))
 
   # Initial residual
   krylov_residual!(V[1],x,A,b,Pl,zl)
@@ -80,9 +127,15 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::GMRESNumericalSetup,b::Abst
     V[1] ./= β
     fill!(H,0.0)
     fill!(g,0.0); g[1] = β
-    idone = init!(ilog,β)
-    while !idone
+    while !done && !restart(solver,j)
+      # Expand Krylov basis if needed
+      if j > m  
+        H, g, c, s = expand_krylov_caches!(ns)
+        m = krylov_cache_length(ns)
+      end
+
       # Arnoldi orthogonalization by Modified Gram-Schmidt
+      fill!(V[j+1],zero(eltype(V[j+1])))
       krylov_mul!(V[j+1],A,V[j],Pr,Pl,zr,zl)
       for i in 1:j
         H[i,j] = dot(V[j+1],V[i])
@@ -105,7 +158,7 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::GMRESNumericalSetup,b::Abst
 
       β  = abs(g[j+1])
       j += 1
-      idone = update!(ilog,β)
+      done = update!(log,β)
     end
     j = j-1
 
@@ -128,8 +181,8 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::GMRESNumericalSetup,b::Abst
       x .+= zr
     end
     krylov_residual!(V[1],x,A,b,Pl,zl)
-    done = update!(log,β)
   end
+
   finalize!(log,β)
   return x
 end
diff --git a/src/LinearSolvers/Krylov/MINRESSolvers.jl b/src/LinearSolvers/Krylov/MINRESSolvers.jl
index fcf4504b..8feb9fb1 100644
--- a/src/LinearSolvers/Krylov/MINRESSolvers.jl
+++ b/src/LinearSolvers/Krylov/MINRESSolvers.jl
@@ -62,6 +62,17 @@ function Gridap.Algebra.numerical_setup!(ns::MINRESNumericalSetup, A::AbstractMa
   ns.A = A
 end
 
+function Gridap.Algebra.numerical_setup!(ns::MINRESNumericalSetup, A::AbstractMatrix, x::AbstractVector)
+  if !isa(ns.Pr_ns,Nothing)
+    numerical_setup!(ns.Pr_ns,A,x)
+  end
+  if !isa(ns.Pl_ns,Nothing)
+    numerical_setup!(ns.Pl_ns,A,x)
+  end
+  ns.A = A
+  return ns
+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
   V, W, zr, zl, H, g, c, s = caches
diff --git a/src/LinearSolvers/LinearSolvers.jl b/src/LinearSolvers/LinearSolvers.jl
index 5eb500ad..863816e3 100644
--- a/src/LinearSolvers/LinearSolvers.jl
+++ b/src/LinearSolvers/LinearSolvers.jl
@@ -43,6 +43,7 @@ include("Krylov/FGMRESSolvers.jl")
 include("Krylov/MINRESSolvers.jl")
 
 include("PETSc/PETScUtils.jl")
+include("PETSc/PETScCaches.jl")
 include("PETSc/ElasticitySolvers.jl")
 include("PETSc/HipmairXuSolvers.jl")
 
@@ -51,7 +52,6 @@ include("JacobiLinearSolvers.jl")
 include("RichardsonSmoothers.jl")
 include("SymGaussSeidelSmoothers.jl")
 include("GMGLinearSolvers.jl")
-include("BlockDiagonalSmoothers.jl")
 include("IterativeLinearSolvers.jl")
 include("SchurComplementSolvers.jl")
 
diff --git a/src/LinearSolvers/PETSc/PETScCaches.jl b/src/LinearSolvers/PETSc/PETScCaches.jl
new file mode 100644
index 00000000..a2845149
--- /dev/null
+++ b/src/LinearSolvers/PETSc/PETScCaches.jl
@@ -0,0 +1,33 @@
+
+"""
+ Notes on this structure: 
+
+ When converting julia vectors/PVectors to PETSc vectors, we purposely create aliasing 
+ of the vector values. This means we can avoid copying data from one to another before solving, 
+ but we need to be careful about it. 
+
+ This structure takes care of this, and makes sure you do not attempt to solve the system 
+ with julia vectors that are not the ones you used to create the solver cache.
+"""
+struct CachedPETScNS{TM,A}
+  ns     :: GridapPETSc.PETScLinearSolverNS{TM}
+  X      :: PETScVector
+  B      :: PETScVector
+  owners :: A
+  function CachedPETScNS(ns::GridapPETSc.PETScLinearSolverNS{TM},x::AbstractVector,b::AbstractVector) where TM
+    X = convert(PETScVector,x)
+    B = convert(PETScVector,b)
+    owners = (x,b)
+
+    A = typeof(owners)
+    new{TM,A}(ns,X,B,owners)
+  end
+end
+
+function Algebra.solve!(x::AbstractVector,ns::CachedPETScNS,b::AbstractVector)
+  @assert x === ns.owners[1]
+  @assert b === ns.owners[2]
+  solve!(ns.X,ns.ns,ns.B)
+  consistent!(x)
+  return x
+end
diff --git a/src/NonlinearSolvers/NewtonRaphsonSolver.jl b/src/NonlinearSolvers/NewtonRaphsonSolver.jl
new file mode 100644
index 00000000..86bd8bf9
--- /dev/null
+++ b/src/NonlinearSolvers/NewtonRaphsonSolver.jl
@@ -0,0 +1,70 @@
+
+# TODO: This should be called NewtonRaphsonSolver, but it would clash with Gridap. 
+struct NewtonSolver <: Algebra.NonlinearSolver
+  ls ::Algebra.LinearSolver
+  log::ConvergenceLog{Float64}
+end
+
+function NewtonSolver(ls;maxiter=100,atol=1e-12,rtol=1.e-6,verbose=0,name="Newton-Raphson")
+  tols = SolverTolerances{Float64}(;maxiter=maxiter,atol=atol,rtol=rtol)
+  log  = ConvergenceLog(name,tols;verbose=verbose)
+  return NewtonSolver(ls,log)
+end
+
+struct NewtonCache
+  A::AbstractMatrix
+  b::AbstractVector
+  dx::AbstractVector
+  ns::NumericalSetup
+end
+
+function Algebra.solve!(x::AbstractVector,nls::NewtonSolver,op::NonlinearOperator,cache::Nothing)
+  b  = residual(op, x)
+  A  = jacobian(op, x)
+  dx = allocate_in_domain(A); fill!(dx,zero(eltype(dx)))
+  ss = symbolic_setup(nls.ls, A)
+  ns = numerical_setup(ss,A,x)
+  _solve_nr!(x,A,b,dx,ns,nls,op)
+  return NewtonCache(A,b,dx,ns)
+end
+
+function Algebra.solve!(x::AbstractVector,nls::NewtonSolver,op::NonlinearOperator,cache::NewtonCache)
+  A,b,dx,ns = cache.A, cache.b, cache.dx, cache.ns
+  residual!(b, op, x)
+  jacobian!(A, op, x)
+  numerical_setup!(ns,A,x)
+  _solve_nr!(x,A,b,dx,ns,nls,op)
+  return cache
+end
+
+function _solve_nr!(x,A,b,dx,ns,nls,op)
+  log = nls.log
+
+  # Check for convergence on the initial residual
+  res = norm(b)
+  done = init!(log,res)
+
+  # Newton-like iterations
+  while !done
+
+    # Solve linearized problem
+    rmul!(b,-1)
+    solve!(dx,ns,b)
+    x .+= dx
+
+    # Check convergence for the current residual
+    residual!(b, op, x)
+    res  = norm(b)
+    done = update!(log,res)
+
+    if !done
+      # Update jacobian and solver
+      jacobian!(A, op, x)
+      numerical_setup!(ns,A,x)
+    end
+
+  end
+
+  finalize!(log,res)
+  return x
+end
diff --git a/src/NonlinearSolvers/NonlinearSolvers.jl b/src/NonlinearSolvers/NonlinearSolvers.jl
new file mode 100644
index 00000000..235769a7
--- /dev/null
+++ b/src/NonlinearSolvers/NonlinearSolvers.jl
@@ -0,0 +1,20 @@
+module NonlinearSolvers
+  using LinearAlgebra
+  using SparseArrays
+  using SparseMatricesCSR
+  using BlockArrays
+  using IterativeSolvers
+
+  using Gridap
+  using Gridap.Helpers, Gridap.Algebra, Gridap.CellData, Gridap.Arrays, Gridap.FESpaces, Gridap.MultiField
+  using PartitionedArrays
+  using GridapDistributed
+
+  using GridapSolvers.SolverInterfaces
+  using GridapSolvers.MultilevelTools
+  using GridapSolvers.SolverInterfaces
+
+  include("NewtonRaphsonSolver.jl")
+  export NewtonSolver
+
+end
\ No newline at end of file
diff --git a/src/PatchBasedSmoothers/seq/PatchTriangulations.jl b/src/PatchBasedSmoothers/seq/PatchTriangulations.jl
index b16abfd3..fead9f41 100644
--- a/src/PatchBasedSmoothers/seq/PatchTriangulations.jl
+++ b/src/PatchBasedSmoothers/seq/PatchTriangulations.jl
@@ -44,12 +44,12 @@ function Gridap.Geometry.Triangulation(PD::PatchDecomposition)
   return PatchTriangulation(trian,PD,patch_cells,nothing,nothing)
 end
 
-function Gridap.Geometry.BoundaryTriangulation(PD::PatchDecomposition{Dc}) where Dc
-  Df       = Dc -1 
+function Gridap.Geometry.BoundaryTriangulation(PD::PatchDecomposition{Dc};tags="boundary") where Dc
+  Df       = Dc-1 
   model    = PD.model
   labeling = get_face_labeling(model)
 
-  is_boundary = get_face_mask(labeling,["boundary"],Df)
+  is_boundary = get_face_mask(labeling,tags,Df)
   patch_faces = get_patch_faces(PD,Df,is_boundary)
   pfaces_to_pcells = get_pfaces_to_pcells(PD,Df,patch_faces)
 
diff --git a/src/SolverInterfaces/ConvergenceLogs.jl b/src/SolverInterfaces/ConvergenceLogs.jl
index 4e29888c..ed877fdd 100644
--- a/src/SolverInterfaces/ConvergenceLogs.jl
+++ b/src/SolverInterfaces/ConvergenceLogs.jl
@@ -5,6 +5,8 @@
   SOLVER_VERBOSE_HIGH = 2
 end
 
+SolverVerboseLevel(verbose::Bool) = (verbose ? SOLVER_VERBOSE_HIGH : SOLVER_VERBOSE_NONE)
+
 mutable struct ConvergenceLog{T<:Real}
   name      :: String
   tols      :: SolverTolerances{T}
@@ -22,8 +24,7 @@ function ConvergenceLog(name :: String,
                         depth   = 0
                         ) where T
   residuals = Vector{T}(undef,tols.maxiter+1)
-  verbose = (isa(verbose,Bool) && verbose) ? SOLVER_VERBOSE_HIGH : verbose
-  verbose = isa(verbose,SolverVerboseLevel) ? verbose : SolverVerboseLevel(verbose)
+  verbose = SolverVerboseLevel(verbose)
   if nested
     depth += 1
   end
@@ -82,6 +83,12 @@ function finalize!(log::ConvergenceLog{T},r::T) where T
   return flag
 end
 
+function print_message(log::ConvergenceLog{T},msg::String) where T
+  if log.verbose > SOLVER_VERBOSE_LOW
+    println(get_tabulation(log),msg)
+  end
+end
+
 function Base.show(io::IO,k::MIME"text/plain",log::ConvergenceLog)
   println(io,"ConvergenceLog[$(log.name)]")
   println(io," > tols: $(summary(log.tols))")
diff --git a/src/SolverInterfaces/GridapExtras.jl b/src/SolverInterfaces/GridapExtras.jl
index a6d5d24f..8d867c55 100644
--- a/src/SolverInterfaces/GridapExtras.jl
+++ b/src/SolverInterfaces/GridapExtras.jl
@@ -1,11 +1,14 @@
 
 # LinearSolvers that depend on the non-linear solution
-"""
-function Gridap.Algebra.numerical_setup!(ns::Gridap.Algebra.LinearSolver,A::AbstractMatrix,x::AbstractVector)
-  numerical_setup!(ns,A)
+
+function Gridap.Algebra.symbolic_setup(ns::Gridap.Algebra.LinearSolver,A::AbstractMatrix,x::AbstractVector)
+  symbolic_setup(ns,A)
 end
 
-function allocate_solver_caches(ns::Gridap.Algebra.LinearSolver,args...;kwargs...)
-  @abstractmethod
+function Gridap.Algebra.numerical_setup(ns::Gridap.Algebra.SymbolicSetup,A::AbstractMatrix,x::AbstractVector)
+  numerical_setup(ns,A)
+end
+
+function Gridap.Algebra.numerical_setup!(ns::Gridap.Algebra.NumericalSetup,A::AbstractMatrix,x::AbstractVector)
+  numerical_setup!(ns,A)
 end
-"""
\ No newline at end of file
diff --git a/src/SolverInterfaces/SolverInterfaces.jl b/src/SolverInterfaces/SolverInterfaces.jl
index 812bf868..24563838 100644
--- a/src/SolverInterfaces/SolverInterfaces.jl
+++ b/src/SolverInterfaces/SolverInterfaces.jl
@@ -14,7 +14,7 @@ include("SolverInfos.jl")
 
 export SolverVerboseLevel, SolverConvergenceFlag
 export SolverTolerances, get_solver_tolerances, set_solver_tolerances!
-export ConvergenceLog, init!, update!, finalize!, reset!
+export ConvergenceLog, init!, update!, finalize!, reset!, print_message
 
 export SolverInfo
 
diff --git a/src/SolverInterfaces/SolverTolerances.jl b/src/SolverInterfaces/SolverTolerances.jl
index 28668d67..1112b5fc 100644
--- a/src/SolverInterfaces/SolverTolerances.jl
+++ b/src/SolverInterfaces/SolverTolerances.jl
@@ -28,22 +28,22 @@ function set_solver_tolerances!(a::SolverTolerances{T};
                                 rtol   = T(1.e-5),
                                 dtol   = T(Inf)) where T
   a.maxiter = maxiter
-  a.atol   = atol
-  a.rtol   = rtol
-  a.dtol   = dtol
+  a.atol = atol
+  a.rtol = rtol
+  a.dtol = dtol
   return a
 end
 
 function finished_flag(tols::SolverTolerances,niter,e_a,e_r)
-  if !finished(tols,niter,e_r,e_a)
+  if !finished(tols,niter,e_a,e_r)
     @warn "finished_flag() called with unfinished solver!"
   end
-  if niter > tols.maxiter
-    return SOLVER_DIVERGED_MAXITER
-  elseif e_r < tols.rtol
+  if e_r < tols.rtol
     return SOLVER_CONVERGED_RTOL
   elseif e_a < tols.atol
     return SOLVER_CONVERGED_ATOL
+  elseif niter >= tols.maxiter
+    return SOLVER_DIVERGED_MAXITER
   else
     return SOLVER_DIVERGED_BREAKDOWN
   end
diff --git a/test/BlockSolvers/BlockDiagonalSolversTests.jl b/test/BlockSolvers/BlockDiagonalSolversTests.jl
new file mode 100644
index 00000000..ad1efac6
--- /dev/null
+++ b/test/BlockSolvers/BlockDiagonalSolversTests.jl
@@ -0,0 +1,56 @@
+using BlockArrays, LinearAlgebra
+using Gridap, Gridap.MultiField, Gridap.Algebra
+using PartitionedArrays, GridapDistributed
+using GridapSolvers
+
+np = (2,2)
+ranks = with_debug() do distribute
+  distribute(LinearIndices((prod(np),)))
+end
+
+model = CartesianDiscreteModel(ranks,np,(0,1,0,1),(8,8))
+
+reffe = ReferenceFE(lagrangian,Float64,1)
+V = FESpace(model,reffe)
+
+mfs = BlockMultiFieldStyle()
+Y = MultiFieldFESpace([V,V];style=mfs)
+
+Ω  = Triangulation(model)
+dΩ = Measure(Ω,4)
+
+sol(x) = sum(x)
+a((u1,u2),(v1,v2)) = ∫(u1⋅v1 + u2⋅v2)*dΩ
+l((v1,v2)) = ∫(sol⋅v1 - sol⋅v2)*dΩ
+
+op = AffineFEOperator(a,l,Y,Y)
+A, b = get_matrix(op), get_vector(op);
+
+
+# 1) From system blocks
+s1  = BlockDiagonalSolver([LUSolver(),LUSolver()])
+ss1 = symbolic_setup(s1,A)
+ns1 = numerical_setup(ss1,A)
+numerical_setup!(ns1,A)
+
+x1  = allocate_in_domain(A); fill!(x1,0.0)
+solve!(x1,ns1,b)
+
+# 2) From matrix blocks
+s2  = BlockDiagonalSolver([A[Block(1,1)],A[Block(2,2)]],[LUSolver(),LUSolver()])
+ss2 = symbolic_setup(s2,A)
+ns2 = numerical_setup(ss2,A)
+numerical_setup!(ns2,A)
+
+x2  = allocate_in_domain(A); fill!(x2,0.0)
+solve!(x2,ns2,b)
+
+# 3) From weakform blocks
+aii = (u,v) -> ∫(u⋅v)*dΩ
+s3  = BlockDiagonalSolver([aii,aii],[V,V],[V,V],[LUSolver(),LUSolver()])
+ss3 = symbolic_setup(s3,A)
+ns3 = numerical_setup(ss3,A)
+numerical_setup!(ns3,A)
+
+x3  = allocate_in_domain(A); fill!(x3,0.0)
+solve!(x3,ns3,b)
diff --git a/test/BlockSolvers/BlockFEOperatorsTests.jl b/test/BlockSolvers/BlockFEOperatorsTests.jl
new file mode 100644
index 00000000..71b94d68
--- /dev/null
+++ b/test/BlockSolvers/BlockFEOperatorsTests.jl
@@ -0,0 +1,57 @@
+using Test
+using BlockArrays, LinearAlgebra
+using Gridap, Gridap.MultiField, Gridap.Algebra
+using PartitionedArrays, GridapDistributed
+using GridapSolvers, GridapSolvers.BlockSolvers
+
+function same_block_array(A,B)
+  map(blocks(A),blocks(B)) do A, B
+    t = map(partition(A),partition(B)) do A, B
+      A ≈ B
+    end
+    reduce(&,t)
+  end |> all
+end
+
+np = (2,2)
+ranks = with_debug() do distribute
+  distribute(LinearIndices((prod(np),)))
+end
+
+model = CartesianDiscreteModel(ranks,np,(0,1,0,1),(8,8))
+
+reffe = ReferenceFE(lagrangian,Float64,1)
+V = FESpace(model,reffe)
+
+mfs = BlockMultiFieldStyle()
+Y = MultiFieldFESpace([V,V];style=mfs)
+
+Ω  = Triangulation(model)
+dΩ = Measure(Ω,4)
+
+u0 = zero(Y)
+sol(x) = sum(x)
+
+# Reference operator
+a_ref((u1,u2),(v1,v2)) = ∫(u1⋅v1 + u2⋅v2)*dΩ
+l_ref((v1,v2)) = ∫(sol⋅v1 + sol⋅v2)*dΩ
+res_ref(u,v) = a_ref(u,v) - l_ref(v)
+jac_ref(u,du,dv) = a_ref(du,dv)
+op_ref = FEOperator(res_ref,jac_ref,Y,Y)
+A_ref = jacobian(op_ref,u0)
+b_ref = residual(op_ref,u0)
+
+# Block operator
+a(u,v) = ∫(u⋅v)*dΩ
+l(v)   = ∫(sol⋅v)*dΩ
+res(u,v) = a(u,v) - l(v)
+jac(u,du,dv) = a(du,dv)
+
+res_blocks = collect(reshape([res,missing,missing,res],(2,2)))
+jac_blocks = collect(reshape([jac,missing,missing,jac],(2,2)))
+op = BlockFEOperator(res_blocks,jac_blocks,Y,Y)
+A  = jacobian(op,u0)
+b  = residual(op,u0)
+
+@test same_block_array(A,A_ref)
+@test same_block_array(b,b_ref)
diff --git a/test/BlockSolvers/BlockTriangularSolversTests.jl b/test/BlockSolvers/BlockTriangularSolversTests.jl
new file mode 100644
index 00000000..2e84624a
--- /dev/null
+++ b/test/BlockSolvers/BlockTriangularSolversTests.jl
@@ -0,0 +1,46 @@
+
+using BlockArrays, LinearAlgebra
+using Gridap, Gridap.MultiField, Gridap.Algebra
+using PartitionedArrays, GridapDistributed
+using GridapSolvers, GridapSolvers.BlockSolvers
+
+np = (2,2)
+ranks = with_debug() do distribute
+  distribute(LinearIndices((prod(np),)))
+end
+
+model = CartesianDiscreteModel(ranks,np,(0,1,0,1),(8,8))
+
+reffe = ReferenceFE(lagrangian,Float64,1)
+V = FESpace(model,reffe)
+
+mfs = BlockMultiFieldStyle()
+Y = MultiFieldFESpace([V,V];style=mfs)
+
+Ω  = Triangulation(model)
+dΩ = Measure(Ω,4)
+
+sol(x) = sum(x)
+a((u1,u2),(v1,v2)) = ∫(u1⋅v1 + u2⋅v2 + u1⋅v2 - u2⋅v1)*dΩ
+l((v1,v2)) = ∫(sol⋅v1 - sol⋅v2)*dΩ
+
+op = AffineFEOperator(a,l,Y,Y)
+A, b = get_matrix(op), get_vector(op);
+
+# Upper
+s1  = BlockTriangularSolver([LUSolver(),LUSolver()];half=:upper)
+ss1 = symbolic_setup(s1,A)
+ns1 = numerical_setup(ss1,A)
+numerical_setup!(ns1,A)
+
+x1 = allocate_in_domain(A); fill!(x1,0.0)
+solve!(x1,ns1,b)
+
+# Lower
+s2  = BlockTriangularSolver([LUSolver(),LUSolver()];half=:lower)
+ss2 = symbolic_setup(s2,A)
+ns2 = numerical_setup(ss2,A)
+numerical_setup!(ns2,A)
+
+x2 = allocate_in_domain(A); fill!(x2,0.0)
+solve!(x2,ns2,b)
diff --git a/test/LinearSolvers/GMGTests.jl b/test/LinearSolvers/GMGTests.jl
index 8cfaebf1..bff9ab83 100644
--- a/test/LinearSolvers/GMGTests.jl
+++ b/test/LinearSolvers/GMGTests.jl
@@ -63,7 +63,8 @@ function gmg_driver(t,parts,mh,spaces,qdegree,smoothers,biform,liform,u)
                         verbose=false,
                         mode=:preconditioner)
   
-  solver = CGSolver(gmg;maxiter=100,atol=1e-10,rtol=1.e-6,verbose=i_am_main(parts))
+  solver = CGSolver(gmg;maxiter=20,atol=1e-14,rtol=1.e-6,verbose=i_am_main(parts))
+  #solver = GMRESSolver(5;Pr=gmg,maxiter=20,atol=1e-14,rtol=1.e-6,verbose=i_am_main(parts))
   ns = numerical_setup(symbolic_setup(solver,A),A)
   toc!(t,"GMG setup")
 
diff --git a/test/LinearSolvers/KrylovSolversTests.jl b/test/LinearSolvers/KrylovSolversTests.jl
index f0857fe2..ca114ff7 100644
--- a/test/LinearSolvers/KrylovSolversTests.jl
+++ b/test/LinearSolvers/KrylovSolversTests.jl
@@ -16,7 +16,7 @@ function test_solver(solver,op,Uh,dΩ)
   A, b = get_matrix(op), get_vector(op);
   ns = numerical_setup(symbolic_setup(solver,A),A)
 
-  x = allocate_in_domain(A)
+  x = allocate_in_domain(A); fill!(x,0.0)
   solve!(x,ns,b)
 
   u  = interpolate(sol,Uh)
@@ -69,10 +69,16 @@ function main(distribute,np)
   test_solver(gmres,op,Uh,dΩ)
 
   # GMRES without preconditioner
-  gmres = LinearSolvers.GMRESSolver(40;rtol=1.e-8,verbose=verbose)
+  gmres = LinearSolvers.GMRESSolver(10;rtol=1.e-8,verbose=verbose)
   test_solver(gmres,op,Uh,dΩ)
 
-  fgmres = LinearSolvers.FGMRESSolver(40,P;rtol=1.e-8,verbose=verbose)
+  gmres = LinearSolvers.GMRESSolver(10;restart=true,rtol=1.e-8,verbose=verbose)
+  test_solver(gmres,op,Uh,dΩ)
+
+  fgmres = LinearSolvers.FGMRESSolver(10,P;rtol=1.e-8,verbose=verbose)
+  test_solver(fgmres,op,Uh,dΩ)
+
+  fgmres = LinearSolvers.FGMRESSolver(10,P;restart=true,rtol=1.e-8,verbose=verbose)
   test_solver(fgmres,op,Uh,dΩ)
 
   pcg = LinearSolvers.CGSolver(P;rtol=1.e-8,verbose=verbose)