diff --git a/src/BlockSolvers/BlockTriangularSolvers.jl b/src/BlockSolvers/BlockTriangularSolvers.jl
index a2afb5fd..0fa9b622 100644
--- a/src/BlockSolvers/BlockTriangularSolvers.jl
+++ b/src/BlockSolvers/BlockTriangularSolvers.jl
@@ -91,7 +91,11 @@ function Gridap.Algebra.numerical_setup!(ns::BlockTriangularSolverNS,mat::Abstra
   solver       = ns.solver
   mat_blocks   = blocks(mat)
   block_caches = map(update_block_cache!,ns.block_caches,solver.blocks,mat_blocks)
-  map(numerical_setup!,ns.block_ns,diag(block_caches))
+  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
 
@@ -102,7 +106,11 @@ function Gridap.Algebra.numerical_setup!(ns::BlockTriangularSolverNS,mat::Abstra
   block_caches = map(CartesianIndices(solver.blocks)) do I
     update_block_cache!(ns.block_caches[I],mat_blocks[I],vec_blocks[I[2]])
   end
-  map(numerical_setup!,ns.block_ns,diag(block_caches),vec_blocks)
+  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
 
diff --git a/src/LinearSolvers/Krylov/CGSolvers.jl b/src/LinearSolvers/Krylov/CGSolvers.jl
index 12a69c3a..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)
diff --git a/src/LinearSolvers/Krylov/FGMRESSolvers.jl b/src/LinearSolvers/Krylov/FGMRESSolvers.jl
index 31adf264..c5de93a1 100644
--- a/src/LinearSolvers/Krylov/FGMRESSolvers.jl
+++ b/src/LinearSolvers/Krylov/FGMRESSolvers.jl
@@ -62,6 +62,16 @@ 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)
diff --git a/src/LinearSolvers/Krylov/GMRESSolvers.jl b/src/LinearSolvers/Krylov/GMRESSolvers.jl
index b2e5de7b..14d603e9 100644
--- a/src/LinearSolvers/Krylov/GMRESSolvers.jl
+++ b/src/LinearSolvers/Krylov/GMRESSolvers.jl
@@ -63,6 +63,18 @@ 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)
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/NonlinearSolvers/NewtonRaphsonSolver.jl b/src/NonlinearSolvers/NewtonRaphsonSolver.jl
index 08e7bdac..a2c17471 100644
--- a/src/NonlinearSolvers/NewtonRaphsonSolver.jl
+++ b/src/NonlinearSolvers/NewtonRaphsonSolver.jl
@@ -60,7 +60,7 @@ function _solve_nr!(x,A,b,dx,ns,nls,op)
     if !done
       # Update jacobian and solver
       jacobian!(A, op, x)
-      numerical_setup!(ns,A)
+      numerical_setup!(ns,A,x)
     end
 
   end
diff --git a/src/NonlinearSolvers/NonlinearSolvers.jl b/src/NonlinearSolvers/NonlinearSolvers.jl
index 0ee37902..235769a7 100644
--- a/src/NonlinearSolvers/NonlinearSolvers.jl
+++ b/src/NonlinearSolvers/NonlinearSolvers.jl
@@ -10,6 +10,7 @@ module NonlinearSolvers
   using PartitionedArrays
   using GridapDistributed
 
+  using GridapSolvers.SolverInterfaces
   using GridapSolvers.MultilevelTools
   using GridapSolvers.SolverInterfaces
 
diff --git a/src/SolverInterfaces/GridapExtras.jl b/src/SolverInterfaces/GridapExtras.jl
index a6d5d24f..88a549ee 100644
--- a/src/SolverInterfaces/GridapExtras.jl
+++ b/src/SolverInterfaces/GridapExtras.jl
@@ -1,11 +1,6 @@
 
 # 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)
 end
-
-function allocate_solver_caches(ns::Gridap.Algebra.LinearSolver,args...;kwargs...)
-  @abstractmethod
-end
-"""
\ No newline at end of file