PETSc preconditioning not done when matrix is not modified

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

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;

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);*/