Index: trunk/SophyaLib/TArray/tarray.cc
===================================================================
--- trunk/SophyaLib/TArray/tarray.cc	(revision 2569)
+++ trunk/SophyaLib/TArray/tarray.cc	(revision 2575)
@@ -482,8 +482,8 @@
 //! 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 ! ");
+TArray<T>& TArray<T>::SetCst(T x)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::SetCst(T )  - Not Allocated Array ! ");
   T * pe;
   sa_size_t j,k;
@@ -528,5 +528,5 @@
   const T * pe;
   T * per;
-  sa_size_t j,k,kr;
+  sa_size_t j,k;
   if (smo && (IsPacked() > 0) && (res.IsPacked() > 0))   {  // regularly spaced elements
     sa_size_t maxx = totsize_;
@@ -538,10 +538,10 @@
     int_4 ax,axr;
     sa_size_t step, stepr;
-    sa_size_t gpas, naxr;
-    GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxr);
-    for(j=0; j<naxr; j++)  {
+    sa_size_t gpas, naxa;
+    GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa);
+    for(j=0; j<naxa; j++)  {
       pe = mNDBlock.Begin()+Offset(ax,j);
       per = res.DataBlock().Begin()+res.Offset(axr,j);
-      for(k=0, kr=0;  k<gpas;  k+=step, kr+=stepr)  per[kr] = pe[k]+x; 
+      for(k=0;  k<gpas;  k+=step, pe+=step, per+=stepr)  *per = *pe+x; 
     }
   }
@@ -557,5 +557,5 @@
 \param x : constant to subtract from the array elements
 \param res : Output array containing the result (res=this+x or res=x-this). 
-\param fginv == true : Invert subtraction argument order (*this = x-(*this)) 
+\param fginv == true : Invert subtraction argument order (res = x-(*this)) 
 */
 template <class T>
@@ -571,5 +571,5 @@
   const T * pe;
   T * per;
-  sa_size_t j,k,kr;
+  sa_size_t j,k;
   if (smo && (IsPacked() > 0) && (res.IsPacked() > 0))   {  // regularly spaced elements
     sa_size_t maxx = totsize_;
@@ -584,13 +584,13 @@
     int_4 ax,axr;
     sa_size_t step, stepr;
-    sa_size_t gpas, naxr;
-    GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxr);
-    for(j=0; j<naxr; j++)  {
+    sa_size_t gpas, naxa;
+    GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa);
+    for(j=0; j<naxa; j++)  {
       pe = mNDBlock.Begin()+Offset(ax,j);
       per = res.DataBlock().Begin()+res.Offset(axr,j);
       if (!fginv) 
-        for(k=0, kr=0;  k<gpas;  k+=step, kr+=stepr)  per[kr] = pe[k]-x; 
-      else
-        for(k=0, kr=0;  k<gpas;  k+=step, kr+=stepr)  per[kr] = x-pe[k]; 
+        for(k=0;  k<gpas;  k+=step, pe+=step, per+=stepr)  *per = *pe-x; 
+       else
+        for(k=0;  k<gpas;  k+=step, pe+=step, per+=stepr)  *per = x-*pe; 
     }
   }
@@ -619,5 +619,5 @@
   const T * pe;
   T * per;
-  sa_size_t j,k,kr;
+  sa_size_t j,k;
   if (smo && (IsPacked() > 0) && (res.IsPacked() > 0))   {  // regularly spaced elements
     sa_size_t maxx = totsize_;
@@ -629,10 +629,10 @@
     int_4 ax,axr;
     sa_size_t step, stepr;
-    sa_size_t gpas, naxr;
-    GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxr);
-    for(j=0; j<naxr; j++)  {
+    sa_size_t gpas, naxa;
+    GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa);
+    for(j=0; j<naxa; j++)  {
       pe = mNDBlock.Begin()+Offset(ax,j);
       per = res.DataBlock().Begin()+res.Offset(axr,j);
-      for(k=0, kr=0;  k<gpas;  k+=step, kr+=stepr)  per[kr] = pe[k]*x; 
+      for(k=0;  k<gpas;  k+=step, pe+=step, per+=stepr)  *per = (*pe)*x; 
     }
   }
@@ -648,5 +648,5 @@
 \param x : Array elements are divied by x
 \param res : Output array containing the result (res=(*this)/x or res=x/(*this)). 
-\param fginv == true : Invert the operation order (res = x/(*this)) 
+\param fginv == true : Invert the division argument order (res = x/(*this)) 
 */
 template <class T>
