gemv.hpp
Go to the documentation of this file.
1 /*!
2  *
3  *
4  * \brief -
5  *
6  * \author O. Krause
7  * \date 2012
8  *
9  *
10  * \par Copyright 1995-2015 Shark Development Team
11  *
12  * <BR><HR>
13  * This file is part of Shark.
14  * <http://image.diku.dk/shark/>
15  *
16  * Shark is free software: you can redistribute it and/or modify
17  * it under the terms of the GNU Lesser General Public License as published
18  * by the Free Software Foundation, either version 3 of the License, or
19  * (at your option) any later version.
20  *
21  * Shark is distributed in the hope that it will be useful,
22  * but WITHOUT ANY WARRANTY; without even the implied warranty of
23  * MatAERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24  * GNU Lesser General Public License for more details.
25  *
26  * You should have received a copy of the GNU Lesser General Public License
27  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
28  *
29  */
30 #ifndef REMORA_KERNELS_DEFAULT_GEMatAV_HPP
31 #define REMORA_KERNELS_DEFAULT_GEMatAV_HPP
32 
33 #include "../../expression_types.hpp" //matrix/vector_expression
34 #include "../../detail/matrix_proxy_classes.hpp" //matrix_row, matrix_transpose
35 #include "../../detail/traits.hpp" //matrix orientations
36 #include "../dot.hpp" //inner product
37 #include "../../assignment.hpp" //plus_assign
38 #include <type_traits> //std::false_type marker for unoptimized
39 
40 namespace remora{namespace bindings {
41 
42 //row major can be further reduced to inner_prod()
43 template<class ResultV, class MatA, class V>
44 void gemv_impl(
45  matrix_expression<MatA, cpu_tag> const& A,
46  vector_expression<V, cpu_tag> const& x,
47  vector_expression<ResultV, cpu_tag>& result,
48  typename ResultV::value_type alpha,
49  row_major
50 ) {
51  typedef typename ResultV::value_type value_type;
52  value_type value;
53  for(std::size_t i = 0; i != A().size1();++i){
54  matrix_row<typename const_expression<MatA>::type > rowA(A(),i);
55  kernels::dot(rowA,x,value);
56  if(value != value_type())//handling of sparse results.
57  result()(i) += alpha * value;
58  }
59 }
60 
61 //column major is implemented by computing a linear combination of matrix-rows
62 template<class ResultV, class MatA, class V>
63 void gemv_impl(
64  matrix_expression<MatA, cpu_tag> const& A,
65  vector_expression<V, cpu_tag> const& x,
66  vector_expression<ResultV, cpu_tag>& result,
67  typename ResultV::value_type alpha,
68  column_major
69 ) {
70  //instead of a matrix column, we have matrix_row
71  typedef matrix_transpose<typename const_expression<MatA>::type > TransA;
72  TransA transA(A());
73  typedef typename V::const_iterator iterator;
74  typedef typename ResultV::value_type value_type;
75  iterator end = x().end();
76  for(iterator it = x().begin(); it != end; ++it) {
77  value_type multiplier = alpha * (*it);
78  matrix_row<TransA> colA(transA,it.index());
79  //FIXME: for sparse result vectors, this might hurt.
80  plus_assign(result,colA,multiplier);
81  }
82 }
83 
84 //unknown orientation is dispatched to row_major
85 template<class ResultV, class MatA, class V>
86 void gemv_impl(
87  matrix_expression<MatA, cpu_tag> const& A,
88  vector_expression<V, cpu_tag> const& x,
89  vector_expression<ResultV, cpu_tag>& result,
90  typename ResultV::value_type alpha,
91  unknown_orientation
92 ) {
93  gemv_impl(A,x,result,alpha,row_major());
94 }
95 
96 // result += alpha * A * x
97 template<class ResultV, class MatA, class V>
98 void gemv(
99  matrix_expression<MatA, cpu_tag> const& A,
100  vector_expression<V, cpu_tag> const& x,
101  vector_expression<ResultV, cpu_tag>& result,
102  typename ResultV::value_type alpha,
103  std::false_type
104 ) {
105  typedef typename MatA::orientation orientation;
106 
107  gemv_impl(A, x, result, alpha, orientation());
108 }
109 
110 }}
111 #endif