diff --git a/Numeric/fullMatrix.cpp b/Numeric/fullMatrix.cpp
index ea9d2a66611cde989e501ed06ab516e2ebe547ad..7d2349efc364c56dfb0674e8b3b0aa9a1622364a 100644
--- a/Numeric/fullMatrix.cpp
+++ b/Numeric/fullMatrix.cpp
@@ -17,6 +17,16 @@
#define F77NAME(x) (x##_)
#endif
+
+
+/*==============================================================================
+ * This file improve the methods of fullVector and fullMatrix using
+ * - BLAS : Basic Linear Algebra Subprograms
+ * - LAPACK : Linear Algebra PACKage
+ *============================================================================*/
+
+
+
#if defined(HAVE_BLAS)
extern "C" {
@@ -130,6 +140,8 @@ void fullMatrix<std::complex<double> >::mult(const fullVector<std::complex<doubl
#endif
+
+
#if defined(HAVE_LAPACK)
extern "C" {
@@ -307,6 +319,12 @@ bool fullMatrix<double>::svd(fullMatrix<double> &V, fullVector<double> &S)
#endif
+
+
+/*==============================================================================
+ * BINDINGS
+ *============================================================================*/
+
#include "Bindings.h"
template<>
@@ -317,6 +335,9 @@ void fullMatrix<double>::registerBindings(binding *b)
"The memory is allocated in one continuous block and stored "
"in column major order (like in fortran).");
methodBinding *cm;
+ cm = cb->setConstructor<fullMatrix<double>,int,int>();
+ cm->setDescription ("A new matrix of size 'nRows' x 'nColumns'");
+ cm->setArgNames("nRows","nColumns",NULL);
cm = cb->addMethod("size1", &fullMatrix<double>::size1);
cm->setDescription("Returns the number of rows in the matrix");
cm = cb->addMethod("size2", &fullMatrix<double>::size2);
@@ -340,9 +361,6 @@ void fullMatrix<double>::registerBindings(binding *b)
cm = cb->addMethod("print", &fullMatrix<double>::print);
cm->setArgNames("name",NULL);
cm->setDescription("print the matrix");
- cm = cb->setConstructor<fullMatrix<double>,int,int>();
- cm->setDescription ("A new matrix of size 'nRows' x 'nColumns'");
- cm->setArgNames("nRows","nColumns",NULL);
cm = cb->addMethod("invertInPlace", &fullMatrix<double>::invertInPlace);
cm->setDescription("invert the matrix and return the determinant");
}
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index fecedd92165b7d0d97ebeecb77002d5475e767be..1eec0f1d4557872ebe89ac1166b2365025b5e394 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -14,32 +14,60 @@
class binding;
template <class scalar> class fullMatrix;
+
+
+/*==============================================================================
+ * class fullVector : An abstract interface for vectors of scalar
+ *============================================================================*/
+
template <class scalar>
class fullVector
{
+
private:
- int _r;
- scalar *_data;
+
+ int _r; // Size of the vector
+ scalar *_data; // Pointer on the first element
friend class fullMatrix<scalar>;
+
public:
- inline const scalar* getDataPtr() const{
- return _data;
- }
+
+ // Constructor and destructor
+
+ fullVector(void) : _r(0),_data(0) {}
fullVector(int r) : _r(r)
{
_data = new scalar[_r];
setAll(0.);
}
- fullVector(void) : _r(0),_data(0) {}
fullVector(const fullVector<scalar> &other) : _r(other._r)
{
_data = new scalar[_r];
for(int i = 0; i < _r; ++i) _data[i] = other._data[i];
}
- ~fullVector() { if(_data) delete [] _data; }
- bool resize(int r)
+ ~fullVector()
+ {
+ if(_data) delete [] _data;
+ }
+
+ // Get information (size, value)
+
+ inline int size() const { return _r; }
+ inline const scalar * getDataPtr() const { return _data; }
+ inline scalar operator () (int i) const { return _data[i]; }
+ inline scalar & operator () (int i) { return _data[i]; }
+
+ // Operations
+
+ inline scalar norm() const
+ {
+ scalar n = 0.;
+ for(int i = 0; i < _r; ++i) n += _data[i] * _data[i];
+ return sqrt(n);
+ }
+ bool resize(int r)
{
- if (_r < r){
+ if (_r < r) {
if (_data) delete[] _data;
_r = r;
_data = new scalar[_r];
@@ -49,26 +77,10 @@ class fullVector
}
void setAsProxy(const fullVector<scalar> &original, int r_start, int r)
{
- if(_data)
- delete [] _data;
+ if(_data) delete [] _data;
_r = r;
_data = original._data + r_start;
}
- inline scalar operator () (int i) const
- {
- return _data[i];
- }
- inline int size() const { return _r; }
- inline scalar & operator () (int i)
- {
- return _data[i];
- }
- inline scalar norm() const
- {
- scalar n = 0.;
- for(int i = 0; i < _r; ++i) n += _data[i] * _data[i];
- return sqrt(n);
- }
inline void scale(const scalar s)
{
if(s == 0.)
@@ -103,6 +115,9 @@ class fullVector
{
for (int i = 0; i < _r; i++) _data[i] *= x._data[i];
}
+
+ // Printing and file treatment
+
void print(const char *name="") const
{
printf("Printing vector %s:\n", name);
@@ -112,46 +127,36 @@ class fullVector
}
printf("\n");
}
- void binarySave (FILE *f) const{
+ void binarySave (FILE *f) const
+ {
fwrite (_data, sizeof(scalar), _r, f);
}
- void binaryLoad (FILE *f){
+ void binaryLoad (FILE *f)
+ {
if(fread (_data, sizeof(scalar), _r, f) != _r) return;
}
};
+
+
+/*==============================================================================
+ * class fullMatrix : An abstract interface for matrix of scalar
+ *============================================================================*/
+
template <class scalar>
class fullMatrix
{
+
private:
- bool _own_data; // should data be freed on delete ?
- int _r, _c;
- scalar *_data;
+
+ bool _own_data; // should data be freed on delete ?
+ int _r, _c; // Size of the matrix
+ scalar *_data; // Pointer on the first element
+
public:
- inline scalar get(int r, int c) const {
- #ifdef _DEBUG
- if (r >= _r || r < 0 || c >= _c || c < 0)
- Msg::Fatal("invalid index to access fullMatrix : %i %i (size = %i %i)",
- r, c, _r, _c);
- #endif
- return (*this)(r, c);
- }
- inline void set(int r, int c, scalar v){
- #ifdef _DEBUG
- if (r >= _r || r < 0 || c >= _c || c < 0)
- Msg::Fatal("invalid index to access fullMatrix : %i %i (size = %i %i)",
- r, c, _r, _c);
- #endif
- (*this)(r, c) = v;
- }
- inline scalar norm() const
- {
- scalar n = 0.;
- for(int i = 0; i < _r; ++i)
- for(int j = 0; j < _c; ++j)
- n += (*this)(i, j) * (*this)(i, j);
- return sqrt(n);
- }
+
+ // Constructor and destructor
+
fullMatrix(scalar *original, int r, int c)
{
_r = r;
@@ -184,10 +189,45 @@ class fullMatrix
for(int i = 0; i < _r * _c; ++i) _data[i] = other._data[i];
}
fullMatrix() : _own_data(false),_r(0), _c(0), _data(0) {}
- ~fullMatrix() { if(_data && _own_data) delete [] _data; }
+ ~fullMatrix()
+ {
+ if(_data && _own_data) delete [] _data;
+ }
+
+ // Get information (size, value)
+
+ inline int size1() const { return _r; }
+ inline int size2() const { return _c; }
+ inline scalar get(int r, int c) const {
+ #ifdef _DEBUG
+ if (r >= _r || r < 0 || c >= _c || c < 0)
+ Msg::Fatal("invalid index to access fullMatrix : %i %i (size = %i %i)",
+ r, c, _r, _c);
+ #endif
+ return (*this)(r, c);
+ }
+
+ // Operations
+
+ inline void set(int r, int c, scalar v){
+ #ifdef _DEBUG
+ if (r >= _r || r < 0 || c >= _c || c < 0)
+ Msg::Fatal("invalid index to access fullMatrix : %i %i (size = %i %i)",
+ r, c, _r, _c);
+ #endif
+ (*this)(r, c) = v;
+ }
+ inline scalar norm() const
+ {
+ scalar n = 0.;
+ for(int i = 0; i < _r; ++i)
+ for(int j = 0; j < _c; ++j)
+ n += (*this)(i, j) * (*this)(i, j);
+ return sqrt(n);
+ }
bool resize(int r, int c, bool resetValue = true) // data will be owned (same as constructor)
{
- if ((r * c > _r * _c) || !_own_data){
+ if ((r * c > _r * _c) || !_own_data) {
_r = r;
_c = c;
if (_own_data && _data) delete[] _data;
@@ -197,7 +237,7 @@ class fullMatrix
setAll(0.);
return true;
}
- else{
+ else {
_r = r;
_c = c;
}
@@ -234,8 +274,6 @@ class fullMatrix
_own_data = false;
_data = original._data;
}
- inline int size1() const { return _r; }
- inline int size2() const { return _c; }
fullMatrix<scalar> & operator = (const fullMatrix<scalar> &other)
{
if(this != &other){
@@ -276,8 +314,7 @@ class fullMatrix
#endif
return _data[i + _r * j];
}
- void copy(const fullMatrix<scalar> &a, int i0, int ni, int j0, int nj,
- int desti0, int destj0)
+ void copy(const fullMatrix<scalar> &a, int i0, int ni, int j0, int nj, int desti0, int destj0)
{
for(int i = i0, desti = desti0; i < i0 + ni; i++, desti++)
for(int j = j0, destj = destj0; j < j0 + nj; j++, destj++)
@@ -293,7 +330,7 @@ class fullMatrix
_own_data = true;
setAll(a);
}
- void mult_naive(const fullMatrix<scalar> &b, fullMatrix<scalar> &c)const
+ void mult_naive(const fullMatrix<scalar> &b, fullMatrix<scalar> &c) const
{
c.scale(0.);
for(int i = 0; i < _r; i++)
@@ -301,8 +338,7 @@ class fullMatrix
for(int k = 0; k < _c; k++)
c._data[i + _r * j] += (*this)(i, k) * b(k, j);
}
- ;
- void mult(const fullMatrix<scalar> &b, fullMatrix<scalar> &c)const
+ void mult(const fullMatrix<scalar> &b, fullMatrix<scalar> &c) const
#if !defined(HAVE_BLAS)
{
mult_naive(b,c);
@@ -313,8 +349,7 @@ class fullMatrix
{
for (int i = 0; i < _r * _c; i++) _data[i] *= a._data[i];
}
- void gemm_naive(const fullMatrix<scalar> &a, const fullMatrix<scalar> &b,
- scalar alpha=1., scalar beta=1.)
+ void gemm_naive(const fullMatrix<scalar> &a, const fullMatrix<scalar> &b, scalar alpha=1., scalar beta=1.)
{
fullMatrix<scalar> temp(a.size1(), b.size2());
a.mult_naive(b, temp);
@@ -322,9 +357,7 @@ class fullMatrix
scale(beta);
add(temp);
}
- ;
- void gemm(const fullMatrix<scalar> &a, const fullMatrix<scalar> &b,
- scalar alpha=1., scalar beta=1.)
+ void gemm(const fullMatrix<scalar> &a, const fullMatrix<scalar> &b, scalar alpha=1., scalar beta=1.)
#if !defined(HAVE_BLAS)
{
gemm_naive(a,b,alpha,beta);
@@ -341,7 +374,7 @@ class fullMatrix
}
void scale(const double s)
#if !defined(HAVE_BLAS)
- {
+ {
if(s == 0.) // this is not really correct nan*0 (or inf*0) is expected to give nan
for(int i = 0; i < _r * _c; ++i) _data[i] = 0.;
else
@@ -435,12 +468,10 @@ class fullMatrix
int ni = size1();
int nj = size2();
fullMatrix<scalar> cof(ni - 1, nj - 1);
- for(int I = 0; I < ni; I++){
- for(int J = 0; J < nj; J++){
+ for(int I = 0; I < ni; I++)
+ for(int J = 0; J < nj; J++)
if(J != j && I != i)
cof(I < i ? I : I - 1, J < j ? J : J - 1) = (*this)(I, J);
- }
- }
return cof;
}
scalar determinant() const
@@ -459,6 +490,9 @@ class fullMatrix
}
#endif
;
+
+ // Printing
+
void print(const std::string name = "") const
{
printf("Printing matrix %s:\n", name.c_str());
@@ -472,6 +506,9 @@ class fullMatrix
printf("\n");
}
}
+
+ // Bindings
+
static void registerBindings(binding *b);
};