tensor.h

Go to the documentation of this file.

// AMB 98/04
// last modified AMB 980428
// KG, 981205: added layout in order to provide an effective means to
// swap indices, Index iterators
// $Id: tensor.h,v 1.33.2.30 2022/11/03 17:28:11 garloff Exp $

#ifndef TBCI_TENSOR_H
#define TBCI_TENSOR_H  

#include "tbci/vector.h"
#include "tbci/index.h"

// Avoid -fguiding-decls
#if !defined(NO_GD) && !defined(AUTO_DECL)
# include "tbci/tensor_gd.h"
#endif

NAMESPACE_TBCI

#ifndef TBCI_DISABLE_EXCEPT
//#include "except.h" is included in basics.h
class TensErr : public NumErr
{
    public:
        TensErr() 
        : NumErr("Error in Tensor class") {}
        TensErr(const char* t, const long i = 0) 
        : NumErr(t, i) {}
        TensErr(const TensErr& te)
        : NumErr(te) {}
        virtual ~TensErr() throw() {}   
};
#endif

#ifdef PRAGMA_I
# pragma interface "tensor.h"
#endif

//----------------------------------------------------------------

template <typename T>
class CTensor
{
  
protected:
  Vector <T> data;              // linear storage
  unsigned rank;                // rank of tensor
  unsigned long noel;           // number of elements
  Index shape;                  // shape of tensor, eg. 2x3x4
  Index layout;                 // memory layout of tensor, eg. 12,4,1
  unsigned long calcsize (void)
  {
    noel = 1;
    for (REGISTER unsigned i = 0; i < rank; i++)
      noel *= shape (i);
    return noel;
  }
  unsigned calclayout ()
  {
    unsigned long fac = 1;
    for (REGISTER int i = rank-1; i >= 0; i--)
      {
        layout(i) = fac;
        fac *= shape(i);
      }
    return fac;
  }

public:

  typedef T value_type;
  typedef T element_type;
  typedef T aligned_value_type TALIGN(MIN_ALIGN2);

  CTensor () { rank = 0; }                      // default constructor
  explicit CTensor (const unsigned dim_rank);   // "square" constructor
  CTensor (const Index & ix);                   // standard constructor
  CTensor (const T & value, const Index & ix);  // standard constructor
  CTensor (const CTensor & ct);                 // copy constructor
  CTensor (vararg va, ...);                     // arbitrary no of arguments
  CTensor (const T, vararg va, ...);            // arb no of arg & ini value 

  ~CTensor () {}

  // mapping and reverse mapping
  unsigned long calc_offs (const Index & ix) const;
  Index calc_indx (const unsigned long i) const;

  // access
  //T & operator ()(const TVector<T> & tv) const
  //      { Vector<T> ix(tv); return data (calc_offs (ix)); }
        T & operator () (const Index & ix)
          { return data (calc_offs (ix)); }
  typename tbci_traits<T>::const_refval_type    
            operator () (const Index & ix) const
                { return data (calc_offs (ix)); }
        T & operator () (vararg va, ...);
  typename tbci_traits<T>::const_refval_type    
            operator () (vararg va, ...) const;

        T & operator () (const unsigned long i)
          { BCHKNR(i>=noel, TensErr, & op(): out of range, i);
            return data(i); }
  typename tbci_traits<T>::const_refval_type    
            operator () (const unsigned long i) const
          { BCHKNR(i>=noel, TensErr, op() c: out of range, i);
            return data(i); }

  const T& getcref (vararg va, ...) const;
     
        T& get_lin_idx (const unsigned long i) //const
          {  return (*this)(i); }
  typename tbci_traits<T>::const_refval_type    
           get_lin_idx (const unsigned long i) const
          {  return (*this)(i); }
  
#ifndef NO_POD
  unsigned lin_read (vararg va, ...);
#endif

