Skip to content
Snippets Groups Projects
Select Git revision
1 result Searching

GF_LaplacexForm.cpp

Blame
  • 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