librheolef-dev 6.7-6 / usr / include / rheolef / vec_expr_v2.h

#ifndef _RHEOLEF_VEC_EXPR_v2_H
#define _RHEOLEF_VEC_EXPR_v2_H
///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/// 
/// =========================================================================
// vec-valued affine expressions:
//
//      vec w = 2*u - v + 1:
//
// author: Pierre.Saramito@imag.fr
//
// date: 14 september 2015
//
// Notes; implementation uses template expressions and SFINAE techiques
// * template expressions allows to eliminate temporaries
// * SFINAE reduces the code size
//
// OVERVIEW:
// 1. linear/affine algebra operators
// 1.1. vec_expr, a type concept for expression types
// 1.2. unary operations
// 1.3. binary operations
// 2. completing the interface of the vec<T,M> class
// 2.1. cstor & assignment
// 2.2. computed assignment
// 2.3. scalar product
//
#include "rheolef/vec.h"
#include "rheolef/dis_inner_product.h"
#include "rheolef/dis_accumulate.h"
#include "rheolef/expr_utilities.h"

namespace rheolef {

// ===================================================================
// 1. linear/affine algebra operators
// ===================================================================
// -------------------------------------------------------------------
// 1.1. vec_expr, a type concept for expression types
// -------------------------------------------------------------------
namespace details {

// Define a trait type for detecting vec expression valid arguments
template<class T>          struct is_vec              : std::false_type {};
template<class T, class M> struct is_vec<vec<T, M> >  : std::true_type  {};

// Define a trait type for detecting vec expression valid arguments
template<class T>          struct is_vec_expr_v2_arg              : std::false_type {};
template<class T, class M> struct is_vec_expr_v2_arg <vec<T, M> > : std::true_type  {};

} // namespace details
// -------------------------------------------
// 1.2. unary operations
// -------------------------------------------
namespace details {

template <class Op, class Expr>
struct vec_expr_v2_unary {

// typedefs:

  typedef typename Expr::size_type      size_type;
  typedef typename Expr::value_type     value_type;
  typedef typename Expr::memory_type    memory_type;
  typedef typename Expr::const_iterator expr_const_iterator;
  typedef typename float_traits <value_type>::type   float_type;

// allocatos:

  vec_expr_v2_unary (const Op& op, const Expr& expr)
   : _op(op), _expr_iter(expr.begin()), _ownership(expr.ownership()) {}

  template <class BinaryOp, class Constant>
  vec_expr_v2_unary (const BinaryOp& binop, const Constant& c, const Expr& expr)
   : _op(binop,c), _expr_iter(expr.begin()), _ownership(expr.ownership()) {}

  template <class BinaryOp, class Constant>
  vec_expr_v2_unary (const BinaryOp& binop, const Expr& expr, const Constant& c)
   : _op(binop,c), _expr_iter(expr.begin()), _ownership(expr.ownership()) {}

// accessors:

  const distributor& ownership() const { return _ownership; }
  size_type size() const { return _ownership.size(); }

// minimal forward iterator interface:

  struct const_iterator {
    typedef std::forward_iterator_tag         iterator_category;
    typedef typename Expr::value_type         value_type;
    typedef value_type&                       reference;
    typedef value_type*                       pointer;
    typedef std::ptrdiff_t                    difference_type;
    const_iterator (Op op, expr_const_iterator expr_iter)
     : _op(op), _expr_iter (expr_iter) {}
    const_iterator& operator++ () { ++_expr_iter; return *this; }
    value_type operator* () const {  return _op (*_expr_iter); }
  protected:
    const Op                                   _op;
    expr_const_iterator                        _expr_iter;
  };
  const_iterator begin() const { return const_iterator (_op, _expr_iter); }
protected:
  const Op                                   _op;
  const expr_const_iterator                  _expr_iter;
  const distributor                          _ownership;
};
template<class Op, class Expr> struct is_vec_expr_v2_arg <vec_expr_v2_unary<Op,Expr> > : std::true_type {};

} // namespace details

#define _RHEOLEF_vec_expr_v2_unary_operator(OP, FUNCTOR) 			\
template <class Expr>								\
inline										\
typename 									\
std::enable_if<									\
  details::is_vec_expr_v2_arg<Expr>::value,					\
  details::vec_expr_v2_unary<							\
    FUNCTOR,									\
    Expr									\
  >										\
>::type										\
operator OP (const Expr& expr)							\
{										\
  typedef details::vec_expr_v2_unary <FUNCTOR, Expr>   expr_t;			\
  return expr_t (FUNCTOR(), expr); 						\
}

_RHEOLEF_vec_expr_v2_unary_operator(+, details::generic_unary_plus<>)
_RHEOLEF_vec_expr_v2_unary_operator(-, details::generic_negate<>)
#undef _RHEOLEF_vec_expr_v2_unary_operator

// -------------------------------------------
// 1.3. binary operations
// -------------------------------------------
namespace details {

template <class Op, class Expr1, class Expr2>
struct vec_expr_v2_binary {

