source: Sophya/trunk/SophyaLib/TArray/tarray.cc@ 2153

//      template array class for numerical types
//                     R. Ansari, C.Magneville   03/2000

#include "machdefs.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "pexceptions.h"
#include "tarray.h"

/*!
  \class SOPHYA::TArray
  \ingroup TArray
  Class for template arrays

  This class implements arrays of dimensions up to
  \ref BASEARRAY_MAXNDIMS "BASEARRAY_MAXNDIMS"
*/

/*! \ingroup TArray
  \typedef sa_size_t
  \brief Array index range and size, defined to be a 4-byte or 8-byte integer 
*/

// -------------------------------------------------------
//                Methodes de la classe
// -------------------------------------------------------

////////////////////////// Les constructeurs / destructeurs

//! Default constructor
template <class T>
TArray<T>::TArray()
  : BaseArray() , mNDBlock() 
{
}

//! Constructor
/*!
  \param ndim : number of dimensions (less or equal to
                \ref BASEARRAY_MAXNDIMS "BASEARRAY_MAXNDIMS")
  \param siz[ndim] : size along each dimension
  \param step : step (same for all dimensions)
 */
template <class T>
TArray<T>::TArray(int_4 ndim, const sa_size_t * siz, sa_size_t step)
  : BaseArray() , mNDBlock(ComputeTotalSize(ndim, siz, step, 1))
{
  string exmsg = "TArray<T>::TArray(int_4, sa_size_t *, sa_size_t)";
  if (!UpdateSizes(ndim, siz, step, 0, exmsg))  throw( ParmError(exmsg) );
}

//! Constructor
/*!
  \param nx,ny,nz,nt,nu : sizes along first, second, third, fourth and fifth dimension
 */
template <class T>
TArray<T>::TArray(sa_size_t nx, sa_size_t ny, sa_size_t nz, sa_size_t nt, sa_size_t nu)
  : BaseArray() , mNDBlock(nx*((ny>0)?ny:1)*((nz>0)?nz:1)*((nt>0)?nt:1)*((nu>0)?nu:1)) 
{
  sa_size_t size[BASEARRAY_MAXNDIMS];
  size[0] = nx;   size[1] = ny;   size[2] = nz;
  size[3] = nt;   size[4] = nu;
  int_4 ndim = 1;
  if ((size[1] > 0) && (size[2] > 0) && (size[3] > 0) && (size[4] > 0) ) ndim = 5;
  else if ((size[1] > 0) && (size[2] > 0) && (size[3] > 0) ) ndim = 4;
  else if ((size[1] > 0) && (size[2] > 0)) ndim = 3;
  else if (size[1] > 0)  ndim = 2;
  else ndim = 1;
  string exmsg = "TArray<T>::TArray(sa_size_t, sa_size_t, sa_size_t, sa_size_t, sa_size_t)";
  if (!UpdateSizes(ndim, size, 1, 0, exmsg))  throw( ParmError(exmsg) );
}

//! Constructor
/*!
  \param ndim : number of dimensions
  \param siz[ndim] : size along each dimension
  \param db : datas are given by this NDataBlock
  \param share : if true, data are shared, if false they are copied
  \param step : step (same for all dimensions) in data block
  \param offset : offset for first element in data block
 */
template <class T>
TArray<T>::TArray(int_4 ndim, const sa_size_t * siz, NDataBlock<T> & db, bool share, sa_size_t step, sa_size_t offset)
  : BaseArray() , mNDBlock(db, share) 
{
  string exmsg = "TArray<T>::TArray(int_4, sa_size_t *,  NDataBlock<T> & ... )";
  if (!UpdateSizes(ndim, siz, step, offset, exmsg))  throw( ParmError(exmsg) );
  if (mNDBlock.Size() < ComputeTotalSize(ndim, siz, step, offset)) {
    exmsg += " DataBlock.Size() < ComputeTotalSize(...) " ;
    throw( ParmError(exmsg) );
  }
}

//! Constructor
/*!
  \param ndim : number of dimensions
  \param siz[ndim] : size along each dimension
  \param values : datas are given by this pointer
  \param share : if true, data are shared, if false they are copied
  \param step : step (same for all dimensions) in data block
  \param offset : offset for first element in data block
  \param br : if not NULL, dats are bridge with other datas
  \sa NDataBlock
 */
template <class T>
TArray<T>::TArray(int_4 ndim, const sa_size_t * siz, T* values, sa_size_t step, sa_size_t offset, Bridge* br)
  : BaseArray() , mNDBlock(ComputeTotalSize(ndim, siz, step, 1), values, br) 
{
  string exmsg = "TArray<T>::TArray(int_4, sa_size_t *, T* ... )";
  if (!UpdateSizes(ndim, siz, step, offset, exmsg))  throw( ParmError(exmsg) );
}

//! Constructor by copy
/*!
  \warning datas are \b SHARED with \b a.
  \sa NDataBlock::NDataBlock(const NDataBlock<T>&)
 */
template <class T>
TArray<T>::TArray(const TArray<T>& a)
  : BaseArray() , mNDBlock(a.mNDBlock) 
{
  string exmsg = "TArray<T>::TArray(const TArray<T>&)";
  if (!UpdateSizes(a, exmsg))  throw( ParmError(exmsg) );
  if (a.mInfo) mInfo = new DVList(*(a.mInfo));
}

