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