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

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

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



// Variables statiques globales
template <class T>
char * TArray<T>::ck_op_msg_[6] = {"???", "Size(int )", "IsPacked(int )", 
                                 "Stride(int )", "ElemCheckBound()", "operator()" };
template <class T>
uint_4 TArray<T>::max_nprt_ = 50;

// Methodes statiques globales
template <class T>
void TArray<T>::SetMaxPrint(uint_4 nprt)
{
  max_nprt_ = nprt;
}


template <class T>
uint_8 TArray<T>::ComputeTotalSize(uint_4 ndim, uint_4 * siz, uint_4 step, uint_8 offset)
{
  uint_8 rs = step;
  for(int k=0; k<ndim; k++) rs *= siz[k];
  return(rs+offset);
}

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

// Les constructeurs 
template <class T>
TArray<T>::TArray()
  : mNDBlock() , mInfo(NULL)
// Default constructor
{
  ndim_ = 0;
  for(int k=0; k<TARRAY_MAXNDIMS; k++) step_[k] = size_[k] = 0;
  totsize_ = 0;
  minstep_ = 0;
  moystep_ = 0;
  offset_ = 0;
  // Default for matrices : Lines = Y , Columns = X
  macoli_ = 0;
  marowi_ = 1;
}

template <class T>
TArray<T>::TArray(uint_4 ndim, uint_4 * siz, uint_4 step)
  : mNDBlock(ComputeTotalSize(ndim, siz, step, 1)) , mInfo(NULL)
{
  string exmsg = "TArray<T>::TArray(uint_4, uint_4 *, uint_4)";
  if (!UpdateSizes(ndim, siz, step, 0, exmsg))  throw( ParmError(exmsg) );
}

template <class T>
TArray<T>::TArray(uint_4 nx, uint_4 ny=1, uint_4 nz=1)
  : mNDBlock(nx*((ny>0)?ny:1)*((nz>0)?nz:1)) , mInfo(NULL)
{
  uint_4 size[TARRAY_MAXNDIMS];
  size[0] = nx;   size[1] = ny;   size[2] = nz;
  int ndim = 1;
  if ((size[1] > 0) && (size[2] > 0)) ndim = 3;
  else if (size[1] > 0)  ndim = 2;
  else ndim = 1;
  string exmsg = "TArray<T>::TArray(uint_4, uint_4, uint_4)";
  if (!UpdateSizes(ndim, size, 1, 0, exmsg))  throw( ParmError(exmsg) );
}

template <class T>
TArray<T>::TArray(uint_4 ndim, uint_4 * siz, NDataBlock<T> & db, bool share, uint_4 step, uint_8 offset)
  : mNDBlock(db, share) , mInfo(NULL)
{
  string exmsg = "TArray<T>::TArray(uint_4, uint_4 *,  NDataBlock<T> & ... )";
  if (!UpdateSizes(ndim, siz, step, offset, exmsg))  throw( ParmError(exmsg) );
  
}

template <class T>
TArray<T>::TArray(uint_4 ndim, uint_4 * siz, T* values, uint_4 step, uint_8 offset, Bridge* br)
  : mNDBlock(ComputeTotalSize(ndim, siz, step, 1), values, br) , mInfo(NULL)
{
  string exmsg = "TArray<T>::TArray(uint_4, uint_4 *, T* ... )";
  if (!UpdateSizes(ndim, siz, step, offset, exmsg))  throw( ParmError(exmsg) );
}

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

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

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

template <class T>
TArray<T>& TArray<T>::operator = (const TArray<T>& a)
{
  if (this != &a) { 
    CloneOrShare(a);
    if (mInfo) delete mInfo;
    if (a.mInfo) mInfo = new DVList(*(a.mInfo));
  }
  return(*this);
}

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;
  if (a.mInfo) mInfo = new DVList(*(a.mInfo));
}

template <class T>
void TArray<T>::ReSize(uint_4 ndim, uint_4 * siz, uint_4 step)
{
  string exmsg = "TArray<T>::ReSize()";
  if (!UpdateSizes(ndim, siz, step, 0, exmsg))  throw( ParmError(exmsg) );
  mNDBlock.ReSize(totsize_);    
}

