Changeset 3332 in Sophya for trunk/SophyaLib/TArray/tarray.cc

Timestamp:

Oct 3, 2007, 3:04:34 PM (18 years ago)

Author:

ansari

Message:

1- Methode TArray::SumX2() renommee en ::SumSq()
2- Methode TArray::Norm2() appelle maintenant SumSq() - sauf pour les
tableaux complexes, ou on calcule Sum[el x conj(el)]=Sum[module^{2]
3- Nettoyage fichier sopemtx.cc (utilisation sa_size_t pour supprimer
les warnings g++}

Reza , 02/10/2007

File:

• trunk/SophyaLib/TArray/tarray.cc (modified) (3 diffs)

trunk/SophyaLib/TArray/tarray.cc

@@ -160 +160 @@
   string exmsg = "TArray<T>::TArray(int_4, sa_size_t *,  NDataBlock<T> & ... )";
   if (!UpdateSizes(ndim, siz, step, offset, exmsg))  throw( ParmError(exmsg) );
-  if (mNDBlock.Size() < ComputeTotalSize(ndim, siz, step, offset)) {
+  if (mNDBlock.Size() < (size_t)ComputeTotalSize(ndim, siz, step, offset)) {
     exmsg += " DataBlock.Size() < ComputeTotalSize(...) " ;
     throw( ParmError(exmsg) );
@@ -1395 +1395 @@
 //! Returns the sum of all array elements squared (Sum_k((*this)(k)*(*this)(k)).
 template <class T>
-T TArray<T>::SumX2() const
-{
-  if (NbDimensions() < 1) 
-    throw RangeCheckError("TArray<T>::SumX2()  - Not Allocated Array ! ");
+T TArray<T>::SumSq() const
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::SumSq()  - Not Allocated Array ! ");
   T ret=0;
   const T * pe;
@@ -1420 +1420 @@
   return ret;
 }
+
+/*!
+\brief Returns the array norm squared, defined as  Sum_k [ el(k)* el(k) ] 
+  For arrays with integer or real data, this method calls SumSq(), which computes
+  the sum of array elements squared. For complex arrays, it computes and returns 
+  the sum of array elements module squared (= Sum_k [el(k)*conj(el(k))]
+*/
+template <class T>
+T TArray<T>::Norm2() const
+{
+  return SumSq();
+}
+
+
+// Fonction auxiliaire pour specialisation de la methode Norm2() pour tableaux complexes
+template <class T> 
+complex<T>  _ComputeComplexNorm_Private_(TArray< complex<T> > const & ca) 
+{
+  if (ca.NbDimensions() < 1) 
+    throw RangeCheckError("TArray< complex<T> >::Norm2()  - Not Allocated Array ! ");
+  complex<T> ret= complex<T>(0., 0.);
+  const complex<T> * pe;
+  sa_size_t j,k;
+  if (ca.AvgStep() > 0)   {  // regularly spaced elements
+    sa_size_t step = ca.AvgStep();
+    sa_size_t maxx = ca.Size()*step;
+    pe = ca.Data();
+    for(k=0; k<maxx; k+=step )  ret += pe[k]*conj(pe[k]);
+  }
+  else {    // Non regular data spacing ...
+    int_4 ka = ca.MaxSizeKA();
+    sa_size_t step = ca.Step(ka);
+    sa_size_t gpas = ca.Size(ka)*step;
+    sa_size_t naxa = ca.Size()/ca.Size(ka);
+    for(j=0; j<naxa; j++)  {
+      pe = ca.DataBlock().Begin()+ca.Offset(ka,j);
+      for(k=0; k<gpas; k+=step)  ret += pe[k]*conj(pe[k]) ;
+    }
+  }
+  return ret;
+
+}
+
+// --- Specialisation de la methode Norm2() pour tableaux complexes ---
+DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
+complex<r_4> TArray< complex<r_4> >::Norm2() const
+{
+  return  _ComputeComplexNorm_Private_(*this);
+}
+DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
+complex<r_8> TArray< complex<r_8> >::Norm2() const
+{
+  return  _ComputeComplexNorm_Private_(*this);
+}
+//-------------------
 
 //! Return the minimum and the maximum values of the array elements