scalar.hpp

Go to the documentation of this file.

#ifndef VIENNACL_SCALAR_HPP_
#define VIENNACL_SCALAR_HPP_

/* =========================================================================
   Copyright (c) 2010-2015, 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 manual)

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

#include <iostream>

#include "viennacl/forwards.h"
#include "viennacl/backend/memory.hpp"
#include "viennacl/meta/result_of.hpp"
#include "viennacl/linalg/scalar_operations.hpp"
#include "viennacl/traits/handle.hpp"

#ifdef VIENNACL_WITH_OPENCL
#include "viennacl/ocl/backend.hpp"
#endif

namespace viennacl
{
template<typename LHS, typename RHS, typename OP>
class scalar_expression
{
  typedef typename LHS::value_type          DummyType; //Visual C++ 2005 does not allow to write LHS::value_type::value_type
public:
  typedef typename viennacl::result_of::cpu_value_type<DummyType>::type    ScalarType;

  scalar_expression(LHS & lhs, RHS & rhs) : lhs_(lhs), rhs_(rhs) {}

  LHS & lhs() const { return lhs_; }
  RHS & rhs() const { return rhs_; }

  operator ScalarType () const
  {
    viennacl::scalar<ScalarType> temp;
    temp = *this;
    return temp;
  }

private:
  LHS & lhs_;
  RHS & rhs_;
};


template<typename LHS, typename RHS>
class scalar_expression<LHS, RHS, op_inner_prod>
{
  //typedef typename LHS::value_type          DummyType; //Visual C++ 2005 does not allow to write LHS::value_type::value_type
public:
  typedef typename viennacl::result_of::cpu_value_type<LHS>::type    ScalarType;

  scalar_expression(LHS & lhs, RHS & rhs) : lhs_(lhs), rhs_(rhs) {}

  LHS & lhs() const { return lhs_; }
  RHS & rhs() const { return rhs_; }

  operator ScalarType () const
  {
    ScalarType result;
    viennacl::linalg::inner_prod_cpu(lhs_, rhs_, result);
    return result;
  }

private:
  LHS & lhs_;
  RHS & rhs_;
};


template<typename LHS, typename RHS>
class scalar_expression<LHS, RHS, op_norm_1>
{
  //typedef typename LHS::value_type          DummyType; //Visual C++ 2005 does not allow to write LHS::value_type::value_type
public:
  typedef typename viennacl::result_of::cpu_value_type<LHS>::type    ScalarType;

  scalar_expression(LHS & lhs, RHS & rhs) : lhs_(lhs), rhs_(rhs) {}

  LHS & lhs() const { return lhs_; }
  RHS & rhs() const { return rhs_; }

  operator ScalarType () const
  {
    ScalarType result;
    viennacl::linalg::norm_1_cpu(lhs_, result);
    return result;
  }

private:
  LHS & lhs_;
  RHS & rhs_;
};

template<typename LHS, typename RHS>
class scalar_expression<LHS, RHS, op_norm_2>
{
  //typedef typename LHS::value_type          DummyType; //Visual C++ 2005 does not allow to write LHS::value_type::value_type
public:
  typedef typename viennacl::result_of::cpu_value_type<LHS>::type    ScalarType;

  scalar_expression(LHS & lhs, RHS & rhs) : lhs_(lhs), rhs_(rhs) {}

  LHS & lhs() const { return lhs_; }
  RHS & rhs() const { return rhs_; }

  operator ScalarType () const
  {
    ScalarType result;
    viennacl::linalg::norm_2_cpu(lhs_, result);
    return result;
  }

private:
  LHS & lhs_;
  RHS & rhs_;
};


template<typename LHS, typename RHS>
class scalar_expression<LHS, RHS, op_norm_inf>
{
  //typedef typename LHS::value_type          DummyType; //Visual C++ 2005 does not allow to write LHS::value_type::value_type
public:
  typedef typename viennacl::result_of::cpu_value_type<LHS>::type    ScalarType;

  scalar_expression(LHS & lhs, RHS & rhs) : lhs_(lhs), rhs_(rhs) {}

  LHS & lhs() const { return lhs_; }
  RHS & rhs() const { return rhs_; }

  operator ScalarType () const
  {
    ScalarType result;
    viennacl::linalg::norm_inf_cpu(lhs_, result);
    return result;
  }

private:
  LHS & lhs_;
  RHS & rhs_;
};

template<typename LHS, typename RHS>
class scalar_expression<LHS, RHS, op_max>
{
  //typedef typename LHS::value_type          DummyType; //Visual C++ 2005 does not allow to write LHS::value_type::value_type
public:
  typedef typename viennacl::result_of::cpu_value_type<LHS>::type    ScalarType;

  scalar_expression(LHS & lhs, RHS & rhs) : lhs_(lhs), rhs_(rhs) {}

  LHS & lhs() const { return lhs_; }
  RHS & rhs() const { return rhs_; }

  operator ScalarType () const
  {
    ScalarType result;
    viennacl::linalg::max_cpu(lhs_, result);
    return result;
  }

private:
  LHS & lhs_;
  RHS & rhs_;
};


template<typename LHS, typename RHS>
class scalar_expression<LHS, RHS, op_min>
{
  //typedef typename LHS::value_type          DummyType; //Visual C++ 2005 does not allow to write LHS::value_type::value_type
public:
  typedef typename viennacl::result_of::cpu_value_type<LHS>::type    ScalarType;

  scalar_expression(LHS & lhs, RHS & rhs) : lhs_(lhs), rhs_(rhs) {}

  LHS & lhs() const { return lhs_; }
  RHS & rhs() const { return rhs_; }

  operator ScalarType () const
  {
    ScalarType result;
    viennacl::linalg::min_cpu(lhs_, result);
    return result;
  }

private:
  LHS & lhs_;
  RHS & rhs_;
};

template<typename LHS, typename RHS>
class scalar_expression<LHS, RHS, op_sum>
{
  //typedef typename LHS::value_type          DummyType; //Visual C++ 2005 does not allow to write LHS::value_type::value_type
public:
  typedef typename viennacl::result_of::cpu_value_type<LHS>::type    ScalarType;

  scalar_expression(LHS & lhs, RHS & rhs) : lhs_(lhs), rhs_(rhs) {}

  LHS & lhs() const { return lhs_; }
  RHS & rhs() const { return rhs_; }

  operator ScalarType () const
  {
    ScalarType result;
    viennacl::linalg::sum_cpu(lhs_, result);
    return result;
  }

private:
  LHS & lhs_;
  RHS & rhs_;
};


template<typename LHS, typename RHS>
class scalar_expression<LHS, RHS, op_norm_frobenius>
{
  //typedef typename LHS::value_type          DummyType; //Visual C++ 2005 does not allow to write LHS::value_type::value_type
public:
  typedef typename viennacl::result_of::cpu_value_type<LHS>::type    ScalarType;

  scalar_expression(LHS & lhs, RHS & rhs) : lhs_(lhs), rhs_(rhs) {}

  LHS & lhs() const { return lhs_; }
  RHS & rhs() const { return rhs_; }

  operator ScalarType () const
  {
    ScalarType result;
    viennacl::linalg::norm_frobenius_cpu(lhs_, result);
    return result;
  }

private:
  LHS & lhs_;
  RHS & rhs_;
};




template<class NumericT>
class scalar
{
  typedef scalar<NumericT>         self_type;
public:
  typedef viennacl::backend::mem_handle                     handle_type;
  typedef vcl_size_t                                        size_type;

  typedef NumericT   value_type;

  scalar() {}

  scalar(NumericT val, viennacl::context ctx = viennacl::context())
  {
    viennacl::backend::memory_create(val_, sizeof(NumericT), ctx, &val);
  }

#ifdef VIENNACL_WITH_OPENCL

  explicit scalar(cl_mem mem, size_type /*size*/)
  {
    val_.switch_active_handle_id(viennacl::OPENCL_MEMORY);
    val_.opencl_handle() = mem;
    val_.opencl_handle().inc();  //prevents that the user-provided memory is deleted once the vector object is destroyed.
  }
