diff --git a/src/BlockSolvers/BlockDiagonalSolvers.jl b/src/BlockSolvers/BlockDiagonalSolvers.jl
index 9b3a30ed..30787edb 100644
--- a/src/BlockSolvers/BlockDiagonalSolvers.jl
+++ b/src/BlockSolvers/BlockDiagonalSolvers.jl
@@ -1,7 +1,7 @@
 
 struct BlockDiagonalSolver{N,A,B} <: Gridap.Algebra.LinearSolver
-  blocks  :: B
-  solvers :: C
+  blocks  :: A
+  solvers :: B
   function BlockDiagonalSolver(
     blocks  :: AbstractVector{<:SolverBlock},
     solvers :: AbstractVector{<:Gridap.Algebra.LinearSolver}
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
index 82f86a90..38bc28ca 100644
--- a/src/BlockSolvers/BlockSolverInterfaces.jl
+++ b/src/BlockSolvers/BlockSolverInterfaces.jl
@@ -19,6 +19,12 @@ struct BiformBlock <: LinearSolverBlock
   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
@@ -26,6 +32,12 @@ struct TriformBlock <: NonlinearSolverBlock
   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
 
 # Instantiate blocks
@@ -48,7 +60,7 @@ function instantiate_block_cache(block::BiformBlock,mat::AbstractMatrix)
 end
 instantiate_block_cache(block::LinearSystemBlock,mat::AbstractMatrix) = mat
 
-function instantiate_block_cache(block::TriformSolverBlock,mat::AbstractMatrix,x::AbstractVector)
+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)
@@ -64,7 +76,7 @@ function update_block_cache!(cache,block::NonlinearSolverBlock,mat::AbstractMatr
   @abstractmethod
 end
 function update_block_cache!(cache,block::LinearSolverBlock,mat::AbstractMatrix,x::AbstractVector)
-  update_block!(cache,block,mat)
+  update_block_cache!(cache,block,mat)
 end
 
 function update_block_cache!(cache,block::TriformBlock,mat::AbstractMatrix,x::AbstractVector)
diff --git a/src/BlockSolvers/BlockSolvers.jl b/src/BlockSolvers/BlockSolvers.jl
index f3b78d65..2a7d0fa9 100644
--- a/src/BlockSolvers/BlockSolvers.jl
+++ b/src/BlockSolvers/BlockSolvers.jl
@@ -13,10 +13,14 @@ module BlockSolvers
   using GridapSolvers.MultilevelTools
   using GridapSolvers.SolverInterfaces
 
+  include("BlockFEOperators.jl")
+
   include("BlockSolverInterfaces.jl")
   include("BlockDiagonalSolvers.jl")
   include("BlockTriangularSolvers.jl")
 
+  export BlockFEOperator
+
   export BlockDiagonalSolver
   export BlockTriangularSolver
 end
diff --git a/src/GridapSolvers.jl b/src/GridapSolvers.jl
index 12870f53..5316d068 100644
--- a/src/GridapSolvers.jl
+++ b/src/GridapSolvers.jl
@@ -2,11 +2,13 @@ module GridapSolvers
 
   include("SolverInterfaces/SolverInterfaces.jl")
   include("MultilevelTools/MultilevelTools.jl")
+  include("BlockSolvers/BlockSolvers.jl")
   include("LinearSolvers/LinearSolvers.jl")
   include("PatchBasedSmoothers/PatchBasedSmoothers.jl")
 
   using GridapSolvers.SolverInterfaces
   using GridapSolvers.MultilevelTools
+  using GridapSolvers.BlockSolvers
   using GridapSolvers.LinearSolvers
   using GridapSolvers.PatchBasedSmoothers
 
@@ -25,6 +27,9 @@ module GridapSolvers
   export RestrictionOperator, ProlongationOperator
   export setup_transfer_operators
 
+  # BlockSolvers
+  export BlockDiagonalSolver
+
   # LinearSolvers
   export JacobiLinearSolver
   export RichardsonSmoother
@@ -37,7 +42,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/Krylov/FGMRESSolvers.jl b/src/LinearSolvers/Krylov/FGMRESSolvers.jl
index 34da4ae4..31adf264 100644
--- a/src/LinearSolvers/Krylov/FGMRESSolvers.jl
+++ b/src/LinearSolvers/Krylov/FGMRESSolvers.jl
@@ -70,7 +70,6 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::FGMRESNumericalSetup,b::Abs
   V, Z, zl, H, g, c, s = caches
 
   fill!(V[1],zero(eltype(V[1])))
-  fill!(zr,zero(eltype(zr)))
   fill!(zl,zero(eltype(zl)))
 
   # Initial residual
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)