Skip to content
Snippets Groups Projects
Select Git revision
  • c6063bc44a0555792b01e66d106ccf76d1c80aed
  • master default protected
  • hierarchical-basis
  • alphashapes
  • bl
  • relaying
  • new_export_boris
  • oras_vs_osm
  • reassign_partitions
  • distributed_fwi
  • rename-classes
  • fix/fortran-api-example-t4
  • robust_partitions
  • reducing_files
  • fix_overlaps
  • 3115-issue-fix
  • 3023-Fillet2D-Update
  • convert_fdivs
  • tmp_jcjc24
  • fixedMeshIF
  • save_edges
  • gmsh_4_14_0
  • gmsh_4_13_1
  • gmsh_4_13_0
  • gmsh_4_12_2
  • gmsh_4_12_1
  • gmsh_4_12_0
  • gmsh_4_11_1
  • 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
41 results

linearSystemGMM.h

Blame
  • linearSystemGMM.h 3.96 KiB
    // Gmsh - Copyright (C) 1997-2019 C. Geuzaine, J.-F. Remacle
    //
    // See the LICENSE.txt file for license information. Please report all
    // issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
    
    #ifndef _LINEAR_SYSTEM_GMM_H_
    #define _LINEAR_SYSTEM_GMM_H_
    
    // Interface to GMM++
    
    #include "GmshConfig.h"
    #include "GmshMessage.h"
    #include "linearSystem.h"
    
    #if defined(HAVE_GMM)
    #include <gmm.h>
    
    template <class scalar> class linearSystemGmm : public linearSystem<scalar> {
    protected:
      std::vector<scalar>
        *_x; // the nonLinearSystemGmm has to access to this vector
      std::vector<scalar> *_b; // idem
      gmm::row_matrix<gmm::wsvector<scalar> > *_a;
    
    private:
      double _prec;
      int _noisy, _gmres;
    
    public:
      linearSystemGmm() : _x(0), _b(0), _a(0), _prec(1.e-8), _noisy(0), _gmres(0) {}
      virtual bool isAllocated() const { return _a != 0; }
      virtual void allocate(int nbRows)
      {
        clear();
        _a = new gmm::row_matrix<gmm::wsvector<scalar> >(nbRows, nbRows);
        _b = new std::vector<scalar>(nbRows);
        _x = new std::vector<scalar>(nbRows);
      }
      virtual ~linearSystemGmm() { clear(); }
      virtual void clear()
      {
        if(_a) {
          delete _a;
          delete _b;
          delete _x;
        }
        _a = 0;
      }
    
      virtual void addToMatrix(int row, int col, const scalar &val)
      {
        if(val != 0.0) (*_a)(row, col) += val;
      }
      virtual void getFromMatrix(int row, int col, scalar &val) const
      {
        val = (*_a)(row, col);
      }
    
      virtual void addToRightHandSide(int row, const scalar &val, int ith = 0)
      {
        if(val != 0.0) (*_b)[row] += val;
      }
    
      virtual void getFromRightHandSide(int row, scalar &val) const
      {
        val = (*_b)[row];
      }
    
      virtual void addToSolution(int row, const scalar &val)
      {
        if(val != 0.0) (*_x)[row] += val;
      }
    
      virtual void getFromSolution(int row, scalar &val) const { val = (*_x)[row]; }
      virtual void zeroMatrix() { gmm::clear(*_a); }
      virtual void zeroRightHandSide()
      {
        for(unsigned int i = 0; i < _b->size(); i++) (*_b)[i] = 0.;
      }
      virtual void zeroSolution()
      {
        for(unsigned int i = 0; i < _x->size(); i++) (*_x)[i] = 0.;
      }
    
      virtual double normInfRightHandSide() const
      {
        double nor = 0.;
        double temp;
        for(unsigned int i = 0; i < _b->size(); i++) {
          temp = abs((*_b)[i]); // this is valid also for complex
          // if(temp<0) temp = -temp;
          if(nor < temp) nor = temp;
        }
        return nor;
      }
    
      void setPrec(double p) { _prec = p; }
      void setNoisy(int n) { _noisy = n; }
      void setGmres(int n) { _gmres = n; }
      virtual int systemSolve()
      {
    #if defined(HAVE_MUMPS)
        gmm::MUMPS_solve(*_a, *_x, *_b);
    #else
        // gmm::ilutp_precond<gmm::row_matrix<gmm::wsvector<scalar> > > P(*_a, 25,
        // 0.);
        gmm::ildltt_precond<gmm::row_matrix<gmm::wsvector<scalar> > > P(*_a, 30,
                                                                        1.e-10);
        gmm::iteration iter(_prec);
        iter.set_noisy(_noisy);
        if(_gmres)
          gmm::gmres(*_a, *_x, *_b, P, 100, iter);
        else
          gmm::cg(*_a, *_x, *_b, P, iter);
    #endif
        return 1;
      }
    };
    
    #else
    
    template <class scalar> class linearSystemGmm : public linearSystem<scalar> {
    public:
      linearSystemGmm()
      {
        Msg::Error("Gmm++ is not available in this version of Gmsh");
      }
      virtual bool isAllocated() const { return false; }
      virtual void allocate(int nbRows) {}
      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, int ith = 0) {}
      virtual void getFromRightHandSide(int row, scalar &val) const {}
      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() { return 0; }
      virtual double normInfRightHandSide() const { return 0.; }
      void setPrec(double p) {}
      virtual void clear() {}
      void setNoisy(int n) {}
      void setGmres(int n) {}
    };
    
    #endif
    
    #endif