#endif

  template<typename T1, typename T2, typename OP>
  scalar(scalar_expression<T1, T2, OP> const & proxy)
  {
    val_.switch_active_handle_id(viennacl::traits::handle(proxy.lhs()).get_active_handle_id());
    viennacl::backend::memory_create(val_, sizeof(NumericT), viennacl::traits::context(proxy));
    *this = proxy;
  }

  //copy constructor
  scalar(const scalar & other)
  {
    if (other.handle().get_active_handle_id() != viennacl::MEMORY_NOT_INITIALIZED)
    {
      //copy value:
      val_.switch_active_handle_id(other.handle().get_active_handle_id());
      viennacl::backend::memory_create(val_, sizeof(NumericT), viennacl::traits::context(other));
      viennacl::backend::memory_copy(other.handle(), val_, 0, 0, sizeof(NumericT));
    }
  }

  operator NumericT() const
  {
    // make sure the scalar contains reasonable data:
    assert( val_.get_active_handle_id() != viennacl::MEMORY_NOT_INITIALIZED && bool("Scalar not initialized, cannot read!"));

    NumericT tmp;
    viennacl::backend::memory_read(val_, 0, sizeof(NumericT), &tmp);
    return tmp;
  }

  self_type & operator= (entry_proxy<NumericT> const & other)
  {
    init_if_necessary(viennacl::traits::context(other));
    viennacl::backend::memory_copy(other.handle(), val_, other.index() * sizeof(NumericT), 0, sizeof(NumericT));
    return *this;
  }

  self_type & operator= (scalar<NumericT> const & other)
  {
    init_if_necessary(viennacl::traits::context(other));
    viennacl::backend::memory_copy(other.handle(), val_, 0, 0, sizeof(NumericT));
    return *this;
  }

  self_type & operator= (float cpu_other)
  {
    init_if_necessary(viennacl::context());

    //copy value:
    NumericT value = static_cast<NumericT>(cpu_other);
    viennacl::backend::memory_write(val_, 0, sizeof(NumericT), &value);
    return *this;
  }

  self_type & operator= (double cpu_other)
  {
    init_if_necessary(viennacl::context());

    NumericT value = static_cast<NumericT>(cpu_other);
    viennacl::backend::memory_write(val_, 0, sizeof(NumericT), &value);
    return *this;
  }

  self_type & operator= (long cpu_other)
  {
    init_if_necessary(viennacl::context());

    NumericT value = static_cast<NumericT>(cpu_other);
    viennacl::backend::memory_write(val_, 0, sizeof(NumericT), &value);
    return *this;
  }

  self_type & operator= (unsigned long cpu_other)
  {
    init_if_necessary(viennacl::context());

    NumericT value = static_cast<NumericT>(cpu_other);
    viennacl::backend::memory_write(val_, 0, sizeof(NumericT), &value);
    return *this;
  }

  self_type & operator= (int cpu_other)
  {
    init_if_necessary(viennacl::context());

    NumericT value = static_cast<NumericT>(cpu_other);
    viennacl::backend::memory_write(val_, 0, sizeof(NumericT), &value);
    return *this;
  }

  self_type & operator= (unsigned int cpu_other)
  {
    init_if_necessary(viennacl::context());

    NumericT value = static_cast<NumericT>(cpu_other);
    viennacl::backend::memory_write(val_, 0, sizeof(NumericT), &value);
    return *this;
  }

  template<typename T1, typename T2>
  self_type & operator= (scalar_expression<T1, T2, op_inner_prod> const & proxy)
  {
    init_if_necessary(viennacl::traits::context(proxy));

    viennacl::linalg::inner_prod_impl(proxy.lhs(), proxy.rhs(), *this);
    return *this;
  }

  template<typename T1, typename T2>
  self_type & operator= (scalar_expression<T1, T2, op_norm_1> const & proxy)
  {
    init_if_necessary(viennacl::traits::context(proxy));

    viennacl::linalg::norm_1_impl(proxy.lhs(), *this);
    return *this;
  }

  template<typename T1, typename T2>
  self_type & operator= (scalar_expression<T1, T2, op_norm_2> const & proxy)
  {
    init_if_necessary(viennacl::traits::context(proxy));

    viennacl::linalg::norm_2_impl(proxy.lhs(), *this);
    return *this;
  }

  template<typename T1, typename T2>
  self_type & operator= (scalar_expression<T1, T2, op_norm_inf> const & proxy)
  {
    init_if_necessary(viennacl::traits::context(proxy));

    viennacl::linalg::norm_inf_impl(proxy.lhs(), *this);
    return *this;
  }

  template<typename T1, typename T2>
  self_type & operator= (scalar_expression<T1, T2, op_max> const & proxy)
  {
    init_if_necessary(viennacl::traits::context(proxy));

    viennacl::linalg::max_impl(proxy.lhs(), *this);
    return *this;
  }

  template<typename T1, typename T2>
  self_type & operator= (scalar_expression<T1, T2, op_min> const & proxy)
  {
    init_if_necessary(viennacl::traits::context(proxy));

    viennacl::linalg::min_impl(proxy.lhs(), *this);
    return *this;
  }

  template<typename T1, typename T2>
  self_type & operator= (scalar_expression<T1, T2, op_sum> const & proxy)
  {
    init_if_necessary(viennacl::traits::context(proxy));

    viennacl::linalg::sum_impl(proxy.lhs(), *this);
    return *this;
  }


  template<typename T1, typename T2>
  self_type & operator= (scalar_expression<T1, T2, op_norm_frobenius> const & proxy)
  {
    init_if_necessary(viennacl::traits::context(proxy));

    viennacl::linalg::norm_frobenius_impl(proxy.lhs(), *this);
    return *this;
  }

  template<typename T1, typename T2>
  self_type & operator= (scalar_expression<T1, T2, op_flip_sign> const & proxy)
  {
    init_if_necessary(viennacl::traits::context(proxy));

    viennacl::linalg::as(*this, proxy.lhs(), NumericT(-1.0), 1, false, true);
    return *this;
  }


  self_type & operator += (scalar<NumericT> const & other)
  {
    assert( val_.get_active_handle_id() != viennacl::MEMORY_NOT_INITIALIZED && bool("Scalar not initialized!"));

    viennacl::linalg::asbs(*this,                                       // s1 =
                           *this, NumericT(1.0), 1, false, false,     //       s1 * 1.0
                           other, NumericT(1.0), 1, false, false);    //     + s2 * 1.0
    return *this;
  }
  self_type & operator += (NumericT other)
  {
    assert( val_.get_active_handle_id() != viennacl::MEMORY_NOT_INITIALIZED && bool("Scalar not initialized!"));

    viennacl::linalg::asbs(*this,                                       // s1 =
                           *this, NumericT(1.0), 1, false, false,     //       s1 * 1.0
                           other, NumericT(1.0), 1, false, false);    //     + s2 * 1.0
    return *this;
  }


  self_type & operator -= (scalar<NumericT> const & other)
  {
    assert( val_.get_active_handle_id() != viennacl::MEMORY_NOT_INITIALIZED && bool("Scalar not initialized!"));

    viennacl::linalg::asbs(*this,                                       // s1 =
                           *this, NumericT(1.0), 1, false, false,     //       s1 * 1.0
                           other, NumericT(-1.0), 1, false, false);   //     + s2 * (-1.0)
    return *this;
  }
  self_type & operator -= (NumericT other)
  {
    assert( val_.get_active_handle_id() != viennacl::MEMORY_NOT_INITIALIZED && bool("Scalar not initialized!"));

    viennacl::linalg::asbs(*this,                                       // s1 =
                           *this, NumericT(1.0), 1, false, false,     //       s1 * 1.0
                           other, NumericT(-1.0), 1, false, false);   //     + s2 * (-1.0)
    return *this;
  }


  self_type & operator *= (scalar<NumericT> const & other)
  {
    assert( val_.get_active_handle_id() != viennacl::MEMORY_NOT_INITIALIZED && bool("Scalar not initialized!"));

    viennacl::linalg::as(*this,                                       // s1 =
                         *this, other, 1, false, false);              //      s1 * s2
    return *this;
  }
  self_type & operator *= (NumericT other)
  {
    assert( val_.get_active_handle_id() != viennacl::MEMORY_NOT_INITIALIZED && bool("Scalar not initialized!"));

    viennacl::linalg::as(*this,                                       // s1 =
                         *this, other, 1, false, false);              //      s1 * s2
    return *this;
  }



  self_type & operator /= (scalar<NumericT> const & other)
  {
    assert( val_.get_active_handle_id() != viennacl::MEMORY_NOT_INITIALIZED && bool("Scalar not initialized!"));

    viennacl::linalg::as(*this,                                       // s1 =
                         *this, other, 1, true, false);              //      s1 / s2
    return *this;
  }
  self_type & operator /= (NumericT other)
  {
    assert( val_.get_active_handle_id() != viennacl::MEMORY_NOT_INITIALIZED && bool("Scalar not initialized!"));

    viennacl::linalg::as(*this,                                       // s1 =
                         *this, other, 1, true, false);              //      s1 / s2
    return *this;
  }



  self_type operator + (scalar<NumericT> const & other)
  {
    assert( val_.get_active_handle_id() != viennacl::MEMORY_NOT_INITIALIZED && bool("Scalar not initialized!"));

    self_type result = 0;

    viennacl::linalg::asbs(result,                                       // result =
                           *this, NumericT(1.0), 1, false, false,      //            *this * 1.0
                           other, NumericT(1.0), 1, false, false);     //          + other * 1.0

    return result;
  }
  template<typename T1, typename T2, typename OP>
  self_type operator + (scalar_expression<T1, T2, OP> const & proxy) const
  {
    assert( val_.get_active_handle_id() != viennacl::MEMORY_NOT_INITIALIZED && bool("Scalar not initialized!"));

    self_type result = proxy;

    viennacl::linalg::asbs(result,                                       // result =
                           *this, NumericT(1.0), 1, false, false,      //            *this * 1.0
                           result, NumericT(1.0), 1, false, false);     //        + result * 1.0

    return result;
  }
  self_type operator + (NumericT other)
  {
    assert( val_.get_active_handle_id() != viennacl::MEMORY_NOT_INITIALIZED && bool("Scalar not initialized!"));

    self_type result = 0;

    viennacl::linalg::asbs(result,                                       // result =
                           *this, NumericT(1.0), 1, false, false,      //            *this * 1.0
                           other, NumericT(1.0), 1, false, false);     //          + other * 1.0

    return result;
  }



  scalar_expression<const self_type, const self_type, op_flip_sign> operator-() const
  {
    return scalar_expression<const self_type, const self_type, op_flip_sign>(*this, *this);
  }


  self_type operator - (scalar<NumericT> const & other) const
  {
    assert( val_.get_active_handle_id() != viennacl::MEMORY_NOT_INITIALIZED && bool("Scalar not initialized!"));

    self_type result = 0;

    viennacl::linalg::asbs(result,                                       // result =
                           *this, NumericT(1.0), 1, false, false,      //            *this * 1.0
                           other, NumericT(-1.0), 1, false, false);    //          + other * (-1.0)

    return result;
  }
  template<typename T1, typename T2, typename OP>
  self_type operator - (scalar_expression<T1, T2, OP> const & proxy) const
  {
    assert( val_.get_active_handle_id() != viennacl::MEMORY_NOT_INITIALIZED && bool("Scalar not initialized!"));

    self_type result = proxy;

    viennacl::linalg::asbs(result,                                       // result =
                           *this, NumericT(1.0), 1 , false, false,    //            *this * 1.0
                           result, NumericT(-1.0), 1, false, false);  //          + result * (-1.0)

    return result;
  }
  scalar<NumericT> operator - (NumericT other) const
  {
    assert( val_.get_active_handle_id() != viennacl::MEMORY_NOT_INITIALIZED && bool("Scalar not initialized!"));

    self_type result = 0;

    viennacl::linalg::asbs(result,                                       // result =
                           *this, NumericT(1.0), 1, false, false,      //            *this * 1.0
                           other, NumericT(-1.0), 1, false, false);    //          + other * (-1.0)

    return result;
  }


  self_type operator * (scalar<NumericT> const & other) const
  {
    assert( val_.get_active_handle_id() != viennacl::MEMORY_NOT_INITIALIZED && bool("Scalar not initialized!"));

    scalar<NumericT> result = 0;

    viennacl::linalg::as(result,                                     // result =
                         *this, other, 1, false, false);              //          *this * other

    return result;
  }
  template<typename T1, typename T2, typename OP>
  self_type operator * (scalar_expression<T1, T2, OP> const & proxy) const
  {
    assert( val_.get_active_handle_id() != viennacl::MEMORY_NOT_INITIALIZED && bool("Scalar not initialized!"));

    self_type result = proxy;

    viennacl::linalg::as(result,                                       // result =
                         *this, result, 1, false, false);              //            *this * proxy

    return result;
  }
  self_type operator * (NumericT other) const
  {
    assert( val_.get_active_handle_id() != viennacl::MEMORY_NOT_INITIALIZED && bool("Scalar not initialized!"));

    scalar<NumericT> result = 0;

    viennacl::linalg::as(result,                                     // result =
                         *this, other, 1, false, false);              //          *this * other

    return result;
  }


  self_type operator / (scalar<NumericT> const & other) const
  {
    assert( val_.get_active_handle_id() != viennacl::MEMORY_NOT_INITIALIZED && bool("Scalar not initialized!"));

    self_type result = 0;

    viennacl::linalg::as(result,                                     // result =
                         *this, other, 1, true, false);              //           *this / other

    return result;
  }
  template<typename T1, typename T2, typename OP>
  self_type operator / (scalar_expression<T1, T2, OP> const & proxy) const
  {
    assert( val_.get_active_handle_id() != viennacl::MEMORY_NOT_INITIALIZED && bool("Scalar not initialized!"));

    self_type result = proxy;

    viennacl::linalg::as(result,                                     // result =
                         *this, result, 1, true, false);              //          *this / proxy

    return result;
  }
  self_type operator / (NumericT other) const
  {
    assert( val_.get_active_handle_id() != viennacl::MEMORY_NOT_INITIALIZED && bool("Scalar not initialized!"));

    self_type result = 0;

    viennacl::linalg::as(result,                                     // result =
                         *this, other, 1, true, false);              //            *this / other

    return result;
  }

  handle_type & handle() { return val_; }

  const handle_type & handle() const { return val_; }

private:

  void init_if_necessary(viennacl::context ctx)
  {
    if (val_.get_active_handle_id() == viennacl::MEMORY_NOT_INITIALIZED)
    {
      viennacl::backend::memory_create(val_, sizeof(NumericT), ctx);
    }
  }

  handle_type val_;
};


//stream operators:
template<class NumericT>
std::ostream & operator<<(std::ostream & s, const scalar<NumericT> & val)
{
  NumericT temp = val;
  s << temp;
  return s;
}

template<class NumericT>
std::istream & operator>>(std::istream & s, const scalar<NumericT> & val)
{
  NumericT temp;
  s >> temp;
  val = temp;
  return s;
}

} //namespace viennacl

#endif