Skip to content
Snippets Groups Projects
Select Git revision
  • 777465d7c2c40c4b7fe87af9a36b12cf3a9ef059
  • master default protected
  • alphashapes
  • quadMeshingTools
  • cygwin_conv_path
  • macos_arm64
  • add-transfiniteautomatic-to-geo
  • patch_releases_4_10
  • HierarchicalHDiv
  • isuruf-master-patch-63355
  • hyperbolic
  • hexdom
  • hxt_update
  • jf
  • 1618-pythonocc-and-gmsh-api-integration
  • octreeSizeField
  • hexbl
  • alignIrregularVertices
  • getEdges
  • patch_releases_4_8
  • isuruf-master-patch-51992
  • gmsh_4_11_0
  • gmsh_4_10_5
  • gmsh_4_10_4
  • gmsh_4_10_3
  • gmsh_4_10_2
  • gmsh_4_10_1
  • gmsh_4_10_0
  • gmsh_4_9_5
  • gmsh_4_9_4
  • gmsh_4_9_3
  • gmsh_4_9_2
  • gmsh_4_9_1
  • gmsh_4_9_0
  • gmsh_4_8_4
  • gmsh_4_8_3
  • gmsh_4_8_2
  • gmsh_4_8_1
  • gmsh_4_8_0
  • gmsh_4_7_1
  • gmsh_4_7_0
41 results

OptHomRun.cpp

Blame
  • Forked from gmsh / gmsh
    Source project has a limited visibility.
    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