  typedef typename Expr1::size_type      size_type;
  typedef typename promote<
    typename Expr1::value_type,
    typename Expr2::value_type>::type   value_type;
  typedef typename Expr1::memory_type   memory_type; // TODO: check Expr2::memory_type
  typedef typename Expr1::const_iterator expr1_const_iterator;
  typedef typename Expr2::const_iterator expr2_const_iterator;
  typedef typename float_traits <value_type>::type   float_type;

// allocators:

  vec_expr_v2_binary (const Op& op, const Expr1& expr1, const Expr2& expr2)
   : _op (op),
     _iter1 (expr1.begin()),
     _iter2 (expr2.begin()),
     _ownership (expr1.ownership())
  {
    check_macro (expr1.size() == expr2.size(), "linear binary vec expression: incompatible sizes "
                << expr1.size() << " and " << expr2.size());
  }

// accessors:

  const distributor& ownership() const { return _ownership; }
  size_type size() const { return _ownership.size(); }

// minimal forward iterator interface:

  struct const_iterator {
    typedef std::forward_iterator_tag         iterator_category;
    typedef typename promote<
      typename Expr1::value_type,
      typename Expr2::value_type>::type       value_type;
    typedef value_type&                       reference;
    typedef value_type*                       pointer;
    typedef std::ptrdiff_t                    difference_type;
    const_iterator (Op op, expr1_const_iterator iter1, expr2_const_iterator iter2)
     : _op(op), _iter1 (iter1), _iter2 (iter2) {}
    const_iterator& operator++ () { ++_iter1; ++_iter2; return *this; }
    value_type operator* () const {  return _op (*_iter1, *_iter2); }
  protected:
    const Op             _op;
    expr1_const_iterator _iter1;
    expr2_const_iterator _iter2;
  };
  const_iterator begin() const { return const_iterator (_op, _iter1, _iter2); }
protected:
  const Op               _op;
  expr1_const_iterator   _iter1;
  expr2_const_iterator   _iter2;
  const distributor      _ownership;
};
template<class Op, class Expr1, class Expr2> struct is_vec_expr_v2_arg <vec_expr_v2_binary<Op,Expr1,Expr2> > : std::true_type {};

template <class Op, class Expr1, class Expr2, class Sfinae = void>
struct vec_expr_v2_binary_traits { /* catch-all case */ };

// vec_expr +- vec_expr
template <class Op, class Expr1, class Expr2>
struct vec_expr_v2_binary_traits <Op, Expr1, Expr2,
  typename std::enable_if<
    details::is_vec_expr_v2_arg<Expr1>::value &&
    details::is_vec_expr_v2_arg<Expr2>::value>::type>
{
  typedef vec_expr_v2_binary <Op,Expr1,Expr2> type;
};
// constant +-* vec_expr
template <class Op, class Expr1, class Expr2>
struct vec_expr_v2_binary_traits <Op, Expr1, Expr2,
  typename std::enable_if<
    details::is_rheolef_arithmetic<Expr1>::value &&
    details::is_vec_expr_v2_arg<Expr2>::value>::type>
{
  typedef generic_binder1st <Op, Expr1>                           fun_t;
  typedef vec_expr_v2_unary <fun_t,Expr2>                         type;
};
// vec_expr +-* constant
template <class Op, class Expr1, class Expr2>
struct vec_expr_v2_binary_traits <Op, Expr1, Expr2,
  typename std::enable_if<
    details::is_vec_expr_v2_arg<Expr1>::value &&
    details::is_rheolef_arithmetic<Expr2>::value>::type>
{
  typedef generic_binder2nd <Op, Expr2>                           fun_t;
  typedef vec_expr_v2_unary <fun_t,Expr1>                         type;
};

} // namespace details

// x+y; x+c ; c+x
#define _RHEOLEF_vec_expr_v2_binary_operator(OP, FUNCTOR) 			\
template <class Expr1, class Expr2>						\
inline										\
typename 									\
std::enable_if<									\
    (details::is_vec_expr_v2_arg<Expr1>::value &&				\
     details::is_vec_expr_v2_arg<Expr2>::value) ||				\
    (details::is_rheolef_arithmetic<Expr1>::value &&				\
     details::is_vec_expr_v2_arg<Expr2>::value) ||				\
    (details::is_vec_expr_v2_arg<Expr1>::value &&				\
     details::is_rheolef_arithmetic<Expr2>::value),				\
  typename 									\
  details::vec_expr_v2_binary_traits<						\
    FUNCTOR,									\
    Expr1, Expr2								\
  >::type									\
>::type										\
operator OP (const Expr1& expr1, const Expr2& expr2)				\
{										\
  typedef typename details::vec_expr_v2_binary_traits <FUNCTOR, Expr1, Expr2>::type  expr_t; \
  return expr_t (FUNCTOR(), expr1, expr2);					\
}
_RHEOLEF_vec_expr_v2_binary_operator(+, details::generic_plus<>)
_RHEOLEF_vec_expr_v2_binary_operator(-, details::generic_minus<>)
#undef _RHEOLEF_vec_expr_v2_binary_operator

// c*x ; x*c
template <class Expr1, class Expr2>
inline
typename
std::enable_if<
    (details::is_rheolef_arithmetic<Expr1>::value &&
     details::is_vec_expr_v2_arg<Expr2>::value) ||
    (details::is_vec_expr_v2_arg<Expr1>::value &&
     details::is_rheolef_arithmetic<Expr2>::value),
  typename
  details::vec_expr_v2_binary_traits<
    details::generic_multiplies<>,
    Expr1, Expr2
  >::type
>::type
operator* (const Expr1& expr1, const Expr2& expr2)
{
  typedef details::generic_multiplies<> fun_t;
  typedef typename details::vec_expr_v2_binary_traits <fun_t, Expr1, Expr2>::type  expr_t;
  return expr_t (fun_t(), expr1, expr2);
}

// x/c
template <class Expr1, class Expr2>
inline
typename
std::enable_if<
    (details::is_vec_expr_v2_arg<Expr1>::value &&
     details::is_rheolef_arithmetic<Expr2>::value),
  typename
  details::vec_expr_v2_binary_traits<
    details::generic_divides<>,
    Expr1, Expr2
  >::type
>::type
operator/ (const Expr1& expr1, const Expr2& expr2)
{
  typedef details::generic_divides<> fun_t;
  typedef typename details::vec_expr_v2_binary_traits <fun_t, Expr1, Expr2>::type  expr_t;
  return expr_t (fun_t(), expr1, expr2);
}
// ===================================================================
// 2. completing the interface of the vec<T,M> class
// ===================================================================
// ---------------------------------------------------------------------------
// 2.1. cstor & assignment
// ---------------------------------------------------------------------------
// vec<double> x = expr;
template<class T, class M>
template<class Expr, class Sfinae>
inline
vec<T,M>::vec (const Expr& expr)
  : disarray<T,M>()
{
    operator= (expr);
}
// x = expr;
template< class T, class M>
template <class Expr, class Sfinae>
inline
vec<T, M>&
vec<T,M>::operator=  (const Expr& expr)
{
  if (disarray<T,M>::dis_size() == 0) {
    resize (expr.ownership());
  } else {
    std::size_t n = disarray<T,M>::size();
    check_macro (n == expr.size(),
      "vec = vec_expression : incompatible size "
         << n << " and " << expr.size());
  }
  details::assign_with_operator (disarray<T,M>::begin(), disarray<T,M>::end(), expr.begin(), details::assign_op());
  return *this;
}
// ---------------------------------------------------------------------------
// 2.2) computed assignment
// ---------------------------------------------------------------------------
// x +-= expr;
#define _RHEOLEF_vec_expr_v2_op_assign(OP, FUNCTOR) 				\
template<class T, class M, class Expr>						\
inline										\
typename std::enable_if<							\
  details::is_vec_expr_v2_arg<Expr>::value,					\
  vec<T,M>&>::type								\
operator OP (vec<T,M>& x, const Expr& expr)					\
{										\
  check_macro (x.size() == expr.size(), "vec " << #OP << " vec_expression : incompatible spaces " \
            << x.size() << " and " << expr.size());				\
  details::assign_with_operator (x.begin(), x.end(), expr.begin(), FUNCTOR()); \
  return x;									\
}
_RHEOLEF_vec_expr_v2_op_assign(+=, details::plus_assign)
_RHEOLEF_vec_expr_v2_op_assign(-=, details::minus_assign)
#undef _RHEOLEF_vec_expr_v2_op_assign

// x -+*/= c
#define _RHEOLEF_vec_expr_v2_op_assign_constant(OP, FUNCTOR)			\
template<class T, class M, class Expr>						\
inline										\
typename std::enable_if<							\
  details::is_rheolef_arithmetic<Expr>::value,						\
  vec<T,M>&>::type								\
operator OP (vec<T,M>& x, const Expr& expr)					\
{										\
  details::assign_with_operator (x.begin(), x.end(), details::iterator_on_constant<Expr>(expr), FUNCTOR()); \
  return x;									\
}
_RHEOLEF_vec_expr_v2_op_assign_constant(+=, details::plus_assign)
_RHEOLEF_vec_expr_v2_op_assign_constant(-=, details::minus_assign)
_RHEOLEF_vec_expr_v2_op_assign_constant(*=, details::multiplies_assign)
_RHEOLEF_vec_expr_v2_op_assign_constant(/=, details::divides_assign)
#undef _RHEOLEF_vec_expr_v2_op_assign_constant

// ---------------------------------------------------------------------------
// 2.3. scalar product
// ---------------------------------------------------------------------------
// dot (x,y)
template <class Expr1, class Expr2>
inline
typename
std::enable_if<
  details::is_vec_expr_v2_arg<Expr1>::value &&
  details::is_vec_expr_v2_arg<Expr2>::value,
  typename promote<
    typename Expr1::float_type,
    typename Expr2::float_type>::type
>::type
dot (const Expr1& expr1, const Expr2& expr2)
{
  typedef typename Expr1::memory_type M;
  return dis_inner_product (expr1.begin(), expr2.begin(), expr1.size(), expr1.ownership().comm(), M());
}
// dot (c,x)
template <class Expr1, class Expr2>
inline
typename
std::enable_if<
  details::is_rheolef_arithmetic<Expr1>::value &&
  details::is_vec_expr_v2_arg<Expr2>::value,
  typename Expr2::float_type
>::type
dot (const Expr1& expr1, const Expr2& expr2)
{
  typedef typename Expr2::memory_type M;
  return expr1*dis_accumulate (expr2.begin(), expr2.size(), expr2.ownership().comm(), M());
}
// dot (x,c) 
template <class Expr1, class Expr2>
inline
typename
std::enable_if<
  details::is_vec_expr_v2_arg<Expr1>::value &&
  details::is_rheolef_arithmetic<Expr2>::value,
  typename Expr1::float_type
>::type
dot (const Expr1& expr1, const Expr2& expr2)
{
  typedef typename Expr1::memory_type M;
  return dis_accumulate (expr1.begin(), expr1.size(), expr1.ownership().comm(), M())*expr2;
}

} // namespace rheolef
#endif // _RHEOLEF_VEC_EXPR_v2_H