Select Git revision
Options.cpp
Forked from
gmsh / gmsh
Source project has a limited visibility.
-
Christophe Geuzaine authoredChristophe Geuzaine authored
linearSystemPETSc.cpp 4.08 KiB
#include "GmshConfig.h"
#if defined(HAVE_PETSC)
#include "linearSystemPETSc.h"
#include "fullMatrix.h"
#include <stdlib.h>
#include "GmshMessage.h"
template <>
void linearSystemPETSc<fullMatrix<PetscScalar> >::addToMatrix(int row, int col, const fullMatrix<PetscScalar> &val)
{
fullMatrix<PetscScalar> &modval = *const_cast<fullMatrix<PetscScalar> *>(&val);
for (int ii = 0; ii < val.size1(); ii++)
for (int jj = 0; jj < ii; jj++) {
PetscScalar buff = modval(ii, jj);
modval(ii, jj) = modval (jj, ii);
modval(jj, ii) = buff;
}
PetscInt i = row, j = col;
_try(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++) {
PetscScalar buff = modval(ii,jj);
modval(ii, jj) = modval (jj,ii);
modval(jj, ii) = buff;
}
}
template<>
void linearSystemPETSc<fullMatrix<PetscScalar> >::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));
}
}
template<>
void linearSystemPETSc<fullMatrix<PetscScalar> >::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
{
PetscScalar *tmp;
_try(VecGetArray(_b, &tmp));
for (int i = 0; i < _blockSize; i++) {
PetscScalar s = tmp[row * _blockSize + i];
#if defined(PETSC_USE_COMPLEX)
val(i,0) = s.real();
#else
val(i,0) = s;
#endif
}
_try(VecRestoreArray(_b, &tmp));
}
template<>
void linearSystemPETSc<fullMatrix<PetscScalar> >::getFromSolution(int row, fullMatrix<PetscScalar> &val) const
{
PetscScalar *tmp;
_try(VecGetArray(_x, &tmp));
for (int i = 0; i < _blockSize; i++) {
PetscScalar s = tmp[row * _blockSize + i];
#if defined(PETSC_USE_COMPLEX)
val(i,0) = s.real();
#else
val(i,0) = s;
#endif
}
_try(VecRestoreArray(_x, &tmp));
}
template<>
void linearSystemPETSc<fullMatrix<PetscScalar> >::allocate(int nbRows)
{
_blockSize = strtol (_parameters["blockSize"].c_str(), NULL, 10);
if (_blockSize == 0)
Msg::Error ("'blockSize' parameters must be set for linearSystemPETScBlock");
clear();
_try(MatCreate(PETSC_COMM_WORLD, &_a));
_try(MatSetSizes(_a, PETSC_DECIDE, PETSC_DECIDE, nbRows * _blockSize, nbRows * _blockSize));
_try(MatSetType(_a, MATSEQBAIJ));
// override the default options with the ones from the option
// database (if any)
_try(MatSetFromOptions(_a));
_try(MatSeqBAIJSetPreallocation(_a, _blockSize, 5, NULL)); //todo preAllocate off-diagonal part
//_try(MatMPIBAIJSetPreallocation(_a, _blockSize, 4, NULL, 0, NULL)); //todo preAllocate off-diagonal part
_try(VecCreate(PETSC_COMM_WORLD, &_x));
_try(VecSetSizes(_x, PETSC_DECIDE, nbRows * _blockSize));
// override the default options with the ones from the option
// database (if any)
_try(VecSetFromOptions(_x));
_try(VecDuplicate(_x, &_b));
_isAllocated = true;
}
#include "Bindings.h"
void linearSystemPETScRegisterBindings(binding *b)
{
// FIXME on complex arithmetic this crashes
#if !defined(PETSC_USE_COMPLEX)
classBinding *cb;
methodBinding *cm;
cb = b->addClass<linearSystemPETSc<PetscScalar> >("linearSystemPETSc");
cb->setDescription("A linear system solver, based on PETSc");
cm = cb->setConstructor<linearSystemPETSc<PetscScalar> >();
cm->setDescription ("A new PETSc<PetscScalar> solver");
cb->setParentClass<linearSystem<PetscScalar> >();
cm->setArgNames(NULL);
cb = b->addClass<linearSystemPETSc<fullMatrix<PetscScalar> > >("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
}
#endif // HAVE_PETSC