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

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


// Classe utilitaires 
Sequence::Sequence(double start, double step)
{
  start_ = start;
  step_ = step;
}

// 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;
  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>::SubArray(uint_4 ndim, uint_4 * siz, uint_4 * pos)
{
  if ( (ndim > ndim_) || (ndim < 1) ) throw(SzMismatchError("TArray<T>::SubArray( ... ) NDim Error") );
  int k;
  for(k=0; k<ndim; k++) 
    if ( (siz[k]+pos[k]) > size_[k] )  
      throw(SzMismatchError("TArray<T>::SubArray( ... ) Size/Pos Error (1)") );
  for(k=ndim; k<ndim_; k++)   
    if ( (pos[k]) >= size_[k] )  
      throw(SzMismatchError("TArray<T>::SubArray( ... ) Size/Pos Error (2)") );
  TArray<T> ra(*this);
  ra.ndim_ = ndim;
  for(k=ndim; k<TARRAY_MAXNDIMS; k++) {
    ra.size_[k] = 1;
    ra.step_[k] = 0;
  }
  ra.offset_ = offset_;
  for(k=0; k<ndim_; k++) ra.offset_ += pos[k]*step_[k]; 
  ra.SetTemp(true);
  return(ra);   
}

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

//   ...... 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 = Data();
  uint_8 step = Step();
  for(uint_8 k=0; k<totsize_; k++ )  pe[k*step] = seq(k);
  return(*this);
}

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

template <class T>
TArray<T>& TArray<T>::operator = (T x)
{
  T * pe = Data();
  uint_8 step = Step();
  for(uint_8 k=0; k<totsize_*step; k += step)  pe[k] = x;
  return(*this);
}

template <class T>
TArray<T>& TArray<T>::operator += (T x)
{
  T * pe = Data();
  uint_8 step = Step();
  for(uint_8 k=0; k<totsize_*step; k += step)  pe[k] += x;
  return(*this);
}

template <class T>
TArray<T>& TArray<T>::operator -= (T x)
{
  T * pe = Data();
  uint_8 step = Step();
  for(uint_8 k=0; k<totsize_*step; k += step)  pe[k] -= x;
  return(*this);
}

template <class T>
TArray<T>& TArray<T>::operator *= (T x)
{
  T * pe = Data();
  uint_8 step = Step();
  for(uint_8 k=0; k<totsize_; k += step)  pe[k] *= x;
  return(*this);
}

template <class T>
TArray<T>& TArray<T>::operator /= (T x)
{
  T * pe = Data();
  uint_8 step = Step();
  for(uint_8 k=0; k<totsize_*step; 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 = Data();
  const T * pea = a.Data();
  uint_8 step = Step();
  uint_8 stepa = a.Step();
  uint_8 k, ka;
  k = ka = 0;
  for(k=0; k<totsize_; k += step)  {
    pe[k] += pea[ka];
    ka += stepa;
  }
  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 = Data();
  const T * pea = a.Data();
  uint_8 step = Step();
  uint_8 stepa = a.Step();
  uint_8 k, ka;
  k = ka = 0;
  for(k=0; k<totsize_; k += step)  {
    pe[k] -= pea[ka];
    ka += stepa;
  }
  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 = Data();
  const T * pea = a.Data();
  uint_8 step = Step();
  uint_8 stepa = a.Step();
  uint_8 k, ka;
  k = ka = 0;
  for(k=0; k<totsize_; k += step)  {
    pe[k] *= pea[ka];
    ka += stepa;
  }
  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 = Data();
  const T * pea = a.Data();
  uint_8 step = Step();
  uint_8 stepa = a.Step();
  uint_8 k, ka;
  k = ka = 0;
  for(k=0; k<totsize_; k += step)  {
    pe[k] /= pea[ka];
    ka += stepa;
  }
  return(*this);
}


// ----------------------------------------------------
//          Application d'une fonction
// ----------------------------------------------------
template <class T>
TArray<T>& TArray<T>::ApplyFunction(Arr_DoubleFunctionOfX f)
{
  T * pe = Data();
  uint_8 step = Step();
  for(uint_8 k=0; k<totsize_*step; k += step)  pe[k] = (T)(f((double)pe[k]));
  return(*this);
}

template <class T>
TArray<T>& TArray<T>::ApplyFunction(Arr_FloatFunctionOfX f)
{
  T * pe = Data();
  uint_8 step = Step();
  for(uint_8 k=0; k<totsize_*step; 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 << "\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;
    for(k3=0; k3<size_[3]; k3++) {
      if (size_[3] > 1) cout << "\n ----- Dimension 4 (T) K3= " << k3;
      for(k2=0; k2<size_[2]; k2++) {
        if (size_[2] > 1) cout << "\n ----- Dimension 3 (Z) K2= " << k2;
        for(k1=0; k1<size_[1]; k1++) {
          if ( (size_[1] > 1) && (size_[0] > 10) ) cout << "----- Dimension 2 (Y) K1= " << k1;
          for(k0=0; k0<size_[0]; k0++) {
            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_ = 1;

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

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


///////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////
#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