#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
// 
// =========================================================================
// AUTHORS: Pierre.Saramito@imag.fr
// DATE:   19 november 1998

namespace rheolef {
/**
@linalgclassfile vec distributed vector

Description
===========
This vector class supports both the `sequential`
and the `distributed` memory model.
In addition, standard linear algebra is supported.

Example
=======

        vec<double> x (100, 3.14);
        vec<double> y (100, 6.28);
        vec<double> z = 2.5*x + y;
        dout << x << endl;

Implementation
==============
@showfromfile
The `vec` class is a template class
with both the floating type and the memory model as parameters.
The implementation bases on the @ref disarray_4 container.

@snippet vec.h verbatim_vec
*/
} // namespace rheolef


# include "rheolef/disarray.h"
# include "rheolef/range.h"

namespace rheolef {

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;
template <class T, class M> class vec_concat_value; // for vec = {x,y};
} // namespace details


// [verbatim_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 range                                            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 (const std::initializer_list<details::vec_concat_value<T,M> >& init_list);

    vec<T,M>& operator= (const std::initializer_list<details::vec_concat_value<T,M> >& init_list);
    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;

    int constraint_process_rank() 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);
};
// [verbatim_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)
// -------------------------------------------
//! @brief norm2(x): see the @ref expression_3 page for the full documentation
template<class T, class M>
inline
T
norm2 (const vec<T,M>& x)
{
  return dot(x,x);
}
//! @brief norm(x): see the @ref expression_3 page for the full documentation
template<class T, class M>
inline
T
norm (const vec<T,M>& x)
{
  return sqrt(norm2(x));
}

} // namespace rheolef
#endif // _RHEO_VEC_H