librheolef-dev 6.5-1+b1 / usr / include / rheolef / vec_expr.h

#ifndef _RHEOLEF_VEC_EXPR_H
#define _RHEOLEF_VEC_EXPR_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
/// 
/// =========================================================================
//
// expressions of vectors without temporaries, based on boost::proto
// and that uses an iterator for the assignment
//
// author: Pierre.Saramito@imag.fr
//
// date: 13 january 2011
//
#include "rheolef/vec.h"
#include "rheolef/dis_inner_product.h"
#include "rheolef/dis_accumulate.h"
#include "rheolef/pretty_name.h"

#include <boost/mpl/bool.hpp>
#include <boost/proto/core.hpp>
#include <boost/proto/debug.hpp>
#include <boost/proto/context.hpp>
#include <boost/proto/transform.hpp>
#include <boost/proto/operators.hpp>
#include <boost/utility/enable_if.hpp>
#include <boost/typeof/std/vector.hpp>
#include <boost/typeof/std/complex.hpp>
#include <boost/type_traits/remove_reference.hpp>

namespace rheolef {
namespace mpl   = boost::mpl;
namespace proto = boost::proto;
using proto::_;

template <class Expr>
struct vec_expr;

// -----------------------------------------------------------
// iterator begin
// -----------------------------------------------------------
template<class Iterator>
struct vec_iterator_wrapper {
    typedef Iterator iterator;
    explicit vec_iterator_wrapper (Iterator iter1) : iter(iter1) {} 
    mutable Iterator iter;
};
struct vec_begin : proto::callable {
    template<class Sig>
    struct result;

    template<class This, class Container>
    struct result<This(Container)> : proto::result_of::as_expr
       <
        vec_iterator_wrapper<
	  typename boost::remove_reference<Container>::type::const_iterator>
       > {};

    template<class Container>
    typename result<vec_begin(const Container&)>::type
    operator ()(const Container& cont) const {
        vec_iterator_wrapper<typename Container::const_iterator> iter (cont.begin());
        return proto::as_expr (iter);
    }
};
// Here is the grammar that replaces vector terminals with their begin iterators
struct vec_grammar_begin
  : proto::or_<
        proto::when <proto::terminal<vec<_, _> >, vec_begin(proto::_value)>
      , proto::when <proto::terminal<_> >
      , proto::when <proto::nary_expr<_, proto::vararg<vec_grammar_begin> > >
    >
{};

// Here is an evaluation context that dereferences iterator terminals.
struct vec_dereference_context {
    // Unless this is an iterator terminal, use the default evaluation context
    template<class Expr, class EnableIf = void>
    struct eval : proto::default_eval<Expr, const vec_dereference_context> {};

    // Dereference iterator terminals.
    template<class Expr>
    struct eval<Expr,
       typename boost::enable_if<
         proto::matches<Expr, proto::terminal<vec_iterator_wrapper<_> > >
        >::type
       > {
        typedef typename proto::result_of::value<Expr>::type IteratorWrapper;
        typedef typename IteratorWrapper::iterator iterator;
        typedef typename std::iterator_traits<iterator>::reference result_type;

        result_type operator() (Expr &expr, const vec_dereference_context&) const {
            return *proto::value(expr).iter;
        }
    };
};
// Here is an evaluation context that increments iterator
// terminals.
struct vec_increment_context {
    // Unless this is an iterator terminal, use the default evaluation context
    template<class Expr, class EnableIf = void>
    struct eval : proto::null_eval<Expr, const vec_increment_context> {};

    // advance iterator terminals.
    template<class Expr> struct eval<Expr,
	typename boost::enable_if<
            proto::matches<Expr, proto::terminal<vec_iterator_wrapper<_> > >
        >::type
      > {
        typedef void result_type;

        result_type operator() (Expr &expr, const vec_increment_context&) const {
            ++proto::value(expr).iter;
        }
    };
};
// -----------------------------------------------------------
// Here is an evaluation context that indexes into a vec
// expression and combines the result.
// -----------------------------------------------------------
struct vec_subscript_context {
    typedef vec<Float>::size_type size_type;
    vec_subscript_context (size_type i) : _i(i) {}
    // Unless this is a vector terminal, use the default evaluation context
    template<class Expr, class EnableIf = void>
    struct eval : proto::default_eval<Expr, const vec_subscript_context> {};
    // Index vector terminals with our subscript.
    template<class Expr>
    struct eval<Expr, typename boost::enable_if<
              proto::matches<Expr, proto::terminal<vec<_, _> > > >::type> {
         typedef typename proto::result_of::value<Expr>::type::value_type result_type;
         result_type operator() (Expr &expr, const vec_subscript_context& ctx) const {
             return proto::value(expr)[ctx._i];
         }
    };
    // index iterator terminals
    template<class Expr>
    struct eval<Expr, typename boost::enable_if<
              proto::matches<Expr, proto::terminal<vec_iterator_wrapper<_> > > >::type > {
        typedef typename proto::result_of::value<Expr>::type IteratorWrapper;
        typedef typename IteratorWrapper::iterator iterator;
        typedef typename std::iterator_traits<iterator>::reference result_type;
        result_type operator() (Expr &expr, const vec_subscript_context& ctx) const {
            return *proto::value(expr).iter[ctx._i];
        }
    };
    mutable size_type _i;
};
// -----------------------------------------------------------
// get size
// -----------------------------------------------------------
struct vec_get_size_context {
    typedef vec<Float>::size_type size_type;
// allocator:
    vec_get_size_context () {}

// eval:
    template<class Expr, class EnableIf = void>
    struct eval : proto::null_eval<Expr, const vec_get_size_context> {};