template <class T>
void TArray<T>::Realloc(uint_4 ndim, uint_4 * siz, uint_4 step, bool force)
{
  string exmsg = "TArray<T>::Realloc()";
  if (!UpdateSizes(ndim, siz, step, 0, exmsg))  throw( ParmError(exmsg) );
  mNDBlock.ReSize(totsize_);    
}

template <class T>
bool TArray<T>::CompareSizes(const TArray<T>& a)
{
  if (ndim_ != a.ndim_)  return(false);
  for(int k=0; k<ndim_; k++) 
    if (size_[k] != a.size_[k])  return(false);
  if ( (macoli_ != a.macoli_) || (marowi_ != a.marowi_) )  return(false);
  return(true);
}

template <class T>
TArray<T>& TArray<T>::CompactDim()
{
  if (ndim_ < 2)  return(*this);
  uint_4 ndim = 0;
  uint_4 size[TARRAY_MAXNDIMS];
  uint_4 step[TARRAY_MAXNDIMS];
  for(int k=0; k<ndim_; k++) {
    if (size_[k] < 2)  continue;
    size[ndim] = size_[k];
    step[ndim] = step_[k];
    ndim++;
  }
  if (ndim == 0)  {
    size[0] = size_[0];
    step[0] = step_[0];
    ndim = 1;
  }
  string exmsg = "TArray<T>::CompactDim() ";
  if (!UpdateSizes(ndim, size, step, offset_, exmsg))  throw( ParmError(exmsg) );
  return(*this);
}


template <class T>
uint_4 TArray<T>::MinStepKA() const
{
  for(int ka=0; ka<ndim_; ka++)
    if (step_[ka] == minstep_) return(ka);
  return(0);
}

template <class T>
uint_4 TArray<T>::MaxSizeKA() const
{
  uint_4 ka = 0;
  uint_4 mx = size_[0];
  for(int k=0; k<ndim_; k++)  
    if (size_[k] > mx) {  ka = k;  mx = size_[k];  }
  return(ka);
}


//  Acces lineaire aux elements ....  Calcul d'offset

template <class T>
uint_8 TArray<T>::Offset(uint_8 ip) const
{
if (ip == 0)  return(offset_);
uint_4 idx[TARRAY_MAXNDIMS];
int k;
uint_8 rest = ip;
for(k=0; k<ndim_; k++)   { idx[k] = rest%size_[k];   rest /= size_[k];  }   
uint_8 off = offset_;
for(k=0; k<ndim_; k++)  off += idx[k]*step_[k];
return (off);
}

