Solver: added matMult function

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){};

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

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
 {