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