// SubArrays 
template <class T>
TArray<T> TArray<T>::operator () (Range rx, Range ry, Range rz, Range rt, Range ru)
{
  uint_4 ndim = 0;
  uint_4 size[TARRAY_MAXNDIMS];
  uint_4 step[TARRAY_MAXNDIMS];
  uint_4 pos[TARRAY_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;
  SubArray(ra, ndim, size, pos, step); 
  ra.SetTemp(true);
  return(ra);
}

//   ...... Operation de calcul sur les tableaux ......
//              ------- Attention --------
//  Boucles normales prenant en compte les steps ....
//  Possibilite de // , vectorisation 
template <class T>
TArray<T>& TArray<T>::operator = (Sequence seq)
{
  T * pe;
  uint_8 j,k;
  if (AvgStep() > 0)   {  // regularly spaced elements
    uint_8 step = AvgStep();
    pe = Data();
    for(k=0; k<totsize_; k++ )  pe[k*step] = seq(k);
  }
  else {    // Non regular data spacing ...
    uint_4 ka = MaxSizeKA();
    uint_8 step = Step(ka);
    uint_8 gpas = Size(ka);
    for(j=0; j<totsize_; j += Size(ka))  {
      pe = mNDBlock.Begin()+Offset(j);
      for(k=0; k<gpas; k++)  pe[k*step] = seq(j+k);
    }
  }
  return(*this);
}

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

template <class T>
TArray<T>& TArray<T>::operator = (T x)
{
  T * pe;
  uint_8 j,k;
  if (AvgStep() > 0)   {  // regularly spaced elements
    uint_8 step = AvgStep();
    uint_8 maxx = totsize_*step;
    pe = Data();
    for(k=0; k<maxx; k+=step )  pe[k] = x;
  }
  else {    // Non regular data spacing ...
    uint_4 ka = MaxSizeKA();
    uint_8 step = Step(ka);
    uint_8 gpas = Size(ka)*step;
    for(j=0; j<totsize_; j += Size(ka))  {
      pe = mNDBlock.Begin()+Offset(j);
      for(k=0; k<gpas; k+=step)  pe[k] = x;
    }
  }
  return(*this);
}

template <class T>
TArray<T>& TArray<T>::operator += (T x)
{
  T * pe;
  uint_8 j,k;
  if (AvgStep() > 0)   {  // regularly spaced elements
    uint_8 step = AvgStep();
    uint_8 maxx = totsize_*step;
    pe = Data();
    for(k=0; k<maxx; k+=step )  pe[k] += x;
  }
  else {    // Non regular data spacing ...
    uint_4 ka = MaxSizeKA();
    uint_8 step = Step(ka);
    uint_8 gpas = Size(ka)*step;
    for(j=0; j<totsize_; j += Size(ka))  {
      pe = mNDBlock.Begin()+Offset(j);
      for(k=0; k<gpas; k+=step)  pe[k] += x;
    }
  }
  return(*this);
}

template <class T>
TArray<T>& TArray<T>::operator -= (T x)
{
  T * pe;
  uint_8 j,k;
  if (AvgStep() > 0)   {  // regularly spaced elements
    uint_8 step = AvgStep();
    uint_8 maxx = totsize_*step;
    pe = Data();
    for(k=0; k<maxx; k+=step )  pe[k] -= x;
  }
  else {    // Non regular data spacing ...
    uint_4 ka = MaxSizeKA();
    uint_8 step = Step(ka);
    uint_8 gpas = Size(ka)*step;
    for(j=0; j<totsize_; j += Size(ka))  {
      pe = mNDBlock.Begin()+Offset(j);
      for(k=0; k<gpas; k+=step)  pe[k] -= x;
    }
  }
  return(*this);
}

template <class T>
TArray<T>& TArray<T>::operator *= (T x)
{
  T * pe;
  uint_8 j,k;
  if (AvgStep() > 0)   {  // regularly spaced elements
    uint_8 step = AvgStep();
    uint_8 maxx = totsize_*step;
    pe = Data();
    for(k=0; k<maxx; k+=step )  pe[k] *= x;
  }
  else {    // Non regular data spacing ...
    uint_4 ka = MaxSizeKA();
    uint_8 step = Step(ka);
    uint_8 gpas = Size(ka)*step;
    for(j=0; j<totsize_; j += Size(ka))  {
      pe = mNDBlock.Begin()+Offset(j);
      for(k=0; k<gpas; k+=step)  pe[k] *= x;
    }
  }
  return(*this);
}

template <class T>
TArray<T>& TArray<T>::operator /= (T x)
{
  T * pe;
  uint_8 j,k;
  if (AvgStep() > 0)   {  // regularly spaced elements
    uint_8 step = AvgStep();
    uint_8 maxx = totsize_*step;
    pe = Data();
    for(k=0; k<maxx; k+=step )  pe[k] /= x;
  }
  else {    // Non regular data spacing ...
    uint_4 ka = MaxSizeKA();
    uint_8 step = Step(ka);
    uint_8 gpas = Size(ka)*step;
    for(j=0; j<totsize_; j += Size(ka))  {
      pe = mNDBlock.Begin()+Offset(j);
      for(k=0; k<gpas; k+=step)  pe[k] /= x;
    }
  }
  return(*this);
}


//  >>>> Operations avec 2nd membre de type tableau
template <class T>
TArray<T>& TArray<T>::operator += (const TArray<T>& a)
{
  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::operator += (const TArray<T>&)")) ;

  T * pe;
  const T * pea;
  uint_8 j,k,ka;
  if ((AvgStep() > 0) && (a.AvgStep() > 0) )   {  // regularly spaced elements
    uint_8 step = AvgStep();
    uint_8 stepa = a.AvgStep();
    uint_8 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 ...
    uint_4 ax = MaxSizeKA();
    uint_8 step = Step(ax);
    uint_8 stepa = a.Step(ax);
    uint_8 gpas = Size(ax)*step;
    for(j=0; j<totsize_; j += Size(ax))  {
      pe = mNDBlock.Begin()+Offset(j);
      pea = a.DataBlock().Begin()+a.Offset(j);
      for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] += pea[ka]; 
    }
  }
  return(*this);
}