//! Constructor by copy
/*!
  \param share : if true, data are shared, if false they are copied
 */
template <class T>
TArray<T>::TArray(const TArray<T>& a, bool share)
  : BaseArray() , mNDBlock(a.mNDBlock, share) 
{
  if (a.NbDimensions() == 0)  return;  
  string exmsg = "TArray<T>::TArray(const TArray<T>&, bool)";
  if (!UpdateSizes(a, exmsg))  throw( ParmError(exmsg) );
  if (a.mInfo) mInfo = new DVList(*(a.mInfo));
}

//! Constructor with size and contents copied (after conversion) from a different type TArray
template <class T>
TArray<T>::TArray(const BaseArray& a)
  : BaseArray() , mNDBlock() 
{
  if (a.NbDimensions() == 0)  return;  
  string exmsg = "TArray<T>::TArray(const BaseArray&)";
  if (!UpdateSizes(a, exmsg))  throw( ParmError(exmsg) );
  mNDBlock.ReSize(totsize_);    
  //  if (a.mInfo) mInfo = new DVList(*(a.mInfo));  - pb protected !
  ConvertAndCopyElt(a);
}

//! Destructor 
template <class T>
TArray<T>::~TArray()
{
}

////////////////////////// Les methodes de copie/share

//! Set array equal to \b a and return *this
/*!
  If the array is already allocated, CopyElt() is called
  for checking that the two arrays have the same size and
  for copying the array element values. For non allocated 
  arrays,  CloneOrShare() is called. The array memory 
  organization is also copied from \b a.
  \warning Datas are copied (cloned) from \b a.
  \sa CopyElt 
  \sa CloneOrShare
  \sa NDataBlock::operator=(const NDataBlock<T>&)
*/
template <class T>
TArray<T>& TArray<T>::Set(const TArray<T>& a)
{
  if (this == &a) return(*this);
  if (a.NbDimensions() < 1) 
    throw RangeCheckError("TArray<T>::Set(a ) - Array a not allocated ! ");
  if (NbDimensions() < 1) CloneOrShare(a);
  else CopyElt(a);
  return(*this);
}

//! Set array elements equal to the \b a array elements, after conversion 
template <class T>
TArray<T>& TArray<T>::SetBA(const BaseArray& a)
{
  if (this == &a) return(*this);
  if (a.NbDimensions() < 1) 
    throw RangeCheckError("TArray<T>::SetBA(a ) - Array a not allocated ! ");
  if (NbDimensions() < 1) {
    string exmsg = "TArray<T>::SetBA(const BaseArray& a)";
    if (!UpdateSizes(a, exmsg))  throw( ParmError(exmsg) );
    mNDBlock.ReSize(totsize_);    
  }
  ConvertAndCopyElt(a);
  return(*this);
}

//! Clone array \b a
template <class T>
void TArray<T>::Clone(const TArray<T>& a)
{
  string exmsg = "TArray<T>::Clone()";
  if (!UpdateSizes(a, exmsg))  throw( ParmError(exmsg) );
  mNDBlock.Clone(a.mNDBlock);
  if (mInfo) {delete mInfo; mInfo = NULL;}
  if (a.mInfo) mInfo = new DVList(*(a.mInfo));
}

//! Clone if \b a is not temporary, share if temporary
/*! \sa NDataBlock::CloneOrShare(const NDataBlock<T>&) */
template <class T>
void TArray<T>::CloneOrShare(const TArray<T>& a)
{
  string exmsg = "TArray<T>::CloneOrShare()";
  if (!UpdateSizes(a, exmsg))  throw( ParmError(exmsg) );
  mNDBlock.CloneOrShare(a.mNDBlock);
  if (mInfo) {delete mInfo; mInfo = NULL;}
  if (a.mInfo) mInfo = new DVList(*(a.mInfo));
}

//! Share data with a
template <class T>
void TArray<T>::Share(const TArray<T>& a)
{
  string exmsg = "TArray<T>::Share()";
  if (!UpdateSizes(a, exmsg))  throw( ParmError(exmsg) );
  mNDBlock.Share(a.mNDBlock);
  if (mInfo) {delete mInfo; mInfo = NULL;}
  if (a.mInfo) mInfo = new DVList(*(a.mInfo));
}


//! Sets or changes the array size 
/*!
  \param ndim : number of dimensions
  \param siz[ndim] : size along each dimension
  \param step : step (same for all dimensions)
 */
template <class T>
void TArray<T>::ReSize(int_4 ndim, sa_size_t * siz, sa_size_t step)
{
  if (arrtype_ != 0) {
    if (ndim != 2) 
       throw( ParmError("TArray<T>::ReSize(ndim!=2,...) for Matrix" ) ); 
    if ((arrtype_ == 2) && (siz[0] > 1) && (siz[1] > 1)) 
       throw( ParmError("TArray<T>::ReSize(,siz[0]>1 && size[1]>1) for Vector" ) ); 
  }
  string exmsg = "TArray<T>::ReSize(int_4 ...)";
  if (!UpdateSizes(ndim, siz, step, 0, exmsg))  throw( ParmError(exmsg) );
  mNDBlock.ReSize(totsize_);    
}

//! Sets or changes the array size.
/*!
  The array size and memory layout are copied from the array \b a.
  \param a : Array used as template for setting the size and memory layout.
 */