    template<class Expr>
    struct eval<Expr, typename boost::enable_if<
            proto::matches<Expr, proto::terminal<vec<_,_> > > >::type> {
        typedef void result_type;
        result_type operator() (Expr &expr, const vec_get_size_context& ctx) const {
            ctx._ownership = proto::value(expr).ownership();
        }
    };
// accessors:
    size_type        size() const { return _ownership.size(); };
    distributor ownership() const { return _ownership; };
// data:
    mutable distributor _ownership;
};
// -----------------------------------------------------------
// check that sizes match
// -----------------------------------------------------------
struct vec_check_size_context {
    typedef vec<Float>::size_type size_type;
    vec_check_size_context (size_type size) : _size(size) {}

    template<class Expr, class EnableIf = void>
    struct eval : proto::null_eval<Expr, const vec_check_size_context> {};

    template<class Expr>
    struct eval<Expr, typename boost::enable_if<
            proto::matches<Expr, proto::terminal<vec<_,_> > > >::type> {
        typedef void result_type;
        result_type operator() (Expr &expr, const vec_check_size_context& ctx) const {
            if (ctx._size != proto::value(expr).size()) {
                error_macro ("incompatible sizes "<<ctx._size<< " and " 
		  << proto::value(expr).size() << " in vec<T> expression "
                  <<typename_macro(Expr));
            }
        }
    };
    size_type _size;
};
// -----------------------------------------------------------
// A grammar which matches all the assignment operators,
// so we can easily disable them.
// -----------------------------------------------------------
struct vec_assign_operators : proto::switch_<struct vec_assign_operators_cases> {};

// Here are the cases used by the switch_ above.
struct vec_assign_operators_cases {
    template<class Tag, int D = 0> struct case_  : proto::not_<_> {};

    template<int D> struct case_< proto::tag::plus_assign, D >         : _ {};
    template<int D> struct case_< proto::tag::minus_assign, D >        : _ {};
    template<int D> struct case_< proto::tag::multiplies_assign, D >   : _ {};
    template<int D> struct case_< proto::tag::divides_assign, D >      : _ {};
    template<int D> struct case_< proto::tag::modulus_assign, D >      : _ {};
    template<int D> struct case_< proto::tag::shift_left_assign, D >   : _ {};
    template<int D> struct case_< proto::tag::shift_right_assign, D >  : _ {};
    template<int D> struct case_< proto::tag::bitwise_and_assign, D >  : _ {};
    template<int D> struct case_< proto::tag::bitwise_or_assign, D >   : _ {};
    template<int D> struct case_< proto::tag::bitwise_xor_assign, D >  : _ {};
};
// An expression conforms to the vec_grammar if it is a terminal or some
// op that is not an assignment op. (Assignment will be handled specially.)
struct vec_grammar
  : proto::or_<
        proto::terminal<_>
      , proto::and_<
            proto::nary_expr<_, proto::vararg<vec_grammar> >
          , proto::not_<vec_assign_operators>
        >
    >
{};

// Expressions in the vec_domain will be wrapped in vec_expr<>
// and must conform to the vec_grammar
struct vec_domain : proto::domain<proto::generator<vec_expr>, vec_grammar> {};

// -------------------------------------------------------------------
// Here is vec_expr, a wrapper for expression types in the vec_domain.
// It mimics the array<T> and vec<T> interface
// -------------------------------------------------------------------
template<class Expr>
struct vec_expr : proto::extends<Expr, vec_expr<Expr>, vec_domain> {

// typedefs :

    typedef proto::extends<Expr, vec_expr<Expr>, vec_domain> base_type;

    typedef vec_get_size_context::size_type size_type;
    typedef typename boost::result_of<vec_grammar_begin(const Expr&)>::type raw_iterator;

// allocators:

