viennacl/vector.hpp

Go to the documentation of this file.

#ifndef VIENNACL_VECTOR_HPP_
#define VIENNACL_VECTOR_HPP_

/* =========================================================================
   Copyright (c) 2010-2013, Institute for Microelectronics,
                            Institute for Analysis and Scientific Computing,
                            TU Wien.
   Portions of this software are copyright by UChicago Argonne, LLC.

                            -----------------
                  ViennaCL - The Vienna Computing Library
                            -----------------

   Project Head:    Karl Rupp                   rupp@iue.tuwien.ac.at
               
   (A list of authors and contributors can be found in the PDF manual)

   License:         MIT (X11), see file LICENSE in the base directory
============================================================================= */

#include "viennacl/forwards.h"
#include "viennacl/backend/memory.hpp"
#include "viennacl/scalar.hpp"
#include "viennacl/tools/tools.hpp"
#include "viennacl/tools/entry_proxy.hpp"
#include "viennacl/linalg/vector_operations.hpp"
#include "viennacl/meta/result_of.hpp"


namespace viennacl
{
  //
  // Initializer types
  //
  template <typename SCALARTYPE>
  class unit_vector
  {
    public:
      typedef vcl_size_t        size_type;
      
      unit_vector(size_type s, size_type ind) : size_(s), index_(ind) 
      {
        assert( (ind < s) && bool("Provided index out of range!") );
      }
      
      size_type size() const { return size_; }
      size_type index() const { return index_; }
      
    private:
      size_type size_;
      size_type index_;
  };

  
  template <typename SCALARTYPE>
  class zero_vector
  {
    public:
      typedef vcl_size_t        size_type;
      typedef SCALARTYPE        const_reference;
      
      zero_vector(size_type s) : size_(s) {}
      
      size_type size() const { return size_; }
      const_reference operator()(size_type /*i*/) const { return 0; }
      const_reference operator[](size_type /*i*/) const { return 0; }
      
    private:
      size_type size_;
  };
  
  
  template <typename SCALARTYPE>
  class scalar_vector
  {
    public:
      typedef vcl_size_t         size_type;
      typedef SCALARTYPE const & const_reference;
      
      scalar_vector(size_type s, SCALARTYPE val) : size_(s), value_(val) {}
      
      size_type size() const { return size_; }
      const_reference operator()(size_type /*i*/) const { return value_; }
      const_reference operator[](size_type /*i*/) const { return value_; }
      
    private:
      size_type size_;
      SCALARTYPE value_;
  };
  
  
  //
  // Vector expression
  //
    
  template <typename LHS, typename RHS, typename OP>
  class vector_expression
  {
      typedef typename result_of::reference_if_nonscalar<LHS>::type     lhs_reference_type;
      typedef typename result_of::reference_if_nonscalar<RHS>::type     rhs_reference_type;
    
    public:
      enum { alignment = 1 };
      
      typedef typename viennacl::tools::VECTOR_EXTRACTOR<LHS, RHS>::ResultType    VectorType;
      typedef vcl_size_t       size_type;
      
      vector_expression(LHS & l, RHS & r) : lhs_(l), rhs_(r) {}
      
      lhs_reference_type lhs() const { return lhs_; }
      rhs_reference_type rhs() const { return rhs_; }
      
      size_type size() const { return viennacl::traits::size(*this); }
      
    private:
      lhs_reference_type lhs_;
      rhs_reference_type rhs_;
  };
  
  template<class SCALARTYPE, unsigned int ALIGNMENT>
  class const_vector_iterator
  {
      typedef const_vector_iterator<SCALARTYPE, ALIGNMENT>    self_type;
    public:
      typedef scalar<SCALARTYPE>            value_type;
      typedef long                          difference_type;
      typedef backend::mem_handle           handle_type;
      
      //const_vector_iterator() {};
      
      const_vector_iterator(vector_base<SCALARTYPE> const & vec,
                            std::size_t index,
                            std::size_t start = 0,
                            vcl_ptrdiff_t stride = 1) : elements_(vec.handle()), index_(index), start_(start), stride_(stride) {};
                            
      const_vector_iterator(handle_type const & elements,
                            std::size_t index,
                            std::size_t start = 0,
                            vcl_ptrdiff_t stride = 1) : elements_(elements), index_(index), start_(start), stride_(stride) {};

      value_type operator*(void) const 
      { 
          value_type result;
          result = const_entry_proxy<SCALARTYPE>(start_ + index_ * stride_, elements_);
          return result;
      }
      self_type operator++(void) { index_ += stride_; return *this; }
      self_type operator++(int) { self_type tmp = *this; ++(*this); return tmp; }
      
      bool operator==(self_type const & other) const { return index_ == other.index_; }
      bool operator!=(self_type const & other) const { return index_ != other.index_; }
      
//        self_type & operator=(self_type const & other)
//        {
//           index_ = other._index;
//           elements_ = other._elements;
//           return *this;
//        }   

      difference_type operator-(self_type const & other) const 
      {
        assert( (other.start_ == start_) && (other.stride_ == stride_) && bool("Iterators are not from the same vector (proxy)!"));
        return static_cast<difference_type>(index_) - static_cast<difference_type>(other.index_); 
      }
      self_type operator+(difference_type diff) const { return self_type(elements_, index_ + diff * stride_, start_, stride_); }
      
      //std::size_t index() const { return index_; }
      std::size_t offset() const { return start_ + index_ * stride_; }
      
      std::size_t stride() const { return stride_; }
      handle_type const & handle() const { return elements_; }

    protected:
      handle_type const & elements_;
      std::size_t index_;  //offset from the beginning of elements_
      std::size_t start_;
      vcl_ptrdiff_t stride_;
  };
  

  template<class SCALARTYPE, unsigned int ALIGNMENT>
  class vector_iterator : public const_vector_iterator<SCALARTYPE, ALIGNMENT>
  {
      typedef const_vector_iterator<SCALARTYPE, ALIGNMENT>  base_type;
      typedef vector_iterator<SCALARTYPE, ALIGNMENT>        self_type;
    public:
      typedef typename base_type::handle_type               handle_type;
      typedef typename base_type::difference_type           difference_type;
      
      vector_iterator() : base_type(), elements_(NULL) {};
      vector_iterator(handle_type & elements,
                      std::size_t index,
                      std::size_t start = 0,
                      vcl_ptrdiff_t stride = 1)  : base_type(elements, index, start, stride), elements_(elements) {};
      vector_iterator(vector_base<SCALARTYPE> & vec,
                      std::size_t index,
                      std::size_t start = 0,
                      vcl_ptrdiff_t stride = 1) : base_type(vec, index, start, stride), elements_(vec.handle()) {};
      //vector_iterator(base_type const & b) : base_type(b) {};

      typename base_type::value_type operator*(void)  
      { 
          typename base_type::value_type result;
          result = entry_proxy<SCALARTYPE>(base_type::start_ + base_type::index_ * base_type::stride_, elements_); 
          return result;
      }
      
      difference_type operator-(self_type const & other) const { difference_type result = base_type::index_; return (result - static_cast<difference_type>(other.index_)); }
      self_type operator+(difference_type diff) const { return self_type(elements_, base_type::index_ + diff * base_type::stride_, base_type::start_, base_type::stride_); }
      
      handle_type       & handle()       { return elements_; }
      handle_type const & handle() const { return base_type::elements_; }
      
      //operator base_type() const
      //{
      //  return base_type(base_type::elements_, base_type::index_, base_type::start_, base_type::stride_);
      //}
    private:
      handle_type & elements_;
  };
  
  
  template<class SCALARTYPE, typename SizeType /* see forwards.h for default type */, typename DistanceType /* see forwards.h for default type */>
  class vector_base
  {
      typedef vector_base<SCALARTYPE>         self_type;
      
    public:
      typedef scalar<typename viennacl::tools::CHECK_SCALAR_TEMPLATE_ARGUMENT<SCALARTYPE>::ResultType>   value_type;
      typedef SCALARTYPE                                        cpu_value_type;
      typedef backend::mem_handle                               handle_type;
      typedef SizeType                                        size_type;
      typedef DistanceType                                     difference_type;
      typedef const_vector_iterator<SCALARTYPE, 1>              const_iterator;
      typedef vector_iterator<SCALARTYPE, 1>                    iterator;
      
      static const size_type alignment = 1;
  
      explicit vector_base() : size_(0), start_(0), stride_(1) { /* Note: One must not call ::init() here because a vector might have been created globally before the backend has become available */ }
  
      explicit vector_base(viennacl::backend::mem_handle & h,
                           size_type vec_size, size_type vec_start, difference_type vec_stride) : size_(vec_size), start_(vec_start), stride_(vec_stride), elements_(h) {}
      
      explicit vector_base(size_type vec_size) : size_(vec_size), start_(0), stride_(1)
      {
        if (size_ > 0)
        {
          std::vector<SCALARTYPE> temp(internal_size());
          viennacl::backend::memory_create(elements_, sizeof(SCALARTYPE)*internal_size(), &(temp[0]));
          pad();
        }
      }
  
      explicit vector_base(size_type vec_size, viennacl::memory_types mem_type) : size_(vec_size), start_(0), stride_(1)
      {
        elements_.switch_active_handle_id(mem_type);
        if (size_ > 0)
        {
          std::vector<SCALARTYPE> temp(internal_size());
          viennacl::backend::memory_create(elements_, sizeof(SCALARTYPE)*internal_size(), &(temp[0]));
          pad();
        }
      }
      
      //
      // operator=
      //
      
  
      self_type & operator=(const self_type & vec)
      {
        assert( ( (vec.size() == size()) || (size() == 0) )
                && bool("Incompatible vector sizes!"));

        if (vec.size() > 0)
        {
          if (size_ == 0)
          {
            size_ = vec.size();
            elements_.switch_active_handle_id(vec.handle().get_active_handle_id());
            viennacl::backend::memory_create(elements_, sizeof(SCALARTYPE)*internal_size());
          }
            
          viennacl::linalg::av(*this,
                               vec, cpu_value_type(1.0), 1, false, false);
        }
        
        return *this;
      }
  
  
      template <typename T, typename S1, typename OP>
      typename viennacl::enable_if< viennacl::is_any_scalar<S1>::value,
                                    self_type & >::type
      operator = (const vector_expression< const vector_base<T>, const S1, OP> & proxy)
      {
        assert( ( (proxy.lhs().size() == size()) || (size() == 0) )
                && bool("Incompatible vector sizes!"));
        
        // initialize the necessary buffer
        if (size() == 0)
        {
          size_ = proxy.lhs().size();
          viennacl::backend::memory_create(elements_, sizeof(SCALARTYPE)*internal_size());
          pad();
        } 
  
        viennacl::linalg::av(*this,
                             proxy.lhs(), proxy.rhs(), 1, (viennacl::is_division<OP>::value ? true : false), false);
        
        return *this;
      }
  
      //v1 = v2 +- v3; 
      template <typename T, typename OP>
      typename viennacl::enable_if< viennacl::is_addition<OP>::value || viennacl::is_subtraction<OP>::value,
                                    self_type &>::type
      operator = (const vector_expression< const vector_base<T>,
                                           const vector_base<T>,
                                            OP> & proxy)
      {
        assert( ( (proxy.lhs().size() == size()) || (size() == 0) )
                && bool("Incompatible vector sizes!"));
        
        if (size() == 0)
        {
          size_ = proxy.lhs().size();
          viennacl::backend::memory_create(elements_, sizeof(SCALARTYPE)*internal_size());
          pad();
        } 
  
        viennacl::linalg::avbv(*this, 
                               proxy.lhs(), SCALARTYPE(1.0), 1, false, false,
                               proxy.rhs(), SCALARTYPE(1.0), 1, false, (viennacl::is_subtraction<OP>::value ? true : false));
        return *this;
      }
      
      template <typename T, typename S2, typename OP2,
                typename OP>
      typename viennacl::enable_if< viennacl::is_any_scalar<S2>::value && (viennacl::is_product<OP2>::value || viennacl::is_division<OP2>::value)
                                    && (viennacl::is_addition<OP>::value || viennacl::is_subtraction<OP>::value),
                                    self_type &>::type
      operator = (const vector_expression< const vector_base<T>,
                                           const vector_expression<const vector_base<T>, const S2, OP2>,
                                           OP> & proxy)
      {
        assert( ( (proxy.lhs().size() == size()) || (size() == 0) )
                && bool("Incompatible vector sizes!"));
        
        if (size() == 0)
        {
          size_ = proxy.lhs().size();
          viennacl::backend::memory_create(elements_, sizeof(SCALARTYPE)*internal_size());
          pad();
        } 
  
        bool flip_sign_2 = (viennacl::is_subtraction<OP>::value ? true : false);
        if (viennacl::is_flip_sign_scalar<S2>::value)
          flip_sign_2 = !flip_sign_2;
        viennacl::linalg::avbv(*this, 
                              proxy.lhs(),         SCALARTYPE(1.0), 1, false                                             , false,
                              proxy.rhs().lhs(), proxy.rhs().rhs(), 1, (viennacl::is_division<OP2>::value ? true : false), flip_sign_2);

        return *this;
      }
  
      template <typename T, typename S1, typename OP1,
                typename OP>
      typename viennacl::enable_if< viennacl::is_any_scalar<S1>::value && (viennacl::is_product<OP1>::value || viennacl::is_division<OP1>::value)
                                    && (viennacl::is_addition<OP>::value || viennacl::is_subtraction<OP>::value),
                                    self_type &>::type
      operator = (const vector_expression< const vector_expression<const vector_base<T>, const S1, OP1>,
                                           const vector_base<T>,
                                           OP> & proxy)
      {
        assert( ( (proxy.size() == size()) || (size() == 0) )
                && bool("Incompatible vector sizes!"));
        
        if (size() == 0)
        {
          size_ = proxy.lhs().size();
          viennacl::backend::memory_create(elements_, sizeof(SCALARTYPE)*internal_size());
          pad();
        } 
  
        viennacl::linalg::avbv(*this, 
                              proxy.lhs().lhs(), proxy.lhs().rhs(), 1, (viennacl::is_division<OP1>::value ? true : false), (viennacl::is_flip_sign_scalar<S1>::value ? true : false),
                              proxy.rhs(),         SCALARTYPE(1.0), 1, false                                             , (viennacl::is_subtraction<OP>::value ? true : false));
        return *this;
      }
      
      template <typename T,
                typename S1, typename OP1,
                typename S2, typename OP2,
                typename OP>
      typename viennacl::enable_if<    viennacl::is_any_scalar<S1>::value && (viennacl::is_product<OP1>::value || viennacl::is_division<OP1>::value)
                                    && viennacl::is_any_scalar<S2>::value && (viennacl::is_product<OP2>::value || viennacl::is_division<OP2>::value)
                                    && (viennacl::is_addition<OP>::value || viennacl::is_subtraction<OP>::value),
                                    self_type &>::type
      operator = (const vector_expression< const vector_expression<const vector_base<T>, const S1, OP1>,
                                           const vector_expression<const vector_base<T>, const S2, OP2>,
                                           OP> & proxy)
      {
        assert( ( (proxy.lhs().size() == size()) || (size() == 0) )
                && bool("Incompatible vector sizes!"));
        
        if (size() == 0)
        {
          size_ = proxy.lhs().size();
          viennacl::backend::memory_create(elements_, sizeof(SCALARTYPE)*internal_size());
          pad();
        } 
  
        bool flip_sign_2 = (viennacl::is_subtraction<OP>::value ? true : false);
        if (viennacl::is_flip_sign_scalar<S2>::value)
          flip_sign_2 = !flip_sign_2;
        viennacl::linalg::avbv(*this, 
                              proxy.lhs().lhs(), proxy.lhs().rhs(), 1, (viennacl::is_division<OP1>::value ? true : false), (viennacl::is_flip_sign_scalar<S1>::value ? true : false),
                              proxy.rhs().lhs(), proxy.rhs().rhs(), 1, (viennacl::is_division<OP2>::value ? true : false), flip_sign_2);

        return *this;
      }
      
      
      //v1 = v2 * / v3; 
      template <typename T, typename OP>
      typename viennacl::enable_if< (viennacl::is_product<OP>::value || viennacl::is_division<OP>::value),
                                    self_type &>::type
      operator = (const vector_expression< const vector_base<T>,
                                           const vector_base<T>,
                                           OP> & proxy)
      {
        assert( ( (proxy.lhs().size() == size()) || (size() == 0) )
                && bool("Incompatible vector sizes!"));
        
        if (size() == 0)
        {
          size_ = proxy.lhs().size();
          viennacl::backend::memory_create(elements_, sizeof(SCALARTYPE)*internal_size());
          pad();
        } 
  
        viennacl::linalg::element_op(*this, proxy);
        return *this;
      }
      
      
      
      // assign vector range or vector slice
      template <typename T>
      self_type &
      operator = (const vector_base<T> & v1)
      {
        assert( ( (v1.size() == size()) || (size() == 0) )
                && bool("Incompatible vector sizes!"));
        
        if (size() == 0)
        {
          size_ = v1.size();
          if (size_ > 0)
          {
            viennacl::backend::memory_create(elements_, sizeof(SCALARTYPE)*internal_size());
            pad();
          }
        } 
        
        viennacl::linalg::av(*this, 
                             v1, SCALARTYPE(1.0), 1, false, false);
        
        return *this;
      }
      
