matrix_proxy.hpp
Go to the documentation of this file.
1 /*!
2  *
3  *
4  * \brief Matrix proxy classes.
5  *
6  *
7  *
8  * \author O. Krause
9  * \date 2013
10  *
11  *
12  * \par Copyright 1995-2015 Shark Development Team
13  *
14  * <BR><HR>
15  * This file is part of Shark.
16  * <http://image.diku.dk/shark/>
17  *
18  * Shark is free software: you can redistribute it and/or modify
19  * it under the terms of the GNU Lesser General Public License as published
20  * by the Free Software Foundation, either version 3 of the License, or
21  * (at your option) any later version.
22  *
23  * Shark is distributed in the hope that it will be useful,
24  * but WITHOUT ANY WARRANTY; without even the implied warranty of
25  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26  * GNU Lesser General Public License for more details.
27  *
28  * You should have received a copy of the GNU Lesser General Public License
29  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
30  *
31  */
32 
33 #ifndef SHARK_LINALG_BLAS_MATRIX_PROXY_HPP
34 #define SHARK_LINALG_BLAS_MATRIX_PROXY_HPP
35 
36 #include "assignment.hpp"
37 
38 namespace shark {
39 namespace blas {
40 
41 ///\brief Wraps another expression as a reference.
42 template<class M>
43 class matrix_reference:public matrix_expression<matrix_reference<M> > {
44 public:
45  typedef typename M::size_type size_type;
46  typedef typename M::difference_type difference_type;
47  typedef typename M::value_type value_type;
48  typedef typename M::scalar_type scalar_type;
49  typedef typename M::const_reference const_reference;
50  typedef typename reference<M>::type reference;
51  typedef typename M::const_pointer const_pointer;
52  typedef typename pointer<M>::type pointer;
53 
54  typedef typename M::index_type index_type;
55  typedef typename M::const_index_pointer const_index_pointer;
57 
60  typedef typename M::orientation orientation;
61  typedef typename M::storage_category storage_category;
62  typedef elementwise_tag evaluation_category;
63 
64 
65  // Construction and destruction
66  matrix_reference(M& m):m_expression(&m) {}
67  template<class E>
69  :m_expression(&other.expression()){}
70 
71  // Accessors
72 
73  M& expression() const {
74  return *m_expression;
75  }
76  M& expression() {
77  return *m_expression;
78  }
79 
80  ///\brief Returns the number of rows of the matrix.
81  size_type size1() const {
82  return expression().size1();
83  }
84  ///\brief Returns the number of columns of the matrix.
85  size_type size2() const {
86  return expression().size2();
87  }
88 
89  // ---------
90  // Dense low level interface
91  // ---------
92 
93  ///\brief Returns the stride in memory between two rows.
94  difference_type stride1()const{
95  return expression().stride1();
96  }
97  ///\brief Returns the stride in memory between two columns.
98  difference_type stride2()const{
99  return expression().stride2();
100  }
101 
102  ///\brief Returns the pointer to the beginning of the matrix storage.
103  ///
104  /// Grants low-level access to the matrix internals. Element order depends on whether the matrix is row_major or column_major.
105  /// to access element (i,j) use storage()[i*stride1()+j*stride2()].
106  pointer storage()const{
107  return expression().storage();
108  }
109 
110  // ---------
111  // Sparse low level interface
112  // ---------
113 
114  /// \brief Number of nonzero elements of the matrix.
115  size_type nnz()const{
116  return expression().nnz();
117  }
118  /// \brief Array of values of the nonzero elements.
119  const_pointer values()const{
120  return expression().values();
121  }
122 
123  /// \brief Array of indices of the nonzero elements.
124  ///
125  /// Note that there is a pair of indices needed:
126  /// When accessing the j-th element in row i you have to write code like this:
127  /// index_type start = outer_indices()[i] //aquire start of the i-th row
128  /// index = inner_indices()[start+j];
129  /// All elements in the row are contained in the range [outer_indices()[i],outer_indices_end()[i])
130  /// there might be gaps between the end of the one line and the start of the next!
131  index_pointer inner_indices()const{
132  return expression().inner_indices();
133  }
134 
135  ///\brief Returns an array containing the start of the rows.
136  ///
137  /// See documentation of inner_indices() for more details
138  index_pointer outer_indices()const{
139  return expression().outer_indices();
140  }
141 
142  ///\brief Returns an array containing the end of the rows.
143  ///
144  /// See documentation of inner_indices() for more details
145  index_pointer outer_indices_end()const{
146  return expression().outer_indices_end();
147  }
148 
149  /// \brief Returns the number of nonzero elements in the i-th row/column.
150  size_type inner_nnz(index_type i) const {
151  return expression().inner_nnz(i);
152  }
153 
154  // ---------
155  // High level interface
156  // ---------
157 
158  // Element access
159  reference operator()(index_type i, index_type j) const {
160  return expression()(i, j);
161  }
162 
163  void set_element(size_type i, size_type j,value_type t){
164  expression().set_element(i,j,t);
165  }
166 
167 
168  // Assignment
169  template<class E>
171  expression() = e();
172  return *this;
173  }
174 
175  // Iterator types
177  typedef row_iterator const_row_iterator;
179  typedef column_iterator const_column_iterator;
180 
181  // Iterators are the iterators of the referenced expression.
182  const_row_iterator row_begin(index_type i) const {
183  return expression().row_begin(i);
184  }
185  const_row_iterator row_end(index_type i) const {
186  return expression().row_end(i);
187  }
188  row_iterator row_begin(index_type i){
189  return expression().row_begin(i);
190  }
191  row_iterator row_end(index_type i){
192  return expression().row_end(i);
193  }
194 
195  const_column_iterator column_begin(index_type j) const {
196  return expression().column_begin(j);
197  }
198  const_column_iterator column_end(index_type j) const {
199  return expression().column_end(j);
200  }
201  column_iterator column_begin(index_type j){
202  return expression().column_begin(j);
203  }
204  column_iterator column_end(index_type j){
205  return expression().column_end(j);
206  }
207 
208  row_iterator set_element(row_iterator pos, index_type index, value_type value) {
209  return expression().set_element(pos,index,value);
210  }
211 
212  row_iterator clear_range(row_iterator start, row_iterator end) {
213  return expression().clear_range(start,end);
214  }
215 
216  row_iterator clear_element(row_iterator elem) {
217  return expression().clear_element(elem);
218  }
219 
220  void clear(){
221  expression().clear();
222  }
223 
224  void reserve(size_type non_zeros) {
225  expression().reserve(non_zeros);
226  }
227 
228  void reserve_row(size_type row, size_type non_zeros) {
229  expression().reserve_row(row,non_zeros);
230  }
231 
232  void reserve_column(size_type column, size_type non_zeros) {
233  expression().reserve_column(column,non_zeros);
234  }
235 
236 
238  expression().swap(m.expression());
239  }
240 
241  friend void swap(matrix_reference& m1, matrix_reference& m2) {
242  m1.swap(m2);
243  }
244 
245  void swap_rows(index_type i, index_type j){
246  expression().swap_rows(i,j);
247  }
248 
249  void swap_columns(index_type i, index_type j){
250  expression().swap_columns(i,j);
251  }
252 private:
253  M* m_expression;
254 };
255 
256 /// \brief Matrix transpose.
257 template<class M>
258 class matrix_transpose: public matrix_expression<matrix_transpose<M> > {
259 public:
260  typedef typename M::size_type size_type;
261  typedef typename M::difference_type difference_type;
262  typedef typename M::value_type value_type;
263  typedef typename M::scalar_type scalar_type;
264  typedef typename M::const_reference const_reference;
265  typedef typename reference<M>::type reference;
266  typedef typename M::const_pointer const_pointer;
267  typedef typename pointer<M>::type pointer;
268 
269  typedef typename M::index_type index_type;
270  typedef typename M::const_index_pointer const_index_pointer;
272 
273  typedef typename closure<M>::type matrix_closure_type;
276  typedef typename M::orientation::transposed_orientation orientation;
277  typedef typename M::storage_category storage_category;
278  typedef elementwise_tag evaluation_category;
279 
280  // Construction and destruction
281  explicit matrix_transpose(matrix_closure_type const& m):
282  m_expression(m) {}
283 
284  //conversion closure->const_closure
285  template<class E>
287  matrix_transpose<E> const& m,
288  typename boost::disable_if<
289  boost::mpl::or_<
290  boost::is_same<matrix_transpose<E>,matrix_transpose>,
291  boost::is_same<matrix_transpose<E>,matrix_closure_type>
292  >
293  >::type* dummy = 0
294  ):m_expression(m.expression()) {}
295 
296  // Expression accessors
297  matrix_closure_type const& expression() const{
298  return m_expression;
299  }
300  matrix_closure_type expression(){
301  return m_expression;
302  }
303 
304  ///\brief Returns the number of rows of the matrix.
305  size_type size1() const {
306  return expression().size2();
307  }
308  ///\brief Returns the number of columns of the matrix.
309  size_type size2() const {
310  return expression().size1();
311  }
312 
313  // ---------
314  // Dense low level interface
315  // ---------
316 
317  ///\brief Returns the stride in memory between two rows.
318  difference_type stride1()const{
319  return expression().stride2();
320  }
321  ///\brief Returns the stride in memory between two columns.
322  difference_type stride2()const{
323  return expression().stride1();
324  }
325 
326  // ---------
327  // Sparse low level interface
328  // ---------
329 
330  /// \brief Number of nonzero elements of the matrix.
331  size_type nnz()const{
332  return expression().nnz();
333  }
334  /// \brief Array of values of the nonzero elements.
335  const_pointer values()const{
336  return expression().values();
337  }
338 
339  /// \brief Array of indices of the nonzero elements.
340  ///
341  /// Note that there is a pair of indices needed:
342  /// When accessing the j-th element in row i you have to write code like this:
343  /// index_type start = outer_indices()[i] //aquire start of the i-th row
344  /// index = inner_indices()[start+j];
345  /// All elements in the row are contained in the range [outer_indices()[i],outer_indices_end()[i])
346  /// there might be gaps between the end of the one line and the start of the next!
347  index_pointer inner_indices()const{
348  return expression().inner_indices();
349  }
350 
351  ///\brief Returns an array containing the start of the rows.
352  ///
353  /// See documentation of inner_indices() for more details.
354  index_pointer outer_indices()const{
355  return expression().outer_indices();
356  }
357 
358  ///\brief Returns an array containing the end of the rows.
359  ///
360  /// See documentation of inner_indices() for more details.
361  index_pointer outer_indices_end()const{
362  return expression().outer_indices_end();
363  }
364 
365  /// \brief Returns the number of nonzero elements in the i-th row/column.
366  size_type inner_nnz(index_type i) const {
367  return expression().inner_nnz(i);
368  }
369 
370  ///\brief Returns the pointer to the beginning of the matrix storage
371  ///
372  /// Grants low-level access to the matrix internals. Element order depends on whether the matrix is row_major or column_major.
373  /// to access element (i,j) use storage()[i*stride1()+j*stride2()].
374  pointer storage()const{
375  return expression().storage();
376  }
377 
378  // ---------
379  // High level interface
380  // ---------
381 
382  // Element access
383  reference operator()(index_type i, index_type j)const{
384  return expression()(j, i);
385  }
386 
387  void set_element(size_type i, size_type j,value_type t){
388  expression().set_element(j,i,t);
389  }
390 
391  typedef typename matrix_closure_type::const_column_iterator const_row_iterator;
392  typedef typename matrix_closure_type::column_iterator row_iterator;
393  typedef typename matrix_closure_type::const_row_iterator const_column_iterator;
394  typedef typename matrix_closure_type::row_iterator column_iterator;
395 
396  //iterators
397  const_row_iterator row_begin(index_type i) const {
398  return expression().column_begin(i);
399  }
400  const_row_iterator row_end(index_type i) const {
401  return expression().column_end(i);
402  }
403  const_column_iterator column_begin(index_type j) const {
404  return expression().row_begin(j);
405  }
406  const_column_iterator column_end(index_type j) const {
407  return expression().row_end(j);
408  }
409 
410  row_iterator row_begin(index_type i) {
411  return expression().column_begin(i);
412  }
413  row_iterator row_end(index_type i) {
414  return expression().column_end(i);
415  }
416  column_iterator column_begin(index_type j) {
417  return expression().row_begin(j);
418  }
419  column_iterator column_end(index_type j) {
420  return expression().row_end(j);
421  }
422 
424 
425  major_iterator set_element(major_iterator pos, index_type index, value_type value){
426  return expression().set_element(pos,index,value);
427  }
428 
429  major_iterator clear_range(major_iterator start, major_iterator end){
430  return expression().clear_range(start,end);
431  }
432 
433  major_iterator clear_element(major_iterator elem){
434  return expression().clear_element(elem);
435  }
436 
437  void clear(){
438  expression().clear();
439  }
440 
441  void reserve(size_type non_zeros) {
442  expression().reserve(non_zeros);
443  }
444 
445  void reserve_row(size_type row, size_type non_zeros) {
446  expression().reserve_row(row,non_zeros);
447  }
448  void reserve_column(size_type column, size_type non_zeros) {
449  expression().reserve_column(column,non_zeros);
450  }
451 
452  // Assignment
453  //we implement it by using the identity A^T op= B <=> A op= B^T where op= is one of =,-=,+=
455  expression() = m.expression();
456  return *this;
457  }
458  template<class E>
461  return *this;
462  }
463 private:
464  matrix_closure_type m_expression;
465 };
466 
467 
468 // (trans m) [i] [j] = m [j] [i]
469 template<class M>
473 }
474 template<class M>
476  return matrix_transpose<M>(m());
477 }
478 
479 template<class M>
481  return trans(static_cast<M&>(m));
482 }
483 
484 template<class M>
485 class matrix_row: public vector_expression<matrix_row<M> > {
486 public:
487  typedef M matrix_type;
488  typedef std::size_t size_type;
489  typedef typename M::difference_type difference_type;
490  typedef typename M::value_type value_type;
491  typedef typename M::scalar_type scalar_type;
492  typedef typename M::const_reference const_reference;
493  typedef typename reference<M>::type reference;
494  typedef typename M::const_pointer const_pointer;
495  typedef typename pointer<M>::type pointer;
496 
497  typedef typename M::index_type index_type;
498  typedef typename M::const_index_pointer const_index_pointer;
500 
501  typedef typename closure<M>::type matrix_closure_type;
502  typedef matrix_row<typename const_expression<M>::type> const_closure_type;
503  typedef matrix_row<M> closure_type;
504  typedef typename M::storage_category storage_category;
505  typedef elementwise_tag evaluation_category;
506 
507  // Construction and destruction
508  matrix_row(matrix_closure_type const& expression, index_type i):m_expression(expression), m_i(i) {
509  SIZE_CHECK (i < expression.size1());
510  }
511 
512  template<class E>
513  matrix_row(matrix_row<E> const& other)
514  :m_expression(other.expression()),m_i(other.index()){}
515 
516  matrix_closure_type const& expression() const {
517  return m_expression;
518  }
519  matrix_closure_type& expression() {
520  return m_expression;
521  }
522 
523  index_type index() const {
524  return m_i;
525  }
526 
527  ///\brief Returns the size of the vector
528  size_type size() const {
529  return expression().size2();
530  }
531 
532  // ---------
533  // Dense low level interface
534  // ---------
535 
536  ///\brief Returns the stride in memory between two elements
537  difference_type stride()const{
538  return expression().stride2();
539  }
540 
541  ///\brief Returns the pointer to the beginning of the vector storage
542  ///
543  /// Grants low-level access to the vector internals.
544  /// to access element i use storage()[i*stride()].
545  pointer storage()const{
546  return expression().storage()+index()*expression().stride1();
547  }
548 
549  // ---------
550  // Sparse low level interface
551  // ---------
552 
553  /// \brief Number of nonzero elements of the vector.
554  size_type nnz()const{
555  return expression().inner_nnz(m_i);
556  }
557  /// \brief Array of values of the nonzero elements.
558  const_pointer values()const{
559  return expression().values()+expression().outer_indices()[m_i];
560  }
561 
562  /// \brief Array of indices of the nonzero elements.
563  index_pointer indices()const{
564  return expression().inner_indices()+expression().outer_indices()[m_i];
565  }
566 
567  // ---------
568  // High level interface
569  // ---------
570 
571  // Element access
572  reference operator()(index_type j) const {
573  return m_expression(m_i, j);
574  }
575  reference operator [](index_type j) const {
576  return (*this)(j);
577  }
578 
579  void set_element(size_type j,value_type t){
580  expression().set_element(m_i,j,t);
581  }
582 
583  // Assignment
584 
585  template<class E>
587  return assign(*this, typename vector_temporary<M>::type(e));
588  }
590  return assign(*this, typename vector_temporary<M>::type(e));
591  }
592 
593  // Iterator types
594  typedef typename M::const_row_iterator const_iterator;
595  typedef typename row_iterator<M>::type iterator;
596 
597  iterator begin() {
598  return expression().row_begin(m_i);
599  }
600  iterator end() {
601  return expression().row_end(m_i);
602  }
603  const_iterator begin()const{
604  return expression().row_begin(m_i);
605  }
606  const_iterator end()const{
607  return expression().row_end(m_i);
608  }
609 
610  iterator set_element(iterator pos, index_type index, value_type value) {
611  return set_element(pos, index, value,
612  typename M::orientation(), typename iterator::iterator_category()
613  );
614  }
615 
616  iterator clear_range(iterator start, iterator end) {
617  return clear_range(start,end,
618  typename M::orientation(), typename iterator::iterator_category()
619  );
620  }
621 
622  iterator clear_element(iterator pos) {
623  return clear_element(pos,
624  typename M::orientation(), typename iterator::iterator_category()
625  );
626  }
627 
628  void clear(){
629  clear_range(begin(),end());
630  }
631 
632  void reserve(size_type non_zeros) {
633  expression().reserve_row(m_i,non_zeros);
634  }
635 
636 private:
637  //we need two implementations of the sparse-interface,
638  //depending on whether M is row or column major.
639 
640  //row major case is trivial
641  template<class Tag>
642  iterator set_element(iterator pos, index_type index, value_type value, row_major, Tag) {
643  return expression().set_element(pos,index,value);
644  }
645  template<class Tag>
646  iterator clear_range(iterator start, iterator end, row_major, Tag) {
647  return expression().clear_range(start,end);
648  }
649  template<class Tag>
650  iterator clear_element(iterator pos, row_major, Tag) {
651  return expression().clear_element(pos);
652  }
653  //dense row major case
654  iterator set_element(iterator pos, index_type index, value_type value,
655  column_major,
656  dense_random_access_iterator_tag
657  ) {
658  RANGE_CHECK(pos.index() == index);
659  *pos = value;
660  return pos;
661  }
662 
663  iterator clear_element(iterator pos,
664  column_major m,
665  dense_random_access_iterator_tag t
666  ) {
667  return set_element(pos,pos.index(),value_type(),m,t);
668  }
669  iterator clear_range(iterator start, iterator end,
670  column_major m,
671  dense_random_access_iterator_tag t
672  ) {
673  for(;start != end; ++start)
674  clear_element(start,m,t);
675  return end;
676  }
677  //todo: sparse column major case.
678 
679  matrix_closure_type m_expression;
680  size_type m_i;
681 };
682 
683 // Projections
684 template<class M>
686  return matrix_row<M> (expression(), i);
687 }
688 template<class M>
690 row(matrix_expression<M> const& expression, typename M::index_type i) {
692 }
693 
694 template<class M>
696  return row(static_cast<M&>(expression), i);
697 }
698 
699 template<class M>
700 class matrix_column: public vector_expression<matrix_column<M> > {
701 public:
702  typedef M matrix_type;
703  typedef std::size_t size_type;
704  typedef typename M::difference_type difference_type;
705  typedef typename M::value_type value_type;
706  typedef typename M::scalar_type scalar_type;
707  typedef typename M::const_reference const_reference;
708  typedef typename reference<M>::type reference;
709  typedef typename M::const_pointer const_pointer;
710  typedef typename pointer<M>::type pointer;
711 
712  typedef typename M::index_type index_type;
713  typedef typename M::const_index_pointer const_index_pointer;
715 
716  typedef typename closure<M>::type matrix_closure_type;
719  typedef typename M::storage_category storage_category;
720  typedef elementwise_tag evaluation_category;
721 
722  // Construction and destruction
723  matrix_column(matrix_closure_type const& expression, index_type j)
724  :m_expression(expression), m_j(j) {
725  SIZE_CHECK (j < expression.size2());
726  }
727 
728  template<class E>
730  :m_expression(other.expression()),m_j(other.index()){}
731 
732  matrix_closure_type const& expression() const {
733  return m_expression;
734  }
735  matrix_closure_type& expression() {
736  return m_expression;
737  }
738 
739  index_type index() const {
740  return m_j;
741  }
742 
743  ///\brief Returns the size of the vector
744  size_type size() const {
745  return expression().size1();
746  }
747 
748  // ---------
749  // Dense low level interface
750  // ---------
751 
752  ///\brief Returns the stride in memory between two elements
753  difference_type stride()const{
754  return expression().stride1();
755  }
756 
757  ///\brief Returns the pointer to the beginning of the vector storage
758  ///
759  /// Grants low-level access to the vector internals.
760  /// to access element i use storage()[i*stride()].
761  pointer storage()const{
762  return expression().storage()+index()*expression().stride2();
763  }
764 
765  // ---------
766  // Sparse low level interface
767  // ---------
768 
769  //~ /// \brief Number of nonzero elements of the vector.
770  //~ size_type nnz()const{
771  //~ return expression().inner_nnz(m_i);
772  //~ }
773  //~ /// \brief Array of values of the nonzero elements.
774  //~ const_pointer values()const{
775  //~ return expression().values()+expression().outer_indices()[m_i];
776  //~ }
777 
778  //~ /// \brief Array of indices of the nonzero elements.
779  //~ index_pointer indices()const{
780  //~ return expression().inner_indices()+expression().outer_indices()[m_i];
781  //~ }
782 
783  // ---------
784  // High level interface
785  // ---------
786 
787  // Element access
788  reference operator()(index_type i) const {
789  return m_expression(i,m_j);
790  }
791  reference operator [](index_type i) const {
792  return (*this)(i);
793  }
794 
795  void set_element(size_type i,value_type t){
796  expression().set_element(i,m_j,t);
797  }
798 
799  // Assignment
800 
801  template<class E>
803  return assign(*this, typename vector_temporary<M>::type(e));
804  }
806  return assign(*this, typename vector_temporary<M>::type(e));
807  }
808 
809  // Iterator types
810  typedef typename M::const_column_iterator const_iterator;
811  typedef typename column_iterator<M>::type iterator;
812 
813  iterator begin() {
814  return expression().column_begin(m_j);
815  }
816  iterator end() {
817  return expression().column_end(m_j);
818  }
819  const_iterator begin()const{
820  return expression().column_begin(m_j);
821  }
822  const_iterator end()const{
823  return expression().column_end(m_j);
824  }
825 
826  iterator set_element(iterator pos, index_type index, value_type value) {
827  return set_element(pos, index, value,
828  typename M::orientation(), typename iterator::iterator_category()
829  );
830  }
831 
832  iterator clear_range(iterator start, iterator end) {
833  return clear_range(start,end,
834  typename M::orientation(), typename iterator::iterator_category()
835  );
836  }
837 
838  iterator clear_element(iterator pos) {
839  return clear_element(pos,
840  typename M::orientation(), typename iterator::iterator_category()
841  );
842  }
843 
844  void clear(){
845  clear_range(begin(),end());
846  }
847 
848  //~ void reserve(size_type non_zeros) {
849  //~ expression().reserve_row(m_i,non_zeros);
850  //~ }
851 
852 private:
853  //we need two implementations of the sparse-interface,
854  //depending on whether M is row or column major.
855 
856  //column major case is trivial
857  template<class Tag>
858  iterator set_element(iterator pos, index_type index, value_type value, column_major, Tag) {
859  return expression().set_element(pos,index,value);
860  }
861  template<class Tag>
862  iterator clear_range(iterator start, iterator end, column_major, Tag) {
863  return expression().clear_range(start,end);
864  }
865  template<class Tag>
866  iterator clear_element(iterator pos, column_major, Tag) {
867  return expression().clear_element(pos);
868  }
869  //dense row major case
870  iterator set_element(iterator pos, index_type index, value_type value,
871  row_major,
872  dense_random_access_iterator_tag
873  ) {
874  RANGE_CHECK(pos.index() == index);
875  *pos = value;
876  return pos;
877  }
878 
879  iterator clear_element(iterator pos,
880  row_major m,
881  dense_random_access_iterator_tag t
882  ) {
883  return set_element(pos,pos.index(),value_type(),m,t);
884  }
885  iterator clear_range(iterator start, iterator end,
886  row_major m,
887  dense_random_access_iterator_tag t
888  ) {
889  for(;start != end; ++start)
890  clear_element(start,m,t);
891  return end;
892  }
893  //todo: sparse column major case.
894 
895  matrix_closure_type m_expression;
896  size_type m_j;
897 };
898 
899 // Projections
900 template<class M>
902  return matrix_column<M> (expression(), j);
903 }
904 template<class M>
906 column(matrix_expression<M> const& expression, typename M::index_type j) {
908 }
909 
910 template<class M>
912  return column(static_cast<M&>(expression), j);
913 }
914 
915 // Matrix based vector range class representing (off-)diagonals of a matrix.
916 template<class M>
917 class matrix_vector_range: public vector_expression<matrix_vector_range<M> > {
918 public:
919  typedef M matrix_type;
920  typedef std::size_t size_type;
921  typedef typename M::difference_type difference_type;
922  typedef typename M::value_type value_type;
923  typedef typename M::scalar_type scalar_type;
924  typedef typename M::const_reference const_reference;
925  typedef typename reference<M>::type reference;
926  typedef typename M::const_pointer const_pointer;
927  typedef typename pointer<M>::type pointer;
928 
929  typedef typename M::index_type index_type;
930  typedef typename M::const_index_pointer const_index_pointer;
932 
933  typedef typename closure<M>::type matrix_closure_type;
936  typedef typename M::storage_category storage_category;
937  typedef elementwise_tag evaluation_category;
938 
939  // Construction and destruction
940  matrix_vector_range(matrix_type& expression, range const&r1, range const&r2):
941  m_expression(expression), m_range1(r1), m_range2(r2) {
942  SIZE_CHECK (m_range1.start() <= expression.size1());
943  SIZE_CHECK (m_range1.start() + m_range1.size () <= expression.size1());
944  SIZE_CHECK (m_range2.start() <= expression.size2());
945  SIZE_CHECK (m_range2.start() + m_range2.size() <= expression.size2());
946  SIZE_CHECK (m_range1.size() == m_range2.size());
947  }
948 
949  template<class E>
951  :m_expression(other.expression()),m_range1(other.range1()),m_range2(other.range2()){}
952 
953  // Accessors
954  size_type start1() const {
955  return m_range1.start();
956  }
957  size_type start2() const {
958  return m_range2.start();
959  }
960 
961  const matrix_closure_type& expression() const {
962  return m_expression;
963  }
964  matrix_closure_type& expression() {
965  return m_expression;
966  }
967 
968  range const& range1() const{
969  return m_range1;
970  }
971  range const& range2() const{
972  return m_range2;
973  }
974 
975  // ---------
976  // Dense low level interface
977  // ---------
978 
979  ///\brief Returns the size of the vector
980  size_type size() const {
981  return m_range1.size();
982  }
983 
984  ///\brief Returns the stride in memory between two elements
985  difference_type stride()const{
986  return expression().stride1()+expression().stride2();
987  }
988 
989  void set_element(size_type i,value_type t){
990  expression().set_element(m_range1(i),m_range2(i),t);
991  }
992 
993 
994  ///\brief Returns the pointer to the beginning of the vector storage
995  ///
996  /// Grants low-level access to the vector internals.
997  /// to access element i use storage()[i*stride()].
998  pointer storage()const{
999  return expression().storage()+start1()*expression().stride1()+start2()*expression().stride1();
1000  }
1001 
1002  // ---------
1003  // High level interface
1004  // ---------
1005 
1006  // Element access
1007  reference operator()(index_type i) const {
1008  return m_expression(m_range1(i), m_range2(i));
1009  }
1010  reference operator [](index_type i) const {
1011  return (*this)(i);
1012  }
1013 
1014  // Assignment
1015 
1016  template<class E>
1018  return assign(*this, typename vector_temporary<M>::type(e));
1019  }
1020 
1021  typedef indexed_iterator<closure_type> iterator;
1022  typedef indexed_iterator<const_closure_type> const_iterator;
1023 
1024  // Element lookup
1025  const_iterator begin()const{
1026  return const_iterator(*this, 0);
1027  }
1028  const_iterator end()const{
1029  return const_iterator(*this, size());
1030  }
1031 
1032  iterator begin() {
1033  return iterator(*this, 0);
1034  }
1035  iterator end() {
1036  return iterator(*this, size());
1037  }
1038 
1039  void reserve(){}
1040  void reserve_row(size_type, size_type) {}
1041  void reserve_column(size_type, size_type ){}
1042 
1043 private:
1044  matrix_closure_type m_expression;
1045  range m_range1;
1046  range m_range2;
1047 };
1048 
1049 ///\brief Returns the diagonal of a constant square matrix as vector.
1050 ///
1051 /// given a matrix
1052 /// A = (1 2 3)
1053 /// (4 5 6)
1054 /// (7 8 9)
1055 ///
1056 /// the diag operation results in
1057 /// diag(A) = (1,5,9)
1058 template<class M>
1060  SIZE_CHECK(mat().size1() == mat().size2());
1061  matrix_vector_range<typename const_expression<M>::type > diagonal(mat(),range(0,mat().size1()),range(0,mat().size1()));
1062  return diagonal;
1063 }
1064 
1065 ///\brief Returns the diagonal of a square matrix as vector.
1066 ///
1067 /// given a matrix
1068 /// A = (1 2 3)
1069 /// (4 5 6)
1070 /// (7 8 9)
1071 ///
1072 /// the diag operation results in
1073 /// diag(A) = (1,5,9)
1074 template<class M>
1076  SIZE_CHECK(mat().size1() == mat().size2());
1077  matrix_vector_range<M> diagonal(mat(),range(0,mat().size1()),range(0,mat().size1()));
1078  return diagonal;
1079 }
1080 
1081 template<class M>
1083  return diag(static_cast<M&>(mat));
1084 }
1085 
1086 // Matrix based range class
1087 template<class M>
1088 class matrix_range:public matrix_expression<matrix_range<M> > {
1089 public:
1090  typedef M matrix_type;
1091  typedef std::size_t size_type;
1092  typedef typename M::difference_type difference_type;
1093  typedef typename M::value_type value_type;
1094  typedef typename M::scalar_type scalar_type;
1095  typedef typename M::const_reference const_reference;
1096  typedef typename reference<M>::type reference;
1097  typedef typename M::const_pointer const_pointer;
1098  typedef typename pointer<M>::type pointer;
1099 
1100  typedef typename M::index_type index_type;
1101  typedef typename M::const_index_pointer const_index_pointer;
1103 
1104  typedef typename closure<M>::type matrix_closure_type;
1105  typedef matrix_range<typename const_expression<M>::type> const_closure_type;
1106  typedef matrix_range<M> closure_type;
1107  typedef typename M::storage_category storage_category;
1108  typedef elementwise_tag evaluation_category;
1109  typedef typename M::orientation orientation;
1110 
1111  // Construction and destruction
1112 
1113  matrix_range(matrix_closure_type expression, range const&r1, range const&r2)
1114  :m_expression(expression), m_range1(r1), m_range2(r2) {
1115  SIZE_CHECK(r1.start() <= expression.size1());
1116  SIZE_CHECK(r1.start() +r1.size() <= expression.size1());
1117  SIZE_CHECK(r2.start() <= expression.size2());
1118  SIZE_CHECK(r2.start() +r2.size() <= expression.size2());
1119  }
1120 
1121  //conversion closure->const_closure
1122  template<class E>
1124  matrix_range<E> const& other,
1125  typename boost::disable_if<
1126  boost::is_same<E,matrix_range>
1127  >::type* dummy = 0
1128  ):m_expression(other.expression())
1129  , m_range1(other.range1())
1130  , m_range2(other.range2()){}
1131 
1132  // Accessors
1133  size_type start1() const {
1134  return m_range1.start();
1135  }
1136  size_type start2() const {
1137  return m_range2.start();
1138  }
1139 
1140  matrix_closure_type expression() const {
1141  return m_expression;
1142  }
1143  matrix_closure_type& expression() {
1144  return m_expression;
1145  }
1146 
1147  range const& range1() const{
1148  return m_range1;
1149  }
1150  range const& range2() const{
1151  return m_range2;
1152  }
1153 
1154  // ---------
1155  // Dense Low level interface
1156  // ---------
1157 
1158  ///\brief Returns the number of rows of the matrix.
1159  size_type size1() const {
1160  return m_range1.size();
1161  }
1162  ///\brief Returns the number of columns of the matrix.
1163  size_type size2() const {
1164  return m_range2.size();
1165  }
1166 
1167  ///\brief Returns the stride in memory between two rows.
1168  difference_type stride1()const{
1169  return expression().stride1();
1170  }
1171  ///\brief Returns the stride in memory between two columns.
1172  difference_type stride2()const{
1173  return expression().stride2();
1174  }
1175 
1176  ///\brief Returns the pointer to the beginning of the matrix storage
1177  ///
1178  /// Grants low-level access to the matrix internals. Element order depends on whether the matrix is row_major or column_major.
1179  /// to access element (i,j) use storage()[i*stride1()+j*stride2()].
1180  pointer storage()const{
1181  return expression().storage()+start1()*stride1()+start2()*stride2();
1182  }
1183 
1184  // ---------
1185  // High level interface
1186  // ---------
1187 
1188 
1189  // Element access
1190  reference operator()(index_type i, index_type j)const{
1191  return m_expression(m_range1(i), m_range2(j));
1192  }
1193 
1194  // Assignment
1195 
1197  return assign(*this, typename matrix_temporary<matrix_range>::type(e));
1198  }
1199  template<class E>
1201  return assign(*this, typename matrix_temporary<E>::type(e));
1202  }
1203 
1204  // Iterator types
1205  typedef subrange_iterator<typename row_iterator<M>::type> row_iterator;
1206  typedef subrange_iterator<typename column_iterator<M>::type> column_iterator;
1207  typedef subrange_iterator<typename M::const_row_iterator> const_row_iterator;
1208  typedef subrange_iterator<typename M::const_column_iterator> const_column_iterator;
1209 
1210  // Element lookup
1211  const_row_iterator row_begin(index_type i) const {
1212  return const_row_iterator(
1213  expression().row_begin(i+start1()),expression().row_end(i+start1()),
1214  start2(),start2()
1215  );
1216  }
1217  row_iterator row_begin(index_type i){
1218  return row_iterator(
1219  expression().row_begin(i+start1()),expression().row_end(i+start1()),
1220  start2(),start2()
1221  );
1222  }
1223  const_row_iterator row_end(index_type i) const {
1224  return const_row_iterator(
1225  expression().row_begin(i+start1()),expression().row_end(i+start1()),
1226  start2()+size2(),start2()
1227  );
1228  }
1229  row_iterator row_end(index_type i){
1230  return row_iterator(
1231  expression().row_begin(i+start1()),expression().row_end(i+start1()),
1232  start2()+size2(),start2()
1233  );
1234  }
1235  const_column_iterator column_begin(index_type j) const {
1236  return const_column_iterator(
1237  expression().column_begin(j+start2()),expression().column_end(j+start2()),
1238  start1(),start1()
1239  );
1240  }
1241  column_iterator column_begin(index_type j) {
1242  return column_iterator(
1243  expression().column_begin(j+start2()),expression().column_end(j+start2()),
1244  start1(),start1()
1245  );
1246  }
1247  const_column_iterator column_end(index_type j) const {
1248  return const_column_iterator(
1249  expression().column_begin(j+start2()),expression().column_end(j+start2()),
1250  start1()+size1(),start1()
1251  );
1252  }
1253  column_iterator column_end(index_type j) {
1254  return column_iterator(
1255  expression().column_begin(j+start2()),expression().column_end(j+start2()),
1256  start1()+size1(),start1()
1257  );
1258  }
1260 
1261  major_iterator set_element(major_iterator pos, index_type index, value_type value) {
1262  return expression().set_element(pos.inner(),index+orientation::index_m(start1(),start2()),value);
1263  }
1264 
1265  major_iterator clear_element(major_iterator elem) {
1266  return major_iterator(expression().clear_element(elem.inner()),orientation::size_m(start1(),start2()));
1267  }
1268 
1269  major_iterator clear_range(major_iterator start, major_iterator end) {
1270  return major_iterator(expression().clear_range(start.inner(),end.inner()),orientation::size_m(start1(),start2()));
1271  }
1272 
1273  void clear(){
1274  for(index_type i = 0; i != orientation::index_M(size1(),size2()); ++i)
1275  clear_range(major_begin(*this,i),major_end(*this,i));
1276  }
1277 
1278  void reserve(size_type){}
1279  void reserve_row(size_type, size_type) {}
1280  void reserve_column(size_type, size_type ){}
1281 private:
1282  matrix_closure_type m_expression;
1283  range m_range1;
1284  range m_range2;
1285 };
1286 
1287 // Simple Projections
1288 template<class M>
1291  std::size_t start1, std::size_t stop1,
1292  std::size_t start2, std::size_t stop2
1293 ) {
1294  RANGE_CHECK(start1 <= stop1);
1295  RANGE_CHECK(start2 <= stop2);
1296  SIZE_CHECK(stop1 <= expression().size1());
1297  SIZE_CHECK(stop2 <= expression().size2());
1298  return matrix_range<M> (expression(), range(start1, stop1), range(start2, stop2));
1299 }
1300 template<class M>
1303  std::size_t start1, std::size_t stop1,
1304  std::size_t start2, std::size_t stop2
1305 ) {
1306  RANGE_CHECK(start1 <= stop1);
1307  RANGE_CHECK(start2 <= stop2);
1308  SIZE_CHECK(stop1 <= expression().size1());
1309  SIZE_CHECK(stop2 <= expression().size2());
1310  return matrix_range<typename const_expression<M>::type> (expression(), range(start1, stop1), range(start2, stop2));
1311 }
1312 
1313 template<class M>
1316  std::size_t start1, std::size_t stop1,
1317  std::size_t start2, std::size_t stop2
1318 ) {
1319  return subrange(static_cast<M&>(expression),start1,stop1,start2,stop2);
1320 }
1321 
1322 template<class M>
1325  std::size_t start, std::size_t stop
1326 ) {
1327  RANGE_CHECK(start <= stop);
1328  SIZE_CHECK(stop <= expression().size1());
1329  return subrange(expression, start, stop, 0,expression().size2());
1330 }
1331 
1332 template<class M>
1335  std::size_t start, std::size_t stop
1336 ) {
1337  RANGE_CHECK(start <= stop);
1338  SIZE_CHECK(stop <= expression().size1());
1339  return subrange(expression, start, stop, 0,expression().size2());
1340 }
1341 
1342 template<class M>
1345  std::size_t start, std::size_t stop
1346 ) {
1347  return rows(static_cast<M&>(expression),start,stop);
1348 }
1349 
1350 template<class M>
1353  typename M::index_type start, typename M::index_type stop
1354 ) {
1355  RANGE_CHECK(start <= stop);
1356  SIZE_CHECK(stop <= expression().size2());
1357  return subrange(expression, 0,expression().size1(), start, stop);
1358 }
1359 
1360 template<class M>
1363  typename M::index_type start, typename M::index_type stop
1364 ) {
1365  RANGE_CHECK(start <= stop);
1366  SIZE_CHECK(stop <= expression().size2());
1367  return subrange(expression, 0,expression().size1(), start, stop);
1368 }
1369 
1370 template<class M>
1373  std::size_t start, std::size_t stop
1374 ) {
1375  return columns(static_cast<M&>(expression),start,stop);
1376 }
1377 
1378 template<class T,class Orientation=row_major>
1379 class dense_matrix_adaptor: public matrix_expression<dense_matrix_adaptor<T,Orientation> > {
1381 public:
1382 
1383  //std::container types
1384  typedef typename Orientation::orientation orientation;
1385  typedef std::size_t size_type;
1386  typedef std::ptrdiff_t difference_type;
1387  typedef typename boost::remove_const<T>::type value_type;
1388  typedef value_type scalar_type;
1389  typedef value_type const& const_reference;
1390  typedef T& reference;
1391  typedef T* pointer;
1392  typedef value_type const* const_pointer;
1393 
1394  typedef std::size_t index_type;
1395  typedef index_type const* const_index_pointer;
1396  typedef index_type* index_pointer;
1397 
1398  //ublas types
1401  typedef dense_tag storage_category;
1402  typedef elementwise_tag evaluation_category;
1403 
1404 
1405  // Construction and destruction
1406 
1408  : m_values(expression.storage())
1409  , m_size1(expression.size1())
1410  , m_size2(expression.size2())
1411  , m_stride1(expression.stride1())
1412  , m_stride2(expression.stride2())
1413  {}
1414 
1415  /// \brief Constructor of a vector proxy from a Dense MatrixExpression
1416  ///
1417  /// Be aware that the expression must live longer than the proxy!
1418  /// \param expression Expression from which to construct the Proxy
1419  template<class E>
1421  : m_values(expression().storage())
1422  , m_size1(expression().size1())
1423  , m_size2(expression().size2())
1424  , m_stride1(expression().stride1())
1425  , m_stride2(expression().stride2())
1426  {
1427  BOOST_STATIC_ASSERT((
1428  boost::is_same<typename E::orientation,orientation>::value
1429  ));
1430  }
1431 
1432  /// \brief Constructor of a vector proxy from a Dense MatrixExpression
1433  ///
1434  /// Be aware that the expression must live longer than the proxy!
1435  /// \param expression Expression from which to construct the Proxy
1436  template<class E>
1438  : m_values(expression().storage())
1439  , m_size1(expression().size1())
1440  , m_size2(expression().size2())
1441  , m_stride1(expression().stride1())
1442  , m_stride2(expression().stride2())
1443  {
1444  BOOST_STATIC_ASSERT(
1445  (boost::is_same<typename E::orientation,orientation>::value)
1446  );
1447  }
1448 
1449  /// \brief Constructor of a vector proxy from a block of memory
1450  /// \param values the block of memory used
1451  /// \param size1 size in 1st direction
1452  /// \param size2 size in 2nd direction
1453  /// \param stride1 distance in 1st direction between elements of the self_type in memory
1454  /// \param stride2 distance in 2nd direction between elements of the self_type in memory
1456  pointer values,
1457  size_type size1, size_type size2,
1458  difference_type stride1 = 0,difference_type stride2 = 0
1459  )
1460  : m_values(values)
1461  , m_size1(size1)
1462  , m_size2(size2)
1463  , m_stride1(stride1)
1464  , m_stride2(stride2)
1465  {
1466  if(!m_stride1)
1467  m_stride1= Orientation::stride1(m_size1,m_size2);
1468  if(!m_stride2)
1469  m_stride2= Orientation::stride2(m_size1,m_size2);
1470  }
1471 
1472  // ---------
1473  // Dense low level interface
1474  // ---------
1475 
1476  /// \brief Return the number of rows of the matrix
1477  size_type size1() const {
1478  return m_size1;
1479  }
1480  /// \brief Return the number of columns of the matrix
1481  size_type size2() const {
1482  return m_size2;
1483  }
1484 
1485  difference_type stride1()const{
1486  return m_stride1;
1487  }
1488  difference_type stride2()const{
1489  return m_stride2;
1490  }
1491 
1492  pointer storage()const{
1493  return m_values;
1494  }
1495 
1496  // ---------
1497  // High level interface
1498  // ---------
1499 
1500  // -------
1501  // ASSIGNING
1502  // -------
1503 
1504  self_type& operator = (self_type const& e) {
1505  SIZE_CHECK(size1() == e().size1());
1506  SIZE_CHECK(size2() == e().size2());
1507  return assign(*this, typename matrix_temporary<self_type>::type(e));
1508  }
1509  template<class E>
1510  self_type& operator = (matrix_expression<E> const& e) {
1511  SIZE_CHECK(size1() == e().size1());
1512  SIZE_CHECK(size2() == e().size2());
1513  return assign(*this, typename matrix_temporary<self_type>::type(e));
1514  }
1515 
1516  // --------------
1517  // Element access
1518  // --------------
1519 
1520  const_reference operator () (index_type i, index_type j) const {
1521  return m_values[i*m_stride1+j*m_stride2];
1522  }
1523  reference operator () (index_type i, index_type j) {
1524  return m_values[i*m_stride1+j*m_stride2];
1525  }
1526  void set_element(size_type i, size_type j,value_type t){
1527  m_values[i*m_stride1+j*m_stride2] = t;
1528  }
1529 
1530 
1531  // --------------
1532  // ITERATORS
1533  // --------------
1534 
1535 
1536  typedef dense_storage_iterator<T> row_iterator;
1537  typedef dense_storage_iterator<T> column_iterator;
1538  typedef dense_storage_iterator<value_type const> const_row_iterator;
1539  typedef dense_storage_iterator<value_type const> const_column_iterator;
1540 
1541  const_row_iterator row_begin(index_type i) const {
1542  return const_row_iterator(m_values+i*stride1(),0,stride2());
1543  }
1544  const_row_iterator row_end(index_type i) const {
1545  return const_row_iterator(m_values+i*stride1()+size2()*stride2(),size2(),stride2());
1546  }
1547  row_iterator row_begin(index_type i){
1548  return row_iterator(m_values+i*stride1(),0,stride2());
1549  }
1550  row_iterator row_end(index_type i){
1551  return row_iterator(m_values+i*stride1()+size2()*stride2(),size2(),stride2());
1552  }
1553 
1554  const_column_iterator column_begin(index_type j) const {
1555  return const_column_iterator(m_values+j*stride2(),0,stride1());
1556  }
1557  const_column_iterator column_end(index_type j) const {
1558  return const_column_iterator(m_values+j*stride2()+size1()*stride1(),size1(),stride1());
1559  }
1560  column_iterator column_begin(index_type j){
1561  return column_iterator(m_values+j*stride2(),0,stride1());
1562  }
1563  column_iterator column_end(index_type j){
1564  return column_iterator(m_values+j*stride2()+size1()*stride1(),size1(),stride1());
1565  }
1566 
1568 
1569  major_iterator set_element(major_iterator pos, index_type index, value_type value) {
1570  RANGE_CHECK(pos.index() == index);
1571  *pos=value;
1572  return pos;
1573  }
1574 
1575  major_iterator clear_element(major_iterator elem) {
1576  *elem = value_type();
1577  return elem+1;
1578  }
1579 
1580  major_iterator clear_range(major_iterator start, major_iterator end) {
1581  std::fill(start,end,value_type());
1582  return end;
1583  }
1584 
1585 
1586  void clear(){
1587  for(index_type i = 0; i != size1(); ++i){
1588  for(index_type j = 0; j != size2(); ++j){
1589  (*this)(i,j) = value_type();
1590  }
1591  }
1592  }
1593 private:
1594  pointer m_values;
1595  size_type m_size1;
1596  size_type m_size2;
1597  difference_type m_stride1;
1598  difference_type m_stride2;
1599 };
1600 
1601 /// \brief Converts a chunk of memory into a matrix of given size.
1602 template <class T>
1604  return dense_matrix_adaptor<T>(data,size1, size2);
1605 }
1606 
1607 /// \brief Converts a 2D C-style array into a matrix of given size.
1608 template <class T, std::size_t M, std::size_t N>
1610  return dense_matrix_adaptor<T>(&(array[0][0]),M,N);
1611 }
1612 
1613 /// \brief Converts a dense vector to a matrix of a given size
1614 template <class V>
1615 typename boost::enable_if<
1616  boost::is_same<typename V::storage_category,dense_tag>,
1618  typename boost::remove_reference<typename V::reference>::type
1619  > >
1620 >::type
1623  std::size_t size1, std::size_t size2
1624 ){
1625  typedef typename boost::remove_reference<typename V::reference>::type ElementType;
1626  return dense_matrix_adaptor<ElementType>(v().storage(), size1, size2);
1627 }
1628 
1629 /// \brief Converts a dense vector to a matrix of a given size
1630 template <class V>
1631 typename boost::enable_if<
1632  boost::is_same<typename V::storage_category,dense_tag>,
1634 >::type
1636  vector_expression<V> const& v,
1637  std::size_t size1, std::size_t size2
1638 ){
1640 }
1641 
1642 template <class E>
1643 typename boost::enable_if<
1644  boost::is_same<typename E::storage_category,dense_tag>,
1646  typename boost::remove_reference<typename E::reference>::type
1647  > >
1648 >::type
1651  std::size_t size1, std::size_t size2
1652 ){
1653  return to_matrix(static_cast<E&>(v),size1,size2);
1654 }
1655 
1656 }}
1657 #endif