template <class T>
void TArray<T>::ReSize(const BaseArray& a)
{
  if (arrtype_ != 0) {
    if (a.NbDimensions() != 2) 
       throw( ParmError("TArray<T>::ReSize(a.NbDimensions()!=2,...) for Matrix" ) ); 
    if ((arrtype_ == 2) && (a.Size(0) > 1) && (a.Size(1) > 1)) 
       throw( ParmError("TArray<T>::ReSize(a.Size(0)>1 && a.Size(1)>1) for Vector" ) ); 
  }
  string exmsg = "TArray<T>::ReSize(const TArray<T>&)";
  if (!UpdateSizes(a, exmsg))  throw( ParmError(exmsg) );
  mNDBlock.ReSize(totsize_);      
}

//! Re-allocate space for array
/*!
  \param ndim : number of dimensions
  \param siz[ndim] : size along each dimension
  \param step : step (same for all dimensions)
  \param force : if true re-allocation is forced, if not it occurs
          only if the required space is greater than the old one.
 */
template <class T>
void TArray<T>::Realloc(int_4 ndim, sa_size_t * siz, sa_size_t step, bool force)
{
  if (arrtype_ != 0) {
    if (ndim != 2) 
       throw( ParmError("TArray<T>::Realloc(ndim!=2,...) for Matrix" ) ); 
    if ((arrtype_ == 2) && (siz[0] > 1) && (siz[1] > 1)) 
       throw( ParmError("TArray<T>::Realloc(,siz[0]>1 && size[1]>1) for Vector" ) ); 
  }
  string exmsg = "TArray<T>::Realloc()";
  if (!UpdateSizes(ndim, siz, step, 0, exmsg))  throw( ParmError(exmsg) );
  mNDBlock.Realloc(totsize_, force);    
}


//! Compact dimensions in one or more is equal to 1.
template <class T>
TArray<T>& TArray<T>::CompactAllDimensions()
{
  CompactAllDim();
  return(*this);
}

//! Compact dimensions if the last one is equal to 1.
template <class T>
TArray<T>& TArray<T>::CompactTrailingDimensions()
{
  CompactTrailingDim();
  return(*this);
}

//! Give value (in \b double) for element at position \b ip..
template <class T>
MuTyV & TArray<T>::ValueAtPosition(sa_size_t ip) const
{
#ifdef SO_BOUNDCHECKING
  if (ip >= totsize_)  throw( ParmError("TArray<T>::ValueAtPosition(sa_size_t ip) Out-of-bound Error") );
#endif
  my_mtv = *(mNDBlock.Begin()+Offset(ip));
  return( my_mtv );
}

//! Return array with elements packed
/*!
  \param force : if true, pack elements in a new array.
                 If false and array is already packed, return
                 an array that share data with the current one.
  \return packed array
 */
template <class T>
TArray<T> TArray<T>::PackElements(bool force) const
{
  if (NbDimensions() < 1) 
    throw RangeCheckError("TArray<T>::PackElements() - Not Allocated Array ! ");
  if ( !force && (AvgStep() == 1) ) {
    TArray<T> ra;
    ra.Share(*this);
    ra.SetTemp(true);
    return(ra);
  }
  else {
    TArray<T> ra(ndim_, size_, 1);
    ra.CopyElt(*this);
    ra.SetTemp(true);
    return(ra);
  }
}

// SubArrays 
// $CHECK$ Reza 03/2000  Doit-on declarer cette methode const ?
//! Extract a sub-array
/*!
  \param rx,ry,rz,rt,ru : range of extraction along dimensions
  \sa Range
 */
template <class T>
TArray<T> TArray<T>::SubArray(Range rx, Range ry, Range rz, Range rt, Range ru) const
{
  if (NbDimensions() < 1) 
    throw RangeCheckError("TArray<T>::operator () (Range, ...) - Not Allocated Array ! ");
  int_4 ndim = 0;
  sa_size_t size[BASEARRAY_MAXNDIMS];
  sa_size_t step[BASEARRAY_MAXNDIMS];
  sa_size_t pos[BASEARRAY_MAXNDIMS];
  size[0] = rx.Size();
  size[1] = ry.Size();
  size[2] = rz.Size();
  size[3] = rt.Size();
  size[4] = ru.Size();

  step[0] = rx.Step();
  step[1] = ry.Step();
  step[2] = rz.Step();
  step[3] = rt.Step();
  step[4] = ru.Step();

  pos[0] = rx.Start();
  pos[1] = ry.Start();
  pos[2] = rz.Start();
  pos[3] = rt.Start();
  pos[4] = ru.Start();

  ndim = ndim_;
  TArray<T> ra;
  UpdateSubArraySizes(ra, ndim, size, pos, step); 
  ra.DataBlock().Share(this->DataBlock());
  ra.SetTemp(true);
  return(ra);
}

//   ...... Operation de calcul sur les tableaux ......
//              ------- Attention --------
//  Boucles normales prenant en compte les steps ....
//  Possibilite de // , vectorisation

//! Fill TArray with Sequence \b seq
/*!
  \param seq : sequence to fill the array
  \sa Sequence
 */
