Skip to content
Snippets Groups Projects
Select Git revision
  • fb3251cd5f312f9512fff4ada6ad60f02b1a0c0a
  • master default
  • cgnsUnstructured
  • partitioning
  • poppler
  • HighOrderBLCurving
  • gmsh_3_0_4
  • gmsh_3_0_3
  • gmsh_3_0_2
  • gmsh_3_0_1
  • gmsh_3_0_0
  • gmsh_2_16_0
  • gmsh_2_15_0
  • gmsh_2_14_1
  • gmsh_2_14_0
  • gmsh_2_13_2
  • gmsh_2_13_1
  • gmsh_2_12_0
  • gmsh_2_11_0
  • gmsh_2_10_1
  • gmsh_2_10_0
  • gmsh_2_9_3
  • gmsh_2_9_2
  • gmsh_2_9_1
  • gmsh_2_9_0
  • gmsh_2_8_6
26 results

GmshMatrix.h

Blame
  • Forked from gmsh / gmsh
    Source project has a limited visibility.
    GmshMatrix.h 4.71 KiB
    #ifndef _GMSH_BOOSTMATRIX_
    #define _GMSH_BOOSTMATRIX_
    
    #include <assert.h>
    
    template <class SCALAR>
    class Gmsh_Vector
    {
    private:
      int r;
    public:
      inline int size() const {return r;}
      SCALAR *data;
      ~Gmsh_Vector() {delete [] data;}
      Gmsh_Vector(int R)
        : r(R)
      {
        data = new SCALAR [r];
        scal(0);
      }
      Gmsh_Vector(const Gmsh_Vector<SCALAR> &other)
        : r(other.r)
      {
        data = new double [r];
        copy ( other.data ) ;
      }
      inline SCALAR operator () (int i) const
      {
        return data[i];
      }
      inline SCALAR & operator () (int i)
      {
        return data[i];
      }
      inline SCALAR operator *(const Gmsh_Vector<SCALAR> & other)
      {
        throw;
      }
      inline void scal ( const SCALAR s)
      {
        for (int i=0;i<r;++i)data[i]*=s;
      }
      inline void copy ( const SCALAR **other)
      {
        for (int i=0;i<r;++i)data[i]=other.data[i];
      }
      inline void lu_solve (const  Gmsh_Vector<SCALAR>& rhs,  Gmsh_Vector<SCALAR> & result)
          {
    	  throw;
          }
    };
    
    template <class SCALAR>
    class Gmsh_Matrix
    {
    private:
      int r,c;
    public:
      inline int size1() const {return r;}
      inline int size2() const {return c;}
      SCALAR *data;
      ~Gmsh_Matrix() {delete [] data;}
      Gmsh_Matrix(int R,int C)
        : r(R),c(C)
      {
        data = new SCALAR [r*c];
        scal(0);
      }
      Gmsh_Matrix(const Gmsh_Matrix<SCALAR> &other)
        : r(other.r),c(other.c)
      {
        data = new double [r*c];
        copy ( other.data ) ;
      }
      inline SCALAR operator () (int i, int j) const
      {
        return data[i+r*j];
      }
      inline SCALAR & operator () (int i, int j)
      {
        return data[i+r*j];
      }
      inline Gmsh_Matrix operator *(const Gmsh_Matrix<SCALAR> & other)
      {
        throw;
      }
      inline void scal ( const SCALAR s)
      {
        for (int i=0;i<r*c;++i)data[i]*=s;
      }
      inline void copy ( const SCALAR **other)
      {
        for (int i=0;i<r*c;++i)data[i]=other.data[i];
      }
      inline void mult(const Gmsh_Matrix<SCALAR> & x, const Gmsh_Matrix<SCALAR> & b)
      {
        throw;
      }
      inline void mult (const Gmsh_Vector<SCALAR> & x, Gmsh_Vector<SCALAR> & b )
      {
        throw;
      }
      inline void least_squares (const Gmsh_Vector<SCALAR> & rhs, Gmsh_Vector<SCALAR> & result)
      {
        throw;
      }
      inline void lu_solve (const Gmsh_Vector<SCALAR> & rhs, Gmsh_Vector<SCALAR> & result)
      {
        throw;
      }
    };
    
    #ifdef HAVE_GSL
    #include <gsl/gsl_linalg.h>
    #include <gsl/gsl_matrix.h>
    #include <gsl/gsl_blas.h>
    class GSL_Vector
    {
    private:
      int r;
    public:
      inline int size() const {return r;}
      gsl_vector *data;
      ~GSL_Vector() {gsl_vector_free (data);}
      GSL_Vector(int R)
        : r(R)
      {
        data = gsl_vector_calloc (r);
      }
      GSL_Vector(const GSL_Vector&other)
        : r(other.r)
      {
        data = gsl_vector_calloc (r);
        gsl_vector_memcpy (data, other.data);
      }
      inline double operator () (int i) const
      {
        return gsl_vector_get (data,i);
      }
      inline double & operator () (int i)
      {
        return *gsl_vector_ptr (data,i);
      }
    };
    
    class GSL_Matrix
    {
    private:
      int r,c;
    public:
      inline int size1() const {return r;}
      inline int size2() const {return c;}
      gsl_matrix *data;
      ~GSL_Matrix() {gsl_matrix_free (data);}
      GSL_Matrix(int R,int C)
        : r(R),c(C)
      {
        data = gsl_matrix_calloc (r, c);
      }
      GSL_Matrix(const GSL_Matrix&other)
        : r(other.r),c(other.c)
      {
        data = gsl_matrix_calloc (r, c);
        gsl_matrix_memcpy (data, other.data);
      }
      inline double operator () (int i, int j) const
      {
        return gsl_matrix_get (data,i,j);
      }
      inline double & operator () (int i, int j)
      {
        return *gsl_matrix_ptr (data,i,j);
      }
      inline void mult(const GSL_Matrix & x, const GSL_Matrix & b)
      {
        gsl_blas_dgemm (CblasNoTrans,CblasNoTrans, 1.0, data, x.data, 1.0, b.data);
      }
    
      inline void least_squares (const GSL_Vector & rhs, GSL_Vector & result)
          {
    	  assert (r > c);
    	  assert (rhs.size() == r);
    	  assert (result.size() == c);
    	  GSL_Matrix *ls     = new GSL_Matrix(c, c);
    	  GSL_Vector *ls_rhs = new GSL_Vector(c);
    	  //GSL_Vector *test = new GSL_Vector(c);
    	  gsl_blas_dgemm (CblasTrans,CblasNoTrans, 1.0, data, data, 1.0, ls->data);
    	  gsl_blas_dgemv (CblasTrans, 1.0, data, rhs.data, 1.0, ls_rhs->data);
    	  ls->lu_solve (*ls_rhs,result);
    	  
    
    	  delete ls;
    	  delete ls_rhs;
          }
      inline void lu_solve (const GSL_Vector & rhs, GSL_Vector & result)
          {
    	  int s;
    	  gsl_permutation * p = gsl_permutation_alloc (size1());
    	  gsl_linalg_LU_decomp ( data, p, &s);
    	  gsl_linalg_LU_solve ( data ,  p, rhs.data, result.data ) ;
    	  gsl_permutation_free (p);
          }
    
      inline void mult (const GSL_Vector & x, GSL_Vector & b )
      {
        gsl_blas_dgemv (CblasNoTrans, 1.0, data, x.data, 1.0, b.data);
      }
    };
    typedef GSL_Matrix Double_Matrix;
    typedef GSL_Vector Double_Vector;
    #else
    typedef Gmsh_Matrix<double> Double_Matrix;
    typedef Gmsh_Vector<double> Double_Vector;
    #endif
    typedef Gmsh_Matrix<int>    Int_Matrix;
    typedef Gmsh_Vector<int>    Int_Vector;
    #endif