Tensor3.h

#ifndef QNANO_NEW_TENSOR3_DEFINED_H
#define QNANO_NEW_TENSOR3_DEFINED_H

#include <complex>
#include <iostream>
#include <cstdlib>

template <typename T> class template_Tensor3{
public:
  //sizes: 
  size_t size[3];
  //data block:
  T *data;
  
  inline size_t get_index(size_t i, size_t j, size_t k)const{
    return (i*size[1]+j)*size[2]+k;
  }

public: 
  inline size_t get_totalsize()const{ return size[0]*size[1]*size[2];  }

protected:

  void allocate(){  data=new T[get_totalsize()];  }
  void deallocate(){ delete[] data; }
public:
  inline size_t get_size(size_t i)const{return size[i];}
  void resize(size_t i, size_t j, size_t k){
    size[0]=i; size[1]=j; size[2]=k;
    deallocate();
    allocate();
  }
  
  const T &operator()(size_t i, size_t j, size_t k)const{
    return data[get_index(i,j,k)];
  }
  T &operator()(size_t i, size_t j, size_t k){
    return data[get_index(i,j,k)];
  }
  const T &get_check(int i, int j, int k){
    if(i<0 || i>=(int)size[0] || j<0 || j>=(int)size[1] || k<0 || k>=(int)size[2] ){
      std::cerr<<"Tensor3: out of bounds!"<<std::endl;
      exit(1);
    }
    return data[get_index(i,j,k)];
  }
 

  //Constructors and copy operators
  void copy(const template_Tensor3 &other){
    resize(other.get_size(0),other.get_size(1),other.get_size(2));
    for(size_t i=0; i<get_totalsize(); i++){
      data[i]=other.data[i];
    }
  } 
  template_Tensor3 &operator=(const template_Tensor3 &other){
    copy(other);
    return *this;
  }
  template_Tensor3(const template_Tensor3 &other){
    copy(other);
  }
  template_Tensor3(){
    size[0]=size[1]=size[2]=0;
    allocate();
  } 
  template_Tensor3(size_t i, size_t j, size_t k){
    size[0]=i; size[1]=j; size[2]=k;
    allocate();
  }
  ~template_Tensor3(){
    deallocate();
  }
  

};


typedef template_Tensor3<std::complex<double> >Tensor3;
typedef template_Tensor3<double> Tensor3d;

#endif