getfem-commits
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Getfem-commits] (no subject)


From: Konstantinos Poulios
Subject: [Getfem-commits] (no subject)
Date: Wed, 16 Jan 2019 08:09:58 -0500 (EST)

branch: integration-point-variables
commit 5c0db7f3ab476c0d8bc3877480f45e1bba426909
Author: Konstantinos Poulios <address@hidden>
Date:   Wed Jan 16 14:09:42 2019 +0100

    Code cleanup and small extension of bgeot::tensor constructors
---
 src/Makefile.am                                 |   2 +-
 src/getfem/bgeot_tensor.h                       |  24 +-
 src/getfem_generic_assembly_compile_and_exec.cc |  71 +--
 src/gmm/gmm_vector.h                            | 570 ++++++++++++------------
 4 files changed, 336 insertions(+), 331 deletions(-)

diff --git a/src/Makefile.am b/src/Makefile.am
index d7e1a61..ee36640 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -209,7 +209,7 @@ SRC =                                               \
        getfem_mat_elem.cc                              \
        getfem_mat_elem_type.cc                         \
        getfem_level_set.cc                                 \
-    getfem_locale.cc                        \
+       getfem_locale.cc                        \
        getfem_mesh_level_set.cc                        \
        getfem_mesh_im_level_set.cc                     \
        getfem_mesh_fem_level_set.cc                    \
diff --git a/src/getfem/bgeot_tensor.h b/src/getfem/bgeot_tensor.h
index 43d779d..815f6eb 100644
--- a/src/getfem/bgeot_tensor.h
+++ b/src/getfem/bgeot_tensor.h
@@ -62,7 +62,7 @@ namespace bgeot {
       } else resize(1);
     }
 
