vector_assign.hpp
Go to the documentation of this file.
1/*!
2 * \brief Assignment kernels for vector expressions
3 *
4 * \author O. Krause
5 * \date 2015
6 *
7 *
8 * \par Copyright 1995-2015 Shark Development Team
9 *
10 * <BR><HR>
11 * This file is part of Shark.
12 * <http://image.diku.dk/shark/>
13 *
14 * Shark is free software: you can redistribute it and/or modify
15 * it under the terms of the GNU Lesser General Public License as published
16 * by the Free Software Foundation, either version 3 of the License, or
17 * (at your option) any later version.
18 *
19 * Shark is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU Lesser General Public License for more details.
23 *
24 * You should have received a copy of the GNU Lesser General Public License
25 * along with Shark. If not, see <http://www.gnu.org/licenses/>.
26 *
27 */
28#ifndef REMORA_KERNELS_DEFAULT_VECTOR_ASSIGN_HPP
29#define REMORA_KERNELS_DEFAULT_VECTOR_ASSIGN_HPP
30
31#include "../../detail/traits.hpp"
32
33namespace remora{namespace bindings{
34
35template<class F, class V>
36void apply(vector_expression<V, cpu_tag>& v, F const& f) {
37 typedef typename V::iterator iterator;
38 iterator end = v().end();
39 for (iterator it = v().begin(); it != end; ++it){
40 *it = f(*it);
41 }
42}
43
44template<class F, class V>
45void assign(vector_expression<V, cpu_tag>& v, typename V::value_type t) {
46 F f;
47 apply(v, [=](typename V::value_type x){return f(x,t);});
48}
49
50/////////////////////////////////////////////////////////
51//direct assignment of two vectors
52////////////////////////////////////////////////////////
53
54// Dense-Dense case
55template< class V, class E>
56void vector_assign(
57 vector_expression<V, cpu_tag>& v, vector_expression<E, cpu_tag> const& e,
58 dense_tag, dense_tag
59) {
60 for(std::size_t i = 0; i != v().size(); ++i){
61 v()(i) = static_cast<typename V::value_type>(e()(i));
62 }
63}
64// Dense-packed case
65template< class V, class E>
66void vector_assign(
67 vector_expression<V, cpu_tag>& v, vector_expression<E, cpu_tag> const& e,
68 dense_tag, packed_tag
69) {
70 typedef typename E::const_iterator EIterator;
71 typedef typename V::value_type value_type;
72 EIterator eiter = e.begin();
73 EIterator eend = e.end();
74 //special case:
75 //right hand side is completely 0
76 if(eiter == eend){
77 v().clear();
78 return;
79 }
80 EIterator viter = v.begin();
81 EIterator vend = v.end();
82
83 //set the first elements to zero
84 for(;viter.index() != eiter.index(); ++viter){
85 *viter= value_type/*zero*/();
86 }
87 //copy contents of right-hand side
88 for(;eiter != eend; ++eiter,++viter){
89 *viter= *eiter;
90 }
91
92 for(;viter!= vend; ++viter){
93 *viter= value_type/*zero*/();
94 }
95}
96
97// packed-packed case
98template< class V, class E>
99void vector_assign(
100 vector_expression<V, cpu_tag>& v, vector_expression<E, cpu_tag> const& e,
101 packed_tag, packed_tag
102) {
103 typedef typename E::const_iterator EIterator;
104 EIterator eiter = e.begin();
105 EIterator eend = e.end();
106 //special case:
107 //right hand side is completely 0
108 if(eiter == eend){
109 v().clear();
110 return;
111 }
112 EIterator viter = v.begin();
113 EIterator vend = v.end();
114
115 //check for compatible layout
116 REMORA_SIZE_CHECK(vend-viter);//empty ranges can't be compatible
117 //check whether the right hand side range is included in the left hand side range
118 REMORA_SIZE_CHECK(viter.index() <= eiter.index());
119 REMORA_SIZE_CHECK(viter.index()+(vend-viter) >= eiter.index()+(eend-eiter));
120
121 //copy contents of right-hand side
122 viter += eiter.index()-viter.index();
123 for(;eiter != eend; ++eiter,++viter){
124 *viter= *eiter;
125 }
126}
127
128//Dense-Sparse case
129template<class V, class E>
130void vector_assign(
131 vector_expression<V, cpu_tag>& v,
132 vector_expression<E, cpu_tag> const& e,
133 dense_tag,
134 sparse_tag
135) {
136 v().clear();
137 typedef typename E::const_iterator iterator;
138 iterator end = e().end();
139 for(iterator it = e().begin(); it != end; ++it){
140 v()(it.index()) = *it;
141 }
142}
143//Sparse-Dense
144template<class V, class E>
145void vector_assign(
146 vector_expression<V, cpu_tag>& v,
147 vector_expression<E, cpu_tag> const& e,
148 sparse_tag,
149 dense_tag
150) {
151 v().clear();
152 v().reserve(e().size());
153 typename V::iterator targetPos = v().begin();
154 REMORA_RANGE_CHECK(targetPos == v().end());//as v is cleared, pos must be equal to end
155 for(std::size_t i = 0; i != e().size(); ++i){
156 targetPos = v().set_element(targetPos,i,e()(i));
157 }
158}
159// Sparse-Sparse case
160template<class V, class E>
161void vector_assign(
162 vector_expression<V, cpu_tag>& v,
163 vector_expression<E, cpu_tag> const& e,
164 sparse_tag,
165 sparse_tag
166) {
167 v().clear();
168 typedef typename E::const_iterator iteratorE;
169 typename V::iterator targetPos = v().begin();
170 REMORA_RANGE_CHECK(targetPos == v().end());//as v is cleared, pos must be equal to end
171 iteratorE end = e().end();
172 for(iteratorE it = e().begin(); it != end; ++it){
173 targetPos = v().set_element(targetPos,it.index(),*it);
174 }
175}
176
177////////////////////////////////////////////
178//assignment with functor
179////////////////////////////////////////////
180
181//dense dense case
182template<class V, class E, class F>
183void vector_assign_functor(
184 vector_expression<V, cpu_tag>& v,
185 vector_expression<E, cpu_tag> const& e,
186 F f,
187 dense_tag, dense_tag
188) {
189 for(std::size_t i = 0; i != v().size(); ++i){
190 v()(i) = f(v()(i),e()(i));
191 }
192}
193
194//dense packed case
195template<class V, class E, class F>
196void vector_assign_functor(
197 vector_expression<V, cpu_tag>& v,
198 vector_expression<E, cpu_tag> const& e,
199 F f,
200 dense_tag, packed_tag
201) {
202 typedef typename E::const_iterator EIterator;
203 typedef typename V::const_iterator VIterator;
204 typedef typename V::value_type value_type;
205 EIterator eiter = e().begin();
206 EIterator eend = e().end();
207 VIterator viter = v().begin();
208 VIterator vend = v().end();
209 //right hand side hasnonzero elements
210 if(eiter != eend){
211 //apply f to the first elements for which the right hand side is 0, unless f is the identity
212 for(;viter.index() != eiter.index() &&!F::right_zero_identity; ++viter){
213 *viter = f(*viter,value_type/*zero*/());
214 }
215 //copy contents of right-hand side
216 for(;eiter != eend; ++eiter,++viter){
217 *viter = f(*viter,*eiter);
218 }
219 }
220 //apply f to the last elements for which the right hand side is 0, unless f is the identity
221 for(;viter!= vend &&!F::right_zero_identity; ++viter){
222 *viter= value_type/*zero*/();
223 }
224}
225
226//packed-packed case
227template<class V, class E, class F>
228void vector_assign_functor(
229 vector_expression<V, cpu_tag>& v,
230 vector_expression<E, cpu_tag> const& e,
231 F f,
232 packed_tag, packed_tag
233) {
234 typedef typename E::const_iterator EIterator;
235 typedef typename V::const_iterator VIterator;
236 typedef typename V::value_type value_type;
237 EIterator eiter = e().begin();
238 EIterator eend = e().end();
239 VIterator viter = v().begin();
240 VIterator vend = v().end();
241
242 //right hand side has nonzero elements
243 if(eiter != eend){
244
245 //check for compatible layout
246 REMORA_SIZE_CHECK(vend-viter);//empty ranges can't be compatible
247 //check whether the right hand side range is included in the left hand side range
248 REMORA_SIZE_CHECK(viter.index() <= eiter.index());
249 REMORA_SIZE_CHECK(viter.index()+(vend-viter) >= eiter.index()+(eend-eiter));
250
251 //apply f to the first elements for which the right hand side is 0, unless f is the identity
252 for(;viter.index() != eiter.index() &&!F::right_zero_identity; ++viter){
253 *viter = f(*viter,value_type/*zero*/());
254 }
255 //copy contents of right-hand side
256 for(;eiter != eend; ++eiter,++viter){
257 *viter = f(*viter,*eiter);
258 }
259 }
260 //apply f to the last elements for which the right hand side is 0, unless f is the identity
261 for(;viter!= vend &&!F::right_zero_identity; ++viter){
262 *viter= f(*viter,value_type/*zero*/());
263 }
264}
265
266//Dense-Sparse case
267template<class V, class E, class F>
268void vector_assign_functor(
269 vector_expression<V, cpu_tag>& v,
270 vector_expression<E, cpu_tag> const& e,
271 F f,
272 dense_tag, sparse_tag
273) {
274 typedef typename E::const_iterator iterator;
275 iterator end = e().end();
276 for(iterator it = e().begin(); it != end; ++it){
277 v()(it.index()) = f(v()(it.index()),*it);
278 }
279}
280
281//sparse-dense case
282template<class V, class E, class F>
283void vector_assign_functor(
284 vector_expression<V, cpu_tag>& v,
285 vector_expression<E, cpu_tag> const& e,
286 F f,
287 sparse_tag tag, dense_tag
288){
289 typedef typename V::size_type size_type;
290 size_type size = e().size();
291 auto it = v().begin();
292 for(size_type i = 0; i != size; ++i){
293 auto val = e()(i);
294 if(it == v().end() || it.index() != i){//insert missing elements
295 it = v().set_element(it,i,val);
296 }else{
297 *it = f(*it, val);
298 ++it;
299 }
300 }
301}
302
303template<class V, class E, class F>
304void vector_assign_functor(
305 vector_expression<V, cpu_tag>& v,
306 vector_expression<E, cpu_tag> const& e,
307 F f,
308 sparse_tag tag, sparse_tag
309){
310 typedef typename V::value_type value_type;
311 typedef typename V::size_type size_type;
312 value_type zero = value_type();
313
314 typename V::iterator it = v().begin();
315 typename E::const_iterator ite = e().begin();
316 typename E::const_iterator ite_end = e().end();
317 while(it != v().end() && ite != ite_end) {
318 size_type it_index = it.index();
319 size_type ite_index = ite.index();
320 if (it_index == ite_index) {
321 *it = f(*it, *ite);
322 ++ ite;
323 ++it;
324 } else if (it_index < ite_index) {
325 *it = f(*it, zero);
326 ++it;
327 } else{//it_index > ite_index so insert new element in v()
328 it = v().set_element(it,ite_index,f(zero, *ite));
329 ++ite;
330 }
331
332 }
333 //apply zero transformation on elements which are not transformed yet
334 for(;it != v().end();++it) {
335 *it = f(*it, zero);
336 }
337 //add missing elements
338 for(;ite != ite_end;++it,++ite) {
339 it = v().set_element(it,ite.index(),zero);
340 *it = f(*it, *ite);
341 }
342}
343
344}}
345#endif