
// Type-secure C++ interface for multi-dimensional arrays.


// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - //
//                                                                     //
// template class `Section_': Adds a data field to `DAD'.              //
//                                                                     //
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - //

template<int R, class T>
class Section_ : public DAD {
public :

  Section_(const Group& _group, T* _data) :
#ifndef RESTRICTED_OVERLOADING
               DAD(R, _group), data(_data) {}
#else
               DAD(R, _group), data(_data), isArray(0) {}
#endif

  T* dat() const {return data ; }

  operator IOArg() const {return IOArg((DAD*) this, data, sizeof(T)) ; }
      // Conversion operator to support `aprintf' etc.

protected :
  Section_(const Group& _group) :
#ifndef RESTRICTED_OVERLOADING
               DAD(R, _group) {}
#else
               DAD(R, _group), isArray(1) {}
#endif
                   // Used by `Array' subclasses.

  T* data ;

#ifdef RESTRICTED_OVERLOADING
  int isArray ;
#endif
} ;


// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - //
//                                                                     //
// template classes `SectionR': Add subscripting operations to         //
//                              `Section_'.                            //
//                                                                     //
// template classes `ArrayR': Add element allocation to `SectionR'     //
//                                                                     //
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - //

// Zero-dimensional sections and arrays

#ifdef RESTRICTED_OVERLOADING
#define Array0 Section0
#endif

template<class T>
class Section0 : public Section_<0, T> {
public :
  // Constructors

  Section0(const Section0<T>& a) : Section_<0, T>(a.grp(), a.dat()) {}

  Section0(const Group& p, T* _data) : Section_<0, T>(p, _data) {}

  // Simple subscripting

  T& operator()() {
    return *(dat()) ;
  }

#ifndef RESTRICTED_OVERLOADING
protected :
  Section0(const Group& p) : Section_<0, T>(p) {}
               // Used by `Array0'.
#endif

#ifndef RESTRICTED_OVERLOADING
} ;

template<class T>
class Array0 : public Section0<T> {
#endif

public :

  Array0(const Group& p = apg) :
#ifndef RESTRICTED_OVERLOADING
             Section0<T>(p) {
#else
             Section_<0, T>(p) {
#endif

    data = new T ;
  }

#ifndef RESTRICTED_OVERLOADING
  ~Array0() {
    delete data ;
  }

private :
  Array0(const Array0<T>& a) : Section0<T>(a) {}  // Prohibit copying.

#else
  ~Array0() {
    if(isArray) delete data ;
  }
#endif
} ;


// One-dimensional sections and arrays

#ifdef RESTRICTED_OVERLOADING
#define Array1 Section1
#endif

template<class T>
class Section1 : public Section_<1, T> {
public :
  // Constructors

  Section1(const Section1<T>& a) :
               Section_<1, T>(a.grp(), a.dat()) {
    rngs [0] = a.rng(0) ;

    strs [0] = a.str(0) ;

    DAD :: ranges  = rngs ;
    DAD :: strides = strs ;
  }

  Section1(const Range& x_1,
           const Group& p,
	   const Stride& str_1,
	   T* _data) : Section_<1, T>(p, _data) {
    rngs [0] = x_1 ;

    strs [0] = str_1 ;

    DAD :: ranges  = rngs ;
    DAD :: strides = strs ;
  }

  // Simple subscripting

  T& operator()(const Location& s_1) const {
    return dat() [str(0).offset(s_1)] ;
  }

  // Section subscripting

  Section1<T> sect(const Range& x_1) const {
    return Section1<T>(x_1, grp(), str(0), dat()) ;
  }

  Section0<T> sect(const Location& s_1) const {
    return Section0<T>(grp() / s_1, dat() + str(0).offset(s_1)) ;
  }

#ifndef RESTRICTED_OVERLOADING
protected :
  Section1(const Group& p) : Section_<1, T>(p) {}  // Used by `Array1'
#endif

protected :
  Range  rngs [1] ;
  Stride strs [1] ;

#ifndef RESTRICTED_OVERLOADING
} ;

template<class T>
class Array1 : public Section1<T> {
#endif

public :