template <class T>
TArray<T>& TArray<T>::operator -= (const TArray<T>& a)
{
  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::operator -= (const TArray<T>&)")) ;

  T * pe;
  const T * pea;
  uint_8 j,k,ka;
  if ((AvgStep() > 0) && (a.AvgStep() > 0) )   {  // regularly spaced elements
    uint_8 step = AvgStep();
    uint_8 stepa = a.AvgStep();
    uint_8 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 ...
    uint_4 ax = MaxSizeKA();
    uint_8 step = Step(ax);
    uint_8 stepa = a.Step(ax);
    uint_8 gpas = Size(ax)*step;
    for(j=0; j<totsize_; j += Size(ax))  {
      pe = mNDBlock.Begin()+Offset(j);
      pea = a.DataBlock().Begin()+a.Offset(j);
      for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] -= pea[ka]; 
    }
  }
  return(*this);
}


template <class T>
TArray<T>& TArray<T>::MultElt(const TArray<T>& a)
{
  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::MultElt(const TArray<T>&)")) ;

  T * pe;
  const T * pea;
  uint_8 j,k,ka;
  if ((AvgStep() > 0) && (a.AvgStep() > 0) )   {  // regularly spaced elements
    uint_8 step = AvgStep();
    uint_8 stepa = a.AvgStep();
    uint_8 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 ...
    uint_4 ax = MaxSizeKA();
    uint_8 step = Step(ax);
    uint_8 stepa = a.Step(ax);
    uint_8 gpas = Size(ax)*step;
    for(j=0; j<totsize_; j += Size(ax))  {
      pe = mNDBlock.Begin()+Offset(j);
      pea = a.DataBlock().Begin()+a.Offset(j);
      for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] *= pea[ka]; 
    }
  }
  return(*this);
}

template <class T>
TArray<T>& TArray<T>::DivElt(const TArray<T>& a)
{
  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::DivElt(const TArray<T>&)")) ;

  T * pe;
  const T * pea;
  uint_8 j,k,ka;
  if ((AvgStep() > 0) && (a.AvgStep() > 0) )   {  // regularly spaced elements
    uint_8 step = AvgStep();
    uint_8 stepa = a.AvgStep();
    uint_8 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 ...
    uint_4 ax = MaxSizeKA();
    uint_8 step = Step(ax);
    uint_8 stepa = a.Step(ax);
    uint_8 gpas = Size(ax)*step;
    for(j=0; j<totsize_; j += Size(ax))  {
      pe = mNDBlock.Begin()+Offset(j);
      pea = a.DataBlock().Begin()+a.Offset(j);
      for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] /= pea[ka]; 
    }
  }
  return(*this);
}


// ----------------------------------------------------
//          Application d'une fonction
// ----------------------------------------------------
template <class T>
TArray<T>& TArray<T>::ApplyFunction(Arr_DoubleFunctionOfX f)
{
  T * pe;
  uint_8 j,k;
  if (AvgStep() > 0)   {  // regularly spaced elements
    uint_8 step = AvgStep();
    uint_8 maxx = totsize_*step;
    pe = Data();
    for(k=0; k<maxx; k+=step )  pe[k] = (T)(f((double)pe[k]));
  }
  else {    // Non regular data spacing ...
    uint_4 ka = MaxSizeKA();
    uint_8 step = Step(ka);
    uint_8 gpas = Size(ka)*step;
    for(j=0; j<totsize_; j += Size(ka))  {
      pe = mNDBlock.Begin()+Offset(j);
      for(k=0; k<gpas; k+=step)  pe[k] = (T)(f((double)pe[k]));
    }
  }
  return(*this);
}

template <class T>
TArray<T>& TArray<T>::ApplyFunction(Arr_FloatFunctionOfX f)
{
  T * pe;
  uint_8 j,k;
  if (AvgStep() > 0)   {  // regularly spaced elements
    uint_8 step = AvgStep();
    uint_8 maxx = totsize_*step;
    pe = Data();
    for(k=0; k<maxx; k+=step )  pe[k] = (T)(f((float)pe[k]));
  }
  else {    // Non regular data spacing ...
    uint_4 ka = MaxSizeKA();
    uint_8 step = Step(ka);
    uint_8 gpas = Size(ka)*step;
    for(j=0; j<totsize_; j += Size(ka))  {
      pe = mNDBlock.Begin()+Offset(j);
      for(k=0; k<gpas; k+=step)  pe[k] = (T)(f((float)pe[k]));
    }
  }
  return(*this);
}

