Proto_BoxData.H

#pragma once
#ifndef _PROTO_RECTMDARRAY_H_
#define _PROTO_RECTMDARRAY_H_
#include "Proto_Timer.H"

#include "Proto_MemType.H"

#include "Proto_Box.H"
#include "Proto_PowerItoI.H"
#include "Proto_Stencil.H"
#include "Proto_Reduction.H"

#include <memory>   //for shared_ptr
#include <iostream>
#include <fstream>  // for parallel pretty printing
#include <iomanip>  // for pretty printing
#include <limits>   // for max / min functions
#include <cstdlib>  // for size_t
#include <cmath>    // floating point abs
#include <map>
#include <unordered_map>
#include "implem/Proto_Stack.H"
#include "Proto_DisjointBoxLayout.H"

// Please do not set this value to false (slow, race condition).
#define DEFAULT_USE_STACK true

namespace Proto
{
   enum BoxDataOp
  {
    Invalid = -1,
    Copy    = 0,
    Add     = 1,
    Assign,
    Subtract, 
    Multiply, 
    Divide, 
    NUM_OP_TYPES
  };
}

 
// This is to facilitate automated testing.
#ifdef PROTO_MEM_CHECK
class memcheck
{
public:
  static int numcopies;
  inline  static void FLUSH_CPY(){numcopies = 0;}
};
#endif
using std::shared_ptr;
namespace Proto 
{

// forward declarations
  class H5Handler;
  //class DisjointBoxLayout;  //used for write(...) hack
  //class ProblemDomain;      //used for write(...) hack
  //template<typename T, unsigned int C, MemType, MEM, unsigned int CTR>
  //class LevelBoxData;       //used for write(...) hack
  
  template<typename T> class Stencil;
  template<typename T, unsigned int C, MemType MEMTYPE, unsigned char D, unsigned char E>
  class LazyStencil;
/*
    template <typename T, unsigned int C, unsigned char D, unsigned char E, MemType MEMTYPE>
    class BoxData;
*/

    template<MemType MEMTYPE>
  inline void null_deleter_boxdata(void* ptr)
  {
  }




class CInterval 
{
public:

    inline CInterval(unsigned int a_c0, unsigned int a_c1,
                     unsigned int a_d0 = 0, unsigned int a_d1 = 0,
                     unsigned int a_e0 = 0, unsigned int a_e1 = 0)
    {
      PROTO_ASSERT(a_c1 >= a_c0,
      "CInterval(...) invalid for a_c0 = %i and a_c1 = %i. \
      lower bound must be less than or equal to high bound.",
      a_c0, a_c1);
      PROTO_ASSERT(a_d1 >= a_d0,
      "CInterval(...) invalid for a_d0 = %i and a_d1 = %i. \
      lower bound must be less than or equal to high bound.",
      a_d0, a_d1);
      PROTO_ASSERT(a_e1 >= a_e0,
      "CInterval(...) invalid for a_e0 = %i and a_e1 = %i. \
      lower bound must be less than or equal to high bound.",
      a_e0, a_e1);
      m_comps[0][0] = a_c0; m_comps[0][1] = a_c1;
      m_comps[1][0] = a_d0; m_comps[1][1] = a_d1;
      m_comps[2][0] = a_e0; m_comps[2][1] = a_e1;
    }

    CInterval(std::initializer_list<
             std::initializer_list<unsigned int>> a_lst)
    {
      int comp = 0;
      for (auto& c : a_lst)
      {
        if (c.size() == 0)
        {
          m_comps[comp][0] = 0;
          m_comps[comp][1] = 0;
          comp++;
          continue;
        }
        int bound = 0;
        for (auto& b : c)
        {
          m_comps[comp][bound] = b;
          bound++;
        }
        comp++;
      }
      PROTO_ASSERT(m_comps[0][1] >= m_comps[0][0],
      "CInterval(...) invalid for a_c0 = %i and a_c1 = %i. \
      lower bound must be less than or equal to high bound.",
      m_comps[0][0], m_comps[0][1]);
      PROTO_ASSERT(m_comps[1][1] >= m_comps[1][0],
      "CInterval(...) invalid for a_d0 = %i and a_d1 = %i. \
      lower bound must be less than or equal to high bound.",
      m_comps[1][0], m_comps[1][1]);
      PROTO_ASSERT(m_comps[2][1] >= m_comps[2][0],
      "CInterval(...) invalid for a_e0 = %i and a_e1 = %i. \
      lower bound must be less than or equal to high bound.",
      m_comps[2][0], m_comps[2][1]);
    }
    

