39 #ifndef GMM_INTERFACE_H__ 
   40 #define GMM_INTERFACE_H__ 
   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;
 
   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;
 
   95     iterator begin_, end_;
 
   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)) {}
 
  104     simple_vector_ref(
const simple_vector_ref<CPT> &cr)
 
  105       : begin_(cr.begin_),end_(cr.end_),origin(cr.origin),size_(cr.size_) {}
 
  107     simple_vector_ref(
void) {}
 
  110     { 
return linalg_traits<V>::access(origin, begin_, end_, i); }
 
  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());
 
  119   template <
typename IT, 
typename ORG, 
typename PT> 
inline 
  120   void set_to_begin(IT &it, ORG o, 
const simple_vector_ref<PT> *,
 
  122     typedef typename linalg_traits<simple_vector_ref<PT> >::V_reference ref_t;
 
  123     set_to_begin(it, o, PT(), ref_t());
 
  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());
 
  132   template <
typename IT, 
typename ORG, 
typename PT> 
inline 
  133   void set_to_end(IT &it, ORG o, 
const simple_vector_ref<PT> *,
 
  135     typedef typename linalg_traits<simple_vector_ref<PT> >::V_reference ref_t;
 
  136     set_to_end(it, o, PT(), ref_t());
 
  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;
 
  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()); 
 
  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());
 
  170     static inline iterator end(this_type &v) {
 
  171       iterator it = v.end_;
 
  172       set_to_end(it, v.origin, pthis_type(), is_reference());
 
  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());
 
  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,
 
  187     { 
return linalg_traits<V>::access(o, it, ite, i); }
 
  188     static reference access(origin_type *o, 
const iterator &it,
 
  190     { 
return linalg_traits<V>::access(o, it, ite, i); }
 
  193   template <
typename PT>
 
  194   std::ostream &operator << (std::ostream &o, 
const simple_vector_ref<PT>& v)
 
  195   { gmm::write(o,v); 
return o; }
 
  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); }
 
  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,
 
  234     static reference access(origin_type *, 
const iterator &it,
 
  242   template <
typename T>
 
  243   inline size_type nnz(
const std::vector<T>& l) { 
return l.size(); }
 
  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;
 
  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;
 
  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) {}
 
  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))) {}
 
  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,
 
  309     static reference access(origin_type *, 
const iterator &it, 
 
  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; }
 
  319   template <
typename IT, 
typename V>
 
  321     typedef  tab_ref_reg_spaced_with_origin<IT, V> this_type;
 
  322     typedef typename linalg_traits<this_type>::porigin_type porigin_type;
 
  326     tab_ref_reg_spaced_with_origin(
void) {}
 
  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))) {}
 
  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,
 
  369     static reference access(origin_type *, 
const iterator &it, 
 
  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; }
 
  379   template <
typename IT, 
typename ITINDEX, 
typename V>
 
  380   struct tab_ref_index_ref_with_origin 
 
  382     typedef tab_ref_index_ref_with_origin<IT, ITINDEX, V> this_type;
 
  383     typedef typename linalg_traits<this_type>::porigin_type porigin_type;
 
  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) {}
 
  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))) {}
 
  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,
 
  431     static reference access(origin_type *, 
const iterator &it,
 
  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; }
 
  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;
 
  450     typedef dense_compressed_iterator<ITER, MIT, PT> iterator;
 
  451     typedef typename std::iterator_traits<PT>::value_type *MPT;
 
  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; }
 
  464     { iterator itt = *
this; 
return (itt += ii); }
 
  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; }
 
  470     ITER operator *()
 const { 
return it+i*N; }
 
  471     ITER operator [](
int ii)
 const { 
return it + (i+ii) * N; }
 
  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); }
 
  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),
 
  489       : it(iter), N(n), nrows(r), ncols(c), i(ii), origin(o) { }
 
  497   template <
