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

get rid of random errors in linearSystemPETScBlock by changing template 

structure (not sure what was wrong though)
parent 11ae31b5
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Common/gmshpy.i
+ 10
− 12
View file @ 78ba2114
......@@ -11,7 +11,6 @@
#include "linearSystem.h"
#include "linearSystemFull.h"
#include "linearSystemPETSc.h"
#include "linearSystemPETSc.cpp" // needed for the specialization
#include "linearSystemCSR.h"
#include "GEntity.h"
#include "GVertex.h"
......@@ -79,15 +78,25 @@ namespace std {
%template(GFaceList) list<GFace*, std::allocator<GFace*> >;
}
%include "GmshConfig.h"
%include "Context.h"
%include "fullMatrix.h"
%template(fullMatrixDouble) fullMatrix<double>;
%template(fullVectorDouble) fullVector<double>;
%include "dofManager.h"
%template(dofManagerDouble) dofManager<double>;
%include "GModel.h"
%include "function.h"
%include "linearSystem.h"
%template(linearSystemDouble) linearSystem<double>;
%template(linearSystemFullMatrixDouble) linearSystem<fullMatrix<double> >;
%include "linearSystemFull.h"
%template(linearSystemFullDouble) linearSystemFull<double> ;
%include "linearSystemPETSc.h"
%template(linearSystemPETScDouble) linearSystemPETSc<double>;
%include "linearSystemCSR.h"
%template(linearSystemTAUCSDouble) linearSystemCSRTaucs<double>;
%include "GEntity.h"
%include "GVertex.h"
%include "GEdge.h"
......@@ -104,7 +113,6 @@ namespace std {
%include "polynomialBasis.h"
%include "Gauss.h"
%include "meshPartitionOptions.h"
%include "linearSystemCSR.h"
%include "elasticitySolver.h"
%include "PView.h"
%include "PViewData.h"
......@@ -115,15 +123,5 @@ namespace std {
%include "SPoint3.h"
%include "SPoint2.h"
%include "GPoint.h"
%template(dofManagerDouble) dofManager<double>;
%template(linearSystemDouble) linearSystem<double>;
%template(linearSystemFullMatrixDouble) linearSystem<fullMatrix<double> >;
%template(linearSystemFullDouble) linearSystemFull<double> ;
%template(linearSystemPETScDouble) linearSystemPETSc<double>;
%template(linearSystemTAUCSDouble) linearSystemCSRTaucs<double>;
%template(fullMatrixDouble) fullMatrix<double>;
%template(fullVectorDouble) fullVector<double>;
%template(linearSystemPETScBlockDouble) linearSystemPETSc<fullMatrix<double> >;
%include "functionPython.h"
Solver/elasticitySolver.h
+ 0
− 1
View file @ 78ba2114
......@@ -86,7 +86,6 @@ class elasticitySolver
void addDirichletBCLua (int dim, int entityId, int component, std::string luaFunctionName, lua_State *L);
void addNeumannBCLua (int dim, int entityId, std::string luaFunctionName, lua_State *L);
void addNeumannBCFct (int dim, int entityId, const function* luaFunction, lua_State *L);
#endif
......
Solver/linearSystemPETSc.cpp
+ 138
− 55
View file @ 78ba2114
......@@ -6,8 +6,15 @@
#include <stdlib.h>
#include "GmshMessage.h"
template <>
void linearSystemPETSc<fullMatrix<PetscScalar> >::addToMatrix(int row, int col, const fullMatrix<PetscScalar> &val)
void linearSystemPETScBlockDouble::_kspCreate() {
KSPCreate(PETSC_COMM_WORLD, &_ksp);
if (this->_parameters.count("petscPrefix"))
KSPAppendOptionsPrefix(_ksp, this->_parameters["petscPrefix"].c_str());
KSPSetFromOptions(_ksp);
_kspAllocated = true;
}
void linearSystemPETScBlockDouble::addToMatrix(int row, int col, const fullMatrix<PetscScalar> &val)
{
if (!_entriesPreAllocated)
preAllocateEntries();
......@@ -19,7 +26,7 @@ void linearSystemPETSc<fullMatrix<PetscScalar> >::addToMatrix(int row, int col,
modval(jj, ii) = buff;
}
PetscInt i = row, j = col;
_try(MatSetValuesBlocked(_a, 1, &i, 1, &j, &modval(0,0), ADD_VALUES));
MatSetValuesBlocked(_a, 1, &i, 1, &j, &modval(0,0), ADD_VALUES);
//transpose back so that the original matrix is not modified
for (int ii = 0; ii < val.size1(); ii++)
for (int jj = 0; jj < ii; jj++) {
......@@ -29,99 +36,175 @@ void linearSystemPETSc<fullMatrix<PetscScalar> >::addToMatrix(int row, int col,
}
}
template<>
void linearSystemPETSc<fullMatrix<PetscScalar> >::addToRightHandSide(int row, const fullMatrix<PetscScalar> &val)
void linearSystemPETScBlockDouble::addToRightHandSide(int row, const fullMatrix<PetscScalar> &val)
{
for (int ii = 0; ii < _blockSize; ii++) {
PetscInt i = row * _blockSize + ii;
PetscScalar v = val(ii, 0);
_try(VecSetValues(_b, 1, &i, &v, ADD_VALUES));
VecSetValues(_b, 1, &i, &v, ADD_VALUES);
}
}
template<>
void linearSystemPETSc<fullMatrix<PetscScalar> >::getFromMatrix(int row, int col, fullMatrix<PetscScalar> &val ) const
void linearSystemPETScBlockDouble::getFromMatrix(int row, int col, fullMatrix<PetscScalar> &val ) const
{
Msg::Error("getFromMatrix not implemented for PETSc");
}
template<>
void linearSystemPETSc<fullMatrix<PetscScalar> >::getFromRightHandSide(int row, fullMatrix<PetscScalar> &val) const
void linearSystemPETScBlockDouble::getFromRightHandSide(int row, fullMatrix<PetscScalar> &val) const
{
#if defined(PETSC_USE_COMPLEX)
PetscScalar *tmp;
_try(VecGetArray(_b, &tmp));
for (int i = 0; i < _blockSize; i++) {
PetscScalar s = tmp[row * _blockSize + i];
val(i,0) = s.real();
}
_try(VecRestoreArray(_b, &tmp));
#else
for (int i = 0; i < _blockSize; i++) {
int ii = row*_blockSize +i;
VecGetValues ( _b, 1, &ii, &val(i,0));
}
#endif
}
template<>
void linearSystemPETSc<fullMatrix<PetscScalar> >::getFromSolution(int row, fullMatrix<PetscScalar> &val) const {
#if defined(PETSC_USE_COMPLEX)
PetscScalar *tmp;
_try(VecGetArray(_x, &tmp));
for (int i = 0; i < _blockSize; i++) {
PetscScalar s = tmp[row * _blockSize + i];
val(i,0) = s.real();
}
_try(VecRestoreArray(_x, &tmp));
#else
void linearSystemPETScBlockDouble::getFromSolution(int row, fullMatrix<PetscScalar> &val) const
{
for (int i = 0; i < _blockSize; i++) {
int ii = row*_blockSize +i;
VecGetValues ( _x, 1, &ii, &val(i,0));
}
#endif
}
template<>
void linearSystemPETSc<fullMatrix<PetscScalar> >::allocate(int nbRows)
void linearSystemPETScBlockDouble::allocate(int nbRows)
{
if (this->_parameters.count("petscOptions"))
_try(PetscOptionsInsertString(this->_parameters["petscOptions"].c_str()));
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();
//printf("allocate linear system petsc \n");
//_try(PetscOptionsInsertString("-ksp_monitor_true_residual -ksp_rtol 1e-10"));
_try(MatCreate(PETSC_COMM_WORLD, &_a));
_try(MatSetSizes(_a,nbRows * _blockSize, nbRows * _blockSize, PETSC_DETERMINE, PETSC_DETERMINE));
MatCreate(PETSC_COMM_WORLD, &_a);
MatSetSizes(_a,nbRows * _blockSize, nbRows * _blockSize, PETSC_DETERMINE, PETSC_DETERMINE);
if (Msg::GetCommSize() > 1) {
_try(MatSetType(_a, MATMPIBAIJ));
MatSetType(_a, MATMPIBAIJ);
} else {
_try(MatSetType(_a, MATSEQBAIJ));
MatSetType(_a, MATSEQBAIJ);
}
if (_parameters.count("petscPrefix"))
_try(MatAppendOptionsPrefix(_a, _parameters["petscPrefix"].c_str()));
_try(MatSetFromOptions(_a));
_try(MatGetOwnershipRange(_a, &_localRowStart, &_localRowEnd));
_try(MatGetSize(_a, &_globalSize, &_localSize));
MatAppendOptionsPrefix(_a, _parameters["petscPrefix"].c_str());
MatSetFromOptions(_a);
MatGetOwnershipRange(_a, &_localRowStart, &_localRowEnd);
MatGetSize(_a, &_globalSize, &_localSize);
_globalSize /= _blockSize;
_localSize /= _blockSize;
_localRowStart /= _blockSize;
_localRowEnd /= _blockSize;
// override the default options with the ones from the option
// database (if any)
_try(VecCreate(PETSC_COMM_WORLD, &_x));
_try(VecSetSizes(_x, nbRows * _blockSize, PETSC_DETERMINE));
VecCreate(PETSC_COMM_WORLD, &_x);
VecSetSizes(_x, nbRows * _blockSize, PETSC_DETERMINE);
// override the default options with the ones from the option
// database (if any)
if (_parameters.count("petscPrefix"))
_try(VecAppendOptionsPrefix(_x, _parameters["petscPrefix"].c_str()));
_try(VecSetFromOptions(_x));
_try(VecDuplicate(_x, &_b));
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;
}
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")
KSPSetOperators(_ksp, _a, _a, SAME_PRECONDITIONER);
else
KSPSetOperators(_ksp, _a, _a, DIFFERENT_NONZERO_PATTERN);
MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY);
VecAssemblyBegin(_b);
VecAssemblyEnd(_b);
KSPSolve(_ksp, _b, _x);
return 1;
}
void linearSystemPETScBlockDouble::preAllocateEntries()
{
if (_entriesPreAllocated) return;
if (!_isAllocated) Msg::Fatal("system must be allocated first");
if (_sparsity.getNbRows() == 0) {
PetscInt prealloc = 300;
PetscTruth 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);
}
}
linearSystemPETScBlockDouble::linearSystemPETScBlockDouble()
{
_entriesPreAllocated = false;
_isAllocated = false;
_kspAllocated = false;
printf("init\n");
}
double linearSystemPETScBlockDouble::normInfRightHandSide() const
{
PetscReal nor;
VecNorm(_b, NORM_INFINITY, &nor);
return nor;
}
#include "Bindings.h"
void linearSystemPETScRegisterBindings(binding *b)
......@@ -136,12 +219,12 @@ void linearSystemPETScRegisterBindings(binding *b)
cm->setDescription ("A new PETSc<PetscScalar> solver");
cb->setParentClass<linearSystem<PetscScalar> >();
cm->setArgNames(NULL);
cb = b->addClass<linearSystemPETSc<fullMatrix<PetscScalar> > >("linearSystemPETScBlock");
/* cb = b->addClass<linearSystemPETScBlockDouble >("linearSystemPETScBlock");
cb->setDescription("A linear system solver, based on PETSc");
cm = cb->setConstructor<linearSystemPETSc<fullMatrix<PetscScalar> > >();
cm->setDescription ("A new PETScBlock<PetscScalar> solver");
cb->setParentClass<linearSystem<fullMatrix<PetscScalar> > >();
#endif // FIXME
cm = cb->setConstructor<linearSystemPETScBlockDouble >();
cm->setDescription ("A new PETScBlockDouble solver");
cb->setParentClass<linearSystem<fullMatrix<PetscScalar> > >();*/
#endif
}
#endif // HAVE_PETSC
Solver/linearSystemPETSc.h
+ 29
− 3
View file @ 78ba2114
......@@ -236,7 +236,7 @@ class linearSystemPETSc : public linearSystem<scalar> {
_try(VecZeroEntries(_b));
}
}
int systemSolve()
virtual int systemSolve()
{
if (!_kspAllocated)
_kspCreate();
......@@ -266,10 +266,34 @@ class linearSystemPETSc : public linearSystem<scalar> {
class binding;
void linearSystemPETScRegisterBindings(binding *b);
class linearSystemPETScBlockDouble : public linearSystem<fullMatrix<double> > {
bool _entriesPreAllocated, _isAllocated, _kspAllocated;
sparsityPattern _sparsity;
Mat _a;
Vec _b, _x;
KSP _ksp;
int _blockSize, _localRowStart, _localRowEnd, _globalSize, _localSize;
public:
void _kspCreate();
virtual void addToMatrix(int row, int col, const fullMatrix<double> &val);
virtual void addToRightHandSide(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();
double normInfRightHandSide() const;
linearSystemPETScBlockDouble();
};
#else
template <class scalar>
/*template <class scalar>
class linearSystemPETSc : public linearSystem<scalar> {
public :
linearSystemPETSc()
......@@ -288,8 +312,10 @@ class linearSystemPETSc : public linearSystem<scalar> {
virtual void zeroRightHandSide() {}
virtual int systemSolve() { return 0; }
virtual double normInfRightHandSide() const{return 0;}
};
};*/
#endif
#endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment