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 SHARK_LINALG_BLAS_KERNELS_VECTOR_ASSIGN_HPP
29 #define SHARK_LINALG_BLAS_KERNELS_VECTOR_ASSIGN_HPP
30 
31 #include "../detail/functional.hpp"
32 #include "../expression_types.hpp"
33 
34 namespace shark {
35 namespace blas {
36 namespace kernels {
37 
38 template<template <class T1, class T2> class F, class V>
39 void assign(vector_expression<V>& v, typename V::value_type t) {
40  typedef F<typename V::iterator::reference, typename V::value_type> Function;
41  Function f;
42  typedef typename V::iterator iterator;
43  iterator end = v().end();
44  for (iterator it = v().begin(); it != end; ++it){
45  f(*it, t);
46  }
47 }
48 
49 /////////////////////////////////////////////////////////
50 //direct assignment of two vectors
51 ////////////////////////////////////////////////////////
52 
53 // Dense-Dense case
54 template< class V, class E>
55 void assign(
57  dense_random_access_iterator_tag, dense_random_access_iterator_tag
58 ) {
59  SIZE_CHECK(v().size() == e().size());
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
65 template< class V, class E>
66 void assign(
68  dense_random_access_iterator_tag, packed_random_access_iterator_tag
69 ) {
70  SIZE_CHECK(v().size() == e().size());
71  typedef typename E::const_iterator EIterator;
72  typedef typename V::scalar_type scalar_type;
73  EIterator eiter = e.begin();
74  EIterator eend = e.end();
75  //special case:
76  //right hand side is completely 0
77  if(eiter == eend){
78  v().clear();
79  return;
80  }
81  EIterator viter = v.begin();
82  EIterator vend = v.end();
83 
84  //set the first elements to zero
85  for(;viter.index() != eiter.index(); ++viter){
86  *viter= scalar_type/*zero*/();
87  }
88  //copy contents of right-hand side
89  for(;eiter != eend; ++eiter,++viter){
90  *viter= *eiter;
91  }
92 
93  for(;viter!= vend; ++viter){
94  *viter= scalar_type/*zero*/();
95  }
96 }
97 
98 // packed-packed case
99 template< class V, class E>
100 void assign(
102  packed_random_access_iterator_tag, packed_random_access_iterator_tag
103 ) {
104  SIZE_CHECK(v().size() == e().size());
105  typedef typename E::const_iterator EIterator;
106  EIterator eiter = e.begin();
107  EIterator eend = e.end();
108  //special case:
109  //right hand side is completely 0
110  if(eiter == eend){
111  v().clear();
112  return;
113  }
114  EIterator viter = v.begin();
115  EIterator vend = v.end();
116 
117  //check for compatible layout
118  SIZE_CHECK(vend-viter);//empty ranges can't be compatible
119  //check whether the right hand side range is included in the left hand side range
120  SIZE_CHECK(viter.index() <= eiter.index());
121  SIZE_CHECK(viter.index()+(vend-viter) >= eiter.index()+(eend-eiter));
122 
123  //copy contents of right-hand side
124  viter += eiter.index()-viter.index();
125  for(;eiter != eend; ++eiter,++viter){
126  *viter= *eiter;
127  }
128 }
129 
130 //Dense-Sparse case
131 template<class V, class E>
132 void assign(
134  vector_expression<E> const& e,
135  dense_random_access_iterator_tag,
136  sparse_bidirectional_iterator_tag
137 ) {
138  v().clear();
139  typedef typename E::const_iterator iterator;
140  iterator end = e().end();
141  for(iterator it = e().begin(); it != end; ++it){
142  v()(it.index()) = *it;
143  }
144 }
145 //Sparse-Dense
146 template<class V, class E>
147 void assign(
149  vector_expression<E> const& e,
150  sparse_bidirectional_iterator_tag,
151  dense_random_access_iterator_tag
152 ) {
153  v().clear();
154  v().reserve(e().size());
155  typename V::iterator targetPos = v().begin();
156  RANGE_CHECK(targetPos == v().end());//as v is cleared, pos must be equal to end
157  for(std::size_t i = 0; i != e().size(); ++i,++targetPos){
158  targetPos = v().set_element(targetPos,i,e()(i));
159  }
160 }
161 // Sparse-Sparse case
162 template<class V, class E>
163 void assign(
165  vector_expression<E> const& e,
166  sparse_bidirectional_iterator_tag,
167  sparse_bidirectional_iterator_tag
168 ) {
169  v().clear();
170  typedef typename E::const_iterator iteratorE;
171  typename V::iterator targetPos = v().begin();
172  RANGE_CHECK(targetPos == v().end());//as v is cleared, pos must be equal to end
173  iteratorE end = e().end();
174  for(iteratorE it = e().begin(); it != end; ++it,++targetPos){
175  targetPos = v().set_element(targetPos,it.index(),*it);
176  }
177 }
178 
179 //dispatcher
180 template< class V, class E>
182  SIZE_CHECK(v().size() == e().size());
183  typedef typename V::const_iterator::iterator_category CategoryV;
184  typedef typename E::const_iterator::iterator_category CategoryE;
185  assign(v, e, CategoryV(),CategoryE());
186 }
187 
188 ////////////////////////////////////////////
189 //assignment with functor
190 ////////////////////////////////////////////
191 
192 //dense dense case
193 template<class V, class E, class F>
194 void assign(
196  vector_expression<E> const& e,
197  F f,
198  dense_random_access_iterator_tag, dense_random_access_iterator_tag
199 ) {
200  SIZE_CHECK(v().size() == e().size());
201  for(std::size_t i = 0; i != v().size(); ++i){
202  f(v()(i),e()(i));
203  }
204 }
205 
206 //dense packed case
207 template<class V, class E, class F>
208 void assign(
210  vector_expression<E> const& e,
211  F f,
212  dense_random_access_iterator_tag, packed_random_access_iterator_tag
213 ) {
214  SIZE_CHECK(v().size() == e().size());
215  typedef typename E::const_iterator EIterator;
216  typedef typename V::const_iterator VIterator;
217  typedef typename V::scalar_type scalar_type;
218  EIterator eiter = e().begin();
219  EIterator eend = e().end();
220  VIterator viter = v().begin();
221  VIterator vend = v().end();
222  //right hand side hasnonzero elements
223  if(eiter != eend){
224  //apply f to the first elements for which the right hand side is 0, unless f is the identity
225  for(;viter.index() != eiter.index() &&!F::right_zero_identity; ++viter){
226  f(*viter,scalar_type/*zero*/());
227  }
228  //copy contents of right-hand side
229  for(;eiter != eend; ++eiter,++viter){
230  f(*viter,*eiter);
231  }
232  }
233  //apply f to the last elements for which the right hand side is 0, unless f is the identity
234  for(;viter!= vend &&!F::right_zero_identity; ++viter){
235  *viter= scalar_type/*zero*/();
236  }
237 }
238 
239 //packed-packed case
240 template<class V, class E, class F>
241 void assign(
243  vector_expression<E> const& e,
244  F f,
245  packed_random_access_iterator_tag, packed_random_access_iterator_tag
246 ) {
247  SIZE_CHECK(v().size() == e().size());
248  typedef typename E::const_iterator EIterator;
249  typedef typename V::const_iterator VIterator;
250  typedef typename V::scalar_type scalar_type;
251  EIterator eiter = e().begin();
252  EIterator eend = e().end();
253  VIterator viter = v().begin();
254  VIterator vend = v().end();
255 
256  //right hand side has nonzero elements
257  if(eiter != eend){
258 
259  //check for compatible layout
260  SIZE_CHECK(vend-viter);//empty ranges can't be compatible
261  //check whether the right hand side range is included in the left hand side range
262  SIZE_CHECK(viter.index() <= eiter.index());
263  SIZE_CHECK(viter.index()+(vend-viter) >= eiter.index()+(eend-eiter));
264 
265  //apply f to the first elements for which the right hand side is 0, unless f is the identity
266  for(;viter.index() != eiter.index() &&!F::right_zero_identity; ++viter){
267  f(*viter,scalar_type/*zero*/());
268  }
269  //copy contents of right-hand side
270  for(;eiter != eend; ++eiter,++viter){
271  f(*viter,*eiter);
272  }
273  }
274  //apply f to the last elements for which the right hand side is 0, unless f is the identity
275  for(;viter!= vend &&!F::right_zero_identity; ++viter){
276  *viter= scalar_type/*zero*/();
277  }
278 }
279 
280 //Dense-Sparse case
281 template<class V, class E, class F>
282 void assign(
284  vector_expression<E> const& e,
285  F f,
286  dense_random_access_iterator_tag, sparse_bidirectional_iterator_tag
287 ) {
288  typedef typename E::const_iterator iterator;
289  iterator end = e().end();
290  for(iterator it = e().begin(); it != end; ++it){
291  f(v()(it.index()),*it);
292  }
293 }
294 
295 //sparse-dense case
296 template<class V, class E, class F>
297 void assign(
299  vector_expression<E> const& e,
300  F f,
301  sparse_bidirectional_iterator_tag tag, dense_random_access_iterator_tag
302 ){
303  typedef typename V::value_type value_type;
304  typedef typename V::size_type size_type;
305  value_type zero = value_type();
306  size_type size = e().size();
307 
308  typename V::iterator it = v().begin();
309  for(size_type i = 0; i != size; ++i,++it){
310  if(it == v().end() || it.index() != i){//insert missing elements
311  it = v().set_element(it,i,zero);
312  }
313  f(*it, e()(i));
314  }
315 }
316 
317 // Sparse-Sparse case has three implementations.
318 //the stupidity of this case is, that we have to assume in the general case v and e share the same
319 //array memory and thus changing v might invalidate the iterators of e.
320 //This is not the same as aliasing of v and e, as v might be for example one matrix row and e another
321 //of the same matrix.
322 //thus we look at the cases where (at least) one of the arguments is a vector-container, which means
323 //that we are not facing the problem of same memory as this would otherwise mean that we are aliasing
324 //in which case the expression is not defined anyways.
325 
326 //called for independent argumeents v and e
327 template<class V, class E, class F>
330  vector_expression<E> const& e,
331  F f
332 ){
333  typedef typename V::value_type value_type;
334  typedef typename V::size_type size_type;
335  value_type zero = value_type();
336  size_type size = v().size();
337 
338  typename V::iterator it = v().begin();
339  typename E::const_iterator ite = e().begin();
340  typename E::const_iterator ite_end = e().end();
341  while(it != v().end() && ite != ite_end) {
342  size_type it_index = it.index();
343  size_type ite_index = ite.index();
344  if (it_index == ite_index) {
345  f(*it, *ite);
346  ++ ite;
347  } else if (it_index < ite_index) {
348  f(*it, zero);
349  } else{//it_index > ite_index so insert new element in v()
350  it = v().set_element(it,ite_index,zero);
351  f(*it, *ite);
352  ++ite;
353  }
354  ++it;
355  }
356  //apply zero transformation on elements which are not transformed yet
357  for(;it != v().end();++it) {
358  f(*it, zero);
359  }
360  //add missing elements
361  for(;ite != ite_end;++it,++ite) {
362  it = v().set_element(it,ite.index(),zero);
363  f(*it, *ite);
364  }
365 }
366 //as long as one argument is not a proxy, we are in the good case.
367 template<class V, class E, class F>
368 void assign(
370  vector_container<E> const& e,
371  F f,
372  sparse_bidirectional_iterator_tag tag, sparse_bidirectional_iterator_tag
373 ){
374  assign_sparse(v,e);
375 }
376 template<class V, class E, class F>
377 void assign(
379  vector_expression<E> const& e,
380  F f,
381  sparse_bidirectional_iterator_tag tag, sparse_bidirectional_iterator_tag
382 ){
383  assign_sparse(v,e,f);
384 }
385 template<class V, class E, class F>
386 void assign(
388  vector_container<E> const& e,
389  F f,
390  sparse_bidirectional_iterator_tag tag, sparse_bidirectional_iterator_tag
391 ){
392  assign_sparse(v,e,f);
393 }
394 
395 //In the general case we have to take one bullet:
396 //either use a temporary, which has copying time and allocation time
397 //or count the non-zero elements of e which might be expensive as well. we decide
398 //to take the first route for now.
399 template<class V, class E, class F>
400 void assign(
402  vector_expression<E> const& e,
403  F f,
404  sparse_bidirectional_iterator_tag tag, sparse_bidirectional_iterator_tag
405 ){
406  typename vector_temporary<V>::type temporary(v());
407  assign_sparse(temporary,e, f);
408  v().clear();
409  kernels::assign(v, temporary);
410 }
411 
412 // Dispatcher
413 template<template <class T1, class T2> class F, class V, class E>
415  SIZE_CHECK(v().size() == e().size());
416  typedef typename V::const_iterator::iterator_category CategoryV;
417  typedef typename E::const_iterator::iterator_category CategoryE;
418  typedef F<typename V::iterator::reference, typename E::value_type> functor_type;
419  assign(v(), e(), functor_type(), CategoryV(),CategoryE());
420 }
421 
422 }}}
423 #endif