pp

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");
 }

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);
 };