template <class T>
TArray<T>& TArray<T>::SetSeq(Sequence const & seq)
{
  if (NbDimensions() < 1) 
    throw RangeCheckError("TArray<T>::SetSeq(Sequence ) - Not Allocated Array ! ");
  
  T * pe;
  sa_size_t j,k;
  int_4 ka;
  if (arrtype_ == 0)  ka = 0;
  else ka = macoli_;
  sa_size_t step = Step(ka);
  sa_size_t gpas = Size(ka);
  sa_size_t naxa = Size()/Size(ka);
  for(j=0; j<naxa; j++)  {
    pe = mNDBlock.Begin()+Offset(ka,j);
/*  
    Appel explicite de l'operateur de conversion 
    suite a la suggestion de M. Reinecke, Reza 31/7/2002 
#if !defined(__GNUG__)
    for(k=0; k<gpas; k++)  pe[k*step] = (T) seq(j*gpas+k);
#else
    // g++ (up to 2.95.1) se melange les pinceaux  s'il y a le cast (T) pour l'instanciation des complexes 
    for(k=0; k<gpas; k++)  pe[k*step] = seq(j*gpas+k);
#endif
   --- Appel explicite de l'operateur de conversion sur l'objet MuTyV 
*/
    for(k=0; k<gpas; k++)  pe[k*step] = seq(j*gpas+k).operator T();
  }
  return(*this);
}

//  >>>> Operations avec 2nd membre de type scalaire 

//! Fill an array with a constant value \b x
template <class T>
TArray<T>& TArray<T>::SetT(T x)
{
  if (NbDimensions() < 1) 
    throw RangeCheckError("TArray<T>::SetT(T )  - Not Allocated Array ! ");
  T * pe;
  sa_size_t j,k;
  if (AvgStep() > 0)   {  // regularly spaced elements
    sa_size_t step = AvgStep();
    sa_size_t maxx = totsize_*step;
    pe = Data();
    for(k=0; k<maxx; k+=step )  pe[k] = x;
  }
  else {    // Non regular data spacing ...
    int_4 ka = MaxSizeKA();
    sa_size_t step = Step(ka);
    sa_size_t gpas = Size(ka)*step;
    sa_size_t naxa = Size()/Size(ka);
    for(j=0; j<naxa; j++)  {
      pe = mNDBlock.Begin()+Offset(ka,j);
      for(k=0; k<gpas; k+=step)  pe[k] = x;
    }
  }
  return(*this);
}

//! Add a constant value \b x to an array
template <class T>
TArray<T>& TArray<T>::Add(T x)
{
  if (NbDimensions() < 1) 
    throw RangeCheckError("TArray<T>::Add(T )  - Not Allocated Array ! ");
  T * pe;
  sa_size_t j,k;
  if (AvgStep() > 0)   {  // regularly spaced elements
    sa_size_t step = AvgStep();
    sa_size_t maxx = totsize_*step;
    pe = Data();
    for(k=0; k<maxx; k+=step )  pe[k] += x;
  }
  else {    // Non regular data spacing ...
    int_4 ka = MaxSizeKA();
    sa_size_t step = Step(ka);
    sa_size_t gpas = Size(ka)*step;
    sa_size_t naxa = Size()/Size(ka);
    for(j=0; j<naxa; j++)  {
      pe = mNDBlock.Begin()+Offset(ka,j);
      for(k=0; k<gpas; k+=step)  pe[k] += x;
    }
  }
  return(*this);
}

//! Substract a constant value \b x to an array
/*! 
Substract a constant from the *this = *this-x
\param fginv == true : Perfoms the inverse subtraction (*this = x-(*this)) 
*/
template <class T>
TArray<T>& TArray<T>::Sub(T x, bool fginv)
{
  if (NbDimensions() < 1) 
    throw RangeCheckError("TArray<T>::Sub(T )  - Not Allocated Array ! ");
  T * pe;
  sa_size_t j,k;
  if (AvgStep() > 0)   {  // regularly spaced elements
    sa_size_t step = AvgStep();
    sa_size_t maxx = totsize_*step;
    pe = Data();
    if (fginv)
      for(k=0; k<maxx; k+=step )  pe[k] = x-pe[k];
    else
      for(k=0; k<maxx; k+=step )  pe[k] -= x;
  }
  else {    // Non regular data spacing ...
    int_4 ka = MaxSizeKA();
    sa_size_t step = Step(ka);
    sa_size_t gpas = Size(ka)*step;
    sa_size_t naxa = Size()/Size(ka);
    for(j=0; j<naxa; j++)  {
      pe = mNDBlock.Begin()+Offset(ka,j);
      if (fginv)
        for(k=0; k<gpas; k+=step)  pe[k] = x-pe[k];
      else
        for(k=0; k<gpas; k+=step)  pe[k] -= x;
    }
  }
  return(*this);
}

//! Multiply an array by a constant value \b x
template <class T>
TArray<T>& TArray<T>::Mul(T x)
{
  if (NbDimensions() < 1) 
    throw RangeCheckError("TArray<T>::Mul(T )  - Not Allocated Array ! ");
  T * pe;
  sa_size_t j,k;
  if (AvgStep() > 0)   {  // regularly spaced elements
    sa_size_t step = AvgStep();
    sa_size_t maxx = totsize_*step;
    pe = Data();
    for(k=0; k<maxx; k+=step )  pe[k] *= x;
  }
  else {    // Non regular data spacing ...
    int_4 ka = MaxSizeKA();
    sa_size_t step = Step(ka);
    sa_size_t gpas = Size(ka)*step;
    sa_size_t naxa = Size()/Size(ka);
    for(j=0; j<naxa; j++)  {
      pe = mNDBlock.Begin()+Offset(ka,j);
      for(k=0; k<gpas; k+=step)  pe[k] *= x;
    }
  }
  return(*this);
}