typename PT1, 
typename PT2, 
int shift = 0>
 
  498   struct cs_vector_ref_iterator {
 
  502     typedef typename std::iterator_traits<PT1>::value_type value_type;
 
  504     typedef typename std::iterator_traits<PT1>::reference  reference;
 
  506     typedef ptrdiff_t     difference_type;
 
  507     typedef std::bidirectional_iterator_tag iterator_category;
 
  508     typedef cs_vector_ref_iterator<PT1, PT2, shift> iterator;
 
  510     cs_vector_ref_iterator(
void) {}
 
  511     cs_vector_ref_iterator(PT1 p1, PT2 p2) : pr(p1), ir(p2) {}
 
  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; }
 
  519     reference operator  *()
 const { 
return *pr; }
 
  520     pointer   operator ->()
 const { 
return pr; }
 
  522     bool operator ==(
const iterator &i)
 const { 
return (i.pr==pr);}
 
  523     bool operator !=(
const iterator &i)
 const { 
return (i.pr!=pr);}
 
  526   template <
typename PT1, 
typename PT2, 
int shift = 0> 
struct cs_vector_ref {
 
  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;
 
  536       : pr(pt1), ir(pt2), n(
nnz), size_(ns) {}
 
  537     cs_vector_ref(
void) {}
 
  539     size_type size(
void)
 const { 
return size_; }
 
  541     const_iterator begin(
void)
 const { 
return const_iterator(pr, ir); }
 
  542     const_iterator end(
void)
 const { 
return const_iterator(pr+n, ir+n); }
 
  545     { 
return linalg_traits<this_type>::access(pr, begin(), end(),i); }
 
  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,
 
  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);
 
  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; }
 
  580   template <
typename PT1, 
typename PT2, 
int shift>
 
  581   inline size_type nnz(
const cs_vector_ref<PT1, PT2, shift>& l) { 
return l.n; }
 
  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;
 
  594     typedef std::random_access_iterator_tag iterator_category;
 
  595     typedef sparse_compressed_iterator<PT1, PT2, PT3, shift> iterator;
 
  601     const value_type *origin;
 
  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; }
 
  610     { iterator itt = *
this; 
return (itt += i); }
 
  612     { iterator itt = *
this; 
return (itt -= i); }
 
  613     difference_type 
operator -(
const iterator &i)
 const { 
return jc - i.jc; }
 
  615     reference operator *()
 const { 
return pr + *jc - shift; }
 
  616     reference operator [](
int ii) { 
return pr + *(jc+ii) - shift; }
 
  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); }
 
  624     sparse_compressed_iterator(
void) {}
 
  625     sparse_compressed_iterator(PT1 p1, PT2 p2, PT3 p3, 
size_type nn,
 
  627       : pr(p1), ir(p2), jc(p3), n(nn), origin(o) { }
 
  631   template <
typename PT1, 
typename PT2, 
typename PT3, 
int shift = 0>
 
  632   struct csc_matrix_ref {
 
  638     typedef typename std::iterator_traits<PT1>::value_type value_type;
 
  640       : pr(pt1), ir(pt2), jc(pt3), nc(ncc), nr(nrr) {}
 
  641     csc_matrix_ref(
void) {}
 
  643     size_type nrows(
void)
 const { 
return nr; }
 
  644     size_type ncols(
void)
 const { 
return nc; }
 
  647       { 
return mat_col(*
this, j)[i]; }
 
  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);
 
  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]; }
 
  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; }
 
  698   template <
typename PT1, 
typename PT2, 
typename PT3, 
int shift = 0>
 
  699   struct csr_matrix_ref {
 
  705     typedef typename std::iterator_traits<PT1>::value_type value_type;
 
  707       : pr(pt1), ir(pt2), jc(pt3), nc(ncc), nr(nrr) {}
 
  708     csr_matrix_ref(
void) {}
 
  710     size_type nrows(
void)
 const { 
return nr; }
 
  711     size_type ncols(
void)
 const { 
return nc; }
 
  714       { 
return mat_row(*
this, i)[j]; }
 
  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>
 
  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);
 
  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]; }
 
  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; }
 
  767   template <
class PT> 
struct array1D_reference {
 
  769     typedef typename std::iterator_traits<PT>::value_type value_type;
 
  773     const value_type &operator[](
size_type i)
 const { 
return *(begin+i); }
 
  774     value_type &operator[](
size_type i) { 
return *(begin+i); }
 