-    void reset(void) { std::fill(begin(), end(), 0); }
+    void reset() { std::fill(begin(), end(), 0); }
 
     inline bool finished(const multi_index &m) {
       if (m.size() == 0)
@@ -83,7 +83,7 @@ namespace bgeot {
       : std::vector<size_type>(4)
     { (*this)[0] = i; (*this)[1] = j; (*this)[2] = k; (*this)[3] = l; }
 
-    multi_index(void) {}
+    multi_index() {}
 
     bool is_equal(const multi_index &m) const {
       if (this->size() != m.size()) return false;
@@ -92,7 +92,7 @@ namespace bgeot {
       return true;
     }
 
-    size_type total_size(void) const {
+    size_type total_size() const {
       size_type s = 1;
       for (size_type k = 0; k < this->size(); ++k) s *= (*this)[k];
       return s;
@@ -194,10 +194,10 @@ namespace bgeot {
       return *(this->begin() + d);
     }
 
-    inline size_type size(void) const { return std::vector<T>::size(); }
+    inline size_type size() const { return std::vector<T>::size(); }
     inline size_type size(size_type i) const { return sizes_[i]; }
-    inline const multi_index &sizes(void) const { return sizes_; }
-    inline size_type order(void) const { return sizes_.size(); }
+    inline const multi_index &sizes() const { return sizes_; }
+    inline size_type order() const { return sizes_.size(); }
 
     void init(const multi_index &c) {
       auto it = c.begin();
@@ -236,7 +236,7 @@ namespace bgeot {
     }
 
     inline void adjust_sizes(const multi_index &mi) { init(mi); }
-    inline void adjust_sizes(void) { init(); }
+    inline void adjust_sizes() { init(); }
     inline void adjust_sizes(size_type i) { init(i); }
     inline void adjust_sizes(size_type i, size_type j) { init(i, j); }
     inline void adjust_sizes(size_type i, size_type j, size_type k)
@@ -294,8 +294,8 @@ namespace bgeot {
         + sizeof(*this) + sizes_.memsize() + coeff.memsize();
     }
 
-    std::vector<T> &as_vector(void) { return *this; }
-    const std::vector<T> &as_vector(void) const { return *this; }
+    std::vector<T> &as_vector() { return *this; }
+    const std::vector<T> &as_vector() const { return *this; }
 
 
     tensor<T>& operator +=(const tensor<T>& w)
@@ -323,9 +323,13 @@ namespace bgeot {
     }
 
     tensor(const multi_index &c) { init(c); }
+    tensor(size_type i, size_type j)
+    { init(multi_index(i, j)); }
+    tensor(size_type i, size_type j, size_type k)
+    { init(multi_index(i, j, k)); }
     tensor(size_type i, size_type j, size_type k, size_type l)
     { init(multi_index(i, j, k, l)); }
-    tensor(void) {}
+    tensor() {}
   };
 
   template<class T> void tensor<T>::mat_transp_reduction
diff --git a/src/getfem_generic_assembly_compile_and_exec.cc 
b/src/getfem_generic_assembly_compile_and_exec.cc
index 5c32f28..4a2e2b0 100644
--- a/src/getfem_generic_assembly_compile_and_exec.cc
+++ b/src/getfem_generic_assembly_compile_and_exec.cc
@@ -4135,7 +4135,7 @@ namespace getfem {
   };
 
   template <class MAT>
-  inline void add_elem_matrix_
+  inline void add_elem_matrix
   (MAT &K, const std::vector<size_type> &dofs1,
    const std::vector<size_type> &dofs2, std::vector<size_type> &/*dofs1_sort*/,
    base_vector &elem, scalar_type threshold, size_type /* N */) {
@@ -4155,7 +4155,7 @@ namespace getfem {
   //   return  int((*the_indto_sort)[aa]) - int((*the_indto_sort)[bb]);
   // }
 
-  inline void add_elem_matrix_
+  inline void add_elem_matrix
   (gmm::col_matrix<gmm::rsvector<scalar_type>> &K,
    const std::vector<size_type> &dofs1, const std::vector<size_type> &dofs2,
    std::vector<size_type> &dofs1_sort,
@@ -4318,12 +4318,12 @@ namespace getfem {
 
         if (pmf1 == pmf2 && cv1 == cv2) {
           if (ifirst1 == ifirst2) {
-            add_elem_matrix_(K, dofs1, dofs1, dofs1_sort, elem, ninf*1E-14, N);
+            add_elem_matrix(K, dofs1, dofs1, dofs1_sort, elem, ninf*1E-14, N);
           } else {
             dofs2.resize(dofs1.size());
             for (size_type i = 0; i < dofs1.size(); ++i)
               dofs2[i] =  dofs1[i] + ifirst2 - ifirst1;
-            add_elem_matrix_(K, dofs1, dofs2, dofs1_sort, elem, ninf*1E-14, N);
+            add_elem_matrix(K, dofs1, dofs2, dofs1_sort, elem, ninf*1E-14, N);
           }
         } else {
           dofs2.assign(s2, ifirst2);
@@ -4345,7 +4345,7 @@ namespace getfem {
           } else
             for (size_type i=0; i < s2; ++i) dofs2[i] += i;
 
-          add_elem_matrix_(K, dofs1, dofs2, dofs1_sort, elem, ninf*1E-14, N);
+          add_elem_matrix(K, dofs1, dofs2, dofs1_sort, elem, ninf*1E-14, N);
         }
       }
       return 0;
@@ -4423,12 +4423,12 @@ namespace getfem {
 
         if (pmf2 == pmf1 && cv1 == cv2) {
           if (I1.first() == I2.first()) {
-            add_elem_matrix_(K, dofs1, dofs1, dofs1_sort, elem, ninf*1E-14, N);
+            add_elem_matrix(K, dofs1, dofs1, dofs1_sort, elem, ninf*1E-14, N);
           } else {
             dofs2.resize(dofs1.size());
             for (size_type i = 0; i < dofs1.size(); ++i)
               dofs2[i] =  dofs1[i] + I2.first() - I1.first();
-            add_elem_matrix_(K, dofs1, dofs2, dofs1_sort, elem, ninf*1E-14, N);
+            add_elem_matrix(K, dofs1, dofs2, dofs1_sort, elem, ninf*1E-14, N);
           }
         } else {
           if (cv2 == size_type(-1)) return 0;
@@ -4437,7 +4437,7 @@ namespace getfem {
           dofs2.resize(ct2.size());
           for (size_type i = 0; i < ct2.size(); ++i)
             dofs2[i] = ct2[i] + I2.first();
-          add_elem_matrix_(K, dofs1, dofs2, dofs1_sort, elem, ninf*1E-14, N);
+          add_elem_matrix(K, dofs1, dofs2, dofs1_sort, elem, ninf*1E-14, N);
         }
       }
       return 0;
@@ -4518,12 +4518,12 @@ namespace getfem {
 
         if (pmf2 == pmf1 && cv1 == cv2) {
           if (I1.first() == I2.first()) {
-            add_elem_matrix_(K, dofs1, dofs1, dofs1_sort, elem, ninf*1E-14, N);
+            add_elem_matrix(K, dofs1, dofs1, dofs1_sort, elem, ninf*1E-14, N);
           } else {
             dofs2.resize(dofs1.size());
             for (size_type i = 0; i < dofs1.size(); ++i)
               dofs2[i] =  dofs1[i] + I2.first() - I1.first();
-            add_elem_matrix_(K, dofs1, dofs2, dofs1_sort, elem, ninf*1E-14, N);
+            add_elem_matrix(K, dofs1, dofs2, dofs1_sort, elem, ninf*1E-14, N);
           }
         } else {
           if (cv2 == size_type(-1)) return 0;
@@ -4536,7 +4536,7 @@ namespace getfem {
             for (size_type q = 0; q < qmult2; ++q)
               *itd++ += *itt + q;
 
-          add_elem_matrix_(K, dofs1, dofs2, dofs1_sort, elem, ninf*1E-14, N);
+          add_elem_matrix(K, dofs1, dofs2, dofs1_sort, elem, ninf*1E-14, N);
         }
       }
       return 0;
@@ -4604,26 +4604,26 @@ namespace getfem {
 
         if (pmf2 == pmf1 && cv1 == cv2) {
           if (i1 == i2) {
-            add_elem_matrix_(K, dofs1, dofs1, dofs1_sort, elem, ninf, N);
+            add_elem_matrix(K, dofs1, dofs1, dofs1_sort, elem, ninf, N);
             for (size_type i = 0; i < ss1; ++i) (dofs1[i])++;
-            add_elem_matrix_(K, dofs1, dofs1, dofs1_sort, elem, ninf, N);
+            add_elem_matrix(K, dofs1, dofs1, dofs1_sort, elem, ninf, N);
           } else {
             dofs2.resize(ss2);
             for (size_type i = 0; i < ss2; ++i) dofs2[i] = i2 + ct1[i];
-            add_elem_matrix_(K, dofs1, dofs2, dofs1_sort, elem, ninf, N);
+            add_elem_matrix(K, dofs1, dofs2, dofs1_sort, elem, ninf, N);
             for (size_type i = 0; i < ss1; ++i) (dofs1[i])++;
             for (size_type i = 0; i < ss2; ++i) (dofs2[i])++;
-            add_elem_matrix_(K, dofs1, dofs2, dofs1_sort, elem, ninf, N);
+            add_elem_matrix(K, dofs1, dofs2, dofs1_sort, elem, ninf, N);
           }
         } else {
           if (cv2 == size_type(-1)) return 0;
           auto &ct2 = pmf2->ind_scalar_basic_dof_of_element(cv2);
           dofs2.resize(ss2);
           for (size_type i = 0; i < ss2; ++i) dofs2[i] = i2 + ct2[i];
-          add_elem_matrix_(K, dofs1, dofs2, dofs1_sort, elem, ninf, N);
+          add_elem_matrix(K, dofs1, dofs2, dofs1_sort, elem, ninf, N);
           for (size_type i = 0; i < ss1; ++i) (dofs1[i])++;
           for (size_type i = 0; i < ss2; ++i) (dofs2[i])++;
-          add_elem_matrix_(K, dofs1, dofs2, dofs1_sort, elem, ninf, N);
+          add_elem_matrix(K, dofs1, dofs2, dofs1_sort, elem, ninf, N);
         }
       }
       return 0;
@@ -4691,34 +4691,34 @@ namespace getfem {
 
         if (pmf2 == pmf1 && cv1 == cv2) {
           if (i1 == i2) {
-            add_elem_matrix_(K, dofs1, dofs1, dofs1_sort, elem, ninf, N);
+            add_elem_matrix(K, dofs1, dofs1, dofs1_sort, elem, ninf, N);
             for (size_type i = 0; i < ss1; ++i) (dofs1[i])++;
-            add_elem_matrix_(K, dofs1, dofs1, dofs1_sort, elem, ninf, N);
+            add_elem_matrix(K, dofs1, dofs1, dofs1_sort, elem, ninf, N);
             for (size_type i = 0; i < ss1; ++i) (dofs1[i])++;
-            add_elem_matrix_(K, dofs1, dofs1, dofs1_sort, elem, ninf, N);
+            add_elem_matrix(K, dofs1, dofs1, dofs1_sort, elem, ninf, N);
           } else {
             dofs2.resize(ss2);
             for (size_type i = 0; i < ss2; ++i) dofs2[i] = i2 + ct1[i];
-            add_elem_matrix_(K, dofs1, dofs2, dofs1_sort, elem, ninf, N);
+            add_elem_matrix(K, dofs1, dofs2, dofs1_sort, elem, ninf, N);
             for (size_type i = 0; i < ss1; ++i) (dofs1[i])++;
             for (size_type i = 0; i < ss2; ++i) (dofs2[i])++;
-            add_elem_matrix_(K, dofs1, dofs2, dofs1_sort, elem, ninf, N);
+            add_elem_matrix(K, dofs1, dofs2, dofs1_sort, elem, ninf, N);
             for (size_type i = 0; i < ss1; ++i) (dofs1[i])++;
             for (size_type i = 0; i < ss2; ++i) (dofs2[i])++;
-            add_elem_matrix_(K, dofs1, dofs2, dofs1_sort, elem, ninf, N);
+            add_elem_matrix(K, dofs1, dofs2, dofs1_sort, elem, ninf, N);
           }
         } else {
           if (cv2 == size_type(-1)) return 0;
           auto &ct2 = pmf2->ind_scalar_basic_dof_of_element(cv2);
           dofs2.resize(ss2);
           for (size_type i = 0; i < ss2; ++i) dofs2[i] = i2 + ct2[i];
-          add_elem_matrix_(K, dofs1, dofs2, dofs1_sort, elem, ninf, N);
+          add_elem_matrix(K, dofs1, dofs2, dofs1_sort, elem, ninf, N);
           for (size_type i = 0; i < ss1; ++i) (dofs1[i])++;
           for (size_type i = 0; i < ss2; ++i) (dofs2[i])++;
-          add_elem_matrix_(K, dofs1, dofs2, dofs1_sort, elem, ninf, N);
+          add_elem_matrix(K, dofs1, dofs2, dofs1_sort, elem, ninf, N);
           for (size_type i = 0; i < ss1; ++i) (dofs1[i])++;
           for (size_type i = 0; i < ss2; ++i) (dofs2[i])++;
-          add_elem_matrix_(K, dofs1, dofs2, dofs1_sort, elem, ninf, N);
+          add_elem_matrix(K, dofs1, dofs2, dofs1_sort, elem, ninf, N);
         }
       }
       return 0;
@@ -6578,7 +6578,7 @@ namespace getfem {
       }
     }
     rmi.node_list[pnode->hash_value].push_back(pnode);
-  }
+  } // ga_compile_node
 
   void ga_compile_function(ga_workspace &workspace,
                            ga_instruction_set &gis, bool scalar) {
@@ -6857,18 +6857,19 @@ namespace getfem {
                   GMM_ASSERT1(root->tensor_proper_size() == 1,
                               "Invalid vector or tensor quantity. An order 2 "
                               "weak form has to be a scalar quantity");
-                  const mesh_fem 
*mf1=workspace.associated_mf(root->name_test1);
-                  const mesh_fem 
*mf2=workspace.associated_mf(root->name_test2);
-                  const mesh_fem **mfg1 = 0, **mfg2 = 0;
+                  const mesh_fem 
+                    *mf1=workspace.associated_mf(root->name_test1),
+                    *mf2=workspace.associated_mf(root->name_test2),
+                    **mfg1 = 0, **mfg2 = 0;
                   const im_data
                     *imd1 = workspace.associated_im_data(root->name_test1),
                     *imd2 = workspace.associated_im_data(root->name_test2);
-                  const std::string &intn1 = root->interpolate_name_test1;
-                  const std::string &intn2 = root->interpolate_name_test2;
+                  const std::string &intn1 = root->interpolate_name_test1,
+                                    &intn2 = root->interpolate_name_test2;
                   bool secondary1 = intn1.size() &&
-                      workspace.secondary_domain_exists(intn1);
+                                    workspace.secondary_domain_exists(intn1);
                   bool secondary2 = intn2.size() &&
-                      workspace.secondary_domain_exists(intn2);
+                                    workspace.secondary_domain_exists(intn2);
                   fem_interpolation_context
                     &ctx1 = intn1.size()
                           ? (secondary1 ? rmi.secondary_domain_infos.ctx
@@ -6962,7 +6963,7 @@ namespace getfem {
                 }
               }
               if (pgai)
-                
gis.all_instructions[rm].instructions.push_back(std::move(pgai));
+                rmi.instructions.push_back(std::move(pgai));
             }
           }
         }
diff --git a/src/gmm/gmm_vector.h b/src/gmm/gmm_vector.h
index e69931d..25481db 100644
--- a/src/gmm/gmm_vector.h
+++ b/src/gmm/gmm_vector.h
@@ -53,7 +53,7 @@ namespace gmm {
 
     V *pm;
     size_type l;
-    
+
     public :
 
     operator T() const { return pm->r(l); }
@@ -92,7 +92,7 @@ namespace gmm {
 
     V *pm;
     size_type l;
-    
+
     public :
 
     operator std::complex<T>() const { return pm->r(l); }
@@ -139,21 +139,21 @@ namespace gmm {
     { return std::complex<T>(*this)* v; }
     std::complex<T> operator /(std::complex<T> v)
     { return std::complex<T>(*this)/ v; }
-  };  
+  };
+
 
-  
   template<typename T, typename V> inline
   bool operator ==(T v, const ref_elt_vector<T, V> &re) { return (v==T(re)); }
   template<typename T, typename V> inline
   bool operator !=(T v, const ref_elt_vector<T, V> &re) { return (v!=T(re)); }
   template<typename T, typename V> inline
-  T &operator +=(T &v, const ref_elt_vector<T, V> &re) 
+  T &operator +=(T &v, const ref_elt_vector<T, V> &re)
   { v += T(re); return v; }
   template<typename T, typename V> inline
   T &operator -=(T &v, const ref_elt_vector<T, V> &re)
   { v -= T(re); return v; }
   template<typename T, typename V> inline
-  T &operator *=(T &v, const ref_elt_vector<T, V> &re) 
+  T &operator *=(T &v, const ref_elt_vector<T, V> &re)
   { v *= T(re); return v; }
   template<typename T, typename V> inline
   T &operator /=(T &v, const ref_elt_vector<T, V> &re)
@@ -224,7 +224,7 @@ namespace gmm {
     size_type i;    // Current index.
     T* p;           // Pointer to the current position.
     dsvector<T> *v; // Pointer to the vector.
-    
+
     typedef T                   value_type;
     typedef value_type*         pointer;
     typedef const value_type*   const_pointer;
@@ -233,20 +233,20 @@ namespace gmm {
     typedef ptrdiff_t           difference_type;
     typedef std::bidirectional_iterator_tag iterator_category;
     typedef dsvector_iterator<T> iterator;
-    
+
     reference operator *() const { return *p; }
     pointer operator->() const { return &(operator*()); }
 
     iterator &operator ++() {
       for (size_type k = (i & 15); k < 15; ++k)
-       { ++p; ++i; if (*p != T(0)) return *this; }
+        { ++p; ++i; if (*p != T(0)) return *this; }
       v->next_pos(*(const_cast<const_pointer *>(&(p))), i);
       return *this;
     }
     iterator operator ++(int) { iterator tmp = *this; ++(*this); return tmp; }
     iterator &operator --() {
       for (size_type k = (i & 15); k > 0; --k)
-       { --p; --i; if (*p != T(0)) return *this; }
+        { --p; --i; if (*p != T(0)) return *this; }
       v->previous_pos(p, i);
       return *this;
     }
@@ -256,10 +256,10 @@ namespace gmm {
     { return (i == it.i && p == it.p && v == it.v); }
     bool operator !=(const iterator &it) const
     { return !(it == *this); }
-    
-    size_type index(void) const { return i; }
 
-    dsvector_iterator(void) : i(size_type(-1)), p(0), v(0) {}
+    size_type index() const { return i; }
+
+    dsvector_iterator() : i(size_type(-1)), p(0), v(0) {}
     dsvector_iterator(dsvector<T> &w) : i(size_type(-1)), p(0), v(&w) {};
   };
 
@@ -268,7 +268,7 @@ namespace gmm {
     size_type i;          // Current index.
     const T* p;           // Pointer to the current position.
     const dsvector<T> *v; // Pointer to the vector.
-    
+
     typedef T                   value_type;
     typedef const value_type*   pointer;
     typedef const value_type&   reference;
@@ -276,19 +276,19 @@ namespace gmm {
     typedef ptrdiff_t           difference_type;
     typedef std::bidirectional_iterator_tag iterator_category;
     typedef dsvector_const_iterator<T> iterator;
-   
+
     reference operator *() const { return *p; }
     pointer operator->() const { return &(operator*()); }
     iterator &operator ++() {
       for (size_type k = (i & 15); k < 15; ++k)
-       { ++p; ++i; if (*p != T(0)) return *this; }
+        { ++p; ++i; if (*p != T(0)) return *this; }
       v->next_pos(p, i);
       return *this;
     }
     iterator operator ++(int) { iterator tmp = *this; ++(*this); return tmp; }
     iterator &operator --() {
       for (size_type k = (i & 15); k > 0; --k)
-       { --p; --i; if (*p != T(0)) return *this; }
+        { --p; --i; if (*p != T(0)) return *this; }
       v->previous_pos(p, i);
       return *this;
     }
@@ -298,17 +298,17 @@ namespace gmm {
     { return (i == it.i && p == it.p && v == it.v); }
     bool operator !=(const iterator &it) const
     { return !(it == *this); }
-    
-    size_type index(void) const { return i; }
 
-    dsvector_const_iterator(void) : i(size_type(-1)), p(0) {}
+    size_type index() const { return i; }
+
+    dsvector_const_iterator() : i(size_type(-1)), p(0) {}
     dsvector_const_iterator(const dsvector_iterator<T> &it)
       : i(it.i), p(it.p), v(it.v) {}
     dsvector_const_iterator(const dsvector<T> &w)
       : i(size_type(-1)), p(0), v(&w) {};
   };
 
-  
+
   /**
      Sparse vector built on distribution sort principle.
      Read and write access have a constant complexity depending only on the
@@ -323,7 +323,7 @@ namespace gmm {
     typedef const T *                  const_pointer;
     typedef void *                     void_pointer;
     typedef const void *               const_void_pointer;
- 
+
   protected:
     size_type    n;         // Potential vector size
     size_type    depth;     // Number of row of pointer arrays
@@ -337,10 +337,10 @@ namespace gmm {
       void_pointer p = root_ptr;
       if (!p) return 0;
       for (size_type k = 0; k < depth; ++k) {
-       p = ((void **)(p))[(i & my_mask) >> my_shift];
-       if (!p) return 0;
-       my_mask = (my_mask >> 4);
-       my_shift -= 4;
+        p = ((void **)(p))[(i & my_mask) >> my_shift];
+        if (!p) return 0;
+        my_mask = (my_mask >> 4);
+        my_shift -= 4;
       }
       GMM_ASSERT1(my_shift == 0, "internal error");
       GMM_ASSERT1(my_mask == 15, "internal error");
@@ -351,32 +351,32 @@ namespace gmm {
       GMM_ASSERT1(i < n, "index " << i << " out of range (size " << n << ")");
       size_type my_mask = mask, my_shift = shift;
       if (!root_ptr) {
-       if (depth) {
-         root_ptr = new void_pointer[16];
-         std::memset(root_ptr, 0, 16*sizeof(void_pointer));
-       } else {
-         root_ptr = new T[16];
-         for (size_type l = 0; l < 16; ++l) ((T *)(root_ptr))[l] = T(0);
-       }
+        if (depth) {
+          root_ptr = new void_pointer[16];
+          std::memset(root_ptr, 0, 16*sizeof(void_pointer));
+        } else {
+          root_ptr = new T[16];
+          for (size_type l = 0; l < 16; ++l) ((T *)(root_ptr))[l] = T(0);
+        }
       }
 
       void_pointer p = root_ptr;
       for (size_type k = 0; k < depth; ++k) {
-       size_type j = (i & my_mask) >> my_shift;
-       void_pointer q = ((void_pointer *)(p))[j];
-       if (!q) {
-         if (k+1 != depth) {
-           q = new void_pointer[16];
-           std::memset(q, 0, 16*sizeof(void_pointer));
-         } else {
-           q = new T[16];
-           for (size_type l = 0; l < 16; ++l) ((T *)(q))[l] = T(0);
-         }
-         ((void_pointer *)(p))[j] = q;
-       }
-       p = q;
-       my_mask = (my_mask >> 4);
-       my_shift -= 4;
+        size_type j = (i & my_mask) >> my_shift;
+        void_pointer q = ((void_pointer *)(p))[j];
+        if (!q) {
+          if (k+1 != depth) {
+            q = new void_pointer[16];
+            std::memset(q, 0, 16*sizeof(void_pointer));
+          } else {
+            q = new T[16];
+            for (size_type l = 0; l < 16; ++l) ((T *)(q))[l] = T(0);
+          }
+          ((void_pointer *)(p))[j] = q;
+        }
+        p = q;
+        my_mask = (my_mask >> 4);
+        my_shift -= 4;
       }
       GMM_ASSERT1(my_shift == 0, "internal error");
       GMM_ASSERT1(my_mask == 15, "internal error " << my_mask);
@@ -392,65 +392,65 @@ namespace gmm {
 
     void rec_del(void_pointer p, size_type my_depth) {
       if (my_depth) {
-       for (size_type k = 0; k < 16; ++k)
-         if (((void_pointer *)(p))[k])
-           rec_del(((void_pointer *)(p))[k], my_depth-1);
-       delete[] ((void_pointer *)(p));
+        for (size_type k = 0; k < 16; ++k)
+          if (((void_pointer *)(p))[k])
+            rec_del(((void_pointer *)(p))[k], my_depth-1);
+        delete[] ((void_pointer *)(p));
       } else {
-       delete[] ((T *)(p));
+        delete[] ((T *)(p));
       }
     }
 
     void rec_clean(void_pointer p, size_type my_depth, double eps) {
       if (my_depth) {
-       for (size_type k = 0; k < 16; ++k)
-         if (((void_pointer *)(p))[k])
-           rec_clean(((void_pointer *)(p))[k], my_depth-1, eps);
+        for (size_type k = 0; k < 16; ++k)
+          if (((void_pointer *)(p))[k])
+            rec_clean(((void_pointer *)(p))[k], my_depth-1, eps);
       } else {
-       for (size_type k = 0; k < 16; ++k)
-         if (gmm::abs(((T *)(p))[k]) <= eps) ((T *)(p))[k] = T(0);
+        for (size_type k = 0; k < 16; ++k)
+          if (gmm::abs(((T *)(p))[k]) <= eps) ((T *)(p))[k] = T(0);
       }
     }
 
     void rec_clean_i(void_pointer p, size_type my_depth, size_type my_mask,
-                    size_type i, size_type base) {
+                     size_type i, size_type base) {
       if (my_depth) {
-       my_mask = (my_mask >> 4);
-       for (size_type k = 0; k < 16; ++k)
-         if (((void_pointer *)(p))[k] && (base + (k+1)*(mask+1)) >= i)
-           rec_clean_i(((void_pointer *)(p))[k], my_depth-1, my_mask,
-                       i, base + k*(my_mask+1));
+        my_mask = (my_mask >> 4);
+        for (size_type k = 0; k < 16; ++k)
+          if (((void_pointer *)(p))[k] && (base + (k+1)*(mask+1)) >= i)
+            rec_clean_i(((void_pointer *)(p))[k], my_depth-1, my_mask,
+                        i, base + k*(my_mask+1));
       } else {
-       for (size_type k = 0; k < 16; ++k)
-         if (base+k > i) ((T *)(p))[k] = T(0);
+        for (size_type k = 0; k < 16; ++k)
+          if (base+k > i) ((T *)(p))[k] = T(0);
       }
     }
- 
-      
+
+
     size_type rec_nnz(void_pointer p, size_type my_depth) const {
       size_type nn = 0;
       if (my_depth) {
-       for (size_type k = 0; k < 16; ++k)
-         if (((void_pointer *)(p))[k])
-           nn += rec_nnz(((void_pointer *)(p))[k], my_depth-1);
+        for (size_type k = 0; k < 16; ++k)
+          if (((void_pointer *)(p))[k])
+            nn += rec_nnz(((void_pointer *)(p))[k], my_depth-1);
       } else {
-       for (size_type k = 0; k < 16; ++k)
-         if (((const T *)(p))[k] != T(0)) nn++;
+        for (size_type k = 0; k < 16; ++k)
+          if (((const T *)(p))[k] != T(0)) nn++;
       }
       return nn;
     }
 
     void copy_rec(void_pointer &p, const_void_pointer q, size_type my_depth) {
       if (my_depth) {
-       p = new void_pointer[16];
-       std::memset(p, 0, 16*sizeof(void_pointer));
-       for (size_type l = 0; l < 16; ++l)
-         if (((const const_void_pointer *)(q))[l])
-           copy_rec(((void_pointer *)(p))[l],
-                    ((const const_void_pointer *)(q))[l], my_depth-1);
+        p = new void_pointer[16];
+        std::memset(p, 0, 16*sizeof(void_pointer));
+        for (size_type l = 0; l < 16; ++l)
+          if (((const const_void_pointer *)(q))[l])
+            copy_rec(((void_pointer *)(p))[l],
+                     ((const const_void_pointer *)(q))[l], my_depth-1);
       } else {
-       p = new T[16];
-       for (size_type l = 0; l < 16; ++l) ((T *)(p))[l] = ((const T *)(q))[l];
+        p = new T[16];
+        for (size_type l = 0; l < 16; ++l) ((T *)(p))[l] = ((const T *)(q))[l];
       }
     }
 
@@ -462,75 +462,75 @@ namespace gmm {
     }
 
     void next_pos_rec(void_pointer p, size_type my_depth, size_type my_mask,
-                     const_pointer &pp, size_type &i, size_type base) const {
+                      const_pointer &pp, size_type &i, size_type base) const {
       size_type ii = i;
       if (my_depth) {
-       my_mask = (my_mask >> 4);
-       for (size_type k = 0; k < 16; ++k)
-         if (((void_pointer *)(p))[k] && (base + (k+1)*(my_mask+1)) >= i) {
-           next_pos_rec(((void_pointer *)(p))[k], my_depth-1, my_mask,
-                        pp, i, base + k*(my_mask+1));
-           if (i != size_type(-1)) return; else i = ii;
-       }
-       i = size_type(-1); pp = 0;
+        my_mask = (my_mask >> 4);
+        for (size_type k = 0; k < 16; ++k)
+          if (((void_pointer *)(p))[k] && (base + (k+1)*(my_mask+1)) >= i) {
+            next_pos_rec(((void_pointer *)(p))[k], my_depth-1, my_mask,
+                         pp, i, base + k*(my_mask+1));
+            if (i != size_type(-1)) return; else i = ii;
+        }
+        i = size_type(-1); pp = 0;
       } else {
-       for (size_type k = 0; k < 16; ++k)
-         if (base+k > i && ((const_pointer)(p))[k] != T(0))
-           { i = base+k; pp = &(((const_pointer)(p))[k]); return; }
-       i = size_type(-1); pp = 0;
+        for (size_type k = 0; k < 16; ++k)
+          if (base+k > i && ((const_pointer)(p))[k] != T(0))
+            { i = base+k; pp = &(((const_pointer)(p))[k]); return; }
+        i = size_type(-1); pp = 0;
       }
     }
 
     void previous_pos_rec(void_pointer p, size_type my_depth, size_type 
my_mask,
-                         const_pointer &pp, size_type &i,
-                         size_type base) const {
+                          const_pointer &pp, size_type &i,
+                          size_type base) const {
       size_type ii = i;
       if (my_depth) {
-       my_mask = (my_mask >> 4);
-       for (size_type k = 15; k != size_type(-1); --k)
-         if (((void_pointer *)(p))[k] && ((base + k*(my_mask+1)) < i)) {
-           previous_pos_rec(((void_pointer *)(p))[k], my_depth-1,
-                            my_mask, pp, i, base + k*(my_mask+1));
-           if (i != size_type(-1)) return; else i = ii;
-       }
-       i = size_type(-1); pp = 0;
+        my_mask = (my_mask >> 4);
+        for (size_type k = 15; k != size_type(-1); --k)
+          if (((void_pointer *)(p))[k] && ((base + k*(my_mask+1)) < i)) {
+            previous_pos_rec(((void_pointer *)(p))[k], my_depth-1,
+                             my_mask, pp, i, base + k*(my_mask+1));
+            if (i != size_type(-1)) return; else i = ii;
+        }
+        i = size_type(-1); pp = 0;
       } else {
-       for (size_type k = 15; k != size_type(-1); --k)
-         if (base+k < i && ((const_pointer)(p))[k] != T(0))
-           { i = base+k; pp = &(((const_pointer)(p))[k]); return; }
-       i = size_type(-1); pp = 0;
+        for (size_type k = 15; k != size_type(-1); --k)
+          if (base+k < i && ((const_pointer)(p))[k] != T(0))
+            { i = base+k; pp = &(((const_pointer)(p))[k]); return; }
+        i = size_type(-1); pp = 0;
       }
     }
-    
-    
+
+
   public:
     void clean(double eps) { if (root_ptr) rec_clean(root_ptr, depth); }
     void resize(size_type n_) {
       if (n_ != n) {
-       n = n_;
-       if (n_ < n) { // Depth unchanged (a choice)
-         if (root_ptr) rec_clean_i(root_ptr, depth, mask, n_, 0);
-       } else {
-         // may change the depth (add some levels)
-         size_type my_depth = 0, my_shift = 0, my_mask = 1; if (n_) --n_;
-         while (n_) { n_ /= 16; ++my_depth; my_shift += 4; my_mask *= 16; }
-         my_mask--; if (my_shift) my_shift -= 4; if (my_depth) --my_depth;
-         if (my_depth > depth || depth == 0) {
-           if (root_ptr) {
-             for (size_type k = depth; k < my_depth; ++k) {
-               void_pointer *q = new void_pointer [16];
-               std::memset(q, 0, 16*sizeof(void_pointer));
-               q[0] = root_ptr; root_ptr = q;
-             }
-           }
-           mask = my_mask; depth = my_depth; shift = my_shift;
-         }
-       }
+        n = n_;
+        if (n_ < n) { // Depth unchanged (a choice)
+          if (root_ptr) rec_clean_i(root_ptr, depth, mask, n_, 0);
+        } else {
+          // may change the depth (add some levels)
+          size_type my_depth = 0, my_shift = 0, my_mask = 1; if (n_) --n_;
+          while (n_) { n_ /= 16; ++my_depth; my_shift += 4; my_mask *= 16; }
+          my_mask--; if (my_shift) my_shift -= 4; if (my_depth) --my_depth;
+          if (my_depth > depth || depth == 0) {
+            if (root_ptr) {
+              for (size_type k = depth; k < my_depth; ++k) {
+                void_pointer *q = new void_pointer [16];
+                std::memset(q, 0, 16*sizeof(void_pointer));
+                q[0] = root_ptr; root_ptr = q;
+              }
+            }
+            mask = my_mask; depth = my_depth; shift = my_shift;
+          }
+        }
       }
     }
-    
-    void clear(void) { if (root_ptr) rec_del(root_ptr, depth); root_ptr = 0; }
-    
+
+    void clear() { if (root_ptr) rec_del(root_ptr, depth); root_ptr = 0; }
+
     void next_pos(const_pointer &pp, size_type &i) const {
       if (!root_ptr || i >= n) { pp = 0, i = size_type(-1); return; }
       next_pos_rec(root_ptr, depth, mask, pp, i, 0);
@@ -542,29 +542,29 @@ namespace gmm {
       previous_pos_rec(root_ptr, depth, mask, pp, i, 0);
     }
 
-    iterator begin(void) {
-      iterator it(*this); 
+    iterator begin() {
+      iterator it(*this);
       if (n && root_ptr) {
-       it.i = 0; it.p = const_cast<T *>(read_access(0));
-       if (!(it.p) || *(it.p) == T(0))
-         next_pos(*(const_cast<const_pointer *>(&(it.p))), it.i);
+        it.i = 0; it.p = const_cast<T *>(read_access(0));
+        if (!(it.p) || *(it.p) == T(0))
+          next_pos(*(const_cast<const_pointer *>(&(it.p))), it.i);
       }
       return it;
     }
 
-    iterator end(void) { return iterator(*this); }
+    iterator end() { return iterator(*this); }
 
-    const_iterator begin(void) const {
+    const_iterator begin() const {
       const_iterator it(*this);
       if (n && root_ptr) {
-       it.i = 0; it.p = read_access(0);
-       if (!(it.p) || *(it.p) == T(0)) next_pos(it.p, it.i);
+        it.i = 0; it.p = read_access(0);
+        if (!(it.p) || *(it.p) == T(0)) next_pos(it.p, it.i);
       }
       return it;
     }
 
-    const_iterator end(void) const { return const_iterator(*this); }
-    
+    const_iterator end() const { return const_iterator(*this); }
+
     inline ref_elt_vector<T, dsvector<T> > operator [](size_type c)
     { return ref_elt_vector<T, dsvector<T> >(this, c); }
 
@@ -580,22 +580,22 @@ namespace gmm {
     { const T *p = read_access(c); if (p) return *p; else return T(0); }
 
     inline T operator [](size_type c) const { return r(c); }
-    
-    size_type nnz(void) const
+
+    size_type nnz() const
     { if (root_ptr) return rec_nnz(root_ptr, depth); else return 0; }
-    size_type size(void) const { return n; }
+    size_type size() const { return n; }
 
     void swap(dsvector<T> &v) {
       std::swap(n, v.n); std::swap(root_ptr, v.root_ptr);
       std::swap(depth, v.depth); std::swap(shift, v.shift);
       std::swap(mask, v.mask);
     }
-    
+
     /* Constructors */
     dsvector(const dsvector<T> &v) { init(0); copy(v); }
     dsvector<T> &operator =(const dsvector<T> &v) { copy(v); return *this; }
     explicit dsvector(size_type l){ init(l); }
-    dsvector(void) { init(0); }
+    dsvector() { init(0); }
     ~dsvector() { if (root_ptr) rec_del(root_ptr, depth); root_ptr = 0; }
   };
 
@@ -621,10 +621,10 @@ namespace gmm {
     { o->clear(); }
     static void do_clear(this_type &v) { v.clear(); }
     static value_type access(const origin_type *o, const const_iterator &,
-                            const const_iterator &, size_type i)
+                             const const_iterator &, size_type i)
     { return (*o)[i]; }
     static reference access(origin_type *o, const iterator &, const iterator &,
-                           size_type i)
+                            size_type i)
     { return (*o)[i]; }
     static void resize(this_type &v, size_type n) { v.resize(n); }
   };
@@ -635,12 +635,12 @@ namespace gmm {
   /******* Optimized operations for dsvector<T> ****************************/
 
   template <typename T> inline void copy(const dsvector<T> &v1,
-                                        dsvector<T> &v2) {
+                                          dsvector<T> &v2) {
     GMM_ASSERT2(v1.size() == v2.size(), "dimensions mismatch");
     v2 = v1;
   }
   template <typename T> inline void copy(const dsvector<T> &v1,
-                                        const dsvector<T> &v2) {
+                                         const dsvector<T> &v2) {
     GMM_ASSERT2(v1.size() == v2.size(), "dimensions mismatch");
     v2 = const_cast<dsvector<T> &>(v1);
   }
@@ -655,23 +655,23 @@ namespace gmm {
   }
   template <typename T> inline
   void copy(const simple_vector_ref<const dsvector<T> *> &v1,
-           dsvector<T> &v2)
+            dsvector<T> &v2)
   { copy(*(v1.origin), v2); }
   template <typename T> inline
   void copy(const simple_vector_ref<dsvector<T> *> &v1, dsvector<T> &v2)
   { copy(*(v1.origin), v2); }
   template <typename T> inline
   void copy(const simple_vector_ref<dsvector<T> *> &v1,
-           const simple_vector_ref<dsvector<T> *> &v2)
+            const simple_vector_ref<dsvector<T> *> &v2)
   { copy(*(v1.origin), v2); }
   template <typename T> inline
   void copy(const simple_vector_ref<const dsvector<T> *> &v1,
-           const simple_vector_ref<dsvector<T> *> &v2)
+            const simple_vector_ref<dsvector<T> *> &v2)
   { copy(*(v1.origin), v2); }
-  
+
   template <typename T>
   inline size_type nnz(const dsvector<T>& l) { return l.nnz(); }
-  
+
   /*************************************************************************/
   /*                                                                       */
   /* Class wsvector: sparse vector optimized for random write operations,  */
@@ -679,7 +679,7 @@ namespace gmm {
   /* Based on std::map                                                     */
   /*                                                                       */
   /*************************************************************************/
-  
+
   template<typename T> struct wsvector_iterator
     : public std::map<size_type, T>::iterator {
     typedef typename std::map<size_type, T>::iterator base_it_type;
@@ -689,12 +689,12 @@ namespace gmm {
     // typedef size_t              size_type;
     typedef ptrdiff_t           difference_type;
     typedef std::bidirectional_iterator_tag iterator_category;
-    
+
     reference operator *() const { return (base_it_type::operator*()).second; }
     pointer operator->() const { return &(operator*()); }
-    size_type index(void) const { return (base_it_type::operator*()).first; }
+    size_type index() const { return (base_it_type::operator*()).first; }
 
-    wsvector_iterator(void) {}
+    wsvector_iterator() {}
     wsvector_iterator(const base_it_type &it) : base_it_type(it) {}
   };
 
@@ -707,12 +707,12 @@ namespace gmm {
     // typedef size_t              size_type;
     typedef ptrdiff_t           difference_type;
     typedef std::bidirectional_iterator_tag iterator_category;
-    
+
     reference operator *() const { return (base_it_type::operator*()).second; }
     pointer operator->() const { return &(operator*()); }
-    size_type index(void) const { return (base_it_type::operator*()).first; }
+    size_type index() const { return (base_it_type::operator*()).first; }
 
-    wsvector_const_iterator(void) {}
+    wsvector_const_iterator() {}
     wsvector_const_iterator(const wsvector_iterator<T> &it)
       : base_it_type(it) {}
     wsvector_const_iterator(const base_it_type &it) : base_it_type(it) {}
@@ -725,7 +725,7 @@ namespace gmm {
   */
   template<typename T> class wsvector : public std::map<size_type, T> {
   public:
-    
+
     typedef typename std::map<int, T>::size_type size_type;
     typedef std::map<size_type, T> base_type;
     typedef typename base_type::iterator iterator;
@@ -733,11 +733,11 @@ namespace gmm {
 
   protected:
     size_type nbl;
-    
+
   public:
     void clean(double eps);
     void resize(size_type);
-    
+
     inline ref_elt_vector<T, wsvector<T> > operator [](size_type c)
     { return ref_elt_vector<T, wsvector<T> >(this, c); }
 
@@ -750,9 +750,9 @@ namespace gmm {
     inline void wa(size_type c, const T &e) {
       GMM_ASSERT2(c < nbl, "out of range");
       if (e != T(0)) {
-       iterator it = this->lower_bound(c);
-       if (it != this->end() && it->first == c) it->second += e;
-       else base_type::operator [](c) = e;
+        iterator it = this->lower_bound(c);
+        if (it != this->end() && it->first == c) it->second += e;
+        else base_type::operator [](c) = e;
       }
     }
 
@@ -764,18 +764,18 @@ namespace gmm {
     }
 
     inline T operator [](size_type c) const { return r(c); }
-    
-    size_type nb_stored(void) const { return base_type::size(); }
-    size_type size(void) const { return nbl; }
+
+    size_type nb_stored() const { return base_type::size(); }
+    size_type size() const { return nbl; }
 
     void swap(wsvector<T> &v)
     { std::swap(nbl, v.nbl); std::map<size_type, T>::swap(v); }
-                                      
+
 
     /* Constructors */
     void init(size_type l) { nbl = l; this->clear(); }
     explicit wsvector(size_type l){ init(l); }
-    wsvector(void) { init(0); }
+    wsvector() { init(0); }
   };
 
   template<typename T>  void wsvector<T>::clean(double eps) {
@@ -815,10 +815,10 @@ namespace gmm {
     { o->clear(); }
     static void do_clear(this_type &v) { v.clear(); }
     static value_type access(const origin_type *o, const const_iterator &,
-                            const const_iterator &, size_type i)
+                             const const_iterator &, size_type i)
     { return (*o)[i]; }
     static reference access(origin_type *o, const iterator &, const iterator &,
-                           size_type i)
+                            size_type i)
     { return (*o)[i]; }
     static void resize(this_type &v, size_type n) { v.resize(n); }
   };
@@ -829,7 +829,7 @@ namespace gmm {
   /******* Optimized BLAS for wsvector<T> **********************************/
 
   template <typename T> inline void copy(const wsvector<T> &v1,
-                                        wsvector<T> &v2) {
+                                         wsvector<T> &v2) {
     GMM_ASSERT2(vect_size(v1) == vect_size(v2), "dimensions mismatch");
     v2 = v1;
   }
@@ -844,7 +844,7 @@ namespace gmm {
   }
   template <typename T> inline
   void copy(const simple_vector_ref<const wsvector<T> *> &v1,
-           wsvector<T> &v2)
+            wsvector<T> &v2)
   { copy(*(v1.origin), v2); }
   template <typename T> inline
   void copy(const simple_vector_ref<wsvector<T> *> &v1, wsvector<T> &v2)
@@ -853,9 +853,9 @@ namespace gmm {
   template <typename T> inline void clean(wsvector<T> &v, double eps) {
     typedef typename number_traits<T>::magnitude_type R;
     typename wsvector<T>::iterator it = v.begin(), ite = v.end(), itc;
-    while (it != ite) 
+    while (it != ite)
       if (gmm::abs((*it).second) <= R(eps))
-       { itc=it; ++it; v.erase(itc); } else ++it; 
+        { itc=it; ++it; v.erase(itc); } else ++it;
   }
 
   template <typename T>
@@ -881,7 +881,7 @@ namespace gmm {
     size_type c; T e;
     /* e is initialized by default to avoid some false warnings of valgrind.
        (from http://valgrind.org/docs/manual/mc-manual.html:
-      
+
        When memory is read into the CPU's floating point registers, the
        relevant V bits are read from memory and they are immediately
        checked. If any are invalid, an uninitialised value error is
@@ -890,7 +890,7 @@ namespace gmm {
        does not have to track the validity status of the floating-point
        registers.
     */
-    elt_rsvector_(void) : e(0) {}
+    elt_rsvector_() : e(0) {}
     elt_rsvector_(size_type cc) : c(cc), e(0) {}
     elt_rsvector_(size_type cc, const T &ee) : c(cc), e(ee) {}
     bool operator < (const elt_rsvector_ &a) const { return c < a.c; }
@@ -921,8 +921,8 @@ namespace gmm {
     bool operator ==(const iterator &i) const { return it == i.it; }
     bool operator !=(const iterator &i) const { return !(i == *this); }
 
-    size_type index(void) const { return it->c; }
-    rsvector_iterator(void) {}
+    size_type index() const { return it->c; }
+    rsvector_iterator() {}
     rsvector_iterator(const IT &i) : it(i) {}
   };
 
@@ -940,7 +940,7 @@ namespace gmm {
 
     reference operator *() const { return it->e; }
     pointer operator->() const { return &(operator*()); }
-    size_type index(void) const { return it->c; }
+    size_type index() const { return it->c; }
 
     iterator &operator ++() { ++it; return *this; }
     iterator operator ++(int) { iterator tmp = *this; ++(*this); return tmp; }
@@ -950,18 +950,18 @@ namespace gmm {
     bool operator ==(const iterator &i) const { return it == i.it; }
     bool operator !=(const iterator &i) const { return !(i == *this); }
 
-    rsvector_const_iterator(void) {}
+    rsvector_const_iterator() {}
     rsvector_const_iterator(const rsvector_iterator<T> &i) : it(i.it) {}
     rsvector_const_iterator(const IT &i) : it(i) {}
   };
 
   /**
      sparse vector built upon std::vector. Read access is fast,
-     but insertion is O(n) 
+     but insertion is O(n)
   */
   template<typename T> class rsvector : public std::vector<elt_rsvector_<T> > {
   public:
-    
+
     typedef std::vector<elt_rsvector_<T> > base_type_;
     typedef typename base_type_::iterator iterator;
     typedef typename base_type_::const_iterator const_iterator;
@@ -969,14 +969,14 @@ namespace gmm {
     typedef T value_type;
 
   protected:
-    size_type nbl;     /* size of the vector.                    */
-    
+    size_type nbl;    // size of the vector
+
   public:
 
     void sup(size_type j);
     void base_resize(size_type n) { base_type_::resize(n); }
     void resize(size_type);
-    
+
     ref_elt_vector<T, rsvector<T> > operator [](size_type c)
     { return ref_elt_vector<T, rsvector<T> >(this, c); }
 
@@ -986,16 +986,16 @@ namespace gmm {
     void swap_indices(size_type i, size_type j);
 
     inline T operator [](size_type c) const { return r(c); }
-    
-    size_type nb_stored(void) const { return base_type_::size(); }
-    size_type size(void) const { return nbl; }
-    void clear(void) { base_type_::resize(0); }
+
+    size_type nb_stored() const { return base_type_::size(); }
+    size_type size() const { return nbl; }
+    void clear() { base_type_::resize(0); }
     void swap(rsvector<T> &v)
     { std::swap(nbl, v.nbl); std::vector<elt_rsvector_<T> >::swap(v); }
 
     /* Constructeurs */
     explicit rsvector(size_type l) : nbl(l) { }
-    rsvector(void) : nbl(0) { }
+    rsvector() : nbl(0) { }
   };
 
   template <typename T>
@@ -1012,18 +1012,18 @@ namespace gmm {
 
       switch (situation) {
       case 1 : a = *iti; a.c = j; it = iti; ++it; ite = this->end();
-              for (; it != ite && it->c <= j; ++it, ++iti) *iti = *it;
-              *iti = a;
-              break;
+               for (; it != ite && it->c <= j; ++it, ++iti) *iti = *it;
+               *iti = a;
+               break;
       case 2 : a = *itj; a.c = i; it = itj; ite = this->begin();
-       if (it != ite) {
-         --it;
-         while (it->c >= i) { *itj = *it;  --itj; if (it==ite) break; --it; }
-       }
-       *itj = a;
-       break;
+        if (it != ite) {
+          --it;
+          while (it->c >= i) { *itj = *it;  --itj; if (it==ite) break; --it; }
+        }
+        *itj = a;
+        break;
       case 3 : std::swap(iti->e, itj->e);
-              break;
+               break;
       }
     }
   }
@@ -1033,8 +1033,8 @@ namespace gmm {
       elt_rsvector_<T> ev(j);
       iterator it = std::lower_bound(this->begin(), this->end(), ev);
       if (it != this->end() && it->c == j) {
-       for (iterator ite = this->end() - 1; it != ite; ++it) *it = *(it+1);
-       base_resize(nb_stored()-1);
+        for (iterator ite = this->end() - 1; it != ite; ++it) *it = *(it+1);
+        base_resize(nb_stored()-1);
       }
     }
   }
@@ -1042,7 +1042,7 @@ namespace gmm {
   template<typename T>  void rsvector<T>::resize(size_type n) {
     if (n < nbl) {
       for (size_type i = 0; i < nb_stored(); ++i)
-       if (base_type_::operator[](i).c >= n) { base_resize(i); break; }
+        if (base_type_::operator[](i).c >= n) { base_resize(i); break; }
     }
     nbl = n;
   }
@@ -1053,24 +1053,24 @@ namespace gmm {
     else {
       elt_rsvector_<T> ev(c, e);
       if (nb_stored() == 0) {
-       base_type_::push_back(ev);
+        base_type_::push_back(ev);
       }
       else {
-       iterator it = std::lower_bound(this->begin(), this->end(), ev);
-       if (it != this->end() && it->c == c) it->e = e;
-       else {
-         size_type ind = it - this->begin(), nb = this->nb_stored();
+        iterator it = std::lower_bound(this->begin(), this->end(), ev);
+        if (it != this->end() && it->c == c) it->e = e;
+        else {
+          size_type ind = it - this->begin(), nb = this->nb_stored();
           if (nb - ind > 1100)
             GMM_WARNING2("Inefficient addition of element in rsvector with "
                          << this->nb_stored() - ind << " non-zero entries");
-         base_type_::push_back(ev);
-         if (ind != nb) {
-           it = this->begin() + ind;
-           iterator ite = this->end(); --ite; iterator itee = ite; 
-           for (; ite != it; --ite) { --itee; *ite = *itee; }
-           *it = ev;
-         }
-       }
+          base_type_::push_back(ev);
+          if (ind != nb) {
+            it = this->begin() + ind;
+            iterator ite = this->end(); --ite; iterator itee = ite;
+            for (; ite != it; --ite) { --itee; *ite = *itee; }
+            *it = ev;
+          }
+        }
       }
     }
   }
@@ -1080,31 +1080,31 @@ namespace gmm {
     if (e != T(0)) {
       elt_rsvector_<T> ev(c, e);
       if (nb_stored() == 0) {
-       base_type_::push_back(ev);
+        base_type_::push_back(ev);
       }
       else {
-       iterator it = std::lower_bound(this->begin(), this->end(), ev);
-       if (it != this->end() && it->c == c) it->e += e;
-       else {
-         size_type ind = it - this->begin(), nb = this->nb_stored();
+        iterator it = std::lower_bound(this->begin(), this->end(), ev);
+        if (it != this->end() && it->c == c) it->e += e;
+        else {
+          size_type ind = it - this->begin(), nb = this->nb_stored();
           if (nb - ind > 1100)
             GMM_WARNING2("Inefficient addition of element in rsvector with "
                          << this->nb_stored() - ind << " non-zero entries");
-         base_type_::push_back(ev);
-         if (ind != nb) {
-           it = this->begin() + ind;
-           iterator ite = this->end(); --ite; iterator itee = ite; 
-           for (; ite != it; --ite) { --itee; *ite = *itee; }
-           *it = ev;
-         }
-       }
+          base_type_::push_back(ev);
+          if (ind != nb) {
+            it = this->begin() + ind;
+            iterator ite = this->end(); --ite; iterator itee = ite;
+            for (; ite != it; --ite) { --itee; *ite = *itee; }
+            *it = ev;
+          }
+        }
       }
     }
   }
-  
+
   template <typename T> T rsvector<T>::r(size_type c) const {
-    GMM_ASSERT2(c < nbl, "out of range. Index " << c 
-               << " for a length of " << nbl);
+    GMM_ASSERT2(c < nbl, "out of range. Index " << c
+                << " for a length of " << nbl);
     if (nb_stored() != 0) {
       elt_rsvector_<T> ev(c);
       const_iterator it = std::lower_bound(this->begin(), this->end(), ev);
@@ -1137,10 +1137,10 @@ namespace gmm {
     { o->clear(); }
     static void do_clear(this_type &v) { v.clear(); }
     static value_type access(const origin_type *o, const const_iterator &,
-                            const const_iterator &, size_type i)
+                             const const_iterator &, size_type i)
     { return (*o)[i]; }
     static reference access(origin_type *o, const iterator &, const iterator &,
-                           size_type i)
+                            size_type i)
     { return (*o)[i]; }
     static void resize(this_type &v, size_type n) { v.resize(n); }
   };
@@ -1151,7 +1151,7 @@ namespace gmm {
   /******* Optimized operations for rsvector<T> ****************************/
 
   template <typename T> inline void copy(const rsvector<T> &v1,
-                                        rsvector<T> &v2) {
+                                          rsvector<T> &v2) {
     GMM_ASSERT2(vect_size(v1) == vect_size(v2), "dimensions mismatch");
     v2 = v1;
   }
@@ -1166,39 +1166,39 @@ namespace gmm {
   }
   template <typename T> inline
   void copy(const simple_vector_ref<const rsvector<T> *> &v1,
-           rsvector<T> &v2)
+            rsvector<T> &v2)
   { copy(*(v1.origin), v2); }
   template <typename T> inline
   void copy(const simple_vector_ref<rsvector<T> *> &v1, rsvector<T> &v2)
   { copy(*(v1.origin), v2); }
 
   template <typename V, typename T> inline void add(const V &v1,
-                                                   rsvector<T> &v2) {
+                                                    rsvector<T> &v2) {
     if ((const void *)(&v1) != (const void *)(&v2)) {
       GMM_ASSERT2(vect_size(v1) == vect_size(v2), "dimensions mismatch");
-       add_rsvector(v1, v2, typename linalg_traits<V>::storage_type());
+        add_rsvector(v1, v2, typename linalg_traits<V>::storage_type());
     }
   }
 
-  template <typename V, typename T> 
+  template <typename V, typename T>
   inline void add_rsvector(const V &v1, rsvector<T> &v2, abstract_dense)
   { add(v1, v2, abstract_dense(), abstract_sparse()); }
 
-  template <typename V, typename T> 
+  template <typename V, typename T>
   inline void add_rsvector(const V &v1, rsvector<T> &v2, abstract_skyline)
   { add(v1, v2, abstract_skyline(), abstract_sparse()); }
 
-  template <typename V, typename T> 
+  template <typename V, typename T>
   void add_rsvector(const V &v1, rsvector<T> &v2, abstract_sparse) {
     add_rsvector(v1, v2, typename linalg_traits<V>::index_sorted());
   }
 
-  template <typename V, typename T> 
+  template <typename V, typename T>
   void add_rsvector(const V &v1, rsvector<T> &v2, linalg_false) {
     add(v1, v2, abstract_sparse(), abstract_sparse());
   }
 
-  template <typename V, typename T> 
+  template <typename V, typename T>
   void add_rsvector(const V &v1, rsvector<T> &v2, linalg_true) {
     typename linalg_traits<V>::const_iterator it1 = vect_const_begin(v1),
       ite1 = vect_const_end(v1);
@@ -1227,16 +1227,16 @@ namespace gmm {
     if ((const void *)(&v1) != (const void *)(&v2)) {
       GMM_ASSERT2(vect_size(v1) == vect_size(v2), "dimensions mismatch");
       if (same_origin(v1, v2))
-       GMM_WARNING2("a conflict is possible in vector copy\n");
+        GMM_WARNING2("a conflict is possible in vector copy\n");
       copy_rsvector(v1, v2, typename linalg_traits<V>::storage_type());
     }
   }
 
-  template <typename V, typename T> 
+  template <typename V, typename T>
   void copy_rsvector(const V &v1, rsvector<T> &v2, abstract_dense)
   { copy_vect(v1, v2, abstract_dense(), abstract_sparse()); }
 
-  template <typename V, typename T> 
+  template <typename V, typename T>
   void copy_rsvector(const V &v1, rsvector<T> &v2, abstract_skyline)
   { copy_vect(v1, v2, abstract_skyline(), abstract_sparse()); }
 
@@ -1244,7 +1244,7 @@ namespace gmm {
   void copy_rsvector(const V &v1, rsvector<T> &v2, abstract_sparse) {
     copy_rsvector(v1, v2, typename linalg_traits<V>::index_sorted());
   }
-  
+
   template <typename V, typename T2>
   void copy_rsvector(const V &v1, rsvector<T2> &v2, linalg_true) {
     typedef typename linalg_traits<V>::value_type T1;
@@ -1271,7 +1271,7 @@ namespace gmm {
     v2.base_resize(nn);
     std::sort(v2.begin(), v2.end());
   }
-  
+
   template <typename T> inline void clean(rsvector<T> &v, double eps) {
     typedef typename number_traits<T>::magnitude_type R;
     typename rsvector<T>::iterator it = v.begin(), ite = v.end();
@@ -1280,7 +1280,7 @@ namespace gmm {
       typename rsvector<T>::iterator itc = it;
       size_type erased = 1;
       for (++it; it != ite; ++it)
-       { *itc = *it; if (gmm::abs((*it).e) <= R(eps)) ++erased; else ++itc; }
+        { *itc = *it; if (gmm::abs((*it).e) <= R(eps)) ++erased; else ++itc; }
       v.base_resize(v.nb_stored() - erased);
     }
   }
@@ -1294,7 +1294,7 @@ namespace gmm {
     clean(*pv, eps);
     svr->begin_ = vect_begin(*pv); svr->end_ = vect_end(*pv);
   }
-  
+
   template <typename T>
   inline size_type nnz(const rsvector<T>& l) { return l.nb_stored(); }
 
@@ -1316,8 +1316,8 @@ namespace gmm {
 
     base_iterator it;
     size_type shift;
-    
-   
+
+
     iterator &operator ++()
     { ++it; ++shift; return *this; }
     iterator &operator --()
@@ -1336,21 +1336,21 @@ namespace gmm {
     { iterator tmp = *this; return (tmp -= i); }
     difference_type operator -(const iterator &i) const
     { return it - i.it; }
-       
+
     reference operator *() const
     { return *it; }
     reference operator [](int ii)
     { return *(it + ii); }
-    
+
     bool operator ==(const iterator &i) const
     { return it == i.it; }
     bool operator !=(const iterator &i) const
     { return !(i == *this); }
     bool operator < (const iterator &i) const
     { return it < i.it; }
-    size_type index(void) const { return shift; }
+    size_type index() const { return shift; }
 
-    slvector_iterator(void) {}
+    slvector_iterator() {}
     slvector_iterator(const base_iterator &iter, size_type s)
       : it(iter), shift(s) {}
   };
@@ -1367,8 +1367,8 @@ namespace gmm {
 
     base_iterator it;
     size_type shift;
-    
-   
+
+
     iterator &operator ++()
     { ++it; ++shift; return *this; }
     iterator &operator --()
@@ -1387,21 +1387,21 @@ namespace gmm {
     { iterator tmp = *this; return (tmp -= i); }
     difference_type operator -(const iterator &i) const
     { return it - i.it; }
-       
+
     value_type operator *() const
     { return *it; }
     value_type operator [](int ii)
     { return *(it + ii); }
-    
+
     bool operator ==(const iterator &i) const
     { return it == i.it; }
     bool operator !=(const iterator &i) const
     { return !(i == *this); }
     bool operator < (const iterator &i) const
     { return it < i.it; }
-    size_type index(void) const { return shift; }
+    size_type index() const { return shift; }
 
-    slvector_const_iterator(void) {}
+    slvector_const_iterator() {}
     slvector_const_iterator(const slvector_iterator<T>& iter)
       : it(iter.it), shift(iter.shift) {}
     slvector_const_iterator(const base_iterator &iter, size_type s)
@@ -1412,7 +1412,7 @@ namespace gmm {
   /** skyline vector.
    */
   template <typename T> class slvector {
-    
+
   public :
     typedef slvector_iterator<T> iterators;
     typedef slvector_const_iterator<T> const_iterators;
@@ -1427,17 +1427,17 @@ namespace gmm {
 
   public :
 
-    size_type size(void) const { return size_; }
-    size_type first(void) const { return shift; }
-    size_type last(void) const { return shift + data.size(); }
+    size_type size() const { return size_; }
+    size_type first() const { return shift; }
+    size_type last() const { return shift + data.size(); }
     ref_elt_vector<T, slvector<T> > operator [](size_type c)
     { return ref_elt_vector<T, slvector<T> >(this, c); }
 
-    typename std::vector<T>::iterator data_begin(void) { return data.begin(); }
-    typename std::vector<T>::iterator data_end(void) { return data.end(); }
-    typename std::vector<T>::const_iterator data_begin(void) const
+    typename std::vector<T>::iterator data_begin() { return data.begin(); }
+    typename std::vector<T>::iterator data_end() { return data.end(); }
+    typename std::vector<T>::const_iterator data_begin() const
       { return data.begin(); }
-    typename std::vector<T>::const_iterator data_end(void) const
+    typename std::vector<T>::const_iterator data_end() const
       { return data.end(); }
 
     void w(size_type c, const T &e);
@@ -1450,7 +1450,7 @@ namespace gmm {
 
     inline T operator [](size_type c) const { return r(c); }
     void resize(size_type);
-    void clear(void) { data.resize(0); shift = 0; }
+    void clear() { data.resize(0); shift = 0; }
     void swap(slvector<T> &v) {
       std::swap(data, v.data);
       std::swap(shift, v.shift);
@@ -1458,7 +1458,7 @@ namespace gmm {
     }
 
 
-    slvector(void) : data(0), shift(0), size_(0) {}
+    slvector() : data(0), shift(0), size_(0) {}
     explicit slvector(size_type l) : data(0), shift(0), size_(l) {}
     slvector(size_type l, size_type d, size_type s)
       : data(d), shift(s), size_(l) {}
@@ -1477,7 +1477,7 @@ namespace gmm {
     size_type s = data.size();
     if (!s) { data.resize(1); shift = c; }
     else if (c < shift) {
-      data.resize(s + shift - c); 
+      data.resize(s + shift - c);
       typename std::vector<T>::iterator it = data.begin(),it2=data.end()-1;
       typename std::vector<T>::iterator it3 = it2 - shift + c;
       for (; it3 >= it; --it3, --it2) *it2 = *it3;
@@ -1496,7 +1496,7 @@ namespace gmm {
     size_type s = data.size();
     if (!s) { data.resize(1, e); shift = c; return; }
     else if (c < shift) {
-      data.resize(s + shift - c); 
+      data.resize(s + shift - c);
       typename std::vector<T>::iterator it = data.begin(),it2=data.end()-1;
       typename std::vector<T>::iterator it3 = it2 - shift + c;
       for (; it3 >= it; --it3, --it2) *it2 = *it3;
@@ -1513,8 +1513,8 @@ namespace gmm {
     }
     data[c - shift] += e;
   }
-  
-  
+
+
   template <typename T> struct linalg_traits<slvector<T> > {
     typedef slvector<T> this_type;
     typedef this_type origin_type;
@@ -1541,10 +1541,10 @@ namespace gmm {
     { o->clear(); }
     static void do_clear(this_type &v) { v.clear(); }
     static value_type access(const origin_type *o, const const_iterator &,
-                            const const_iterator &, size_type i)
+                             const const_iterator &, size_type i)
     { return (*o)[i]; }
     static reference access(origin_type *o, const iterator &, const iterator &,
-                           size_type i)
+                            size_type i)
     { return (*o)[i]; }
     static void resize(this_type &v, size_type n) { v.resize(n); }
   };



reply via email to

[Prev in Thread] Current Thread [Next in Thread]