From 8b57824cb621e986a3763a711652fbc6c9c66b1c Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sun, 21 Dec 2008 23:07:08 +0000
Subject: [PATCH] fix compile
---
Numeric/GmshMatrix.h | 251 ++++++++++++++++++++++---------------------
1 file changed, 126 insertions(+), 125 deletions(-)
diff --git a/Numeric/GmshMatrix.h b/Numeric/GmshMatrix.h
index 4b60e129cb..7922f6a041 100644
--- a/Numeric/GmshMatrix.h
+++ b/Numeric/GmshMatrix.h
@@ -9,46 +9,48 @@
#include <math.h>
#include "GmshMessage.h"
+template <class SCALAR> class Gmsh_Matrix;
+
template <class SCALAR>
class Gmsh_Vector
{
private:
- int r;
- SCALAR *data;
- friend class Gmsh_Matrix;
+ int _r;
+ SCALAR *_data;
+ friend class Gmsh_Matrix<SCALAR>;
public:
- inline int size() const { return r; }
- Gmsh_Vector(int R) : r(R)
+ Gmsh_Vector(int r) : _r(r)
{
- data = new SCALAR[r];
+ _data = new SCALAR[_r];
scale(0.);
}
- Gmsh_Vector(const Gmsh_Vector<SCALAR> &other) : r(other.r)
+ Gmsh_Vector(const Gmsh_Vector<SCALAR> &other) : _r(other._r)
{
- data = new SCALAR[r];
- for(int i = 0; i < r; ++i) data[i] = other.data[i];
+ _data = new SCALAR[_r];
+ for(int i = 0; i < _r; ++i) _data[i] = other._data[i];
}
- ~Gmsh_Vector() { delete [] data; }
+ ~Gmsh_Vector() { if(_data) delete [] _data; }
inline SCALAR operator () (int i) const
{
- return data[i];
+ return _data[i];
}
+ inline int size() const { return _r; }
inline SCALAR & operator () (int i)
{
- return data[i];
+ return _data[i];
}
inline SCALAR norm()
{
SCALAR n = 0.;
- for(int i = 0; i < r; ++i) n += data[i] * data[i];
+ for(int i = 0; i < _r; ++i) n += _data[i] * _data[i];
return sqrt(n);
}
inline void scale(const SCALAR s)
{
if(s == 0.)
- for(int i = 0; i < r; ++i) data[i] = 0.;
+ for(int i = 0; i < _r; ++i) _data[i] = 0.;
else
- for(int i = 0; i < r; ++i) data[i] *= s;
+ for(int i = 0; i < _r; ++i) _data[i] *= s;
}
};
@@ -56,13 +58,13 @@ template <class SCALAR>
class Gmsh_Matrix
{
private:
- int r, c;
- SCALAR *data;
+ int _r, _c;
+ SCALAR *_data;
void _back_substitution(int *indx, SCALAR *b)
{
int i, ii = -1, ip, j;
SCALAR sum;
- for(i = 0; i < c; i++){
+ for(i = 0; i < _c; i++){
ip = indx[i];
sum = b[ip];
b[ip] = b[i];
@@ -71,38 +73,38 @@ class Gmsh_Matrix
else if(sum) ii = i;
b[i] = sum;
}
- for(i = c - 1; i >= 0; i--){
+ for(i = _c - 1; i >= 0; i--){
sum = b[i];
- for(j = i + 1; j < c; j++) sum -= (*this)(i, j) * b[j];
+ for(j = i + 1; j < _c; j++) sum -= (*this)(i, j) * b[j];
b[i] = sum / (*this)(i, i);
}
}
bool _lu_decomposition(int *indx , SCALAR &determinant)
{
- if(r != c)
- Msg::Fatal("Gmsh_Matrix::_lu_decomposition : cannot lu decompose a non-square matrix");
+ if(_r != _c)
+ Msg::Fatal("Cannot compute lu factorization of non-square matrix");
int i, imax, j, k;
SCALAR big, dum, sum, temp;
- SCALAR *vv = new SCALAR[c];
+ SCALAR *vv = new SCALAR[_c];
determinant = 1.0;
- for(i = 0; i < c; i++){
- big = 0.0;
- for(j = 0; j < c; j++)
- if((temp=fabs((*this)(i, j))) > big) big = temp;
- if(big == 0.0) {
+ for(i = 0; i < _c; i++){
+ big = 0.;
+ for(j = 0; j < _c; j++)
+ if((temp = fabs((*this)(i, j))) > big) big = temp;
+ if(big == 0.) {
return false;
big = 1.e-12;
}
vv[i] = 1.0 / big;
}
- for(j = 0; j < c; j++){
+ for(j = 0; j < _c; j++){
for(i = 0; i < j; i++){
- sum=(*this)(i, j);
+ sum = (*this)(i, j);
for(k = 0; k < i; k++) sum -= (*this)(i, k)*(*this)(k, j);
(*this)(i, j) = sum;
}
- big = 0.0;
- for(i = j; i < c; i++){
+ big = 0.;
+ for(i = j; i < _c; i++){
sum = (*this)(i, j);
for(k = 0; k < j; k++)
sum -= (*this)(i, k) * (*this)(k, j);
@@ -113,7 +115,7 @@ class Gmsh_Matrix
}
}
if(j != imax){
- for(k = 0; k < c; k++){
+ for(k = 0; k < _c; k++){
dum = (*this)(imax, k);
(*this)(imax, k) = (*this)(j, k);
(*this)(j, k) = dum;
@@ -122,59 +124,59 @@ class Gmsh_Matrix
vv[imax] = vv[j];
}
indx[j] = imax;
- if((*this)(j, j) == 0.0) (*this)(j, j) = 1.e-20;
- if(j != c){
+ if((*this)(j, j) == 0.) (*this)(j, j) = 1.e-20;
+ if(j != _c){
dum = 1.0 / ((*this)(j, j));
- for(i = j + 1; i < c; i++) (*this)(i, j) *= dum;
+ for(i = j + 1; i < _c; i++) (*this)(i, j) *= dum;
}
}
delete [] vv;
return true;
}
public:
- inline int size1() const { return r; }
- inline int size2() const { return c; }
- Gmsh_Matrix(int R, int C) : r(R), c(C)
+ Gmsh_Matrix(int r, int c) : _r(r), _c(c)
{
- data = new SCALAR[r * c];
+ _data = new SCALAR[_r * _c];
scale(0.);
}
- Gmsh_Matrix(const Gmsh_Matrix<SCALAR> &other) : r(other.r), c(other.c)
+ Gmsh_Matrix(const Gmsh_Matrix<SCALAR> &other) : _r(other._r), _c(other._c)
{
- data = new SCALAR[r * c];
+ _data = new SCALAR[_r * _c];
memcpy(other);
}
- ~Gmsh_Matrix() { delete [] data; }
- Gmsh_Matrix & operator=(const Gmsh_Matrix<SCALAR> &other)
+ Gmsh_Matrix() : _r(0), _c(0), _data(0) {}
+ ~Gmsh_Matrix() { if(_data) delete [] _data; }
+ inline int size1() const { return _r; }
+ inline int size2() const { return _c; }
+ Gmsh_Matrix<SCALAR> & operator=(const Gmsh_Matrix<SCALAR> &other)
{
if(this != &other){
- r = other.r;
- c = other.c;
- data = new SCALAR[r * c];
+ _r = other._r;
+ _c = other._c;
+ _data = new SCALAR[_r * _c];
memcpy(other);
}
return *this;
}
- Gmsh_Matrix() : r(0), c(0), data(0) {}
- void memcpy(const Gmsh_Matrix &other)
+ void memcpy(const Gmsh_Matrix<SCALAR> &other)
{
- for(int i = 0; i < r * c; ++i) data[i] = other.data[i];
+ for(int i = 0; i < _r * _c; ++i) _data[i] = other._data[i];
}
inline SCALAR operator () (int i, int j) const
{
- return data[i + r * j];
+ return _data[i + _r * j];
}
inline SCALAR & operator () (int i, int j)
{
- return data[i + r * j];
+ return _data[i + _r * j];
}
inline void mult(const Gmsh_Matrix<SCALAR> &b, Gmsh_Matrix<SCALAR> &c)
{
c.scale(0.);
- for(int i = 0; i < r; i++)
+ for(int i = 0; i < _r; i++)
for(int j = 0; j < b.size2(); j++)
- for(int k = 0; k < c; k++)
- c.data[i + r * j] += (*this)(i, k) * b(k, j);
+ for(int k = 0; k < _c; k++)
+ c._data[i + _r * j] += (*this)(i, k) * b(k, j);
}
inline void blas_dgemm(const Gmsh_Matrix<SCALAR> &b, Gmsh_Matrix<SCALAR> &c,
const SCALAR alpha=1.0, const SCALAR beta=1.0)
@@ -187,19 +189,19 @@ class Gmsh_Matrix
}
inline void set_all(const SCALAR &m)
{
- for(int i = 0; i < r * c; i++) data[i] = m;
+ for(int i = 0; i < _r * _c; i++) _data[i] = m;
}
inline void lu_solve(const Gmsh_Vector<SCALAR> &rhs, Gmsh_Vector<SCALAR> &result)
{
- int *indx = new int [c];
+ int *indx = new int[_c];
SCALAR d;
if(!_lu_decomposition(indx, d))
- Msg::Fatal("Singular matrix in Gmsh_Matrix::_lu_decomposition");
- for(int i = 0; i < c; i++) result[i] = rhs[i];
- _back_substitution(indx, result.data);
+ Msg::Fatal("LU fatorization failed (singular matrix)");
+ for(int i = 0; i < _c; i++) result(i) = rhs(i);
+ _back_substitution(indx, result._data);
delete [] indx;
}
- Gmsh_Matrix cofactor(int i, int j) const
+ Gmsh_Matrix<SCALAR> cofactor(int i, int j) const
{
int ni = size1();
int nj = size2();
@@ -215,39 +217,39 @@ class Gmsh_Matrix
inline void mult(const Gmsh_Vector<SCALAR> &x, Gmsh_Vector<SCALAR> &y)
{
y.scale(0.);
- for(int i = 0; i < r; i++)
- for(int j = 0; j < c; j++)
- y.data[i] += (*this)(i, j) * x(j);
+ for(int i = 0; i < _r; i++)
+ for(int j = 0; j < _c; j++)
+ y._data[i] += (*this)(i, j) * x(j);
}
SCALAR determinant() const
{
- Gmsh_Matrix copy = *this;
+ Gmsh_Matrix<SCALAR> copy = *this;
SCALAR factor = 1.0;
- int *indx = new int[c];
- if(!copy._lu_decomposition(indx, factor)) return 0.0;
+ int *indx = new int[_c];
+ if(!copy._lu_decomposition(indx, factor)) return 0.;
SCALAR det = factor;
- for(int i = 0; i < c; i++) det *= copy(i, i);
+ for(int i = 0; i < _c; i++) det *= copy(i, i);
delete [] indx;
return det;
}
- inline Gmsh_Matrix touchSubmatrix(int i0, int ni, int j0, int nj)
+ inline Gmsh_Matrix<SCALAR> touchSubmatrix(int i0, int ni, int j0, int nj)
{
Msg::Fatal("Gmsh_Matrix::touchSubmatrix is not implemented");
- Gmsh_Matrix subm(ni, nj);
+ Gmsh_Matrix<SCALAR> subm(ni, nj);
return subm;
}
inline void scale(const double s)
{
if(s == 0.)
- for(int i = 0; i < r * c; ++i) data[i] = 0.;
+ for(int i = 0; i < _r * _c; ++i) _data[i] = 0.;
else
- for(int i = 0; i < r * c; ++i) data[i] *= s;
+ for(int i = 0; i < _r * _c; ++i) _data[i] *= s;
}
inline void add(const double &a)
{
- for(int i = 0; i < r * c; ++i) data[i] += a;
+ for(int i = 0; i < _r * _c; ++i) _data[i] += a;
}
- inline void add(const Gmsh_Matrix &m)
+ inline void add(const Gmsh_Matrix<SCALAR> &m)
{
for(int i = 0; i < size1(); i++)
for(int j = 0; j < size2(); j++)
@@ -262,105 +264,107 @@ class Gmsh_Matrix
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
+class GSL_Matrix;
+
class GSL_Vector
{
private:
- int r;
- gsl_vector *data;
+ int _r;
+ gsl_vector *_data;
friend class GSL_Matrix;
public:
- inline int size() const { return r; }
- GSL_Vector(int R) : r(R)
+ GSL_Vector(int r) : _r(r)
{
- data = gsl_vector_calloc(r);
+ _data = gsl_vector_calloc(_r);
}
- GSL_Vector(const GSL_Vector &other) : r(other.r)
+ GSL_Vector(const GSL_Vector &other) : _r(other._r)
{
- data = gsl_vector_calloc(r);
- gsl_vector_memcpy(data, other.data);
+ _data = gsl_vector_calloc(_r);
+ gsl_vector_memcpy(_data, other._data);
}
- ~GSL_Vector() { gsl_vector_free(data); }
+ ~GSL_Vector() { gsl_vector_free(_data); }
+ inline int size() const { return _r; }
inline double operator () (int i) const
{
- return gsl_vector_get(data, i);
+ return gsl_vector_get(_data, i);
}
inline double & operator () (int i)
{
- return *gsl_vector_ptr(data, i);
+ return *gsl_vector_ptr(_data, i);
}
inline double norm()
{
- return gsl_blas_dnrm2(data);
+ return gsl_blas_dnrm2(_data);
}
inline void scale(const double s)
{
- if(s == 0.) gsl_vector_set_zero(data);
- else gsl_vector_scale(data, s);
+ if(s == 0.) gsl_vector_set_zero(_data);
+ else gsl_vector_scale(_data, s);
}
};
class GSL_Matrix
{
private:
- gsl_matrix_view view;
- gsl_matrix *data;
+ gsl_matrix_view _view;
+ gsl_matrix *_data;
inline const gsl_matrix_view _see_submatrix(int i0, int ni, int j0, int nj) const
{
- return gsl_matrix_submatrix(data, i0, j0, ni, nj);
+ return gsl_matrix_submatrix(_data, i0, j0, ni, nj);
}
public:
- inline int size1() const { return data->size1; }
- inline int size2() const { return data->size2; }
- GSL_Matrix(gsl_matrix_view _data) : view(_data), data(&view.matrix) {}
- GSL_Matrix(int r, int c) { data = gsl_matrix_calloc(r, c); }
- GSL_Matrix() : data(0) {}
- GSL_Matrix(const GSL_Matrix &other) : data(0)
- {
- if(data) gsl_matrix_free(data);
- data = gsl_matrix_calloc(other.data->size1, other.data->size2);
- gsl_matrix_memcpy(data, other.data);
- }
- ~GSL_Matrix() { if(data && data->owner == 1) gsl_matrix_free(data); }
+ GSL_Matrix(gsl_matrix_view view) : _view(view), _data(&_view.matrix) {}
+ GSL_Matrix(int r, int c) { _data = gsl_matrix_calloc(r, c); }
+ GSL_Matrix(const GSL_Matrix &other) : _data(0)
+ {
+ if(_data) gsl_matrix_free(_data);
+ _data = gsl_matrix_calloc(other._data->size1, other._data->size2);
+ gsl_matrix_memcpy(_data, other._data);
+ }
+ GSL_Matrix() : _data(0) {}
+ ~GSL_Matrix() { if(_data && _data->owner == 1) gsl_matrix_free(_data); }
+ inline int size1() const { return _data->size1; }
+ inline int size2() const { return _data->size2; }
GSL_Matrix & operator = (const GSL_Matrix&other)
{
if(&other != this){
- if(data) gsl_matrix_free(data);
- data = gsl_matrix_calloc(other.data->size1, other.data->size2);
- gsl_matrix_memcpy(data, other.data);
+ if(_data) gsl_matrix_free(_data);
+ _data = gsl_matrix_calloc(other._data->size1, other._data->size2);
+ gsl_matrix_memcpy(_data, other._data);
}
return *this;
}
void memcpy(const GSL_Matrix &other)
{
- gsl_matrix_memcpy(data, other.data);
+ gsl_matrix_memcpy(_data, other._data);
}
inline double operator () (int i, int j) const
{
- return gsl_matrix_get(data, i, j);
+ return gsl_matrix_get(_data, i, j);
}
inline double & operator () (int i, int j)
{
- return *gsl_matrix_ptr(data, i, j);
+ return *gsl_matrix_ptr(_data, i, j);
}
inline void mult(const GSL_Matrix &b, GSL_Matrix &c)
{
- gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1., data, b.data, 0., c.data);
+ gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1., _data, b._data, 0., c._data);
}
inline void blas_dgemm(const GSL_Matrix &x, GSL_Matrix &b,
const double alpha = 1.0, const double beta = 1.0)
{
- gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, alpha, x.data, b.data, beta, data);
+ gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, alpha, x._data, b._data, beta, _data);
}
inline void set_all(const double &m)
{
- gsl_matrix_set_all(data, m);
+ gsl_matrix_set_all(_data, m);
}
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_linalg_LU_decomp(_data, p, &s);
+ gsl_linalg_LU_solve(_data, p, rhs._data, result._data);
gsl_permutation_free(p);
}
GSL_Matrix cofactor(int i, int j) const
@@ -388,33 +392,33 @@ class GSL_Matrix
}
inline void mult(const GSL_Vector &x, GSL_Vector &b)
{
- gsl_blas_dgemv(CblasNoTrans, 1., data, x.data, 0., b.data);
+ gsl_blas_dgemv(CblasNoTrans, 1., _data, x._data, 0., b._data);
}
double determinant() const
{
GSL_Matrix copy = *this;
gsl_permutation *p = gsl_permutation_alloc(size1());
int s;
- gsl_linalg_LU_decomp(copy.data, p, &s);
+ gsl_linalg_LU_decomp(copy._data, p, &s);
gsl_permutation_free(p);
- return gsl_linalg_LU_det(copy.data, s);
+ return gsl_linalg_LU_det(copy._data, s);
}
inline gsl_matrix_view touchSubmatrix(int i0, int ni, int j0, int nj)
{
- return gsl_matrix_submatrix(data, i0, j0, ni, nj);
+ return gsl_matrix_submatrix(_data, i0, j0, ni, nj);
}
inline void scale(const double s)
{
- if(s == 0.) gsl_matrix_set_zero(data);
- else gsl_matrix_scale(data, s);
+ if(s == 0.) gsl_matrix_set_zero(_data);
+ else gsl_matrix_scale(_data, s);
}
inline void add(const double &a)
{
- gsl_matrix_add_constant(data, a);
+ gsl_matrix_add_constant(_data, a);
}
inline void add(const GSL_Matrix &m)
{
- gsl_matrix_add(data, m.data);
+ gsl_matrix_add(_data, m._data);
}
};
@@ -428,7 +432,4 @@ typedef Gmsh_Vector<double> Double_Vector;
#endif
-typedef Gmsh_Matrix<int> Int_Matrix;
-typedef Gmsh_Vector<int> Int_Vector;
-
#endif
--
GitLab