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
36namespace remora{namespace bindings{
37
38// Dense case
39template<class E1, class E2, class result_type>
40void 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
54template<class E1, class E2, class result_type>
55void 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
87template<class E1, class E2, class result_type>
88void 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.
103template<class E1, class E2, class result_type>
104void 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