    inline unsigned int low(unsigned int a_comp) const
    {
      PROTO_ASSERT((a_comp < 3),
      "CInterval::low(unsigned int a_comp) invalid for a_comp = %i.\
      a_comp must be in [0,3).",
      a_comp);
      return m_comps[a_comp][0];
    };
    

    inline unsigned int high(unsigned int a_comp) const
    {
      PROTO_ASSERT((a_comp < 3),
      "CInterval::high(unsigned int a_comp) invalid for a_comp = %i.\
      a_comp must be in [0,3).",
      a_comp);
      return m_comps[a_comp][1];
    };
    

    inline bool contains(unsigned int a_index,
                         unsigned int a_comp) const
    {
        PROTO_ASSERT((a_comp < 3),
        "CInterval::contains(unsigned int a_comp) \
        invalid for a_comp = %i. a_comp must be in [0,3).",
        a_comp);
        return ((a_index >= m_comps[a_comp][0])
             && (a_index <= m_comps[a_comp][1]));
    }


    inline unsigned int size(unsigned int a_comp) const
    {
        PROTO_ASSERT((a_comp < 3),
        "CInterval::size(unsigned int a_comp) \
        invalid for a_comp = %i. a_comp must be in [0,3).",
        a_comp);
        return (m_comps[a_comp][1] - m_comps[a_comp][0] + 1);
    }

    inline bool operator==(CInterval& a_rhs)
    {
        bool equal = true;
        equal &= (m_comps[0][0] == a_rhs.m_comps[0][0]);
        equal &= (m_comps[1][0] == a_rhs.m_comps[1][0]);
        equal &= (m_comps[2][0] == a_rhs.m_comps[2][0]);
        equal &= (m_comps[0][1] == a_rhs.m_comps[0][1]);
        equal &= (m_comps[1][1] == a_rhs.m_comps[1][1]);
        equal &= (m_comps[2][1] == a_rhs.m_comps[2][1]);
        return equal;
    }

    inline void print() const
    {
      std::cout << "{ ";
      std::cout << "{" << m_comps[0][0] << ", " << m_comps[0][1] << "}";
      std::cout << ", {" << m_comps[1][0] << ", " << m_comps[1][1] << "}";
      std::cout << ", {" << m_comps[2][0] << ", " << m_comps[2][1] << "}";
      std::cout << " }";
    }

private:
    unsigned int m_comps[3][2]; 
};

inline std::ostream& operator<<(std::ostream& a_os,
                                const CInterval& a_int)
{
      a_os << "{ ";
      a_os << "{" << a_int.low(0) << ", " << a_int.high(0) << "}";
      a_os << ", {" << a_int.low(1) << ", " << a_int.high(1) << "}";
      a_os << ", {" << a_int.low(2) << ", " << a_int.high(2) << "}";
      a_os << " }";
      return a_os;
}


  template<typename T, unsigned int C=1 ,MemType MEMTYPE=MEMTYPE_DEFAULT, unsigned char D=1, unsigned char E=1>
  class Var
  {
  public:


#ifdef PROTO_CUDA //__CUDA_ARCH__
    __device__
    inline T& getValDevice(unsigned int a_c, unsigned char a_d = 0, unsigned char a_e = 0)
    {
      int idx = threadIdx.x + blockIdx.x*blockDim.x;
      int idy = blockIdx.y;
#if DIM == 3
      int idz = blockIdx.z;
      return (m_ptrs[a_c + C*a_d + C*D*a_e])[idx + boxDimX * (idy + idz * boxDimY)];
#else
      return (m_ptrs[a_c + C*a_d + C*D*a_e])[idx + idy * boxDimX];
#endif
    }

#ifdef PROTO_HIP
    __host__
    inline T& getValDevice(unsigned int a_c, unsigned char a_d = 0, unsigned char a_e = 0)
    {
      return *(m_ptrs[a_c + C*a_d + C*D*a_e]);
    }
#endif
#else
    inline T& getValDevice(unsigned int a_c, unsigned char a_d = 0, unsigned char a_e = 0)
    {
      return *(m_ptrs[a_c + C*a_d + C*D*a_e]);
    }
#endif

    CUDA_DECORATION
    __attribute__((always_inline)) T& operator()(unsigned int a_c, unsigned char a_d = 0, unsigned char a_e = 0)  
    {
#ifdef PROTO_CUDA
      if(MEMTYPE==MemType::DEVICE)
    return getValDevice(a_c,a_d,a_e);
      else
#endif
    return *(m_ptrs[a_c + C*a_d + C*D*a_e]);
    }


#ifdef PROTO_CUDA //__CUDA_ARCH__
    __device__
    inline const T& getValDevice(unsigned int a_c, unsigned char a_d = 0, unsigned char a_e = 0) const
    {
      int idx = threadIdx.x + blockIdx.x*blockDim.x;
      int idy = blockIdx.y;
#if DIM == 3
      int idz = blockIdx.z;
      return (m_ptrs[a_c + C*a_d + C*D*a_e])[idx + boxDimX * (idy + idz * boxDimY)];
#else
      return (m_ptrs[a_c + C*a_d + C*D*a_e])[idx + idy * boxDimX];
#endif
    }

#ifdef PROTO_HIP
    __host__
    inline const T& getValDevice(unsigned int a_c, unsigned char a_d = 0, unsigned char a_e = 0) const
    {
      return *(m_ptrs[a_c + C*a_d + C*D*a_e]);
    }
#endif
#else
    inline const T& getValDevice(unsigned int a_c, unsigned char a_d = 0, unsigned char a_e = 0) const
    {
      return *(m_ptrs[a_c + C*a_d + C*D*a_e]);
    }
#endif
    
    CUDA_DECORATION
    __attribute__((always_inline)) const T& operator()(unsigned int a_c, unsigned char a_d = 0, unsigned char a_e = 0) const 
    {
#ifdef PROTO_CUDA
      if(MEMTYPE==MemType::DEVICE)
    return getValDevice(a_c,a_d,a_e);
      else
#endif
        return *(m_ptrs[a_c + C*a_d + C*D*a_e]);
    }
 
    CUDA_DECORATION
    inline Var&  operator+=(unsigned int& a_increment) 
    {
      for (int ii = 0; ii < C*D*E; ii++)
      {
        m_ptrs[ii] += a_increment;
      }
      return *this;
    }

    // test
    CUDA_DECORATION
    inline Var&  operator+=(Point& a_p)
    {
#if DIM == 3
      unsigned int shift = a_p[0] + (a_p[1] + a_p[2]*boxDimY)*boxDimX;
#else
      unsigned int shift = a_p[0] + a_p[1]*boxDimX;
#endif
      return *this+=shift;
    }

    CUDA_DECORATION
    inline Var&  operator+=(const Point& a_p)
    {
#if DIM == 3
      unsigned int shift = a_p[0] + (a_p[1] + a_p[2]*boxDimY)*boxDimX;
#else
      unsigned int shift = a_p[0] + a_p[1]*boxDimX;
#endif
      return *this+=shift;
    }
   
    CUDA_DECORATION
    inline Var& operator++() 
    {
      for (int ii = 0; ii < C*D*E; ii++)
      {
        ++m_ptrs[ii];
      }
      return *this;
    }

    CUDA_DECORATION
    inline Var& operator--() 
    {
      for (int ii = 0; ii < C*D*E; ii++)
      {
        --m_ptrs[ii];
      }
      return *this;
    }

    unsigned int boxDimX;
    unsigned int boxDimY;
    unsigned int subBoxDimX;
    unsigned int subBoxDimY;
    T* m_ptrs[C*D*E];
  };


  template <class T=double, unsigned int C=1, MemType MEMTYPE=MEMTYPE_DEFAULT, unsigned char D=1, unsigned char E=1>
  class BoxData 
  {
  public:

    typedef Var<T,C,MEMTYPE,D,E> reference;

  
      BoxData(const BoxData<T,C,MEMTYPE,D,E>& a_src) = delete;
      BoxData& operator=(const BoxData<T,C,MEMTYPE,D,E>& a_src) = delete;


    BoxData();


    explicit BoxData(const Box& a_box, bool a_stackAllocation=DEFAULT_USE_STACK);
    

    BoxData(const Box& a_box, T a_init);
    


      BoxData(BoxData<T,C,MEMTYPE,D,E>&& a_src);
    

    void define(const Box& a_box, bool a_stackAllocation=DEFAULT_USE_STACK);



