Changeset 2575 in Sophya


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

Location:
trunk/SophyaLib/TArray
Files:
6 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

  • trunk/SophyaLib/TArray/tarray.h

    r2569 r2575  
    125125  //! Fill TArray with Sequence \b seq
    126126  inline  TArray<T>&  operator = (Sequence const & seq)    { return SetSeq(seq); }
     127
    127128// A = x (tous les elements a x)
    128   virtual TArray<T>&  SetT(T x);
     129  virtual TArray<T>&  SetCst(T x);
     130  //! Fill an array with a constant value \b x ( alias for \b SetCst() method )
     131  inline  TArray<T>&  SetT(T x) { return SetCst(x); }
    129132  //! Fill TArray with all elements equal to \b x
    130133  inline  TArray<T>&  operator = (T x)             { return SetT(x); }
    … …  
    138141
    139142// A += -= *= /= x (ajoute, soustrait, ... x a tous les elements)
     143  // Methodes Add/Sub/Mul/Div() sont la pour compatibilite avant V=2 (1.818)
     144  // Faut-il les garder ? Reza, Juillet 2004
    140145  inline  TArray<T>&  Add(T x)                     { return AddCst(x, *this); }
    141146  inline  TArray<T>&  Sub(T x, bool fginv=false)   { return SubCst(x, *this, fginv); }
    … …  
    158163
    159164// A += -=  (ajoute, soustrait element par element les deux tableaux )
    160   virtual TArray<T>&  AddElt(const TArray<T>& a);
    161   //! Operator TArray += TArray
    162   inline  TArray<T>&  operator += (const TArray<T>& a)  { return AddElt(a); }
    163   virtual TArray<T>&  SubElt(const TArray<T>& a, bool fginv=false);
    164   //! Operator TArray -= TArray
    165   inline  TArray<T>&  operator -= (const TArray<T>& a)  { return SubElt(a); }
     165  virtual TArray<T>&  AddElt(const TArray<T>& a, TArray<T>& res) const ;
     166  virtual TArray<T>&  SubElt(const TArray<T>& a, TArray<T>& res, bool fginv=false) const ;
    166167// Multiplication, division element par element les deux tableaux
    167   virtual TArray<T>&  MulElt(const TArray<T>& a);
    168   virtual TArray<T>&  DivElt(const TArray<T>& a, bool fginv=false, bool divzero=false);
     168  virtual TArray<T>&  MulElt(const TArray<T>& a, TArray<T>& res) const ;
     169  virtual TArray<T>&  DivElt(const TArray<T>& a, TArray<T>& res, bool fginv=false, bool divzero=false) const ;
     170
     171  //! Operator TArray += TArray (element by element addition in place)
     172  inline  TArray<T>&  operator += (const TArray<T>& a)  { return AddElt(a, *this); }
     173  //! Operator TArray -= TArray (element by element subtraction in place)
     174  inline  TArray<T>&  operator -= (const TArray<T>& a)  { return SubElt(a, *this); }
     175
    169176// Recopie des valeurs, element par element
    170177  virtual TArray<T>&  CopyElt(const TArray<T>& a);
    171178// Recopie des valeurs avec conversion prealable, element par element
    172179  virtual TArray<T>&  ConvertAndCopyElt(const BaseArray& a);
     180
     181// Calcul du produit scalaire ( Somme_i (*this)(i)*a(i) )
     182  virtual T ScalarProduct(const TArray<T>& a) const ;
     183// Norme(^2)
     184  //! Returns the squarred of the array norm, defined as  Sum_k (*this)(k)*(*this)(k)
     185  inline T Norm2() const { return ScalarProduct(*this); }
    173186
    174187// Somme et produit des elements
    … …  
    283296inline TArray<T> operator + (const TArray<T>& a,const TArray<T>& b)
    284297    { TArray<T> result; result.SetTemp(true);
    285     if (b.IsTemp())  { result.Share(b); result.AddElt(a); }
    286     else { result.CloneOrShare(a); result.AddElt(b); }
    287     return result; }
     298    a.AddElt(b, result);    return result; }
    288299
    289300/*! \ingroup TArray \fn operator-(const TArray<T>&,const TArray<T>&)
    … …  
    292303inline TArray<T> operator - (const TArray<T>& a,const TArray<T>& b)
    293304    { TArray<T> result; result.SetTemp(true);
    294     if (b.IsTemp())  { result.Share(b); result.SubElt(a, true); }
    295     else { result.CloneOrShare(a); result.SubElt(b); }
    296     return result; }
     305    a.SubElt(b, result);    return result; }
    297306
    298307

  • trunk/SophyaLib/TArray/tmatrix.cc

    r2574 r2575  
    1 // $Id: tmatrix.cc,v 1.25 2004-07-29 08:40:49 cmv Exp $
     1// $Id: tmatrix.cc,v 1.26 2004-07-29 12:31:16 ansari Exp $
    22//                         C.Magneville          04/99
    33#include "machdefs.h"
    … …  
    6161  \param c : number of columns
    6262  \param mm : define the memory mapping type
     63  \param fzero : if \b true , set matrix elements to zero
    6364  \sa ReSize
    6465 */
    6566template <class T>
    66 TMatrix<T>::TMatrix(sa_size_t r,sa_size_t c, short mm)
     67TMatrix<T>::TMatrix(sa_size_t r,sa_size_t c, short mm, bool fzero)
    6768// Construit une matrice de r lignes et c colonnes.
    6869  :  TArray<T>()
    … …  
    7172    throw ParmError("TMatrix<T>::TMatrix(sa_size_t r,sa_size_t c) NRows or NCols = 0");
    7273  arrtype_ = 1;   // Type = Matrix
    73   ReSize(r, c, mm);
     74  ReSize(r, c, mm, fzero);
    7475}
    7576
    … …  
    204205         (SameMemoryMapping,CMemoryMapping
    205206         ,FortranMemoryMapping,DefaultMemoryMapping)
    206  */
    207 template <class T>
    208 void TMatrix<T>::ReSize(sa_size_t r, sa_size_t c, short mm)
     207  \param fzero : if \b true , set matrix elements to zero
     208 */
     209template <class T>
     210void TMatrix<T>::ReSize(sa_size_t r, sa_size_t c, short mm, bool fzero)
    209211{
    210212  if(r==0||c==0)
    … …  
    223225    size[0] = r;  size[1] = c;
    224226  }
    225   TArray<T>::ReSize(2, size, 1);
     227  TArray<T>::ReSize(2, size, 1, fzero);
    226228  UpdateMemoryMapping(mm);
    227229}

  • trunk/SophyaLib/TArray/tmatrix.h

    r2564 r2575  
    1515  // Creation / destruction
    1616  TMatrix();
    17   TMatrix(sa_size_t r,sa_size_t c, short mm=BaseArray::AutoMemoryMapping);
     17  TMatrix(sa_size_t r,sa_size_t c, short mm=BaseArray::AutoMemoryMapping, bool fzero=true);
    1818  TMatrix(const TMatrix<T>& a);
    1919  TMatrix(const TMatrix<T>& a, bool share);
    … …  
    4848  inline sa_size_t NCol()  const {return size_[macoli_]; } // back-compat Peida
    4949
    50   void ReSize(sa_size_t r,sa_size_t c, short mm=BaseArray::SameMemoryMapping);
     50  void ReSize(sa_size_t r,sa_size_t c, short mm=BaseArray::SameMemoryMapping, bool fzero=true);
    5151  //! a synonym (alias) for method ReSize(sa_size_t, sa_size_t, short)
    52   inline void SetSize(sa_size_t r,sa_size_t c, short mm=BaseArray::SameMemoryMapping)
    53                  { ReSize(r, c, mm); }
     52  inline void SetSize(sa_size_t r,sa_size_t c, short mm=BaseArray::SameMemoryMapping, bool fzero=true)
     53                 { ReSize(r, c, mm, fzero); }
    5454  // Reallocation de place
    5555  void Realloc(sa_size_t r,sa_size_t c, short mm=BaseArray::SameMemoryMapping, bool force=false);
    … …  
    103103  //  operations avec matrices
    104104  //! += : add a matrix
    105   inline  TMatrix<T>&  operator += (const TMatrix<T>& a)  { AddElt(a); return(*this); }
     105  inline  TMatrix<T>&  operator += (const TMatrix<T>& a)  { AddElt(a,*this); return(*this); }
    106106  //! -= : substract a matrix
    107   inline  TMatrix<T>&  operator -= (const TMatrix<T>& a)  { SubElt(a); return(*this); }
     107  inline  TMatrix<T>&  operator -= (const TMatrix<T>& a)  { SubElt(a,*this); return(*this); }
     108
    108109  TMatrix<T>  Multiply(const TMatrix<T>& b, short mm=BaseArray::SameMemoryMapping) const;
    109   //! *= : matrix product : C = (*this)*B
    110   inline  TMatrix<T>&  operator *= (const TMatrix<T>& b)
    111           { this->Set(Multiply(b));  return(*this); }
     110  //A supprimer ? Reza Juillet 2004 ! *= : matrix product : C = (*this)*B
     111  //  inline  TMatrix<T>&  operator *= (const TMatrix<T>& b)
     112  //       { this->Set(Multiply(b));  return(*this); }
    112113
    113114  // I/O print, ...
    … …  
    201202  the original matrix elements.  */
    202203template <class T> inline TMatrix<T> operator - (const TMatrix<T>& a)
    203     {TMatrix<T> result; result.CloneOrShare(a); result.SetTemp(true);
    204     result.NegateElt(); return result;}
     204    {TMatrix<T> result; result.SetTemp(true);
     205    a.NegateElt(result); return result;}
    205206
    206207
    … …  
    215216inline TMatrix<T> operator + (const TMatrix<T>& a,const TMatrix<T>& b)
    216217    {TMatrix<T> result; result.SetTemp(true);
    217     if (b.IsTemp())  { result.Share(b); result.AddElt(a); }
    218     else { result.CloneOrShare(a); result.AddElt(b); }
    219     return result; }
     218    a.AddElt(b, result);     return result; }
    220219
    221220
    … …  
    225224inline TMatrix<T> operator - (const TMatrix<T>& a,const TMatrix<T>& b)
    226225    {TMatrix<T> result; result.SetTemp(true);
    227     if (b.IsTemp())  { result.Share(b); result.SubElt(a, true); }
    228     else { result.CloneOrShare(a); result.SubElt(b); }
    229     return result; }
     226    a.SubElt(b, result);     return result; }
     227   
    230228
    231229// Surcharge d'operateurs C = A * B

  • trunk/SophyaLib/TArray/tvector.cc

    r2267 r2575  
    1 // $Id: tvector.cc,v 1.15 2002-11-15 09:35:44 ansari Exp $
     1// $Id: tvector.cc,v 1.16 2004-07-29 12:31:16 ansari Exp $
    22//                         C.Magneville          04/99
    33#include "machdefs.h"
    … …  
    5252  \param lcv : line or column vector ?
    5353  \param mm : memory mapping type
     54  \param fzero : if \b true , set vector elements to zero
    5455  \sa SelectVectorType
    5556 */
    5657template <class T>
    57 TVector<T>::TVector(sa_size_t n, short lcv, short mm)
     58TVector<T>::TVector(sa_size_t n, short lcv, short mm, bool fzero)
    5859// Constructeur
    59   : TMatrix<T>(1,1,mm)
     60  : TMatrix<T>(1,1,mm,false)
    6061{
    6162  arrtype_ = 2;   // Type = Vector
    6263  lcv = SelectVectorType(lcv);
    63   ReSize(n,lcv);
     64  ReSize(n,lcv,fzero);
    6465}
    6566
    … …  
    145146 */
    146147template <class T>
    147 void TVector<T>::ReSize(sa_size_t n, short lcv)
     148void TVector<T>::ReSize(sa_size_t n, short lcv, bool fzero)
    148149{
    149150  if( n == 0 )
    … …  
    154155  if (lcv == ColumnVector) { r = n;  c = 1; }
    155156  else { c = n; r = 1; }
    156   TMatrix<T>::ReSize(r,c);
     157  TMatrix<T>::ReSize(r,c,BaseArray::SameMemoryMapping,fzero);
    157158  veceli_ = (lcv ==  ColumnVector ) ?  marowi_ : macoli_;
    158159}
    … …  
    191192  TVector sv( mtx(rr, cr) , true, GetVectorType(), GetMemoryMapping() );
    192193  return(sv); 
    193 }
    194 
    195 //! Return the norm \f$ \sum{V(i)^2} \f$
    196 template <class T>
    197 T TVector<T>::Norm2() const
    198 {
    199   T ret = 0;
    200   for(sa_size_t k=0; k<Size(); k++)   ret += (*this)(k)*(*this)(k);
    201   return ret;
    202194}
    203195

  • trunk/SophyaLib/TArray/tvector.h

    r2564 r2575  
    1414  // Creation / destruction
    1515  TVector();
    16   TVector(sa_size_t n, short lcv=AutoVectorType, short mm=AutoMemoryMapping);
     16  TVector(sa_size_t n, short lcv=AutoVectorType, short mm=AutoMemoryMapping, bool fzero=true);
    1717  TVector(const TVector<T>& v);
    1818  TVector(const TVector<T>& v, bool share);
    … …  
    3939   
    4040  // Gestion taille/Remplissage
    41   void ReSize(sa_size_t n, short lcv=SameVectorType );
     41  void ReSize(sa_size_t n, short lcv=SameVectorType, bool fzero=true);
    4242  //! a synonym (alias) for method ReSize(sa_size_t, short)
    43   inline void SetSize(sa_size_t n, short lcv=SameVectorType )
    44               { ReSize(n, lcv); }
     43  inline void SetSize(sa_size_t n, short lcv=SameVectorType, bool fzero=true)
     44              { ReSize(n, lcv, fzero); }
    4545  void Realloc(sa_size_t n, short lcv=SameVectorType, bool force=false);
    4646
    … …  
    7777  //  operations avec matrices
    7878  //! += : add a vector in place
    79   inline  TVector<T>&  operator += (const TVector<T>& a)  { AddElt(a); return(*this); }
     79  inline  TVector<T>&  operator += (const TVector<T>& a)  { AddElt(a,*this); return(*this); }
    8080  //! += : substract a vector in place
    81   inline  TVector<T>&  operator -= (const TVector<T>& a)  { SubElt(a); return(*this); }
    82 
    83   // Norme(^2)
    84   T Norm2() const ;
     81  inline  TVector<T>&  operator -= (const TVector<T>& a)  { SubElt(a,*this); return(*this); }
    8582
    8683  virtual string InfoString() const;
Note: See TracChangeset for help on using the changeset viewer.