diff --git a/Solver/linearSystemPETSc.cpp b/Solver/linearSystemPETSc.cpp
index 0bba26d01357dc87f5f43e6cde7cb2d859c1c7d3..461199370999d792171770bbf730a761ffd7fde9 100644
--- a/Solver/linearSystemPETSc.cpp
+++ b/Solver/linearSystemPETSc.cpp
@@ -61,6 +61,7 @@ void linearSystemPETScBlockDouble::addToMatrix(int row, int col,
       modval(jj, ii) = buff;
     }
   #endif
+  _matrixModified=true;
 }
 
 void linearSystemPETScBlockDouble::addToRightHandSide(int row,
@@ -199,14 +200,17 @@ int linearSystemPETScBlockDouble::systemSolve()
 {
   if (!_kspAllocated)
     _kspCreate();
-  if (_parameters["matrix_reuse"] == "same_sparsity")
-    KSPSetOperators(_ksp, _a, _a, SAME_NONZERO_PATTERN);
-  else if (_parameters["matrix_reuse"] == "same_matrix")
+  if (!_matrixModified || _parameters["matrix_reuse"] == "same_matrix")
     KSPSetOperators(_ksp, _a, _a, SAME_PRECONDITIONER);
+  else if (_parameters["matrix_reuse"] == "same_sparsity")
+    KSPSetOperators(_ksp, _a, _a, SAME_NONZERO_PATTERN);
   else
     KSPSetOperators(_ksp, _a, _a, DIFFERENT_NONZERO_PATTERN);
-  MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY);
-  MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY);
+  if (_matrixModified && _parameters["matrix_reuse"]!="same_matrix"){
+    MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY);
+    MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY);
+  }
+  _matrixModified=false;
   VecAssemblyBegin(_b);
   VecAssemblyEnd(_b);
   KSPSolve(_ksp, _b, _x);
@@ -289,6 +293,7 @@ linearSystemPETScBlockDouble::linearSystemPETScBlockDouble(bool sequential)
   _entriesPreAllocated = false;
   _isAllocated = false;
   _kspAllocated = false;
+  _matrixModified=true;
   _sequential = sequential;
 }
 
diff --git a/Solver/linearSystemPETSc.h b/Solver/linearSystemPETSc.h
index 1cd061baff0e4e2d6a5b059e05a7e9a972a234b4..ad773c6c81182f2af7a7f6afe988777ac54fc80e 100644
--- a/Solver/linearSystemPETSc.h
+++ b/Solver/linearSystemPETSc.h
@@ -48,7 +48,7 @@ class linearSystemPETSc : public linearSystem<scalar> {
   protected:
   MPI_Comm _comm;
   int _blockSize; // for block Matrix
-  bool _isAllocated, _kspAllocated, _entriesPreAllocated;
+  bool _isAllocated, _kspAllocated, _entriesPreAllocated, _matrixModified;
   Mat _a;
   Vec _b, _x;
   KSP _ksp;
@@ -82,7 +82,7 @@ class linearSystemPETSc : public linearSystem<scalar> {
 };
 
 class linearSystemPETScBlockDouble : public linearSystem<fullMatrix<double> > {
-  bool _entriesPreAllocated, _isAllocated, _kspAllocated;
+  bool _entriesPreAllocated, _isAllocated, _kspAllocated, _matrixModified;
   sparsityPattern _sparsity;
   Mat _a;
   Vec _b, _x;
diff --git a/Solver/linearSystemPETSc.hpp b/Solver/linearSystemPETSc.hpp
index 655642f0723fd7226b0bdeb95d75bc8a027d0e63..f5216ebdb3a16d1ae4aee8b8f3dfddb7ad9812c1 100644
--- a/Solver/linearSystemPETSc.hpp
+++ b/Solver/linearSystemPETSc.hpp
@@ -40,6 +40,7 @@ linearSystemPETSc<scalar>::linearSystemPETSc(MPI_Comm com)
   _isAllocated = false;
   _blockSize = 0;
   _kspAllocated = false;
+  _matrixModified=true;
 }
 
 template <class scalar>
@@ -212,6 +213,7 @@ void linearSystemPETSc<scalar>::addToMatrix(int row, int col, const scalar &val)
   PetscInt i = row, j = col;
   PetscScalar s = val;
   _try(MatSetValues(_a, 1, &i, 1, &j, &s, ADD_VALUES));
+  _matrixModified=true;
 }
 
 template <class scalar>
@@ -272,14 +274,17 @@ int linearSystemPETSc<scalar>::systemSolve()
 {
   if (!_kspAllocated)
     _kspCreate();
-  if (linearSystem<scalar>::_parameters["matrix_reuse"] == "same_sparsity")
-    _try(KSPSetOperators(_ksp, _a, _a, SAME_NONZERO_PATTERN));
-  else if (linearSystem<scalar>::_parameters["matrix_reuse"] == "same_matrix")
+  if (!_matrixModified || linearSystem<scalar>::_parameters["matrix_reuse"] == "same_matrix")
     _try(KSPSetOperators(_ksp, _a, _a, SAME_PRECONDITIONER));
+  else if (linearSystem<scalar>::_parameters["matrix_reuse"] == "same_sparsity")
+    _try(KSPSetOperators(_ksp, _a, _a, SAME_NONZERO_PATTERN));
   else
     _try(KSPSetOperators(_ksp, _a, _a, DIFFERENT_NONZERO_PATTERN));
-  _try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
-  _try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
+  if (_matrixModified && _parameters["matrix_reuse"]!="same_matrix"){
+    _try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
+    _try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
+  }
+  _matrixModified=false;
   /*MatInfo info;
     MatGetInfo(_a, MAT_LOCAL, &info);
     printf("mallocs %.0f    nz_allocated %.0f    nz_used %.0f    nz_unneeded %.0f\n", info.mallocs, info.nz_allocated, info.nz_used, info.nz_unneeded);*/