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