GetFEM  5.5
gmm_vector.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2002-2026 Yves Renard
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program. If not, see https://www.gnu.org/licenses/.
19 
20  As a special exception, you may use this file as it is a part of a free
21  software library without restriction. Specifically, if other files
22  instantiate templates or use macros or inline functions from this file,
23  or you compile this file and link it with other files to produce an
24  executable, this file does not by itself cause the resulting executable
25  to be covered by the GNU Lesser General Public License. This exception
26  does not however invalidate any other reasons why the executable file
27  might be covered by the GNU Lesser General Public License.
28 
29 ===========================================================================*/
30 /**@file gmm_vector.h
31  @author Yves Renard <[email protected]>
32  @date October 13, 2002.
33  @brief Declaration of the vector types (gmm::rsvector, gmm::wsvector,
34  gmm::slvector ,..)
35 */
36 #ifndef GMM_VECTOR_H__
37 #define GMM_VECTOR_H__
38 
39 #include <map>
40 #include "gmm_interface.h"
41 
42 namespace gmm {
43 
44  /*************************************************************************/
45  /* */
46  /* Class ref_elt_vector: reference on a vector component. */
47  /* */
48  /*************************************************************************/
49 
50 
51  template<typename T, typename V>
52  class ref_elt_vector {
53 
54  V *pm;
55  size_type l;
56 
57  public :
58 
59  operator T() const { return pm->r(l); }
60  ref_elt_vector(V *p, size_type ll) : pm(p), l(ll) {}
61  inline bool operator ==(T v) const { return ((*pm).r(l) == v); }
62  inline bool operator !=(T v) const { return ((*pm).r(l) != v); }
63  inline bool operator ==(std::complex<T> v) const
64  { return ((*pm).r(l) == v); }
65  inline bool operator !=(std::complex<T> v) const
66  { return ((*pm).r(l) != v); }
67  inline ref_elt_vector &operator +=(T v)
68  { (*pm).wa(l, v); return *this; }
69  inline ref_elt_vector &operator -=(T v)
70  { (*pm).wa(l, -v); return *this; }
71  inline ref_elt_vector &operator /=(T v)
72  { (*pm).w(l,(*pm).r(l) / v); return *this; }
73  inline ref_elt_vector &operator *=(T v)
74  { (*pm).w(l,(*pm).r(l) * v); return *this; }
75  inline ref_elt_vector &operator =(const ref_elt_vector &re)
76  { *this = T(re); return *this; }
77  inline ref_elt_vector &operator =(T v)
78  { (*pm).w(l,v); return *this; }
79  T operator +() { return T(*this); }
80  T operator -() { return -T(*this); }
81  T operator +(T v) { return T(*this)+ v; }
82  T operator -(T v) { return T(*this)- v; }
83  T operator *(T v) { return T(*this)* v; }
84  T operator /(T v) { return T(*this)/ v; }
85  std::complex<T> operator +(std::complex<T> v) { return T(*this)+ v; }
86  std::complex<T> operator -(std::complex<T> v) { return T(*this)- v; }
87  std::complex<T> operator *(std::complex<T> v) { return T(*this)* v; }
88  std::complex<T> operator /(std::complex<T> v) { return T(*this)/ v; }
89  };
90 
91 
92  template<typename T, typename V>
93  class ref_elt_vector<std::complex<T>,V> {
94 
95  V *pm;
96  size_type l;
97 
98  public :
99 
100  operator std::complex<T>() const { return pm->r(l); }
101  ref_elt_vector(V *p, size_type ll) : pm(p), l(ll) {}
102  ref_elt_vector(const ref_elt_vector &re) : pm(re.pm), l(re.l) {}
103  inline bool operator ==(std::complex<T> v) const
104  { return ((*pm).r(l) == v); }
105  inline bool operator !=(std::complex<T> v) const
106  { return ((*pm).r(l) != v); }
107  inline bool operator ==(T v) const { return ((*pm).r(l) == v); }
108  inline bool operator !=(T v) const { return ((*pm).r(l) != v); }
109  inline ref_elt_vector &operator +=(std::complex<T> v)
110  { (*pm).w(l,(*pm).r(l) + v); return *this; }
111  inline ref_elt_vector &operator -=(std::complex<T> v)
112  { (*pm).w(l,(*pm).r(l) - v); return *this; }
113  inline ref_elt_vector &operator /=(std::complex<T> v)
114  { (*pm).w(l,(*pm).r(l) / v); return *this; }
115  inline ref_elt_vector &operator *=(std::complex<T> v)
116  { (*pm).w(l,(*pm).r(l) * v); return *this; }
117  inline ref_elt_vector &operator =(const ref_elt_vector &re)
118  { *this = std::complex<T>(re); return *this; }
119  inline ref_elt_vector &operator =(std::complex<T> v)
120  { (*pm).w(l,v); return *this; }
121  inline ref_elt_vector &operator =(T v)
122  { (*pm).w(l,std::complex<T>(v)); return *this; }
123  inline ref_elt_vector &operator +=(T v)
124  { (*pm).w(l,(*pm).r(l) + v); return *this; }
125  inline ref_elt_vector &operator -=(T v)
126  { (*pm).w(l,(*pm).r(l) - v); return *this; }
127  inline ref_elt_vector &operator /=(T v)
128  { (*pm).w(l,(*pm).r(l) / v); return *this; }
129  inline ref_elt_vector &operator *=(T v)
130  { (*pm).w(l,(*pm).r(l) * v); return *this; }
131  std::complex<T> operator +() { return std::complex<T>(*this); }
132  std::complex<T> operator -() { return -std::complex<T>(*this); }
133  std::complex<T> operator +(T v) { return std::complex<T>(*this)+ v; }
134  std::complex<T> operator -(T v) { return std::complex<T>(*this)- v; }
135  std::complex<T> operator *(T v) { return std::complex<T>(*this)* v; }
136  std::complex<T> operator /(T v) { return std::complex<T>(*this)/ v; }
137  std::complex<T> operator +(std::complex<T> v)
138  { return std::complex<T>(*this)+ v; }
139  std::complex<T> operator -(std::complex<T> v)
140  { return std::complex<T>(*this)- v; }
141  std::complex<T> operator *(std::complex<T> v)
142  { return std::complex<T>(*this)* v; }
143  std::complex<T> operator /(std::complex<T> v)
144  { return std::complex<T>(*this)/ v; }
145  };
146 
147 
148  template<typename T, typename V> inline
149  bool operator ==(T v, const ref_elt_vector<T, V> &re) { return (v==T(re)); }
150  template<typename T, typename V> inline
151  bool operator !=(T v, const ref_elt_vector<T, V> &re) { return (v!=T(re)); }
152  template<typename T, typename V> inline
153  T &operator +=(T &v, const ref_elt_vector<T, V> &re)
154  { v += T(re); return v; }
155  template<typename T, typename V> inline
156  T &operator -=(T &v, const ref_elt_vector<T, V> &re)
157  { v -= T(re); return v; }
158  template<typename T, typename V> inline
159  T &operator *=(T &v, const ref_elt_vector<T, V> &re)
160  { v *= T(re); return v; }
161  template<typename T, typename V> inline
162  T &operator /=(T &v, const ref_elt_vector<T, V> &re)
163  { v /= T(re); return v; }
164  template<typename T, typename V> inline
165  T operator +(T v, const ref_elt_vector<T, V> &re) { return v+ T(re); }
166  template<typename T, typename V> inline
167  T operator -(T v, const ref_elt_vector<T, V> &re) { return v- T(re); }
168  template<typename T, typename V> inline
169  T operator *(T v, const ref_elt_vector<T, V> &re) { return v* T(re); }
170  template<typename T, typename V> inline
171  T operator /(T v, const ref_elt_vector<T, V> &re) { return v/ T(re); }
172  template<typename T, typename V> inline
173  std::complex<T> operator +(std::complex<T> v, const ref_elt_vector<T, V> &re)
174  { return v+ T(re); }
175  template<typename T, typename V> inline
176  std::complex<T> operator -(std::complex<T> v, const ref_elt_vector<T, V> &re)
177  { return v- T(re); }
178  template<typename T, typename V> inline
179  std::complex<T> operator *(std::complex<T> v, const ref_elt_vector<T, V> &re)
180  { return v* T(re); }
181  template<typename T, typename V> inline
182  std::complex<T> operator /(std::complex<T> v, const ref_elt_vector<T, V> &re)
183  { return v/ T(re); }
184  template<typename T, typename V> inline
185  std::complex<T> operator +(T v, const ref_elt_vector<std::complex<T>, V> &re)
186  { return v+ std::complex<T>(re); }
187  template<typename T, typename V> inline
188  std::complex<T> operator -(T v, const ref_elt_vector<std::complex<T>, V> &re)
189  { return v- std::complex<T>(re); }
190  template<typename T, typename V> inline
191  std::complex<T> operator *(T v, const ref_elt_vector<std::complex<T>, V> &re)
192  { return v* std::complex<T>(re); }
193  template<typename T, typename V> inline
194  std::complex<T> operator /(T v, const ref_elt_vector<std::complex<T>, V> &re)
195  { return v/ std::complex<T>(re); }
196  template<typename T, typename V> inline
197  typename number_traits<T>::magnitude_type
198  abs(const ref_elt_vector<T, V> &re) { return gmm::abs(T(re)); }
199  template<typename T, typename V> inline
200  T sqr(const ref_elt_vector<T, V> &re) { return gmm::sqr(T(re)); }
201  template<typename T, typename V> inline
202  typename number_traits<T>::magnitude_type
203  abs_sqr(const ref_elt_vector<T, V> &re) { return gmm::abs_sqr(T(re)); }
204  template<typename T, typename V> inline
205  T conj(const ref_elt_vector<T, V> &re) { return gmm::conj(T(re)); }
206  template<typename T, typename V> std::ostream &operator <<
207  (std::ostream &o, const ref_elt_vector<T, V> &re) { o << T(re); return o; }
208  template<typename T, typename V> inline
209  typename number_traits<T>::magnitude_type
210  real(const ref_elt_vector<T, V> &re) { return gmm::real(T(re)); }
211  template<typename T, typename V> inline
212  typename number_traits<T>::magnitude_type
213  imag(const ref_elt_vector<T, V> &re) { return gmm::imag(T(re)); }
214 
215 
216  /*************************************************************************/
217  /* */
218  /* Class dsvector: sparse vector optimized for random write operations */
219  /* with constant complexity for read and write operations. */
220  /* Based on distribution sort principle. */
221  /* Cheap for densely populated vectors. */
222  /* */
223  /*************************************************************************/
224 
225  template<typename T> class dsvector;
226 
227  template<typename T>
228  struct dsvector_iterator {
229  size_type i; // Current index.
230  T* p; // Pointer to the current position.
231  dsvector<T> *v; // Pointer to the vector.
232 
233  typedef T value_type;
234  typedef value_type* pointer;
235  typedef const value_type* const_pointer;
236  typedef value_type& reference;
237  // typedef size_t size_type;
238  typedef ptrdiff_t difference_type;
239  typedef std::bidirectional_iterator_tag iterator_category;
240  typedef dsvector_iterator<T> iterator;
241 
242  reference operator *() const { return *p; }
243  pointer operator->() const { return &(operator*()); }
244 
245  iterator &operator ++() {
246  for (size_type k = (i & 15); k < 15; ++k)
247  { ++p; ++i; if (*p != T(0)) return *this; }
248  v->next_pos(*(const_cast<const_pointer *>(&(p))), i);
249  return *this;
250  }
251  iterator operator ++(int) { iterator tmp = *this; ++(*this); return tmp; }
252  iterator &operator --() {
253  for (size_type k = (i & 15); k > 0; --k)
254  { --p; --i; if (*p != T(0)) return *this; }
255  v->previous_pos(p, i);
256  return *this;
257  }
258  iterator operator --(int) { iterator tmp = *this; --(*this); return tmp; }
259 
260  bool operator ==(const iterator &it) const
261  { return (i == it.i && p == it.p && v == it.v); }
262  bool operator !=(const iterator &it) const
263  { return !(it == *this); }
264 
265  size_type index() const { return i; }
266 
267  dsvector_iterator() : i(size_type(-1)), p(0), v(0) {}
268  dsvector_iterator(dsvector<T> &w) : i(size_type(-1)), p(0), v(&w) {};
269  };
270 
271 
272  template<typename T>
273  struct dsvector_const_iterator {
274  size_type i; // Current index.
275  const T* p; // Pointer to the current position.
276  const dsvector<T> *v; // Pointer to the vector.
277 
278  typedef T value_type;
279  typedef const value_type* pointer;
280  typedef const value_type& reference;
281  // typedef size_t size_type;
282  typedef ptrdiff_t difference_type;
283  typedef std::bidirectional_iterator_tag iterator_category;
284  typedef dsvector_const_iterator<T> iterator;
285 
286  reference operator *() const { return *p; }
287  pointer operator->() const { return &(operator*()); }
288  iterator &operator ++() {
289  for (size_type k = (i & 15); k < 15; ++k)
290  { ++p; ++i; if (*p != T(0)) return *this; }
291  v->next_pos(p, i);
292  return *this;
293  }
294  iterator operator ++(int) { iterator tmp = *this; ++(*this); return tmp; }
295  iterator &operator --() {
296  for (size_type k = (i & 15); k > 0; --k)
297  { --p; --i; if (*p != T(0)) return *this; }
298  v->previous_pos(p, i);
299  return *this;
300  }
301  iterator operator --(int) { iterator tmp = *this; --(*this); return tmp; }
302 
303  bool operator ==(const iterator &it) const
304  { return (i == it.i && p == it.p && v == it.v); }
305  bool operator !=(const iterator &it) const
306  { return !(it == *this); }
307 
308  size_type index() const { return i; }
309 
310  dsvector_const_iterator() : i(size_type(-1)), p(0) {}
311  dsvector_const_iterator(const dsvector_iterator<T> &it)
312  : i(it.i), p(it.p), v(it.v) {}
313  dsvector_const_iterator(const dsvector<T> &w)
314  : i(size_type(-1)), p(0), v(&w) {};
315  };
316 
317 
318  /**
319  Sparse vector built on distribution sort principle.
320  Read and write access have a constant complexity depending only on the
321  vector size.
322  */
323  template<typename T>
324  class dsvector {
325 
326  typedef dsvector_iterator<T> iterator;
327  typedef dsvector_const_iterator<T> const_iterator;
328  typedef dsvector<T> this_type;
329  typedef T * pointer;
330  typedef const T * const_pointer;
331  typedef void * void_pointer;
332  typedef const void * const_void_pointer;
333 
334  protected:
335  size_type n; // Potential vector size
336  size_type depth; // Number of row of pointer arrays
337  size_type mask; // Mask for the first pointer array
338  size_type shift; // Shift for the first pointer array
339  void_pointer root_ptr; // Root pointer
340 
341  const T *read_access(size_type i) const {
342  GMM_ASSERT1(i < n, "index out of range");
343  size_type my_mask = mask, my_shift = shift;
344  void_pointer p = root_ptr;
345  if (!p) return 0;
346  for (size_type k = 0; k < depth; ++k) {
347  p = ((void **)(p))[(i & my_mask) >> my_shift];
348  if (!p) return 0;
349  my_mask = (my_mask >> 4);
350  my_shift -= 4;
351  }
352  GMM_ASSERT1(my_shift == 0, "internal error");
353  GMM_ASSERT1(my_mask == 15, "internal error");
354  return &(((const T *)(p))[i & 15]);
355  }
356 
357  T *write_access(size_type i) {
358  GMM_ASSERT1(i < n, "index " << i << " out of range (size " << n << ")");
359  size_type my_mask = mask, my_shift = shift;
360  if (!root_ptr) {
361  if (depth) {
362  root_ptr = new void_pointer[16];
363  std::memset(root_ptr, 0, 16*sizeof(void_pointer));
364  } else {
365  root_ptr = new T[16];
366  for (size_type l = 0; l < 16; ++l) ((T *)(root_ptr))[l] = T(0);
367  }
368  }
369 
370  void_pointer p = root_ptr;
371  for (size_type k = 0; k < depth; ++k) {
372  size_type j = (i & my_mask) >> my_shift;
373  void_pointer q = ((void_pointer *)(p))[j];
374  if (!q) {
375  if (k+1 != depth) {
376  q = new void_pointer[16];
377  std::memset(q, 0, 16*sizeof(void_pointer));
378  } else {
379  q = new T[16];
380  for (size_type l = 0; l < 16; ++l) ((T *)(q))[l] = T(0);
381  }
382  ((void_pointer *)(p))[j] = q;
383  }
384  p = q;
385  my_mask = (my_mask >> 4);
386  my_shift -= 4;
387  }
388  GMM_ASSERT1(my_shift == 0, "internal error");
389  GMM_ASSERT1(my_mask == 15, "internal error " << my_mask);
390  return &(((T *)(p))[i & 15]);
391  }
392 
393  void init(size_type n_) {
394  n = n_; depth = 0; shift = 0; mask = 1; if (n_) --n_;
395  while (n_) { n_ /= 16; ++depth; shift += 4; mask *= 16; }
396  mask--; if (shift) shift -= 4; if (depth) --depth;
397  root_ptr = 0;
398  }
399 
400  void rec_del(void_pointer p, size_type my_depth) {
401  if (my_depth) {
402  for (size_type k = 0; k < 16; ++k)
403  if (((void_pointer *)(p))[k])
404  rec_del(((void_pointer *)(p))[k], my_depth-1);
405  delete[] ((void_pointer *)(p));
406  } else {
407  delete[] ((T *)(p));
408  }
409  }
410 
411  void rec_clean(void_pointer p, size_type my_depth, double eps) {
412  if (my_depth) {
413  for (size_type k = 0; k < 16; ++k)
414  if (((void_pointer *)(p))[k])
415  rec_clean(((void_pointer *)(p))[k], my_depth-1, eps);
416  } else {
417  for (size_type k = 0; k < 16; ++k)
418  if (gmm::abs(((T *)(p))[k]) <= eps) ((T *)(p))[k] = T(0);
419  }
420  }
421 
422  void rec_clean_i(void_pointer p, size_type my_depth, size_type my_mask,
423  size_type i, size_type base) {
424  if (my_depth) {
425  my_mask = (my_mask >> 4);
426  for (size_type k = 0; k < 16; ++k)
427  if (((void_pointer *)(p))[k] && (base + (k+1)*(mask+1)) >= i)
428  rec_clean_i(((void_pointer *)(p))[k], my_depth-1, my_mask,
429  i, base + k*(my_mask+1));
430  } else {
431  for (size_type k = 0; k < 16; ++k)
432  if (base+k > i) ((T *)(p))[k] = T(0);
433  }
434  }
435 
436 
437  size_type rec_nnz(void_pointer p, size_type my_depth) const {
438  size_type nn = 0;
439  if (my_depth) {
440  for (size_type k = 0; k < 16; ++k)
441  if (((void_pointer *)(p))[k])
442  nn += rec_nnz(((void_pointer *)(p))[k], my_depth-1);
443  } else {
444  for (size_type k = 0; k < 16; ++k)
445  if (((const T *)(p))[k] != T(0)) nn++;
446  }
447  return nn;
448  }
449 
450  void copy_rec(void_pointer &p, const_void_pointer q, size_type my_depth) {
451  if (my_depth) {
452  p = new void_pointer[16];
453  std::memset(p, 0, 16*sizeof(void_pointer));
454  for (size_type l = 0; l < 16; ++l)
455  if (((const const_void_pointer *)(q))[l])
456  copy_rec(((void_pointer *)(p))[l],
457  ((const const_void_pointer *)(q))[l], my_depth-1);
458  } else {
459  p = new T[16];
460  for (size_type l = 0; l < 16; ++l) ((T *)(p))[l] = ((const T *)(q))[l];
461  }
462  }
463 
464  void copy(const dsvector<T> &v) {
465  if (root_ptr) rec_del(root_ptr, depth);
466  root_ptr = 0;
467  mask = v.mask; depth = v.depth; n = v.n; shift = v.shift;
468  if (v.root_ptr) copy_rec(root_ptr, v.root_ptr, depth);
469  }
470 
471  void next_pos_rec(void_pointer p, size_type my_depth, size_type my_mask,
472  const_pointer &pp, size_type &i, size_type base) const {
473  size_type ii = i;
474  if (my_depth) {
475  my_mask = (my_mask >> 4);
476  for (size_type k = 0; k < 16; ++k)
477  if (((void_pointer *)(p))[k] && (base + (k+1)*(my_mask+1)) >= i) {
478  next_pos_rec(((void_pointer *)(p))[k], my_depth-1, my_mask,
479  pp, i, base + k*(my_mask+1));
480  if (i != size_type(-1)) return; else i = ii;
481  }
482  i = size_type(-1); pp = 0;
483  } else {
484  for (size_type k = 0; k < 16; ++k)
485  if (base+k > i && ((const_pointer)(p))[k] != T(0))
486  { i = base+k; pp = &(((const_pointer)(p))[k]); return; }
487  i = size_type(-1); pp = 0;
488  }
489  }
490 
491  void previous_pos_rec(void_pointer p, size_type my_depth, size_type my_mask,
492  const_pointer &pp, size_type &i,
493  size_type base) const {
494  size_type ii = i;
495  if (my_depth) {
496  my_mask = (my_mask >> 4);
497  for (size_type k = 15; k != size_type(-1); --k)
498  if (((void_pointer *)(p))[k] && ((base + k*(my_mask+1)) < i)) {
499  previous_pos_rec(((void_pointer *)(p))[k], my_depth-1,
500  my_mask, pp, i, base + k*(my_mask+1));
501  if (i != size_type(-1)) return; else i = ii;
502  }
503  i = size_type(-1); pp = 0;
504  } else {
505  for (size_type k = 15; k != size_type(-1); --k)
506  if (base+k < i && ((const_pointer)(p))[k] != T(0))
507  { i = base+k; pp = &(((const_pointer)(p))[k]); return; }
508  i = size_type(-1); pp = 0;
509  }
510  }
511 
512 
513  public:
514  void clean(double eps) { if (root_ptr) rec_clean(root_ptr, depth); }
515  void resize(size_type n_) {
516  if (n_ != n) {
517  n = n_;
518  if (n_ < n) { // Depth unchanged (a choice)
519  if (root_ptr) rec_clean_i(root_ptr, depth, mask, n_, 0);
520  } else {
521  // may change the depth (add some levels)
522  size_type my_depth = 0, my_shift = 0, my_mask = 1; if (n_) --n_;
523  while (n_) { n_ /= 16; ++my_depth; my_shift += 4; my_mask *= 16; }
524  my_mask--; if (my_shift) my_shift -= 4; if (my_depth) --my_depth;
525  if (my_depth > depth || depth == 0) {
526  if (root_ptr) {
527  for (size_type k = depth; k < my_depth; ++k) {
528  void_pointer *q = new void_pointer [16];
529  std::memset(q, 0, 16*sizeof(void_pointer));
530  q[0] = root_ptr; root_ptr = q;
531  }
532  }
533  mask = my_mask; depth = my_depth; shift = my_shift;
534  }
535  }
536  }
537  }
538 
539  void clear() {
540  if (root_ptr)
541  rec_del(root_ptr, depth);
542  root_ptr = 0;
543  }
544 
545  void next_pos(const_pointer &pp, size_type &i) const {
546  if (!root_ptr || i >= n) { pp = 0, i = size_type(-1); return; }
547  next_pos_rec(root_ptr, depth, mask, pp, i, 0);
548  }
549 
550  void previous_pos(const_pointer &pp, size_type &i) const {
551  if (!root_ptr) { pp = 0, i = size_type(-1); return; }
552  if (i == size_type(-1)) { i = n; }
553  previous_pos_rec(root_ptr, depth, mask, pp, i, 0);
554  }
555 
556  iterator begin() {
557  iterator it(*this);
558  if (n && root_ptr) {
559  it.i = 0; it.p = const_cast<T *>(read_access(0));
560  if (!(it.p) || *(it.p) == T(0))
561  next_pos(*(const_cast<const_pointer *>(&(it.p))), it.i);
562  }
563  return it;
564  }
565 
566  iterator end() { return iterator(*this); }
567 
568  const_iterator begin() const {
569  const_iterator it(*this);
570  if (n && root_ptr) {
571  it.i = 0; it.p = read_access(0);
572  if (!(it.p) || *(it.p) == T(0)) next_pos(it.p, it.i);
573  }
574  return it;
575  }
576 
577  const_iterator end() const { return const_iterator(*this); }
578 
579  inline ref_elt_vector<T, dsvector<T> > operator [](size_type c)
580  { return ref_elt_vector<T, dsvector<T> >(this, c); }
581 
582  inline void w(size_type c, const T &e) {
583  if (e == T(0)) { if (read_access(c)) *(write_access(c)) = e; }
584  else *(write_access(c)) = e;
585  }
586 
587  inline void wa(size_type c, const T &e)
588  { if (e != T(0)) { *(write_access(c)) += e; } }
589 
590  inline T r(size_type c) const
591  { const T *p = read_access(c); if (p) return *p; else return T(0); }
592 
593  inline T operator [](size_type c) const { return r(c); }
594 
595  size_type nnz() const
596  { if (root_ptr) return rec_nnz(root_ptr, depth); else return 0; }
597  size_type size() const { return n; }
598 
599  void swap(dsvector<T> &v) {
600  std::swap(n, v.n); std::swap(root_ptr, v.root_ptr);
601  std::swap(depth, v.depth); std::swap(shift, v.shift);
602  std::swap(mask, v.mask);
603  }
604 
605  /* Constructors */
606  dsvector(const dsvector<T> &v) { init(0); copy(v); }
607  dsvector<T> &operator =(const dsvector<T> &v) { copy(v); return *this; }
608  explicit dsvector(size_type l){ init(l); }
609  dsvector() { init(0); }
610  ~dsvector() { if (root_ptr) rec_del(root_ptr, depth); root_ptr = 0; }
611  };
612 
613 
614  template <typename T>
615  struct linalg_traits<dsvector<T>> {
616  typedef dsvector<T> this_type;
617  typedef this_type origin_type;
618  typedef linalg_false is_reference;
619  typedef abstract_vector linalg_type;
620  typedef T value_type;
621  typedef ref_elt_vector<T, dsvector<T> > reference;
622  typedef dsvector_iterator<T> iterator;
623  typedef dsvector_const_iterator<T> const_iterator;
624  typedef abstract_sparse storage_type;
625  typedef linalg_true index_sorted;
626  static size_type size(const this_type &v) { return v.size(); }
627  static iterator begin(this_type &v) { return v.begin(); }
628  static const_iterator begin(const this_type &v) { return v.begin(); }
629  static iterator end(this_type &v) { return v.end(); }
630  static const_iterator end(const this_type &v) { return v.end(); }
631  static origin_type* origin(this_type &v) { return &v; }
632  static const origin_type* origin(const this_type &v) { return &v; }
633  static void clear(origin_type* o, const iterator &, const iterator &)
634  { o->clear(); }
635  static void do_clear(this_type &v) { v.clear(); }
636  static value_type access(const origin_type *o, const const_iterator &,
637  const const_iterator &, size_type i)
638  { return (*o)[i]; }
639  static reference access(origin_type *o, const iterator &, const iterator &,
640  size_type i)
641  { return (*o)[i]; }
642  static void resize(this_type &v, size_type n) { v.resize(n); }
643  };
644 
645  template<typename T>
646  std::ostream &operator <<(std::ostream &o, const dsvector<T>& v)
647  { gmm::write(o,v); return o; }
648 
649  /******* Optimized operations for dsvector<T> ****************************/
650 
651  template <typename T> inline
652  void copy(const dsvector<T> &v1, dsvector<T> &v2) {
653  GMM_ASSERT2(v1.size() == v2.size(), "dimensions mismatch");
654  v2 = v1;
655  }
656 
657  template <typename T> inline
658  void copy(const dsvector<T> &v1, const dsvector<T> &v2) {
659  GMM_ASSERT2(v1.size() == v2.size(), "dimensions mismatch");
660  v2 = const_cast<dsvector<T> &>(v1);
661  }
662 
663  template <typename T> inline
664  void copy(const dsvector<T> &v1,
665  const simple_vector_ref<dsvector<T> *> &v2) {
666  simple_vector_ref<dsvector<T> *>
667  *svr = const_cast<simple_vector_ref<dsvector<T> *> *>(&v2);
668  dsvector<T>
669  *pv = const_cast<dsvector<T> *>((v2.origin));
670  GMM_ASSERT2(vect_size(v1) == vect_size(v2), "dimensions mismatch");
671  *pv = v1; svr->begin_ = vect_begin(*pv); svr->end_ = vect_end(*pv);
672  }
673 
674  template <typename T> inline
675  void copy(const simple_vector_ref<const dsvector<T> *> &v1,
676  dsvector<T> &v2)
677  { copy(*(v1.origin), v2); }
678 
679  template <typename T> inline
680  void copy(const simple_vector_ref<dsvector<T> *> &v1, dsvector<T> &v2)
681  { copy(*(v1.origin), v2); }
682 
683  template <typename T> inline
684  void copy(const simple_vector_ref<dsvector<T> *> &v1,
685  const simple_vector_ref<dsvector<T> *> &v2)
686  { copy(*(v1.origin), v2); }
687 
688  template <typename T> inline
689  void copy(const simple_vector_ref<const dsvector<T> *> &v1,
690  const simple_vector_ref<dsvector<T> *> &v2)
691  { copy(*(v1.origin), v2); }
692 
693  template <typename T> inline
694  size_type nnz(const dsvector<T>& l) { return l.nnz(); }
695 
696 
697  /*************************************************************************/
698  /* */
699  /* Class wsvector: sparse vector optimized for random write operations, */
700  /* with log(n) complexity for read and write operations. */
701  /* Based on std::map */
702  /* */
703  /*************************************************************************/
704 
705  template<typename T>
706  struct wsvector_iterator : public std::map<size_type, T>::iterator {
707 
708  typedef typename std::map<size_type, T>::iterator base_it_type;
709  typedef T value_type;
710  typedef value_type* pointer;
711  typedef value_type& reference;
712  // typedef size_t size_type;
713  typedef ptrdiff_t difference_type;
714  typedef std::bidirectional_iterator_tag iterator_category;
715 
716  reference operator *() const { return (base_it_type::operator*()).second; }
717  pointer operator->() const { return &(operator*()); }
718  size_type index() const { return (base_it_type::operator*()).first; }
719 
720  wsvector_iterator() {}
721  wsvector_iterator(const base_it_type &it) : base_it_type(it) {}
722  };
723 
724  template<typename T>
725  struct wsvector_const_iterator
726  : public std::map<size_type, T>::const_iterator {
727 
728  typedef typename std::map<size_type, T>::const_iterator base_it_type;
729  typedef T value_type;
730  typedef const value_type* pointer;
731  typedef const value_type& reference;
732  // typedef size_t size_type;
733  typedef ptrdiff_t difference_type;
734  typedef std::bidirectional_iterator_tag iterator_category;
735 
736  reference operator *() const { return (base_it_type::operator*()).second; }
737  pointer operator->() const { return &(operator*()); }
738  size_type index() const { return (base_it_type::operator*()).first; }
739 
740  wsvector_const_iterator() {}
741  wsvector_const_iterator(const wsvector_iterator<T> &it)
742  : base_it_type(typename std::map<size_type, T>::iterator(it)) {}
743  wsvector_const_iterator(const base_it_type &it) : base_it_type(it) {}
744  };
745 
746 
747  /**
748  sparse vector built upon std::map.
749  Read and write access are quite fast (log n)
750  */
751  template<typename T>
752  class wsvector : public std::map<size_type, T> {
753 
754  public:
755 
756  typedef typename std::map<int, T>::size_type size_type;
757  typedef std::map<size_type, T> base_type;
758  typedef typename base_type::iterator iterator;
759  typedef typename base_type::const_iterator const_iterator;
760 
761  protected:
762  size_type nbl;
763 
764  public:
765  void clean(double eps);
766  void resize(size_type);
767 
768  inline ref_elt_vector<T, wsvector<T> > operator [](size_type c)
769  { return ref_elt_vector<T, wsvector<T> >(this, c); }
770 
771  inline void w(size_type c, const T &e) {
772  GMM_ASSERT2(c < nbl, "out of range");
773  if (e == T(0)) { this->erase(c); }
774  else base_type::operator [](c) = e;
775  }
776 
777  inline void wa(size_type c, const T &e) {
778  GMM_ASSERT2(c < nbl, "out of range");
779  if (e != T(0)) {
780  iterator it = this->lower_bound(c);
781  if (it != this->end() && it->first == c) it->second += e;
782  else base_type::operator [](c) = e;
783  }
784  }
785 
786  inline T r(size_type c) const {
787  GMM_ASSERT2(c < nbl, "out of range");
788  const_iterator it = this->lower_bound(c);
789  if (it != this->end() && c == it->first) return it->second;
790  else return T(0);
791  }
792 
793  inline T operator [](size_type c) const { return r(c); }
794 
795  size_type nb_stored() const { return base_type::size(); }
796  size_type size() const { return nbl; }
797 
798  void swap(wsvector<T> &v)
799  { std::swap(nbl, v.nbl); std::map<size_type, T>::swap(v); }
800 
801 
802  /* Constructors */
803  void init(size_type l) { nbl = l; this->clear(); }
804  explicit wsvector(size_type l){ init(l); }
805  wsvector() { init(0); }
806  };
807 
808  template<typename T>
809  void wsvector<T>::clean(double eps) {
810  iterator it = this->begin(), itf = it, ite = this->end();
811  while (it != ite) {
812  ++itf; if (gmm::abs(it->second) <= eps) this->erase(it); it = itf;
813  }
814  }
815 
816  template<typename T>
817  void wsvector<T>::resize(size_type n) {
818  if (n < nbl) {
819  iterator it = this->begin(), itf = it, ite = this->end();
820  while (it != ite) { ++itf; if (it->first >= n) this->erase(it); it=itf; }
821  }
822  nbl = n;
823  }
824 
825  template <typename T>
826  struct linalg_traits<wsvector<T> > {
827 
828  typedef wsvector<T> this_type;
829  typedef this_type origin_type;
830  typedef linalg_false is_reference;
831  typedef abstract_vector linalg_type;
832  typedef T value_type;
833  typedef ref_elt_vector<T, wsvector<T> > reference;
834  typedef wsvector_iterator<T> iterator;
835  typedef wsvector_const_iterator<T> const_iterator;
836  typedef abstract_sparse storage_type;
837  typedef linalg_true index_sorted;
838  static size_type size(const this_type &v) { return v.size(); }
839  static iterator begin(this_type &v) { return v.begin(); }
840  static const_iterator begin(const this_type &v) { return v.begin(); }
841  static iterator end(this_type &v) { return v.end(); }
842  static const_iterator end(const this_type &v) { return v.end(); }
843  static origin_type* origin(this_type &v) { return &v; }
844  static const origin_type* origin(const this_type &v) { return &v; }
845  static void clear(origin_type* o, const iterator &, const iterator &)
846  { o->clear(); }
847  static void do_clear(this_type &v) { v.clear(); }
848  static value_type access(const origin_type *o, const const_iterator &,
849  const const_iterator &, size_type i)
850  { return (*o)[i]; }
851  static reference access(origin_type *o, const iterator &, const iterator &,
852  size_type i)
853  { return (*o)[i]; }
854  static void resize(this_type &v, size_type n) { v.resize(n); }
855  };
856 
857  template<typename T>
858  std::ostream &operator <<(std::ostream &o, const wsvector<T>& v)
859  { gmm::write(o,v); return o; }
860 
861 
862  /******* Optimized BLAS for wsvector<T> **********************************/
863 
864  template <typename T> inline void copy(const wsvector<T> &v1,
865  wsvector<T> &v2) {
866  GMM_ASSERT2(vect_size(v1) == vect_size(v2), "dimensions mismatch");
867  v2 = v1;
868  }
869  template <typename T> inline
870  void copy(const wsvector<T> &v1, const simple_vector_ref<wsvector<T> *> &v2){
871  simple_vector_ref<wsvector<T> *>
872  *svr = const_cast<simple_vector_ref<wsvector<T> *> *>(&v2);
873  wsvector<T>
874  *pv = const_cast<wsvector<T> *>(v2.origin);
875  GMM_ASSERT2(vect_size(v1) == vect_size(v2), "dimensions mismatch");
876  *pv = v1; svr->begin_ = vect_begin(*pv); svr->end_ = vect_end(*pv);
877  }
878  template <typename T> inline
879  void copy(const simple_vector_ref<const wsvector<T> *> &v1,
880  wsvector<T> &v2)
881  { copy(*(v1.origin), v2); }
882  template <typename T> inline
883  void copy(const simple_vector_ref<wsvector<T> *> &v1, wsvector<T> &v2)
884  { copy(*(v1.origin), v2); }
885 
886  template <typename T> inline void clean(wsvector<T> &v, double eps) {
887  typedef typename number_traits<T>::magnitude_type R;
888  typename wsvector<T>::iterator it = v.begin(), ite = v.end(), itc;
889  while (it != ite)
890  if (gmm::abs((*it).second) <= R(eps))
891  { itc=it; ++it; v.erase(itc); } else ++it;
892  }
893 
894  template <typename T>
895  inline void clean(const simple_vector_ref<wsvector<T> *> &l, double eps) {
896  simple_vector_ref<wsvector<T> *>
897  *svr = const_cast<simple_vector_ref<wsvector<T> *> *>(&l);
898  wsvector<T>
899  *pv = const_cast<wsvector<T> *>((l.origin));
900  clean(*pv, eps);
901  svr->begin_ = vect_begin(*pv); svr->end_ = vect_end(*pv);
902  }
903 
904  template <typename T>
905  inline size_type nnz(const wsvector<T>& l) { return l.nb_stored(); }
906 
907  /*************************************************************************/
908  /* */
909  /* rsvector: sparse vector optimized for linear algebra operations. */
910  /* */
911  /*************************************************************************/
912 
913  template<typename T> struct elt_rsvector_ {
914  size_type c; T e;
915  /* e is initialized by default to avoid some false warnings of valgrind.
916  (from http://valgrind.org/docs/manual/mc-manual.html:
917 
918  When memory is read into the CPU's floating point registers, the
919  relevant V bits are read from memory and they are immediately
920  checked. If any are invalid, an uninitialised value error is
921  emitted. This precludes using the floating-point registers to copy
922  possibly-uninitialised memory, but simplifies Valgrind in that it
923  does not have to track the validity status of the floating-point
924  registers.
925  */
926  elt_rsvector_() : e(0) {}
927  elt_rsvector_(size_type cc) : c(cc), e(0) {}
928  elt_rsvector_(size_type cc, const T &ee) : c(cc), e(ee) {}
929  bool operator < (const elt_rsvector_ &a) const { return c < a.c; }
930  bool operator == (const elt_rsvector_ &a) const { return c == a.c; }
931  bool operator != (const elt_rsvector_ &a) const { return c != a.c; }
932  };
933 
934  template<typename T> struct rsvector_iterator {
935  typedef typename std::vector<elt_rsvector_<T> >::iterator IT;
936  typedef T value_type;
937  typedef value_type* pointer;
938  typedef value_type& reference;
939  typedef size_t size_type;
940  typedef ptrdiff_t difference_type;
941  typedef std::bidirectional_iterator_tag iterator_category;
942  typedef rsvector_iterator<T> iterator;
943 
944  IT it;
945 
946  reference operator *() const { return it->e; }
947  pointer operator->() const { return &(operator*()); }
948 
949  iterator &operator ++() { ++it; return *this; }
950  iterator operator ++(int) { iterator tmp = *this; ++(*this); return tmp; }
951  iterator &operator --() { --it; return *this; }
952  iterator operator --(int) { iterator tmp = *this; --(*this); return tmp; }
953 
954  bool operator ==(const iterator &i) const { return it == i.it; }
955  bool operator !=(const iterator &i) const { return !(i == *this); }
956 
957  size_type index() const { return it->c; }
958  rsvector_iterator() {}
959  rsvector_iterator(const IT &i) : it(i) {}
960  };
961 
962  template<typename T> struct rsvector_const_iterator {
963  typedef typename std::vector<elt_rsvector_<T> >::const_iterator IT;
964  typedef T value_type;
965  typedef const value_type* pointer;
966  typedef const value_type& reference;
967  typedef size_t size_type;
968  typedef ptrdiff_t difference_type;
969  typedef std::forward_iterator_tag iterator_category;
970  typedef rsvector_const_iterator<T> iterator;
971 
972  IT it;
973 
974  reference operator *() const { return it->e; }
975  pointer operator->() const { return &(operator*()); }
976  size_type index() const { return it->c; }
977 
978  iterator &operator ++() { ++it; return *this; }
979  iterator operator ++(int) { iterator tmp = *this; ++(*this); return tmp; }
980  iterator &operator --() { --it; return *this; }
981  iterator operator --(int) { iterator tmp = *this; --(*this); return tmp; }
982 
983  bool operator ==(const iterator &i) const { return it == i.it; }
984  bool operator !=(const iterator &i) const { return !(i == *this); }
985 
986  rsvector_const_iterator() {}
987  rsvector_const_iterator(const rsvector_iterator<T> &i) : it(i.it) {}
988  rsvector_const_iterator(const IT &i) : it(i) {}
989  };
990 
991  /**
992  sparse vector built upon std::vector. Read access is fast,
993  but insertion is O(n)
994  */
995  template<typename T> class rsvector : public std::vector<elt_rsvector_<T> > {
996  public:
997 
998  typedef std::vector<elt_rsvector_<T> > base_type_;
999  typedef typename base_type_::iterator iterator;
1000  typedef typename base_type_::const_iterator const_iterator;
1001  typedef typename base_type_::size_type size_type;
1002  typedef T value_type;
1003 
1004  protected:
1005  size_type nbl; // size of the vector
1006 
1007  public:
1008 
1009  void sup(size_type j);
1010  void base_resize(size_type n) { base_type_::resize(n); }
1011  void resize(size_type);
1012 
1013  ref_elt_vector<T, rsvector<T> > operator [](size_type c)
1014  { return ref_elt_vector<T, rsvector<T> >(this, c); }
1015 
1016  void w(size_type c, const T &e);
1017  void wa(size_type c, const T &e);
1018  T r(size_type c) const;
1019  void swap_indices(size_type i, size_type j);
1020 
1021  inline T operator [](size_type c) const { return r(c); }
1022 
1023  size_type nb_stored() const { return base_type_::size(); }
1024  size_type size() const { return nbl; }
1025  void clear() { base_type_::resize(0); }
1026  void swap(rsvector<T> &v)
1027  { std::swap(nbl, v.nbl); std::vector<elt_rsvector_<T> >::swap(v); }
1028 
1029  /* Constructeurs */
1030  explicit rsvector(size_type l) : nbl(l) { }
1031  rsvector() : nbl(0) { }
1032  };
1033 
1034  template <typename T>
1035  void rsvector<T>::swap_indices(size_type i, size_type j) {
1036  if (i > j) std::swap(i, j);
1037  if (i != j) {
1038  int situation = 0;
1039  elt_rsvector_<T> ei(i), ej(j), a;
1040  iterator it, ite, iti, itj;
1041  iti = std::lower_bound(this->begin(), this->end(), ei);
1042  if (iti != this->end() && iti->c == i) situation += 1;
1043  itj = std::lower_bound(this->begin(), this->end(), ej);
1044  if (itj != this->end() && itj->c == j) situation += 2;
1045 
1046  switch (situation) {
1047  case 1 : a = *iti; a.c = j; it = iti; ++it; ite = this->end();
1048  for (; it != ite && it->c <= j; ++it, ++iti) *iti = *it;
1049  *iti = a;
1050  break;
1051  case 2 : a = *itj; a.c = i; it = itj; ite = this->begin();
1052  if (it != ite) {
1053  --it;
1054  while (it->c >= i) { *itj = *it; --itj; if (it==ite) break; --it; }
1055  }
1056  *itj = a;
1057  break;
1058  case 3 : std::swap(iti->e, itj->e);
1059  break;
1060  }
1061  }
1062  }
1063 
1064  template <typename T> void rsvector<T>::sup(size_type j) {
1065  if (nb_stored() != 0) {
1066  elt_rsvector_<T> ev(j);
1067  iterator it = std::lower_bound(this->begin(), this->end(), ev);
1068  if (it != this->end() && it->c == j) {
1069  for (iterator ite = this->end() - 1; it != ite; ++it) *it = *(it+1);
1070  base_resize(nb_stored()-1);
1071  }
1072  }
1073  }
1074 
1075  template<typename T> void rsvector<T>::resize(size_type n) {
1076  if (n < nbl) {
1077  for (size_type i = 0; i < nb_stored(); ++i)
1078  if (base_type_::operator[](i).c >= n) { base_resize(i); break; }
1079  }
1080  nbl = n;
1081  }
1082 
1083  template <typename T> void rsvector<T>::w(size_type c, const T &e) {
1084  GMM_ASSERT2(c < nbl, "out of range");
1085  if (e == T(0)) sup(c);
1086  else {
1087  elt_rsvector_<T> ev(c, e);
1088  if (nb_stored() == 0) {
1089  base_type_::push_back(ev);
1090  }
1091  else {
1092  iterator it = std::lower_bound(this->begin(), this->end(), ev);
1093  if (it != this->end() && it->c == c) it->e = e;
1094  else {
1095  size_type ind = it - this->begin(), nb = this->nb_stored();
1096  if (nb - ind > 1100)
1097  GMM_WARNING2("Inefficient addition of element in rsvector with "
1098  << this->nb_stored() - ind << " non-zero entries");
1099  base_type_::push_back(ev);
1100  if (ind != nb) {
1101  it = this->begin() + ind;
1102  iterator ite = this->end(); --ite; iterator itee = ite;
1103  for (; ite != it; --ite) { --itee; *ite = *itee; }
1104  *it = ev;
1105  }
1106  }
1107  }
1108  }
1109  }
1110 
1111  template <typename T> void rsvector<T>::wa(size_type c, const T &e) {
1112  GMM_ASSERT2(c < nbl, "out of range");
1113  if (e != T(0)) {
1114  elt_rsvector_<T> ev(c, e);
1115  if (nb_stored() == 0) {
1116  base_type_::push_back(ev);
1117  }
1118  else {
1119  iterator it = std::lower_bound(this->begin(), this->end(), ev);
1120  if (it != this->end() && it->c == c) it->e += e;
1121  else {
1122  size_type ind = it - this->begin(), nb = this->nb_stored();
1123  if (nb - ind > 1100)
1124  GMM_WARNING2("Inefficient addition of element in rsvector with "
1125  << this->nb_stored() - ind << " non-zero entries");
1126  base_type_::push_back(ev);
1127  if (ind != nb) {
1128  it = this->begin() + ind;
1129  iterator ite = this->end(); --ite; iterator itee = ite;
1130  for (; ite != it; --ite) { --itee; *ite = *itee; }
1131  *it = ev;
1132  }
1133  }
1134  }
1135  }
1136  }
1137 
1138  template <typename T> T rsvector<T>::r(size_type c) const {
1139  GMM_ASSERT2(c < nbl, "out of range. Index " << c
1140  << " for a length of " << nbl);
1141  if (nb_stored() != 0) {
1142  elt_rsvector_<T> ev(c);
1143  const_iterator it = std::lower_bound(this->begin(), this->end(), ev);
1144  if (it != this->end() && it->c == c) return it->e;
1145  }
1146  return T(0);
1147  }
1148 
1149  template <typename T> struct linalg_traits<rsvector<T> > {
1150  typedef rsvector<T> this_type;
1151  typedef this_type origin_type;
1152  typedef linalg_false is_reference;
1153  typedef abstract_vector linalg_type;
1154  typedef T value_type;
1155  typedef ref_elt_vector<T, rsvector<T> > reference;
1156  typedef rsvector_iterator<T> iterator;
1157  typedef rsvector_const_iterator<T> const_iterator;
1158  typedef abstract_sparse storage_type;
1159  typedef linalg_true index_sorted;
1160  static size_type size(const this_type &v) { return v.size(); }
1161  static iterator begin(this_type &v) { return iterator(v.begin()); }
1162  static const_iterator begin(const this_type &v)
1163  { return const_iterator(v.begin()); }
1164  static iterator end(this_type &v) { return iterator(v.end()); }
1165  static const_iterator end(const this_type &v)
1166  { return const_iterator(v.end()); }
1167  static origin_type* origin(this_type &v) { return &v; }
1168  static const origin_type* origin(const this_type &v) { return &v; }
1169  static void clear(origin_type* o, const iterator &, const iterator &)
1170  { o->clear(); }
1171  static void do_clear(this_type &v) { v.clear(); }
1172  static value_type access(const origin_type *o, const const_iterator &,
1173  const const_iterator &, size_type i)
1174  { return (*o)[i]; }
1175  static reference access(origin_type *o, const iterator &, const iterator &,
1176  size_type i)
1177  { return (*o)[i]; }
1178  static void resize(this_type &v, size_type n) { v.resize(n); }
1179  };
1180 
1181  template<typename T> std::ostream &operator <<
1182  (std::ostream &o, const rsvector<T>& v) { gmm::write(o,v); return o; }
1183 
1184  /******* Optimized operations for rsvector<T> ****************************/
1185 
1186  template <typename T> inline void copy(const rsvector<T> &v1,
1187  rsvector<T> &v2) {
1188  GMM_ASSERT2(vect_size(v1) == vect_size(v2), "dimensions mismatch");
1189  v2 = v1;
1190  }
1191  template <typename T> inline
1192  void copy(const rsvector<T> &v1, const simple_vector_ref<rsvector<T> *> &v2){
1193  simple_vector_ref<rsvector<T> *>
1194  *svr = const_cast<simple_vector_ref<rsvector<T> *> *>(&v2);
1195  rsvector<T>
1196  *pv = const_cast<rsvector<T> *>((v2.origin));
1197  GMM_ASSERT2(vect_size(v1) == vect_size(v2), "dimensions mismatch");
1198  *pv = v1; svr->begin_ = vect_begin(*pv); svr->end_ = vect_end(*pv);
1199  }
1200  template <typename T> inline
1201  void copy(const simple_vector_ref<const rsvector<T> *> &v1,
1202  rsvector<T> &v2)
1203  { copy(*(v1.origin), v2); }
1204  template <typename T> inline
1205  void copy(const simple_vector_ref<rsvector<T> *> &v1, rsvector<T> &v2)
1206  { copy(*(v1.origin), v2); }
1207 
1208  template <typename V, typename T> inline void add(const V &v1,
1209  rsvector<T> &v2) {
1210  if ((const void *)(&v1) != (const void *)(&v2)) {
1211  GMM_ASSERT2(vect_size(v1) == vect_size(v2), "dimensions mismatch");
1212  add_rsvector(v1, v2, typename linalg_traits<V>::storage_type());
1213  }
1214  }
1215 
1216  template <typename V, typename T>
1217  inline void add_rsvector(const V &v1, rsvector<T> &v2, abstract_dense)
1218  { add(v1, v2, abstract_dense(), abstract_sparse()); }
1219 
1220  template <typename V, typename T>
1221  inline void add_rsvector(const V &v1, rsvector<T> &v2, abstract_skyline)
1222  { add(v1, v2, abstract_skyline(), abstract_sparse()); }
1223 
1224  template <typename V, typename T>
1225  void add_rsvector(const V &v1, rsvector<T> &v2, abstract_sparse) {
1226  add_rsvector(v1, v2, typename linalg_traits<V>::index_sorted());
1227  }
1228 
1229  template <typename V, typename T>
1230  void add_rsvector(const V &v1, rsvector<T> &v2, linalg_false) {
1231  add(v1, v2, abstract_sparse(), abstract_sparse());
1232  }
1233 
1234  template <typename V, typename T>
1235  void add_rsvector(const V &v1, rsvector<T> &v2, linalg_true) {
1236  typename linalg_traits<V>::const_iterator it1 = vect_const_begin(v1),
1237  ite1 = vect_const_end(v1);
1238  typename rsvector<T>::iterator it2 = v2.begin(), ite2 = v2.end(), it3;
1239  size_type nbc = 0, old_nbc = v2.nb_stored();
1240  for (; it1 != ite1 && it2 != ite2 ; ++nbc)
1241  if (it1.index() == it2->c) { ++it1; ++it2; }
1242  else if (it1.index() < it2->c) ++it1; else ++it2;
1243  for (; it1 != ite1; ++it1) ++nbc;
1244  for (; it2 != ite2; ++it2) ++nbc;
1245 
1246  v2.base_resize(nbc);
1247  it3 = v2.begin() + old_nbc;
1248  it2 = v2.end();
1249  ite2 = v2.begin();
1250  it1 = vect_end(v1);
1251  ite1 = vect_const_begin(v1);
1252  while (it1 != ite1 && it2 != ite2 && it3 != ite2){
1253  --it3;
1254  --it1;
1255  --it2;
1256  if (it3->c > it1.index()) {
1257  *it2 = *it3;
1258  ++it1;
1259  }
1260  else if (it3->c == it1.index()) {
1261  *it2=*it3;
1262  it2->e+=*it1;
1263  }
1264  else {
1265  it2->c = it1.index();
1266  it2->e = *it1; ++it3;
1267  }
1268  }
1269  while (it1 != ite1 && it2 != ite2) {
1270  --it1;
1271  --it2;
1272  it2->c = it1.index();
1273  it2->e = *it1;
1274  }
1275  }
1276 
1277  template <typename V, typename T> void copy(const V &v1, rsvector<T> &v2) {
1278  if ((const void *)(&v1) != (const void *)(&v2)) {
1279  GMM_ASSERT2(vect_size(v1) == vect_size(v2), "dimensions mismatch");
1280  if (same_origin(v1, v2))
1281  GMM_WARNING2("a conflict is possible in vector copy\n");
1282  copy_rsvector(v1, v2, typename linalg_traits<V>::storage_type());
1283  }
1284  }
1285 
1286  template <typename V, typename T>
1287  void copy_rsvector(const V &v1, rsvector<T> &v2, abstract_dense)
1288  { copy_vect(v1, v2, abstract_dense(), abstract_sparse()); }
1289 
1290  template <typename V, typename T>
1291  void copy_rsvector(const V &v1, rsvector<T> &v2, abstract_skyline)
1292  { copy_vect(v1, v2, abstract_skyline(), abstract_sparse()); }
1293 
1294  template <typename V, typename T>
1295  void copy_rsvector(const V &v1, rsvector<T> &v2, abstract_sparse) {
1296  copy_rsvector(v1, v2, typename linalg_traits<V>::index_sorted());
1297  }
1298 
1299  template <typename V, typename T2>
1300  void copy_rsvector(const V &v1, rsvector<T2> &v2, linalg_true) {
1301  typedef typename linalg_traits<V>::value_type T1;
1302  typename linalg_traits<V>::const_iterator it = vect_const_begin(v1),
1303  ite = vect_const_end(v1);
1304  v2.base_resize(nnz(v1));
1305  typename rsvector<T2>::iterator it2 = v2.begin();
1306  size_type nn = 0;
1307  for (; it != ite; ++it)
1308  if ((*it) != T1(0)) { it2->c = it.index(); it2->e = *it; ++it2; ++nn; }
1309  v2.base_resize(nn);
1310  }
1311 
1312  template <typename V, typename T2>
1313  void copy_rsvector(const V &v1, rsvector<T2> &v2, linalg_false) {
1314  typedef typename linalg_traits<V>::value_type T1;
1315  typename linalg_traits<V>::const_iterator it = vect_const_begin(v1),
1316  ite = vect_const_end(v1);
1317  v2.base_resize(nnz(v1));
1318  typename rsvector<T2>::iterator it2 = v2.begin();
1319  size_type nn = 0;
1320  for (; it != ite; ++it)
1321  if ((*it) != T1(0)) { it2->c = it.index(); it2->e = *it; ++it2; ++nn; }
1322  v2.base_resize(nn);
1323  std::sort(v2.begin(), v2.end());
1324  }
1325 
1326  template <typename T> inline void clean(rsvector<T> &v, double eps) {
1327  typedef typename number_traits<T>::magnitude_type R;
1328  typename rsvector<T>::iterator it = v.begin(), ite = v.end();
1329  for (; it != ite; ++it) if (gmm::abs((*it).e) <= eps) break;
1330  if (it != ite) {
1331  typename rsvector<T>::iterator itc = it;
1332  size_type erased = 1;
1333  for (++it; it != ite; ++it)
1334  { *itc = *it; if (gmm::abs((*it).e) <= R(eps)) ++erased; else ++itc; }
1335  v.base_resize(v.nb_stored() - erased);
1336  }
1337  }
1338 
1339  template <typename T>
1340  inline void clean(const simple_vector_ref<rsvector<T> *> &l, double eps) {
1341  simple_vector_ref<rsvector<T> *>
1342  *svr = const_cast<simple_vector_ref<rsvector<T> *> *>(&l);
1343  rsvector<T>
1344  *pv = const_cast<rsvector<T> *>((l.origin));
1345  clean(*pv, eps);
1346  svr->begin_ = vect_begin(*pv); svr->end_ = vect_end(*pv);
1347  }
1348 
1349  template <typename T>
1350  inline size_type nnz(const rsvector<T>& l) { return l.nb_stored(); }
1351 
1352  /*************************************************************************/
1353  /* */
1354  /* Class slvector: 'sky-line' vector. */
1355  /* */
1356  /*************************************************************************/
1357 
1358  template<typename T> struct slvector_iterator {
1359  typedef T value_type;
1360  typedef T *pointer;
1361  typedef T &reference;
1362  typedef ptrdiff_t difference_type;
1363  typedef std::random_access_iterator_tag iterator_category;
1364  typedef size_t size_type;
1365  typedef slvector_iterator<T> iterator;
1366  typedef typename std::vector<T>::iterator base_iterator;
1367 
1368  base_iterator it;
1369  size_type shift;
1370 
1371 
1372  iterator &operator ++()
1373  { ++it; ++shift; return *this; }
1374  iterator &operator --()
1375  { --it; --shift; return *this; }
1376  iterator operator ++(int)
1377  { iterator tmp = *this; ++(*(this)); return tmp; }
1378  iterator operator --(int)
1379  { iterator tmp = *this; --(*(this)); return tmp; }
1380  iterator &operator +=(difference_type i)
1381  { it += i; shift += i; return *this; }
1382  iterator &operator -=(difference_type i)
1383  { it -= i; shift -= i; return *this; }
1384  iterator operator +(difference_type i) const
1385  { iterator tmp = *this; return (tmp += i); }
1386  iterator operator -(difference_type i) const
1387  { iterator tmp = *this; return (tmp -= i); }
1388  difference_type operator -(const iterator &i) const
1389  { return it - i.it; }
1390 
1391  reference operator *() const
1392  { return *it; }
1393  reference operator [](int ii)
1394  { return *(it + ii); }
1395 
1396  bool operator ==(const iterator &i) const
1397  { return it == i.it; }
1398  bool operator !=(const iterator &i) const
1399  { return !(i == *this); }
1400  bool operator < (const iterator &i) const
1401  { return it < i.it; }
1402  bool operator > (const iterator &i) const
1403  { return it > i.it; }
1404  bool operator >=(const iterator &i) const
1405  { return it >= i.it; }
1406  size_type index() const { return shift; }
1407 
1408  slvector_iterator() {}
1409  slvector_iterator(const base_iterator &iter, size_type s)
1410  : it(iter), shift(s) {}
1411  };
1412 
1413  template<typename T> struct slvector_const_iterator {
1414  typedef T value_type;
1415  typedef const T *pointer;
1416  typedef value_type reference;
1417  typedef ptrdiff_t difference_type;
1418  typedef std::random_access_iterator_tag iterator_category;
1419  typedef size_t size_type;
1420  typedef slvector_const_iterator<T> iterator;
1421  typedef typename std::vector<T>::const_iterator base_iterator;
1422 
1423  base_iterator it;
1424  size_type shift;
1425 
1426 
1427  iterator &operator ++()
1428  { ++it; ++shift; return *this; }
1429  iterator &operator --()
1430  { --it; --shift; return *this; }
1431  iterator operator ++(int)
1432  { iterator tmp = *this; ++(*(this)); return tmp; }
1433  iterator operator --(int)
1434  { iterator tmp = *this; --(*(this)); return tmp; }
1435  iterator &operator +=(difference_type i)
1436  { it += i; shift += i; return *this; }
1437  iterator &operator -=(difference_type i)
1438  { it -= i; shift -= i; return *this; }
1439  iterator operator +(difference_type i) const
1440  { iterator tmp = *this; return (tmp += i); }
1441  iterator operator -(difference_type i) const
1442  { iterator tmp = *this; return (tmp -= i); }
1443  difference_type operator -(const iterator &i) const
1444  { return it - i.it; }
1445 
1446  value_type operator *() const
1447  { return *it; }
1448  value_type operator [](int ii)
1449  { return *(it + ii); }
1450 
1451  bool operator ==(const iterator &i) const
1452  { return it == i.it; }
1453  bool operator !=(const iterator &i) const
1454  { return !(i == *this); }
1455  bool operator < (const iterator &i) const
1456  { return it < i.it; }
1457  bool operator > (const iterator &i) const
1458  { return it > i.it; }
1459  bool operator >=(const iterator &i) const
1460  { return it >= i.it; }
1461  size_type index() const { return shift; }
1462 
1463  slvector_const_iterator() {}
1464  slvector_const_iterator(const slvector_iterator<T>& iter)
1465  : it(iter.it), shift(iter.shift) {}
1466  slvector_const_iterator(const base_iterator &iter, size_type s)
1467  : it(iter), shift(s) {}
1468  };
1469 
1470 
1471  /** skyline vector.
1472  */
1473  template <typename T> class slvector {
1474 
1475  public :
1476  typedef slvector_iterator<T> iterators;
1477  typedef slvector_const_iterator<T> const_iterators;
1478  typedef typename std::vector<T>::size_type size_type;
1479  typedef T value_type;
1480 
1481  protected :
1482  std::vector<T> data;
1483  size_type shift;
1484  size_type size_;
1485 
1486 
1487  public :
1488 
1489  size_type size() const { return size_; }
1490  size_type first() const { return shift; }
1491  size_type last() const { return shift + data.size(); }
1492  ref_elt_vector<T, slvector<T> > operator [](size_type c)
1493  { return ref_elt_vector<T, slvector<T> >(this, c); }
1494 
1495  typename std::vector<T>::iterator data_begin() { return data.begin(); }
1496  typename std::vector<T>::iterator data_end() { return data.end(); }
1497  typename std::vector<T>::const_iterator data_begin() const
1498  { return data.begin(); }
1499  typename std::vector<T>::const_iterator data_end() const
1500  { return data.end(); }
1501 
1502  void w(size_type c, const T &e);
1503  void wa(size_type c, const T &e);
1504  T r(size_type c) const {
1505  GMM_ASSERT2(c < size_, "out of range");
1506  if (c < shift || c >= shift + data.size()) return T(0);
1507  return data[c - shift];
1508  }
1509 
1510  inline T operator [](size_type c) const { return r(c); }
1511  void resize(size_type);
1512  void clear() { data.resize(0); shift = 0; }
1513  void swap(slvector<T> &v) {
1514  std::swap(data, v.data);
1515  std::swap(shift, v.shift);
1516  std::swap(size_, v.size_);
1517  }
1518 
1519 
1520  slvector() : data(0), shift(0), size_(0) {}
1521  explicit slvector(size_type l) : data(0), shift(0), size_(l) {}
1522  slvector(size_type l, size_type d, size_type s)
1523  : data(d), shift(s), size_(l) {}
1524 
1525  };
1526 
1527  template<typename T> void slvector<T>::resize(size_type n) {
1528  if (n < last()) {
1529  if (shift >= n) clear(); else { data.resize(n-shift); }
1530  }
1531  size_ = n;
1532  }
1533 
1534  template<typename T> void slvector<T>::w(size_type c, const T &e) {
1535  GMM_ASSERT2(c < size_, "out of range");
1536  size_type s = data.size();
1537  if (!s) { data.resize(1); shift = c; }
1538  else if (c < shift) {
1539  data.resize(s + shift - c);
1540  typename std::vector<T>::iterator it = data.begin(),it2=data.end()-1;
1541  typename std::vector<T>::iterator it3 = it2 - shift + c;
1542  for (; it3 >= it; --it3, --it2) *it2 = *it3;
1543  std::fill(it, it + shift - c, T(0));
1544  shift = c;
1545  }
1546  else if (c >= shift + s) {
1547  data.resize(c - shift + 1, T(0));
1548  // std::fill(data.begin() + s, data.end(), T(0));
1549  }
1550  data[c - shift] = e;
1551  }
1552 
1553  template<typename T> void slvector<T>::wa(size_type c, const T &e) {
1554  GMM_ASSERT2(c < size_, "out of range");
1555  size_type s = data.size();
1556  if (!s) { data.resize(1, e); shift = c; return; }
1557  else if (c < shift) {
1558  data.resize(s + shift - c);
1559  typename std::vector<T>::iterator it = data.begin(),it2=data.end()-1;
1560  typename std::vector<T>::iterator it3 = it2 - shift + c;
1561  for (; it3 >= it; --it3, --it2) *it2 = *it3;
1562  std::fill(it, it + shift - c, T(0));
1563  shift = c;
1564  data[c - shift] = e;
1565  return;
1566  }
1567  else if (c >= shift + s) {
1568  data.resize(c - shift + 1, T(0));
1569  data[c - shift] = e;
1570  return;
1571  // std::fill(data.begin() + s, data.end(), T(0));
1572  }
1573  data[c - shift] += e;
1574  }
1575 
1576 
1577  template <typename T> struct linalg_traits<slvector<T> > {
1578  typedef slvector<T> this_type;
1579  typedef this_type origin_type;
1580  typedef linalg_false is_reference;
1581  typedef abstract_vector linalg_type;
1582  typedef T value_type;
1583  typedef ref_elt_vector<T, slvector<T> > reference;
1584  typedef slvector_iterator<T> iterator;
1585  typedef slvector_const_iterator<T> const_iterator;
1586  typedef abstract_skyline storage_type;
1587  typedef linalg_true index_sorted;
1588  static size_type size(const this_type &v) { return v.size(); }
1589  static iterator begin(this_type &v)
1590  { return iterator(v.data_begin(), v.first()); }
1591  static const_iterator begin(const this_type &v)
1592  { return const_iterator(v.data_begin(), v.first()); }
1593  static iterator end(this_type &v)
1594  { return iterator(v.data_end(), v.last()); }
1595  static const_iterator end(const this_type &v)
1596  { return const_iterator(v.data_end(), v.last()); }
1597  static origin_type* origin(this_type &v) { return &v; }
1598  static const origin_type* origin(const this_type &v) { return &v; }
1599  static void clear(origin_type* o, const iterator &, const iterator &)
1600  { o->clear(); }
1601  static void do_clear(this_type &v) { v.clear(); }
1602  static value_type access(const origin_type *o, const const_iterator &,
1603  const const_iterator &, size_type i)
1604  { return (*o)[i]; }
1605  static reference access(origin_type *o, const iterator &, const iterator &,
1606  size_type i)
1607  { return (*o)[i]; }
1608  static void resize(this_type &v, size_type n) { v.resize(n); }
1609  };
1610 
1611  template<typename T> std::ostream &operator <<
1612  (std::ostream &o, const slvector<T>& v) { gmm::write(o,v); return o; }
1613 
1614  template <typename T>
1615  inline size_type nnz(const slvector<T>& l) { return l.last() - l.first(); }
1616 
1617 }
1618 
1619 namespace std {
1620  template <typename T> void swap(gmm::wsvector<T> &v, gmm::wsvector<T> &w)
1621  { v.swap(w);}
1622  template <typename T> void swap(gmm::rsvector<T> &v, gmm::rsvector<T> &w)
1623  { v.swap(w);}
1624  template <typename T> void swap(gmm::slvector<T> &v, gmm::slvector<T> &w)
1625  { v.swap(w);}
1626 }
1627 
1628 
1629 
1630 #endif /* GMM_VECTOR_H__ */
gmm::dsvector
Sparse vector built on distribution sort principle.
Definition: gmm_vector.h:324
gmm::rsvector
sparse vector built upon std::vector.
Definition: gmm_vector.h:995
gmm::slvector
skyline vector.
Definition: gmm_vector.h:1473
gmm::wsvector
sparse vector built upon std::map.
Definition: gmm_vector.h:752
gmm::nnz
size_type nnz(const L &l)
count the number of non-zero entries of a vector or matrix.
Definition: gmm_blas.h:68
gmm::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:976
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:58
gmm::resize
void resize(V &v, size_type n)
*‍/
Definition: gmm_blas.h:209
gmm::clean
void clean(L &l, double threshold)
Clean a vector or matrix (replace near-zero entries with zeroes).
gmm::add
void add(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:1275
gmm_interface.h
gmm interface for STL vectors.
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:48