//! Divide an array by a constant value \b x
/*! 
Divide the array by a constant *this = *this/x
\param fginv == true : Perfoms the inverse division (*this = x/(*this)) 
*/
template <class T>
TArray<T>& TArray<T>::Div(T x, bool fginv)
{
  if (NbDimensions() < 1) 
    throw RangeCheckError("TArray<T>::Div(T )  - Not Allocated Array ! ");
  if (!fginv && (x == (T) 0) ) 
    throw MathExc("TArray<T>::Div(T )  - Divide by zero ! ");
  T * pe;
  sa_size_t j,k;
  if (AvgStep() > 0)   {  // regularly spaced elements
    sa_size_t step = AvgStep();
    sa_size_t maxx = totsize_*step;
    pe = Data();
    if (fginv)
      for(k=0; k<maxx; k+=step )  pe[k] = x/pe[k];
    else 
      for(k=0; k<maxx; k+=step )  pe[k] /= x;
  }
  else {    // Non regular data spacing ...
    int_4 ka = MaxSizeKA();
    sa_size_t step = Step(ka);
    sa_size_t gpas = Size(ka)*step;
    sa_size_t naxa = Size()/Size(ka);
    for(j=0; j<naxa; j++)  {
      pe = mNDBlock.Begin()+Offset(ka,j);
      if (fginv)
        for(k=0; k<gpas; k+=step)  pe[k] = x/pe[k];
      else 
        for(k=0; k<gpas; k+=step)  pe[k] /= x;
    }
  }
  return(*this);
}


//! Replace array elements values by their opposite ( a(i) -> -a(i) )
template <class T>
TArray<T>& TArray<T>::NegateElt()
{
  if (NbDimensions() < 1) 
    throw RangeCheckError("TArray<T>::NegateElt()  - Not Allocated Array ! ");
  T * pe;
  sa_size_t j,k;
  if (AvgStep() > 0)   {  // regularly spaced elements
    sa_size_t step = AvgStep();
    sa_size_t maxx = totsize_*step;
    pe = Data();
    for(k=0; k<maxx; k+=step )  pe[k] = -pe[k];
  }
  else {    // Non regular data spacing ...
    int_4 ka = MaxSizeKA();
    sa_size_t step = Step(ka);
    sa_size_t gpas = Size(ka)*step;
    sa_size_t naxa = Size()/Size(ka);
    for(j=0; j<naxa; j++)  {
      pe = mNDBlock.Begin()+Offset(ka,j);
      for(k=0; k<gpas; k+=step)  pe[k] = -pe[k];
    }
  }
  return(*this);
}

//  >>>> Operations avec 2nd membre de type tableau
//! Add two TArrays
template <class T>
TArray<T>& TArray<T>::AddElt(const TArray<T>& a)
{
  if (NbDimensions() < 1) 
    throw RangeCheckError("TArray<T>::AddElt(const TArray<T>& )  - Not Allocated Array ! ");
  bool smo;
  if (!CompareSizes(a, smo)) 
    throw(SzMismatchError("TArray<T>::AddElt(const TArray<T>&) SizeMismatch")) ;

  T * pe;
  const T * pea;
  sa_size_t j,k,ka;
  if (smo && (AvgStep() > 0) && (a.AvgStep() > 0))   {  // regularly spaced elements
    sa_size_t step = AvgStep();
    sa_size_t stepa = a.AvgStep();
    sa_size_t maxx = totsize_*step;
    pe = Data();
    pea = a.Data();
    for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] += pea[ka] ;
  }
  else {    // Non regular data spacing ...
    int_4 ax,axa;
    sa_size_t step, stepa;
    sa_size_t gpas, naxa;
    GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa);
    for(j=0; j<naxa; j++)  {
      pe = mNDBlock.Begin()+Offset(ax,j);
      pea = a.DataBlock().Begin()+a.Offset(axa,j);
      for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] += pea[ka]; 
    }
  }
  return(*this);
}

//! Substract two TArrays
/*! 
Substract two TArrays *this = *this-a
\param fginv == true : Perfoms the inverse subtraction (*this = a-(*this)) 
*/
template <class T>
TArray<T>& TArray<T>::SubElt(const TArray<T>& a, bool fginv)
{
  if (NbDimensions() < 1) 
    throw RangeCheckError("TArray<T>::SubElt(const TArray<T>& )  - Not Allocated Array ! ");
  bool smo;
  if (!CompareSizes(a, smo)) 
    throw(SzMismatchError("TArray<T>::SubElt(const TArray<T>&) SizeMismatch")) ;

  T * pe;
  const T * pea;
  sa_size_t j,k,ka;
  if (smo && (AvgStep() > 0) && (a.AvgStep() > 0) )   {  // regularly spaced elements
    sa_size_t step = AvgStep();
    sa_size_t stepa = a.AvgStep();
    sa_size_t maxx = totsize_*step;
    pe = Data();
    pea = a.Data();
    if (fginv) 
      for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] = pea[ka]-pe[k] ;
    else 
      for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] -= pea[ka] ;
  }
  else {    // Non regular data spacing ...
    int_4 ax,axa;
    sa_size_t step, stepa;
    sa_size_t gpas, naxa;
    GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa);
    for(j=0; j<naxa; j++)  {
      pe = mNDBlock.Begin()+Offset(ax,j);
      pea = a.DataBlock().Begin()+a.Offset(axa,j);
      if (fginv) 
        for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] = pea[ka]-pe[k] ;
      else
        for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] -= pea[ka]; 
    }
  }
  return(*this);
}


