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  * MERCHANTABILITY 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 SHARK_LINALG_BLAS_KERNELS_DEFAULT_GEMV_HPP
31 #define SHARK_LINALG_BLAS_KERNELS_DEFAULT_GEMV_HPP
32 
33 #include "../../matrix_proxy.hpp"
34 #include "../../vector_expression.hpp"
35 #include <boost/mpl/bool.hpp>
36 
37 namespace shark {namespace blas {namespace bindings {
38 
39 //row major can be further reduced to inner_prod()
40 template<class ResultV, class M, class V>
41 void gemv_impl(
42  matrix_expression<M> const& A,
43  vector_expression<V> const& x,
44  vector_expression<ResultV>& result,
45  typename ResultV::value_type alpha,
46  row_major
47 ) {
48  typedef typename ResultV::value_type value_type;
49  for(std::size_t i = 0; i != A().size1();++i){
50  value_type value = inner_prod(row(A,i),x);
51  if(value != value_type())//handling of sparse results.
52  result()(i) += alpha* value;
53  }
54 }
55 
56 //column major is implemented by computing a linear combination of matrix-rows
57 template<class ResultV, class M, class V>
58 void gemv_impl(
59  matrix_expression<M> const& A,
60  vector_expression<V> const& x,
61  vector_expression<ResultV>& result,
62  typename ResultV::value_type alpha,
63  column_major
64 ) {
65  typedef typename V::const_iterator iterator;
66  typedef typename ResultV::value_type value_type;
67  iterator end = x().end();
68  for(iterator it = x().begin(); it != end; ++it) {
69  value_type multiplier = alpha * (*it);
70  //fixme: for sparse result vectors, this might hurt.
71  noalias(result)+= multiplier * column(A(), it.index());
72  }
73 }
74 
75 //unknown orientation is dispatched to row_major
76 template<class ResultV, class M, class V>
77 void gemv_impl(
78  matrix_expression<M> const& A,
79  vector_expression<V> const& x,
80  vector_expression<ResultV>& result,
81  typename ResultV::value_type alpha,
82  unknown_orientation
83 ) {
84  gemv_impl(A,x,result,alpha,row_major());
85 }
86 
87 // result += alpha * A * x
88 template<class ResultV, class M, class V>
89 void gemv(
90  matrix_expression<M> const& A,
91  vector_expression<V> const& x,
92  vector_expression<ResultV>& result,
93  typename ResultV::value_type alpha,
94  boost::mpl::false_
95 ) {
96  SIZE_CHECK(A().size1()==result().size());
97  SIZE_CHECK(A().size2()==x().size());
98  typedef typename M::orientation orientation;
99 
100  gemv_impl(A, x, result, alpha, orientation());
101 }
102 
103 }}}
104 #endif