  // swap indices
  CTensor<T>& transpose (const unsigned i, const unsigned j)
    {
      BCHK (i>=rank, TensErr, swap_idx: idx1 out of range, i, *this);
      BCHK (j>=rank, TensErr, swap_idx: idx2 out of range, j, *this);
      SWAP (shape(i), shape(j)); SWAP (layout(i), layout(j));
      return *this;
    }

  // assignment
  CTensor <T>& fill (const T & value)
    { data.fill (value); return *this; }

  // resizing equation operator!
  CTensor <T> & operator = (const CTensor <T> & ct)
  {
    rank = ct.rank; shape.resize(ct.shape);
    layout.resize (ct.layout);
    noel = ct.noel; data.resize(ct.data);
    return *this;
  }
  CTensor <T>& operator = (const T & value)
    { return this->fill (value); }

  //misc
  CTensor<T>& resize (const Index& ix)
    {
      shape.resize(ix); rank = ix.size();
      layout.resize(rank); calclayout ();
      calcsize(); data.resize (noel);
      return *this;
    }
  CTensor<T>& resize (const T& val, const Index& ix)
    {
      shape.resize(ix); rank = ix.size();
      layout.resize(rank); calclayout ();
      calcsize(); data.resize (val, noel);
      return *this;
    }
  CTensor<T>& resize (const CTensor<T>& ct)
    {
      shape.resize(ct.shape); layout.resize (ct.layout);
      rank = ct.rank; noel = ct.noel;
      data.resize (ct.data);
      return *this;
    }

  T max (unsigned long& pos) const { return data.max (pos); }
  T min (unsigned long& pos) const { return data.min (pos); }
  
  T max () const { return data.max (); }
  T min () const { return data.min (); }

  T trace () const { return 0; } // ????????????????????????????

  // compare
  bool operator == (const CTensor < T > &ct)
    {
      if (shape != ct.shape) return false;
      if (layout == ct.layout) {
        if (data == ct.data) return true; else return false; }
      // layout is different, now we are in trouble
      Index ix(0, rank);
      while (ix(0) < shape(0))
        {
          if ((*this)(ix) != ct(ix)) return false;
          ix.next_idx (shape);
        }
      return true;
    }

  bool operator != (const CTensor < T > &ct)
    { return !(this->operator == (ct)); }
  
  // info
  unsigned long lin_size (void) const
    { return noel; }
  Index index_size (void) const
    { return shape; }
  unsigned rank_size (void) const
    { return rank; }
  
  // i/o operations
  friend STD__ ostream & operator << FGD (STD__ ostream&, const CTensor <T>&);
  friend STD__ istream & operator >> FGD (STD__ istream&, CTensor <T>&);

};

// class CTensor definitions

template < typename T >
inline CTensor < T >::CTensor (const unsigned dim_rank) : 
  rank (dim_rank), shape (dim_rank, dim_rank), layout (dim_rank)
  // "square" constructor
{
  calcsize (); calclayout ();
  data.resize (noel);
}

template < typename T >
inline CTensor < T >::CTensor (const Index & ix) : 
  rank (ix.size()), shape (ix), layout (ix.size())
  // standard constructor
{
  calcsize (); calclayout ();
  data.resize (noel);
}

template < typename T >
inline CTensor < T >::CTensor (const T & value, const Index & ix) : 
   rank (ix.size()), shape (ix), layout (ix.size())
  // standard constructor        
{
  //rank = shape.size ();
  calcsize (); calclayout ();
  data.resize (value, noel);
}

// copy constructor
template < typename T >
inline CTensor < T >::CTensor (const CTensor & ct) 
  : data(ct.data), rank(ct.rank), noel(ct.noel), 
    shape(ct.shape), layout(ct.layout)
{}

// arb. arg. & ini
template < typename T >
inline CTensor < T >::CTensor (const T value, vararg va, ...) :
  rank (va), shape ((unsigned)va), layout ((unsigned)va)
{
  va_list vl;
  va_start (vl, va);
  for (unsigned i = 0; i < (unsigned)va; ++i) 
        shape (i) = va_arg (vl, unsigned);
  va_end (vl);
  calcsize (); calclayout ();
  data.resize (value, noel);
}