//! Multiply two TArrays (elements by elements)
template <class T>
TArray<T>& TArray<T>::MulElt(const TArray<T>& a)
{
  if (NbDimensions() < 1) 
    throw RangeCheckError("TArray<T>::MulElt(const TArray<T>& )  - Not Allocated Array ! ");
  bool smo;
  if (!CompareSizes(a, smo)) 
    throw(SzMismatchError("TArray<T>::MulElt(const TArray<T>&) SizeMismatch")) ;

  T * pe;
  const T * pea;
  sa_size_t j,k,ka;
  if (smo && (AvgStep() > 0) && (a.AvgStep() > 0) )   {  // regularly spaced elements
    sa_size_t step = AvgStep();
    sa_size_t stepa = a.AvgStep();
    sa_size_t maxx = totsize_*step;
    pe = Data();
    pea = a.Data();
    for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] *= pea[ka] ;
  }
  else {    // Non regular data spacing ...
    int_4 ax,axa;
    sa_size_t step, stepa;
    sa_size_t gpas, naxa;
    GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa);
    for(j=0; j<naxa; j++)  {
      pe = mNDBlock.Begin()+Offset(axa,j);
      pea = a.DataBlock().Begin()+a.Offset(axa,j);
      for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] *= pea[ka]; 
    }
  }
  return(*this);
}


//! Divide two TArrays (elements by elements)
/*! 
Divide two TArrays *this = (*this)/a
\param fginv == true : Perfoms the inverse division (*this = a/(*this)) 
\param divzero == true : if a(i)==0. result is set to zero: (*this)(i)==0.
*/
template <class T>
TArray<T>& TArray<T>::DivElt(const TArray<T>& a, bool fginv, bool divzero)
{
  if (NbDimensions() < 1) 
    throw RangeCheckError("TArray<T>::DivElt(const TArray<T>& )  - Not Allocated Array ! ");
  bool smo;
  if (!CompareSizes(a, smo)) 
    throw(SzMismatchError("TArray<T>::DivElt(const TArray<T>&) SizeMismatch")) ;

  T * pe;
  const T * pea;
  sa_size_t j,k,ka;
  if (smo && (AvgStep() > 0) && (a.AvgStep() > 0) )   {  // regularly spaced elements
    sa_size_t step = AvgStep();
    sa_size_t stepa = a.AvgStep();
    sa_size_t maxx = totsize_*step;
    pe = Data();
    pea = a.Data();
    if(divzero) {
      if (fginv) 
        for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )
          {if(pe[k]==(T)0) pe[k] = (T)0; else pe[k] = pea[ka]/pe[k];}
      else 
        for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )
          {if(pea[k]==(T)0) pe[k] = (T)0; else pe[k] /= pea[ka] ;}
    } else {
      if (fginv) 
        for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] = pea[ka]/pe[k];
      else 
        for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] /= pea[ka] ;
    }
  }
  else {    // Non regular data spacing ...
    int_4 ax,axa;
    sa_size_t step, stepa;
    sa_size_t gpas, naxa;
    GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa);
    for(j=0; j<naxa; j++)  {
      pe = mNDBlock.Begin()+Offset(ax,j);
      pea = a.DataBlock().Begin()+a.Offset(axa,j);
      if(divzero) {
        if (fginv) 
          for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)
            {if(pe[k]==(T)0) pe[k] = (T)0; else pe[k] = pea[ka]/pe[k];} 
        else
          for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)
            {if(pea[k]==(T)0) pe[k] = (T)0; else pe[k] /= pea[ka];} 
      } else {
        if (fginv) 
          for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] = pea[ka]/pe[k]; 
        else
          for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] /= pea[ka]; 
      }
    }
  }
  return(*this);
}

//! Copy elements of \b a
template <class T>
TArray<T>& TArray<T>::CopyElt(const TArray<T>& a)
{
  if (NbDimensions() < 1) 
    throw RangeCheckError("TArray<T>::CopyElt(const TArray<T>& )  - Not Allocated Array ! ");
  bool smo;
  if (!CompareSizes(a, smo)) 
    throw(SzMismatchError("TArray<T>::CopyElt(const TArray<T>&) SizeMismatch")) ;

  T * pe;
  const T * pea;
  sa_size_t j,k,ka;
  if (smo && (AvgStep() > 0) && (a.AvgStep() > 0) )   {  // regularly spaced elements
    sa_size_t step = AvgStep();
    sa_size_t stepa = a.AvgStep();
    sa_size_t maxx = totsize_*step;
    pe = Data();
    pea = a.Data();
    for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] = pea[ka] ;
  }
  else {    // Non regular data spacing ...
    int_4 ax,axa;
    sa_size_t step, stepa;
    sa_size_t gpas, naxa;
    GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa);
    for(j=0; j<naxa; j++)  {
      pe = mNDBlock.Begin()+Offset(ax,j);
      pea = a.DataBlock().Begin()+a.Offset(axa,j);
      for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] = pea[ka]; 
    }
  }
  return(*this);
}