    void defineMalloc(const Box& a_box);
   


    void defineStack(const Box& a_box);
    

    template<unsigned int ncomp>
    void define(BoxData<T,ncomp,MEMTYPE,1,1>& a_input,
                unsigned int      & a_comp)
    {
      T* dataptr = a_input.dataPtr(a_comp);
      m_data   = a_input.getData();
      m_box    = a_input.box();
      m_rawPtr = dataptr;
    }           
   


    BoxData(const T* a_ptr, const Box& a_box, int a_ncomp = C);



    void define(const T* a_ptr, const Box& a_box, int a_ncomp = C);



    BoxData(LazyStencil<T,C,MEMTYPE,D,E>&& a_op);
   


    BoxData(shared_ptr<T> a_data,const T* a_ptr,const Box& a_box);


    ~BoxData(); 
  
   


    template<BoxDataOp>
    void operatorT(const BoxData<T,C,MEMTYPE,D,E>& a_src);


    template<BoxDataOp>
    void operatorT(const T a_scale);


    BoxData& operator=(BoxData<T,C,MEMTYPE,D,E>&& a_src);


    void copyTo(BoxData<T,C,MEMTYPE,D,E>& a_dest) const;


    void copyTo(BoxData<T,C,MEMTYPE,D,E>&   a_dest,
                const Box&           a_srcBox,
                const Point&        a_destShift = Point::Zeros()) const;
    

    void copyTo(BoxData<T,C,MEMTYPE,D,E>& a_dest,
                const Box&    a_srcBox,
                const Box&    a_destBox) const;


    template <unsigned int Cdest,
              unsigned char Ddest,
              unsigned char Edest>
    void copyTo(BoxData<T,Cdest,MEMTYPE,Ddest,Edest>& a_dest,
                const Box&    a_srcBox,
                CInterval    a_srcComps,
                const Point& a_destShift,
                CInterval    a_destComps) const;


    template<unsigned int Csrc>
    void copy(const BoxData<T,Csrc,MEMTYPE,D,E>& a_dsrc,
              const Box&                  a_srcBox,
              unsigned int               a_srcComp,
              const Box&                  a_destBox,
              unsigned int               a_destComp,
              unsigned int               a_numcomp);
    
  

    

    CUDA_DECORATION
    inline T& operator()(const Point& a_pt,
                        unsigned int  a_c = 0,
                        unsigned char a_d = 0,
                        unsigned char a_e = 0) const
    {
      PROTO_ASSERT(m_box.contains(a_pt),
                   "BoxData::operator(Point a_pt, uint a_c, uchar a_d, uchar a_e) invalid for a_pt not in this->box()");
      PROTO_ASSERT((a_c < C),
                   "BoxData::operator(Point a_pt, uint a_c, uchar a_d, uchar a_e) invalid for a_c = %i. a_c must be in [0,%i)",a_c,C);
      PROTO_ASSERT((a_d < D),
                   "BoxData::operator(Point a_pt, uint a_c, uchar a_d, uchar a_e) invalid for a_d = %i. a_d must be in [0,%i)",a_d,D);
      PROTO_ASSERT( (a_e < E),
                   "BoxData::operator(Point a_pt, uint a_c, uchar a_d, uchar a_e) invalid for a_e = %i. a_e must be in [0,%i)",a_e,E);
      return m_rawPtr[index(a_pt,a_c,a_d,a_e)];
    }
  

    inline T* operator[](unsigned int a_index); 
    


    inline const T* operator[](unsigned int a_index) const; 
   


    inline const T* data(const Point& a_p,
                         unsigned int a_c = 0,
                         unsigned int a_d = 0,
                         unsigned int a_e = 0) const
    {
      size_t ioff = index(a_p, a_c, a_d, a_e);
      return &(m_rawPtr[ioff]);
    }


    inline T* data() 
    {
       return m_rawPtr;
    }


    inline const T* data() const
    {
       return m_rawPtr;
    }


    inline T* data(const Point& a_p,
                   unsigned int a_c = 0,
                   unsigned int a_d = 0,
                   unsigned int a_e = 0)
    {
      size_t ioff  = index(a_p, a_c, a_d, a_e);
      return &(m_rawPtr[ioff]);
    }
    

    inline const T* dataPtr(unsigned int a_c = 0,
                            unsigned int a_d = 0,
                            unsigned int a_e = 0) const
    {
      Point pt = m_box.low();
      size_t ioff = index(pt, a_c, a_d, a_e);
      return &(m_rawPtr[ioff]);
    }


