Select Git revision
meshGRegionDelaunayInsertion.cpp
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