// arb. arg.
template < typename T >
inline CTensor <T>::CTensor (vararg va, ...) :
  rank (va), shape ((unsigned)va), layout ((unsigned)va)
{
  va_list vl;
  va_start (vl, va);
  for (unsigned i = 0; i < rank; ++i) 
        shape (i) = va_arg (vl, unsigned);
  va_end (vl);
  calcsize (); calclayout ();
  data.resize (noel);
}

#ifndef NO_POD
template < typename T >
unsigned CTensor < T >::lin_read (vararg va, ...)
{
  va_list vl;
  va_start (vl, va);
  unsigned i;
#if !defined(__clang__) || !defined(CPLX)
  for (i = 0; i < (unsigned)va; ++i)
        data (i) = va_arg (vl, T);
#else
  throw TensErr("vararg not supported for cplx in clang");
#warning no vararg support with cplx numbers and clang
#endif
  va_end (vl);
  return i;
}
#endif

//Calculates eff. index
template < typename T >
inline unsigned long CTensor < T >::calc_offs (const Index & ix) const
{
  BCHK (rank != ix.size(), TensErr, Rank of Tensor and Index size dont match, ix.size(), 0);
  unsigned long offs (0);
  for (REGISTER unsigned i = 0; i < rank; i++)
    {
      offs += ix (i) * layout (i);
    }
  BCHK (offs >= noel, TensErr, Too large index, offs, offs);
  return offs;
}

#if 0
template < typename T >
inline Index CTensor < T >::calc_indx (const unsigned long i) const
{
  Index ix (rank);
  unsigned long mod = i;
  for (unsigned r = rank; r > 0; r--)
    {
      ix (r - 1) = mod % shape (r - 1);
      mod /= shape (r - 1);
    }
  BCHK (mod > 0, TensErr, Mapping error, mod, ix);
  return ix;
}
#endif

/* This function is inefficient by design! */
template < typename T >
inline Index CTensor < T >::calc_indx (const unsigned long i) const
{
  Index ix (rank); ix.fill (0);
  unsigned long mod = i;
  BCHK (i >= noel, TensErr, Linear index too large, i, ix);
  // tricky stuff to start here ...
  unsigned long maxdiv = noel;
  while (mod)
    {
      // find largest divisor
      unsigned long div = 0; unsigned ixnr = 0;
      for (unsigned i = 0; i < rank; i++)
        if (layout(i) > div && layout(i) < maxdiv && shape(i) > 1)
          { div = layout(i); ixnr = i; }
      maxdiv = div;
      ix(ixnr) = mod / div; mod %= div;
    }
  return ix;
}


// formatted output
template < typename T >
STD__ ostream & operator << (STD__ ostream & os, const CTensor < T > & ct)
{
  if (ct.rank > 1)
    {
      //os << "\n";
      Index ix (0, ct.rank);
      for (unsigned long i = 0; i < ct.noel; i++)
        {
          os << ct(ix) << "\t"; ix.next_idx (ct.shape);
          for (unsigned j = ct.rank - 1; j > 0; j--)
            if (ix(j) == 0) os << "\n"; else break;
        }
    }
  else
    os << ct.data;
  return os;
}

template <typename T>
STD__ istream& operator >> (STD__ istream & is, CTensor<T>& ct)
{
  Index ix (0, ct.rank);
  for (unsigned long i = 0; i < ct.noel; i++)
    {
      is >> ct(ix); ix.next_idx (ct.shape);
    }
  return is;
}

//access
#if 0
template <typename T>
T& CTensor <T>::operator () (vararg va, ...) //const
{
  Index ix ((unsigned)va);
  va_list vl;
  va_start (vl, va);
  for (unsigned i = 0; i < (unsigned)va; ++i) 
        ix (i) = va_arg (vl, int);
  va_end (vl);
  return (*this) (ix);
}