      self_type & operator = (unit_vector<SCALARTYPE> const & v)
      {
        assert( ( (v.size() == size()) || (size() == 0) )
                && bool("Incompatible vector sizes!"));
        
        if (size() == 0)
        {
          size_ = v.size();
          if (size_ > 0)
            viennacl::backend::memory_create(elements_, sizeof(SCALARTYPE)*internal_size());
        }
        
        if (size_ > 0)
        {
          clear();
          this->operator()(v.index()) = SCALARTYPE(1);
        }
        
        return *this;
      }
      
      self_type & operator = (zero_vector<SCALARTYPE> const & v)
      {
        assert( ( (v.size() == size()) || (size() == 0) )
                && bool("Incompatible vector sizes!"));
        
        if (size() == 0)
        {
          size_ = v.size();
          if (size_ > 0)
            viennacl::backend::memory_create(elements_, sizeof(SCALARTYPE)*internal_size());
        }
        
        if (size_ > 0)
          clear();
        
        return *this;
      }
  
      self_type & operator = (scalar_vector<SCALARTYPE> const & v)
      {
        assert( ( (v.size() == size()) || (size() == 0) )
                && bool("Incompatible vector sizes!"));
        
        if (size() == 0)
        {
          size_ = v.size();
          if (size_ > 0)
            viennacl::backend::memory_create(elements_, sizeof(SCALARTYPE)*internal_size());
        }
        
        if (size_ > 0)
          viennacl::linalg::vector_assign(*this, v[0]);
        
        return *this;
      }
  
      
      
  
      //Note: The following operator overloads are defined in matrix_operations.hpp, compressed_matrix_operations.hpp and coordinate_matrix_operations.hpp
      //This is certainly not the nicest approach and will most likely by changed in the future, but it works :-)
      