    inline T* dataPtr(unsigned int a_c = 0,
                      unsigned int a_d = 0,
                      unsigned int a_e = 0)
    {
      Point pt = m_box.low();
      size_t ioff = index(pt, a_c, a_d, a_e);
      return &(m_rawPtr[ioff]);
    }


    inline shared_ptr<T> getData() const
    {
      return m_data; 
    }


    CUDA_DECORATION
    inline size_t index(const Point a_pt,
                        unsigned int a_c = 0,
                        unsigned int a_d = 0,
                        unsigned int a_e = 0) const
    {
      PROTO_ASSERT(m_box.contains(a_pt),
                   "BoxData::operator(Point a_pt, uint a_c, uchar a_d, uchar a_e)\
                   invalid for a_pt not in this->box()");
      PROTO_ASSERT((a_c < C),
                   "BoxData::operator(Point a_pt, uint a_c, uchar a_d, uchar a_e)\
                    invalid for a_c = %i. a_c must be in [0,%i)",a_c,C);
      PROTO_ASSERT((a_d < D),
                   "BoxData::operator(Point a_pt, uint a_c, uchar a_d, uchar a_e)\
                    invalid for a_d = %i. a_d must be in [0,%i)",a_d,D);
      PROTO_ASSERT((a_e < E),
                   "BoxData::operator(Point a_pt, uint a_c, uchar a_d, uchar a_e)\
                    invalid for a_e = %i. a_e must be in [0,%i)",a_e,E);
      size_t m = m_box.size();
      size_t k = m_box.index(a_pt);
      return k+m*a_c+a_d*m*C+a_e*m*C*D;
    }
  

    // Get Domain Box
    inline Box box() const {return m_box;};
    
    // Get Size
    inline std::size_t size() const {return m_box.size()*C*D*E;};


    inline bool defined() const {return bool(m_data);};
   


    inline Var<T,C,MEMTYPE,D,E> var(const Point& a_pt);
    

    inline Var<T,C,MEMTYPE,D,E> var(const Point& a_pt) const;
 
    



    BoxData<T,C,MEMTYPE,D,E>& operator+=(const BoxData<T,C,MEMTYPE,D,E>& a_rhs);
    

    BoxData<T,C,MEMTYPE,D,E>& operator-=(const BoxData<T,C,MEMTYPE,D,E>& a_rhs);


    BoxData<T,C,MEMTYPE,D,E>& operator*=(const BoxData<T,C,MEMTYPE,D,E>& a_rhs);


    BoxData<T,C,MEMTYPE,D,E>& operator/=(const BoxData<T,C,MEMTYPE,D,E>& a_rhs);


    BoxData<T,C,MEMTYPE,D,E>& operator+=(T a_scale);


    BoxData<T,C,MEMTYPE,D,E>& operator-=(T scale);


    BoxData<T,C,MEMTYPE,D,E>& operator*=(T scale);


    BoxData<T,C,MEMTYPE,D,E>& operator/=(T scale);
  




    void setVal(const T& a_val);
  

    void setVal(const T&    a_val,
                const Box& a_box);


    void setVal(const T& a_val,
                const Box& a_box,
                int a_c,
                int a_d = 0,
                int a_e = 0);


    void absMax(Reduction<T,Abs>& a_Rxn) const;


    T absMax() const;


    T absMax(int a_c, int a_d = 0, int a_e = 0) const;


    T min() const;


    T min(int a_c, int a_d = 0, int a_e = 0) const;


    T max() const;


    T max(int a_c, int a_d = 0, int a_e = 0) const;


    T sum() const;


    T sum(int a_c, int a_d = 0, int a_e = 0) const;


    inline void shift(const Point a_pt)
    {
      m_box = m_box.shift(a_pt);
    };
  

    void linearOut(void*     a_buf,
                   const Box& a_box,
                   CInterval a_comps) const;

    void linearOut(void*                  a_buf,
                   const ::Proto::Box&      a_bx,
                   unsigned int a_startcomp,
                   unsigned int a_numcomps) const;



    void linearIn(void*     a_buf,
                  const Box& a_box,
                  CInterval a_comps);

    void linearIn(void*                  a_buf,
                  const ::Proto::Box&      a_bx,
                  unsigned int a_startcomp,
                  unsigned int a_numcomps);
 

   