  Array1(const Range& x_1, const Group& p = apg) :
#ifndef RESTRICTED_OVERLOADING
             Section1<T>(p) {
#else
             Section_<1, T>(p) {
#endif
    rngs [0] = x_1 ;

    int size = 1 ; 
    strs [0] = x_1.stride(&size) ;

    DAD :: ranges  = rngs ;
    DAD :: strides = strs ;

    data = new T [size] ;
  }

#ifndef RESTRICTED_OVERLOADING
  ~Array1() {
    delete [] data ;
  }

private :
  Array1(const Array1<T>& a) : Section1<T>(a) {}  // Prohibit copying.

#else
  ~Array1() {
    if(isArray) delete [] data ;
  }
#endif
} ;


// Two-dimensional sections and arrays

#ifdef RESTRICTED_OVERLOADING
#define Array2 Section2
#endif

template<class T>
class Section2 : public Section_<2, T> {
public :
  // Constructors

  Section2(const Section2<T>& a) :
               Section_<2, T>(a.grp(), a.dat()) {
    rngs [0] = a.rng(0) ;
    rngs [1] = a.rng(1) ;

    strs [0] = a.str(0) ;
    strs [1] = a.str(1) ;

    DAD :: ranges  = rngs ;
    DAD :: strides = strs ;
  }

  Section2(const Range& x_1, const Range& x_2, const Group& p,
	   const Stride& str_1, const Stride& str_2, T* _data) :
               Section_<2, T>(p, _data) {
    rngs [0] = x_1 ;
    rngs [1] = x_2 ;

    strs [0] = str_1 ;
    strs [1] = str_2 ;

    DAD :: ranges  = rngs ;
    DAD :: strides = strs ;
  }

  // Simple subscripting

  T& operator()(const Location& s_1, const Location& s_2) const {
    return dat() [str(0).offset(s_1) + str(1).offset(s_2)] ;
  }

  // Section subscripting

  Section2<T> sect(const Range& x_1, const Range& x_2) const {
    return Section2<T>(x_1, x_2, grp(),
		       str(0), str(1), dat()) ;
  }

  Section1<T> sect(const Range& x_1, const Location& s_2) const {
    return Section1<T>(x_1, grp() / s_2,
		       str(0), dat() + str(1).offset(s_2)) ;
  }

  Section1<T> sect(const Location& s_1, const Range& x_2) const {
    return Section1<T>(x_2, grp() / s_1,
		       str(1), dat() + str(0).offset(s_1)) ;
  }

  Section0<T> sect(const Location& s_1, const Location& s_2) const {
    return Section0<T>(grp() / s_1 / s_2,
		       dat() + str(0).offset(s_1) + str(1).offset(s_2)) ;
  }

#ifndef RESTRICTED_OVERLOADING
protected :
  Section2(const Group& p) : Section_<2, T>(p) {}  // Used by `Array2'
#endif

protected :
  Range  rngs [2] ;
  Stride strs [2] ;

#ifndef RESTRICTED_OVERLOADING
} ;

template<class T>
class Array2 : public Section2<T> {
#endif

public :

  Array2(const Range& x_1, const Range& x_2, const Group& p = apg) :
#ifndef RESTRICTED_OVERLOADING
             Section2<T>(p) {
#else
             Section_<2, T>(p) {
#endif
    rngs [0] = x_1 ;
    rngs [1] = x_2 ;

    int size = 1 ; 
    strs [0] = x_1.stride(&size) ;
    strs [1] = x_2.stride(&size) ;

    DAD :: ranges  = rngs ;
    DAD :: strides = strs ;

    data = new T [size] ;
  }

#ifndef RESTRICTED_OVERLOADING
  ~Array2() {
    delete [] data ;
  }

private :
  Array2(const Array2<T>& a) : Section2<T>(a) {}  // Prohibit copying.

#else
  ~Array2() {
    if(isArray) delete [] data ;
  }
#endif
} ;


// Three-dimensional sections and arrays

#ifdef RESTRICTED_OVERLOADING
#define Array3 Section3
#endif

template<class T>
class Section3 : public Section_<3, T> {
public :
  // Constructors

  Section3(const Section3<T>& a) :
               Section_<3, T>(a.grp(), a.dat()) {
    rngs [0] = a.rng(0) ;
    rngs [1] = a.rng(1) ;
    rngs [2] = a.rng(2) ;

    strs [0] = a.str(0) ;
    strs [1] = a.str(1) ;
    strs [2] = a.str(2) ;

    DAD :: ranges  = rngs ;
    DAD :: strides = strs ;
  }