      //matrix<>
      template <typename F>
      self_type & operator=(const viennacl::vector_expression< const matrix_base<SCALARTYPE, F>, const vector_base<SCALARTYPE>, viennacl::op_prod> & proxy)
      {
        assert(viennacl::traits::size1(proxy.lhs()) == size() && bool("Size check failed for v1 = A * v2: size1(A) != size(v1)"));
        
        // check for the special case x = A * x
        if (viennacl::traits::handle(proxy.rhs()) == viennacl::traits::handle(*this))
        {
          viennacl::vector<SCALARTYPE> result(viennacl::traits::size1(proxy.lhs()));
          viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), result);
          *this = result;
        }
        else
        {
          viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), *this);
        }
        return *this;
      }
  
      
      //transposed_matrix_proxy:
      template <typename F>
      self_type & operator=(const vector_expression< const matrix_expression< const matrix_base<SCALARTYPE, F>, const matrix_base<SCALARTYPE, F>, op_trans >,
                                                     const vector_base<SCALARTYPE>,
                                                     op_prod> & proxy)
      {
        assert(viennacl::traits::size1(proxy.lhs()) == size() && bool("Size check failed in v1 = trans(A) * v2: size2(A) != size(v1)"));
  
        // check for the special case x = trans(A) * x
        if (viennacl::traits::handle(proxy.rhs()) == viennacl::traits::handle(*this))
        {
          viennacl::vector<SCALARTYPE> result(viennacl::traits::size1(proxy.lhs()));
          viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), result);
          *this = result;
        }
        else
        {
          viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), *this);
        }
        return *this;
      }
  

  
      //read-write access to an element of the vector
      entry_proxy<SCALARTYPE> operator()(size_type index)
      {
        assert( (size() > 0)  && bool("Cannot apply operator() to vector of size zero!"));
        assert( index < size() && bool("Index out of bounds!") );
        
        return entry_proxy<SCALARTYPE>(start_ + stride_ * index, elements_);
      }
  
      entry_proxy<SCALARTYPE> operator[](size_type index)
      {
        assert( (size() > 0)  && bool("Cannot apply operator() to vector of size zero!"));
        assert( index < size() && bool("Index out of bounds!") );
        
        return entry_proxy<SCALARTYPE>(start_ + stride_ * index, elements_);
      }
  
  
      const_entry_proxy<SCALARTYPE> operator()(size_type index) const
      {
        assert( (size() > 0)  && bool("Cannot apply operator() to vector of size zero!"));
        assert( index < size() && bool("Index out of bounds!") );
        
        return const_entry_proxy<SCALARTYPE>(start_ + stride_ * index, elements_);
      }
      
      const_entry_proxy<SCALARTYPE> operator[](size_type index) const
      {
        assert( (size() > 0)  && bool("Cannot apply operator() to vector of size zero!"));
        assert( index < size() && bool("Index out of bounds!") );
        
        return const_entry_proxy<SCALARTYPE>(start_ + stride_ * index, elements_);
      }
      
      //
      // Operator overloads with implicit conversion (thus cannot be made global without introducing additional headache)
      //
      
      self_type & operator += (const self_type & vec)
      {
        assert(vec.size() == size() && bool("Incompatible vector sizes!"));
  
        if (size() > 0)
          viennacl::linalg::avbv(*this, 
                                  *this, SCALARTYPE(1.0), 1, false, false,
                                  vec,   SCALARTYPE(1.0), 1, false, false);
        return *this;
      }
      
      self_type & operator -= (const self_type & vec)
      {
        assert(vec.size() == size() && bool("Incompatible vector sizes!"));
  
        if (size() > 0)
          viennacl::linalg::avbv(*this, 
                                  *this, SCALARTYPE(1.0),  1, false, false,
                                  vec,   SCALARTYPE(-1.0), 1, false, false);
        return *this;
      }
      
      self_type & operator *= (SCALARTYPE val)
      {
        if (size() > 0)
          viennacl::linalg::av(*this,
                                *this, val, 1, false, false);
        return *this;
      }
      
      self_type & operator /= (SCALARTYPE val)
      {
        if (size() > 0)
          viennacl::linalg::av(*this,
                               *this, val, 1, true, false);
        return *this;
      }
      
      
      vector_expression< const self_type, const SCALARTYPE, op_prod> 
      operator * (SCALARTYPE value) const
      {
        return vector_expression< const self_type, const SCALARTYPE, op_prod>(*this, value);
      }
  
  
      vector_expression< const self_type, const SCALARTYPE, op_prod> 
      operator / (SCALARTYPE value) const
      {
        return vector_expression< const self_type, const SCALARTYPE, op_prod>(*this, SCALARTYPE(1.0) / value);
      }
      
      
      vector_expression<const self_type, const SCALARTYPE, op_prod> operator-() const
      {
        return vector_expression<const self_type, const SCALARTYPE, op_prod>(*this, SCALARTYPE(-1.0));
      }
      
      //
      //
      
      iterator begin()
      {
        return iterator(*this, 0, start_, stride_);
      }
  
      iterator end()
      {
        return iterator(*this, size(), start_, stride_);
      }
  
      const_iterator begin() const
      {
        return const_iterator(*this, 0, start_, stride_);
      }
  
      const_iterator end() const
      {
        return const_iterator(*this, size(), start_, stride_);
      }
  
      self_type & swap(self_type & other)
      {
        viennacl::linalg::vector_swap(*this, other);
        return *this;
      };
      
      
      size_type size() const { return size_; }
      
      size_type internal_size() const { return viennacl::tools::roundUpToNextMultiple<size_type>(size_, 1); }

      size_type start() const { return start_; }
      
      size_type stride() const { return stride_; }
      
      
      bool empty() const { return size_ == 0; }
      
      const handle_type & handle() const { return elements_; }
  
      handle_type & handle() { return elements_; }
      
      void clear()
      {
        viennacl::linalg::vector_assign(*this, cpu_value_type(0.0));
      }
      
      viennacl::memory_types memory_domain() const
      {
        return elements_.get_active_handle_id();
      }
      
    protected:
      
      void set_handle(viennacl::backend::mem_handle const & h)
      {
        elements_ = h;
      }
        
      self_type & fast_swap(self_type & other) 
      { 
        assert(this->size_ == other.size_ && bool("Vector size mismatch")); 
        this->elements_.swap(other.elements_); 
        return *this; 
      }
      
      void pad()
      {
        if (internal_size() != size())
        {
          std::vector<SCALARTYPE> pad(internal_size() - size());
          viennacl::backend::memory_write(elements_, 0, sizeof(SCALARTYPE) * pad.size(), &(pad[0]));
        }
      }
      
      void switch_memory_domain(viennacl::memory_types new_domain)
      {
        viennacl::backend::switch_memory_domain<SCALARTYPE>(elements_, new_domain);
      }
      
      //TODO: Think about implementing the following public member functions
      //void insert_element(unsigned int i, SCALARTYPE val){}
      //void erase_element(unsigned int i){}
      
      //enlarge or reduce allocated memory and set unused memory to zero
      void resize(size_type new_size, bool preserve = true)
      {
        assert(new_size > 0 && bool("Positive size required when resizing vector!"));
        
        if (new_size != size_)
        {
          std::size_t new_internal_size = viennacl::tools::roundUpToNextMultiple<std::size_t>(new_size, 1);
        
          std::vector<SCALARTYPE> temp(size_);
          if (preserve && size_ > 0)
            fast_copy(*this, temp);
          temp.resize(new_size);  //drop all entries above new_size
          temp.resize(new_internal_size); //enlarge to fit new internal size
          
          if (new_internal_size != internal_size())
          {
            viennacl::backend::memory_create(elements_, sizeof(SCALARTYPE)*new_internal_size);
          }
          
          fast_copy(temp, *this);
          size_ = new_size;
        }
        
      }
      
      
    private:
      size_type       size_;
      size_type       start_;
      difference_type stride_;
      handle_type elements_;
  }; //vector_base
  
  

  // forward definition in forwards.h!
  template<class SCALARTYPE, unsigned int ALIGNMENT>
  class vector : public vector_base<SCALARTYPE>
  {
    typedef vector<SCALARTYPE, ALIGNMENT>         self_type;
    typedef vector_base<SCALARTYPE>               base_type;
    
  public:
    typedef typename base_type::size_type                  size_type;

    explicit vector() : base_type() { /* Note: One must not call ::init() here because the vector might have been created globally before the backend has become available */ }

    explicit vector(size_type vec_size) : base_type(vec_size) {}

#ifdef VIENNACL_WITH_OPENCL

    explicit vector(cl_mem existing_mem, size_type vec_size) : base_type(vec_size)
    {
      viennacl::backend::mem_handle h;
      h.switch_active_handle_id(viennacl::OPENCL_MEMORY);
      h.opencl_handle() = existing_mem;
      h.opencl_handle().inc();  //prevents that the user-provided memory is deleted once the vector object is destroyed.
      h.raw_size(sizeof(SCALARTYPE) * vec_size);
      
      base_type::set_handle(h);
    }
#endif
    
    template <typename LHS, typename RHS, typename OP>
    vector(vector_expression<const LHS, const RHS, OP> const & proxy) : base_type(proxy.size(), viennacl::traits::active_handle_id(proxy))
    {
      self_type::operator=(proxy);
    }

    vector(const base_type & v) : base_type(v.size(), v.handle().get_active_handle_id())
    {
      if (v.size() > 0)
        base_type::operator=(v);
    }

    vector(const self_type & v) : base_type(v.size(), v.handle().get_active_handle_id())
    {
      if (v.size() > 0)
        base_type::operator=(v);
    }
    
    vector(unit_vector<SCALARTYPE> const & v) : base_type(v.size())
    {
      if (v.size() > 0)
        this->operator()(v.index()) = 1;
    }
    
    vector(zero_vector<SCALARTYPE> const & v) : base_type(v.size()) 
    {
      if (v.size() > 0)
        viennacl::linalg::vector_assign(*this, SCALARTYPE(0.0));
    }

    vector(scalar_vector<SCALARTYPE> const & v) : base_type(v.size())
    {
      if (v.size() > 0)
        viennacl::linalg::vector_assign(*this, v[0]);
    }

    using base_type::operator=;
    

    // Note: Dense operations are already handled in vector_base

    //
    // Sparse matrices
    //
    template <typename SparseMatrixType>
    typename viennacl::enable_if< viennacl::is_any_sparse_matrix<SparseMatrixType>::value,
                                  self_type & >::type
    operator=(const viennacl::vector_expression< const SparseMatrixType,
                                                 const base_type,
                                                 viennacl::op_prod> & proxy) ;
    
    //
    // circulant_matrix<>
    //
    template <unsigned int MAT_ALIGNMENT>
    self_type & operator=(const vector_expression< const circulant_matrix<SCALARTYPE, MAT_ALIGNMENT>,
                                                   const base_type,
                                                   op_prod> & proxy);

    template <unsigned int MAT_ALIGNMENT>
    self_type & operator+=(const vector_expression< const circulant_matrix<SCALARTYPE, MAT_ALIGNMENT>,
                                                    const base_type,
                                                    op_prod> & proxy);
                                              
    template <unsigned int MAT_ALIGNMENT>
    self_type & operator-=(const vector_expression< const circulant_matrix<SCALARTYPE, MAT_ALIGNMENT>,
                                                    const base_type,
                                                    op_prod> & proxy);

    template <unsigned int MAT_ALIGNMENT>
    self_type operator+(const vector_expression< const circulant_matrix<SCALARTYPE, MAT_ALIGNMENT>,
                                                 const base_type,
                                                 op_prod> & proxy);

    template <unsigned int MAT_ALIGNMENT>
    self_type operator-(const vector_expression< const circulant_matrix<SCALARTYPE, MAT_ALIGNMENT>,
                                                 const base_type,
                                                 op_prod> & proxy);


    //
    // hankel_matrix<>
    //
    template <unsigned int MAT_ALIGNMENT>
    self_type & operator=(const vector_expression< const hankel_matrix<SCALARTYPE, MAT_ALIGNMENT>,
                                                   const base_type,
                                                   op_prod> & proxy);

    template <unsigned int MAT_ALIGNMENT>
    self_type & operator+=(const vector_expression< const hankel_matrix<SCALARTYPE, MAT_ALIGNMENT>,
                                                    const base_type,
                                                    op_prod> & proxy);
                                              
    template <unsigned int MAT_ALIGNMENT>
    self_type & operator-=(const vector_expression< const hankel_matrix<SCALARTYPE, MAT_ALIGNMENT>,
                                                    const base_type,
                                                    op_prod> & proxy);

    template <unsigned int MAT_ALIGNMENT>
    self_type operator+(const vector_expression< const hankel_matrix<SCALARTYPE, MAT_ALIGNMENT>,
                                                 const base_type,
                                                 op_prod> & proxy);

    template <unsigned int MAT_ALIGNMENT>
    self_type operator-(const vector_expression< const hankel_matrix<SCALARTYPE, MAT_ALIGNMENT>,
                                                 const base_type,
                                                 op_prod> & proxy);

    //
    // toeplitz_matrix<>
    //
    template <unsigned int MAT_ALIGNMENT>
    self_type & operator=(const vector_expression< const toeplitz_matrix<SCALARTYPE, MAT_ALIGNMENT>,
                                                   const base_type,
                                                   op_prod> & proxy);

    template <unsigned int MAT_ALIGNMENT>
    self_type & operator+=(const vector_expression< const toeplitz_matrix<SCALARTYPE, MAT_ALIGNMENT>,
                                                    const base_type,
                                                    op_prod> & proxy);
                                              
    template <unsigned int MAT_ALIGNMENT>
    self_type & operator-=(const vector_expression< const toeplitz_matrix<SCALARTYPE, MAT_ALIGNMENT>,
                                                    const base_type,
                                                    op_prod> & proxy);

    template <unsigned int MAT_ALIGNMENT>
    self_type operator+(const vector_expression< const toeplitz_matrix<SCALARTYPE, MAT_ALIGNMENT>,
                                                  const base_type,
                                                  op_prod> & proxy);

    template <unsigned int MAT_ALIGNMENT>
    self_type operator-(const vector_expression< const toeplitz_matrix<SCALARTYPE, MAT_ALIGNMENT>,
                                                  const base_type,
                                                  op_prod> & proxy);

    
    //
    // vandermonde_matrix<>
    //
    template <unsigned int MAT_ALIGNMENT>
    self_type & operator=(const vector_expression< const vandermonde_matrix<SCALARTYPE, MAT_ALIGNMENT>,
                                                    const base_type,
                                                    op_prod> & proxy);

    template <unsigned int MAT_ALIGNMENT>
    self_type & operator+=(const vector_expression< const vandermonde_matrix<SCALARTYPE, MAT_ALIGNMENT>,
                                                    const base_type,
                                                    op_prod> & proxy);
                                              
    template <unsigned int MAT_ALIGNMENT>
    self_type & operator-=(const vector_expression< const vandermonde_matrix<SCALARTYPE, MAT_ALIGNMENT>,
                                                    const base_type,
                                                    op_prod> & proxy);

    template <unsigned int MAT_ALIGNMENT>
    self_type operator+(const vector_expression< const vandermonde_matrix<SCALARTYPE, MAT_ALIGNMENT>,
                                                 const base_type,
                                                 op_prod> & proxy);

    template <unsigned int MAT_ALIGNMENT>
    self_type operator-(const vector_expression< const vandermonde_matrix<SCALARTYPE, MAT_ALIGNMENT>,
                                                 const base_type,
                                                 op_prod> & proxy);
    
    
    vector_expression<const vector_base<SCALARTYPE>, const SCALARTYPE, op_prod> operator-() const
    {
      return vector_expression<const vector_base<SCALARTYPE>, const SCALARTYPE, op_prod>(*this, SCALARTYPE(-1.0));
    }
    

    //enlarge or reduce allocated memory and set unused memory to zero
    void resize(size_type new_size, bool preserve = true)
    {
      base_type::resize(new_size, preserve);
    }
    

    self_type & fast_swap(self_type & other) 
    { 
      base_type::fast_swap(other);
      return *this;
    }
    
    void switch_memory_domain(viennacl::memory_types new_domain)
    {
      base_type::switch_memory_domain(new_domain);
    }
    
  }; //vector
  

  //
  //
  

  template <typename SCALARTYPE, unsigned int ALIGNMENT, typename CPU_ITERATOR>
  void fast_copy(const const_vector_iterator<SCALARTYPE, ALIGNMENT> & gpu_begin,
                  const const_vector_iterator<SCALARTYPE, ALIGNMENT> & gpu_end,
                  CPU_ITERATOR cpu_begin )
  {
    if (gpu_begin != gpu_end)
    {
      if (gpu_begin.stride() == 1)
      {
        viennacl::backend::memory_read(gpu_begin.handle(), 
                                      sizeof(SCALARTYPE)*gpu_begin.offset(),
                                      sizeof(SCALARTYPE)*gpu_begin.stride() * (gpu_end - gpu_begin),
                                      &(*cpu_begin));
      }
      else
      {
        std::size_t gpu_size = (gpu_end - gpu_begin);
        std::vector<SCALARTYPE> temp_buffer(gpu_begin.stride() * gpu_size);
        viennacl::backend::memory_read(gpu_begin.handle(), sizeof(SCALARTYPE)*gpu_begin.offset(), sizeof(SCALARTYPE)*temp_buffer.size(), &(temp_buffer[0]));

        for (std::size_t i=0; i<gpu_size; ++i)
        {
          (&(*cpu_begin))[i] = temp_buffer[i * gpu_begin.stride()];
        }
      }
    }
  }

  template <typename NumericT, typename CPUVECTOR>
  void fast_copy(vector_base<NumericT> const & gpu_vec, CPUVECTOR & cpu_vec )
  {
    viennacl::fast_copy(gpu_vec.begin(), gpu_vec.end(), cpu_vec.begin());
  }

  
  template <typename SCALARTYPE, unsigned int ALIGNMENT, typename CPU_ITERATOR>
  void copy(const const_vector_iterator<SCALARTYPE, ALIGNMENT> & gpu_begin,
            const const_vector_iterator<SCALARTYPE, ALIGNMENT> & gpu_end,
            CPU_ITERATOR cpu_begin )
  {
    assert(gpu_end - gpu_begin >= 0 && bool("Iterators incompatible"));
    if (gpu_end - gpu_begin != 0)
    {
      std::vector<SCALARTYPE> temp_buffer(gpu_end - gpu_begin);
      fast_copy(gpu_begin, gpu_end, temp_buffer.begin());
      
      //now copy entries to cpu_vec:
      std::copy(temp_buffer.begin(), temp_buffer.end(), cpu_begin);
    }
  }

  template <typename SCALARTYPE, unsigned int ALIGNMENT, typename CPU_ITERATOR>
  void copy(const vector_iterator<SCALARTYPE, ALIGNMENT> & gpu_begin,
            const vector_iterator<SCALARTYPE, ALIGNMENT> & gpu_end,
            CPU_ITERATOR cpu_begin )

  {
    viennacl::copy(const_vector_iterator<SCALARTYPE, ALIGNMENT>(gpu_begin),
                    const_vector_iterator<SCALARTYPE, ALIGNMENT>(gpu_end),
                    cpu_begin);
  }
  
  template <typename NumericT, typename CPUVECTOR>
  void copy(vector_base<NumericT> const & gpu_vec, CPUVECTOR & cpu_vec )
  {
    viennacl::copy(gpu_vec.begin(), gpu_vec.end(), cpu_vec.begin());
  }



  #ifdef VIENNACL_WITH_EIGEN
  template <unsigned int ALIGNMENT>
  void copy(vector<float, ALIGNMENT> const & gpu_vec,
            Eigen::VectorXf & eigen_vec)
  {
    viennacl::fast_copy(gpu_vec.begin(), gpu_vec.end(), &(eigen_vec[0]));
  }
  
  template <unsigned int ALIGNMENT>
  void copy(vector<double, ALIGNMENT> & gpu_vec,
            Eigen::VectorXd & eigen_vec)
  {
    viennacl::fast_copy(gpu_vec.begin(), gpu_vec.end(), &(eigen_vec[0]));
  }
  #endif


  //
  //

  template <typename CPU_ITERATOR, typename SCALARTYPE, unsigned int ALIGNMENT>
  void fast_copy(CPU_ITERATOR const & cpu_begin,
                  CPU_ITERATOR const & cpu_end,
                  vector_iterator<SCALARTYPE, ALIGNMENT> gpu_begin)
  {
    if (cpu_end - cpu_begin > 0)
    {
      if (gpu_begin.stride() == 1)
      {
        viennacl::backend::memory_write(gpu_begin.handle(),
                                        sizeof(SCALARTYPE)*gpu_begin.offset(),
                                        sizeof(SCALARTYPE)*gpu_begin.stride() * (cpu_end - cpu_begin), &(*cpu_begin));
      }
      else //writing to slice:
      {
        std::size_t cpu_size = (cpu_end - cpu_begin);
        std::vector<SCALARTYPE> temp_buffer(gpu_begin.stride() * cpu_size);
        
        viennacl::backend::memory_read(gpu_begin.handle(), sizeof(SCALARTYPE)*gpu_begin.offset(), sizeof(SCALARTYPE)*temp_buffer.size(), &(temp_buffer[0]));

        for (std::size_t i=0; i<cpu_size; ++i)
          temp_buffer[i * gpu_begin.stride()] = (&(*cpu_begin))[i];
        
        viennacl::backend::memory_write(gpu_begin.handle(), sizeof(SCALARTYPE)*gpu_begin.offset(), sizeof(SCALARTYPE)*temp_buffer.size(), &(temp_buffer[0]));
      }
    }
  }


  template <typename CPUVECTOR, typename NumericT>
  void fast_copy(const CPUVECTOR & cpu_vec, vector_base<NumericT> & gpu_vec)
  {
    viennacl::fast_copy(cpu_vec.begin(), cpu_vec.end(), gpu_vec.begin());
  }
  
  //from cpu to gpu. Safe assumption: cpu_vector does not necessarily occupy a linear memory segment, but is not larger than the allocated memory on the GPU
  template <typename SCALARTYPE, unsigned int ALIGNMENT, typename CPU_ITERATOR>
  void copy(CPU_ITERATOR const & cpu_begin,
            CPU_ITERATOR const & cpu_end,
            vector_iterator<SCALARTYPE, ALIGNMENT> gpu_begin)
  {
    assert(cpu_end - cpu_begin > 0 && bool("Iterators incompatible"));
    if (cpu_begin != cpu_end)
    {
      //we require that the size of the gpu_vector is larger or equal to the cpu-size
      std::vector<SCALARTYPE> temp_buffer(cpu_end - cpu_begin);
      std::copy(cpu_begin, cpu_end, temp_buffer.begin());
      viennacl::fast_copy(temp_buffer.begin(), temp_buffer.end(), gpu_begin);
    }
  }

  // for things like copy(std_vec.begin(), std_vec.end(), vcl_vec.begin() + 1);

  template <typename CPUVECTOR, typename T>
  void copy(const CPUVECTOR & cpu_vec, vector_base<T> & gpu_vec)
  {
    viennacl::copy(cpu_vec.begin(), cpu_vec.end(), gpu_vec.begin());
  }


  #ifdef VIENNACL_WITH_EIGEN
  template <unsigned int ALIGNMENT>
  void copy(Eigen::VectorXf const & eigen_vec,
            vector<float, ALIGNMENT> & gpu_vec)
  {
    std::vector<float> entries(eigen_vec.size());
    for (std::size_t i = 0; i<entries.size(); ++i)
      entries[i] = eigen_vec(i);
    viennacl::fast_copy(entries.begin(), entries.end(), gpu_vec.begin());
  }
  
  template <unsigned int ALIGNMENT>
  void copy(Eigen::VectorXd const & eigen_vec,
            vector<double, ALIGNMENT> & gpu_vec)
  {
    std::vector<double> entries(eigen_vec.size());
    for (std::size_t i = 0; i<entries.size(); ++i)
      entries[i] = eigen_vec(i);
    viennacl::fast_copy(entries.begin(), entries.end(), gpu_vec.begin());
  }
  #endif
  


  //
  //
  template <typename SCALARTYPE, unsigned int ALIGNMENT_SRC, unsigned int ALIGNMENT_DEST>
  void copy(const_vector_iterator<SCALARTYPE, ALIGNMENT_SRC> const & gpu_src_begin,
            const_vector_iterator<SCALARTYPE, ALIGNMENT_SRC> const & gpu_src_end,
            vector_iterator<SCALARTYPE, ALIGNMENT_DEST> gpu_dest_begin)
  {
    assert(gpu_src_end - gpu_src_begin >= 0);
    assert(gpu_src_begin.stride() == 1 && bool("ViennaCL ERROR: copy() for GPU->GPU not implemented for slices! Use operator= instead for the moment."));

    if (gpu_src_begin.stride() == 1 && gpu_dest_begin.stride() == 1)
    {
      if (gpu_src_begin != gpu_src_end)
        viennacl::backend::memory_copy(gpu_src_begin.handle(), gpu_dest_begin.handle(),
                                        sizeof(SCALARTYPE) * gpu_src_begin.offset(),
                                        sizeof(SCALARTYPE) * gpu_dest_begin.offset(),
                                        sizeof(SCALARTYPE) * (gpu_src_end.offset() - gpu_src_begin.offset()));
    }
    else
    {
      assert( false && bool("not implemented yet"));
    }
  }

  template <typename SCALARTYPE, unsigned int ALIGNMENT_SRC, unsigned int ALIGNMENT_DEST>
  void copy(vector_iterator<SCALARTYPE, ALIGNMENT_SRC> const & gpu_src_begin,
            vector_iterator<SCALARTYPE, ALIGNMENT_SRC> const & gpu_src_end,
            vector_iterator<SCALARTYPE, ALIGNMENT_DEST> gpu_dest_begin)
  {
    viennacl::copy(static_cast<const_vector_iterator<SCALARTYPE, ALIGNMENT_SRC> >(gpu_src_begin),
                    static_cast<const_vector_iterator<SCALARTYPE, ALIGNMENT_SRC> >(gpu_src_end),
                    gpu_dest_begin);
  }

  template <typename SCALARTYPE, unsigned int ALIGNMENT_SRC, unsigned int ALIGNMENT_DEST>
  void copy(vector<SCALARTYPE, ALIGNMENT_SRC> const & gpu_src_vec,
            vector<SCALARTYPE, ALIGNMENT_DEST> & gpu_dest_vec )
  {
    viennacl::copy(gpu_src_vec.begin(), gpu_src_vec.end(), gpu_dest_vec.begin());
  } 


  
  
  

  //global functions for handling vectors:
  template <typename T>
  std::ostream & operator<<(std::ostream & s, vector_base<T> const & val)
  {
    std::vector<T> tmp(val.size());
    viennacl::copy(val.begin(), val.end(), tmp.begin());
    std::cout << "[" << val.size() << "](";
    for (typename std::vector<T>::size_type i=0; i<val.size(); ++i)
    {
      if (i > 0)
        s << ",";
      s << tmp[i];
    }
    std::cout << ")";
    return s;
  }

  template <typename LHS, typename RHS, typename OP>
  std::ostream & operator<<(std::ostream & s, vector_expression<LHS, RHS, OP> const & proxy)
  {
    typedef typename viennacl::result_of::cpu_value_type<typename LHS::value_type>::type ScalarType;
    viennacl::vector<ScalarType> result = proxy;
    s << result;
    return s;
  }
  
  template <typename T>
  void swap(vector_base<T> & vec1, vector_base<T> & vec2)
  {
    viennacl::linalg::vector_swap(vec1, vec2);
  }
  
  template <typename SCALARTYPE, unsigned int ALIGNMENT>
  vector<SCALARTYPE, ALIGNMENT> & fast_swap(vector<SCALARTYPE, ALIGNMENT> & v1,
                                            vector<SCALARTYPE, ALIGNMENT> & v2) 
  { 
    return v1.fast_swap(v2);
  }       
  
  
  
  
  
  //
  //
  //
  //
  
  //
  // operator +=
  //
  
                                              
  template <typename T>
  vector_base<T> & operator += (vector_base<T> & v1, const vector_base<T> & v2)
  {
    
    assert(v1.size() == v2.size() && bool("Incompatible vector sizes!"));

    if (v1.size() > 0)
      viennacl::linalg::avbv(v1, 
                             v1, T(1.0), 1, false, false,
                             v2, T(1.0), 1, false, false);
    return v1;
  }
  
  template <typename T, typename S2, typename OP>
  typename viennacl::enable_if< viennacl::is_any_scalar<S2>::value,
                                vector_base<T> &>::type
  operator += (vector_base<T> & v1, 
                const vector_expression< const vector_base<T>, const S2, OP> & proxy)
  {
    assert(proxy.lhs().size() == v1.size() && bool("Incompatible vector sizes!"));

    if (v1.size() > 0)
      viennacl::linalg::avbv(v1, 
                              v1,               T(1.0), 1, false,                                             false,
                              proxy.lhs(), proxy.rhs(), 1, (viennacl::is_division<OP>::value ? true : false), (viennacl::is_flip_sign_scalar<S2>::value ? true : false) );
    return v1;
  }

  
  template <typename T, typename OP>
  typename viennacl::enable_if< (viennacl::is_addition<OP>::value || viennacl::is_subtraction<OP>::value),
                                vector_base<T> &>::type
  operator += (vector_base<T> & v1, 
               const vector_expression< const vector_base<T>, const vector_base<T>, OP> & proxy)
  {
    assert(proxy.lhs().size() == v1.size() && bool("Incompatible vector sizes!"));

    if (v1.size() > 0)
      viennacl::linalg::avbv_v(v1, 
                                proxy.lhs(), T(1.0), 1, false, false,
                                proxy.rhs(), T(1.0), 1, false, (viennacl::is_subtraction<OP>::value ? true : false) );
    return v1;
  }
  
  template < typename T, 
             typename S3, typename OP3,
             typename OP>
  typename viennacl::enable_if< viennacl::is_any_scalar<S3>::value && (viennacl::is_product<OP3>::value || viennacl::is_division<OP3>::value)
                                && (viennacl::is_addition<OP>::value || viennacl::is_subtraction<OP>::value),
                                vector_base<T> &>::type
  operator += (vector_base<T> & v1,
               const vector_expression< const vector_base<T>,
                                        const vector_expression<const vector_base<T>, const S3, OP3>,
                                        OP> & proxy)
  {
    assert(proxy.lhs().size() == v1.size() && bool("Incompatible vector sizes!"));

    if (v1.size() > 0)
    {
      bool flip_sign_3 = (viennacl::is_subtraction<OP>::value ? true : false);
      if (viennacl::is_flip_sign_scalar<S3>::value)
        flip_sign_3 = !flip_sign_3;
      viennacl::linalg::avbv_v(v1, 
                                proxy.lhs(),                  T(1.0), 1, false                                             , false,
                                proxy.rhs().lhs(), proxy.rhs().rhs(), 1, (viennacl::is_division<OP3>::value ? true : false), flip_sign_3 );
    }
    return v1;
  }

  template <typename T, typename S2, typename OP2,
            typename OP>
  typename viennacl::enable_if< viennacl::is_any_scalar<S2>::value && (viennacl::is_product<OP2>::value || viennacl::is_division<OP2>::value)
                                && (viennacl::is_addition<OP>::value || viennacl::is_subtraction<OP>::value),
                                vector_base<T> &>::type
  operator += (vector_base<T> & v1,
               const vector_expression< const vector_expression<const vector_base<T>, const S2, OP2>,
                                        const vector_base<T>,
                                        OP> & proxy)
  {
    assert(proxy.size() == v1.size() && bool("Incompatible vector sizes!"));

    if (v1.size() > 0)
      viennacl::linalg::avbv_v(v1, 
                                proxy.lhs().lhs(), proxy.lhs().rhs(), 1, (viennacl::is_division<OP2>::value ? true : false), (viennacl::is_flip_sign_scalar<S2>::value ? true : false),
                                proxy.rhs(),                  T(1.0), 1, false                                             , (viennacl::is_subtraction<OP>::value ? true : false) );
    return v1;
  }
  
  template <typename T, 
            typename S2, typename OP2,
            typename S3, typename OP3,
            typename OP>
  typename viennacl::enable_if<    viennacl::is_any_scalar<S2>::value && (viennacl::is_product<OP2>::value || viennacl::is_division<OP2>::value)
                                && viennacl::is_any_scalar<S3>::value && (viennacl::is_product<OP3>::value || viennacl::is_division<OP3>::value)
                                && (viennacl::is_addition<OP>::value || viennacl::is_subtraction<OP>::value),
                                vector_base<T> &>::type
  operator += (vector_base<T> & v1,
               const vector_expression< const vector_expression<const vector_base<T>, const S2, OP2>,
                                        const vector_expression<const vector_base<T>, const S3, OP3>,
                                        OP> & proxy)
  {
    assert(proxy.lhs().size() == v1.size() && bool("Incompatible vector sizes!"));

    if (v1.size() > 0)
    {
      bool flip_sign_3 = (viennacl::is_subtraction<OP>::value ? true : false);
      if (viennacl::is_flip_sign_scalar<S3>::value)
        flip_sign_3 = !flip_sign_3;
      viennacl::linalg::avbv_v(v1, 
                                proxy.lhs().lhs(), proxy.lhs().rhs(), 1, (viennacl::is_division<OP2>::value ? true : false), (viennacl::is_flip_sign_scalar<S2>::value ? true : false),
                                proxy.rhs().lhs(), proxy.rhs().rhs(), 1, (viennacl::is_division<OP3>::value ? true : false), flip_sign_3 );
    }
    return v1;
  }
  
  
  //
  // operator -=
  //
  
  template <typename T>
  vector_base<T> & operator -= (vector_base<T> & v1, const vector_base<T> & vec)
  {
    assert(vec.size() == v1.size() && bool("Incompatible vector sizes!"));

    if (v1.size() > 0)
      viennacl::linalg::avbv(v1, 
                             v1,  T( 1.0), 1, false, false,
                             vec, T(-1.0), 1, false, false);
    return v1;
  }

  
  template <typename T, typename S2, typename OP>
  typename viennacl::enable_if< viennacl::is_any_scalar<S2>::value,
                                vector_base<T> &>::type
  operator -= (vector_base<T> & v1, 
               const vector_expression< const vector_base<T>, const S2, OP> & proxy)
  {
    assert(proxy.lhs().size() == v1.size() && bool("Incompatible vector sizes!"));

    if (v1.size() > 0)
      viennacl::linalg::avbv(v1, 
                             v1,               T(1.0), 1, false,                                             false,
                             proxy.lhs(), proxy.rhs(), 1, (viennacl::is_division<OP>::value ? true : false), (viennacl::is_flip_sign_scalar<S2>::value ? false : true));
    return v1;
  }
  
  template <typename T, typename OP>
  typename viennacl::enable_if< (viennacl::is_addition<OP>::value || viennacl::is_subtraction<OP>::value),
                                vector_base<T> &>::type
  operator -= (vector_base<T> & v1, 
               const vector_expression< const vector_base<T>, const vector_base<T>, OP> & proxy)
  {
    assert(proxy.lhs().size() == v1.size() && bool("Incompatible vector sizes!"));

    if (v1.size() > 0)
      viennacl::linalg::avbv_v(v1, 
                               proxy.lhs(), T(1.0), 1, false, true,
                               proxy.rhs(), T(1.0), 1, false, (viennacl::is_subtraction<OP>::value ? false : true) );
    return v1;
  }
  
  template <typename T,
            typename S3, typename OP3,
            typename OP>
  typename viennacl::enable_if< viennacl::is_any_scalar<S3>::value && (viennacl::is_product<OP3>::value || viennacl::is_division<OP3>::value)
                                && (viennacl::is_addition<OP>::value || viennacl::is_subtraction<OP>::value),
                                vector_base<T> &>::type
  operator -= (vector_base<T> & v1, 
               const vector_expression< const vector_base<T>,
                                        const vector_expression<const vector_base<T>, const S3, OP3>,
                                        OP> & proxy)
  {
    assert(proxy.lhs().size() == v1.size() && bool("Incompatible vector sizes!"));

    if (v1.size() > 0)
    {
      bool flip_sign_3 = (viennacl::is_subtraction<OP>::value ? false : true);
      if (viennacl::is_flip_sign_scalar<S3>::value)
        flip_sign_3 = !flip_sign_3;
      viennacl::linalg::avbv_v(v1, 
                                proxy.lhs(),                  T(1.0), 1, false                                             , true,
                                proxy.rhs().lhs(), proxy.rhs().rhs(), 1, (viennacl::is_division<OP3>::value ? true : false), flip_sign_3);
    }
    return v1;
  }

  template <typename T, typename S2, typename OP2,
            typename OP>
  typename viennacl::enable_if< viennacl::is_any_scalar<S2>::value && (viennacl::is_product<OP2>::value || viennacl::is_division<OP2>::value)
                                && (viennacl::is_addition<OP>::value || viennacl::is_subtraction<OP>::value),
                                vector_base<T> &>::type
  operator -= (vector_base<T> & v1, 
               const vector_expression< const vector_expression<const vector_base<T>, const S2, OP2>,
                                        const vector_base<T>,
                                        OP> & proxy)
  {
    assert(proxy.size() == v1.size() && bool("Incompatible vector sizes!"));

    if (v1.size() > 0)
      viennacl::linalg::avbv_v(v1, 
                                proxy.lhs().lhs(), proxy.lhs().rhs(), 1, (viennacl::is_division<OP2>::value ? true : false), (viennacl::is_flip_sign_scalar<S2>::value ? false : true),
                                proxy.rhs(),                  T(1.0), 1, false                                             , (viennacl::is_subtraction<OP>::value ? false : true) );
    return v1;
  }
  
  template <typename T,
            typename S2, typename OP2,
            typename S3, typename OP3,
            typename OP>
  typename viennacl::enable_if<    viennacl::is_any_scalar<S2>::value && (viennacl::is_product<OP2>::value || viennacl::is_division<OP2>::value)
                                && viennacl::is_any_scalar<S3>::value && (viennacl::is_product<OP3>::value || viennacl::is_division<OP3>::value)
                                && (viennacl::is_addition<OP>::value || viennacl::is_subtraction<OP>::value),
                                vector_base<T> &>::type
  operator -= (vector_base<T> & v1, 
               const vector_expression< const vector_expression<const vector_base<T>, const S2, OP2>,
                                        const vector_expression<const vector_base<T>, const S3, OP3>,
                                        OP> & proxy)
  {
    assert(proxy.lhs().size() == v1.size() && bool("Incompatible vector sizes!"));

    if (v1.size() > 0)
    {
      bool flip_sign_3 = (viennacl::is_subtraction<OP>::value ? false : true);
      if (viennacl::is_flip_sign_scalar<S3>::value)
        flip_sign_3 = !flip_sign_3;
      viennacl::linalg::avbv_v(v1, 
                                proxy.lhs().lhs(), proxy.lhs().rhs(), 1, (viennacl::is_division<OP2>::value ? true : false), (viennacl::is_flip_sign_scalar<S2>::value ? false : true),
                                proxy.rhs().lhs(), proxy.rhs().rhs(), 1, (viennacl::is_division<OP3>::value ? true : false), flip_sign_3);
    }
    return v1;
  }
  
  
  //
  // operator *=
  //

  template <typename T, typename S1>
  typename viennacl::enable_if< viennacl::is_any_scalar<S1>::value,
                                vector_base<T> & 
                              >::type
  operator *= (vector_base<T> & v1, S1 const & gpu_val)
  {
    if (v1.size() > 0)
      viennacl::linalg::av(v1,
                           v1, gpu_val, 1, false, (viennacl::is_flip_sign_scalar<S1>::value ? true : false));
    return v1;
  }

  
  //
  // operator /=
  //
    

  template <typename T, typename S1>
  typename viennacl::enable_if< viennacl::is_any_scalar<S1>::value,
                                vector_base<T> & 
                              >::type
  operator /= (vector_base<T> & v1, S1 const & gpu_val)
  {
    if (v1.size() > 0)
      viennacl::linalg::av(v1,
                           v1, gpu_val, 1, true, (viennacl::is_flip_sign_scalar<S1>::value ? true : false));
    return v1;
  }
  
  
  //
  // operator +
  //
  
  //addition and subtraction of two vector_expressions:
  template <typename LHS1, typename RHS1, typename OP1,
            typename LHS2, typename RHS2, typename OP2>
  typename vector_expression< LHS1, RHS1, OP1>::VectorType
  operator + (vector_expression< LHS1, RHS1, OP1> const & proxy1,
              vector_expression< LHS2, RHS2, OP2> const & proxy2)
  {
    assert(proxy1.size() == proxy2.size() && bool("Incompatible vector sizes!"));
    typename vector_expression< LHS1, RHS1, OP1>::VectorType result = proxy1;
    result += proxy2;
    return result;
  }
  
  
  template <typename LHS, typename RHS, typename OP, typename T>
  viennacl::vector<T>
  operator + (vector_expression<LHS, RHS, OP> const & proxy,
              vector_base<T> const & vec)
  {
    assert(proxy.size() == vec.size() && bool("Incompatible vector sizes!"));
    viennacl::vector<T> result = proxy;
    result += vec;
    return result;
  }

  template <typename T, typename LHS, typename RHS, typename OP>
  viennacl::vector<T>
  operator + (vector_base<T> const & vec,
              vector_expression<LHS, RHS, OP> const & proxy)
  {
    assert(proxy.size() == vec.size() && bool("Incompatible vector sizes!"));
    viennacl::vector<T> result = vec;
    result += proxy;
    return result;
  }

  
  template <typename T, typename S1, typename OP1>
  typename viennacl::enable_if< viennacl::is_any_scalar<S1>::value,
                                vector_expression<const vector_expression<const vector_base<T>, const S1, OP1>,
                                                  const vector_base<T>,
                                                  op_add>
                              >::type
  operator + (vector_expression<const vector_base<T>, const S1, OP1> const & proxy,
              vector_base<T> const & vec)
  {
    return vector_expression<const vector_expression<const vector_base<T>, const S1, OP1>,
                             const vector_base<T>,
                             op_add>(proxy, vec);
  }
  
  template <typename T,
            typename S1, typename OP1,
            typename S2, typename OP2>
  typename viennacl::enable_if< viennacl::is_any_scalar<S1>::value && viennacl::is_any_scalar<S2>::value,
                                vector_expression<const vector_expression<const vector_base<T>, const S1, OP1>,
                                                  const vector_expression<const vector_base<T>, const S2, OP2>,
                                                  op_add>
                              >::type
  operator + (vector_expression<const vector_base<T>, const S1, OP1> const & lhs,
              vector_expression<const vector_base<T>, const S2, OP2> const & rhs)
  {
    return vector_expression<const vector_expression<const vector_base<T>, const S1, OP1>,
                             const vector_expression<const vector_base<T>, const S2, OP2>,
                             op_add>(lhs, rhs);
  }
  
  
  template <typename T>
  vector_expression< const vector_base<T>, const vector_base<T>, op_add>      
  operator+(const vector_base<T> & v1, const vector_base<T> & v2)
  {
    return vector_expression< const vector_base<T>, const vector_base<T>, op_add>(v1, v2);
  }
  
  template <typename T, typename S2, typename OP2>
  typename viennacl::enable_if< viennacl::is_any_scalar<S2>::value,
                                vector_expression< const vector_base<T>,
                                                    const vector_expression< const vector_base<T>, const S2, OP2>,
                                                    op_add>      
                              >::type
  operator+(const vector_base<T> & v1, 
            const vector_expression< const vector_base<T>,
                                     const S2,
                                     OP2> & proxy)
  {
    return vector_expression< const vector_base<T>,
                              const vector_expression< const vector_base<T>, const S2, OP2>,
                              op_add>(v1, proxy);
  }


  
  
  //
  // operator -
  //
  
  template <typename LHS1, typename RHS1, typename OP1,
            typename LHS2, typename RHS2, typename OP2>
  typename vector_expression< LHS1, RHS1, OP1>::VectorType
  operator - (vector_expression< LHS1, RHS1, OP1> const & proxy1,
              vector_expression< LHS2, RHS2, OP2> const & proxy2)
  {
    assert(proxy1.size() == proxy2.size() && bool("Incompatible vector sizes!"));
    typename vector_expression< LHS1, RHS1, OP1>::VectorType result = proxy1;
    result -= proxy2;
    return result;
  }
  
  
  template <typename LHS, typename RHS, typename OP, typename T>
  viennacl::vector<T>
  operator - (vector_expression< LHS, RHS, OP> const & proxy,
              vector_base<T> const & vec)
  {
    assert(proxy.size() == vec.size() && bool("Incompatible vector sizes!"));
    viennacl::vector<T> result = proxy;
    result -= vec;
    return result;
  }

  template <typename T, typename LHS, typename RHS, typename OP>
  viennacl::vector<T>
  operator - (vector_base<T> const & vec,
              vector_expression< LHS, RHS, OP> const & proxy)
  {
    assert(proxy.size() == vec.size() && bool("Incompatible vector sizes!"));
    viennacl::vector<T> result = vec;
    result -= proxy;
    return result;
  }
  
  template <typename T, typename S1, typename OP1>
  typename viennacl::enable_if< viennacl::is_any_scalar<S1>::value,
                                vector_expression<const vector_expression<const vector_base<T>, const S1, OP1>,
                                                  const vector_base<T>,
                                                  op_sub>
                              >::type
  operator - (vector_expression<const vector_base<T>, const S1, OP1> const & proxy,
              vector_base<T> const & vec)
  {
    return vector_expression<const vector_expression<const vector_base<T>, const S1, OP1>,
                              const vector_base<T>,
                              op_sub>(proxy, vec);
  }
  
  template <typename T,
            typename S1, typename OP1,
            typename S2, typename OP2>
  typename viennacl::enable_if<    viennacl::is_any_scalar<S1>::value
                                && viennacl::is_any_scalar<S2>::value,
                                vector_expression<const vector_expression<const vector_base<T>, const S1, OP1>,
                                                  const vector_expression<const vector_base<T>, const S2, OP2>,
                                                  op_sub>
                              >::type
  operator - (vector_expression<const vector_base<T>, const S1, OP1> const & lhs,
              vector_expression<const vector_base<T>, const S2, OP2> const & rhs)
  {
    return vector_expression<const vector_expression<const vector_base<T>, const S1, OP1>,
                             const vector_expression<const vector_base<T>, const S2, OP2>,
                             op_sub>(lhs, rhs);
  }

  template <typename T>
  vector_expression< const vector_base<T>, const vector_base<T>, op_sub>      
  operator-(const vector_base<T> & v1, const vector_base<T> & v2)
  {
    return vector_expression< const vector_base<T>, const vector_base<T>, op_sub>(v1, v2);
  }


  template <typename T, typename S2, typename OP2>
  typename viennacl::enable_if< viennacl::is_any_scalar<S2>::value,
                                vector_expression< const vector_base<T>,
                                                   const vector_expression< const vector_base<T>, const S2, OP2>,
                                                   op_sub>      
                              >::type
  operator-(const vector_base<T> & v1, 
            const vector_expression< const vector_base<T>,
                                     const S2,
                                     OP2> & proxy)
  {
    return vector_expression< const vector_base<T>,
                              const vector_expression< const vector_base<T>, const S2, OP2>,
                              op_sub>(v1, proxy);
  }

  
  //
  // operator *
  //
  
  
  template <typename S1, typename T>
  typename viennacl::enable_if< viennacl::is_any_scalar<S1>::value,
                                vector_expression< const vector_base<T>, const S1, op_prod> >::type 
  operator * (S1 const & value, vector_base<T> const & vec)
  {
    return vector_expression< const vector_base<T>, const S1, op_prod>(vec, value);
  }


  template <typename LHS, typename RHS, typename OP, typename S1>
  typename viennacl::enable_if< viennacl::is_any_scalar<S1>::value,
                                viennacl::vector<typename viennacl::result_of::cpu_value_type<RHS>::type> >::type
  operator * (vector_expression< LHS, RHS, OP> const & proxy,
              S1 const & val)
  {
    viennacl::vector<typename viennacl::result_of::cpu_value_type<RHS>::type> result = proxy;
    result *= val;
    return result;
  }


  template <typename S1, typename LHS, typename RHS, typename OP>
  typename viennacl::enable_if< viennacl::is_any_scalar<S1>::value,
                                viennacl::vector<typename viennacl::result_of::cpu_value_type<RHS>::type> >::type
  operator * (S1 const & val,
              vector_expression< LHS, RHS, OP> const & proxy)
  {
    viennacl::vector<typename viennacl::result_of::cpu_value_type<RHS>::type> result = proxy;
    result *= val;
    return result;
  }
  
  template <typename T, typename S1>
  typename viennacl::enable_if< viennacl::is_any_scalar<S1>::value,
                                vector_expression< const vector_base<T>, const S1, op_prod> >::type
  operator * (vector_base<T> const & v1, S1 const & s1)
  {
    return vector_expression< const vector_base<T>, const S1, op_prod>(v1, s1);
  }
  
  //
  // operator /
  //  
  
  template <typename S1, typename LHS, typename RHS, typename OP>
  typename viennacl::enable_if< viennacl::is_any_scalar<S1>::value,
                                viennacl::vector<typename viennacl::result_of::cpu_value_type<RHS>::type> >::type
  operator / (vector_expression< LHS, RHS, OP> const & proxy,
              S1 const & val)
  {
    viennacl::vector<typename viennacl::result_of::cpu_value_type<RHS>::type> result = proxy;
    result /= val;
    return result;
  }


  template <typename T, typename S1>
  typename viennacl::enable_if< viennacl::is_any_scalar<S1>::value,
                                vector_expression< const vector_base<T>, const S1, op_div> >::type
  operator / (vector_base<T> const & v1, S1 const & s1)
  {
    return vector_expression< const vector_base<T>, const S1, op_div>(v1, s1);
  }
  
}

#endif