    inline bool contains(CInterval a_interval) const
    {
      bool ret = true;
      ret &= (a_interval.high(0) < C);
      ret &= (a_interval.high(1) < D);
      ret &= (a_interval.high(2) < E);
      return ret;
    }
    

    template<unsigned int CC, unsigned char DD, unsigned char EE>
    inline bool isAlias(const BoxData<T, CC, MEMTYPE, DD, EE>& a_src) const
    {
      return (m_data.get() == a_src.m_data.get());
    }


    void print() const;


    void printData(int a_prec = 2) const;
    

    void printData(const Box& a_box, int a_prec = 2) const;


    void printData(const Box& a_box, int a_c, int a_d, int a_e, int a_prec = 2) const;

    // Friends
    
    template<class TT>
    friend class Stencil;
    
    template<class TT, unsigned int CC, unsigned char DD, unsigned char EE,MemType MM>
      friend BoxData<TT,CC,MM,DD,EE> alias(BoxData<TT,CC,MM,DD,EE>& a_original,
                                                const Point&              a_shift);

    template<class TT, unsigned int CC, unsigned char DD, unsigned char EE,MemType MM>
    friend const BoxData<TT,CC,MM,DD,EE> alias(const BoxData<TT,CC,MM,DD,EE>& a_original,
                                                    const Point&                a_shift);

   class Interval;
   class IntVect;
   size_t size(const Box&      a_box,
                      const size_t a_comps) const
   {
     return a_box.size() * (sizeof(T) * a_comps);
   }



   size_t charsize(const ::Proto::Box&      a_bx,
                          unsigned int a_startcomp,
                          unsigned int a_numcomps) const {

     // Box chbox(a_bx.low(), a_bx.high());

     //const size_t interv =  a_startcomp + a_numcomps + a_startcomp + 1;
     //return size(chbox, interv);
     return a_bx.size()*sizeof(T)*a_numcomps;
  }

  static int preAllocatable()
  {
    return 0; // static preAllocatable
  }

  static int memTypeAllocation()
  {
    return 1;
  }

  private:

    bool                     m_stackAlloc;
    Box                      m_box; 
    ::std::shared_ptr<T>    m_data; 
    T*                      m_rawPtr; 
  };




  template<class T, unsigned int C = 1, unsigned char D = 1, unsigned char E = 1, MemType MEMTYPE=MEMTYPE_DEFAULT>
  BoxData<T,C,MEMTYPE,D,E> alias(BoxData<T,C,MEMTYPE,D,E>& a_original, const Point& shift=Point::Zeros());


  template<class T, unsigned int C = 1, unsigned char D = 1, unsigned char E = 1, MemType MEMTYPE=MEMTYPE_DEFAULT>
  const BoxData<T,C,MEMTYPE,D,E> alias(const BoxData<T,C,MEMTYPE,D,E>& a_original, const Point& shift=Point::Zeros());


    template<typename T, unsigned int C, MemType MEMTYPE=MEMTYPE_DEFAULT, unsigned char D=1, unsigned char E=1>
    BoxData<T,1,MEMTYPE,1,1> slice(const BoxData<T,C,MEMTYPE,D,E>& a_src,
                                   unsigned int a_c,
                                   unsigned int a_d = 0,
                                   unsigned int a_e = 0);


  template<typename T, unsigned int C, unsigned char CC, MemType MEMTYPE=MEMTYPE_DEFAULT>
  BoxData<T,CC,MEMTYPE,1,1> slice(const BoxData<T,C,MEMTYPE,1,1>& a_src, unsigned int a_nstart);





    template<typename T, unsigned int C=1, unsigned char D=1, unsigned char E=1, MemType MEMTYPE=MEMTYPE_DEFAULT,
             typename Func, typename... Srcs>
  inline BoxData<T,C,MEMTYPE,D,E> forall(const Func& a_F, Srcs&&... a_srcs);

  template<typename T, unsigned int C=1, unsigned char D=1, unsigned char E=1, MemType MEMTYPE=MEMTYPE_DEFAULT,
           typename Func, typename... Srcs>
  inline BoxData<T,C,MEMTYPE,D,E> forallOp(unsigned long long int a_num_flops_point,   
                                           const char* a_timername, 
                                           const Func& a_F, Srcs&&... a_srcs);