//! Converts and Copy elements of \b a
template <class T>
TArray<T>& TArray<T>::ConvertAndCopyElt(const BaseArray& a)
{
  if (NbDimensions() < 1) 
    throw RangeCheckError("TArray<T>::ConvertAndCopyElt(const TArray<T>& )  - Not Allocated Array ! ");
  bool smo;
  if (!CompareSizes(a, smo)) 
    throw(SzMismatchError("TArray<T>::ConvertAndCopyElt(const TArray<T>&) SizeMismatch")) ;

  T * pe;
  sa_size_t j,k,ka;
  sa_size_t offa;
  // Non regular data spacing ...
  int_4 ax,axa;
  sa_size_t step, stepa;
  sa_size_t gpas, naxa;
  GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa);
  for(j=0; j<naxa; j++)  {
    pe = mNDBlock.Begin()+Offset(ax,j);
    offa = a.Offset(axa,j);
/*  
    Appel explicite de l'operateur de conversion 
    suite a la suggestion de M. Reinecke, Reza 31/7/2002 
#if !defined(__GNUG__)
    for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] = (T)a.ValueAtPosition(offa+ka); 
#else
    // g++ (up to 2.95.1) se melange les pinceaux  s'il y a le cast (T) pour l'instanciation des complexes 
    for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] = a.ValueAtPosition(offa+ka); 
#endif
   --- Appel explicite de l'operateur de conversion sur l'objet MuTyV 
*/
    for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  
      pe[k] = a.ValueAtPosition(offa+ka).operator T(); 
  }
  return(*this);
}


// Somme et produit des elements
//! Sum all elements
template <class T>
T TArray<T>::Sum() const
{
  if (NbDimensions() < 1) 
    throw RangeCheckError("TArray<T>::Sum()  - Not Allocated Array ! ");
  T ret=0;
  const T * pe;
  sa_size_t j,k;
  if (AvgStep() > 0)   {  // regularly spaced elements
    sa_size_t step = AvgStep();
    sa_size_t maxx = totsize_*step;
    pe = Data();
    for(k=0; k<maxx; k+=step )  ret += pe[k];
  }
  else {    // Non regular data spacing ...
    int_4 ka = MaxSizeKA();
    sa_size_t step = Step(ka);
    sa_size_t gpas = Size(ka)*step;
    sa_size_t naxa = Size()/Size(ka);
    for(j=0; j<naxa; j++)  {
      pe = mNDBlock.Begin()+Offset(ka,j);
      for(k=0; k<gpas; k+=step)  ret += pe[k] ;
    }
  }
  return ret;
}

//! Multiply all elements
template <class T>
T TArray<T>::Product() const
{
  if (NbDimensions() < 1) 
    throw RangeCheckError("TArray<T>::Product()  - Not Allocated Array ! ");
  T ret=(T)1;
  const T * pe;
  sa_size_t j,k;
  if (AvgStep() > 0)   {  // regularly spaced elements
    sa_size_t step = AvgStep();
    sa_size_t maxx = totsize_*step;
    pe = Data();
    for(k=0; k<maxx; k+=step )  ret *= pe[k];
  }
  else {    // Non regular data spacing ...
    int_4 ka = MaxSizeKA();
    sa_size_t step = Step(ka);
    sa_size_t gpas = Size(ka)*step;
    sa_size_t naxa = Size()/Size(ka);
    for(j=0; j<naxa; j++)  {
      pe = mNDBlock.Begin()+Offset(ka,j);
      for(k=0; k<gpas; k+=step)  ret *= pe[k] ;
    }
  }
  return ret;
}

//! Returns the sum of all elements squared (Sum
template <class T>
T TArray<T>::SumX2() const
{
  if (NbDimensions() < 1) 
    throw RangeCheckError("TArray<T>::SumX2()  - Not Allocated Array ! ");
  T ret=0;
  const T * pe;
  sa_size_t j,k;
  if (AvgStep() > 0)   {  // regularly spaced elements
    sa_size_t step = AvgStep();
    sa_size_t maxx = totsize_*step;
    pe = Data();
    for(k=0; k<maxx; k+=step )  ret += pe[k]*pe[k];
  }
  else {    // Non regular data spacing ...
    int_4 ka = MaxSizeKA();
    sa_size_t step = Step(ka);
    sa_size_t gpas = Size(ka)*step;
    sa_size_t naxa = Size()/Size(ka);
    for(j=0; j<naxa; j++)  {
      pe = mNDBlock.Begin()+Offset(ka,j);
      for(k=0; k<gpas; k+=step)  ret += pe[k]*pe[k] ;
    }
  }
  return ret;
}