  776     array1D_reference(PT begin_, 
size_type s) : begin(begin_), end(begin_+s) {}
 
  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;
 
  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,
 
  805     static reference access(origin_type *, 
const iterator &it,
 
  809     { GMM_ASSERT1(
false, 
"Not resizable vector"); }
 
  812   template<
typename PT> std::ostream &
operator <<
 
  813   (std::ostream &o, 
const array1D_reference<PT>& v)
 
  814   { gmm::write(o,v); 
return o; }
 
  816   template <
class PT> 
struct array2D_col_reference {
 
  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;
 
  822     typedef typename const_pointer<PT>::pointer const_iterator;
 
  828       GMM_ASSERT2(l < nbl && c < nbc, 
"out of range");
 
  829       return *(begin_ + c*nbl+l);
 
  832       GMM_ASSERT2(l < nbl && c < nbc, 
"out of range");
 
  833       return *(begin_ + c*nbl+l);
 
  838       GMM_ASSERT2(n*m == nbl*nbc, 
"dimensions mismatch");
 
  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; }
 
  847     inline size_type nrows(
void)
 const { 
return nbl; }
 
  848     inline size_type ncols(
void)
 const { 
return nbc; }
 
  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; }
 
  856       : begin_(begin__), nbl(nrows_), nbc(ncols_) {}
 
  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);
 
  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(),
 
  915     static const_col_iterator col_begin(
const this_type &m) {
 
  916       return const_col_iterator(m.begin(), m.nrows(), m.nrows(),
 
  919     static const_col_iterator col_end(
const this_type &m) {
 
  920       return const_col_iterator(m.begin(), m.nrows(),m.nrows(),m.ncols(),
 
  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]; }
 
  936   template<
typename PT> std::ostream &
operator <<
 
  937     (std::ostream &o, 
const array2D_col_reference<PT>& m)
 
  938   { gmm::write(o,m); 
return o; }
 
  942   template <
class PT> 
struct array2D_row_reference {
 
  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;
 
  948     typedef typename const_pointer<PT>::pointer const_iterator;
 
  954       GMM_ASSERT2(l < nbl && c < nbc, 
"out of range");
 
  955       return *(begin_ + l*nbc+c);
 
  958       GMM_ASSERT2(l < nbl && c < nbc, 
"out of range");
 
  959       return *(begin_ + l*nbc+c);
 
  964       GMM_ASSERT2(n*m == nbl*nbc, 
"dimensions mismatch");
 
  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; }
 
  973     inline size_type nrows(
void)
 const { 
return nbl; }
 
  974     inline size_type ncols(
void)
 const { 
return nbc; }
 
  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; }
 
  982       : begin_(begin__), nbl(nrows_), nbc(ncols_) {}
 
  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);
 
 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(),
 
 1041     static const_row_iterator row_begin(
const this_type &m) {
 
 1042       return const_row_iterator(m.begin(), m.ncols(), m.ncols(), m.nrows(),
 
 1045     static const_row_iterator row_end(
const this_type &m) {
 
 1046       return const_row_iterator(m.begin(), m.ncols(), m.ncols(), m.nrows(),
 
 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]; }
 
 1059     { v.reshape(m, n); }
 
 1062   template<
typename PT> std::ostream &
operator <<
 
 1063     (std::ostream &o, 
const array2D_row_reference<PT>& m)
 
 1064   { gmm::write(o,m); 
return o; }
 
indexed array reference (given a container X, and a set of indexes I, this class provides a pseudo-co...
 
provide a "strided" view a of container
 
Basic linear algebra functions.
 
size_type nnz(const L &l)
count the number of non-zero entries of a vector or matrix.
 
void reshape(M &v, size_type m, size_type n)
*/
 
void fill(L &l, typename gmm::linalg_traits< L >::value_type x)
*/
 
void clear(L &l)
clear (fill with zeros) a vector or matrix.
 
void resize(V &v, size_type n)
*/
 
rational_fraction< T > operator-(const polynomial< T > &P, const rational_fraction< T > &Q)
Subtract Q from P.
 
rational_fraction< T > operator+(const polynomial< T > &P, const rational_fraction< T > &Q)
Add Q to P.
 
size_t size_type
used as the common size type in the library