array3.h

#ifndef ARRAY3_H
#define ARRAY3_H

#include "array1.h"
#include <algorithm>
#include <cassert>
#include <iostream>
#include <vector>

template<class T, class ArrayT=std::vector<T> >
struct Array3
{
   // STL-friendly typedefs

   typedef typename ArrayT::iterator iterator;
   typedef typename ArrayT::const_iterator const_iterator;
   typedef typename ArrayT::size_type size_type;
   typedef long difference_type;
   typedef T& reference;
   typedef const T& const_reference;
   typedef T value_type;
   typedef T* pointer;
   typedef const T* const_pointer;
   typedef typename ArrayT::reverse_iterator reverse_iterator;
   typedef typename ArrayT::const_reverse_iterator const_reverse_iterator;

   // the actual representation

   int ni, nj, nk;
   ArrayT a;

   // the interface

   Array3(void)
      : ni(0), nj(0), nk(0)
   {}

   Array3(int ni_, int nj_, int nk_)
      : ni(ni_), nj(nj_), nk(nk_), a(ni_*nj_*nk_)
   { assert(ni_>=0 && nj_>=0 && nk_>=0); }

   Array3(int ni_, int nj_, int nk_, ArrayT& a_)
      : ni(ni_), nj(nj_), nk(nk_), a(a_)
   { assert(ni_>=0 && nj_>=0 && nk_>=0); }

   Array3(int ni_, int nj_, int nk_, const T& value)
      : ni(ni_), nj(nj_), nk(nk_), a(ni_*nj_*nk_, value)
   { assert(ni_>=0 && nj_>=0 && nk_>=0); }

   Array3(int ni_, int nj_, int nk_, const T& value, size_type max_n_)
      : ni(ni_), nj(nj_), nk(nk_), a(ni_*nj_*nk_, value, max_n_)
   { assert(ni_>=0 && nj_>=0 && nk_>=0); }

   Array3(int ni_, int nj_, int nk_, T* data_)
      : ni(ni_), nj(nj_), nk(nk_), a(ni_*nj_*nk_, data_)
   { assert(ni_>=0 && nj_>=0 && nk_>=0); }

   Array3(int ni_, int nj_, int nk_, T* data_, size_type max_n_)
      : ni(ni_), nj(nj_), nk(nk_), a(ni_*nj_*nk_, data_, max_n_)
   { assert(ni_>=0 && nj_>=0 && nk_>=0); }

   ~Array3(void)
   {
#ifndef NDEBUG
      ni=nj=0;
#endif
   }

   const T& operator()(int i, int j, int k) const
   {
      assert(i>=0 && i<ni && j>=0 && j<nj && k>=0 && k<nk);
      return a[i+ni*(j+nj*k)];
   }

   T& operator()(int i, int j, int k)
   {
     if (i < 0 || i >= ni || j < 0 || j >= nj || k < 0 || k >= nk) {
       std::cout << "Illegal indices" << std::endl;
       std::cout << "IJK: " << i << " " << j << " " << k << std::endl;
       std::cout << "Dims: " << ni << " " << nj << " " << nk << std::endl;
     }
      assert(i>=0 && i<ni && j>=0 && j<nj && k>=0 && k<nk);
      return a[i+ni*(j+nj*k)];
   }

   bool operator==(const Array3<T>& x) const
   { return ni==x.ni && nj==x.nj && nk==x.nk && a==x.a; }

   bool operator!=(const Array3<T>& x) const
   { return ni!=x.ni || nj!=x.nj || nk!=x.nk || a!=x.a; }

   bool operator<(const Array3<T>& x) const
   {
      if(ni<x.ni) return true; else if(ni>x.ni) return false;
      if(nj<x.nj) return true; else if(nj>x.nj) return false;
      if(nk<x.nk) return true; else if(nk>x.nk) return false;
      return a<x.a;
   }

   bool operator>(const Array3<T>& x) const
   {
      if(ni>x.ni) return true; else if(ni<x.ni) return false;
      if(nj>x.nj) return true; else if(nj<x.nj) return false;
      if(nk>x.nk) return true; else if(nk<x.nk) return false;
      return a>x.a;
   }

   bool operator<=(const Array3<T>& x) const
   {
      if(ni<x.ni) return true; else if(ni>x.ni) return false;
      if(nj<x.nj) return true; else if(nj>x.nj) return false;
      if(nk<x.nk) return true; else if(nk>x.nk) return false;
      return a<=x.a;
   }