template <typename T>
typename tbci_traits<T>::const_refval_type
CTensor <T>::operator () (vararg va, ...) const
{
  Index ix ((unsigned)va);
  va_list vl;
  va_start (vl, va);
  for (unsigned i = 0; i < (unsigned)va; ++i) 
        ix (i) = va_arg (vl, int);
  va_end (vl);
  return (*this) (ix);
}

#else

template <typename T>
T& CTensor <T>::operator () (vararg va, ...) //const
{
  BCHK((unsigned)va != rank, TensErr, vararg no must be eq to rank, (unsigned)va, data(0));
  unsigned lin_idx = 0;
  va_list vl;
  va_start (vl, va);
  for (unsigned i = 0; i < rank; ++i) 
        lin_idx += va_arg (vl, int) * layout(i);
  va_end (vl);
  return (data (lin_idx));
}

template <typename T>
typename tbci_traits<T>::const_refval_type
CTensor <T>::operator () (vararg va, ...) const
{
  BCHK((unsigned)va != rank, TensErr, vararg no must be eq to rank, (unsigned)va, data(0));
  unsigned lin_idx = 0;
  va_list vl;
  va_start (vl, va);
  for (unsigned i = 0; i < rank; ++i) 
        lin_idx += va_arg (vl, int) * layout(i);
  va_end (vl);
  return (data (lin_idx));
}
#endif


template <typename T>
const T& CTensor <T>::getcref (vararg va, ...) const
{
  BCHK((unsigned)va != rank, TensErr, vararg no must be eq to rank, (unsigned)va, data.getcref(0));
  unsigned lin_idx = 0;
  va_list vl;
  va_start (vl, va);
  for (unsigned i = 0; i < rank; ++i) 
        lin_idx += va_arg (vl, int) * layout(i);
  va_end (vl);
  return data.getcref (lin_idx);
}


//-----------------------------------------------------------------------------


template<typename T>
class Tensor : public CTensor<T>
{
  
protected:
  // additional data structure
  
public:
  
  typedef T value_type;
  typedef T element_type;
  typedef T aligned_value_type TALIGN(MIN_ALIGN2);
  // constructors
  Tensor () : CTensor <T> () {}
  Tensor (const unsigned dim_rank) : CTensor <T> (dim_rank) {}
  Tensor (const Index & ix) : CTensor <T> (ix) {}
  Tensor (const T & value, const Index & ix) : CTensor <T> (value, ix) {}
  Tensor (const Tensor <T>& ct) : CTensor <T> (ct) {}
  Tensor (vararg va, ...);
  Tensor (const T, vararg va, ...);
  ~Tensor () {}
  
  // basic arithmetics
  Tensor<T>& operator += (const Tensor<T> &);
  Tensor<T>& operator -= (const Tensor<T> &);

  Tensor<T> operator + (const Tensor<T> &) const;
  Tensor<T> operator - (const Tensor<T> &) const;

  Tensor<T>& operator *= (const T m) { this->data *= m; return *this; }
  Tensor<T>& operator /= (const T m) { this->data /= m; return *this; }
  Tensor<T>& operator += (const T m) { this->data += m; return *this; }
  Tensor<T>& operator -= (const T m) { this->data -= m; return *this; }

  Tensor<T> operator * (const T m) const
  {   
        Tensor<T> r (this->shape); 
        for (unsigned long i = 0; i < this->noel; ++i)
                r.data(i) = this->data(i) * m;
        return r;
  }


  //for friend NOINST Tensor<T> operator * FGD (const T, const Tensor<T> &);
  Tensor<T> mult (const T) const;

  Tensor<T> operator / (const T m) const
  {   
        Tensor<T> r (this->shape); 
        for (unsigned long i = 0; i < this->noel; ++i) 
                r.data(i) = this->data(i) / m;
        return r;
  }

