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


Ignore:
Timestamp:
Jul 29, 2004, 2:31:16 PM (21 years ago)
Author:
ansari
Message:

1/ Remplacement des methodes Add/Sub/Mul/DivElt(a) par

Add/Sub/Mul/DivElt(TArray a, TArray res)

2/ Operateurs += -= A+B A-B TArray et TMatrix/TVecteur modifies en consequence
3/ Ajout methode TArray::ScalarProduct()
4/ Methode TArray::SetT renomme en SetCst() SetT garde en alias
5/ Ajout parametre bool fzero (mise a zero) ajoute ds constructeur et

ReSize() de TMatrix et TVecteur.

Reza 29/07/2004

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaLib/TArray/tarray.cc

    r2569 r2575  
    482482//! Fill an array with a constant value \b x
    483483template <class T>
    484 TArray<T>& TArray<T>::SetT(T x)
    485 {
    486   if (NbDimensions() < 1)
    487     throw RangeCheckError("TArray<T>::SetT(T )  - Not Allocated Array ! ");
     484TArray<T>& TArray<T>::SetCst(T x)
     485{
     486  if (NbDimensions() < 1)
     487    throw RangeCheckError("TArray<T>::SetCst(T )  - Not Allocated Array ! ");
    488488  T * pe;
    489489  sa_size_t j,k;
    … …  
    528528  const T * pe;
    529529  T * per;
    530   sa_size_t j,k,kr;
     530  sa_size_t j,k;
    531531  if (smo && (IsPacked() > 0) && (res.IsPacked() > 0))   {  // regularly spaced elements
    532532    sa_size_t maxx = totsize_;
    … …  
    538538    int_4 ax,axr;
    539539    sa_size_t step, stepr;
    540     sa_size_t gpas, naxr;
    541     GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxr);
    542     for(j=0; j<naxr; j++)  {
     540    sa_size_t gpas, naxa;
     541    GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa);
     542    for(j=0; j<naxa; j++)  {
    543543      pe = mNDBlock.Begin()+Offset(ax,j);
    544544      per = res.DataBlock().Begin()+res.Offset(axr,j);
    545       for(k=0, kr=0;  k<gpas;  k+=step, kr+=stepr)  per[kr] = pe[k]+x;
     545      for(k=0;  k<gpas;  k+=step, pe+=step, per+=stepr)  *per = *pe+x;
    546546    }
    547547  }
    … …  
    557557\param x : constant to subtract from the array elements
    558558\param res : Output array containing the result (res=this+x or res=x-this).
    559 \param fginv == true : Invert subtraction argument order (*this = x-(*this))
     559\param fginv == true : Invert subtraction argument order (res = x-(*this))
    560560*/
    561561template <class T>
    … …  
    571571  const T * pe;
    572572  T * per;
    573   sa_size_t j,k,kr;
     573  sa_size_t j,k;
    574574  if (smo && (IsPacked() > 0) && (res.IsPacked() > 0))   {  // regularly spaced elements
    575575    sa_size_t maxx = totsize_;
    … …  
    584584    int_4 ax,axr;
    585585    sa_size_t step, stepr;
    586     sa_size_t gpas, naxr;
    587     GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxr);
    588     for(j=0; j<naxr; j++)  {
     586    sa_size_t gpas, naxa;
     587    GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa);
     588    for(j=0; j<naxa; j++)  {
    589589      pe = mNDBlock.Begin()+Offset(ax,j);
    590590      per = res.DataBlock().Begin()+res.Offset(axr,j);
    591591      if (!fginv)
    592         for(k=0, kr=0;  k<gpas;  k+=step, kr+=stepr)  per[kr] = pe[k]-x;
    593       else
    594         for(k=0, kr=0;  k<gpas;  k+=step, kr+=stepr)  per[kr] = x-pe[k];
     592        for(k=0;  k<gpas;  k+=step, pe+=step, per+=stepr)  *per = *pe-x;
     593       else
     594        for(k=0;  k<gpas;  k+=step, pe+=step, per+=stepr)  *per = x-*pe;
    595595    }
    596596  }
    … …  
    619619  const T * pe;
    620620  T * per;
    621   sa_size_t j,k,kr;
     621  sa_size_t j,k;
    622622  if (smo && (IsPacked() > 0) && (res.IsPacked() > 0))   {  // regularly spaced elements
    623623    sa_size_t maxx = totsize_;
    … …  
    629629    int_4 ax,axr;
    630630    sa_size_t step, stepr;
    631     sa_size_t gpas, naxr;
    632     GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxr);
    633     for(j=0; j<naxr; j++)  {
     631    sa_size_t gpas, naxa;
     632    GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa);
     633    for(j=0; j<naxa; j++)  {
    634634      pe = mNDBlock.Begin()+Offset(ax,j);
    635635      per = res.DataBlock().Begin()+res.Offset(axr,j);
    636       for(k=0, kr=0;  k<gpas;  k+=step, kr+=stepr)  per[kr] = pe[k]*x;
     636      for(k=0;  k<gpas;  k+=step, pe+=step, per+=stepr)  *per = (*pe)*x;
    637637    }
    638638  }
    … …  
    648648\param x : Array elements are divied by x
    649649\param res : Output array containing the result (res=(*this)/x or res=x/(*this)).
    650 \param fginv == true : Invert the operation order (res = x/(*this))
     650\param fginv == true : Invert the division argument order (res = x/(*this))
    651651*/
    652652template <class T>
    … …  
    677677    int_4 ax,axr;
    678678    sa_size_t step, stepr;
    679     sa_size_t gpas, naxr;
    680     GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxr);
    681     for(j=0; j<naxr; j++)  {
     679    sa_size_t gpas, naxa;
     680    GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa);
     681    for(j=0; j<naxa; j++)  {
    682682      pe = mNDBlock.Begin()+Offset(ax,j);
    683683      per = res.DataBlock().Begin()+res.Offset(axr,j);
    684684      if (!fginv)
    685         for(k=0, kr=0;  k<gpas;  k+=step, kr+=stepr)  per[kr] = pe[k]/x;
    686       else
    687         for(k=0, kr=0;  k<gpas;  k+=step, kr+=stepr)  per[kr] = x/pe[k];
     685        for(k=0;  k<gpas;  k+=step, pe+=step, per+=stepr)  *per = (*pe)/x;
     686       else
     687        for(k=0;  k<gpas;  k+=step, pe+=step, per+=stepr)  *per = x/(*pe);
    688688    }
    689689  }
    … …  
    710710  const T * pe;
    711711  T * per;
    712   sa_size_t j,k,kr;
     712  sa_size_t j,k;
    713713  if (smo && (IsPacked() > 0) && (res.IsPacked() > 0))   {  // regularly spaced elements
    714714    sa_size_t maxx = totsize_;
    … …  
    720720    int_4 ax,axr;
    721721    sa_size_t step, stepr;
    722     sa_size_t gpas, naxr;
    723     GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxr);
    724     for(j=0; j<naxr; j++)  {
     722    sa_size_t gpas, naxa;
     723    GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa);
     724    for(j=0; j<naxa; j++)  {
    725725      pe = mNDBlock.Begin()+Offset(ax,j);
    726726      per = res.DataBlock().Begin()+res.Offset(axr,j);
    727       for(k=0, kr=0;  k<gpas;  k+=step, kr+=stepr)  per[kr] = -pe[k];
     727      for(k=0;  k<gpas;  k+=step, pe+=step, per+=stepr)  *per = -(*pe);
    728728    }
    729729  }
    … …  
    732732
    733733//  >>>> Operations avec 2nd membre de type tableau
    734 //! Add two TArrays
    735 template <class T>
    736 TArray<T>& TArray<T>::AddElt(const TArray<T>& a)
    737 {
    738   if (NbDimensions() < 1)
    739     throw RangeCheckError("TArray<T>::AddElt(const TArray<T>& )  - Not Allocated Array ! ");
     734
     735//! Two TArrays element by element addition
     736/*!
     737Perform element by element addition of the source array (this) and the \b a array
     738and store the result in \b res (res = *this+a). The source and argument arrays (this, a)
     739should have the same sizes.
     740If not initially allocated, the output array \b res is automatically
     741resized as a packed array with the same sizes as the source (this) array.
     742Returns a reference to the output array \b res.
     743\param a : Array to be added to the source array.
     744\param res : Output array containing the result (res=this+a).
     745*/
     746template <class T>
     747TArray<T>& TArray<T>::AddElt(const TArray<T>& a, TArray<T>& res) const
     748{
     749  if (NbDimensions() < 1)
     750    throw RangeCheckError("TArray<T>::AddElt(...) - Not allocated source array ! ");
     751  bool smoa;
     752  if (!CompareSizes(a, smoa))
     753    throw(SzMismatchError("TArray<T>::AddElt(...) SizeMismatch(this,a)")) ;
     754  if (res.NbDimensions() < 1)  res.SetSize(*this, true, false);
     755  bool smor;
     756  if (!CompareSizes(res, smor))
     757    throw(SzMismatchError("TArray<T>::AddElt(...) SizeMismatch(this,res) ")) ;
     758
     759  bool smora;
     760  a.CompareSizes(res, smora);
     761
     762  bool smo = smoa && smor;  // The three arrays have same memory organisation
     763
     764  const T * pe;
     765  const T * pea;
     766  T * per;
     767  sa_size_t j,k;
     768  if (smo && IsPacked() && a.IsPacked() && res.IsPacked() )   {  // All packed arrays
     769    sa_size_t maxx = totsize_;
     770    pe = Data();
     771    pea = a.Data();
     772    per = res.Data();
     773    for(k=0; k<maxx;  k++, pe++, pea++, per++)  *per = *pe + *pea ;
     774  }
     775  else {    // Non regular data spacing ...
     776    int_4 ax,axa,axr;
     777    sa_size_t step, stepa;
     778    sa_size_t gpas, naxa;
     779    sa_size_t stepr, stgpas;
     780    if ( !smo && smora ) {   // same mem-org for a,res , different from this
     781      a.GetOpeParams(*this, smo, axa, ax, stepa, step, gpas, naxa);
     782      a.GetOpeParams(res, smo, axa, axr, stepa, stepr, gpas, naxa);
     783      stgpas = stepa;
     784    }
     785    else { // same mem-org for all, or same (this,a) OR same(this,res)
     786      GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa);
     787      GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa);
     788      stgpas = step;
     789    }
     790    for(j=0; j<naxa; j++)  {
     791      pe = mNDBlock.Begin()+Offset(ax,j);
     792      pea = a.DataBlock().Begin()+a.Offset(axa,j);
     793      per = res.DataBlock().Begin()+res.Offset(axr,j);
     794      for(k=0; k<gpas; k+=stgpas, pe+=step, pea+=stepa, per+=stepr) *per = *pe + *pea ;
     795    }
     796  }
     797  return(res);
     798}
     799
     800//! Two TArrays element by element subtraction
     801/*!
     802Perform element by element subtraction of the source array (this) and the \b a array
     803and the store result in \b res (res = *this-a or res=a-(*this)).
     804The source and argument arrays (this, a) should have the same sizes.
     805If not initially allocated, the output array \b res is automatically
     806resized as a packed array with the same sizes as the source (this) array.
     807Returns a reference to the output array \b res.
     808\param a : Array to be added to the source array.
     809\param res : Output array containing the result (res=*this+x).
     810\param fginv == true : Invert subtraction argument order (res = a-(*this)) 
     811*/
     812
     813template <class T>
     814TArray<T>& TArray<T>::SubElt(const TArray<T>& a, TArray<T>& res, bool fginv) const
     815{
     816  if (NbDimensions() < 1)
     817    throw RangeCheckError("TArray<T>::SubElt(...) - Not allocated source array ! ");
     818  bool smoa;
     819  if (!CompareSizes(a, smoa))
     820    throw(SzMismatchError("TArray<T>::SubElt(...) SizeMismatch(this,a)")) ;
     821  if (res.NbDimensions() < 1)  res.SetSize(*this, true, false);
     822  bool smor;
     823  if (!CompareSizes(res, smor))
     824    throw(SzMismatchError("TArray<T>::SubElt(...) SizeMismatch(this,res) ")) ;
     825
     826  bool smora;
     827  a.CompareSizes(res, smora);
     828
     829  bool smo = smoa && smor;  // The three arrays have same memory organisation
     830
     831  const T * pe;
     832  const T * pea;
     833  T * per;
     834  sa_size_t j,k;
     835  if (smo && IsPacked() && a.IsPacked() && res.IsPacked() )   {  // All packed arrays
     836    sa_size_t maxx = totsize_;
     837    pe = Data();
     838    pea = a.Data();
     839    per = res.Data();
     840    if (!fginv)
     841      for(k=0; k<maxx;  k++, pe++, pea++, per++)  *per = *pe - *pea ;
     842    else
     843      for(k=0; k<maxx;  k++, pe++, pea++, per++)  *per = *pea - *pe ;
     844  }
     845  else {    // Non regular data spacing ...
     846    int_4 ax,axa,axr;
     847    sa_size_t step, stepa;
     848    sa_size_t gpas, naxa;
     849    sa_size_t stepr, stgpas;
     850    if ( !smo && smora ) {   // same mem-org for a,res , different from this
     851      a.GetOpeParams(*this, smo, axa, ax, stepa, step, gpas, naxa);
     852      a.GetOpeParams(res, smo, axa, axr, stepa, stepr, gpas, naxa);
     853      stgpas = stepa;
     854    }
     855    else { // same mem-org for all, or same (this,a) OR same(this,res)
     856      GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa);
     857      GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa);
     858      stgpas = step;
     859    }
     860    for(j=0; j<naxa; j++)  {
     861      pe = mNDBlock.Begin()+Offset(ax,j);
     862      pea = a.DataBlock().Begin()+a.Offset(axa,j);
     863      per = res.DataBlock().Begin()+res.Offset(axr,j);
     864      if (!fginv)
     865        for(k=0; k<gpas; k+=stgpas, pe+=step, pea+=stepa, per+=stepr) *per = *pe - *pea ;
     866      else
     867        for(k=0; k<gpas; k+=stgpas, pe+=step, pea+=stepa, per+=stepr) *per = *pea - *pea ;
     868    }
     869  }
     870  return(res);
     871}
     872
     873
     874//! Two TArrays element by element multiplication
     875/*!
     876Perform element by element multiplication of the source array (this) and the \b a array
     877and store the result in \b res (res = *this*a). The source and argument arrays (this, a)
     878should have the same sizes.
     879If not initially allocated, the output array \b res is automatically
     880resized as a packed array with the same sizes as the source (this) array.
     881Returns a reference to the output array \b res.
     882\param a : Array to be added to the source array.
     883\param res : Output array containing the result (res=(*this)*a).
     884*/
     885template <class T>
     886TArray<T>& TArray<T>::MulElt(const TArray<T>& a, TArray<T>& res) const
     887{
     888  if (NbDimensions() < 1)
     889    throw RangeCheckError("TArray<T>::MulElt(...) - Not allocated source array ! ");
     890  bool smoa;
     891  if (!CompareSizes(a, smoa))
     892    throw(SzMismatchError("TArray<T>::MulElt(...) SizeMismatch(this,a)")) ;
     893  if (res.NbDimensions() < 1)  res.SetSize(*this, true, false);
     894  bool smor;
     895  if (!CompareSizes(res, smor))
     896    throw(SzMismatchError("TArray<T>::MulElt(...) SizeMismatch(this,res) ")) ;
     897
     898  bool smora;
     899  a.CompareSizes(res, smora);
     900
     901  bool smo = smoa && smor;  // The three arrays have same memory organisation
     902
     903  const T * pe;
     904  const T * pea;
     905  T * per;
     906  sa_size_t j,k;
     907  if (smo && IsPacked() && a.IsPacked() && res.IsPacked() )   {  // All packed arrays
     908    sa_size_t maxx = totsize_;
     909    pe = Data();
     910    pea = a.Data();
     911    per = res.Data();
     912    for(k=0; k<maxx;  k++, pe++, pea++, per++ )  *per = (*pe) * (*pea) ;
     913  }
     914  else {    // Non regular data spacing ...
     915    int_4 ax,axa,axr;
     916    sa_size_t step, stepa;
     917    sa_size_t gpas, naxa;
     918    sa_size_t stepr, stgpas;
     919    if ( !smo && smora ) {   // same mem-org for a,res , different from this
     920      a.GetOpeParams(*this, smo, axa, ax, stepa, step, gpas, naxa);
     921      a.GetOpeParams(res, smo, axa, axr, stepa, stepr, gpas, naxa);
     922      stgpas = stepa;
     923    }
     924    else { // same mem-org for all, or same (this,a) OR same(this,res)
     925      GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa);
     926      GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa);
     927      stgpas = step;
     928    }
     929    for(j=0; j<naxa; j++)  {
     930      pe = mNDBlock.Begin()+Offset(ax,j);
     931      pea = a.DataBlock().Begin()+a.Offset(axa,j);
     932      per = res.DataBlock().Begin()+res.Offset(axr,j);
     933      for(k=0; k<gpas; k+=stgpas, pe+=step, pea+=stepa, per+=stepr) *per = (*pe) * (*pea);
     934    }
     935  }
     936  return(res);
     937}
     938
     939
     940//! Two TArrays element by element division
     941/*!
     942Perform element by element division of the source array (this) and the \b a array
     943and store the result in \b res (res = *this/a). The source and argument arrays (this, a)
     944should have the same sizes.
     945If not initially allocated, the output array \b res is automatically
     946resized as a packed array with the same sizes as the source (this) array.
     947Returns a reference to the output array \b res.
     948\param a : Array to be added to the source array.
     949\param res : Output array containing the result (res=*this/a).
     950\param fginv == true : Inverts the division argument order (res = a/(*this))
     951\param divzero == true : Result is set to zero (res(i)=0) if the operation's
     952second argument is equal to zero (a(i)/(*this)(i)==0)
     953*/
     954template <class T>
     955TArray<T>& TArray<T>::DivElt(const TArray<T>& a, TArray<T>& res, bool fginv, bool divzero) const
     956{
     957  if (NbDimensions() < 1)
     958    throw RangeCheckError("TArray<T>::DivElt(...) - Not allocated source array ! ");
     959  bool smoa;
     960  if (!CompareSizes(a, smoa))
     961    throw(SzMismatchError("TArray<T>::DivElt(...) SizeMismatch(this,a)")) ;
     962  if (res.NbDimensions() < 1)  res.SetSize(*this, true, false);
     963  bool smor;
     964  if (!CompareSizes(res, smor))
     965    throw(SzMismatchError("TArray<T>::DivElt(...) SizeMismatch(this,res) ")) ;
     966
     967  bool smora;
     968  a.CompareSizes(res, smora);
     969
     970  bool smo = smoa && smor;  // The three arrays have same memory organisation
     971
     972  const T * pe;
     973  const T * pea;
     974  T * per;
     975  sa_size_t j,k;
     976  if (smo && IsPacked() && a.IsPacked() && res.IsPacked() )   {  // All packed arrays
     977    sa_size_t maxx = totsize_;
     978    pe = Data();
     979    pea = a.Data();
     980    per = res.Data();
     981    if(divzero) {
     982      if (!fginv)
     983        for(k=0; k<maxx;  k++, pe++, pea++, per++ ) 
     984          if (*pea==(T)0) *per = (T)0;  else   *per = *pe / *pea ;
     985      else
     986        for(k=0; k<maxx;  k++, pe++, pea++, per++ ) 
     987          if (*pe==(T)0)  *per = (T)0;  else   *per = *pea / *pe ;
     988    }
     989    else {
     990      if (!fginv)
     991        for(k=0; k<maxx;  k++, pe++, pea++, per++ )  *per = *pe / *pea ;
     992      else
     993        for(k=0; k<maxx;  k++, pe++, pea++, per++ )  *per = *pea / *pe ;
     994    }
     995  }
     996  else {    // Non regular data spacing ...
     997    int_4 ax,axa,axr;
     998    sa_size_t step, stepa;
     999    sa_size_t gpas, naxa;
     1000    sa_size_t stepr, stgpas;
     1001    if ( !smo && smora ) {   // same mem-org for a,res , different from this
     1002      a.GetOpeParams(*this, smo, axa, ax, stepa, step, gpas, naxa);
     1003      a.GetOpeParams(res, smo, axa, axr, stepa, stepr, gpas, naxa);
     1004      stgpas = stepa;
     1005      }
     1006    else { // same mem-org for all, or same (this,a) OR same(this,res)
     1007      GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa);
     1008      GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa);
     1009      stgpas = step;
     1010     }
     1011    // DBG cout << "DBG-A-DIVELT naxa=" << naxa << " gpas= " << gpas
     1012    //   << " step=" << step << " stepa=" << stepa << " stepr=" << stepr
     1013    //   << " ax= " << ax << " axa= " << axa << " axr= " << axr << endl;
     1014    for(j=0; j<naxa; j++)  {
     1015      pe = mNDBlock.Begin()+Offset(ax,j);
     1016      pea = a.DataBlock().Begin()+a.Offset(axa,j);
     1017      per = res.DataBlock().Begin()+res.Offset(axr,j);
     1018      if(divzero) {
     1019        if (!fginv)
     1020          for(k=0; k<gpas; k+=stgpas, pe+=step, pea+=stepa, per+=stepr) 
     1021            if (*pea==(T)0) *per = (T)0;  else   *per = *pe / *pea ;
     1022        else
     1023          for(k=0; k<gpas; k+=stgpas, pe+=step, pea+=stepa, per+=stepr) 
     1024            if (*pe==(T)0)  *per = (T)0;  else   *per = *pea / *pe ;
     1025      }
     1026      else {
     1027        if (!fginv)
     1028          for(k=0; k<gpas; k+=stgpas, pe+=step, pea+=stepa, per+=stepr) 
     1029            *per = *pe / *pea ;
     1030        else
     1031          for(k=0; k<gpas; k+=stgpas, pe+=step, pea+=stepa, per+=stepr) 
     1032            *per = *pea / *pe ;
     1033      }
     1034    }
     1035  }
     1036  return(res);
     1037}
     1038
     1039
     1040//! Copy elements of \b a
     1041template <class T>
     1042TArray<T>& TArray<T>::CopyElt(const TArray<T>& a)
     1043{
     1044  if (NbDimensions() < 1)
     1045    throw RangeCheckError("TArray<T>::CopyElt(const TArray<T>& )  - Not Allocated Array ! ");
    7401046  bool smo;
    7411047  if (!CompareSizes(a, smo))
    742     throw(SzMismatchError("TArray<T>::AddElt(const TArray<T>&) SizeMismatch")) ;
     1048    throw(SzMismatchError("TArray<T>::CopyElt(const TArray<T>&) SizeMismatch")) ;
    7431049
    7441050  T * pe;
    7451051  const T * pea;
    746   sa_size_t j,k,ka;
    747   if (smo && (AvgStep() > 0) && (a.AvgStep() > 0))   {  // regularly spaced elements
     1052  sa_size_t j,k;
     1053  if (smo && (AvgStep() > 0) && (a.AvgStep() > 0) )   {  // regularly spaced elements
    7481054    sa_size_t step = AvgStep();
    7491055    sa_size_t stepa = a.AvgStep();
    … …  
    7511057    pe = Data();
    7521058    pea = a.Data();
    753     for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] += pea[ka] ;
     1059    for(k=0;  k<maxx;  k+=step, pe+=step, pea+=stepa )  *pe = *pea ;
    7541060  }
    7551061  else {    // Non regular data spacing ...
    … …  
    7611067      pe = mNDBlock.Begin()+Offset(ax,j);
    7621068      pea = a.DataBlock().Begin()+a.Offset(axa,j);
    763       for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] += pea[ka];
    764     }
    765   }
    766   return(*this);
    767 }
    768 
    769 //! Substract two TArrays
    770 /*!
    771 Substract two TArrays *this = *this-a
    772 \param fginv == true : Perfoms the inverse subtraction (*this = a-(*this))
    773 */
    774 template <class T>
    775 TArray<T>& TArray<T>::SubElt(const TArray<T>& a, bool fginv)
    776 {
    777   if (NbDimensions() < 1)
    778     throw RangeCheckError("TArray<T>::SubElt(const TArray<T>& )  - Not Allocated Array ! ");
    779   bool smo;
    780   if (!CompareSizes(a, smo))
    781     throw(SzMismatchError("TArray<T>::SubElt(const TArray<T>&) SizeMismatch")) ;
    782 
    783   T * pe;
    784   const T * pea;
    785   sa_size_t j,k,ka;
    786   if (smo && (AvgStep() > 0) && (a.AvgStep() > 0) )   {  // regularly spaced elements
    787     sa_size_t step = AvgStep();
    788     sa_size_t stepa = a.AvgStep();
    789     sa_size_t maxx = totsize_*step;
    790     pe = Data();
    791     pea = a.Data();
    792     if (fginv)
    793       for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] = pea[ka]-pe[k] ;
    794     else
    795       for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] -= pea[ka] ;
    796   }
    797   else {    // Non regular data spacing ...
    798     int_4 ax,axa;
    799     sa_size_t step, stepa;
    800     sa_size_t gpas, naxa;
    801     GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa);
    802     for(j=0; j<naxa; j++)  {
    803       pe = mNDBlock.Begin()+Offset(ax,j);
    804       pea = a.DataBlock().Begin()+a.Offset(axa,j);
    805       if (fginv)
    806         for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] = pea[ka]-pe[k] ;
    807       else
    808         for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] -= pea[ka];
    809     }
    810   }
    811   return(*this);
    812 }
    813 
    814 
    815 //! Multiply two TArrays (elements by elements)
    816 template <class T>
    817 TArray<T>& TArray<T>::MulElt(const TArray<T>& a)
    818 {
    819   if (NbDimensions() < 1)
    820     throw RangeCheckError("TArray<T>::MulElt(const TArray<T>& )  - Not Allocated Array ! ");
    821   bool smo;
    822   if (!CompareSizes(a, smo))
    823     throw(SzMismatchError("TArray<T>::MulElt(const TArray<T>&) SizeMismatch")) ;
    824 
    825   T * pe;
    826   const T * pea;
    827   sa_size_t j,k,ka;
    828   if (smo && (AvgStep() > 0) && (a.AvgStep() > 0) )   {  // regularly spaced elements
    829     sa_size_t step = AvgStep();
    830     sa_size_t stepa = a.AvgStep();
    831     sa_size_t maxx = totsize_*step;
    832     pe = Data();
    833     pea = a.Data();
    834     for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] *= pea[ka] ;
    835   }
    836   else {    // Non regular data spacing ...
    837     int_4 ax,axa;
    838     sa_size_t step, stepa;
    839     sa_size_t gpas, naxa;
    840     GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa);
    841     for(j=0; j<naxa; j++)  {
    842       pe = mNDBlock.Begin()+Offset(axa,j);
    843       pea = a.DataBlock().Begin()+a.Offset(axa,j);
    844       for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] *= pea[ka];
    845     }
    846   }
    847   return(*this);
    848 }
    849 
    850 
    851 //! Divide two TArrays (elements by elements)
    852 /*!
    853 Divide two TArrays *this = (*this)/a
    854 \param fginv == true : Perfoms the inverse division (*this = a/(*this))
    855 \param divzero == true : if a(i)==0. result is set to zero: (*this)(i)==0.
    856 */
    857 template <class T>
    858 TArray<T>& TArray<T>::DivElt(const TArray<T>& a, bool fginv, bool divzero)
    859 {
    860   if (NbDimensions() < 1)
    861     throw RangeCheckError("TArray<T>::DivElt(const TArray<T>& )  - Not Allocated Array ! ");
    862   bool smo;
    863   if (!CompareSizes(a, smo))
    864     throw(SzMismatchError("TArray<T>::DivElt(const TArray<T>&) SizeMismatch")) ;
    865 
    866   T * pe;
    867   const T * pea;
    868   sa_size_t j,k,ka;
    869   if (smo && (AvgStep() > 0) && (a.AvgStep() > 0) )   {  // regularly spaced elements
    870     sa_size_t step = AvgStep();
    871     sa_size_t stepa = a.AvgStep();
    872     sa_size_t maxx = totsize_*step;
    873     pe = Data();
    874     pea = a.Data();
    875     if(divzero) {
    876       if (fginv)
    877         for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )
    878           {if(pe[k]==(T)0) pe[k] = (T)0; else pe[k] = pea[ka]/pe[k];}
    879       else
    880         for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )
    881           {if(pea[k]==(T)0) pe[k] = (T)0; else pe[k] /= pea[ka] ;}
    882     } else {
    883       if (fginv)
    884         for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] = pea[ka]/pe[k];
    885       else
    886         for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] /= pea[ka] ;
    887     }
    888   }
    889   else {    // Non regular data spacing ...
    890     int_4 ax,axa;
    891     sa_size_t step, stepa;
    892     sa_size_t gpas, naxa;
    893     GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa);
    894     for(j=0; j<naxa; j++)  {
    895       pe = mNDBlock.Begin()+Offset(ax,j);
    896       pea = a.DataBlock().Begin()+a.Offset(axa,j);
    897       if(divzero) {
    898         if (fginv)
    899           for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)
    900             {if(pe[k]==(T)0) pe[k] = (T)0; else pe[k] = pea[ka]/pe[k];}
    901         else
    902           for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)
    903             {if(pea[k]==(T)0) pe[k] = (T)0; else pe[k] /= pea[ka];}
    904       } else {
    905         if (fginv)
    906           for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] = pea[ka]/pe[k];
    907         else
    908           for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] /= pea[ka];
    909       }
    910     }
    911   }
    912   return(*this);
    913 }
    914 
    915 //! Copy elements of \b a
    916 template <class T>
    917 TArray<T>& TArray<T>::CopyElt(const TArray<T>& a)
    918 {
    919   if (NbDimensions() < 1)
    920     throw RangeCheckError("TArray<T>::CopyElt(const TArray<T>& )  - Not Allocated Array ! ");
    921   bool smo;
    922   if (!CompareSizes(a, smo))
    923     throw(SzMismatchError("TArray<T>::CopyElt(const TArray<T>&) SizeMismatch")) ;
    924 
    925   T * pe;
    926   const T * pea;
    927   sa_size_t j,k,ka;
    928   if (smo && (AvgStep() > 0) && (a.AvgStep() > 0) )   {  // regularly spaced elements
    929     sa_size_t step = AvgStep();
    930     sa_size_t stepa = a.AvgStep();
    931     sa_size_t maxx = totsize_*step;
    932     pe = Data();
    933     pea = a.Data();
    934     for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] = pea[ka] ;
    935   }
    936   else {    // Non regular data spacing ...
    937     int_4 ax,axa;
    938     sa_size_t step, stepa;
    939     sa_size_t gpas, naxa;
    940     GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa);
    941     for(j=0; j<naxa; j++)  {
    942       pe = mNDBlock.Begin()+Offset(ax,j);
    943       pea = a.DataBlock().Begin()+a.Offset(axa,j);
    944       for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] = pea[ka];
     1069      for(k=0;  k<gpas;  k+=step, pe+=step, pea+=stepa)  *pe = *pea;
    9451070    }
    9461071  }
    … …  
    9861111}
    9871112
     1113//! Return the the scalar product of the two arrays (Sum_k[(*this)(k)*a(k)])
     1114template <class T>
     1115T TArray<T>::ScalarProduct(const TArray<T>& a)  const
     1116{
     1117  if (NbDimensions() < 1)
     1118    throw RangeCheckError("TArray<T>::ScalarProduct(...)  - Not allocated source array ");
     1119  bool smo;
     1120  if (!CompareSizes(a, smo))
     1121    throw(SzMismatchError("TArray<T>::ScalarProduct(...) SizeMismatch(this,a) ")) ;
     1122
     1123  T res = (T)(0);
     1124  const T * pe;
     1125  const T * pea;
     1126  sa_size_t j,k;
     1127  if (smo && (IsPacked() > 0) && (a.IsPacked() > 0))   {  // regularly spaced elements
     1128    sa_size_t maxx = totsize_;
     1129    pe = Data();
     1130    pea = a.Data();
     1131    for(k=0; k<maxx; k++, pe++, pea++)  res += (*pe)*(*pea);
     1132  }
     1133  else {    // Non regular data spacing ...
     1134    int_4 ax,axa;
     1135    sa_size_t step, stepa;
     1136    sa_size_t gpas, naxa;
     1137    GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa);
     1138    for(j=0; j<naxa; j++)  {
     1139      pe = mNDBlock.Begin()+Offset(ax,j);
     1140      pea = a.DataBlock().Begin()+a.Offset(axa,j);
     1141      for(k=0;  k<gpas; k+=step, pe+=step, pea+=stepa)  res += (*pe)*(*pea);
     1142    }
     1143  }
     1144  return(res);
     1145}
     1146
    9881147
    9891148// Somme et produit des elements
    990 //! Sum all elements
     1149//! Returns the sum of all array elements
    9911150template <class T>
    9921151T TArray<T>::Sum() const
    … …  
    10161175}
    10171176
    1018 //! Multiply all elements
     1177//! Return the product of all elements
    10191178template <class T>
    10201179T TArray<T>::Product() const
    … …  
    10441203}
    10451204
    1046 //! Returns the sum of all elements squared (Sum
     1205//! Returns the sum of all array elements squared (Sum_k((*this)(k)*(*this)(k)).
    10471206template <class T>
    10481207T TArray<T>::SumX2() const
Note: See TracChangeset for help on using the changeset viewer.