From 55211e9b7627ce14f263887ade9004545ad89b89 Mon Sep 17 00:00:00 2001
From: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date: Wed, 24 Jul 2013 08:24:31 +0000
Subject: [PATCH] petsc : do not use MatAssembled, keep our own flag instead so
that we can control its behaviour on empty matrices
---
Solver/linearSystemPETSc.h | 7 ++++--
Solver/linearSystemPETSc.hpp | 47 +++++++++++++++++++-----------------
2 files changed, 30 insertions(+), 24 deletions(-)
diff --git a/Solver/linearSystemPETSc.h b/Solver/linearSystemPETSc.h
index c38597cd02..ea8579a35b 100644
--- a/Solver/linearSystemPETSc.h
+++ b/Solver/linearSystemPETSc.h
@@ -51,7 +51,7 @@ typedef struct _p_Vec* Vec;
typedef struct _p_KSP* KSP;
#endif
-//support old PETSc version, try to avoid using PETSC_VERSION somewhere else
+//support old PETSc version, try to avoid using PETSC_VERSION in other places
#if PETSC_VERSION_RELEASE != 0 && (PETSC_VERSION_MAJOR < 3 || (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR < 2))
#define KSPDestroy(k) KSPDestroy(*(k))
#define MatDestroy(m) MatDestroy(*(m))
@@ -66,13 +66,16 @@ template <class scalar>
class linearSystemPETSc : public linearSystem<scalar> {
protected:
int _blockSize; // for block Matrix
- bool _isAllocated, _kspAllocated, _entriesPreAllocated, _matrixModified;
+ bool _isAllocated, _kspAllocated, _entriesPreAllocated;
+ bool _matrixChangedSinceLastSolve;
+ bool _valuesNotAssembled; //cannot use MatAssembled since MatAssembled return false for an empty matrix
Mat _a;
Vec _b, _x;
KSP _ksp;
int _localRowStart, _localRowEnd, _localSize, _globalSize;
sparsityPattern _sparsity;
void _kspCreate();
+ void _assembleMatrixIfNeeded();
#ifndef SWIG
MPI_Comm _comm;
#endif
diff --git a/Solver/linearSystemPETSc.hpp b/Solver/linearSystemPETSc.hpp
index 0e42d65d9a..81024fadf3 100644
--- a/Solver/linearSystemPETSc.hpp
+++ b/Solver/linearSystemPETSc.hpp
@@ -36,7 +36,8 @@ linearSystemPETSc<scalar>::linearSystemPETSc(MPI_Comm com)
_isAllocated = false;
_blockSize = 0;
_kspAllocated = false;
- _matrixModified=true;
+ _matrixChangedSinceLastSolve = true;
+ _valuesNotAssembled = false;
}
template <class scalar>
@@ -46,7 +47,8 @@ linearSystemPETSc<scalar>::linearSystemPETSc()
_isAllocated = false;
_blockSize = 0;
_kspAllocated = false;
- _matrixModified=true;
+ _matrixChangedSinceLastSolve = true;
+ _valuesNotAssembled = false;
}
template <class scalar>
@@ -57,6 +59,18 @@ linearSystemPETSc<scalar>::~linearSystemPETSc()
_try(KSPDestroy(&_ksp));
}
+template <class scalar>
+void linearSystemPETSc<scalar>::_assembleMatrixIfNeeded()
+{
+ // avoid to assemble the matrix when not needed since it destroy preallocation (e.g. in zeroMatrix)
+ if (_valuesNotAssembled) {
+ _try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
+ _try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
+ _matrixChangedSinceLastSolve = true;
+ _valuesNotAssembled = false;
+ }
+}
+
template <class scalar>
void linearSystemPETSc<scalar>::insertInSparsityPattern (int i, int j)
{
@@ -159,8 +173,7 @@ void linearSystemPETSc<scalar>::allocate(int nbRows)
template <class scalar>
void linearSystemPETSc<scalar>::print()
{
- _try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
- _try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
+ _assembleMatrixIfNeeded();
_try(VecAssemblyBegin(_b));
_try(VecAssemblyEnd(_b));
@@ -241,7 +254,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;
+ _valuesNotAssembled = true;
}
template <class scalar>
@@ -280,14 +293,8 @@ void linearSystemPETSc<scalar>::zeroMatrix()
}
}
if (_isAllocated && _entriesPreAllocated) {
- PetscBool assembled;
- _try(MatAssembled(_a, &assembled));
- if (!assembled) {
- _try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
- _try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
- }
+ _assembleMatrixIfNeeded();
_try(MatZeroEntries(_a));
- _matrixModified = true;
}
}
@@ -317,15 +324,14 @@ int linearSystemPETSc<scalar>::systemSolve()
{
if (!_kspAllocated)
_kspCreate();
- if (!_matrixModified || linearSystem<scalar>::_parameters["matrix_reuse"] == "same_matrix")
+ _assembleMatrixIfNeeded();
+ if (!_matrixChangedSinceLastSolve || 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));
- _matrixModified=false;
+ _matrixChangedSinceLastSolve = 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);*/
@@ -342,8 +348,7 @@ int linearSystemPETSc<scalar>::systemSolve()
/*template <class scalar>
std::vector<scalar> linearSystemPETSc<scalar>::getData()
{
- _try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
- _try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
+ _assembleMatrixIfNeeded();
PetscScalar *v;
_try(MatGetArray(_a,&v));
MatInfo info;
@@ -364,8 +369,7 @@ std::vector<scalar> linearSystemPETSc<scalar>::getData()
template <class scalar>
std::vector<int> linearSystemPETSc<scalar>::getRowPointers()
{
- _try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
- _try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
+ _assembleMatrixIfNeeded();
const PetscInt *rows;
const PetscInt *columns;
PetscInt n;
@@ -381,8 +385,7 @@ std::vector<int> linearSystemPETSc<scalar>::getRowPointers()
template <class scalar>
std::vector<int> linearSystemPETSc<scalar>::getColumnsIndices()
{
- _try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
- _try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
+ _assembleMatrixIfNeeded();
const PetscInt *rows;
const PetscInt *columns;
PetscInt n;
--
GitLab