//! Return the minimum and the maximum values of the array elements 
/*!
  This method generates an exception (\c MathExc) if called for complex arrays
*/
template <class T>
void TArray<T>::MinMax(T& min, T& max) const
{
  const T * pe;
  sa_size_t j,k;
  int_4 ka = MaxSizeKA();
  sa_size_t step = Step(ka);
  sa_size_t gpas = Size(ka)*step;
  sa_size_t naxa = Size()/Size(ka);
  min = (*this)[0];
  max = (*this)[0];
  for(j=0; j<naxa; j++)  {
    pe = mNDBlock.Begin()+Offset(ka,j);
    for(k=0; k<gpas; k+=step) {
      if (pe[k]<min)  min = pe[k];
      else if (pe[k]>max)  max = pe[k];      
    }
  }
  return;
}

void TArray< complex<r_4> >::MinMax(complex<r_4>& min, complex<r_4>& max) const
{
  throw MathExc("TArray< complex<r_4> >::MinMax(...) - No order in complex");
}

void TArray< complex<r_8> >::MinMax(complex<r_8>& min, complex<r_8>& max) const
{
  throw MathExc("TArray< complex<r_4> >::MinMax(...) - No order in complex");
}

// ----------------------------------------------------
//       Impression, etc ...
// ----------------------------------------------------

//! Return a string that contain the type \b T of the array
template <class T>
string TArray<T>::InfoString() const
{
  string rs = "TArray<" ;
  rs += typeid(T).name();
  rs += "> ";
  return(rs);
}

//! Print array
/*!
  \param os : output stream
  \param maxprt : maximum numer of print
  \param si : if true,  display attached DvList
  \param ascd : if true, suppresses the display of line numbers,
  suitable for ascii dump format.
  \sa SetMaxPrint
  \sa WriteASCII
 */
template <class T>
void TArray<T>::Print(ostream& os, sa_size_t maxprt, bool si, bool ascd) const
{
  if (maxprt < 0)  maxprt = max_nprt_;
  sa_size_t npr = 0;
  Show(os, si);
  if (ndim_ < 1) return;
  sa_size_t k0,k1,k2,k3,k4;
  for(k4=0; k4<size_[4]; k4++) {
    if ((size_[4] > 1) && ascd) 
      cout << "\n ----- Dimension 5 (U) K4= " << k4 << endl;
    for(k3=0; k3<size_[3]; k3++) {
      if ((size_[3] > 1) && ascd) 
        cout << "\n ----- Dimension 4 (T) K3= " << k3 << endl;
      for(k2=0; k2<size_[2]; k2++) {
        if ((size_[2] > 1) & ascd) 
          cout << "\n ----- Dimension 3 (Z) K2= " << k2 << endl;
        for(k1=0; k1<size_[1]; k1++) {
          if ( (size_[1] > 1) && (size_[0] > 10) && ascd) 
            cout << "----- Dimension 2 (Y) K1= " << k1 << endl;
          for(k0=0; k0<size_[0]; k0++) {
            if(k0 > 0) os << " ";  
            os << Elem(k0, k1, k2, k3, k4);     npr++;
            if (npr >= (sa_size_t) maxprt) {
              if (npr < totsize_)  os << "\n     .... " << endl; return;
            }
          }
          os << endl;
        }
      }
    }
  }
  os <<  endl;
}

//! Fill the array, decoding the ASCII input stream
/*!
  \param is : input stream (ASCII)
  \param nr : Number of non empty (or comment) lines in stream (return value)
  \param nc : Number of columns (= ntot/nlines) (return value)
  \return Number of decoded elements
 */
template <class T>
sa_size_t TArray<T>::ReadASCII(istream& is, sa_size_t & nr, sa_size_t & nc)
{
  EnumeratedSequence es;
  sa_size_t n = es.FillFromFile(is, nr, nc);
  if ( (n < 1) || (nr < 1) || (nc < 1) ) return(n); 
  if (!IsAllocated()) {
    sa_size_t sz[2];
    if (arrtype_ == 2) { // C'est un vecteur 
      sz[0] = sz[1] = 1;
      sz[veceli_] = n;
    }
    else {
      sz[RowsKA()] = nr;
      sz[ColsKA()] = nc;
    }
    ReSize(2, sz);
  }
  SetSeq(es);
  cout << "TArray<T>::ReadASCII()/Info: " << n << " elements read from stream " 
       << " (Row,Col= " << nr << "," << nc << ")" << endl;
  return(n);
}

//! Writes the array content to the output stream, (in ASCII) 
/*!
  \param os : output stream (ASCII)
  \sa Print
 */
template <class T>
void TArray<T>::WriteASCII(ostream& os) const
{
  Print(os, Size(), false, true);
}



///////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////
#ifdef __CXX_PRAGMA_TEMPLATES__
/*
#pragma define_template TArray<uint_1>
#pragma define_template TArray<int_2>
#pragma define_template TArray<uint_4>
*/
#pragma define_template TArray<uint_2>
#pragma define_template TArray<uint_8>
#pragma define_template TArray<int_4>
#pragma define_template TArray<int_8>
#pragma define_template TArray<r_4>
#pragma define_template TArray<r_8>
#pragma define_template TArray< complex<r_4> > 
#pragma define_template TArray< complex<r_8> > 
#endif

#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
/*
template class TArray<uint_1>; 
template class TArray<int_2>;  
template class TArray<uint_4>;
*/
template class TArray<uint_2>;
template class TArray<uint_8>;
template class TArray<int_4>;
template class TArray<int_8>;
template class TArray<r_4>;
template class TArray<r_8>;
template class TArray< complex<r_4> >;
template class TArray< complex<r_8> >;
#endif