   bool operator>=(const Array3<T>& x) const
   {
      if(ni>x.ni) return true; else if(ni<x.ni) return false;
      if(nj>x.nj) return true; else if(nj<x.nj) return false;
      if(nk>x.nk) return true; else if(nk<x.nk) return false;
      return a>=x.a;
   }

   void assign(const T& value)
   { a.assign(value); }

   void assign(int ni_, int nj_, int nk_, const T& value)
   {
      a.assign(ni_*nj_*nk_, value);
      ni=ni_;
      nj=nj_;
      nk=nk_;
   }
    
   void assign(int ni_, int nj_, int nk_, const T* copydata)
   {
      a.assign(ni_*nj_*nk_, copydata);
      ni=ni_;
      nj=nj_;
      nk=nk_;
   }
    
   const T& at(int i, int j, int k) const
   {
      assert(i>=0 && i<ni && j>=0 && j<nj && k>=0 && k<nk);
      return a[i+ni*(j+nj*k)];
   }

   T& at(int i, int j, int k)
   {
      assert(i>=0 && i<ni && j>=0 && j<nj && k>=0 && k<nk);
      return a[i+ni*(j+nj*k)];
   }

   const T& back(void) const
   { 
      assert(a.size());
      return a.back();
   }

   T& back(void)
   {
      assert(a.size());
      return a.back();
   }

   const_iterator begin(void) const
   { return a.begin(); }

   iterator begin(void)
   { return a.begin(); }

   size_type capacity(void) const
   { return a.capacity(); }

   void clear(void)
   {
      a.clear();
      ni=nj=nk=0;
   }

   bool empty(void) const
   { return a.empty(); }

   const_iterator end(void) const
   { return a.end(); }

   iterator end(void)
   { return a.end(); }

   void fill(int ni_, int nj_, int nk_, const T& value)
   {
      a.fill(ni_*nj_*nk_, value);
      ni=ni_;
      nj=nj_;
      nk=nk_;
   }
    
   const T& front(void) const
   {
      assert(a.size());
      return a.front();
   }

   T& front(void)
   {
      assert(a.size());
      return a.front();
   }

   size_type max_size(void) const
   { return a.max_size(); }

   reverse_iterator rbegin(void)
   { return reverse_iterator(end()); }

   const_reverse_iterator rbegin(void) const
   { return const_reverse_iterator(end()); }

   reverse_iterator rend(void)
   { return reverse_iterator(begin()); }

   const_reverse_iterator rend(void) const
   { return const_reverse_iterator(begin()); }

   void reserve(int reserve_ni, int reserve_nj, int reserve_nk)
   { a.reserve(reserve_ni*reserve_nj*reserve_nk); }

   void resize(int ni_, int nj_, int nk_)
   {
      assert(ni_>=0 && nj_>=0 && nk_>=0);
      a.resize(ni_*nj_*nk_);
      ni=ni_;
      nj=nj_;
      nk=nk_;
   }

   void resize(int ni_, int nj_, int nk_, const T& value)
   {
      assert(ni_>=0 && nj_>=0 && nk_>=0);
      a.resize(ni_*nj_*nk_, value);
      ni=ni_;
      nj=nj_;
      nk=nk_;
   }

   void set_zero(void)
   { a.set_zero(); }

   size_type size(void) const
   { return a.size(); }

   void swap(Array3<T>& x)
   {
      std::swap(ni, x.ni);
      std::swap(nj, x.nj);
      std::swap(nk, x.nk);
      a.swap(x.a);
   }

   void trim(void)
   { a.trim(); }
};

// some common arrays

typedef Array3<double, Array1<double> > Array3d;
typedef Array3<float, Array1<float> > Array3f;
typedef Array3<long long, Array1<long long> > Array3ll;
typedef Array3<unsigned long long, Array1<unsigned long long> > Array3ull;
typedef Array3<int, Array1<int> > Array3i;
typedef Array3<unsigned int, Array1<unsigned int> > Array3ui;
typedef Array3<short, Array1<short> > Array3s;
typedef Array3<unsigned short, Array1<unsigned short> > Array3us;
typedef Array3<char, Array1<char> > Array3c;
typedef Array3<unsigned char, Array1<unsigned char> > Array3uc;

#endif