From 148d583e12a2b76738bdbe123da8fdd9175a3012 Mon Sep 17 00:00:00 2001
From: Sebastien Blaise <sebastien.blaise@uclouvain.be>
Date: Tue, 12 Jun 2012 18:23:10 +0000
Subject: [PATCH] PETSc preconditioning not done when matrix is not modified

---
 Solver/linearSystemPETSc.cpp | 15 ++++++++++-----
 Solver/linearSystemPETSc.h   |  4 ++--
 Solver/linearSystemPETSc.hpp | 15 ++++++++++-----
 3 files changed, 22 insertions(+), 12 deletions(-)

diff --git a/Solver/linearSystemPETSc.cpp b/Solver/linearSystemPETSc.cpp
index 0bba26d013..4611993709 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 1cd061baff..ad773c6c81 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 655642f072..f5216ebdb3 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);*/
-- 
GitLab