  template<typename T, unsigned int C=1, unsigned char D=1, unsigned char E=1, MemType MEMTYPE=MEMTYPE_DEFAULT,
           typename Func, typename... Srcs>
  inline BoxData<T,C,MEMTYPE,D,E> forall(const Func& a_F, Box  a_box, Srcs&&... a_srcs); 

  template<typename T, unsigned int C=1, unsigned char D=1, unsigned char E=1, MemType MEMTYPE=MEMTYPE_DEFAULT,
           typename Func, typename... Srcs>
  inline BoxData<T,C,MEMTYPE,D,E> forallOp(unsigned long long int a_num_flops_point,   
                                   const char* a_timername, 
                                   const Func& a_F, Box  a_box, Srcs&&... a_srcs); 


  template<typename T, unsigned int C=1, unsigned char D=1, unsigned char E=1, MemType MEMTYPE=MEMTYPE_DEFAULT,
           typename Func, typename... Srcs>
  inline BoxData<T,C,MEMTYPE,D,E> forall_p(const Func& a_F, Srcs&&... a_srcs);


  template<typename T, unsigned int C=1, unsigned char D=1, unsigned char E=1, MemType MEMTYPE=MEMTYPE_DEFAULT,
           typename Func, typename... Srcs>
  inline BoxData<T,C,MEMTYPE,D,E> forallOp_p(unsigned long long int a_num_flops_point,   
                                     const char* a_timername, 
                                     const Func& a_F, Srcs&&... a_srcs);


  template<typename T, unsigned int C=1, unsigned char D=1, unsigned char E=1, MemType MEMTYPE=MEMTYPE_DEFAULT,
           typename Func, typename... Srcs>
  inline BoxData<T,C,MEMTYPE,D,E> forall_p(const Func& a_F, Box a_box, Srcs&&... a_srcs); 

  template<typename T, unsigned int C=1, unsigned char D=1, unsigned char E=1, MemType MEMTYPE=MEMTYPE_DEFAULT,
           typename Func, typename... Srcs>
  inline BoxData<T,C,MEMTYPE,D,E> forallOp_p(unsigned long long int a_num_flops_point,   
                                     const char* a_timername, 
                                     const Func& a_F, Box a_box, Srcs&&... a_srcs); 


  template<typename Func, typename... Srcs>
  inline void forallInPlace(const Func& a_F, Srcs&&... a_srcs);


  template<typename Func, typename... Srcs>
  inline void forallInPlaceOp(unsigned long long int a_num_flops_point,   
                              const char* a_timername, 
                              const Func& a_F, Srcs&&... a_srcs);


  template<typename Func, typename... Srcs>
  inline void forallInPlace(const Func& a_F, Box a_box, Srcs&&... a_srcs); 

  template<typename Func, typename... Srcs>
  inline void forallInPlaceOp(unsigned long long int a_num_flops_point,   
                              const char* a_timername, 
                              const Func& a_F, Box a_box, Srcs&&... a_srcs); 



  template<typename Func, typename... Srcs>
  inline void forallInPlace_p(const Func& a_F, Srcs&&... a_srcs);

  template<typename Func, typename... Srcs>
  inline void forallInPlaceOp_p(unsigned long long int a_num_flops_point,   
                                const char* a_timername, 
                                const Func& a_F, Srcs&&... a_srcs);


  template<typename Func, typename... Srcs>
  inline void forallInPlace_p(const Func& a_F,  Box a_box, Srcs&&... a_srcs); 

  template<typename Func, typename... Srcs>
  inline void forallInPlaceOp_p(unsigned long long int a_num_flops_point,   
                                const char* a_timername, 
                                const Func& a_F,  Box a_box, Srcs&&... a_srcs); 

  // End of BoxData Doxygen Module
//========================================================================
//CUDA STUFF
//========================================================================
#ifdef PROTO_CUDA
  template<typename Func, typename... Srcs>
  inline void protoForall(const Func& a_F, Box a_box, Srcs&&... a_srcs); 



  template<typename Func, typename... Srcs>
  inline void protoForallStream(protoStream_t& a_stream, const Func& a_F, Box a_box, Srcs&&... a_srcs); 

  template<typename Func, typename... Srcs>
  inline void protoForall_p(const Func& a_F, Box a_box, Srcs&&... a_srcs); 
#endif

#include "implem/Proto_BoxDataImplem.H"
} //end Proto namespace
#endif //end include guard