@@ -677,13 +677,13 @@
     int_4 ax,axr;
     sa_size_t step, stepr;
-    sa_size_t gpas, naxr;
-    GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxr);
-    for(j=0; j<naxr; j++)  {
+    sa_size_t gpas, naxa;
+    GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa);
+    for(j=0; j<naxa; j++)  {
       pe = mNDBlock.Begin()+Offset(ax,j);
       per = res.DataBlock().Begin()+res.Offset(axr,j);
       if (!fginv) 
-        for(k=0, kr=0;  k<gpas;  k+=step, kr+=stepr)  per[kr] = pe[k]/x; 
-      else
-        for(k=0, kr=0;  k<gpas;  k+=step, kr+=stepr)  per[kr] = x/pe[k]; 
+        for(k=0;  k<gpas;  k+=step, pe+=step, per+=stepr)  *per = (*pe)/x; 
+       else
+        for(k=0;  k<gpas;  k+=step, pe+=step, per+=stepr)  *per = x/(*pe); 
     }
   }
@@ -710,5 +710,5 @@
   const T * pe;
   T * per;
-  sa_size_t j,k,kr;
+  sa_size_t j,k;
   if (smo && (IsPacked() > 0) && (res.IsPacked() > 0))   {  // regularly spaced elements
     sa_size_t maxx = totsize_;
@@ -720,10 +720,10 @@
     int_4 ax,axr;
     sa_size_t step, stepr;
-    sa_size_t gpas, naxr;
-    GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxr);
-    for(j=0; j<naxr; j++)  {
+    sa_size_t gpas, naxa;
+    GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa);
+    for(j=0; j<naxa; j++)  {
       pe = mNDBlock.Begin()+Offset(ax,j);
       per = res.DataBlock().Begin()+res.Offset(axr,j);
-      for(k=0, kr=0;  k<gpas;  k+=step, kr+=stepr)  per[kr] = -pe[k]; 
+      for(k=0;  k<gpas;  k+=step, pe+=step, per+=stepr)  *per = -(*pe); 
     }
   }
@@ -732,18 +732,324 @@
 
 //  >>>> 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 ! ");
+
+//! Two TArrays element by element addition
+/*! 
+Perform element by element addition of the source array (this) and the \b a array 
+and store the result in \b res (res = *this+a). The source and argument arrays (this, a)
+should have the same sizes.
+If not initially allocated, the output array \b res is automatically 
+resized as a packed array with the same sizes as the source (this) array.
+Returns a reference to the output array \b res. 
+\param a : Array to be added to the source array.
+\param res : Output array containing the result (res=this+a). 
+*/
+template <class T>
+TArray<T>& TArray<T>::AddElt(const TArray<T>& a, TArray<T>& res) const 
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::AddElt(...) - Not allocated source array ! ");
+  bool smoa;
+  if (!CompareSizes(a, smoa)) 
+    throw(SzMismatchError("TArray<T>::AddElt(...) SizeMismatch(this,a)")) ;
+  if (res.NbDimensions() < 1)  res.SetSize(*this, true, false);
+  bool smor;
+  if (!CompareSizes(res, smor)) 
+    throw(SzMismatchError("TArray<T>::AddElt(...) SizeMismatch(this,res) ")) ;
+
+  bool smora;
+  a.CompareSizes(res, smora);
+
+  bool smo = smoa && smor;  // The three arrays have same memory organisation
+
+  const T * pe;
+  const T * pea;
+  T * per;
+  sa_size_t j,k;
+  if (smo && IsPacked() && a.IsPacked() && res.IsPacked() )   {  // All packed arrays
+    sa_size_t maxx = totsize_;
+    pe = Data();
+    pea = a.Data();
+    per = res.Data();
+    for(k=0; k<maxx;  k++, pe++, pea++, per++)  *per = *pe + *pea ;
+  }
+  else {    // Non regular data spacing ...
+    int_4 ax,axa,axr;
+    sa_size_t step, stepa;
+    sa_size_t gpas, naxa;
+    sa_size_t stepr, stgpas;
+    if ( !smo && smora ) {   // same mem-org for a,res , different from this
+      a.GetOpeParams(*this, smo, axa, ax, stepa, step, gpas, naxa);
+      a.GetOpeParams(res, smo, axa, axr, stepa, stepr, gpas, naxa);
+      stgpas = stepa;
+    }
+    else { // same mem-org for all, or same (this,a) OR same(this,res)
+      GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa);
+      GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa);
+      stgpas = step;
+    }
+    for(j=0; j<naxa; j++)  {
+      pe = mNDBlock.Begin()+Offset(ax,j);
+      pea = a.DataBlock().Begin()+a.Offset(axa,j);
+      per = res.DataBlock().Begin()+res.Offset(axr,j);
+      for(k=0; k<gpas; k+=stgpas, pe+=step, pea+=stepa, per+=stepr) *per = *pe + *pea ; 
+    }
+  }
+  return(res);
+}
+
+//! Two TArrays element by element subtraction
+/*! 
+Perform element by element subtraction of the source array (this) and the \b a array 
+and the store result in \b res (res = *this-a or res=a-(*this)). 
+The source and argument arrays (this, a) should have the same sizes.
+If not initially allocated, the output array \b res is automatically 
+resized as a packed array with the same sizes as the source (this) array.
+Returns a reference to the output array \b res. 
+\param a : Array to be added to the source array.
+\param res : Output array containing the result (res=*this+x). 
+\param fginv == true : Invert subtraction argument order (res = a-(*this))  
+*/
+
+template <class T>
+TArray<T>& TArray<T>::SubElt(const TArray<T>& a, TArray<T>& res, bool fginv) const
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::SubElt(...) - Not allocated source array ! ");
+  bool smoa;
+  if (!CompareSizes(a, smoa)) 
+    throw(SzMismatchError("TArray<T>::SubElt(...) SizeMismatch(this,a)")) ;
+  if (res.NbDimensions() < 1)  res.SetSize(*this, true, false);
+  bool smor;
+  if (!CompareSizes(res, smor)) 
+    throw(SzMismatchError("TArray<T>::SubElt(...) SizeMismatch(this,res) ")) ;
+
+  bool smora;
+  a.CompareSizes(res, smora);
+
+  bool smo = smoa && smor;  // The three arrays have same memory organisation
+
+  const T * pe;
+  const T * pea;
+  T * per;
+  sa_size_t j,k;
+  if (smo && IsPacked() && a.IsPacked() && res.IsPacked() )   {  // All packed arrays
+    sa_size_t maxx = totsize_;
+    pe = Data();
+    pea = a.Data();
+    per = res.Data();
+    if (!fginv) 
+      for(k=0; k<maxx;  k++, pe++, pea++, per++)  *per = *pe - *pea ;
+    else 
+      for(k=0; k<maxx;  k++, pe++, pea++, per++)  *per = *pea - *pe ;
+  }
+  else {    // Non regular data spacing ...
+    int_4 ax,axa,axr;
+    sa_size_t step, stepa;
+    sa_size_t gpas, naxa;
+    sa_size_t stepr, stgpas;
+    if ( !smo && smora ) {   // same mem-org for a,res , different from this
+      a.GetOpeParams(*this, smo, axa, ax, stepa, step, gpas, naxa);
+      a.GetOpeParams(res, smo, axa, axr, stepa, stepr, gpas, naxa);
+      stgpas = stepa;
+    }
+    else { // same mem-org for all, or same (this,a) OR same(this,res)
+      GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa);
+      GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa);
+      stgpas = step;
+    }
+    for(j=0; j<naxa; j++)  {
+      pe = mNDBlock.Begin()+Offset(ax,j);
+      pea = a.DataBlock().Begin()+a.Offset(axa,j);
+      per = res.DataBlock().Begin()+res.Offset(axr,j);
+      if (!fginv) 
+        for(k=0; k<gpas; k+=stgpas, pe+=step, pea+=stepa, per+=stepr) *per = *pe - *pea ; 
+      else 
+        for(k=0; k<gpas; k+=stgpas, pe+=step, pea+=stepa, per+=stepr) *per = *pea - *pea ; 
+    }
+  }
+  return(res);
+}
+
+
+//! Two TArrays element by element multiplication
+/*! 
+Perform element by element multiplication of the source array (this) and the \b a array 
+and store the result in \b res (res = *this*a). The source and argument arrays (this, a)
+should have the same sizes.
+If not initially allocated, the output array \b res is automatically 
+resized as a packed array with the same sizes as the source (this) array.
+Returns a reference to the output array \b res. 
+\param a : Array to be added to the source array.
+\param res : Output array containing the result (res=(*this)*a). 
+*/
+template <class T>
+TArray<T>& TArray<T>::MulElt(const TArray<T>& a, TArray<T>& res) const 
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::MulElt(...) - Not allocated source array ! ");
+  bool smoa;
+  if (!CompareSizes(a, smoa)) 
+    throw(SzMismatchError("TArray<T>::MulElt(...) SizeMismatch(this,a)")) ;
+  if (res.NbDimensions() < 1)  res.SetSize(*this, true, false);
+  bool smor;
+  if (!CompareSizes(res, smor)) 
+    throw(SzMismatchError("TArray<T>::MulElt(...) SizeMismatch(this,res) ")) ;
+
+  bool smora;
+  a.CompareSizes(res, smora);
+
+  bool smo = smoa && smor;  // The three arrays have same memory organisation
+
+  const T * pe;
+  const T * pea;
+  T * per;
+  sa_size_t j,k;
+  if (smo && IsPacked() && a.IsPacked() && res.IsPacked() )   {  // All packed arrays
+    sa_size_t maxx = totsize_;
+    pe = Data();
+    pea = a.Data();
+    per = res.Data();
+    for(k=0; k<maxx;  k++, pe++, pea++, per++ )  *per = (*pe) * (*pea) ;
+  }
+  else {    // Non regular data spacing ...
+    int_4 ax,axa,axr;
+    sa_size_t step, stepa;
+    sa_size_t gpas, naxa;
+    sa_size_t stepr, stgpas;
+    if ( !smo && smora ) {   // same mem-org for a,res , different from this
+      a.GetOpeParams(*this, smo, axa, ax, stepa, step, gpas, naxa);
+      a.GetOpeParams(res, smo, axa, axr, stepa, stepr, gpas, naxa);
+      stgpas = stepa;
+    }
+    else { // same mem-org for all, or same (this,a) OR same(this,res)
+      GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa);
+      GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa);
+      stgpas = step;
+    }
+    for(j=0; j<naxa; j++)  {
+      pe = mNDBlock.Begin()+Offset(ax,j);
+      pea = a.DataBlock().Begin()+a.Offset(axa,j);
+      per = res.DataBlock().Begin()+res.Offset(axr,j);
+      for(k=0; k<gpas; k+=stgpas, pe+=step, pea+=stepa, per+=stepr) *per = (*pe) * (*pea); 
+    }
+  }
+  return(res);
+}
+
+
+//! Two TArrays element by element division
+/*! 
+Perform element by element division of the source array (this) and the \b a array 
+and store the result in \b res (res = *this/a). The source and argument arrays (this, a)
+should have the same sizes.
+If not initially allocated, the output array \b res is automatically 
+resized as a packed array with the same sizes as the source (this) array.
+Returns a reference to the output array \b res. 
+\param a : Array to be added to the source array.
+\param res : Output array containing the result (res=*this/a). 
+\param fginv == true : Inverts the division argument order (res = a/(*this)) 
+\param divzero == true : Result is set to zero (res(i)=0) if the operation's 
+second argument is equal to zero (a(i)/(*this)(i)==0) 
+*/
+template <class T>
+TArray<T>& TArray<T>::DivElt(const TArray<T>& a, TArray<T>& res, bool fginv, bool divzero) const 
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::DivElt(...) - Not allocated source array ! ");
+  bool smoa;
+  if (!CompareSizes(a, smoa)) 
+    throw(SzMismatchError("TArray<T>::DivElt(...) SizeMismatch(this,a)")) ;
+  if (res.NbDimensions() < 1)  res.SetSize(*this, true, false);
+  bool smor;
+  if (!CompareSizes(res, smor)) 
+    throw(SzMismatchError("TArray<T>::DivElt(...) SizeMismatch(this,res) ")) ;
+
+  bool smora;
+  a.CompareSizes(res, smora);
+
+  bool smo = smoa && smor;  // The three arrays have same memory organisation
+
+  const T * pe;
+  const T * pea;
+  T * per;
+  sa_size_t j,k;
+  if (smo && IsPacked() && a.IsPacked() && res.IsPacked() )   {  // All packed arrays
+    sa_size_t maxx = totsize_;
+    pe = Data();
+    pea = a.Data();
+    per = res.Data();
+    if(divzero) {
+      if (!fginv) 
+        for(k=0; k<maxx;  k++, pe++, pea++, per++ )  
+          if (*pea==(T)0) *per = (T)0;  else   *per = *pe / *pea ;
+      else 
+        for(k=0; k<maxx;  k++, pe++, pea++, per++ )  
+          if (*pe==(T)0)  *per = (T)0;  else   *per = *pea / *pe ;
+    }
+    else {
+      if (!fginv) 
+        for(k=0; k<maxx;  k++, pe++, pea++, per++ )  *per = *pe / *pea ;
+      else
+        for(k=0; k<maxx;  k++, pe++, pea++, per++ )  *per = *pea / *pe ;
+    }
+  }
+  else {    // Non regular data spacing ...
+    int_4 ax,axa,axr;
+    sa_size_t step, stepa;
+    sa_size_t gpas, naxa;
+    sa_size_t stepr, stgpas;
+    if ( !smo && smora ) {   // same mem-org for a,res , different from this
+      a.GetOpeParams(*this, smo, axa, ax, stepa, step, gpas, naxa);
+      a.GetOpeParams(res, smo, axa, axr, stepa, stepr, gpas, naxa);
+      stgpas = stepa;
+      }
+    else { // same mem-org for all, or same (this,a) OR same(this,res)
+      GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa);
+      GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa);
+      stgpas = step;
+     }
+    // DBG cout << "DBG-A-DIVELT naxa=" << naxa << " gpas= " << gpas 
+    //	 << " step=" << step << " stepa=" << stepa << " stepr=" << stepr 
+    //	 << " ax= " << ax << " axa= " << axa << " axr= " << axr << endl;
+    for(j=0; j<naxa; j++)  {
+      pe = mNDBlock.Begin()+Offset(ax,j);
+      pea = a.DataBlock().Begin()+a.Offset(axa,j);
+      per = res.DataBlock().Begin()+res.Offset(axr,j);
+      if(divzero) {
+        if (!fginv) 
+          for(k=0; k<gpas; k+=stgpas, pe+=step, pea+=stepa, per+=stepr)  
+            if (*pea==(T)0) *per = (T)0;  else   *per = *pe / *pea ;
+	else
+          for(k=0; k<gpas; k+=stgpas, pe+=step, pea+=stepa, per+=stepr)  
+            if (*pe==(T)0)  *per = (T)0;  else   *per = *pea / *pe ;
+      }
+      else {
+        if (!fginv) 
+	  for(k=0; k<gpas; k+=stgpas, pe+=step, pea+=stepa, per+=stepr)  
+	    *per = *pe / *pea ;
+	else
+          for(k=0; k<gpas; k+=stgpas, pe+=step, pea+=stepa, per+=stepr)  
+	    *per = *pea / *pe ;
+      }
+    }
+  }
+  return(res);
+}
+
+
+//! 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>::AddElt(const TArray<T>&) SizeMismatch")) ;
+    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 j,k;
+  if (smo && (AvgStep() > 0) && (a.AvgStep() > 0) )   {  // regularly spaced elements
     sa_size_t step = AvgStep();
     sa_size_t stepa = a.AvgStep();
@@ -751,5 +1057,5 @@
     pe = Data();
     pea = a.Data();
-    for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] += pea[ka] ;
+    for(k=0;  k<maxx;  k+=step, pe+=step, pea+=stepa )  *pe = *pea ;
   }
   else {    // Non regular data spacing ...
@@ -761,186 +1067,5 @@
       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]; 
+      for(k=0;  k<gpas;  k+=step, pe+=step, pea+=stepa)  *pe = *pea; 
     }
   }
@@ -986,7 +1111,41 @@
 }
 
+//! Return the the scalar product of the two arrays (Sum_k[(*this)(k)*a(k)])
+template <class T>
+T TArray<T>::ScalarProduct(const TArray<T>& a)  const
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::ScalarProduct(...)  - Not allocated source array ");
+  bool smo;
+  if (!CompareSizes(a, smo)) 
+    throw(SzMismatchError("TArray<T>::ScalarProduct(...) SizeMismatch(this,a) ")) ;
+
+  T res = (T)(0);
+  const T * pe;
+  const T * pea;
+  sa_size_t j,k;
+  if (smo && (IsPacked() > 0) && (a.IsPacked() > 0))   {  // regularly spaced elements
+    sa_size_t maxx = totsize_;
+    pe = Data();
+    pea = a.Data();
+    for(k=0; k<maxx; k++, pe++, pea++)  res += (*pe)*(*pea);
+  }
+  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;  k<gpas; k+=step, pe+=step, pea+=stepa)  res += (*pe)*(*pea);
+    }
+  }
+  return(res);
+}
+
 
 // Somme et produit des elements
-//! Sum all elements
+//! Returns the sum of all array elements 
 template <class T>
 T TArray<T>::Sum() const
@@ -1016,5 +1175,5 @@
 }
 
-//! Multiply all elements
+//! Return the product of all elements
 template <class T>
 T TArray<T>::Product() const
@@ -1044,5 +1203,5 @@
 }
 
-//! Returns the sum of all elements squared (Sum
+//! Returns the sum of all array elements squared (Sum_k((*this)(k)*(*this)(k)).
 template <class T>
 T TArray<T>::SumX2() const