    explicit vec_expr (const Expr& expr)
      : base_type(expr)
    {}

// accessors:

    distributor ownership() const {
        const vec_get_size_context get_size;
        proto::eval (*this, get_size);
        return get_size.ownership();
    }
    size_type size() const {
        return ownership().size();
    }
    size_type dis_size() const {
        return ownership().dis_size();
    }
    struct const_iterator_begin : raw_iterator {

       typedef std::input_iterator_tag        iterator_category;
       typedef typename proto::result_of::eval<raw_iterator,vec_dereference_context>::type value_type;
       typedef value_type reference;
       typedef value_type* pointer;
       typedef std::ptrdiff_t difference_type;

#if BOOST_VERSION < 104601
       const_iterator_begin (const vec_expr<Expr>& expr)
	 : raw_iterator(vec_grammar_begin() (expr))
       {}
#else // BOOST_VERSION < 104601
       const_iterator_begin (const vec_expr<Expr>& expr)
	 : raw_iterator(vec_grammar_begin() (expr.proto_expr_))
	   // ICI: il faut recuperer base_type(expr) de type Expr au lieu vec_expr<Expr>
       {}
#endif // BOOST_VERSION < 104601
       const_iterator_begin& operator++ () { 	
         const vec_increment_context inc = {};
         proto::eval (*this, inc);
	 return *this;
       }
       const_iterator_begin operator++ (int) { 
         const_iterator_begin tmp = *this;
	 operator++();
	 return tmp;
       }
       reference dereference () const {
         const vec_dereference_context deref = {};
         typedef typename proto::result_of::eval<raw_iterator,vec_dereference_context>::type my_ref;
         reference value = proto::eval(*this, deref);
         return value;
       }
       reference operator* () const { return dereference(); }
       reference operator* ()       { return dereference(); }
    };
    typedef const_iterator_begin const_iterator;
    const_iterator begin() const { return const_iterator(*this); }

// random access:

    typedef typename proto::result_of::eval<const Expr, const vec_subscript_context>::type value_type;
    typedef value_type scalar_type;
    value_type operator[] (size_type i) const {
        const vec_subscript_context get_subscript(i);
        return proto::eval(*this, get_subscript);
    }
private:
    // hide this:
    using proto::extends<Expr, vec_expr<Expr>, vec_domain>::operator[];
};
// Define a trait type for detecting vec terminals
template<class T>
struct is_vec : mpl::false_ {};

template<class T, class M>
struct is_vec<vec<T, M> > : mpl::true_ {};

namespace vec_detail {
    template<class ForwardIterator, class Expr, class Op>
    void evaluate (ForwardIterator begin, ForwardIterator end, const Expr& expr, Op op) {
        const vec_increment_context   inc   = {};
        const vec_dereference_context deref = {};
        typename boost::result_of<vec_grammar_begin(const Expr&)>::type expr_begin
	  = vec_grammar_begin() (expr);
#ifndef TRY_VEC_EXPR
        for (; begin != end; ++begin) {
            op (*begin, proto::eval(expr_begin, deref));
            proto::eval (expr_begin, inc);
        }
#endif // TRY_VEC_EXPR
    }
    struct assign_op {
        template<class T, class U>
        void operator() (T &t, U const &u) const { t = u; }
    };
} // namespace vec_detail

// -------------------------------------------
// x = expr; x += expr; x -= expr;
// -------------------------------------------
template<class T, class M>
template<class Expr>
inline
vec<T,M>& 
vec<T,M>::operator= (const vec_expr<Expr>& expr) {
    // get and check sizes:
    const vec_get_size_context get_size;
    proto::eval (proto::as_expr<vec_domain>(expr), get_size);
    size_type   expr_size      = get_size.size();
    if (array<T,M>::dis_size() == 0) {
        distributor expr_ownership = get_size.ownership();
	resize (expr_ownership);
    }
    const vec_check_size_context check_size (array<T,M>::size());
    proto::eval (proto::as_expr<vec_domain>(expr), check_size); // error if the sizes don't match
    // perform all computations here:
    vec_detail::evaluate (array<T,M>::begin(), array<T,M>::end(), proto::as_expr<vec_domain>(expr), vec_detail::assign_op());
    return *this;
}
template<class T, class M>
template<class Expr>
inline
vec<T,M>::vec (const vec_expr<Expr>& expr)
  : array<T,M>()
{
    operator= (expr);
}
template<class T, class M>
inline
T
norm2 (const vec<T,M>& x)
{
  return dot(x,x);
}
template<class T, class M>
inline
T
norm (const vec<T,M>& x)
{
  return sqrt(norm2(x));
}

} // namespace rheolef
#endif // _RHEOLEF_VEC_EXPR_H