  Tensor<T> operator - () const
  {
        Tensor<T> r (this->shape);
        for (unsigned long i=0; i < this->noel; ++i) 
                r.data(i) = - this->data(i);
        return r;
  }

  // tensor arithmetics
  Tensor<T> contract (const unsigned, const unsigned) const;    // contraction over 2 indices
  // direct tensor product
  Tensor<T> drctmul (const Tensor<T>&) const;
  friend Tensor<T> FRIEND_TBCI2__ dctmul FGD (const Tensor<T>&, const Tensor<T>&);
  // multiplication with contraction over 2 indices
  Tensor<T> cntrmul (const Tensor<T>&, const unsigned, const unsigned) const;
  friend Tensor<T> FRIEND_TBCI2__ ctrmul FGD (const Tensor<T>&, const Tensor<T>&, const unsigned, const unsigned);
  // multiplication w/ metric results in differently ordered indices comp. to cntrmul
  friend Tensor<T> FRIEND_TBCI2__ metrmul FGD (const Tensor<T>&, const Tensor<T>&, const unsigned, const unsigned);
  
};


//implementations of class Tensor


// arb. arg. & ini
template < typename T >
inline Tensor < T >::Tensor (const T value, vararg va, ...)
  //: CTensor<T>::rank (va.num), CTensor<T>::layout (va.num), CTensor<T>::shape (va.num)
{
  this->rank = va;
  this->layout.resize ((unsigned)va); 
  this->shape.resize  ((unsigned)va);
  va_list vl;
  va_start (vl, va);
  for (unsigned i = 0; i < (unsigned)va; ++i) 
        this->shape (i) = va_arg (vl, unsigned);
  va_end (vl);
  this->calcsize (); this->calclayout ();
  this->data.resize (value, this->noel);
}

// arb. arg.
template < typename T >
inline Tensor <T>::Tensor (vararg va, ...)
  //: CTensor<T>::rank (va.num), CTensor<T>::layout (va.num), CTensor<T>::shape (va.num)
{
  this->rank = va;
  this->layout.resize ((unsigned)va); 
  this->shape.resize  ((unsigned)va);
  va_list vl;
  va_start (vl, va);
  for (unsigned i = 0; i < this->rank; ++i) 
        this->shape (i) = va_arg (vl, unsigned);
  va_end (vl);
  this->calcsize (); this->calclayout ();
  this->data.resize (this->noel);
}


template <typename T>
Tensor<T> Tensor <T>::contract (const unsigned i1, const unsigned i2) const
{
  // Calculate shape
  Index ix (idx_remove2 (this->shape, i1, i2));
  Tensor < T > t (0, ix);

  BCHK (this->shape (i1) != this->shape (i2), TensErr, Contraction over diff. sized indices, this->shape (i2), t);
  
  Index bx (0, ix.size());
  while (bx(0) < ix(0)) {
      unsigned long l = t.calc_offs (bx);
      for (unsigned m = 0; m < this->shape (i1); ++m)
        t.data (l) += (*this) (idx_fill_in2 (bx, i1, m, i2, m));
      bx.next_idx (ix);
  }
  return t;
}

#if 1
template <typename T>
Tensor<T> Tensor<T>::cntrmul (const Tensor<T>& t2, const unsigned i1, const unsigned i2) const
{
  BCHK(this->shape (i1) != t2.shape (i2), TensErr, CntrMul: Diff. sizes, t2.shape(i2), 0);
  Index sh2 (this->shape); sh2.append (t2.shape);
  Index ix (idx_remove2 (sh2, i1, i2 + this->rank));
  
  Tensor <T> result (0, ix);
  Index ixr (0, result.rank);
  
  while (ixr(0) < ix(0)) {
      unsigned long i = result.calc_offs (ixr);
      for (unsigned m = 0; m < this->shape (i1); ++m) {
          Index ix1 (idx_fill_in1 (ixr.slice(0, this->rank-1), i1, m));
          Index ix2 (idx_fill_in1 (ixr.slice(this->rank-1, result.rank), i2, m));
          result.data(i) += (*this) (ix1) * t2 (ix2);
      }
      ixr.next_idx (ix);
  }
  return result;
}

