diff --git a/src/PETScArrays.jl b/src/PETScArrays.jl
index 1f455ba..6fb6e53 100644
--- a/src/PETScArrays.jl
+++ b/src/PETScArrays.jl
@@ -306,6 +306,13 @@ function Base.copy(a::PETScMatrix)
   Init(v)
 end
 
+function Base.copy!(a::PETScMatrix,b::PETScMatrix)
+  if a !== b
+    @check_error_code PETSC.MatCopy(b.mat[],a.mat[],PETSC.SAME_NONZERO_PATTERN)
+  end
+  a
+end
+
 function Base.copy!(a::PETScMatrix,b::AbstractMatrix)
   _copy!(a.mat[],b)
 end
@@ -363,9 +370,6 @@ function _copy!(petscmat::Mat,mat::AbstractSparseMatrix)
   @check_error_code PETSC.MatAssemblyEnd(petscmat, PETSC.MAT_FINAL_ASSEMBLY)
 end
 
-
-
-
 function Base.convert(::Type{PETScMatrix},a::PETScMatrix)
   a
 end
diff --git a/src/PETScLinearSolvers.jl b/src/PETScLinearSolvers.jl
index 291ad40..76e18e3 100644
--- a/src/PETScLinearSolvers.jl
+++ b/src/PETScLinearSolvers.jl
@@ -98,10 +98,37 @@ function Algebra.solve!(x::PVector,ns::PETScLinearSolverNS,b::PVector)
   x
 end
 
+# NOTE: 
+# We previously threw away PETSc's matrix `ns.B`, and re-set the KSP object (commented code). 
+# This is not only unnecessary, but can also cause some issues with MUMPS. 
+# I am not completely sure why, but here are some notes on the issue: 
+#   - It is matrix-dependent, and only happens in parallel (nprocs > 1). 
+#   - It has to do with the re-use of the symmetric permutation created by MUMPS to 
+#     find the pivots.
+# I think it probably re-orders the matrix internally, and does not re-order it again when we swap it
+# using `KSPSetOperators`. So when accessing the new (non-permuted) matrix using the old permutation, 
+# it throws a stack overflow error.
+# In fact, when updating the SNES setups in the nonlinear solvers, we do not re-set the matrix,
+# but use `copy!` instead.
+#function Algebra.numerical_setup!(ns::PETScLinearSolverNS,A::AbstractMatrix)
+#  # ns.A = A
+#  # ns.B = convert(PETScMatrix,A)
+#  # @check_error_code PETSC.KSPSetOperators(ns.ksp[],ns.B.mat[],ns.B.mat[])
+#  @assert nnz(ns.A) == nnz(A) # This is weak, but it might catch some errors
+#  copy!(ns.B,A)
+#  @check_error_code PETSC.KSPSetUp(ns.ksp[])
+#  ns
+#end
+
 function Algebra.numerical_setup!(ns::PETScLinearSolverNS,A::AbstractMatrix)
-  ns.A = A
-  ns.B = convert(PETScMatrix,A)
-  @check_error_code PETSC.KSPSetOperators(ns.ksp[],ns.B.mat[],ns.B.mat[])
-  @check_error_code PETSC.KSPSetUp(ns.ksp[])
+  if ns.A === A
+    copy!(ns.B,A)
+    @check_error_code PETSC.KSPSetUp(ns.ksp[])
+  else
+    ns.A = A
+    ns.B = convert(PETScMatrix,A)
+    @check_error_code PETSC.KSPSetOperators(ns.ksp[],ns.B.mat[],ns.B.mat[])
+    @check_error_code PETSC.KSPSetUp(ns.ksp[])
+  end
   ns
-end
+end
\ No newline at end of file