28 #ifndef SHARK_LINALG_BLAS_MATRIX_EXPRESSION_HPP 29 #define SHARK_LINALG_BLAS_MATRIX_EXPRESSION_HPP 31 #include <boost/type_traits/is_convertible.hpp> 38 template<
class E1,
class E2>
40 typedef scalar_multiply1<typename E1::value_type, typename E2::value_type> functor_type;
42 typedef typename E1::const_closure_type lhs_closure_type;
43 typedef typename E2::const_closure_type rhs_closure_type;
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;
55 typedef typename E1::index_type index_type;
68 :m_lhs(e1), m_rhs(e2){}
79 const lhs_closure_type &
lhs()
const {
82 const rhs_closure_type &
rhs()
const {
86 const_reference
operator()(index_type i, index_type j)
const {
87 return m_lhs(i) * m_rhs(j);
97 functor_type(m_lhs(i))
100 const_row_iterator
row_end(index_type i)
const {
102 functor_type(m_lhs(i))
108 functor_type(m_rhs(i))
113 functor_type(m_rhs(i))
117 lhs_closure_type m_lhs;
118 rhs_closure_type m_rhs;
122 template<
class E1,
class E2>
134 typedef V expression_type;
136 typedef typename V::const_iterator const_subiterator_type;
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;
148 typedef typename V::index_type index_type;
152 typedef self_type
const const_closure_type;
153 typedef const_closure_type closure_type;
160 m_vector(e), m_rows(rows) {}
167 return m_vector.size();
176 const_reference
operator() (index_type i, index_type j)
const {
190 return m_vector.begin();
192 const_row_iterator
row_end(std::size_t i)
const {
194 return m_vector.end();
206 expression_closure_type m_vector;
219 template<
class Vector>
234 typedef std::size_t size_type;
235 typedef std::ptrdiff_t difference_type;
236 typedef T value_type;
240 typedef value_type scalar_type;
243 typedef std::size_t index_type;
255 m_size1(0), m_size2(0), m_value() {}
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) {}
283 const_row_iterator
row_end(std::size_t i)
const {
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;
317 typedef typename E::const_closure_type expression_closure_type;
319 typedef typename functor_type::result_type value_type;
320 typedef typename E::scalar_type scalar_type;
325 typedef typename E::size_type size_type;
326 typedef typename E::difference_type difference_type;
328 typedef typename E::index_type index_type;
333 typedef const_closure_type closure_type;
340 m_expression(e()), m_scalar(scalar){}
344 return m_expression.size1();
347 return m_expression.size2();
351 const_reference
operator()(index_type i, index_type j)
const {
352 return m_scalar * m_expression(i, j);
358 m_expression.assign_to(X,alpha*m_scalar);
362 m_expression.plus_assign_to(X,alpha*m_scalar);
367 m_expression.minus_assign_to(X,alpha*m_scalar);
379 const_row_iterator
row_end(index_type i)
const {
391 expression_closure_type m_expression;
392 scalar_type m_scalar;
404 template<
class E,
class F>
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;
413 typedef typename expression_type::const_closure_type expression_closure_type;
415 typedef F functor_type;
416 typedef typename functor_type::result_type value_type;
417 typedef value_type scalar_type;
422 typedef typename expression_type::size_type size_type;
423 typedef typename expression_type::difference_type difference_type;
425 typedef typename E::index_type index_type;
429 typedef self_type
const const_closure_type;
430 typedef const_closure_type closure_type;
437 m_expression(e()), m_functor(functor) {}
441 return m_expression.size1();
444 return m_expression.size2();
448 const_reference
operator()(index_type i, index_type j)
const {
449 return m_functor(m_expression(i, j));
461 const_row_iterator
row_end(index_type i)
const {
473 expression_closure_type m_expression;
474 functor_type m_functor;
477 template<
class E,
class T>
478 typename boost::enable_if<
479 boost::is_convertible<T, typename E::scalar_type >,
486 template<
class T,
class E>
487 typename boost::enable_if<
488 boost::is_convertible<T, typename E::scalar_type >,
500 #define SHARK_UNARY_MATRIX_TRANSFORMATION(name, F)\ 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());\ 523 #undef SHARK_UNARY_MATRIX_TRANSFORMATION 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> > \ 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)); \ 545 #undef SHARK_MATRIX_SCALAR_TRANSFORMATION 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> > \ 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)); \ 560 #undef SHARK_MATRIX_SCALAR_TRANSFORMATION_2 562 template<
class E1,
class E2>
565 typedef scalar_binary_plus<
566 typename E1::value_type,
567 typename E2::value_type
570 typedef typename E1::const_closure_type lhs_closure_type;
571 typedef typename E2::const_closure_type rhs_closure_type;
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;
582 typedef typename E1::index_type index_type;
587 typedef const_closure_type closure_type;
594 lhs_closure_type
const& e1,
595 rhs_closure_type
const& e2
596 ): m_lhs (e1), m_rhs (e2){}
600 return m_lhs.size1 ();
603 return m_lhs.size2 ();
607 return m_lhs(i, j) + m_rhs(i,j);
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;
636 typedef binary_transform_iterator<
637 typename E1::const_row_iterator,
638 typename E2::const_row_iterator,
641 typedef binary_transform_iterator<
642 typename E1::const_column_iterator,
643 typename E2::const_column_iterator,
651 m_lhs.row_begin(i),m_lhs.row_end(i),
652 m_rhs.row_begin(i),m_rhs.row_end(i)
657 m_lhs.row_end(i),m_lhs.row_end(i),
658 m_rhs.row_end(i),m_rhs.row_end(i)
664 m_lhs.column_begin(j),m_lhs.column_end(j),
665 m_rhs.column_begin(j),m_rhs.column_end(j)
670 m_lhs.column_begin(j),m_lhs.column_end(j),
671 m_rhs.column_begin(j),m_rhs.column_end(j)
676 lhs_closure_type m_lhs;
677 rhs_closure_type m_rhs;
678 functor_type m_functor;
682 template<
class E1,
class E2>
691 template<
class E1,
class E2>
700 template<
class E,
class T>
701 typename boost::enable_if<
702 boost::is_convertible<T, typename E::value_type>,
712 template<
class T,
class E>
713 typename boost::enable_if<
714 boost::is_convertible<T, typename E::value_type>,
724 template<
class E,
class T>
725 typename boost::enable_if<
726 boost::is_convertible<T, typename E::value_type> ,
736 template<
class E,
class T>
737 typename boost::enable_if<
738 boost::is_convertible<T, typename E::value_type>,
747 template<
class E1,
class E2,
class F>
755 typedef typename E1::const_closure_type lhs_closure_type;
756 typedef typename E2::const_closure_type rhs_closure_type;
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;
767 typedef typename E1::index_type index_type;
772 typedef const_closure_type closure_type;
777 typedef F functor_type;
783 ): m_lhs (e1()), m_rhs (e2()),m_functor(functor) {}
787 return m_lhs.size1 ();
790 return m_lhs.size2 ();
794 return m_functor( m_lhs (i, j), m_rhs(i,j));
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;
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
816 m_lhs.row_begin(i),m_lhs.row_end(i),
817 m_rhs.row_begin(i),m_rhs.row_end(i)
822 m_lhs.row_end(i),m_lhs.row_end(i),
823 m_rhs.row_end(i),m_rhs.row_end(i)
829 m_lhs.column_begin(j),m_lhs.column_end(j),
830 m_rhs.column_begin(j),m_rhs.column_end(j)
835 m_lhs.column_begin(j),m_lhs.column_end(j),
836 m_rhs.column_begin(j),m_rhs.column_end(j)
841 lhs_closure_type m_lhs;
842 rhs_closure_type m_rhs;
843 functor_type m_functor;
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());\ 858 #undef SHARK_BINARY_MATRIX_EXPRESSION 860 template<
class E1,
class E2>
861 matrix_binary<E1, E2,
862 scalar_binary_safe_divide<typename E1::value_type, typename E2::value_type>
867 typename promote_traits<
868 typename E1::value_type,
869 typename E2::value_type
870 >::promote_type defaultValue
872 typedef scalar_binary_safe_divide<typename E1::value_type, typename E2::value_type> functor_type;
876 template<
class MatA,
class VecV>
880 typedef typename MatA::const_closure_type matrix_closure_type;
881 typedef typename VecV::const_closure_type vector_closure_type;
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;
895 typedef const_closure_type closure_type;
907 matrix_closure_type
const&
matrix,
908 vector_closure_type
const&
vector 909 ):m_matrix(matrix), m_vector(vector) {}
912 return m_matrix.size1();
915 matrix_closure_type
const&
matrix()
const {
918 vector_closure_type
const&
vector()
const {
926 plus_assign_to(x,alpha);
939 matrix_closure_type m_matrix;
940 vector_closure_type m_vector;
945 template<
class MatA,
class VecV>
951 template<
class MatA,
class VecV>
958 template<
class MatA,
class MatB>
961 typedef typename MatA::const_closure_type matrix_closure_typeA;
962 typedef typename MatB::const_closure_type matrix_closure_typeB;
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;
976 typedef const_closure_type closure_type;
991 matrix_closure_typeA
const& matrixA,
992 matrix_closure_typeB
const& matrixB
993 ):m_matrixA(matrixA), m_matrixB(matrixB) {}
996 return m_matrixA.size1();
999 return m_matrixB.size2();
1010 template<
class MatX>
1013 plus_assign_to(X,alpha);
1015 template<
class MatX>
1020 template<
class MatX>
1026 matrix_closure_typeA m_matrixA;
1027 matrix_closure_typeB m_matrixB;
1031 template<
class MatA,
class MatB>
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){
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){
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());
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());
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());
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){
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));
1091 template<
class MatA>
1092 typename MatA::value_type sum_impl(MatA
const& matA, unknown_orientation){
1093 return sum_impl(matA,row_major());
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());
1104 template<
class MatA>
1105 typename MatA::value_type sum_impl(MatA
const& matA){
1106 return sum_impl(matA,
typename MatA::orientation());
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){
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){
1127 template<
class MatA>
1128 typename MatA::value_type max_impl(MatA
const& matA, unknown_orientation){
1129 return max_impl(matA,row_major());
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);
1139 template<
class MatA>
1140 typename MatA::value_type max_impl(MatA
const& matA){
1141 return max_impl(matA,
typename MatA::orientation());
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){
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){
1162 template<
class MatA>
1163 typename MatA::value_type min_impl(MatA
const& matA, unknown_orientation){
1164 return min_impl(matA,row_major());
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);
1174 template<
class MatA>
1175 typename MatA::value_type min_impl(MatA
const& matA){
1176 return min_impl(matA,
typename MatA::orientation());
1182 template<
class MatA>
1183 typename vector_temporary_type<
1184 typename MatA::value_type,
1185 dense_random_access_iterator_tag
1188 typename vector_temporary_type<
1189 typename MatA::value_type,
1190 dense_random_access_iterator_tag
1191 >::type result(A().
size2(),0.0);
1196 template<
class MatA>
1197 typename vector_temporary_type<
1198 typename MatA::value_type,
1199 dense_random_access_iterator_tag
1203 typename vector_temporary_type<
1204 typename MatA::value_type,
1205 dense_random_access_iterator_tag
1206 >::type result(A().
size1(),0);
1212 template<
class MatA>
1217 template<
class MatA>
1222 template<
class MatA>
1231 template<
class E1,
class E2>
1232 typename promote_traits <typename E1::value_type,typename E2::value_type>::promote_type
1242 typename matrix_norm_1<E>::result_type
1244 return matrix_norm_1<E>::apply(
eval_block(e));
1248 typename real_traits<typename E::value_type>::type
1255 typename matrix_norm_inf<E>::result_type
1257 return matrix_norm_inf<E>::apply(
eval_block(e));
1280 template <
class MatrixT >
1285 typename MatrixT::value_type t(m()(0, 0));
1286 for (
unsigned i = 1; i < m().size1(); ++i)