#else

template <typename T>
Tensor<T> Tensor<T>::cntrmul (const Tensor<T>& t2, const unsigned i1, const unsigned i2) const
{
   BCHK(this->shape (i1) != t2.shape (i2), TensErr, CrtMul: Diff. sizes, t2.shape(i2), 0);
   Index ix (this->rank + t2.rank - 2); // rank of result tensor
   unsigned c = 0;              // counter
   {
     for (unsigned i = 0; i < this->rank; i++) 
        if (i != i1) 
            ix (c++) = this->shape (i);
   }                            // index from 1rst
   {
     for (unsigned i = 0; i < t2.rank; i++) i
        if (i != i2) 
            ix (c++) = t2.shape (i);
   }                            // & from 2nd tensor excl. sum. index
  
   Tensor <T> result (0, ix);
   Index ixr (result.rank);
   Index ix1 (this->rank);
   Index ix2 (t2.rank);
   for (unsigned long i = 0; i < result.noel; ++i) {
       ixr = result.calc_indx (i);
       for (unsigned m = 0; m < this->shape (i1); ++m) {
           c = 0;                       // reset counter
           for (unsigned k = 0; k < this->rank; ++k) {
               if (k != i1) 
                        ix1 (k) = ixr (c++);
               else
                        ix1 (k) = m;
           }
           for (unsigned l = 0; l < t2.rank; ++l) {
               if (l != i2) 
                        ix2 (l) = ixr (c++);
               else
                        ix2 (l) = m;
           }
           result.data(i) += (*this)(ix1) * t2(ix2);
       }
   }
   return result;
}

#endif

// Same as cntrmul, but friend instead of member
template <typename T>
Tensor<T> ctrmul (const Tensor<T>& t1, const Tensor<T>& t2, const unsigned i1, const unsigned i2)
{
  BCHK(t1.shape (i1) != t2.shape (i2), TensErr, CtrMul: Diff. sizes, t2.shape(i2), 0);
  Index sh2 (t1.shape); sh2.append (t2.shape);
  Index ix (idx_remove2 (sh2, i1, i2 + t1.rank));
  Tensor < T > result (0, ix);

  Index ixr (0, result.rank);
  while (ixr(0) < ix(0)) {
      unsigned long i = result.calc_offs (ixr);
      for (unsigned m = 0; m < t1.shape (i1); ++m) {
          result.data(i) += t1 (idx_fill_in1 (ixr.slice(0,t1.rank-1), i1, m))
            * t2 (idx_fill_in1 (ixr.slice(t1.rank-1,result.rank), i2, m));
      }
      ixr.next_idx (ix);
  }
  return result;
}

template <typename T>
Tensor<T> metrmul (const Tensor<T>& metr, const Tensor<T>& t, const unsigned i1, const unsigned i2)
{
  BCHK(metr.shape (i1) != t.shape (i2), TensErr, MetrMul: Diff. sizes, t.shape(i2), 0);
  BCHK(metr.rank != 2, TensErr, Metric w/ rank != 2 ?, metr.rank, t);
  BCHKNR(metr.shape(0) != metr.shape(1), TensErr, Non-quadr. metric ?, metr.shape(i1));

  Index ix (t.shape);
  Tensor < T > result (0, ix);

  unsigned long io = 1 - i1;            // other index in metric

  Index ixr (0, result.rank);
  Index mix (2);
  while (ixr(0) < ix(0)) {
      unsigned long i = result.calc_offs (ixr);
      for (unsigned m = 0; m < metr.shape (i1); ++m) {
          mix (io) = ixr (i2); mix (i1) = m;
          result.data (i) += metr (mix) * t (idx_set1 (ixr, i2, m));
      }
      ixr.next_idx (ix);
  }
  return result;
}

