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

#ifndef _RHEO_VEC_H
#define _RHEO_VEC_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
/// 
/// =========================================================================

# include "rheolef/disarray.h"
# include "boost/numeric/ublas/storage.hpp"

namespace rheolef {

// for vec[range]:
using boost::numeric::ublas::basic_range;
using boost::numeric::ublas::range;
template <class T, class M> class vec_range;
template <class T, class M> class vec_range_const;

namespace details {
template <class Expr> struct is_vec;
template <class Expr> struct is_vec_expr_v2_arg;
} // namespace details

// for vec = {x,y};
template <class T, class M> class vec_concat_value;

/*Class:vec
NAME: @code{vec} - vector in distributed environment (@PACKAGE@-@VERSION@)
SYNOPSYS:       
  STL-like vector container for a sequential or
  distributed memory machine model.
  Additional operation fom classical algebra.
EXAMPLE:
   A sample usage of the class is:
   @example
     int main(int argc, char**argv) @{
        environment distributed(argc, argv);
        vec<double> x(100, 3.14);
        dout << x << endl;
     @}
   @end example
IMPLEMENTATION NOTE:
   Implementation use disarray<T,M>.
SEE ALSO: "disarray"(1)
AUTHORS: Pierre.Saramito@imag.fr
DATE:   19 november 1998
End:
*/

//<vec:
template <class T, class M = rheo_default_memory_model>
class vec : public disarray<T, M> {
public:

// typedef:

    typedef disarray<T, M> base;
    typedef T                                                value_type;
    typedef typename base::size_type                         size_type;
    typedef std::ptrdiff_t                                   difference_type;
    typedef basic_range<size_type, difference_type>          range_type;
    typedef typename base::reference                         reference;
    typedef typename base::const_reference                   const_reference;
    typedef typename base::iterator                          iterator;
    typedef typename base::const_iterator                    const_iterator;
    typedef typename float_traits <value_type>::type         float_type;

// allocator/deallocator:

    vec (const vec<T,M>&);
    vec<T,M>& operator= (const vec<T,M>& x);

    vec (const distributor& ownership,
	const T&  init_val = std::numeric_limits<T>::max());

    vec(size_type dis_size = 0,
	const T&  init_val = std::numeric_limits<T>::max());

    void resize (
        const distributor& ownership,
        const T&  init_val = std::numeric_limits<T>::max());

    void resize (
        size_type size = 0,
        const T&  init_val = std::numeric_limits<T>::max());

// accessors:

    const_reference operator[] (size_type i) const;
    reference       operator[] (size_type i);

    T min () const;
    T max () const;
    T max_abs () const;

// range:

    vec(const vec_range<T,M>& vr);
    vec(const vec_range_const<T,M>& vr);
    vec<T,M>& operator= (const vec_range<T,M>& vr);
    vec<T,M>& operator= (const vec_range_const<T,M>& vr);

    vec_range_const<T,M> operator[] (const range_type& r) const;
    vec_range<T,M>       operator[] (const range_type& r);

// assignment to a constant:

    vec<T,M>& operator= (const int& expr);
    vec<T,M>& operator= (const T& expr);

// expression template:

    template <class Expr, 
              class Sfinae
	          = typename std::enable_if<
			     details::is_vec_expr_v2_arg<Expr>::value
			&& ! details::is_vec<Expr>::value
		    >::type>
    vec (const Expr& expr);