TArray< complex<r_4> >& TArray< complex<r_4> >::ApplyFunction(Arr_DoubleFunctionOfX f)
{
  throw( NotAvailableOperation("TArray< complex<r_4> >::ApplyFunction() - Unsupported operation ") );
}

TArray< complex<r_4> >& TArray< complex<r_4> >::ApplyFunction(Arr_FloatFunctionOfX f)
{
  throw(NotAvailableOperation( "TArray< complex<r_4> >::ApplyFunction() - Unsupported operation ") );
}

TArray< complex<r_8> >& TArray< complex<r_8> >::ApplyFunction(Arr_DoubleFunctionOfX f)
{
  throw(NotAvailableOperation("TArray< complex<r_8> >::ApplyFunction() - Unsupported operation ") );
}

TArray< complex<r_8> >& TArray< complex<r_8> >::ApplyFunction(Arr_FloatFunctionOfX f)
{
  throw(NotAvailableOperation("TArray< complex<r_8> >::ApplyFunction() - Unsupported operation ") );
}

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

template <class T>
void TArray<T>::Show(ostream& os, bool si) const
{
  os << " TArray< " << typeid(T).name() << " NDim= " << ndim_
     << "(" << typeid(*this).name() << " )" << endl;
  os << "TotSize=" << totsize_ << " Size(X*Y*...)= " ;
  for(int k=0; k<ndim_; k++) { 
    os << size_[k];
    if (k<ndim_-1)  os << " x ";
  }
  os << endl;
  os <<  " Step(X Y ...)= " ;
  for(int k=0; k<ndim_; k++)     os << step_[k] << "  " ;
  os << " Offset= " << offset_  << "\n " << endl;
  if (si && (mInfo != NULL)) os << (*mInfo) << endl;
 
}

template <class T>
void TArray<T>::Print(ostream& os, int_4 maxprt, bool si) const
{
  if (maxprt < 0)  maxprt = max_nprt_;
  uint_4 npr = 0;
  Show(os, si);
  int 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 >= maxprt) {
              if (npr < totsize_)  os << "\n     .... " << endl; return;
            }
          }
          os << endl;
        }
      }
    }
  }
  os << "\n" << endl;
}

template <class T>
DVList& TArray<T>::Info()
{
if (mInfo == NULL)  mInfo = new DVList;
return(*mInfo);
}


template <class T>
bool TArray<T>::UpdateSizes(uint_4 ndim, const uint_4 * siz, uint_4 step, uint_8 offset, string & exmsg)
{
  if (ndim >= TARRAY_MAXNDIMS) {
    exmsg += " NDim Error";  return false;
  }
  if (step < 1) {
    exmsg += " Step(=0) Error";  return false;
  }

  minstep_ = moystep_ = step;

  // Flagging bad updates ...
  ndim_ = 0;

  totsize_ = 1;
  int k;
  for(k=0; k<TARRAY_MAXNDIMS; k++) {
    size_[k] = 1;
    step_[k] = 0;
  }
  for(k=0; k<ndim; k++) {
    size_[k] = siz[k] ;
    step_[k] = totsize_*step;
    totsize_ *= size_[k];
  }
  if (totsize_ < 1) {
    exmsg += " Size Error";  return false;
  }
  offset_ = offset;
  // Default for matrices : Lines = Y , Columns = X
  macoli_ = 0;
  marowi_ = 1;
  // Update OK 
  ndim_ = ndim;
  return true;
}

