dot.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 REMORA_KERNELS_DEFAULT_DOT_HPP
31 #define REMORA_KERNELS_DEFAULT_DOT_HPP
32 
33 #include "../../expression_types.hpp"//vector_expression
34 #include "../../detail/traits.hpp"//storage tags
35 
36 namespace remora{namespace bindings{
37 
38 // Dense case
39 template<class E1, class E2, class result_type>
40 void dot(
41  vector_expression<E1, cpu_tag> const& v1,
42  vector_expression<E2, cpu_tag> const& v2,
43  result_type& result,
44  dense_tag,
45  dense_tag
46 ) {
47  std::size_t size = v1().size();
48  result = result_type();
49  for(std::size_t i = 0; i != size; ++i){
50  result += v1()(i) * v2()(i);
51  }
52 }
53 // Sparse case
54 template<class E1, class E2, class result_type>
55 void dot(
56  vector_expression<E1, cpu_tag> const& v1,
57  vector_expression<E2, cpu_tag> const& v2,
58  result_type& result,
59  sparse_tag,
60  sparse_tag
61 ) {
62  typename E1::const_iterator iter1=v1().begin();
63  typename E1::const_iterator end1=v1().end();
64  typename E2::const_iterator iter2=v2().begin();
65  typename E2::const_iterator end2=v2().end();
66  result = result_type();
67  //be aware of empty vectors!
68  while(iter1 != end1 && iter2 != end2)
69  {
70  std::size_t index1=iter1.index();
71  std::size_t index2=iter2.index();
72  if(index1==index2){
73  result += *iter1 * *iter2;
74  ++iter1;
75  ++iter2;
76  }
77  else if(index1> index2){
78  ++iter2;
79  }
80  else {
81  ++iter1;
82  }
83  }
84 }
85 
86 // Dense-Sparse case
87 template<class E1, class E2, class result_type>
88 void dot(
89  vector_expression<E1, cpu_tag> const& v1,
90  vector_expression<E2, cpu_tag> const& v2,
91  result_type& result,
92  dense_tag,
93  sparse_tag
94 ) {
95  typename E2::const_iterator iter2=v2().begin();
96  typename E2::const_iterator end2=v2().end();
97  result = result_type();
98  for(;iter2 != end2;++iter2){
99  result += v1()(iter2.index()) * *iter2;
100  }
101 }
102 //Sparse-Dense case is reduced to Dense-Sparse using symmetry.
103 template<class E1, class E2, class result_type>
104 void dot(
105  vector_expression<E1, cpu_tag> const& v1,
106  vector_expression<E2, cpu_tag> const& v2,
107  result_type& result,
108  sparse_tag t1,
109  dense_tag t2
110 ) {
111  //use commutativity!
112  dot(v2,v1,result,t2,t1);
113 }
114 
115 }}
116 #endif