Skip to content
Snippets Groups Projects
Commit 999753cd authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

linearSystemPETSc : implement linearSystemPETScBlockDouble as a template 

spetialization of linearSystemPETSc to avoid code duplication
parent 55211e9b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Solver/linearSystemPETSc.cpp
+ 34
− 246
View file @ 999753cd
......@@ -18,28 +18,18 @@ template class linearSystemPETSc<double>;
template class linearSystemPETSc<std::complex<double> >;
#endif
void linearSystemPETScBlockDouble::_kspCreate()
template<>
int linearSystemPETSc<fullMatrix<double> >::_getBlockSizeFromParameters() const
{
KSPCreate(_sequential ? PETSC_COMM_SELF : PETSC_COMM_WORLD, &_ksp);
if (this->_parameters.count("petscPrefix"))
KSPAppendOptionsPrefix(_ksp, this->_parameters["petscPrefix"].c_str());
KSPSetFromOptions(_ksp);
_kspAllocated = true;
if ( _parameters.find("blockSize") == _parameters.end())
Msg::Error ("'blockSize' parameters must be set for linearSystemPETScBlock");
int blockSize = strtol(_parameters.find("blockSize")->second.c_str(), NULL, 10);
return blockSize;
}
linearSystemPETScBlockDouble::~linearSystemPETScBlockDouble()
{
if (_isAllocated) {
MatDestroy(&_a);
VecDestroy(&_b);
VecDestroy(&_x);
}
if (_kspAllocated) {
KSPDestroy(&_ksp);
}
}
void linearSystemPETScBlockDouble::addToMatrix(int row, int col,
template<>
void linearSystemPETSc<fullMatrix<double> >::addToMatrix(int row, int col,
const fullMatrix<double> &val)
{
if (!_entriesPreAllocated)
......@@ -73,30 +63,30 @@ void linearSystemPETScBlockDouble::addToMatrix(int row, int col,
modval(jj, ii) = buff;
}
#endif
_matrixModified=true;
_valuesNotAssembled = true;
}
void linearSystemPETScBlockDouble::addToRightHandSide(int row,
template<>
void linearSystemPETSc<fullMatrix<double> >::addToRightHandSide(int row,
const fullMatrix<double> &val)
{
for (int ii = 0; ii < _blockSize; ii++) {
PetscInt i = row * _blockSize + ii;
int blockSize;
_try(MatGetBlockSize(_a, &blockSize));
for (int ii = 0; ii < blockSize; ii++) {
PetscInt i = row * blockSize + ii;
PetscScalar v = val(ii, 0);
VecSetValues(_b, 1, &i, &v, ADD_VALUES);
}
}
void linearSystemPETScBlockDouble::getFromMatrix(int row, int col,
fullMatrix<double> &val ) const
{
Msg::Error("getFromMatrix not implemented for PETSc");
}
void linearSystemPETScBlockDouble::getFromRightHandSide(int row,
template<>
void linearSystemPETSc<fullMatrix<double > >::getFromRightHandSide(int row,
fullMatrix<double> &val) const
{
for (int i = 0; i < _blockSize; i++) {
int ii = row*_blockSize +i;
int blockSize;
_try(MatGetBlockSize(_a, &blockSize));
for (int i = 0; i < blockSize; i++) {
int ii = row*blockSize +i;
#ifdef PETSC_USE_COMPLEX
PetscScalar s;
VecGetValues ( _b, 1, &ii, &s);
......@@ -107,19 +97,25 @@ void linearSystemPETScBlockDouble::getFromRightHandSide(int row,
}
}
void linearSystemPETScBlockDouble::addToSolution(int row, const fullMatrix<double> &val)
template<>
void linearSystemPETSc<fullMatrix<double> >::addToSolution(int row, const fullMatrix<double> &val)
{
for (int ii = 0; ii < _blockSize; ii++) {
PetscInt i = row * _blockSize + ii;
int blockSize;
_try(MatGetBlockSize(_a, &blockSize));
for (int ii = 0; ii < blockSize; ii++) {
PetscInt i = row * blockSize + ii;
PetscScalar v = val(ii, 0);
VecSetValues(_x, 1, &i, &v, ADD_VALUES);
}
}
void linearSystemPETScBlockDouble::getFromSolution(int row, fullMatrix<double> &val) const
template<>
void linearSystemPETSc<fullMatrix<double> >::getFromSolution(int row, fullMatrix<double> &val) const
{
for (int i = 0; i < _blockSize; i++) {
int ii = row*_blockSize +i;
int blockSize;
_try(MatGetBlockSize(_a, &blockSize));
for (int i = 0; i < blockSize; i++) {
int ii = row*blockSize +i;
#ifdef PETSC_USE_COMPLEX
PetscScalar s;
VecGetValues ( _x, 1, &ii, &s);
......@@ -130,213 +126,5 @@ void linearSystemPETScBlockDouble::getFromSolution(int row, fullMatrix<double> &
}
}
void linearSystemPETScBlockDouble::allocate(int nbRows)
{
MPI_Comm comm = _sequential ? PETSC_COMM_SELF: PETSC_COMM_WORLD;
if (this->_parameters.count("petscOptions"))
PetscOptionsInsertString(this->_parameters["petscOptions"].c_str());
_blockSize = strtol (_parameters["blockSize"].c_str(), NULL, 10);
if (_blockSize == 0)
Msg::Error ("'blockSize' parameters must be set for linearSystemPETScBlock");
clear();
MatCreate(comm, &_a);
MatSetSizes(_a,nbRows * _blockSize, nbRows * _blockSize, PETSC_DETERMINE, PETSC_DETERMINE);
if (Msg::GetCommSize() > 1 && !_sequential) {
MatSetType(_a, MATMPIBAIJ);
}
else {
MatSetType(_a, MATSEQBAIJ);
}
if (_parameters.count("petscPrefix"))
MatAppendOptionsPrefix(_a, _parameters["petscPrefix"].c_str());
MatSetFromOptions(_a);
//since PETSc 3.3 GetOwnershipRange and MatGetSize() cannot be called before SetPreallocation
_localSize = nbRows;
_localRowStart = 0;
_localRowEnd = nbRows;
_globalSize = _localSize;
#ifdef HAVE_MPI
if (!_sequential) {
_localRowStart = 0;
if (Msg::GetCommRank() != 0) {
MPI_Status status;
MPI_Recv((void*)&_localRowStart, 1, MPI_INT, Msg::GetCommRank() - 1, 1, MPI_COMM_WORLD, &status);
}
_localRowEnd = _localRowStart + nbRows;
if (Msg::GetCommRank() != Msg::GetCommSize() - 1) {
MPI_Send((void*)&_localRowEnd, 1, MPI_INT, Msg::GetCommRank() + 1, 1, MPI_COMM_WORLD);
}
MPI_Allreduce((void*)&_localSize, (void*)&_globalSize, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
}
#endif
// override the default options with the ones from the option
// database (if any)
VecCreate(comm, &_x);
VecSetSizes(_x, nbRows * _blockSize, PETSC_DETERMINE);
// override the default options with the ones from the option
// database (if any)
if (_parameters.count("petscPrefix"))
VecAppendOptionsPrefix(_x, _parameters["petscPrefix"].c_str());
VecSetFromOptions(_x);
VecDuplicate(_x, &_b);
_isAllocated = true;
}
bool linearSystemPETScBlockDouble::isAllocated() const
{
return _isAllocated;
}
void linearSystemPETScBlockDouble::clear()
{
if(_isAllocated){
MatDestroy(&_a);
VecDestroy(&_x);
VecDestroy(&_b);
}
_isAllocated = false;
}
void linearSystemPETScBlockDouble::print()
{
_try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
_try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
_try(VecAssemblyBegin(_b));
_try(VecAssemblyEnd(_b));
if(Msg::GetCommRank()==0)
printf("a :\n");
MatView(_a, PETSC_VIEWER_STDOUT_WORLD);
if(Msg::GetCommRank()==0)
printf("b :\n");
VecView(_b, PETSC_VIEWER_STDOUT_WORLD);
if(Msg::GetCommRank()==0)
printf("x :\n");
VecView(_x, PETSC_VIEWER_STDOUT_WORLD);
}
int linearSystemPETScBlockDouble::systemSolve()
{
if (!_kspAllocated)
_kspCreate();
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);
_matrixModified=false;
VecAssemblyBegin(_b);
VecAssemblyEnd(_b);
KSPSolve(_ksp, _b, _x);
return 1;
}
void linearSystemPETScBlockDouble::insertInSparsityPattern (int i, int j)
{
i -= _localRowStart;
if (i<0 || i>= _localSize) return;
_sparsity.insertEntry (i,j);
}
void linearSystemPETScBlockDouble::preAllocateEntries()
{
if (_entriesPreAllocated) return;
if (!_isAllocated) Msg::Fatal("system must be allocated first");
if (_sparsity.getNbRows() == 0) {
PetscInt prealloc = 300;
PetscBool set;
PetscOptionsGetInt(PETSC_NULL, "-petsc_prealloc", &prealloc, &set);
if (_blockSize == 0) {
MatSeqAIJSetPreallocation(_a, prealloc, PETSC_NULL);
} else {
MatSeqBAIJSetPreallocation(_a, _blockSize, 5, PETSC_NULL);
}
} else {
std::vector<int> nByRowDiag (_localSize), nByRowOffDiag (_localSize);
for (int i = 0; i < _localSize; i++) {
int n;
const int *r = _sparsity.getRow(i, n);
for (int j = 0; j < n; j++) {
if (r[j] >= _localRowStart && r[j] < _localRowEnd)
nByRowDiag[i] ++;
else
nByRowOffDiag[i] ++;
}
}
if (_blockSize == 0) {
MatSeqAIJSetPreallocation(_a, 0, &nByRowDiag[0]);
MatMPIAIJSetPreallocation(_a, 0, &nByRowDiag[0], 0, &nByRowOffDiag[0]);
} else {
MatSeqBAIJSetPreallocation(_a, _blockSize, 0, &nByRowDiag[0]);
MatMPIBAIJSetPreallocation(_a, _blockSize, 0, &nByRowDiag[0], 0, &nByRowOffDiag[0]);
}
_sparsity.clear();
}
_entriesPreAllocated = true;
}
void linearSystemPETScBlockDouble::zeroMatrix()
{
if (_isAllocated && _entriesPreAllocated) {
MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY);
MatZeroEntries(_a);
}
}
void linearSystemPETScBlockDouble::zeroRightHandSide()
{
if (_isAllocated) {
VecAssemblyBegin(_b);
VecAssemblyEnd(_b);
VecZeroEntries(_b);
}
}
void linearSystemPETScBlockDouble::zeroSolution()
{
if (_isAllocated) {
VecAssemblyBegin(_x);
VecAssemblyEnd(_x);
VecZeroEntries(_x);
}
}
linearSystemPETScBlockDouble::linearSystemPETScBlockDouble(bool sequential)
{
_entriesPreAllocated = false;
_isAllocated = false;
_kspAllocated = false;
_matrixModified=true;
_sequential = sequential;
}
double linearSystemPETScBlockDouble::normInfRightHandSide() const
{
PetscReal nor;
VecAssemblyBegin(_b);
VecAssemblyEnd(_b);
VecNorm(_b, NORM_INFINITY, &nor);
return nor;
}
void linearSystemPETScBlockDouble::printMatlab(const char *filename) const
{
_try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
_try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
_try(VecAssemblyBegin(_b));
_try(VecAssemblyEnd(_b));
PetscViewer viewer;
PetscViewerASCIIOpen(PETSC_COMM_WORLD, filename, &viewer);
PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
MatView(_a, viewer);
PetscViewerDestroy(&viewer);
return;
}
template class linearSystemPETSc<fullMatrix<double> >;
#endif // HAVE_PETSC
Solver/linearSystemPETSc.h
+ 33
− 64
View file @ 999753cd
......@@ -65,7 +65,6 @@ typedef struct _p_KSP* KSP;
template <class scalar>
class linearSystemPETSc : public linearSystem<scalar> {
protected:
int _blockSize; // for block Matrix
bool _isAllocated, _kspAllocated, _entriesPreAllocated;
bool _matrixChangedSinceLastSolve;
bool _valuesNotAssembled; //cannot use MatAssembled since MatAssembled return false for an empty matrix
......@@ -79,25 +78,27 @@ class linearSystemPETSc : public linearSystem<scalar> {
#ifndef SWIG
MPI_Comm _comm;
#endif
int _getBlockSizeFromParameters() const;
public:
virtual ~linearSystemPETSc();
~linearSystemPETSc();
void insertInSparsityPattern (int i, int j);
virtual bool isAllocated() const { return _isAllocated; }
virtual void preAllocateEntries();
virtual void allocate(int nbRows);
bool isAllocated() const { return _isAllocated; }
void preAllocateEntries();
void allocate(int nbRows);
void print();
virtual void clear();
virtual void getFromMatrix(int row, int col, scalar &val) const;
virtual void addToRightHandSide(int row, const scalar &val);
virtual void getFromRightHandSide(int row, scalar &val) const;
virtual double normInfRightHandSide() const;
virtual void addToMatrix(int row, int col, const scalar &val);
virtual void addToSolution(int row, const scalar &val);
virtual void getFromSolution(int row, scalar &val) const;
virtual void zeroMatrix();
virtual void zeroRightHandSide();
virtual void zeroSolution();
virtual int systemSolve();
void clear();
void getFromMatrix(int row, int col, scalar &val) const;
void addToRightHandSide(int row, const scalar &val);
void getFromRightHandSide(int row, scalar &val) const;
double normInfRightHandSide() const;
void addToMatrix(int row, int col, const scalar &val);
void addToSolution(int row, const scalar &val);
void getFromSolution(int row, scalar &val) const;
void zeroMatrix();
void zeroRightHandSide();
void zeroSolution();
void printMatlab(const char *filename) const;
int systemSolve();
Mat &getMatrix(){ return _a; }
//std::vector<scalar> getData();
//std::vector<int> getRowPointers();
......@@ -108,39 +109,6 @@ class linearSystemPETSc : public linearSystem<scalar> {
#endif
linearSystemPETSc();
};
class linearSystemPETScBlockDouble : public linearSystem<fullMatrix<double> > {
bool _entriesPreAllocated, _isAllocated, _kspAllocated, _matrixModified;
sparsityPattern _sparsity;
Mat _a;
Vec _b, _x;
KSP _ksp;
int _blockSize, _localRowStart, _localRowEnd, _globalSize, _localSize;
bool _sequential;
public:
void _kspCreate();
void print();
void printMatlab(const char *filename) const;
virtual void addToMatrix(int row, int col, const fullMatrix<double> &val);
virtual void addToRightHandSide(int row, const fullMatrix<double> &val);
virtual void addToSolution(int row, const fullMatrix<double> &val);
virtual void getFromMatrix(int row, int col, fullMatrix<double> &val ) const;
virtual void getFromRightHandSide(int row, fullMatrix<double> &val) const;
virtual void getFromSolution(int row, fullMatrix<double> &val) const;
void allocate(int nbRows);
bool isAllocated() const;
int systemSolve();
void preAllocateEntries();
void clear();
void zeroMatrix();
void zeroRightHandSide();
void zeroSolution();
double normInfRightHandSide() const;
void insertInSparsityPattern (int i, int j);
linearSystemPETScBlockDouble(bool sequential = false);
~linearSystemPETScBlockDouble();
};
#else
template <class scalar>
......@@ -150,20 +118,21 @@ class linearSystemPETSc : public linearSystem<scalar> {
{
Msg::Error("PETSc is not available in this version of Gmsh");
}
virtual bool isAllocated() const { return false; }
virtual void allocate(int nbRows) {}
virtual void clear(){}
virtual void addToMatrix(int row, int col, const scalar &val) {}
virtual void getFromMatrix(int row, int col, scalar &val) const {}
virtual void addToRightHandSide(int row, const scalar &val) {}
virtual void addToSolution(int row, const scalar &val) {}
virtual void getFromRightHandSide(int row, scalar &val) const {}
virtual void getFromSolution(int row, scalar &val) const {}
virtual void zeroMatrix() {}
virtual void zeroRightHandSide() {}
virtual void zeroSolution() {}
virtual int systemSolve() { return 0; }
virtual double normInfRightHandSide() const{return 0;}
bool isAllocated() const { return false; }
void allocate(int nbRows) {}
void clear(){}
void addToMatrix(int row, int col, const scalar &val) {}
void getFromMatrix(int row, int col, scalar &val) const {}
void addToRightHandSide(int row, const scalar &val) {}
void addToSolution(int row, const scalar &val) {}
void getFromRightHandSide(int row, scalar &val) const {}
void getFromSolution(int row, scalar &val) const {}
void zeroMatrix() {}
void zeroRightHandSide() {}
void zeroSolution() {}
void printMatlab(const char *filename) const{};
int systemSolve() { return 0; }
double normInfRightHandSide() const{return 0;}
};
#endif
#endif
Solver/linearSystemPETSc.hpp
+ 47
− 19
View file @ 999753cd
......@@ -9,6 +9,12 @@ static void _try(int ierr)
CHKERRABORT(PETSC_COMM_WORLD, ierr);
}
template <class scalar>
int linearSystemPETSc<scalar>::_getBlockSizeFromParameters() const
{
return 1;
}
template <class scalar>
void linearSystemPETSc<scalar>::_kspCreate()
{
......@@ -34,10 +40,10 @@ linearSystemPETSc<scalar>::linearSystemPETSc(MPI_Comm com)
{
_comm = com;
_isAllocated = false;
_blockSize = 0;
_kspAllocated = false;
_matrixChangedSinceLastSolve = true;
_valuesNotAssembled = false;
_entriesPreAllocated = false;
}
template <class scalar>
......@@ -45,10 +51,10 @@ linearSystemPETSc<scalar>::linearSystemPETSc()
{
_comm = PETSC_COMM_WORLD;
_isAllocated = false;
_blockSize = 0;
_kspAllocated = false;
_matrixChangedSinceLastSolve = true;
_valuesNotAssembled = false;
_entriesPreAllocated = false;
}
template <class scalar>
......@@ -84,17 +90,15 @@ void linearSystemPETSc<scalar>::preAllocateEntries()
{
if (_entriesPreAllocated) return;
if (!_isAllocated) Msg::Fatal("system must be allocated first");
int blockSize = _getBlockSizeFromParameters();
std::vector<int> nByRowDiag (_localSize), nByRowOffDiag (_localSize);
if (_sparsity.getNbRows() == 0) {
PetscInt prealloc = 300;
PetscBool set;
PetscOptionsGetInt(PETSC_NULL, "-petsc_prealloc", &prealloc, &set);
if (_blockSize == 0) {
_try(MatSeqAIJSetPreallocation(_a, prealloc, PETSC_NULL));
} else {
_try(MatSeqBAIJSetPreallocation(_a, _blockSize, 5, PETSC_NULL));
}
nByRowDiag.resize(0);
nByRowDiag.resize(_localSize, prealloc);
} else {
std::vector<int> nByRowDiag (_localSize), nByRowOffDiag (_localSize);
for (int i = 0; i < _localSize; i++) {
int n;
const int *r = _sparsity.getRow(i, n);
......@@ -105,24 +109,34 @@ void linearSystemPETSc<scalar>::preAllocateEntries()
nByRowOffDiag[i] ++;
}
}
if (_blockSize == 0) {
_try(MatSeqAIJSetPreallocation(_a, 0, &nByRowDiag[0]));
_try(MatMPIAIJSetPreallocation(_a, 0, &nByRowDiag[0], 0, &nByRowOffDiag[0]));
} else {
_try(MatSeqBAIJSetPreallocation(_a, _blockSize, 0, &nByRowDiag[0]));
_try(MatMPIBAIJSetPreallocation(_a, _blockSize, 0, &nByRowDiag[0], 0, &nByRowOffDiag[0]));
}
_sparsity.clear();
}
_try(MatXAIJSetPreallocation(_a, blockSize, &nByRowDiag[0], &nByRowOffDiag[0], NULL, NULL));
_entriesPreAllocated = true;
}
template <class scalar>
void linearSystemPETSc<scalar>::allocate(int nbRows)
{
#ifdef HAVE_MPI
PetscMPIInt commSize;
MPI_Comm_size(_comm,&commSize);
#endif
int blockSize = _getBlockSizeFromParameters();
clear();
_try(MatCreate(_comm, &_a));
_try(MatSetSizes(_a, nbRows, nbRows, PETSC_DETERMINE, PETSC_DETERMINE));
_try(MatSetSizes(_a, blockSize * nbRows, blockSize * nbRows, PETSC_DETERMINE, PETSC_DETERMINE));
if (blockSize > 1) {
#ifdef HAVE_MPI
if (commSize > 1) {
MatSetType(_a, MATMPIBAIJ);
}
else
#endif
{
MatSetType(_a, MATSEQBAIJ);
}
}
// override the default options with the ones from the option
// database (if any)
if (this->_parameters.count("petscOptions"))
......@@ -133,8 +147,6 @@ void linearSystemPETSc<scalar>::allocate(int nbRows)
//since PETSc 3.3 GetOwnershipRange and MatGetSize cannot be called before MatXXXSetPreallocation
_localSize = nbRows;
#ifdef HAVE_MPI
PetscMPIInt commSize;
MPI_Comm_size(_comm,&commSize);
if (commSize>1){
_localRowStart = 0;
if (Msg::GetCommRank() != 0) {
......@@ -159,7 +171,7 @@ void linearSystemPETSc<scalar>::allocate(int nbRows)
#endif
// preallocation option must be set after other options
_try(VecCreate(_comm, &_x));
_try(VecSetSizes(_x, nbRows, PETSC_DETERMINE));
_try(VecSetSizes(_x, blockSize * nbRows, PETSC_DETERMINE));
// override the default options with the ones from the option
// database (if any)
if (this->_parameters.count("petscPrefix"))
......@@ -345,6 +357,22 @@ int linearSystemPETSc<scalar>::systemSolve()
return 1;
}
template <class scalar>
void linearSystemPETSc<scalar>::printMatlab(const char *filename) const
{
_try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
_try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
_try(VecAssemblyBegin(_b));
_try(VecAssemblyEnd(_b));
PetscViewer viewer;
PetscViewerASCIIOpen(PETSC_COMM_WORLD, filename, &viewer);
PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
MatView(_a, viewer);
PetscViewerDestroy(&viewer);
return;
}
/*template <class scalar>
std::vector<scalar> linearSystemPETSc<scalar>::getData()
{
......
wrappers/gmshpy/gmshSolver.i
+ 1
− 0
View file @ 999753cd
......@@ -35,6 +35,7 @@
#if defined(HAVE_PETSC)
%include "linearSystemPETSc.h"
%template(linearSystemPETScDouble) linearSystemPETSc<double>;
%template(linearSystemPETScBlockDouble) linearSystemPETSc<fullMatrix<double> >;
#endif
%include "eigenSolver.h"
#endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment