Index: trunk/SophyaLib/TArray/tarray.cc
===================================================================
--- trunk/SophyaLib/TArray/tarray.cc	(revision 787)
+++ trunk/SophyaLib/TArray/tarray.cc	(revision 804)
@@ -18,5 +18,5 @@
 template <class T>
 TArray<T>::TArray()
-  : mNDBlock() , BaseArray()
+  : BaseArray() , mNDBlock() 
 // Default constructor
 {
@@ -24,6 +24,6 @@
 
 template <class T>
-TArray<T>::TArray(uint_4 ndim, uint_4 * siz, uint_4 step)
-  : mNDBlock(ComputeTotalSize(ndim, siz, step, 1)) , BaseArray()
+TArray<T>::TArray(uint_4 ndim, const uint_4 * siz, uint_4 step)
+  : BaseArray() , mNDBlock(ComputeTotalSize(ndim, siz, step, 1))
 {
   string exmsg = "TArray<T>::TArray(uint_4, uint_4 *, uint_4)";
@@ -32,11 +32,14 @@
 
 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)) , BaseArray()
+TArray<T>::TArray(uint_4 nx, uint_4 ny, uint_4 nz, uint_4 nt, uint_4 nu)
+  : BaseArray() , mNDBlock(nx*((ny>0)?ny:1)*((nz>0)?nz:1)*((nt>0)?nt:1)*((nu>0)?nu:1)) 
 {
   uint_4 size[BASEARRAY_MAXNDIMS];
   size[0] = nx;   size[1] = ny;   size[2] = nz;
+  size[3] = nt;   size[4] = nu;
   int ndim = 1;
-  if ((size[1] > 0) && (size[2] > 0)) ndim = 3;
+  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;
@@ -46,6 +49,6 @@
 
 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) , BaseArray()
+TArray<T>::TArray(uint_4 ndim, const uint_4 * siz, NDataBlock<T> & db, bool share, uint_4 step, uint_8 offset)
+  : BaseArray() , mNDBlock(db, share) 
 {
   string exmsg = "TArray<T>::TArray(uint_4, uint_4 *,  NDataBlock<T> & ... )";
@@ -55,6 +58,6 @@
 
 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) , BaseArray()