template <class T>
bool TArray<T>::UpdateSizes(uint_4 ndim, const uint_4 * siz, const uint_4 * step, uint_8 offset, string & exmsg)
{
  if (ndim >= TARRAY_MAXNDIMS) {
    exmsg += " NDim Error";  return false;
  }

  // Flagging bad updates ...
  ndim_ = 0;

  totsize_ = 1;
  int k;
  for(k=0; k<TARRAY_MAXNDIMS; k++) {
    size_[k] = 1;
    step_[k] = 0;
  }
  uint_4 minstep = step[0];
  for(k=0; k<ndim; k++) {
    size_[k] = siz[k] ;
    step_[k] = step[k];
    totsize_ *= size_[k];
    if (step_[k] < minstep)  minstep = step_[k];
  }
  if (minstep < 1) {
    exmsg += " Step(=0) Error";  return false;
  }
  if (totsize_ < 1) {
    exmsg += " Size Error";  return false;
  }
  uint_8 plast = 0;
  for(k=0; k<ndim; k++)   plast += (siz[k]-1)*step[k];
  if (plast == minstep*totsize_ )  moystep_ = minstep;
  else moystep_ = 0;
  minstep_ = minstep;
  offset_ = offset;
  // Default for matrices : Lines = Y , Columns = X
  macoli_ = 0;
  marowi_ = 1;
  // Update OK 
  ndim_ = ndim;
  return true;
}

template <class T>
bool TArray<T>::UpdateSizes(const TArray<T>& a, string & exmsg)
{
  if (a.ndim_ >= TARRAY_MAXNDIMS) {
    exmsg += " NDim Error";  return false;
  }

  // Flagging bad updates ...
  ndim_ = 0;

  totsize_ = 1;
  int k;
  for(k=0; k<TARRAY_MAXNDIMS; k++) {
    size_[k] = 1;
    step_[k] = 0;
  }
  uint_4 minstep = a.step_[0];
  for(k=0; k<a.ndim_; k++) {
    size_[k] = a.size_[k] ;
    step_[k] = a.step_[k];
    totsize_ *= size_[k];
    if (step_[k] < minstep)  minstep = step_[k];
  }
  if (minstep < 1) {
    exmsg += " Step(=0) Error";  return false;
  }
  if (totsize_ < 1) {
    exmsg += " Size Error";  return false;
  }

  minstep_ = a.minstep_;
  moystep_ = a.moystep_;
  offset_ = a.offset_;
  macoli_ = a.macoli_;
  marowi_ = a.marowi_;
  // Update OK 
  ndim_ = a.ndim_;
  return true;
}

template <class T>
void TArray<T>::CloneOrShare(const TArray<T>& a)
{
  string exmsg = "TArray<T>::CloneOrShare()";
  if (!UpdateSizes(a.ndim_, a.size_, a.step_, a.offset_, exmsg))  throw( ParmError(exmsg) );
  mNDBlock.CloneOrShare(a.mNDBlock);
}

template <class T>
void TArray<T>::Share(const TArray<T>& a)
{
  string exmsg = "TArray<T>::Share()";
  if (!UpdateSizes(a.ndim_, a.size_, a.step_, a.offset_, exmsg))  throw( ParmError(exmsg) );
  mNDBlock.Share(a.mNDBlock);
}


template <class T>
void TArray<T>::SubArray(TArray<T> & ra, uint_4 ndim, uint_4 * siz, uint_4 * pos, uint_4 * step)
{
  if ( (ndim > ndim_) || (ndim < 1) ) throw(SzMismatchError("TArray<T>::SubArray( ... ) NDim Error") );
  int k;
  for(k=0; k<ndim; k++) 
    if ( (siz[k]*step[k]+pos[k]) > size_[k] )  
      throw(SzMismatchError("TArray<T>::SubArray( ... ) Size/Pos Error") );
  (ra.mNDBlock).Share(mNDBlock);
  uint_8 offset = offset_;
  for(k=0; k<ndim_; k++) { 
    offset += pos[k]*step_[k]; 
    step[k] *= step_[k];
  }
  string exm = "TArray<T>::SubArray() ";
  if (!ra.UpdateSizes(ndim, siz, step, offset,  exm))
     throw( ParmError(exm) );
  (ra.mNDBlock).Share(mNDBlock);
  return;
}

///////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////
#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template TArray<uint_1>
#pragma define_template TArray<uint_2>
#pragma define_template TArray<int_2>
#pragma define_template TArray<int_4>
#pragma define_template TArray<int_8>
#pragma define_template TArray<uint_4>
#pragma define_template TArray<uint_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<uint_2>;
template class TArray<int_2>;
template class TArray<int_4>;
template class TArray<int_8>;
template class TArray<uint_4>;
template class TArray<uint_8>;
template class TArray<r_4>;
template class TArray<r_8>;
template class TArray< complex<r_4> >;
template class TArray< complex<r_8> >;
#endif