diff --git a/Solver/linearSystem.h b/Solver/linearSystem.h
index b008e7e37f577546d7f8f9159e2788045beb3951..2cb4dfa4f7ccbf42fa3c1f8d8dd4deb6c75a891e 100644
--- a/Solver/linearSystem.h
+++ b/Solver/linearSystem.h
@@ -24,6 +24,9 @@ class linearSystemBase {
   virtual void zeroRightHandSide() = 0;
   virtual void zeroSolution() = 0;
   virtual int systemSolve() = 0;
+  // x = A*b
+  virtual int matMult() { return 0; }
+
   void setParameter (std::string key, std::string value);
   std::string getParameter(std::string key) const;
   virtual void insertInSparsityPattern(int _row, int _col){};
diff --git a/Solver/linearSystemPETSc.h b/Solver/linearSystemPETSc.h
index 7db6f53ce630dd7233f41b15e4c421218f026948..e7e1ddc3a24f2150935e1ef39904b315d12440ca 100644
--- a/Solver/linearSystemPETSc.h
+++ b/Solver/linearSystemPETSc.h
@@ -115,6 +115,7 @@ class linearSystemPETSc : public linearSystem<scalar> {
   void printMatlab(const char *filename) const;
   virtual int systemSolve();
   Mat getMatrix(){ return _a; }
+  virtual int matMult();
   //std::vector<scalar> getData();
   //std::vector<int> getRowPointers();
   //std::vector<int> getColumnsIndices();
@@ -149,6 +150,7 @@ class linearSystemPETSc : public linearSystem<scalar> {
   void printMatlab(const char *filename) const{};
   virtual int systemSolve() { return 0; }
   double normInfRightHandSide() const{return 0;}
+  virtual int matMult() { return 0; }
 };
 #endif
 #endif
diff --git a/Solver/linearSystemPETSc.hpp b/Solver/linearSystemPETSc.hpp
index 4db9597b2bc2f4725df64539528a5ba0e7e4a742..d835b73ecbd5abbc571a08a1f73558c5159687c3 100644
--- a/Solver/linearSystemPETSc.hpp
+++ b/Solver/linearSystemPETSc.hpp
@@ -393,6 +393,18 @@ int linearSystemPETSc<scalar>::systemSolve()
   return 1;
 }
 
+template <class scalar>
+int linearSystemPETSc<scalar>::matMult()
+{
+  _try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
+  _try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
+  _try(VecAssemblyBegin(_b));
+  _try(VecAssemblyEnd(_b));
+  _try(MatMult(_a,_b,_x));
+  return 1;
+}
+
+
 template <class scalar>
 void linearSystemPETSc<scalar>::printMatlab(const char *filename) const
 {