+TArray<T>::TArray(uint_4 ndim, const uint_4 * siz, T* values, uint_4 step, uint_8 offset, Bridge* br)
+  : BaseArray() , mNDBlock(ComputeTotalSize(ndim, siz, step, 1), values, br) 
 {
   string exmsg = "TArray<T>::TArray(uint_4, uint_4 *, T* ... )";
@@ -64,5 +67,5 @@
 template <class T>
 TArray<T>::TArray(const TArray<T>& a)
-  : mNDBlock(a.mNDBlock) , BaseArray()
+  : BaseArray() , mNDBlock(a.mNDBlock) 
 {
   string exmsg = "TArray<T>::TArray(const TArray<T>&)";
@@ -73,5 +76,5 @@
 template <class T>
 TArray<T>::TArray(const TArray<T>& a, bool share)
-  : mNDBlock(a.mNDBlock, share) , BaseArray()
+  : BaseArray() , mNDBlock(a.mNDBlock, share) 
 {
   string exmsg = "TArray<T>::TArray(const TArray<T>&, bool)";
@@ -87,5 +90,5 @@
 
 template <class T>
-TArray<T>& TArray<T>::operator = (const TArray<T>& a)
+TArray<T>& TArray<T>::Set(const TArray<T>& a)
 {
   if (this != &a) { 
@@ -171,8 +174,29 @@
 
 
+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(*this, true);
+    ra.SetTemp(true);
+    return(ra);
+  }
+  else {
+    TArray<T> ra(ndim_, size_, 1);
+    ra.CopyElt(*this);
+    ra.SetTemp(true);
+    return(ra);
+  }
+}
+
 // SubArrays 
-template <class T>
-TArray<T> TArray<T>::operator () (Range rx, Range ry, Range rz, Range rt, Range ru)
-{
+// $CHECK$ Reza 03/2000  Doit-on declarer cette methode const ?
+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 ! ");
   uint_4 ndim = 0;
   uint_4 size[BASEARRAY_MAXNDIMS];
@@ -199,5 +223,5 @@
   ndim = ndim_;
   TArray<T> ra;
-  SubArray(ra, ndim, size, pos, step); 
+  UpdateSubArraySizes(ra, ndim, size, pos, step); 
   ra.DataBlock().Share(this->DataBlock());
   ra.SetTemp(true);
@@ -210,6 +234,8 @@
 //  Possibilite de // , vectorisation 
 template <class T>
-TArray<T>& TArray<T>::operator = (Sequence seq)
-{
+TArray<T>& TArray<T>::Set(Sequence seq)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::Set(Sequence ) - Not Allocated Array ! ");
   T * pe;
   uint_8 j,k;
@@ -234,6 +260,8 @@
 
 template <class T>
-TArray<T>& TArray<T>::operator = (T x)
-{
+TArray<T>& TArray<T>::Set(T x)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::Set(T )  - Not Allocated Array ! ");
   T * pe;
   uint_8 j,k;
@@ -257,6 +285,8 @@
 
 template <class T>
-TArray<T>& TArray<T>::operator += (T x)
-{
+TArray<T>& TArray<T>::Add(T x)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::Add(T )  - Not Allocated Array ! ");
   T * pe;
   uint_8 j,k;
@@ -280,6 +310,8 @@
 
 template <class T>
-TArray<T>& TArray<T>::operator -= (T x)
-{
+TArray<T>& TArray<T>::Sub(T x)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::Sub(T )  - Not Allocated Array ! ");
   T * pe;
   uint_8 j,k;
@@ -303,6 +335,8 @@
 
 template <class T>
-TArray<T>& TArray<T>::operator *= (T x)
-{
+TArray<T>& TArray<T>::Mul(T x)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::Mul(T )  - Not Allocated Array ! ");
   T * pe;
   uint_8 j,k;
@@ -326,6 +360,8 @@
 
 template <class T>
-TArray<T>& TArray<T>::operator /= (T x)
-{
+TArray<T>& TArray<T>::Div(T x)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::Div(T )  - Not Allocated Array ! ");
   T * pe;
   uint_8 j,k;
@@ -349,9 +385,62 @@
 
 
+template <class T>
+TArray<T>& TArray<T>::SubInv(T x)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::SubInv(T )  - Not Allocated Array ! ");
+  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-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] = x-pe[k];
+    }
+  }
+  return(*this);
+}
+
+template <class T>
+TArray<T>& TArray<T>::DivInv(T x)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::DivInv(T )  - Not Allocated Array ! ");
+  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/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] = x/pe[k];
+    }
+  }
+  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>&)")) ;
+TArray<T>& TArray<T>::AddElt(const TArray<T>& a)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::AddElt(const TArray<T>& )  - Not Allocated Array ! ");
+  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::AddElt(const TArray<T>&)")) ;
 
   T * pe;
@@ -381,7 +470,9 @@
 
 template <class T>
-TArray<T>& TArray<T>::operator -= (const TArray<T>& a)
-{
-  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::operator -= (const TArray<T>&)")) ;
+TArray<T>& TArray<T>::SubElt(const TArray<T>& a)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::SubElt(const TArray<T>& )  - Not Allocated Array ! ");
+  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::SubElt(const TArray<T>&)")) ;
 
   T * pe;
@@ -412,7 +503,9 @@
 
 template <class T>
-TArray<T>& TArray<T>::MultElt(const TArray<T>& a)
-{
-  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::MultElt(const TArray<T>&)")) ;
+TArray<T>& TArray<T>::MulElt(const TArray<T>& a)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::MulElt(const TArray<T>& )  - Not Allocated Array ! ");
+  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::MulElt(const TArray<T>&)")) ;
 
   T * pe;
@@ -441,7 +534,10 @@
 }
 
+
 template <class T>
 TArray<T>& TArray<T>::DivElt(const TArray<T>& a)
 {
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::DivElt(const TArray<T>& )  - Not Allocated Array ! ");
   if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::DivElt(const TArray<T>&)")) ;
 
@@ -471,4 +567,91 @@
 }
 
+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 ! ");
+  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);
+}
+
+
+// Somme et produit des 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;
+  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 )  ret += 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)  ret += pe[k] ;
+    }
+  }
+  return ret;
+}
+
+template <class T>
+T TArray<T>::Product() const
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::Product()  - Not Allocated Array ! ");
+  T ret=0;
+  const 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 )  ret *= 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)  ret *= pe[k] ;
+    }
+  }
+  return ret;
+}
+
+
 
 // ----------------------------------------------------
@@ -487,7 +670,7 @@
 {
   if (maxprt < 0)  maxprt = max_nprt_;
-  uint_4 npr = 0;
+  int_4 npr = 0;
   Show(os, si);
-  int k0,k1,k2,k3,k4;
+  uint_4 k0,k1,k2,k3,k4;
   for(k4=0; k4<size_[4]; k4++) {
     if (size_[4] > 1) cout << "\n ----- Dimension 5 (U) K4= " << k4 << endl;
@@ -535,11 +718,13 @@
 ///////////////////////////////////////////////////////////////
 #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_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>
@@ -549,11 +734,13 @@
 
 #if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
-template class TArray<uint_1>;
+/*
+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_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>;
