matrix_expression.hpp
Go to the documentation of this file.
1 /*!
2  * \brief Expression templates for expressions involving matrices
3  *
4  * \author O. Krause
5  * \date 2013
6  *
7  *
8  * \par Copyright 1995-2015 Shark Development Team
9  *
10  * <BR><HR>
11  * This file is part of Shark.
12  * <http://image.diku.dk/shark/>
13  *
14  * Shark is free software: you can redistribute it and/or modify
15  * it under the terms of the GNU Lesser General Public License as published
16  * by the Free Software Foundation, either version 3 of the License, or
17  * (at your option) any later version.
18  *
19  * Shark is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU Lesser General Public License for more details.
23  *
24  * You should have received a copy of the GNU Lesser General Public License
25  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
26  *
27  */
28 #ifndef SHARK_LINALG_BLAS_MATRIX_EXPRESSION_HPP
29 #define SHARK_LINALG_BLAS_MATRIX_EXPRESSION_HPP
30 
31 #include <boost/type_traits/is_convertible.hpp>
32 #include "matrix_proxy.hpp"
33 #include "operation.hpp"
34 
35 namespace shark {
36 namespace blas {
37 
38 template<class E1, class E2>
39 class outer_product:public matrix_expression<outer_product<E1, E2> > {
40  typedef scalar_multiply1<typename E1::value_type, typename E2::value_type> functor_type;
41 public:
42  typedef typename E1::const_closure_type lhs_closure_type;
43  typedef typename E2::const_closure_type rhs_closure_type;
44 
45 
46  typedef std::size_t size_type;
47  typedef std::ptrdiff_t difference_type;
48  typedef typename functor_type::result_type value_type;
49  typedef value_type scalar_type;
50  typedef value_type const_reference;
51  typedef const_reference reference;
52  typedef value_type const* const_pointer;
53  typedef const_pointer pointer;
54 
55  typedef typename E1::index_type index_type;
56  typedef typename E1::const_index_pointer const_index_pointer;
58 
61  typedef unknown_orientation orientation;
62  typedef unknown_storage_tag storage_category;
63  typedef typename evaluation_restrict_traits<E1,E2>::type evaluation_category;
64 
65  // Construction and destruction
66 
67  outer_product(lhs_closure_type const& e1, rhs_closure_type const& e2)
68  :m_lhs(e1), m_rhs(e2){}
69 
70  // Accessors
71  size_type size1() const {
72  return m_lhs.size();
73  }
74  size_type size2() const {
75  return m_rhs.size();
76  }
77 
78  // Expression accessors
79  const lhs_closure_type &lhs() const {
80  return m_lhs;
81  }
82  const rhs_closure_type &rhs() const {
83  return m_rhs;
84  }
85  // Element access
86  const_reference operator()(index_type i, index_type j) const {
87  return m_lhs(i) * m_rhs(j);
88  }
89 
90  typedef transform_iterator<typename E2::const_iterator,functor_type> const_row_iterator;
91  typedef transform_iterator<typename E1::const_iterator,functor_type> const_column_iterator;
92  typedef const_row_iterator row_iterator;
93  typedef const_column_iterator column_iterator;
94 
95  const_row_iterator row_begin(index_type i) const {
96  return const_row_iterator(m_rhs.begin(),
97  functor_type(m_lhs(i))
98  );
99  }
100  const_row_iterator row_end(index_type i) const {
101  return const_row_iterator(m_rhs.end(),
102  functor_type(m_lhs(i))
103  );
104  }
105 
106  const_column_iterator column_begin(index_type i) const {
107  return const_column_iterator(m_lhs.begin(),
108  functor_type(m_rhs(i))
109  );
110  }
111  const_column_iterator column_end(index_type i) const {
112  return const_column_iterator(m_lhs.end(),
113  functor_type(m_rhs(i))
114  );
115  }
116 private:
117  lhs_closure_type m_lhs;
118  rhs_closure_type m_rhs;
119 };
120 
121 // (outer_prod (v1, v2)) [i] [j] = v1 [i] * v2 [j]
122 template<class E1, class E2>
125  vector_expression<E1> const& e1,
126  vector_expression<E2> const& e2
127 ) {
128  return outer_product<E1, E2>(e1(), e2());
129 }
130 
131 template<class V>
132 class vector_repeater:public blas::matrix_expression<vector_repeater<V> > {
133 private:
134  typedef V expression_type;
136  typedef typename V::const_iterator const_subiterator_type;
137 public:
138  typedef typename V::const_closure_type expression_closure_type;
139  typedef typename V::size_type size_type;
140  typedef typename V::difference_type difference_type;
141  typedef typename V::value_type value_type;
142  typedef value_type scalar_type;
143  typedef value_type const_reference;
144  typedef const_reference reference;
145  typedef value_type const* pointer;
146  typedef value_type const* const_pointer;
147 
148  typedef typename V::index_type index_type;
149  typedef typename V::const_index_pointer const_index_pointer;
151 
152  typedef self_type const const_closure_type;
153  typedef const_closure_type closure_type;
154  typedef blas::row_major orientation;
155  typedef blas::unknown_storage_tag storage_category;
156  typedef typename V::evaluation_category evaluation_category;
157 
158  // Construction and destruction
159  explicit vector_repeater (expression_type const& e, std::size_t rows):
160  m_vector(e), m_rows(rows) {}
161 
162  // Accessors
163  size_type size1() const {
164  return m_rows;
165  }
166  size_type size2() const {
167  return m_vector.size();
168  }
169 
170  // Expression accessors
171  const expression_closure_type &expression () const {
172  return m_vector;
173  }
174 
175  // Element access
176  const_reference operator() (index_type i, index_type j) const {
177  return m_vector(j);
178  }
179 
180  // Iterator types
181  typedef typename V::const_iterator const_row_iterator;
182  typedef const_row_iterator row_iterator;
183  typedef blas::constant_iterator<value_type> const_column_iterator;
184  typedef const_column_iterator column_iterator;
185 
186  // Element lookup
187 
188  const_row_iterator row_begin(std::size_t i) const {
189  RANGE_CHECK( i < size1());
190  return m_vector.begin();
191  }
192  const_row_iterator row_end(std::size_t i) const {
193  RANGE_CHECK( i < size1());
194  return m_vector.end();
195  }
196 
197  const_column_iterator column_begin(std::size_t j) const {
198  RANGE_CHECK( j < size2());
199  return const_column_iterator(0,m_vector(j));
200  }
201  const_column_iterator column_end(std::size_t j) const {
202  RANGE_CHECK( j < size2());
203  return const_column_iterator(size1(),m_vector(j));
204  }
205 private:
206  expression_closure_type m_vector;
207  std::size_t m_rows;
208 };
209 
210 ///\brief Creates a matrix from a vector by repeating the vector in every row of the matrix
211 ///
212 ///example: vector = (1,2,3)
213 ///repeat(vector,3) results in
214 ///(1,2,3)
215 ///(1,2,3)
216 ///(1,2,3)
217 ///@param vector the vector which is to be repeated as the rows of the resulting matrix
218 ///@param rows the number of rows of the matrix
219 template<class Vector>
221  return vector_repeater<Vector>(vector(),rows);
222 }
223 
224 /** \brief A matrix with all values of type \c T equal to the same value
225  *
226  * \tparam T the type of object stored in the matrix (like double, float, complex, etc...)
227  */
228 template<class T>
230  public matrix_container<scalar_matrix<T> > {
231 
232  typedef scalar_matrix<T> self_type;
233 public:
234  typedef std::size_t size_type;
235  typedef std::ptrdiff_t difference_type;
236  typedef T value_type;
237  typedef const T& const_reference;
238  typedef const_reference reference;
239  typedef value_type const* const_pointer;
240  typedef value_type scalar_type;
241  typedef const_pointer pointer;
242 
243  typedef std::size_t index_type;
244  typedef index_type const* const_index_pointer;
245  typedef index_type index_pointer;
246 
247  typedef self_type const_closure_type;
248  typedef self_type closure_type;
249  typedef dense_tag storage_category;
250  typedef unknown_orientation orientation;
251  typedef elementwise_tag evaluation_category;
252 
253  // Construction and destruction
255  m_size1(0), m_size2(0), m_value() {}
256  scalar_matrix(size_type size1, size_type size2, const value_type& value = value_type(1)):
257  m_size1(size1), m_size2(size2), m_value(value) {}
259  m_size1(m.m_size1), m_size2(m.m_size2), m_value(m.m_value) {}
260 
261  // Accessors
262  size_type size1() const {
263  return m_size1;
264  }
265  size_type size2() const {
266  return m_size2;
267  }
268 
269  // Element access
270  const_reference operator()(size_type /*i*/, size_type /*j*/) const {
271  return m_value;
272  }
273 
274  //Iterators
275  typedef constant_iterator<value_type> const_row_iterator;
276  typedef constant_iterator<value_type> const_column_iterator;
277  typedef const_row_iterator row_iterator;
278  typedef const_column_iterator column_iterator;
279 
280  const_row_iterator row_begin(std::size_t i) const {
281  return const_row_iterator(0, m_value);
282  }
283  const_row_iterator row_end(std::size_t i) const {
284  return const_row_iterator(size2(), m_value);
285  }
286 
287  const_row_iterator column_begin(std::size_t j) const {
288  return const_row_iterator(0, m_value);
289  }
290  const_row_iterator column_end(std::size_t j) const {
291  return const_row_iterator(size1(), m_value);
292  }
293 private:
294  size_type m_size1;
295  size_type m_size2;
296  value_type m_value;
297 };
298 
299 ///brief repeats a single element to form a matrix of size rows x columns
300 ///
301 ///@param scalar the value which is repeated
302 ///@param rows the number of rows of the resulting vector
303 ///@param columns the number of columns of the resulting vector
304 template<class T>
305 typename boost::enable_if<boost::is_arithmetic<T>, scalar_matrix<T> >::type
306 repeat(T scalar, std::size_t rows, std::size_t columns){
307  return scalar_matrix<T>(rows, columns, scalar);
308 }
309 
310 template<class E>
311 class matrix_scalar_multiply:public blas::matrix_expression<matrix_scalar_multiply<E> > {
312 private:
313  typedef typename E::const_row_iterator const_subrow_iterator_type;
314  typedef typename E::const_column_iterator const_subcolumn_iterator_type;
315  typedef scalar_multiply1<typename E::value_type, typename E::scalar_type> functor_type;
316 public:
317  typedef typename E::const_closure_type expression_closure_type;
318 
319  typedef typename functor_type::result_type value_type;
320  typedef typename E::scalar_type scalar_type;
321  typedef value_type const_reference;
322  typedef const_reference reference;
323  typedef value_type const *const_pointer;
324  typedef value_type *pointer;
325  typedef typename E::size_type size_type;
326  typedef typename E::difference_type difference_type;
327 
328  typedef typename E::index_type index_type;
329  typedef typename E::const_index_pointer const_index_pointer;
331 
332  typedef matrix_scalar_multiply const_closure_type;
333  typedef const_closure_type closure_type;
334  typedef typename E::orientation orientation;
335  typedef blas::unknown_storage_tag storage_category;
336  typedef typename E::evaluation_category evaluation_category;
337 
338  // Construction and destruction
339  matrix_scalar_multiply(blas::matrix_expression<E> const &e, scalar_type scalar):
340  m_expression(e()), m_scalar(scalar){}
341 
342  // Accessors
343  size_type size1() const {
344  return m_expression.size1();
345  }
346  size_type size2() const {
347  return m_expression.size2();
348  }
349 
350  // Element access
351  const_reference operator()(index_type i, index_type j) const {
352  return m_scalar * m_expression(i, j);
353  }
354 
355  //computation kernels
356  template<class MatX>
357  void assign_to(matrix_expression<MatX>& X, scalar_type alpha = scalar_type(1) )const{
358  m_expression.assign_to(X,alpha*m_scalar);
359  }
360  template<class MatX>
361  void plus_assign_to(matrix_expression<MatX>& X, scalar_type alpha = scalar_type(1) )const{
362  m_expression.plus_assign_to(X,alpha*m_scalar);
363  }
364 
365  template<class MatX>
366  void minus_assign_to(matrix_expression<MatX>& X, scalar_type alpha = scalar_type(1) )const{
367  m_expression.minus_assign_to(X,alpha*m_scalar);
368  }
369 
370  // Iterator types
371  typedef transform_iterator<typename E::const_row_iterator, functor_type> const_row_iterator;
372  typedef transform_iterator<typename E::const_column_iterator, functor_type> const_column_iterator;
373  typedef const_row_iterator row_iterator;
374  typedef const_column_iterator column_iterator;
375 
376  const_row_iterator row_begin(index_type i) const {
377  return const_row_iterator(m_expression.row_begin(i),functor_type(m_scalar));
378  }
379  const_row_iterator row_end(index_type i) const {
380  return const_row_iterator(m_expression.row_end(i),functor_type(m_scalar));
381  }
382 
383  const_column_iterator column_begin(index_type i) const {
384  return const_row_iterator(m_expression.column_begin(i),functor_type(m_scalar));
385  }
386  const_column_iterator column_end(index_type i) const {
387  return const_row_iterator(m_expression.column_end(i),functor_type(m_scalar));
388  }
389 
390 private:
391  expression_closure_type m_expression;
392  scalar_type m_scalar;
393 };
394 
395 
396 ///\brief class which allows for matrix transformations
397 ///
398 ///transforms a matrix expression e of type E using a Functiof f of type F as an elementwise transformation f(e(i,j))
399 ///This transformation needs to leave f constant, meaning that applying f(x), f(y), f(z) yields the same results independent of the
400 ///order of application. Also F must provide a type F::result_type indicating the result type of the functor.
401 ///F must further provide a boolean flag F::zero_identity which indicates that f(0) = 0. This is needed for correct usage with sparse
402 ///arguments - if f(0) != 0 this expression will be dense!
403 ///todo: desification is not implemented
404 template<class E, class F>
405 class matrix_unary:public blas::matrix_expression<matrix_unary<E, F> > {
406 private:
407  typedef matrix_unary<E, F> self_type;
408  typedef E expression_type;
409  typedef typename expression_type::const_row_iterator const_subrow_iterator_type;
410  typedef typename expression_type::const_column_iterator const_subcolumn_iterator_type;
411 
412 public:
413  typedef typename expression_type::const_closure_type expression_closure_type;
414 
415  typedef F functor_type;
416  typedef typename functor_type::result_type value_type;
417  typedef value_type scalar_type;
418  typedef value_type const_reference;
419  typedef const_reference reference;
420  typedef value_type const *const_pointer;
421  typedef value_type *pointer;
422  typedef typename expression_type::size_type size_type;
423  typedef typename expression_type::difference_type difference_type;
424 
425  typedef typename E::index_type index_type;
426  typedef typename E::const_index_pointer const_index_pointer;
428 
429  typedef self_type const const_closure_type;
430  typedef const_closure_type closure_type;
431  typedef typename E::orientation orientation;
432  typedef blas::unknown_storage_tag storage_category;
433  typedef typename E::evaluation_category evaluation_category;
434 
435  // Construction and destruction
436  matrix_unary(blas::matrix_expression<E> const &e, F const &functor):
437  m_expression(e()), m_functor(functor) {}
438 
439  // Accessors
440  size_type size1() const {
441  return m_expression.size1();
442  }
443  size_type size2() const {
444  return m_expression.size2();
445  }
446 
447  // Element access
448  const_reference operator()(index_type i, index_type j) const {
449  return m_functor(m_expression(i, j));
450  }
451 
452  // Iterator types
453  typedef transform_iterator<typename E::const_row_iterator, F> const_row_iterator;
454  typedef transform_iterator<typename E::const_column_iterator, F> const_column_iterator;
455  typedef const_row_iterator row_iterator;
456  typedef const_column_iterator column_iterator;
457 
458  const_row_iterator row_begin(index_type i) const {
459  return const_row_iterator(m_expression.row_begin(i),m_functor);
460  }
461  const_row_iterator row_end(index_type i) const {
462  return const_row_iterator(m_expression.row_end(i),m_functor);
463  }
464 
465  const_column_iterator column_begin(index_type i) const {
466  return const_row_iterator(m_expression.column_begin(i),m_functor);
467  }
468  const_column_iterator column_end(index_type i) const {
469  return const_row_iterator(m_expression.column_end(i),m_functor);
470  }
471 
472 private:
473  expression_closure_type m_expression;
474  functor_type m_functor;
475 };
476 
477 template<class E, class T>
478 typename boost::enable_if<
479  boost::is_convertible<T, typename E::scalar_type >,
481 >::type
482 operator* (matrix_expression<E> const& e, T scalar){
483  return matrix_scalar_multiply<E>(e(), typename E::scalar_type(scalar));
484 }
485 
486 template<class T, class E>
487 typename boost::enable_if<
488  boost::is_convertible<T, typename E::scalar_type >,
490 >::type
491 operator* (T scalar, matrix_expression<E> const& e){
492  return matrix_scalar_multiply<E>(e(), typename E::scalar_type(scalar));
493 }
494 
495 template<class E>
497  return matrix_scalar_multiply<E>(e(), typename E::scalar_type(-1));
498 }
499 
500 #define SHARK_UNARY_MATRIX_TRANSFORMATION(name, F)\
501 template<class E>\
502 matrix_unary<E,F<typename E::value_type> >\
503 name(matrix_expression<E> const& e){\
504  typedef F<typename E::value_type> functor_type;\
505  return matrix_unary<E, functor_type>(e, functor_type());\
506 }
507 SHARK_UNARY_MATRIX_TRANSFORMATION(conj, scalar_conj)
508 SHARK_UNARY_MATRIX_TRANSFORMATION(real, scalar_real)
509 SHARK_UNARY_MATRIX_TRANSFORMATION(imag, scalar_imag)
510 SHARK_UNARY_MATRIX_TRANSFORMATION(abs, scalar_abs)
511 SHARK_UNARY_MATRIX_TRANSFORMATION(log, scalar_log)
512 SHARK_UNARY_MATRIX_TRANSFORMATION(exp, scalar_exp)
513 SHARK_UNARY_MATRIX_TRANSFORMATION(sin, scalar_sin)
514 SHARK_UNARY_MATRIX_TRANSFORMATION(cos, scalar_cos)
515 SHARK_UNARY_MATRIX_TRANSFORMATION(tanh,scalar_tanh)
516 SHARK_UNARY_MATRIX_TRANSFORMATION(atanh,scalar_atanh)
518 SHARK_UNARY_MATRIX_TRANSFORMATION(abs_sqr, scalar_abs_sqr)
519 SHARK_UNARY_MATRIX_TRANSFORMATION(sqrt, scalar_sqrt)
522 SHARK_UNARY_MATRIX_TRANSFORMATION(elem_inv, scalar_inverse)
523 #undef SHARK_UNARY_MATRIX_TRANSFORMATION
524 
525 #define SHARK_MATRIX_SCALAR_TRANSFORMATION(name, F)\
526 template<class E, class T> \
527 typename boost::enable_if< \
528  boost::is_convertible<T, typename E::value_type >,\
529  matrix_unary<E,F<typename E::value_type,T> > \
530 >::type \
531 name (matrix_expression<E> const& e, T scalar){ \
532  typedef F<typename E::value_type, T> functor_type; \
533  return matrix_unary<E, functor_type>(e, functor_type(scalar)); \
534 }
535 SHARK_MATRIX_SCALAR_TRANSFORMATION(operator/, scalar_divide)
536 SHARK_MATRIX_SCALAR_TRANSFORMATION(operator<, scalar_less_than)
537 SHARK_MATRIX_SCALAR_TRANSFORMATION(operator<=, scalar_less_equal_than)
538 SHARK_MATRIX_SCALAR_TRANSFORMATION(operator>, scalar_bigger_than)
539 SHARK_MATRIX_SCALAR_TRANSFORMATION(operator>=, scalar_bigger_equal_than)
540 SHARK_MATRIX_SCALAR_TRANSFORMATION(operator==, scalar_equal)
541 SHARK_MATRIX_SCALAR_TRANSFORMATION(operator!=, scalar_not_equal)
544 SHARK_MATRIX_SCALAR_TRANSFORMATION(pow, scalar_pow)
545 #undef SHARK_MATRIX_SCALAR_TRANSFORMATION
546 
547 // operations of the form op(t,v)[i,j] = op(t,v[i,j])
548 #define SHARK_MATRIX_SCALAR_TRANSFORMATION_2(name, F)\
549 template<class T, class E> \
550 typename boost::enable_if< \
551  boost::is_convertible<T, typename E::value_type >,\
552  matrix_unary<E,F<typename E::value_type,T> > \
553 >::type \
554 name (T scalar, matrix_expression<E> const& e){ \
555  typedef F<typename E::value_type, T> functor_type; \
556  return matrix_unary<E, functor_type>(e, functor_type(scalar)); \
557 }
560 #undef SHARK_MATRIX_SCALAR_TRANSFORMATION_2
561 
562 template<class E1, class E2>
563 class matrix_addition: public blas::matrix_expression<matrix_addition<E1, E2> > {
564 private:
565  typedef scalar_binary_plus<
566  typename E1::value_type,
567  typename E2::value_type
568  > functor_type;
569 public:
570  typedef typename E1::const_closure_type lhs_closure_type;
571  typedef typename E2::const_closure_type rhs_closure_type;
572 
573  typedef typename E1::size_type size_type;
574  typedef typename E1::difference_type difference_type;
575  typedef typename functor_type::result_type value_type;
576  typedef value_type scalar_type;
577  typedef value_type const_reference;
578  typedef const_reference reference;
579  typedef value_type const* const_pointer;
580  typedef const_pointer pointer;
581 
582  typedef typename E1::index_type index_type;
583  typedef typename E1::const_index_pointer const_index_pointer;
585 
586  typedef const matrix_addition<E1, E2> const_closure_type;
587  typedef const_closure_type closure_type;
588  typedef typename E1::orientation orientation;
589  typedef blas::unknown_storage_tag storage_category;
590  typedef typename evaluation_restrict_traits<E1,E2>::type evaluation_category;
591 
592  // Construction
594  lhs_closure_type const& e1,
595  rhs_closure_type const& e2
596  ): m_lhs (e1), m_rhs (e2){}
597 
598  // Accessors
599  size_type size1 () const {
600  return m_lhs.size1 ();
601  }
602  size_type size2 () const {
603  return m_lhs.size2 ();
604  }
605 
606  const_reference operator () (index_type i, index_type j) const {
607  return m_lhs(i, j) + m_rhs(i,j);
608  }
609 
610  //computation kernels
611  template<class MatX>
612  void assign_to(matrix_expression<MatX>& X, scalar_type alpha = scalar_type(1) )const{
613  assign(X,alpha * m_lhs);
614  plus_assign(X,alpha * m_rhs);
615  }
616  template<class MatX>
617  void plus_assign_to(matrix_expression<MatX>& X, scalar_type alpha = scalar_type(1) )const{
618  plus_assign(X,alpha * m_lhs);
619  plus_assign(X,alpha * m_rhs);
620  }
621 
622  template<class MatX>
623  void minus_assign_to(matrix_expression<MatX>& X, scalar_type alpha = scalar_type(1) )const{
624  minus_assign(X,alpha * m_lhs);
625  minus_assign(X,alpha * m_rhs);
626  }
627 
628  // Iterator types
629 private:
630  typedef typename E1::const_row_iterator const_row_iterator1_type;
631  typedef typename E1::const_column_iterator const_row_column_iterator_type;
632  typedef typename E2::const_row_iterator const_column_iterator1_type;
633  typedef typename E2::const_column_iterator const_column_iterator2_type;
634 
635 public:
636  typedef binary_transform_iterator<
637  typename E1::const_row_iterator,
638  typename E2::const_row_iterator,
639  functor_type
641  typedef binary_transform_iterator<
642  typename E1::const_column_iterator,
643  typename E2::const_column_iterator,
644  functor_type
648 
649  const_row_iterator row_begin(std::size_t i) const {
650  return const_row_iterator (functor_type(),
651  m_lhs.row_begin(i),m_lhs.row_end(i),
652  m_rhs.row_begin(i),m_rhs.row_end(i)
653  );
654  }
655  const_row_iterator row_end(std::size_t i) const {
656  return const_row_iterator (functor_type(),
657  m_lhs.row_end(i),m_lhs.row_end(i),
658  m_rhs.row_end(i),m_rhs.row_end(i)
659  );
660  }
661 
662  const_column_iterator column_begin(std::size_t j) const {
663  return const_column_iterator (functor_type(),
664  m_lhs.column_begin(j),m_lhs.column_end(j),
665  m_rhs.column_begin(j),m_rhs.column_end(j)
666  );
667  }
668  const_column_iterator column_end(std::size_t j) const {
669  return const_column_iterator (functor_type(),
670  m_lhs.column_begin(j),m_lhs.column_end(j),
671  m_rhs.column_begin(j),m_rhs.column_end(j)
672  );
673  }
674 
675 private:
676  lhs_closure_type m_lhs;
677  rhs_closure_type m_rhs;
678  functor_type m_functor;
679 };
680 
681 ///\brief Adds two Matrices
682 template<class E1, class E2>
684  matrix_expression<E1> const& e1,
685  matrix_expression<E2> const& e2
686 ){
687  return matrix_addition<E1, E2>(e1(),e2());
688 }
689 
690 ///\brief Subtracts two Matrices
691 template<class E1, class E2>
693  matrix_expression<E1> const& e1,
694  matrix_expression<E2> const& e2
695 ){
697 }
698 
699 ///\brief Adds a matrix plus a scalar which is interpreted as a constant matrix
700 template<class E, class T>
701 typename boost::enable_if<
702  boost::is_convertible<T, typename E::value_type>,
704 >::type operator+ (
705  matrix_expression<E> const& e,
706  T t
707 ){
708  return e + scalar_matrix<T>(e().size1(),e().size2(),t);
709 }
710 
711 ///\brief Adds a matrix plus a scalar which is interpreted as a constant matrix
712 template<class T, class E>
713 typename boost::enable_if<
714  boost::is_convertible<T, typename E::value_type>,
716 >::type operator+ (
717  T t,
718  matrix_expression<E> const& e
719 ){
720  return e + scalar_matrix<T>(e().size1(),e().size2(),t);
721 }
722 
723 ///\brief Subtracts a scalar which is interpreted as a constant matrix from a matrix.
724 template<class E, class T>
725 typename boost::enable_if<
726  boost::is_convertible<T, typename E::value_type> ,
728 >::type operator- (
729  matrix_expression<E> const& e,
730  T t
731 ){
732  return e - scalar_matrix<T>(e().size1(),e().size2(),t);
733 }
734 
735 ///\brief Subtracts a matrix from a scalar which is interpreted as a constant matrix
736 template<class E, class T>
737 typename boost::enable_if<
738  boost::is_convertible<T, typename E::value_type>,
740 >::type operator- (
741  T t,
742  matrix_expression<E> const& e
743 ){
744  return scalar_matrix<T>(e().size1(),e().size2(),t) - e;
745 }
746 
747 template<class E1, class E2, class F>
749  public blas::matrix_expression<matrix_binary<E1, E2, F> > {
750 private:
751  typedef matrix_binary<E1, E2, F> self_type;
752 public:
753  typedef E1 lhs_type;
754  typedef E2 rhs_type;
755  typedef typename E1::const_closure_type lhs_closure_type;
756  typedef typename E2::const_closure_type rhs_closure_type;
757 
758  typedef typename E1::size_type size_type;
759  typedef typename E1::difference_type difference_type;
760  typedef typename F::result_type value_type;
761  typedef value_type scalar_type;
762  typedef value_type const_reference;
763  typedef const_reference reference;
764  typedef value_type const* const_pointer;
765  typedef const_pointer pointer;
766 
767  typedef typename E1::index_type index_type;
768  typedef typename E1::const_index_pointer const_index_pointer;
770 
771  typedef const matrix_binary<E1, E2, F> const_closure_type;
772  typedef const_closure_type closure_type;
773  typedef typename E1::orientation orientation;
774  typedef blas::unknown_storage_tag storage_category;
775  typedef typename evaluation_restrict_traits<E1,E2>::type evaluation_category;
776 
777  typedef F functor_type;
778 
779  // Construction and destruction
780 
782  matrix_expression<E1> const&e1, matrix_expression<E2> const& e2, functor_type functor
783  ): m_lhs (e1()), m_rhs (e2()),m_functor(functor) {}
784 
785  // Accessors
786  size_type size1 () const {
787  return m_lhs.size1 ();
788  }
789  size_type size2 () const {
790  return m_lhs.size2 ();
791  }
792 
793  const_reference operator () (index_type i, index_type j) const {
794  return m_functor( m_lhs (i, j), m_rhs(i,j));
795  }
796 
797  // Iterator types
798 private:
799  typedef typename E1::const_row_iterator const_row_iterator1_type;
800  typedef typename E1::const_column_iterator const_row_column_iterator_type;
801  typedef typename E2::const_row_iterator const_column_iterator1_type;
802  typedef typename E2::const_column_iterator const_column_iterator2_type;
803 
804 public:
805  typedef binary_transform_iterator<
806  typename E1::const_row_iterator,typename E2::const_row_iterator,functor_type
808  typedef binary_transform_iterator<
809  typename E1::const_column_iterator,typename E2::const_column_iterator,functor_type
813 
814  const_row_iterator row_begin(std::size_t i) const {
815  return const_row_iterator (m_functor,
816  m_lhs.row_begin(i),m_lhs.row_end(i),
817  m_rhs.row_begin(i),m_rhs.row_end(i)
818  );
819  }
820  const_row_iterator row_end(std::size_t i) const {
821  return const_row_iterator (m_functor,
822  m_lhs.row_end(i),m_lhs.row_end(i),
823  m_rhs.row_end(i),m_rhs.row_end(i)
824  );
825  }
826 
827  const_column_iterator column_begin(std::size_t j) const {
828  return const_column_iterator (m_functor,
829  m_lhs.column_begin(j),m_lhs.column_end(j),
830  m_rhs.column_begin(j),m_rhs.column_end(j)
831  );
832  }
833  const_column_iterator column_end(std::size_t j) const {
834  return const_column_iterator (m_functor,
835  m_lhs.column_begin(j),m_lhs.column_end(j),
836  m_rhs.column_begin(j),m_rhs.column_end(j)
837  );
838  }
839 
840 private:
841  lhs_closure_type m_lhs;
842  rhs_closure_type m_rhs;
843  functor_type m_functor;
844 };
845 
846 #define SHARK_BINARY_MATRIX_EXPRESSION(name, F)\
847 template<class E1, class E2>\
848 matrix_binary<E1, E2, F<typename E1::value_type, typename E2::value_type> >\
849 name(matrix_expression<E1> const& e1, matrix_expression<E2> const& e2){\
850  typedef F<typename E1::value_type, typename E2::value_type> functor_type;\
851  return matrix_binary<E1, E2, functor_type>(e1,e2, functor_type());\
852 }
853 SHARK_BINARY_MATRIX_EXPRESSION(operator*, scalar_binary_multiply)
854 SHARK_BINARY_MATRIX_EXPRESSION(element_prod, scalar_binary_multiply)
855 SHARK_BINARY_MATRIX_EXPRESSION(operator/, scalar_binary_divide)
856 SHARK_BINARY_MATRIX_EXPRESSION(pow,scalar_binary_pow)
857 SHARK_BINARY_MATRIX_EXPRESSION(element_div, scalar_binary_divide)
858 #undef SHARK_BINARY_MATRIX_EXPRESSION
859 
860 template<class E1, class E2>
861 matrix_binary<E1, E2,
862  scalar_binary_safe_divide<typename E1::value_type, typename E2::value_type>
863 >
865  matrix_expression<E1> const& e1,
866  matrix_expression<E2> const& e2,
867  typename promote_traits<
868  typename E1::value_type,
869  typename E2::value_type
870  >::promote_type defaultValue
871 ){
872  typedef scalar_binary_safe_divide<typename E1::value_type, typename E2::value_type> functor_type;
873  return matrix_binary<E1, E2, functor_type>(e1,e2, functor_type(defaultValue));
874 }
875 
876 template<class MatA, class VecV>
878  public vector_expression<matrix_vector_prod<MatA, VecV> > {
879 public:
880  typedef typename MatA::const_closure_type matrix_closure_type;
881  typedef typename VecV::const_closure_type vector_closure_type;
882 public:
883  typedef typename promote_traits<
884  typename MatA::scalar_type,
885  typename VecV::scalar_type
886  >::promote_type scalar_type;
887  typedef scalar_type value_type;
888  typedef typename MatA::size_type size_type;
889  typedef typename MatA::difference_type difference_type;
890  typedef typename MatA::index_type index_type;
891  typedef typename MatA::index_pointer index_pointer;
892  typedef typename MatA::const_index_pointer const_index_pointer;
893 
894  typedef matrix_vector_prod<MatA, VecV> const_closure_type;
895  typedef const_closure_type closure_type;
896  typedef unknown_storage_tag storage_category;
897  typedef blockwise_tag evaluation_category;
898 
899 
900  //FIXME: This workaround is required to be able to generate
901  // temporary vectors
902  typedef typename MatA::const_row_iterator const_iterator;
903  typedef const_iterator iterator;
904 
905  // Construction and destruction
907  matrix_closure_type const& matrix,
908  vector_closure_type const& vector
909  ):m_matrix(matrix), m_vector(vector) {}
910 
911  size_type size() const {
912  return m_matrix.size1();
913  }
914 
915  matrix_closure_type const& matrix() const {
916  return m_matrix;
917  }
918  vector_closure_type const& vector() const {
919  return m_vector;
920  }
921 
922  //computation kernels
923  template<class VecX>
924  void assign_to(vector_expression<VecX>& x, scalar_type alpha = scalar_type(1) )const{
925  x().clear();
926  plus_assign_to(x,alpha);
927  }
928  template<class VecX>
929  void plus_assign_to(vector_expression<VecX>& x, scalar_type alpha = scalar_type(1) )const{
930  kernels::gemv(eval_block(m_matrix), eval_block(m_vector), x, alpha);
931  }
932 
933  template<class VecX>
934  void minus_assign_to(vector_expression<VecX>& x, scalar_type alpha = scalar_type(1) )const{
935  kernels::gemv(eval_block(m_matrix), eval_block(m_vector), x, -alpha);
936  }
937 
938 private:
939  matrix_closure_type m_matrix;
940  vector_closure_type m_vector;
941 };
942 
943 
944 /// \brief computes the matrix-vector product x+=Av
945 template<class MatA, class VecV>
947  return matrix_vector_prod<MatA,VecV>(A(),v());
948 }
949 
950 /// \brief computes the matrix-vector product x+=v^TA
951 template<class MatA, class VecV>
953  //compute it using the identity (v^TA)^T= A^Tv
955 }
956 
957 //matrix-matrix prod
958 template<class MatA, class MatB>
959 class matrix_matrix_prod: public matrix_expression<matrix_matrix_prod<MatA, MatB> > {
960 public:
961  typedef typename MatA::const_closure_type matrix_closure_typeA;
962  typedef typename MatB::const_closure_type matrix_closure_typeB;
963 public:
964  typedef typename promote_traits<
965  typename MatA::scalar_type,
966  typename MatB::scalar_type
967  >::promote_type scalar_type;
968  typedef scalar_type value_type;
969  typedef typename MatA::size_type size_type;
970  typedef typename MatA::difference_type difference_type;
971  typedef typename MatA::index_type index_type;
972  typedef typename MatA::index_pointer index_pointer;
973  typedef typename MatA::const_index_pointer const_index_pointer;
974 
975  typedef matrix_matrix_prod<MatA, MatB> const_closure_type;
976  typedef const_closure_type closure_type;
977  typedef unknown_storage_tag storage_category;
978  typedef blockwise_tag evaluation_category;
979  typedef unknown_orientation orientation;
980 
981 
982  //FIXME: This workaround is required to be able to generate
983  // temporary matrices
984  typedef typename MatA::const_row_iterator const_row_iterator;
985  typedef typename MatA::const_column_iterator const_column_iterator;
986  typedef const_row_iterator row_iterator;
987  typedef const_column_iterator column_iterator;
988 
989  // Construction and destruction
991  matrix_closure_typeA const& matrixA,
992  matrix_closure_typeB const& matrixB
993  ):m_matrixA(matrixA), m_matrixB(matrixB) {}
994 
995  size_type size1() const {
996  return m_matrixA.size1();
997  }
998  size_type size2() const {
999  return m_matrixB.size2();
1000  }
1001 
1002  matrix_closure_typeA const& matrixA() const {
1003  return m_matrixA;
1004  }
1005  matrix_closure_typeB const& matrixB() const {
1006  return m_matrixB;
1007  }
1008 
1009  //computation kernels
1010  template<class MatX>
1011  void assign_to(matrix_expression<MatX>& X, scalar_type alpha = scalar_type(1) )const{
1012  X().clear();
1013  plus_assign_to(X,alpha);
1014  }
1015  template<class MatX>
1016  void plus_assign_to(matrix_expression<MatX>& X, scalar_type alpha = scalar_type(1) )const{
1017  kernels::gemm(eval_block(m_matrixA), eval_block(m_matrixB), X, alpha);
1018  }
1019 
1020  template<class MatX>
1021  void minus_assign_to(matrix_expression<MatX>& X, scalar_type alpha = scalar_type(1) )const{
1022  kernels::gemm(eval_block(m_matrixA), eval_block(m_matrixB), X, -alpha);
1023  }
1024 
1025 private:
1026  matrix_closure_typeA m_matrixA;
1027  matrix_closure_typeB m_matrixB;
1028 };
1029 
1030 /// \brief computes the matrix-matrix product X+=AB
1031 template<class MatA, class MatB>
1033  matrix_expression<MatA> const& A,
1034  matrix_expression<MatB> const& B
1035 ) {
1036  return matrix_matrix_prod<MatA,MatB>(A(),B());
1037 }
1038 
1039 namespace detail{
1040 
1041 template<class MatA,class VecB>
1042 void sum_rows_impl(MatA const& matA, VecB& vecB, column_major){
1043  for(std::size_t i = 0; i != matA.size2(); ++i){
1044  vecB(i)+=sum(column(matA,i));
1045  }
1046 }
1047 
1048 template<class MatA,class VecB>
1049 void sum_rows_impl(MatA const& matA, VecB& vecB, row_major){
1050  for(std::size_t i = 0; i != matA.size1(); ++i){
1051  noalias(vecB) += row(matA,i);
1052  }
1053 }
1054 
1055 template<class MatA,class VecB>
1056 void sum_rows_impl(MatA const& matA, VecB& vecB, unknown_orientation){
1057  sum_rows_impl(matA,vecB,row_major());
1058 }
1059 
1060 //dispatcher for triangular matrix
1061 template<class MatA,class VecB,class Orientation,class Triangular>
1062 void sum_rows_impl(MatA const& matA, VecB& vecB, packed<Orientation,Triangular>){
1063  sum_rows_impl(matA,vecB,Orientation());
1064 }
1065 
1066 
1067 //dispatcher
1068 template<class MatA,class VecB>
1069 void sum_rows_impl(MatA const& matA, VecB& vecB){
1070  sum_rows_impl(matA,vecB,typename MatA::orientation());
1071 }
1072 
1073 template<class MatA>
1074 typename MatA::value_type sum_impl(MatA const& matA, column_major){
1075  typename MatA::value_type totalSum = 0;
1076  for(std::size_t j = 0; j != matA.size2(); ++j){
1077  totalSum += sum(column(matA,j));
1078  }
1079  return totalSum;
1080 }
1081 
1082 template<class MatA>
1083 typename MatA::value_type sum_impl(MatA const& matA, row_major){
1084  typename MatA::value_type totalSum = 0;
1085  for(std::size_t i = 0; i != matA.size1(); ++i){
1086  totalSum += sum(row(matA,i));
1087  }
1088  return totalSum;
1089 }
1090 
1091 template<class MatA>
1092 typename MatA::value_type sum_impl(MatA const& matA, unknown_orientation){
1093  return sum_impl(matA,row_major());
1094 }
1095 
1096 
1097 //dispatcher for triangular matrix
1098 template<class MatA,class Orientation,class Triangular>
1099 typename MatA::value_type sum_impl(MatA const& matA, packed<Orientation,Triangular>){
1100  return sum_impl(matA,Orientation());
1101 }
1102 
1103 //dispatcher
1104 template<class MatA>
1105 typename MatA::value_type sum_impl(MatA const& matA){
1106  return sum_impl(matA,typename MatA::orientation());
1107 }
1108 
1109 template<class MatA>
1110 typename MatA::value_type max_impl(MatA const& matA, column_major){
1111  typename MatA::value_type maximum = 0;
1112  for(std::size_t j = 0; j != matA.size2(); ++j){
1113  maximum = std::max(maximum, max(column(matA,j)));
1114  }
1115  return maximum;
1116 }
1117 
1118 template<class MatA>
1119 typename MatA::value_type max_impl(MatA const& matA, row_major){
1120  typename MatA::value_type maximum = 0;
1121  for(std::size_t i = 0; i != matA.size1(); ++i){
1122  maximum= std::max(maximum, max(row(matA,i)));
1123  }
1124  return maximum;
1125 }
1126 
1127 template<class MatA>
1128 typename MatA::value_type max_impl(MatA const& matA, unknown_orientation){
1129  return max_impl(matA,row_major());
1130 }
1131 
1132 //dispatcher for triangular matrix
1133 template<class MatA,class Orientation,class Triangular>
1134 typename MatA::value_type max_impl(MatA const& matA, packed<Orientation,Triangular>){
1135  return std::max(max_impl(matA,Orientation()),0.0);
1136 }
1137 
1138 //dispatcher
1139 template<class MatA>
1140 typename MatA::value_type max_impl(MatA const& matA){
1141  return max_impl(matA,typename MatA::orientation());
1142 }
1143 
1144 template<class MatA>
1145 typename MatA::value_type min_impl(MatA const& matA, column_major){
1146  typename MatA::value_type minimum = 0;
1147  for(std::size_t j = 0; j != matA.size2(); ++j){
1148  minimum= std::min(minimum, min(column(matA,j)));
1149  }
1150  return minimum;
1151 }
1152 
1153 template<class MatA>
1154 typename MatA::value_type min_impl(MatA const& matA, row_major){
1155  typename MatA::value_type minimum = 0;
1156  for(std::size_t i = 0; i != matA.size1(); ++i){
1157  minimum= std::min(minimum, min(row(matA,i)));
1158  }
1159  return minimum;
1160 }
1161 
1162 template<class MatA>
1163 typename MatA::value_type min_impl(MatA const& matA, unknown_orientation){
1164  return min_impl(matA,row_major());
1165 }
1166 
1167 //dispatcher for triangular matrix
1168 template<class MatA,class Orientation,class Triangular>
1169 typename MatA::value_type min_impl(MatA const& matA, packed<Orientation,Triangular>){
1170  return std::min(min_impl(matA,Orientation()),0.0);
1171 }
1172 
1173 //dispatcher
1174 template<class MatA>
1175 typename MatA::value_type min_impl(MatA const& matA){
1176  return min_impl(matA,typename MatA::orientation());
1177 }
1178 
1179 }//end detail
1180 
1181 //dispatcher
1182 template<class MatA>
1183 typename vector_temporary_type<
1184  typename MatA::value_type,
1185  dense_random_access_iterator_tag
1186 >::type
1188  typename vector_temporary_type<
1189  typename MatA::value_type,
1190  dense_random_access_iterator_tag
1191  >::type result(A().size2(),0.0);
1192  detail::sum_rows_impl(eval_block(A),result);
1193  return result;
1194 }
1195 
1196 template<class MatA>
1197 typename vector_temporary_type<
1198  typename MatA::value_type,
1199  dense_random_access_iterator_tag
1200 >::type
1202  //implemented using sum_rows_impl by transposing A
1203  typename vector_temporary_type<
1204  typename MatA::value_type,
1205  dense_random_access_iterator_tag
1206  >::type result(A().size1(),0);
1207  detail::sum_rows_impl(eval_block(trans(A)),result);
1208  return result;
1209 }
1210 
1211 
1212 template<class MatA>
1213 typename MatA::value_type sum(matrix_expression<MatA> const& A){
1214  return detail::sum_impl(eval_block(A));
1215 }
1216 
1217 template<class MatA>
1218 typename MatA::value_type max(matrix_expression<MatA> const& A){
1219  return detail::max_impl(eval_block(A));
1220 }
1221 
1222 template<class MatA>
1223 typename MatA::value_type min(matrix_expression<MatA> const& A){
1224  return detail::min_impl(eval_block(A));
1225 }
1226 
1227 /// \brief Returns the frobenius inner-product between matrices exprssions 1 and e2.
1228 ///
1229 ///The frobenius inner product is defined as \f$ <A,B>_F=\sum_{ij} A_ij*B_{ij} \f$. It induces the
1230 /// Frobenius norm \f$ ||A||_F = \sqrt{<A,A>_F} \f$
1231 template<class E1, class E2>
1232 typename promote_traits <typename E1::value_type,typename E2::value_type>::promote_type
1234  matrix_expression<E1> const& e1,
1235  matrix_expression<E2> const& e2
1236 ) {
1237  return sum(eval_block(e1)*eval_block(e2));
1238 }
1239 
1240 
1241 template<class E>
1242 typename matrix_norm_1<E>::result_type
1244  return matrix_norm_1<E>::apply(eval_block(e));
1245 }
1246 
1247 template<class E>
1248 typename real_traits<typename E::value_type>::type
1250  using std::sqrt;
1251  return sqrt(sum(abs_sqr(eval_block(e))));
1252 }
1253 
1254 template<class E>
1255 typename matrix_norm_inf<E>::result_type
1257  return matrix_norm_inf<E>::apply(eval_block(e));
1258 }
1259 
1260 /*!
1261  * \brief Evaluates the sum of the values at the diagonal of
1262  * matrix "v".
1263  *
1264  * Example:
1265  * \f[
1266  * \left(
1267  * \begin{array}{*{4}{c}}
1268  * {\bf 1} & 5 & 9 & 13\\
1269  * 2 & {\bf 6} & 10 & 14\\
1270  * 3 & 7 & {\bf 11} & 15\\
1271  * 4 & 8 & 12 & {\bf 16}\\
1272  * \end{array}
1273  * \right)
1274  * \longrightarrow 1 + 6 + 11 + 16 = 34
1275  * \f]
1276  *
1277  * \param m square matrix
1278  * \return the sum of the values at the diagonal of \em m
1279  */
1280 template < class MatrixT >
1281 typename MatrixT::value_type trace(matrix_expression<MatrixT> const& m)
1282 {
1283  SIZE_CHECK(m().size1() == m().size2());
1284 
1285  typename MatrixT::value_type t(m()(0, 0));
1286  for (unsigned i = 1; i < m().size1(); ++i)
1287  t += m()(i, i);
1288  return t;
1289 }
1290 
1291 /** \brief An identity matrix with values of type \c T
1292  *
1293  * Elements or cordinates \f$(i,i)\f$ are equal to 1 (one) and all others to 0 (zero).
1294  */
1295 template<class T>
1296 class identity_matrix: public diagonal_matrix<scalar_vector<T> > {
1298 public:
1300  identity_matrix(std::size_t size):base_type(scalar_vector<T>(size,T(1))){}
1301 };
1302 
1303 }
1304 }
1305 
1306 #endif