diff --git a/Solver/linearSystem.h b/Solver/linearSystem.h
index ef701944257310de9f61c03006c44e734471864a..34fbc579302aab28e7528eb4d0dd4b51ac55e933 100644
--- a/Solver/linearSystem.h
+++ b/Solver/linearSystem.h
@@ -26,6 +26,7 @@ class linearSystemBase {
   void setParameter (std::string key, std::string value);
   virtual void insertInSparsityPattern(int _row, int _col){};
   virtual double normInfRightHandSide() const = 0;
+  virtual void createMatrix() {};
 };
 
 template <class scalar>
diff --git a/Solver/linearSystemPETSc.h b/Solver/linearSystemPETSc.h
index ad773c6c81182f2af7a7f6afe988777ac54fc80e..d97e326d826231df5906c94d33656087c7f914ce 100644
--- a/Solver/linearSystemPETSc.h
+++ b/Solver/linearSystemPETSc.h
@@ -79,6 +79,7 @@ class linearSystemPETSc : public linearSystem<scalar> {
   std::vector<scalar> getData();
   std::vector<int> getRowPointers();
   std::vector<int> getColumnsIndices();
+  virtual void createMatrix();
 };
 
 class linearSystemPETScBlockDouble : public linearSystem<fullMatrix<double> > {
diff --git a/Solver/linearSystemPETSc.hpp b/Solver/linearSystemPETSc.hpp
index 4505e6ae5e01e944aa219f8acabbd328f66952a3..ca682349e7d1ab4572a868a62f2c0873d13679b3 100644
--- a/Solver/linearSystemPETSc.hpp
+++ b/Solver/linearSystemPETSc.hpp
@@ -131,6 +131,29 @@ void linearSystemPETSc<scalar>::allocate(int nbRows)
   _entriesPreAllocated = false;
 }
 
+template<class scalar>
+void linearSystemPETSc<scalar>::createMatrix(){
+  if (isAllocated())
+    #if (PETSC_VERSION_RELEASE == 0 || ((PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2)))
+    _try(MatDestroy(&_a));
+    #else
+    _try(MatDestroy(_a));
+    #endif
+  _try(MatCreate(_comm, &_a));
+  _try(MatSetSizes(_a, _localSize, _localSize, _globalSize, _globalSize));
+  // override the default options with the ones from the option
+  // database (if any)
+  if (this->_parameters.count("petscOptions"))
+    _try(PetscOptionsInsertString(this->_parameters["petscOptions"].c_str()));
+  if (this->_parameters.count("petscPrefix"))
+    _try(MatAppendOptionsPrefix(_a, this->_parameters["petscPrefix"].c_str()));
+  _try(MatSetFromOptions(_a));
+
+  _entriesPreAllocated = false;
+  _matrixModified = true;
+  _isAllocated = true;
+};
+
 template <class scalar>
 void linearSystemPETSc<scalar>::print()
 {