template <typename T>
Tensor<T> Tensor<T>::drctmul (const Tensor<T> &t1) const
{
  Index ix (t1.shape);          // rank of result tensor
  ix.append (this->shape);
  Tensor < T > result (0, ix);
  Index ixr (0, result.rank);
  while (ixr(0) < ix(0)) {
      result (ixr) = (*this) (ixr.slice (0, this->rank)) 
        * t1 (ixr.slice (this->rank, result.rank));
      ixr.next_idx (ix);
  }
  return result;
}

// Same as drctmul, but friend
template <typename T>
Tensor<T> dctmul (const Tensor<T> &t0, const Tensor<T> &t1)
{
  Index ix (t1.shape);          // rank of result tensor
  ix.append (t0.shape);
  Tensor < T > result (0, ix);
  Index ixr (0, result.rank);
  while (ixr(0) < ix(0)) {
      result (ixr) = t0 (ixr.slice (0, t0.rank)) 
        * t1 (ixr.slice (t0.rank, result.rank));
      ixr.next_idx (ix);
  }
  return result;
}

template <typename T>
Tensor<T>& Tensor<T>::operator += (const Tensor<T>& t)
{
    BCHK(this->shape!=t.shape, TensErr, Diff Tensor shape in op +=, t.noel, *this);
    if (this->layout == t.layout)
      for (unsigned long i = 0; i < this->noel; i++) 
        this->data(i) += t.data(i);
    else {
      Index ix (0, this->rank);
      while (ix(0) < this->shape(0)) {
            (*this)(ix) += t(ix);
            ix.next_idx(this->shape);
      }
    }
    return *this;
}

template <typename T>
Tensor<T>& Tensor<T>::operator -= (const Tensor<T>& t)
{
    BCHK(this->shape!=t.shape, TensErr, Diff Tensor shape in op -=, t.noel, *this);
    if (this->layout == t.layout)
      for (unsigned long i = 0; i < this->noel; i++) 
        this->data(i) -= t.data(i);
    else {
      Index ix (0, this->rank);
      while (ix(0) < this->shape(0)) {
           (*this)(ix) -= t(ix);
           ix.next_idx(this->shape);
      }
    }
    return *this;
}

template <typename T>
Tensor<T> Tensor<T>::operator + (const Tensor<T>& t) const
{
    Tensor<T> r (this->shape);
    BCHK(this->shape!=t.shape, TensErr, Diff Tensor shape in op +, t.noel, r);
    if (this->layout == t.layout)
      for (unsigned long i = 0; i < this->noel; i++) 
                r.data(i) = this->data(i) + t.data(i);
    else {
      Index ix (0, this->rank);
      while (ix(0) < this->shape(0)) {
            r(ix) = (*this)(ix) + t(ix);
            ix.next_idx(this->shape);
      }
    }
    return r;
}

template <typename T>
Tensor<T> Tensor<T>::operator - (const Tensor<T>& t) const
{
    Tensor<T> r (this->shape);
    BCHK(this->shape!=t.shape, TensErr, Diff Tensor shape in op -, t.noel, r);
    if (this->layout == t.layout)
      for (unsigned long i = 0; i < this->noel; i++) 
                r.data(i) = this->data(i) - t.data(i);
    else {
      Index ix (0, this->rank);
      while (ix(0) < this->shape(0)) {
              r(ix) = (*this)(ix) - t(ix);
              ix.next_idx(this->shape);
      }
    }
    return r;
}


template <typename T>
Tensor<T> Tensor<T>::mult (const T m) const
{
        Tensor<T> res (this->shape);
        for (unsigned long i = 0; i < this->noel; ++i)
                res.data(i) = m * this->data(i);
        return res;
}

INST(template <typename T> class Tensor friend Tensor<T> operator * (const T, const Tensor<T> &);)
template <typename T>
inline Tensor<T> operator * (const T m, const Tensor<T> &t)
{ return t.mult (m); } 


NAMESPACE_END

#endif /* TBCI_TENSOR_H */