    template <class Expr,
              class Sfinae 
	          = typename std::enable_if<
			     details::is_vec_expr_v2_arg<Expr>::value
			&& ! details::is_vec<Expr>::value
		    >::type>
    vec<T, M>& operator=  (const Expr& expr);

// initializer list (c++ 2011):

#ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST
    vec (const std::initializer_list<vec_concat_value<T,M> >& init_list);
    vec<T,M>& operator= (const std::initializer_list<vec_concat_value<T,M> >& init_list);
#endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST
};
//>vec:

// ----------------------------------------------------------------------------
// inlined
// ----------------------------------------------------------------------------
template <class T, class M>
inline
vec<T,M>::vec (const vec<T,M>& x)
  : disarray<T,M>(x)
{
}
template <class T, class M>
inline
vec<T,M>&
vec<T,M>::operator= (const vec<T,M>& x)
{
  disarray<T,M>::operator= (x);
  return *this;
}
template <class T, class M>
inline
vec<T,M>::vec (
    const distributor& ownership,
    const T&  init_val)
  : disarray<T,M>(ownership,init_val)
{
}
template <class T, class M>
inline
vec<T,M>::vec (
    size_type dis_size,
    const T&  init_val)
  : disarray<T,M>(dis_size,init_val)
{
}
template <class T, class M>
inline
void
vec<T,M>::resize (
        const distributor& ownership,
        const T&  init_val)
{
    base::resize (ownership, init_val);
}
template <class T, class M>
inline
void
vec<T,M>::resize (
        size_type dis_size,
        const T&  init_val)
{
    base::resize (dis_size, init_val);
}
// TODO: group cstors vec(int) & vec(T) via Sfinae
template <class T, class M>
inline
vec<T,M>&
vec<T,M>::operator= (const int& expr)
{
    std::fill (disarray<T,M>::begin(), disarray<T,M>::end(), expr);
    return *this;
}
template <class T, class M>
inline
vec<T,M>&
vec<T,M>::operator= (const T& expr)
{
    std::fill (disarray<T,M>::begin(), disarray<T,M>::end(), expr);
    return *this;
}
template <class T, class M>
inline
vec<T,M>&
vec<T,M>::operator= (const vec_range_const<T,M>& vr)
{
    distributor ownership (distributor::decide, vr._u.comm(), vr._r.size());
    resize (ownership);
    std::copy (vr.begin(), vr.end(), base::begin());
    return *this;
}
template <class T, class M>
inline
vec<T,M>&
vec<T,M>::operator= (const vec_range<T,M>& vr)
{
    operator= (vec_range_const<T,M>(vr));
    return *this;
}
template <class T, class M>
inline
vec<T,M>::vec(const vec_range<T,M>& vr)
  : disarray<T,M>()
{
    operator= (vr);
}
template <class T, class M>
inline
vec<T,M>::vec(const vec_range_const<T,M>& vr)
  : disarray<T,M>()
{
    operator= (vr);
}
template <class T, class M>
inline
typename vec<T,M>::const_reference
vec<T,M>::operator[] (size_type i) const
{
    return base::operator[] (i);
}
template <class T, class M>
inline
typename vec<T,M>::reference
vec<T,M>::operator[] (size_type i)
{
    return base::operator[] (i);
}
template <class T, class M>
inline
vec_range<T,M>
vec<T,M>::operator[] (const range_type& r)
{
    return vec_range<T,M> (*this, r);
}
template <class T, class M>
inline
vec_range_const<T,M>
vec<T,M>::operator[] (const range_type& r) const
{
    return vec_range_const<T,M> (*this, r);
}
template <class T, class M>
T
vec<T,M>::min () const
{
    T val = std::numeric_limits<T>::max();
    for (const_iterator iter = base::begin(), last = base::end(); iter != last; iter++) {
      val = std::min(val, *iter);
    }
#ifdef _RHEOLEF_HAVE_MPI
    val = mpi::all_reduce (base::comm(), val, mpi::minimum<T>());
#endif // _RHEOLEF_HAVE_MPI
    return val;
}
template <class T, class M>
T
vec<T,M>::max () const
{
    T val = std::numeric_limits<T>::min();
    for (const_iterator iter = base::begin(), last = base::end(); iter != last; iter++) {
      val = std::max(val, *iter);
    }
#ifdef _RHEOLEF_HAVE_MPI
    val = mpi::all_reduce (base::comm(), val, mpi::maximum<T>());
#endif // _RHEOLEF_HAVE_MPI
    return val;
}
template <class T, class M>
T
vec<T,M>::max_abs () const
{
    T val = 0;
    for (const_iterator iter = base::begin(), last = base::end(); iter != last; iter++) {
      val = std::max(val, abs(*iter));
    }
#ifdef _RHEOLEF_HAVE_MPI
    val = mpi::all_reduce (base::comm(), val, mpi::maximum<T>());
#endif // _RHEOLEF_HAVE_MPI
    return val;
}
template <class T>
inline
idiststream&
operator >> (idiststream& ips,  vec<T,sequential>& x)
{ 
    return x.get_values(ips);
}
template <class T, class M> 
inline
odiststream&
operator << (odiststream& ods, const vec<T,M>& x)
{
    iorheo::flag_type format = iorheo::flags(ods.os()) & iorheo::format_field;
    if (format [iorheo::matlab] || format [iorheo::sparse_matlab]) {
      return x.data().put_matlab (ods);
    }
    // default is raw output
    return x.put_values(ods);
}
#ifdef _RHEOLEF_HAVE_MPI
template <class T>
inline
idiststream&
operator >> (idiststream& ips,  vec<T,distributed>& x)
{ 
    return x.get_values(ips); 
}
#ifdef TO_CLEAN
template <class T> 
inline
odiststream&
operator << (odiststream& ods, const vec<T,distributed>& x)
{
    iorheo::flag_type format = iorheo::flags(ods.os()) & iorheo::format_field;
    if (format [iorheo::matlab] || format [iorheo::sparse_matlab]) {
      return x.put_matlab (ods);
    }
    // default is raw output
    return x.put_values(ods);
}
#endif // TO_CLEAN
#endif // _RHEOLEF_HAVE_MPI
// -------------------------------------------
// norm(x) ; dot(x,y)
// -------------------------------------------
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 // _RHEO_VEC_H