GetFEM  5.4.3
gmm_interface.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 
32 
33 /**@file gmm_interface.h
34  @author Yves Renard <Yves.Renard@insa-lyon.fr>
35  @date October 13, 2002.
36  @brief gmm interface for STL vectors.
37 */
38 
39 #ifndef GMM_INTERFACE_H__
40 #define GMM_INTERFACE_H__
41 
42 #include "gmm_blas.h"
43 #include "gmm_sub_index.h"
44 
45 namespace gmm {
46 
47  /* ********************************************************************* */
48  /* */
49  /* What is needed for a Vector type : */
50  /* Vector v(n) defines a vector with n components. */
51  /* v[i] allows to access to the ith component of v. */
52  /* linalg_traits<Vector> should be filled with appropriate definitions */
53  /* */
54  /* for a dense vector : the minimum is two random iterators (begin and */
55  /* end) and a pointer to a valid origin. */
56  /* for a sparse vector : the minimum is two forward iterators, with */
57  /* a method it.index() which gives the index of */
58  /* a non zero element, an interface object */
59  /* should describe the method to add new non */
60  /* zero element, and a pointer to a valid */
61  /* origin. */
62  /* */
63  /* What is needed for a Matrix type : */
64  /* Matrix m(n, m) defines a matrix with n rows and m columns. */
65  /* m(i, j) allows to access to the element at row i and column j. */
66  /* linalg_traits<Matrix> should be filled with appropriate definitions */
67  /* */
68  /* What is needed for an iterator on dense vector */
69  /* to be standard random access iterator */
70  /* */
71  /* What is needed for an iterator on a sparse vector */
72  /* to be a standard bidirectional iterator */
73  /* elt should be sorted with increasing indices. */
74  /* it.index() gives the index of the non-zero element. */
75  /* */
76  /* Remark : If original iterators are not convenient, they could be */
77  /* redefined and interfaced in linalg_traits<Vector> without changing */
78  /* the original Vector type. */
79  /* */
80  /* ********************************************************************* */
81 
82  /* ********************************************************************* */
83  /* Simple references on vectors */
84  /* ********************************************************************* */
85 
86  template <typename PT> struct simple_vector_ref {
87  typedef simple_vector_ref<PT> this_type;
88  typedef typename std::iterator_traits<PT>::value_type V;
89  typedef V * CPT;
90  typedef typename std::iterator_traits<PT>::reference ref_V;
91  typedef typename linalg_traits<this_type>::iterator iterator;
92  typedef typename linalg_traits<this_type>::reference reference;
93  typedef typename linalg_traits<this_type>::porigin_type porigin_type;
94 
95  iterator begin_, end_;
96  porigin_type origin;
97  size_type size_;
98 
99  simple_vector_ref(ref_V v) : begin_(vect_begin(const_cast<V&>(v))),
100  end_(vect_end(const_cast<V&>(v))),
101  origin(linalg_origin(const_cast<V&>(v))),
102  size_(vect_size(v)) {}
103 
104  simple_vector_ref(const simple_vector_ref<CPT> &cr)
105  : begin_(cr.begin_),end_(cr.end_),origin(cr.origin),size_(cr.size_) {}
106 
107  simple_vector_ref(void) {}
108 
109  reference operator[](size_type i) const
110  { return linalg_traits<V>::access(origin, begin_, end_, i); }
111  };
112 
113  template <typename IT, typename ORG, typename PT> inline
114  void set_to_begin(IT &it, ORG o, simple_vector_ref<PT> *,linalg_modifiable) {
115  typedef typename linalg_traits<simple_vector_ref<PT> >::V_reference ref_t;
116  set_to_begin(it, o, PT(), ref_t());
117  }
118 
119  template <typename IT, typename ORG, typename PT> inline
120  void set_to_begin(IT &it, ORG o, const simple_vector_ref<PT> *,
121  linalg_modifiable) {
122  typedef typename linalg_traits<simple_vector_ref<PT> >::V_reference ref_t;
123  set_to_begin(it, o, PT(), ref_t());
124  }
125 
126  template <typename IT, typename ORG, typename PT> inline
127  void set_to_end(IT &it, ORG o, simple_vector_ref<PT> *, linalg_modifiable) {
128  typedef typename linalg_traits<simple_vector_ref<PT> >::V_reference ref_t;
129  set_to_end(it, o, PT(), ref_t());
130  }
131 
132  template <typename IT, typename ORG, typename PT> inline
133  void set_to_end(IT &it, ORG o, const simple_vector_ref<PT> *,
134  linalg_modifiable) {
135  typedef typename linalg_traits<simple_vector_ref<PT> >::V_reference ref_t;
136  set_to_end(it, o, PT(), ref_t());
137  }
138 
139 
140  template <typename PT> struct linalg_traits<simple_vector_ref<PT> > {
141  typedef simple_vector_ref<PT> this_type;
142  typedef this_type *pthis_type;
143  typedef typename std::iterator_traits<PT>::value_type V;
144  typedef typename linalg_traits<V>::origin_type origin_type;
145  typedef V *pV;
146  typedef typename linalg_traits<V>::is_reference V_reference;
147  typedef typename which_reference<PT>::is_reference is_reference;
148  typedef abstract_vector linalg_type;
149  typedef typename linalg_traits<V>::value_type value_type;
150  typedef typename select_ref<value_type, typename
151  linalg_traits<V>::reference, PT>::ref_type reference;
152  typedef typename select_ref<const origin_type *, origin_type *,
153  PT>::ref_type porigin_type;
154  typedef typename select_ref<typename linalg_traits<V>::const_iterator,
155  typename linalg_traits<V>::iterator, PT>::ref_type iterator;
156  typedef typename linalg_traits<V>::const_iterator const_iterator;
157  typedef typename linalg_traits<V>::storage_type storage_type;
158  typedef linalg_true index_sorted;
159  static size_type size(const this_type &v) { return v.size_; }
160  static inline iterator begin(this_type &v) {
161  iterator it = v.begin_;
162  set_to_begin(it, v.origin, pthis_type(), is_reference());
163  return it;
164  }
165  static inline const_iterator begin(const this_type &v) {
166  const_iterator it = v.begin_;
167  set_to_begin(it, v.origin, pthis_type(), is_reference());
168  return it;
169  }
170  static inline iterator end(this_type &v) {
171  iterator it = v.end_;
172  set_to_end(it, v.origin, pthis_type(), is_reference());
173  return it;
174  }
175  static inline const_iterator end(const this_type &v) {
176  const_iterator it = v.end_;
177  set_to_end(it, v.origin, pthis_type(), is_reference());
178  return it;
179  }
180  static origin_type* origin(this_type &v) { return v.origin; }
181  static const origin_type* origin(const this_type &v) { return v.origin; }
182  static void clear(origin_type* o, const iterator &it, const iterator &ite)
183  { linalg_traits<V>::clear(o, it, ite); }
184  static void do_clear(this_type &v) { clear(v.origin, v.begin_, v.end_); }
185  static value_type access(const origin_type *o, const const_iterator &it,
186  const const_iterator &ite, size_type i)
187  { return linalg_traits<V>::access(o, it, ite, i); }
188  static reference access(origin_type *o, const iterator &it,
189  const iterator &ite, size_type i)
190  { return linalg_traits<V>::access(o, it, ite, i); }
191  };
192 
193  template <typename PT>
194  std::ostream &operator << (std::ostream &o, const simple_vector_ref<PT>& v)
195  { gmm::write(o,v); return o; }
196 
197  template <typename T, typename alloc>
198  simple_vector_ref<const std::vector<T,alloc> *>
199  vref(const std::vector<T, alloc> &vv)
200  { return simple_vector_ref<const std::vector<T,alloc> *>(vv); }
201 
202 
203  /* ********************************************************************* */
204  /* */
205  /* Traits for S.T.L. object */
206  /* */
207  /* ********************************************************************* */
208 
209  template <typename T, typename alloc>
210  struct linalg_traits<std::vector<T, alloc> > {
211  typedef std::vector<T, alloc> this_type;
212  typedef this_type origin_type;
213  typedef linalg_false is_reference;
214  typedef abstract_vector linalg_type;
215  typedef T value_type;
216  typedef T& reference;
217  typedef typename this_type::iterator iterator;
218  typedef typename this_type::const_iterator const_iterator;
219  typedef abstract_dense storage_type;
220  typedef linalg_true index_sorted;
221  static size_type size(const this_type &v) { return v.size(); }
222  static iterator begin(this_type &v) { return v.begin(); }
223  static const_iterator begin(const this_type &v) { return v.begin(); }
224  static iterator end(this_type &v) { return v.end(); }
225  static const_iterator end(const this_type &v) { return v.end(); }
226  static origin_type* origin(this_type &v) { return &v; }
227  static const origin_type* origin(const this_type &v) { return &v; }
228  static void clear(origin_type*, const iterator &it, const iterator &ite)
229  { std::fill(it, ite, value_type(0)); }
230  static void do_clear(this_type &v) { std::fill(v.begin(), v.end(), T(0)); }
231  static value_type access(const origin_type *, const const_iterator &it,
232  const const_iterator &, size_type i)
233  { return it[i]; }
234  static reference access(origin_type *, const iterator &it,
235  const iterator &, size_type i)
236  { return it[i]; }
237  static void resize(this_type &v, size_type n) { v.resize(n); }
238  };
239 
240 
241 
242  template <typename T>
243  inline size_type nnz(const std::vector<T>& l) { return l.size(); }
244 
245  /* ********************************************************************* */
246  /* */
247  /* Traits for ref objects */
248  /* */
249  /* ********************************************************************* */
250 
251  template <typename IT, typename V>
252  struct tab_ref_with_origin : public gmm::tab_ref<IT> {
253  typedef tab_ref_with_origin<IT, V> this_type;
254  // next line replaced by the 4 following lines in order to please aCC
255  //typedef typename linalg_traits<this_type>::porigin_type porigin_type;
256  typedef typename linalg_traits<V>::origin_type origin_type;
257  typedef typename std::iterator_traits<IT>::pointer PT;
258  typedef typename select_ref<const origin_type *, origin_type *,
259  PT>::ref_type porigin_type;
260 
261 
262  porigin_type origin;
263 
264  tab_ref_with_origin(void) {}
265  template <class PT> tab_ref_with_origin(const IT &b, const IT &e, PT p)
266  : gmm::tab_ref<IT>(b,e), origin(porigin_type(p)) {}
267  tab_ref_with_origin(const IT &b, const IT &e, porigin_type p)
268  : gmm::tab_ref<IT>(b,e), origin(p) {}
269 
270  tab_ref_with_origin(const V &v, const sub_interval &si)
271  : gmm::tab_ref<IT>(vect_begin(const_cast<V&>(v))+si.min,
272  vect_begin(const_cast<V&>(v))+si.max),
273  origin(linalg_origin(const_cast<V&>(v))) {}
274  tab_ref_with_origin(V &v, const sub_interval &si)
275  : gmm::tab_ref<IT>(vect_begin(const_cast<V&>(v))+si.min,
276  vect_begin(const_cast<V&>(v))+si.max),
277  origin(linalg_origin(const_cast<V&>(v))) {}
278  };
279 
280  template <typename IT, typename V>
281  struct linalg_traits<tab_ref_with_origin<IT, V> > {
282  typedef typename std::iterator_traits<IT>::pointer PT;
283  typedef typename linalg_traits<V>::origin_type origin_type;
284  typedef tab_ref_with_origin<IT, V> this_type;
285  typedef typename which_reference<PT>::is_reference is_reference;
286  typedef abstract_vector linalg_type;
287  typedef typename select_ref<const origin_type *, origin_type *,
288  PT>::ref_type porigin_type;
289  typedef typename std::iterator_traits<IT>::value_type value_type;
290  typedef typename std::iterator_traits<IT>::reference reference;
291  typedef typename this_type::iterator iterator;
292  typedef typename this_type::iterator const_iterator;
293  typedef abstract_dense storage_type;
294  typedef linalg_true index_sorted;
295  static size_type size(const this_type &v) { return v.size(); }
296  static iterator begin(this_type &v) { return v.begin(); }
297  static const_iterator begin(const this_type &v) { return v.begin(); }
298  static iterator end(this_type &v) { return v.end(); }
299  static const_iterator end(const this_type &v) { return v.end(); }
300  static origin_type* origin(this_type &v) { return v.origin; }
301  static const origin_type* origin(const this_type &v) { return v.origin; }
302  static void clear(origin_type*, const iterator &it, const iterator &ite)
303  { std::fill(it, ite, value_type(0)); }
304  static inline void do_clear(this_type &v)
305  { std::fill(v.begin(), v.end(), value_type(0)); }
306  static value_type access(const origin_type *, const const_iterator &it,
307  const const_iterator &, size_type i)
308  { return it[i]; }
309  static reference access(origin_type *, const iterator &it,
310  const iterator &, size_type i)
311  { return it[i]; }
312  };
313 
314  template <typename IT, typename V> std::ostream &operator <<
315  (std::ostream &o, const tab_ref_with_origin<IT, V>& m)
316  { gmm::write(o,m); return o; }
317 
318 
319  template <typename IT, typename V>
320  struct tab_ref_reg_spaced_with_origin : public gmm::tab_ref_reg_spaced<IT> {
321  typedef tab_ref_reg_spaced_with_origin<IT, V> this_type;
322  typedef typename linalg_traits<this_type>::porigin_type porigin_type;
323 
324  porigin_type origin;
325 
326  tab_ref_reg_spaced_with_origin(void) {}
327  tab_ref_reg_spaced_with_origin(const IT &b, size_type n, size_type s,
328  const porigin_type p)
329  : gmm::tab_ref_reg_spaced<IT>(b,n,s), origin(p) {}
330  tab_ref_reg_spaced_with_origin(const V &v, const sub_slice &si)
331  : gmm::tab_ref_reg_spaced<IT>(vect_begin(const_cast<V&>(v)) + si.min,
332  si.N, (si.max - si.min)/si.N),
333  origin(linalg_origin(const_cast<V&>(v))) {}
334  tab_ref_reg_spaced_with_origin(V &v, const sub_slice &si)
335  : gmm::tab_ref_reg_spaced<IT>(vect_begin(const_cast<V&>(v)) + si.min,
336  si.N, (si.max - si.min)/si.N),
337  origin(linalg_origin(const_cast<V&>(v))) {}
338  };
339 
340  template <typename IT, typename V>
341  struct linalg_traits<tab_ref_reg_spaced_with_origin<IT, V> > {
342  typedef typename std::iterator_traits<IT>::pointer PT;
343  typedef tab_ref_reg_spaced_with_origin<IT, V> this_type;
344  typedef typename linalg_traits<V>::origin_type origin_type;
345  typedef typename select_ref<const origin_type *, origin_type *,
346  PT>::ref_type porigin_type;
347  typedef typename which_reference<PT>::is_reference is_reference;
348  typedef abstract_vector linalg_type;
349  typedef typename std::iterator_traits<IT>::value_type value_type;
350  typedef typename std::iterator_traits<IT>::reference reference;
351  typedef typename this_type::iterator iterator;
352  typedef typename this_type::iterator const_iterator;
353  typedef abstract_dense storage_type;
354  typedef linalg_true index_sorted;
355  static size_type size(const this_type &v) { return v.size(); }
356  static iterator begin(this_type &v) { return v.begin(); }
357  static const_iterator begin(const this_type &v) { return v.begin(); }
358  static iterator end(this_type &v) { return v.end(); }
359  static const_iterator end(const this_type &v) { return v.end(); }
360  static origin_type* origin(this_type &v) { return v.origin; }
361  static const origin_type* origin(const this_type &v) { return v.origin; }
362  static void clear(origin_type*, const iterator &it, const iterator &ite)
363  { std::fill(it, ite, value_type(0)); }
364  static void do_clear(this_type &v)
365  { std::fill(v.begin(), v.end(), value_type(0)); }
366  static value_type access(const origin_type *, const const_iterator &it,
367  const const_iterator &, size_type i)
368  { return it[i]; }
369  static reference access(origin_type *, const iterator &it,
370  const iterator &, size_type i)
371  { return it[i]; }
372  };
373 
374  template <typename IT, typename V> std::ostream &operator <<
375  (std::ostream &o, const tab_ref_reg_spaced_with_origin<IT, V>& m)
376  { gmm::write(o,m); return o; }
377 
378 
379  template <typename IT, typename ITINDEX, typename V>
380  struct tab_ref_index_ref_with_origin
381  : public gmm::tab_ref_index_ref<IT, ITINDEX> {
382  typedef tab_ref_index_ref_with_origin<IT, ITINDEX, V> this_type;
383  typedef typename linalg_traits<this_type>::porigin_type porigin_type;
384 
385  porigin_type origin;
386 
387  tab_ref_index_ref_with_origin(void) {}
388  tab_ref_index_ref_with_origin(const IT &b, const ITINDEX &bi,
389  const ITINDEX &ei, porigin_type p)
390  : gmm::tab_ref_index_ref<IT, ITINDEX>(b, bi, ei), origin(p) {}
391 
392  tab_ref_index_ref_with_origin(const V &v, const sub_index &si)
393  : gmm::tab_ref_index_ref<IT, ITINDEX>(vect_begin(const_cast<V&>(v)),
394  si.begin(), si.end()),
395  origin(linalg_origin(const_cast<V&>(v))) {}
396  tab_ref_index_ref_with_origin(V &v, const sub_index &si)
397  : gmm::tab_ref_index_ref<IT, ITINDEX>(vect_begin(const_cast<V&>(v)),
398  si.begin(), si.end()),
399  origin(linalg_origin(const_cast<V&>(v))) {}
400  };
401 
402  template <typename IT, typename ITINDEX, typename V>
403  struct linalg_traits<tab_ref_index_ref_with_origin<IT, ITINDEX, V> > {
404  typedef typename std::iterator_traits<IT>::pointer PT;
405  typedef tab_ref_index_ref_with_origin<IT, ITINDEX, V> this_type;
406  typedef typename linalg_traits<V>::origin_type origin_type;
407  typedef typename select_ref<const origin_type *, origin_type *,
408  PT>::ref_type porigin_type;
409  typedef typename which_reference<PT>::is_reference is_reference;
410  typedef abstract_vector linalg_type;
411  typedef typename std::iterator_traits<IT>::value_type value_type;
412  typedef typename std::iterator_traits<IT>::reference reference;
413  typedef typename this_type::iterator iterator;
414  typedef typename this_type::iterator const_iterator;
415  typedef abstract_dense storage_type;
416  typedef linalg_true index_sorted;
417  static size_type size(const this_type &v) { return v.size(); }
418  static iterator begin(this_type &v) { return v.begin(); }
419  static const_iterator begin(const this_type &v) { return v.begin(); }
420  static iterator end(this_type &v) { return v.end(); }
421  static const_iterator end(const this_type &v) { return v.end(); }
422  static origin_type* origin(this_type &v) { return v.origin; }
423  static const origin_type* origin(const this_type &v) { return v.origin; }
424  static void clear(origin_type*, const iterator &it, const iterator &ite)
425  { std::fill(it, ite, value_type(0)); }
426  static void do_clear(this_type &v)
427  { std::fill(v.begin(), v.end(), value_type(0)); }
428  static value_type access(const origin_type *, const const_iterator &it,
429  const const_iterator &, size_type i)
430  { return it[i]; }
431  static reference access(origin_type *, const iterator &it,
432  const iterator &, size_type i)
433  { return it[i]; }
434  };
435 
436  template <typename IT, typename ITINDEX, typename V>
437  std::ostream &operator <<
438  (std::ostream &o, const tab_ref_index_ref_with_origin<IT, ITINDEX, V>& m)
439  { gmm::write(o,m); return o; }
440 
441 
442  template<typename ITER, typename MIT, typename PT>
443  struct dense_compressed_iterator {
444  typedef ITER value_type;
445  typedef ITER *pointer;
446  typedef ITER &reference;
447  typedef ptrdiff_t difference_type;
448  typedef std::random_access_iterator_tag iterator_category;
449  typedef size_t size_type;
450  typedef dense_compressed_iterator<ITER, MIT, PT> iterator;
451  typedef typename std::iterator_traits<PT>::value_type *MPT;
452 
453  ITER it;
454  size_type N, nrows, ncols, i;
455  PT origin;
456 
457  iterator operator ++(int) { iterator tmp = *this; i++; return tmp; }
458  iterator operator --(int) { iterator tmp = *this; i--; return tmp; }
459  iterator &operator ++() { ++i; return *this; }
460  iterator &operator --() { --i; return *this; }
461  iterator &operator +=(difference_type ii) { i += ii; return *this; }
462  iterator &operator -=(difference_type ii) { i -= ii; return *this; }
463  iterator operator +(difference_type ii) const
464  { iterator itt = *this; return (itt += ii); }
465  iterator operator -(difference_type ii) const
466  { iterator itt = *this; return (itt -= ii); }
467  difference_type operator -(const iterator &ii) const
468  { return (N ? (it - ii.it) / N : 0) + i - ii.i; }
469 
470  ITER operator *() const { return it+i*N; }
471  ITER operator [](int ii) const { return it + (i+ii) * N; }
472 
473  bool operator ==(const iterator &ii) const
474  { return (*this - ii) == difference_type(0); }
475  bool operator !=(const iterator &ii) const { return !(ii == *this); }
476  bool operator < (const iterator &ii) const
477  { return (*this - ii) < difference_type(0); }
478  bool operator > (const iterator &ii) const
479  { return (*this - ii) > difference_type(0); }
480  bool operator >=(const iterator &ii) const
481  { return (*this - ii) >= difference_type(0); }
482 
483  dense_compressed_iterator(void) {}
484  dense_compressed_iterator(const dense_compressed_iterator<MIT,MIT,MPT> &ii)
485  : it(ii.it), N(ii.N), nrows(ii.nrows), ncols(ii.ncols), i(ii.i),
486  origin(ii.origin) {}
487  dense_compressed_iterator(const ITER &iter, size_type n, size_type r,
488  size_type c, size_type ii, PT o)
489  : it(iter), N(n), nrows(r), ncols(c), i(ii), origin(o) { }
490 
491  };
492 
493  /* ******************************************************************** */
494  /* Read only reference on a compressed sparse vector */
495  /* ******************************************************************** */
496 
497  template <typename PT1, typename PT2, int shift = 0>
498  struct cs_vector_ref_iterator {
499  PT1 pr;
500  PT2 ir;
501 
502  typedef typename std::iterator_traits<PT1>::value_type value_type;
503  typedef PT1 pointer;
504  typedef typename std::iterator_traits<PT1>::reference reference;
505  typedef size_t size_type;
506  typedef ptrdiff_t difference_type;
507  typedef std::bidirectional_iterator_tag iterator_category;
508  typedef cs_vector_ref_iterator<PT1, PT2, shift> iterator;
509 
510  cs_vector_ref_iterator(void) {}
511  cs_vector_ref_iterator(PT1 p1, PT2 p2) : pr(p1), ir(p2) {}
512 
513  inline size_type index(void) const { return (*ir) - shift; }
514  iterator &operator ++() { ++pr; ++ir; return *this; }
515  iterator operator ++(int) { iterator tmp = *this; ++(*this); return tmp; }
516  iterator &operator --() { --pr; --ir; return *this; }
517  iterator operator --(int) { iterator tmp = *this; --(*this); return tmp; }
518 
519  reference operator *() const { return *pr; }
520  pointer operator ->() const { return pr; }
521 
522  bool operator ==(const iterator &i) const { return (i.pr==pr);}
523  bool operator !=(const iterator &i) const { return (i.pr!=pr);}
524  };
525 
526  template <typename PT1, typename PT2, int shift = 0> struct cs_vector_ref {
527  PT1 pr;
528  PT2 ir;
529  size_type n, size_;
530 
531  typedef cs_vector_ref<PT1, PT2, shift> this_type;
532  typedef typename std::iterator_traits<PT1>::value_type value_type;
533  typedef typename linalg_traits<this_type>::const_iterator const_iterator;
534 
535  cs_vector_ref(PT1 pt1, PT2 pt2, size_type nnz, size_type ns)
536  : pr(pt1), ir(pt2), n(nnz), size_(ns) {}
537  cs_vector_ref(void) {}
538 
539  size_type size(void) const { return size_; }
540 
541  const_iterator begin(void) const { return const_iterator(pr, ir); }
542  const_iterator end(void) const { return const_iterator(pr+n, ir+n); }
543 
544  value_type operator[](size_type i) const
545  { return linalg_traits<this_type>::access(pr, begin(), end(),i); }
546  };
547 
548  template <typename PT1, typename PT2, int shift>
549  struct linalg_traits<cs_vector_ref<PT1, PT2, shift> > {
550  typedef cs_vector_ref<PT1, PT2, shift> this_type;
551  typedef linalg_const is_reference;
552  typedef abstract_vector linalg_type;
553  typedef typename std::iterator_traits<PT1>::value_type value_type;
554  typedef value_type origin_type;
555  typedef typename std::iterator_traits<PT1>::value_type reference;
556  typedef cs_vector_ref_iterator<typename const_pointer<PT1>::pointer,
557  typename const_pointer<PT2>::pointer, shift> const_iterator;
558  typedef abstract_null_type iterator;
559  typedef abstract_sparse storage_type;
560  typedef linalg_true index_sorted;
561  static size_type size(const this_type &v) { return v.size(); }
562  static iterator begin(this_type &v) { return v.begin(); }
563  static const_iterator begin(const this_type &v) { return v.begin(); }
564  static iterator end(this_type &v) { return v.end(); }
565  static const_iterator end(const this_type &v) { return v.end(); }
566  static const origin_type* origin(const this_type &v) { return v.pr; }
567  static value_type access(const origin_type *, const const_iterator &b,
568  const const_iterator &e, size_type i) {
569  if (b.ir == e.ir) return value_type(0);
570  PT2 p = std::lower_bound(b.ir, e.ir, i+shift);
571  return (*p == i+shift && p != e.ir) ? b.pr[p-b.ir] : value_type(0);
572  }
573  };
574 
575  template <typename PT1, typename PT2, int shift>
576  std::ostream &operator <<
577  (std::ostream &o, const cs_vector_ref<PT1, PT2, shift>& m)
578  { gmm::write(o,m); return o; }
579 
580  template <typename PT1, typename PT2, int shift>
581  inline size_type nnz(const cs_vector_ref<PT1, PT2, shift>& l) { return l.n; }
582 
583  /* ******************************************************************** */
584  /* Read only reference on a compressed sparse column matrix */
585  /* ******************************************************************** */
586 
587  template <typename PT1, typename PT2, typename PT3, int shift = 0>
588  struct sparse_compressed_iterator {
589  typedef typename std::iterator_traits<PT1>::value_type value_type;
590  typedef const value_type *pointer;
591  typedef const value_type &reference;
592  typedef ptrdiff_t difference_type;
593  typedef size_t size_type;
594  typedef std::random_access_iterator_tag iterator_category;
595  typedef sparse_compressed_iterator<PT1, PT2, PT3, shift> iterator;
596 
597  PT1 pr;
598  PT2 ir;
599  PT3 jc;
600  size_type n;
601  const value_type *origin;
602 
603  iterator operator ++(int) { iterator tmp = *this; jc++; return tmp; }
604  iterator operator --(int) { iterator tmp = *this; jc--; return tmp; }
605  iterator &operator ++() { jc++; return *this; }
606  iterator &operator --() { jc--; return *this; }
607  iterator &operator +=(difference_type i) { jc += i; return *this; }
608  iterator &operator -=(difference_type i) { jc -= i; return *this; }
609  iterator operator +(difference_type i) const
610  { iterator itt = *this; return (itt += i); }
611  iterator operator -(difference_type i) const
612  { iterator itt = *this; return (itt -= i); }
613  difference_type operator -(const iterator &i) const { return jc - i.jc; }
614 
615  reference operator *() const { return pr + *jc - shift; }
616  reference operator [](int ii) { return pr + *(jc+ii) - shift; }
617 
618  bool operator ==(const iterator &i) const { return (jc == i.jc); }
619  bool operator !=(const iterator &i) const { return !(i == *this); }
620  bool operator < (const iterator &i) const { return (jc < i.jc); }
621  bool operator > (const iterator &i) const { return (jc > i.jc); }
622  bool operator >=(const iterator &i) const { return (jc >= i.jc); }
623 
624  sparse_compressed_iterator(void) {}
625  sparse_compressed_iterator(PT1 p1, PT2 p2, PT3 p3, size_type nn,
626  const value_type *o)
627  : pr(p1), ir(p2), jc(p3), n(nn), origin(o) { }
628 
629  };
630 
631  template <typename PT1, typename PT2, typename PT3, int shift = 0>
632  struct csc_matrix_ref {
633  PT1 pr; // values.
634  PT2 ir; // row indexes.
635  PT3 jc; // column repartition on pr and ir.
636  size_type nc, nr;
637 
638  typedef typename std::iterator_traits<PT1>::value_type value_type;
639  csc_matrix_ref(PT1 pt1, PT2 pt2, PT3 pt3, size_type nrr, size_type ncc)
640  : pr(pt1), ir(pt2), jc(pt3), nc(ncc), nr(nrr) {}
641  csc_matrix_ref(void) {}
642 
643  size_type nrows(void) const { return nr; }
644  size_type ncols(void) const { return nc; }
645 
646  value_type operator()(size_type i, size_type j) const
647  { return mat_col(*this, j)[i]; }
648  };
649 
650  template <typename PT1, typename PT2, typename PT3, int shift>
651  struct linalg_traits<csc_matrix_ref<PT1, PT2, PT3, shift> > {
652  typedef csc_matrix_ref<PT1, PT2, PT3, shift> this_type;
653  typedef linalg_const is_reference;
654  typedef abstract_matrix linalg_type;
655  typedef typename std::iterator_traits<PT1>::value_type value_type;
656  typedef typename std::iterator_traits<PT1>::value_type reference;
657  typedef value_type origin_type;
658  typedef abstract_sparse storage_type;
659  typedef abstract_null_type sub_row_type;
660  typedef abstract_null_type const_sub_row_type;
661  typedef abstract_null_type row_iterator;
662  typedef abstract_null_type const_row_iterator;
663  typedef abstract_null_type sub_col_type;
664  typedef cs_vector_ref<typename const_pointer<PT1>::pointer,
665  typename const_pointer<PT2>::pointer, shift> const_sub_col_type;
666  typedef sparse_compressed_iterator<typename const_pointer<PT1>::pointer,
667  typename const_pointer<PT2>::pointer,
668  typename const_pointer<PT3>::pointer,
669  shift> const_col_iterator;
670  typedef abstract_null_type col_iterator;
671  typedef col_major sub_orientation;
672  typedef linalg_true index_sorted;
673  static size_type nrows(const this_type &m) { return m.nrows(); }
674  static size_type ncols(const this_type &m) { return m.ncols(); }
675  static const_col_iterator col_begin(const this_type &m)
676  { return const_col_iterator(m.pr, m.ir, m.jc, m.nr, m.pr); }
677  static const_col_iterator col_end(const this_type &m)
678  { return const_col_iterator(m.pr, m.ir, m.jc + m.nc, m.nr, m.pr); }
679  static const_sub_col_type col(const const_col_iterator &it) {
680  return const_sub_col_type(it.pr + *(it.jc) - shift,
681  it.ir + *(it.jc) - shift, *(it.jc + 1) - *(it.jc), it.n);
682  }
683  static const origin_type* origin(const this_type &m) { return m.pr; }
684  static value_type access(const const_col_iterator &itcol, size_type j)
685  { return col(itcol)[j]; }
686  };
687 
688 
689  template <typename PT1, typename PT2, typename PT3, int shift>
690  std::ostream &operator <<
691  (std::ostream &o, const csc_matrix_ref<PT1, PT2, PT3, shift>& m)
692  { gmm::write(o,m); return o; }
693 
694  /* ******************************************************************** */
695  /* Read only reference on a compressed sparse row matrix */
696  /* ******************************************************************** */
697 
698  template <typename PT1, typename PT2, typename PT3, int shift = 0>
699  struct csr_matrix_ref {
700  PT1 pr; // values.
701  PT2 ir; // column indexes.
702  PT3 jc; // row repartition on pr and ir.
703  size_type nc, nr;
704 
705  typedef typename std::iterator_traits<PT1>::value_type value_type;
706  csr_matrix_ref(PT1 pt1, PT2 pt2, PT3 pt3, size_type nrr, size_type ncc)
707  : pr(pt1), ir(pt2), jc(pt3), nc(ncc), nr(nrr) {}
708  csr_matrix_ref(void) {}
709 
710  size_type nrows(void) const { return nr; }
711  size_type ncols(void) const { return nc; }
712 
713  value_type operator()(size_type i, size_type j) const
714  { return mat_row(*this, i)[j]; }
715  };
716 
717  template <typename PT1, typename PT2, typename PT3, int shift>
718  struct linalg_traits<csr_matrix_ref<PT1, PT2, PT3, shift> > {
719  typedef csr_matrix_ref<PT1, PT2, PT3, shift> this_type;
720  typedef linalg_const is_reference;
721  typedef abstract_matrix linalg_type;
722  typedef typename std::iterator_traits<PT1>::value_type value_type;
723  typedef typename std::iterator_traits<PT1>::value_type reference;
724  typedef value_type origin_type;
725  typedef abstract_sparse storage_type;
726  typedef abstract_null_type sub_col_type;
727  typedef abstract_null_type const_sub_col_type;
728  typedef abstract_null_type col_iterator;
729  typedef abstract_null_type const_col_iterator;
730  typedef abstract_null_type sub_row_type;
731  typedef cs_vector_ref<typename const_pointer<PT1>::pointer,
732  typename const_pointer<PT2>::pointer, shift>
733  const_sub_row_type;
734  typedef sparse_compressed_iterator<typename const_pointer<PT1>::pointer,
735  typename const_pointer<PT2>::pointer,
736  typename const_pointer<PT3>::pointer,
737  shift> const_row_iterator;
738  typedef abstract_null_type row_iterator;
739  typedef row_major sub_orientation;
740  typedef linalg_true index_sorted;
741  static size_type nrows(const this_type &m) { return m.nrows(); }
742  static size_type ncols(const this_type &m) { return m.ncols(); }
743  static const_row_iterator row_begin(const this_type &m)
744  { return const_row_iterator(m.pr, m.ir, m.jc, m.nc, m.pr); }
745  static const_row_iterator row_end(const this_type &m)
746  { return const_row_iterator(m.pr, m.ir, m.jc + m.nr, m.nc, m.pr); }
747  static const_sub_row_type row(const const_row_iterator &it) {
748  return const_sub_row_type(it.pr + *(it.jc) - shift,
749  it.ir + *(it.jc) - shift, *(it.jc + 1) - *(it.jc), it.n);
750  }
751  static const origin_type* origin(const this_type &m) { return m.pr; }
752  static value_type access(const const_row_iterator &itrow, size_type j)
753  { return row(itrow)[j]; }
754  };
755 
756  template <typename PT1, typename PT2, typename PT3, int shift>
757  std::ostream &operator <<
758  (std::ostream &o, const csr_matrix_ref<PT1, PT2, PT3, shift>& m)
759  { gmm::write(o,m); return o; }
760 
761  /* ********************************************************************* */
762  /* */
763  /* Simple interface for C arrays */
764  /* */
765  /* ********************************************************************* */
766 
767  template <class PT> struct array1D_reference {
768 
769  typedef typename std::iterator_traits<PT>::value_type value_type;
770 
771  PT begin, end;
772 
773  const value_type &operator[](size_type i) const { return *(begin+i); }
774  value_type &operator[](size_type i) { return *(begin+i); }
775 
776  array1D_reference(PT begin_, size_type s) : begin(begin_), end(begin_+s) {}
777  };
778 
779  template <typename PT>
780  struct linalg_traits<array1D_reference<PT> > {
781  typedef array1D_reference<PT> this_type;
782  typedef this_type origin_type;
783  typedef typename which_reference<PT>::is_reference is_reference;
784  typedef abstract_vector linalg_type;
785  typedef typename std::iterator_traits<PT>::value_type value_type;
786  typedef typename std::iterator_traits<PT>::reference reference;
787  typedef PT iterator;
788  typedef PT const_iterator;
789  typedef abstract_dense storage_type;
790  typedef linalg_true index_sorted;
791  static size_type size(const this_type &v) { return v.end - v.begin; }
792  static iterator begin(this_type &v) { return v.begin; }
793  static const_iterator begin(const this_type &v) { return v.begin; }
794  static iterator end(this_type &v) { return v.end; }
795  static const_iterator end(const this_type &v) { return v.end; }
796  static origin_type* origin(this_type &v) { return &v; }
797  static const origin_type* origin(const this_type &v) { return &v; }
798  static void clear(origin_type*, const iterator &it, const iterator &ite)
799  { std::fill(it, ite, value_type(0)); }
800  static void do_clear(this_type &v)
801  { std::fill(v.begin, v.end, value_type(0)); }
802  static value_type access(const origin_type *, const const_iterator &it,
803  const const_iterator &, size_type i)
804  { return it[i]; }
805  static reference access(origin_type *, const iterator &it,
806  const iterator &, size_type i)
807  { return it[i]; }
808  static void resize(this_type &, size_type )
809  { GMM_ASSERT1(false, "Not resizable vector"); }
810  };
811 
812  template<typename PT> std::ostream &operator <<
813  (std::ostream &o, const array1D_reference<PT>& v)
814  { gmm::write(o,v); return o; }
815 
816  template <class PT> struct array2D_col_reference {
817 
818  typedef typename std::iterator_traits<PT>::value_type T;
819  typedef typename std::iterator_traits<PT>::reference reference;
820  typedef typename const_reference<reference>::reference const_reference;
821  typedef PT iterator;
822  typedef typename const_pointer<PT>::pointer const_iterator;
823 
824  PT begin_;
825  size_type nbl, nbc;
826 
827  inline const_reference operator ()(size_type l, size_type c) const {
828  GMM_ASSERT2(l < nbl && c < nbc, "out of range");
829  return *(begin_ + c*nbl+l);
830  }
831  inline reference operator ()(size_type l, size_type c) {
832  GMM_ASSERT2(l < nbl && c < nbc, "out of range");
833  return *(begin_ + c*nbl+l);
834  }
835 
836  void resize(size_type, size_type);
837  void reshape(size_type m, size_type n) {
838  GMM_ASSERT2(n*m == nbl*nbc, "dimensions mismatch");
839  nbl = m; nbc = n;
840  }
841 
842  void fill(T a, T b = T(0)) {
843  std::fill(begin_, begin_+nbc*nbl, b);
844  iterator p = begin_, e = begin_+nbc*nbl;
845  while (p < e) { *p = a; p += nbl+1; }
846  }
847  inline size_type nrows(void) const { return nbl; }
848  inline size_type ncols(void) const { return nbc; }
849 
850  iterator begin(void) { return begin_; }
851  const_iterator begin(void) const { return begin_; }
852  iterator end(void) { return begin_+nbl*nbc; }
853  const_iterator end(void) const { return begin_+nbl*nbc; }
854 
855  array2D_col_reference(PT begin__, size_type nrows_, size_type ncols_)
856  : begin_(begin__), nbl(nrows_), nbc(ncols_) {}
857  };
858 
859  template <typename PT> struct linalg_traits<array2D_col_reference<PT> > {
860  typedef array2D_col_reference<PT> this_type;
861  typedef this_type origin_type;
862  typedef typename which_reference<PT>::is_reference is_reference;
863  typedef abstract_matrix linalg_type;
864  typedef typename std::iterator_traits<PT>::value_type value_type;
865  typedef typename std::iterator_traits<PT>::reference reference;
866  typedef abstract_dense storage_type;
867  typedef tab_ref_reg_spaced_with_origin<typename this_type::iterator,
868  this_type> sub_row_type;
869  typedef tab_ref_reg_spaced_with_origin<typename this_type::const_iterator,
870  this_type> const_sub_row_type;
871  typedef dense_compressed_iterator<typename this_type::iterator,
872  typename this_type::iterator,
873  this_type *> row_iterator;
874  typedef dense_compressed_iterator<typename this_type::const_iterator,
875  typename this_type::iterator,
876  const this_type *> const_row_iterator;
877  typedef tab_ref_with_origin<typename this_type::iterator,
878  this_type> sub_col_type;
879  typedef tab_ref_with_origin<typename this_type::const_iterator,
880  this_type> const_sub_col_type;
881  typedef dense_compressed_iterator<typename this_type::iterator,
882  typename this_type::iterator,
883  this_type *> col_iterator;
884  typedef dense_compressed_iterator<typename this_type::const_iterator,
885  typename this_type::iterator,
886  const this_type *> const_col_iterator;
887  typedef col_and_row sub_orientation;
888  typedef linalg_true index_sorted;
889  static size_type nrows(const this_type &m) { return m.nrows(); }
890  static size_type ncols(const this_type &m) { return m.ncols(); }
891  static const_sub_row_type row(const const_row_iterator &it)
892  { return const_sub_row_type(*it, it.nrows, it.ncols, it.origin); }
893  static const_sub_col_type col(const const_col_iterator &it)
894  { return const_sub_col_type(*it, *it + it.nrows, it.origin); }
895  static sub_row_type row(const row_iterator &it)
896  { return sub_row_type(*it, it.nrows, it.ncols, it.origin); }
897  static sub_col_type col(const col_iterator &it)
898  { return sub_col_type(*it, *it + it.nrows, it.origin); }
899  static row_iterator row_begin(this_type &m)
900  { return row_iterator(m.begin(), 1, m.nrows(), m.ncols(), 0, &m); }
901  static row_iterator row_end(this_type &m)
902  { return row_iterator(m.begin(), 1, m.nrows(), m.ncols(), m.nrows(), &m); }
903  static const_row_iterator row_begin(const this_type &m)
904  { return const_row_iterator(m.begin(), 1, m.nrows(), m.ncols(), 0, &m); }
905  static const_row_iterator row_end(const this_type &m) {
906  return const_row_iterator(m.begin(), 1, m.nrows(),
907  m.ncols(), m.nrows(), &m);
908  }
909  static col_iterator col_begin(this_type &m)
910  { return col_iterator(m.begin(), m.nrows(), m.nrows(), m.ncols(), 0, &m); }
911  static col_iterator col_end(this_type &m) {
912  return col_iterator(m.begin(), m.nrows(), m.nrows(), m.ncols(),
913  m.ncols(), &m);
914  }
915  static const_col_iterator col_begin(const this_type &m) {
916  return const_col_iterator(m.begin(), m.nrows(), m.nrows(),
917  m.ncols(), 0, &m);
918  }
919  static const_col_iterator col_end(const this_type &m) {
920  return const_col_iterator(m.begin(), m.nrows(),m.nrows(),m.ncols(),
921  m.ncols(), &m);
922  }
923  static origin_type* origin(this_type &m) { return &m; }
924  static const origin_type* origin(const this_type &m) { return &m; }
925  static void do_clear(this_type &m) { m.fill(value_type(0)); }
926  static value_type access(const const_col_iterator &itcol, size_type j)
927  { return (*itcol)[j]; }
928  static reference access(const col_iterator &itcol, size_type j)
929  { return (*itcol)[j]; }
930  static void resize(this_type &v, size_type m, size_type n)
931  { v.resize(m,n); }
932  static void reshape(this_type &v, size_type m, size_type n)
933  { v.reshape(m, n); }
934  };
935 
936  template<typename PT> std::ostream &operator <<
937  (std::ostream &o, const array2D_col_reference<PT>& m)
938  { gmm::write(o,m); return o; }
939 
940 
941 
942  template <class PT> struct array2D_row_reference {
943 
944  typedef typename std::iterator_traits<PT>::value_type T;
945  typedef typename std::iterator_traits<PT>::reference reference;
946  typedef typename const_reference<reference>::reference const_reference;
947  typedef PT iterator;
948  typedef typename const_pointer<PT>::pointer const_iterator;
949 
950  PT begin_;
951  size_type nbl, nbc;
952 
953  inline const_reference operator ()(size_type l, size_type c) const {
954  GMM_ASSERT2(l < nbl && c < nbc, "out of range");
955  return *(begin_ + l*nbc+c);
956  }
957  inline reference operator ()(size_type l, size_type c) {
958  GMM_ASSERT2(l < nbl && c < nbc, "out of range");
959  return *(begin_ + l*nbc+c);
960  }
961 
962  void resize(size_type, size_type);
963  void reshape(size_type m, size_type n) {
964  GMM_ASSERT2(n*m == nbl*nbc, "dimensions mismatch");
965  nbl = m; nbc = n;
966  }
967 
968  void fill(T a, T b = T(0)) {
969  std::fill(begin_, begin_+nbc*nbl, b);
970  iterator p = begin_, e = begin_+nbc*nbl;
971  while (p < e) { *p = a; p += nbc+1; }
972  }
973  inline size_type nrows(void) const { return nbl; }
974  inline size_type ncols(void) const { return nbc; }
975 
976  iterator begin(void) { return begin_; }
977  const_iterator begin(void) const { return begin_; }
978  iterator end(void) { return begin_+nbl*nbc; }
979  const_iterator end(void) const { return begin_+nbl*nbc; }
980 
981  array2D_row_reference(PT begin__, size_type nrows_, size_type ncols_)
982  : begin_(begin__), nbl(nrows_), nbc(ncols_) {}
983  };
984 
985  template <typename PT> struct linalg_traits<array2D_row_reference<PT> > {
986  typedef array2D_row_reference<PT> this_type;
987  typedef this_type origin_type;
988  typedef typename which_reference<PT>::is_reference is_reference;
989  typedef abstract_matrix linalg_type;
990  typedef typename std::iterator_traits<PT>::value_type value_type;
991  typedef typename std::iterator_traits<PT>::reference reference;
992  typedef abstract_dense storage_type;
993  typedef tab_ref_reg_spaced_with_origin<typename this_type::iterator,
994  this_type> sub_col_type;
995  typedef tab_ref_reg_spaced_with_origin<typename this_type::const_iterator,
996  this_type> const_sub_col_type;
997  typedef dense_compressed_iterator<typename this_type::iterator,
998  typename this_type::iterator,
999  this_type *> col_iterator;
1000  typedef dense_compressed_iterator<typename this_type::const_iterator,
1001  typename this_type::iterator,
1002  const this_type *> const_col_iterator;
1003  typedef tab_ref_with_origin<typename this_type::iterator,
1004  this_type> sub_row_type;
1005  typedef tab_ref_with_origin<typename this_type::const_iterator,
1006  this_type> const_sub_row_type;
1007  typedef dense_compressed_iterator<typename this_type::iterator,
1008  typename this_type::iterator,
1009  this_type *> row_iterator;
1010  typedef dense_compressed_iterator<typename this_type::const_iterator,
1011  typename this_type::iterator,
1012  const this_type *> const_row_iterator;
1013  typedef col_and_row sub_orientation;
1014  typedef linalg_true index_sorted;
1015  static size_type ncols(const this_type &m) { return m.ncols(); }
1016  static size_type nrows(const this_type &m) { return m.nrows(); }
1017  static const_sub_col_type col(const const_col_iterator &it)
1018  { return const_sub_col_type(*it, it.ncols, it.nrows, it.origin); }
1019  static const_sub_row_type row(const const_row_iterator &it)
1020  { return const_sub_row_type(*it, *it + it.ncols, it.origin); }
1021  static sub_col_type col(const col_iterator &it)
1022  { return sub_col_type(*it, *it, it.ncols, it.nrows, it.origin); }
1023  static sub_row_type row(const row_iterator &it)
1024  { return sub_row_type(*it, *it + it.ncols, it.origin); }
1025  static col_iterator col_begin(this_type &m)
1026  { return col_iterator(m.begin(), 1, m.ncols(), m.nrows(), 0, &m); }
1027  static col_iterator col_end(this_type &m)
1028  { return col_iterator(m.begin(), 1, m.ncols(), m.nrows(), m.ncols(), &m); }
1029  static const_col_iterator col_begin(const this_type &m)
1030  { return const_col_iterator(m.begin(), 1, m.ncols(), m.nrows(), 0, &m); }
1031  static const_col_iterator col_end(const this_type &m) {
1032  return const_col_iterator(m.begin(), 1, m.ncols(),
1033  m.nrows(), m.ncols(), &m);
1034  }
1035  static row_iterator row_begin(this_type &m)
1036  { return row_iterator(m.begin(), m.ncols(), m.ncols(), m.nrows(), 0, &m); }
1037  static row_iterator row_end(this_type &m) {
1038  return row_iterator(m.begin(), m.ncols(), m.ncols(), m.nrows(),
1039  m.nrows(), &m);
1040  }
1041  static const_row_iterator row_begin(const this_type &m) {
1042  return const_row_iterator(m.begin(), m.ncols(), m.ncols(), m.nrows(),
1043  0, &m);
1044  }
1045  static const_row_iterator row_end(const this_type &m) {
1046  return const_row_iterator(m.begin(), m.ncols(), m.ncols(), m.nrows(),
1047  m.nrows(), &m);
1048  }
1049  static origin_type* origin(this_type &m) { return &m; }
1050  static const origin_type* origin(const this_type &m) { return &m; }
1051  static void do_clear(this_type &m) { m.fill(value_type(0)); }
1052  static value_type access(const const_row_iterator &itrow, size_type j)
1053  { return (*itrow)[j]; }
1054  static reference access(const row_iterator &itrow, size_type j)
1055  { return (*itrow)[j]; }
1056  static void resize(this_type &v, size_type m, size_type n)
1057  { v.resize(m,n); }
1058  static void reshape(this_type &v, size_type m, size_type n)
1059  { v.reshape(m, n); }
1060  };
1061 
1062  template<typename PT> std::ostream &operator <<
1063  (std::ostream &o, const array2D_row_reference<PT>& m)
1064  { gmm::write(o,m); return o; }
1065 
1066 
1067 
1068 
1069 
1070 
1071 }
1072 
1073 
1074 #endif // GMM_INTERFACE_H__
indexed array reference (given a container X, and a set of indexes I, this class provides a pseudo-co...
Definition: gmm_ref.h:304
provide a "strided" view a of container
Definition: gmm_ref.h:426
Basic linear algebra functions.
size_type nnz(const L &l)
count the number of non-zero entries of a vector or matrix.
Definition: gmm_blas.h:69
void reshape(M &v, size_type m, size_type n)
*‍/
Definition: gmm_blas.h:251
void fill(L &l, typename gmm::linalg_traits< L >::value_type x)
*‍/
Definition: gmm_blas.h:104
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
void resize(V &v, size_type n)
*‍/
Definition: gmm_blas.h:210
sub-indices.
rational_fraction< T > operator-(const polynomial< T > &P, const rational_fraction< T > &Q)
Subtract Q from P.
Definition: bgeot_poly.h:756
rational_fraction< T > operator+(const polynomial< T > &P, const rational_fraction< T > &Q)
Add Q to P.
Definition: bgeot_poly.h:749
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49