From b424a0f225a897a17e99a2103849d65a112ba1e6 Mon Sep 17 00:00:00 2001
From: Sebastien Blaise <sebastien.blaise@uclouvain.be>
Date: Fri, 17 Jul 2015 10:28:32 +0000
Subject: [PATCH] Solver: added matMult function

---
 Solver/linearSystem.h        |  3 +++
 Solver/linearSystemPETSc.h   |  2 ++
 Solver/linearSystemPETSc.hpp | 12 ++++++++++++
 3 files changed, 17 insertions(+)

diff --git a/Solver/linearSystem.h b/Solver/linearSystem.h
index b008e7e37f..2cb4dfa4f7 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 7db6f53ce6..e7e1ddc3a2 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 4db9597b2b..d835b73ecb 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
 {
-- 
GitLab