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

//      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"
*/

// -------------------------------------------------------
//                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) );
  
}

//! 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) 
{
  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() 
{
  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
/*!
  \warning Datas are copied (cloned) from \b a.
  \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));
}


//! Resize array
/*!
  \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()";
  if (!UpdateSizes(ndim, siz, step, 0, 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.ReSize(totsize_);    
}


//! 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);
#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
  }
  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);
#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
  }
  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
  \sa SetMaxPrint
 */
template <class T>
void TArray<T>::Print(ostream& os, int_4 maxprt, bool si) 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) cout << "\n ----- Dimension 5 (U) K4= " << k4 << endl;
    for(k3=0; k3<size_[3]; k3++) {
      if (size_[3] > 1) cout << "\n ----- Dimension 4 (T) K3= " << k3 << endl;
      for(k2=0; k2<size_[2]; k2++) {
        if (size_[2] > 1) cout << "\n ----- Dimension 3 (Z) K2= " << k2 << endl;
        for(k1=0; k1<size_[1]; k1++) {
          if ( (size_[1] > 1) && (size_[0] > 10) ) 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;
}



///////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////
#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_8>
*/
#pragma define_template TArray<uint_2>
#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_8>;
*/
template class TArray<uint_2>;
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