  Section3(const Range& x_1, const Range& x_2, const Range& x_3,
           const Group& p,
	   const Stride& str_1, const Stride& str_2, const Stride& str_3,
	   T* _data) : Section_<3, T>(p, _data) {
    rngs [0] = x_1 ;
    rngs [1] = x_2 ;
    rngs [2] = x_3 ;

    strs [0] = str_1 ;
    strs [1] = str_2 ;
    strs [2] = str_3 ;

    DAD :: ranges  = rngs ;
    DAD :: strides = strs ;
  }

  // Simple subscripting

  T& operator()(const Location& s_1, const Location& s_2, Location& s_3) const {
    return dat() [str(0).offset(s_1) + str(1).offset(s_2) +
                  str(2).offset(s_3)] ;
  }

  // Section subscripting

  Section3<T> sect(const Range& x_1, const Range& x_2,
                   const Range& x_3) const {
    return Section3<T>(x_1, x_2, x_3, grp(),
		       str(0), str(1), str(2), dat()) ;
  }

  Section2<T> sect(const Range& x_1, const Range& x_2,
                   const Location& s_3) const {
    return Section2<T>(x_1, x_2, grp() / s_3, 
		       str(0), str(1), dat() + str(2).offset(s_3)) ;
  }

  Section2<T> sect(const Range& x_1, const Location&
                   s_2, const Range& x_3) const {
    return Section2<T>(x_1, x_3, grp() / s_2,
		       str(0), str(2), dat() + str(1).offset(s_2)) ;
  }

  Section2<T> sect(const Location& s_1, const Range& x_2,
                   const Range& x_3) const {
    return Section2<T>(x_2, x_3, grp() / s_1,
		       str(1), str(2), dat() + str(0).offset(s_1)) ;
  }

  Section1<T> sect(const Range& x_1, const Location& s_2,
                   const Location& s_3) const {
    return Section1<T>(x_1, grp() / s_2 / s_3,
		       str(0),
                       dat() + str(1).offset(s_2) + str(2).offset(s_3)) ;
  }

  Section1<T> sect(const Location& s_1, const Range& x_2,
                   const Location& s_3) const {
    return Section1<T>(x_2, grp() / s_1 / s_3,
		       str(1),
                       dat() + str(0).offset(s_1) + str(2).offset(s_3)) ;
  }

  Section1<T> sect(const Location& s_1, const Location& s_2,
                   const Range& x_3) const {
    return Section1<T>(x_3, grp() / s_1 / s_2,
		       str(1),
                       dat() + str(0).offset(s_1) + str(1).offset(s_2)) ;
  }

  Section0<T> sect(const Location& s_1, const Location& s_2,
                   const Location& s_3) const {
    return Section0<T>(grp() / s_1 / s_2 / s_3,
		       dat() + str(0).offset(s_1) + str(1).offset(s_2) +
                               str(2).offset(s_3)) ;
  }

#ifndef RESTRICTED_OVERLOADING
protected :
  Section3(const Group& p) : Section_<3, T>(p) {}  // Used by `Array3'
#endif

protected :
  Range  rngs [3] ;
  Stride strs [3] ;

#ifndef RESTRICTED_OVERLOADING
} ;

template<class T>
class Array3 : public Section3<T> {
#endif

public :

  Array3(const Range& x_1, const Range& x_2, const Range& x_3,
         const Group& p = apg) :
#ifndef RESTRICTED_OVERLOADING
             Section3<T>(p) {
#else
             Section_<3, T>(p) {
#endif

    rngs [0] = x_1 ;
    rngs [1] = x_2 ;
    rngs [2] = x_3 ;

    int size = 1 ; 
    strs [0] = x_1.stride(&size) ;
    strs [1] = x_2.stride(&size) ;
    strs [2] = x_3.stride(&size) ;

    DAD :: ranges  = rngs ;
    DAD :: strides = strs ;

    data = new T [size] ;
  }

#ifndef RESTRICTED_OVERLOADING
  ~Array3() {
    delete [] data ;
  }

private :
  Array3(const Array3<T>& a) : Section3<T>(a) {}  // Prohibit copying.

#else
  ~Array3() {
    if(isArray) delete [] data ;
  }
#endif
} ;


// ----------------------------------------------------------------------
// Things to do
//
//      Find a better way to deal with ``restricted overloading''?
//
//      (Present solution is to jam `ArrayR' members into the `SectionR'
//      class and use a run-time flag, `isArray', to distinguish the
//      two kinds of object.  The `ArrayR' symbol becomes a macro that
//      expands to `SectionR'.  This is inefficient and will most likely
//      have weird effects in unexpected places.)

