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