Changeset 3332 in Sophya


Ignore:
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[module2]
3- Nettoyage fichier sopemtx.cc (utilisation sa_size_t pour supprimer
les warnings g++

Reza , 02/10/2007

Location:
trunk/SophyaLib/TArray
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaLib/TArray/sopemtx.cc

    r3012 r3332  
    3636  };
    3737  TMatrixRC();
    38   TMatrixRC(TMatrix<T>& m, TRCKind kind, uint_4 index=0);
     38  TMatrixRC(TMatrix<T>& m, TRCKind kind, sa_size_t index=0);
    3939  virtual ~TMatrixRC() {}
    4040
    4141  // Acces aux rangees et colonnes de matrices
    42   static TMatrixRC<T> Row(TMatrix<T> & m, uint_4 r);
    43   static TMatrixRC<T> Col(TMatrix<T> & m, uint_4 c);
     42  static TMatrixRC<T> Row(TMatrix<T> & m, sa_size_t r);
     43  static TMatrixRC<T> Col(TMatrix<T> & m, sa_size_t c);
    4444  static TMatrixRC<T> Diag(TMatrix<T> & m);
    4545
    46   // ---- A virer   $CHECK$ Reza 03/2000
    47   //  int_4 Next();
    48   //  int_4 Prev();
    4946  int_4 SetCol(int_4 c);
    5047  int_4 SetRow(int_4 r);
    5148  int_4 SetDiag();
    5249
    53   static uint_4 Step(const TMatrix<T>& m, TRCKind rckind);
    54   static T* Org(const TMatrix<T>&, TRCKind rckind, uint_4 ind=0);
     50  static sa_size_t Step(const TMatrix<T>& m, TRCKind rckind);
     51  static T* Org(const TMatrix<T>&, TRCKind rckind, sa_size_t ind=0);
    5552
    5653  //! Return the kind of TMatrix (line,column,diagonal)
    5754  TRCKind Kind() const { return kind; }
    58   uint_4 NElts() const;
    59   T& operator()(uint_4 i);
    60   T  operator()(uint_4 i) const;
     55  sa_size_t NElts() const;
     56  T& operator()(sa_size_t i);
     57  T  operator()(sa_size_t i) const;
    6158
    6259
    6360   TMatrixRC<T>& operator = (const TMatrixRC<T>& rc);
    6461
    65 // ---- A virer   $CHECK$ Reza 03/2000
    66 //   TVector<T> GetVect() const;
    67 //   TMatrixRC<T>& operator += (const TMatrixRC<T>& rc);
    68 //   TMatrixRC<T>& operator -= (const TMatrixRC<T>& rc);
    6962
    7063  TMatrixRC<T>& operator *= (T x);
    … …  
    7568 
    7669
    77   TMatrixRC<T>& LinComb(T a, T b, const TMatrixRC& rc, uint_4 first=0);
    78   TMatrixRC<T>& LinComb(T b, const TMatrixRC<T>& rc, uint_4 first=0);
    79 
    80   uint_4 IMaxAbs(uint_4 first=0);
     70  TMatrixRC<T>& LinComb(T a, T b, const TMatrixRC& rc, sa_size_t first=0);
     71  TMatrixRC<T>& LinComb(T b, const TMatrixRC<T>& rc, sa_size_t first=0);
     72
     73  sa_size_t IMaxAbs(sa_size_t first=0);
    8174  void Print(ostream & os) const ;
    8275
    … …  
    9386  //! Define Absolute value for int_8
    9487  inline static double Abs_Value(int_8 v)  {return (v>0)? (double) v: (double) -v;}
    95   //! Define Absolute value for uint_4
     88  //! Define Absolute value for sa_size_t
    9689  inline static double Abs_Value(uint_4 v) {return (double) v;}
    9790  //! Define Absolute value for uint_8
    … …  
    112105  T*          data;    //!< pointer to the beginnig of interesting datas
    113106  int_4       index;   //!< index of the line/column
    114   uint_4      step;    //!< step of the line/column
     107  sa_size_t      step;    //!< step of the line/column
    115108  TRCKind     kind;    //!< type: line, column or diagonal
    116109};
    … …  
    128121    throw(SzMismatchError("TMatrixRC::operator * type mismatch\n"));
    129122  T sum = 0;
    130   for(uint_4 i=0; i<a.NElts(); i++) sum += a(i)*b(i);
     123  for(sa_size_t i=0; i<a.NElts(); i++) sum += a(i)*b(i);
    131124  return sum;
    132125  }
    … …  
    138131 */
    139132template <class T>
    140 inline uint_4 TMatrixRC<T>::Step(const TMatrix<T>& m, TRCKind rckind)
     133inline sa_size_t TMatrixRC<T>::Step(const TMatrix<T>& m, TRCKind rckind)
    141134  { switch (rckind) { case TmatrixRow  : return m.Step(m.ColsKA());
    142135                      case TmatrixCol  : return m.Step(m.RowsKA());
    … …  
    153146 */
    154147template <class T>
    155 inline T* TMatrixRC<T>::Org(const TMatrix<T>& m, TRCKind rckind, uint_4 index)
     148inline T* TMatrixRC<T>::Org(const TMatrix<T>& m, TRCKind rckind, sa_size_t index)
    156149  { switch (rckind) { case TmatrixRow  : return const_cast<T *>(&(m(index,0)));
    157150                      case TmatrixCol  : return const_cast<T *>(&(m(0,index)));
    … …  
    160153
    161154//! return number of elements for a TMatrixRC
    162 template <class T> inline uint_4 TMatrixRC<T>::NElts() const
     155template <class T> inline sa_size_t TMatrixRC<T>::NElts() const
    163156  { if (!matrix) return 0;
    164157    switch (kind) { case TmatrixRow  : return matrix->NCols();
    … …  
    169162//! access of element \b i
    170163template <class T>
    171 inline T& TMatrixRC<T>::operator()(uint_4 i) {return data[i*step];}
     164inline T& TMatrixRC<T>::operator()(sa_size_t i) {return data[i*step];}
    172165//! access of element \b i
    173166template <class T>
    174 inline T  TMatrixRC<T>::operator()(uint_4 i) const {return data[i*step];}
     167inline T  TMatrixRC<T>::operator()(sa_size_t i) const {return data[i*step];}
    175168
    176169////////////////////////////////////////////////////////////////
    … …  
    190183  \param ind : number of the line or column
    191184*/
    192 template <class T> TMatrixRC<T>::TMatrixRC(TMatrix<T>& m,TRCKind rckind,uint_4 ind)
     185template <class T> TMatrixRC<T>::TMatrixRC(TMatrix<T>& m,TRCKind rckind,sa_size_t ind)
    193186: matrix(&m), data(Org(m,rckind,ind)),
    194187  index(ind), step(Step(m,rckind)), kind(rckind)
    … …  
    203196//! Return TMatrixRC for line \b r of matrix \b m
    204197template <class T>
    205 TMatrixRC<T> TMatrixRC<T>::Row(TMatrix<T> & m, uint_4 r)
     198TMatrixRC<T> TMatrixRC<T>::Row(TMatrix<T> & m, sa_size_t r)
    206199{
    207200TMatrixRC<T> rc(m, TmatrixRow, r);
    … …  
    211204//! Return TMatrixRC for column \b r of matrix \b m
    212205template <class T>
    213 TMatrixRC<T> TMatrixRC<T>::Col(TMatrix<T> & m, uint_4 c)
     206TMatrixRC<T> TMatrixRC<T>::Col(TMatrix<T> & m, sa_size_t c)
    214207{
    215208TMatrixRC<T> rc(m, TmatrixCol, c);
    … …  
    226219
    227220
    228 // ---- A virer   $CHECK$ Reza 03/2000
    229 // template <class T> int_4 TMatrixRC<T>::Next()
    230 // {
    231 // if (!matrix || kind==TmatrixDiag) return -1;
    232 // index++;
    233 // if(kind == TmatrixRow) {
    234 //   if(index > (int_4)matrix->NRows()) {index = (int_4)matrix->NRows(); return -1;}
    235 //   data += matrix->NCols();
    236 // } else {
    237 //   if (index > (int_4)matrix->NCols()) {index = (int_4)matrix->NCols(); return -1;}
    238 //   data++;
    239 // }
    240 // return index;
    241 // }
    242 
    243 // template <class T> int_4 TMatrixRC<T>::Prev()
    244 // {
    245 // if (!matrix || kind == TmatrixDiag) return -1;
    246 // index--;
    247 // if(index < 0) {index = 0; return -1;}
    248 // if(kind == TmatrixRow) data -= matrix->NCols();
    249 // else data--;
    250 // return index;
    251 // }
    252221
    253222//! Set column \b c for this TMatrixRC
    … …  
    299268}
    300269
    301 // ---- A virer   $CHECK$ Reza 03/2000
    302 // template <class T> TVector<T> TMatrixRC<T>::GetVect() const
    303 // {
    304 // TVector<T> v(NElts());
    305 // for (uint_4 i=0; i<NElts(); i++) v(i) = (*this)(i);
    306 // return v;
    307 // }
    308 
    309 // template <class T> TMatrixRC<T>& TMatrixRC<T>::operator += (const TMatrixRC<T>& rc)
    310 // {
    311 // if ( NElts() != rc.NElts() )
    312 //   throw(SzMismatchError("TMatrixRC::operator+= size mismatch\n"));
    313 // if ( kind != rc.kind )
    314 //   throw(SzMismatchError("TMatrixRC::operator+= type mismatch\n"));
    315 // for (uint_4 i=0; i<NElts(); i++) (*this)(i) += rc(i);
    316 // return *this;
    317 // }
    318 
    319 // template <class T> TMatrixRC<T>& TMatrixRC<T>::operator -= (const TMatrixRC<T>& rc)
    320 // {
    321 // if( NElts() != rc.NElts() )
    322 //   throw(SzMismatchError("TMatrixRC::operator-= size mismatch\n"));
    323 // if( kind != rc.kind )
    324 //   throw(SzMismatchError("TMatrixRC::operator-= type mismatch\n"));
    325 // for(uint_4 i=0; i<NElts(); i++) (*this)(i) -= rc(i);
    326 // return *this;
    327 // }
    328270
    329271//! Operator to multiply by constant \b x
    330272template <class T> TMatrixRC<T>& TMatrixRC<T>::operator *= (T x)
    331273{
    332 for(uint_4 i=0; i<NElts(); i++) (*this)(i) *= x;
     274for(sa_size_t i=0; i<NElts(); i++) (*this)(i) *= x;
    333275return *this;
    334276}
    … …  
    337279template <class T> TMatrixRC<T>& TMatrixRC<T>::operator /= (T x)
    338280{
    339 for(uint_4 i=0; i<NElts(); i++) (*this)(i) /= x;
     281for(sa_size_t i=0; i<NElts(); i++) (*this)(i) /= x;
    340282return *this;
    341283}
    342284
    343285
    344 // ---- A virer   $CHECK$ Reza 03/2000
    345 // template <class T> TMatrixRC<T>& TMatrixRC<T>::operator -= (T x)
    346 // {
    347 // for(uint_4 i=0; i<NElts(); i++) (*this)(i) -= x;
    348 // return *this;
    349 // }
    350 
    351 // template <class T> TMatrixRC<T>& TMatrixRC<T>::operator += (T x)
    352 // {
    353 // for(uint_4 i=0; i<NElts(); i++) (*this)(i) += x;
    354 // return *this;
    355 // }
    356286
    357287//! Linear combination
    … …  
    361291 */
    362292template <class T>
    363 TMatrixRC<T>& TMatrixRC<T>::LinComb(T a, T b, const TMatrixRC<T>& rc, uint_4 first)
     293TMatrixRC<T>& TMatrixRC<T>::LinComb(T a, T b, const TMatrixRC<T>& rc, sa_size_t first)
    364294{
    365295if ( NElts() != rc.NElts() )
    … …  
    367297if ( kind != rc.kind )
    368298  throw(SzMismatchError("TMatrixRC::LinComb type mismatch\n"));
    369 for(uint_4 i=first; i<NElts(); i++) (*this)(i) = (*this)(i)*a + rc(i)*b;
     299for(sa_size_t i=first; i<NElts(); i++) (*this)(i) = (*this)(i)*a + rc(i)*b;
    370300return *this;
    371301}
    … …  
    376306 */
    377307template <class T>
    378 TMatrixRC<T>& TMatrixRC<T>::LinComb(T b, const TMatrixRC<T>& rc, uint_4 first)
     308TMatrixRC<T>& TMatrixRC<T>::LinComb(T b, const TMatrixRC<T>& rc, sa_size_t first)
    379309{
    380310if ( NElts() != rc.NElts() )
    … …  
    382312if ( kind != rc.kind )
    383313  throw(SzMismatchError("TMatrixRC::LinComb type mismatch\n"));
    384 for(uint_4 i=first; i<NElts(); i++) (*this)(i) += rc(i)*b;
     314for(sa_size_t i=first; i<NElts(); i++) (*this)(i) += rc(i)*b;
    385315return *this;
    386316}
    387317
    388318//! Find maximum absolute value in TMatrixRC, search begin at \b first
    389 template <class T> uint_4 TMatrixRC<T>::IMaxAbs(uint_4 first)
     319template <class T> sa_size_t TMatrixRC<T>::IMaxAbs(sa_size_t first)
    390320{
    391321if (first>NElts())
    392322  throw(SzMismatchError("TMatrixRC::IMaxAbs size mismatch\n"));
    393 uint_4 imax = first;
     323sa_size_t imax = first;
    394324double vmax = Abs_Value((*this)(first));
    395 for(uint_4 i=first+1; i<NElts(); i++) {
     325for(sa_size_t i=first+1; i<NElts(); i++) {
    396326  double v = Abs_Value((*this)(i));
    397327  if(v > vmax) {vmax = v; imax = i;}
    … …  
    406336  os << " TMatrixRC<T>::Print(ostream & os) " << NElts() << " Kind="
    407337     << kind << " Index=" << index << " Step= " << step << endl;
    408   for(uint_4 i=0; i<NElts(); i++) {
     338  for(sa_size_t i=0; i<NElts(); i++) {
    409339    os << (*this)(i);
    410340    if (kind == TmatrixCol)  os << endl;
    … …  
    423353  throw(SzMismatchError("TMatrixRC::Swap type mismatch\n"));
    424354if(rc1.data == rc2.data) return;
    425 for(uint_4 i=0; i<rc1.NElts(); i++)
     355for(sa_size_t i=0; i<rc1.NElts(); i++)
    426356  {T tmp = rc1(i); rc1(i) = rc2(i); rc2(i) = tmp;}
    427357}
    … …  
    468398// {TMatrix A(a); TMatrix B(b); return (T) TMatrix::GausPiv(A,B);}
    469399{
    470 uint_4 n = a.NRows();
     400sa_size_t n = a.NRows();
    471401if(n!=b.NRows())
    472402  throw(SzMismatchError("TMatrix::GausPiv size mismatch\n"));
    … …  
    484414} else if(GausPivScaling==PIV_ROW_SCALE) {
    485415  double nrm = 0.;
    486   for(uint_4 iii=0; iii<a.NRows(); iii++) {
    487     uint_4 jjj; double vmax = -1.;
     416  for(sa_size_t iii=0; iii<a.NRows(); iii++) {
     417    sa_size_t jjj; double vmax = -1.;
    488418    for(jjj=0; jjj<a.NCols(); jjj++) {
    489419      double v = TMatrixRC<T>::Abs_Value(a(iii,jjj));
    … …  
    505435} else {
    506436  double vmin=MAXDOUBLE, vmax=0;
    507   for(uint_4 iii=0; iii<a.NRows(); iii++)
    508     for(uint_4 jjj=0; jjj<a.NCols(); jjj++) {
     437  for(sa_size_t iii=0; iii<a.NRows(); iii++)
     438    for(sa_size_t jjj=0; jjj<a.NCols(); jjj++) {
    509439      double v = TMatrixRC<T>::Abs_Value(a(iii,jjj));
    510440      if(v>vmax) vmax = v;
    … …  
    529459TMatrixRC<T> pivRowb(b,TMatrixRC<T>::TmatrixRow);
    530460
    531 for(uint_4 k=0; k<n-1; k++) {
    532   uint_4 iPiv = TMatrixRC<T>::Col(a, k).IMaxAbs(k);
     461for(sa_size_t k=0; k<n-1; k++) {
     462  sa_size_t iPiv = TMatrixRC<T>::Col(a, k).IMaxAbs(k);
    533463  if(iPiv != k) {
    534464    TMatrixRC<T> aIPiv(TMatrixRC<T>::Row(a,iPiv));
    … …  
    544474  pivRowa.SetRow(k); // to avoid constructors
    545475  pivRowb.SetRow(k);
    546   for (uint_4 i=k+1; i<n; i++) {
     476  for (sa_size_t i=k+1; i<n; i++) {
    547477    T r = -a(i,k)/pivot;
    548478    TMatrixRC<T>::Row(a, i).LinComb(r, pivRowa); // + rapide que -= r * pivRowa
    … …  
    553483
    554484// on remonte
    555 for(uint_4 kk=n-1; kk>0; kk--) {
     485for(sa_size_t kk=n-1; kk>0; kk--) {
    556486  T pivot = a(kk,kk);
    557487  if( TMatrixRC<T>::Abs_Value(pivot) <= 1.e-50 ) return (T) 0;
    558488  pivRowa.SetRow(kk); // to avoid constructors
    559489  pivRowb.SetRow(kk);
    560   for(uint_4 jj=0; jj<kk; jj++) {
     490  for(sa_size_t jj=0; jj<kk; jj++) {
    561491    T r = -a(jj,kk)/pivot;
    562492    TMatrixRC<T>::Row(a, jj).LinComb(r, pivRowa);
    … …  
    565495}
    566496
    567 for(uint_4 l=0; l<n; l++) {
     497for(sa_size_t l=0; l<n; l++) {
    568498  if( TMatrixRC<T>::Abs_Value(a(l,l)) <= 1.e-50 ) return (T) 0;
    569499  TMatrixRC<T>::Row(b, l) /= a(l,l);
    … …  
    616546template <class T>
    617547r_8 LinFitter<T>::LinFit(const TVector<T>& x, const TVector<T>& y,
    618                        uint_4 nf, T (*f)(uint_4,T), TVector<T>& c)
    619 {
    620 uint_4 n = x.NElts();
     548                       sa_size_t nf, T (*f)(sa_size_t,T), TVector<T>& c)
     549{
     550sa_size_t n = x.NElts();
    621551if (n != y.NElts())
    622552  throw SzMismatchError("LinFitter<T>::LinFit(x,y,nf,f,c)/Error x.NElts() <> y.Nelts() ");
    … …  
    624554TMatrix<T> fx(nf,n);
    625555
    626 for(uint_4 i=0; i<nf; i++)
    627   for(uint_4 j=0; j<n; j++) fx(i,j) = f(i,x(j));
     556for(sa_size_t i=0; i<nf; i++)
     557  for(sa_size_t j=0; j<n; j++) fx(i,j) = f(i,x(j));
    628558
    629559return LinFit(fx,y,c);
    … …  
    643573r_8 LinFitter<T>::LinFit(const TMatrix<T>& fx, const TVector<T>& y, TVector<T>& c)
    644574{
    645 uint_4 n = y.NElts();
     575sa_size_t n = y.NElts();
    646576if (n != fx.NCol())
    647577  throw SzMismatchError("LinFitter<T>::LinFit(fx, y, c)/Error y.NElts() <> fx.Nelts() ");
    648578
    649 uint_4 nf = fx.NRows();
     579sa_size_t nf = fx.NRows();
    650580
    651581TMatrix<T> a(nf,nf);
    652582
    653 for(uint_4 j=0; j<nf; j++)
    654   for(uint_4 k=j; k<nf; k++)
     583for(sa_size_t j=0; j<nf; j++)
     584  for(sa_size_t k=j; k<nf; k++)
    655585     a(j,k) = a(k,j) = TMatrixRC<T>::Row(const_cast<TMatrix<T> &>(fx), j)
    656586                     * TMatrixRC<T>::Row(const_cast<TMatrix<T> &>(fx), k);
    … …  
    663593r_8 xi2 = 0., ax;
    664594T x;
    665 for(uint_4 k=0; k<n; k++) {
     595for(sa_size_t k=0; k<n; k++) {
    666596  x = (T) 0;
    667   for(uint_4 i=0; i<nf; i++) x += c(i) * fx(i,k);
     597  for(sa_size_t i=0; i<nf; i++) x += c(i) * fx(i,k);
    668598  x -= y(k);
    669599  ax = TMatrixRC<T>::Abs_Value(x);
    … …  
    690620template <class T>
    691621r_8 LinFitter<T>::LinFit(const TVector<T>& x, const TVector<T>& y,
    692              const TVector<T>& errY2, uint_4 nf, T (*f)(uint_4,T),
     622             const TVector<T>& errY2, sa_size_t nf, T (*f)(sa_size_t,T),
    693623                                  TVector<T>& c, TVector<T>& errC)
    694624{
    695 uint_4 n = x.NElts();
     625sa_size_t n = x.NElts();
    696626if (n != y.NElts())
    697627  throw SzMismatchError("LinFitter<T>::LinFit(x,y,errY2,nf,f,c,errC)/Error x.NElts() <> y.Nelts() ");
    698628
    699629TMatrix<T> fx(nf, n);
    700 for(uint_4 i=0; i<nf; i++)
    701   for(uint_4 j=0; j<n; j++)
     630for(sa_size_t i=0; i<nf; i++)
     631  for(sa_size_t j=0; j<n; j++)
    702632    fx(i,j) = f(i,x(j));
    703633
    … …  
    723653           const TVector<T>& errY2,TVector<T>& c, TVector<T>& errC)
    724654{
    725 uint_4 n = y.NElts();
     655sa_size_t n = y.NElts();
    726656if( n != errY2.NElts())
    727657  throw SzMismatchError("LinFitter<T>::LinFit(fx,y,errY2,c,errC)/Error y.NElts() <> errY2.Nelts() ");
    728658if( n != fx.NCol())
    729659  throw SzMismatchError("LinFitter<T>::LinFit(fx,y,errY2,c,errC)/Error y.NElts() <> fx.NCols() ");
    730 uint_4 nf = fx.NRows();
     660sa_size_t nf = fx.NRows();
    731661
    732662TMatrix<T> a(nf,nf);
    … …  
    735665errC.Realloc(nf);
    736666
    737 for(uint_4 j=0; j<nf; j++)
    738   for(uint_4 k=j; k<nf; k++) {
     667for(sa_size_t j=0; j<nf; j++)
     668  for(sa_size_t k=j; k<nf; k++) {
    739669    T x=0;
    740670    // Matrice a inverser
    741     for(uint_4 l=0; l<n; l++) x += fx(j,l)*fx(k,l)/errY2(l);
     671    for(sa_size_t l=0; l<n; l++) x += fx(j,l)*fx(k,l)/errY2(l);
    742672    a(j,k) = a(k,j) = x;
    743673  }
    744674 
    745675TMatrix<T> d(nf,nf+1);
    746 for(uint_4 k=0; k<nf; k++) {
     676for(sa_size_t k=0; k<nf; k++) {
    747677  T x = (T) 0;
    748678  // Second membre 1ere colonne
    749   for(uint_4 l=0; l<n; l++) x += fx(k,l)*y(l)/errY2(l);             
     679  for(sa_size_t l=0; l<n; l++) x += fx(k,l)*y(l)/errY2(l);             
    750680  d(k,0) = x;
    751681  // Reste second membre = Id.
    752   for(uint_4 m=1; m<=nf; m++) d(k,m) = (k==m)?1:0;
     682  for(sa_size_t m=1; m<=nf; m++) d(k,m) = (k==m)?1:0;
    753683}
    754684
    … …  
    757687
    758688
    759 for(uint_4 l=0; l<nf; l++) {
     689for(sa_size_t l=0; l<nf; l++) {
    760690  c(l) = d(l,0);        // Parametres = 1ere colonne
    761691  errC(l) = d(l,l+1);   // Erreurs = diag inverse.
    … …  
    764694r_8 xi2 = 0., ax;
    765695T x;
    766 for(uint_4 jj=0; jj<n; jj++) {
     696for(sa_size_t jj=0; jj<n; jj++) {
    767697  x = (T) 0;
    768   for(uint_4 ii=0; ii<nf; ii++) x += c(ii) * fx(ii,jj);
     698  for(sa_size_t ii=0; ii<nf; ii++) x += c(ii) * fx(ii,jj);
    769699    x -= y(jj);
    770700    ax = TMatrixRC<T>::Abs_Value(x);

  • trunk/SophyaLib/TArray/sopemtx.h

    r2917 r3332  
    167167//! Linear fitting
    168168r_8 LinFit(const TVector<T>& x, const TVector<T>& y,
    169            uint_4 nf, T (*f)(uint_4,T), TVector<T>& c);
     169           sa_size_t nf, T (*f)(sa_size_t,T), TVector<T>& c);
    170170
    171171//! Linear fitting
    … …  
    174174//! Linear fitting with errors
    175175r_8 LinFit(const TVector<T>& x, const TVector<T>& y, const TVector<T>& errY2,
    176            uint_4 nf,T (*f)(uint_4,T), TVector<T>& c, TVector<T>& errC);
     176           sa_size_t nf,T (*f)(sa_size_t,T), TVector<T>& c, TVector<T>& errC);
    177177
    178178//! Linear fitting with errors

  • trunk/SophyaLib/TArray/tarray.cc

    r3234 r3332  
    160160  string exmsg = "TArray<T>::TArray(int_4, sa_size_t *,  NDataBlock<T> & ... )";
    161161  if (!UpdateSizes(ndim, siz, step, offset, exmsg))  throw( ParmError(exmsg) );
    162   if (mNDBlock.Size() < ComputeTotalSize(ndim, siz, step, offset)) {
     162  if (mNDBlock.Size() < (size_t)ComputeTotalSize(ndim, siz, step, offset)) {
    163163    exmsg += " DataBlock.Size() < ComputeTotalSize(...) " ;
    164164    throw( ParmError(exmsg) );
    … …  
    13951395//! Returns the sum of all array elements squared (Sum_k((*this)(k)*(*this)(k)).
    13961396template <class T>
    1397 T TArray<T>::SumX2() const
    1398 {
    1399   if (NbDimensions() < 1)
    1400     throw RangeCheckError("TArray<T>::SumX2()  - Not Allocated Array ! ");
     1397T TArray<T>::SumSq() const
     1398{
     1399  if (NbDimensions() < 1)
     1400    throw RangeCheckError("TArray<T>::SumSq()  - Not Allocated Array ! ");
    14011401  T ret=0;
    14021402  const T * pe;
    … …  
    14201420  return ret;
    14211421}
     1422
     1423/*!
     1424\brief Returns the array norm squared, defined as  Sum_k [ el(k)* el(k) ]
     1425  For arrays with integer or real data, this method calls SumSq(), which computes
     1426  the sum of array elements squared. For complex arrays, it computes and returns
     1427  the sum of array elements module squared (= Sum_k [el(k)*conj(el(k))]
     1428*/
     1429template <class T>
     1430T TArray<T>::Norm2() const
     1431{
     1432  return SumSq();
     1433}
     1434
     1435
     1436// Fonction auxiliaire pour specialisation de la methode Norm2() pour tableaux complexes
     1437template <class T>
     1438complex<T>  _ComputeComplexNorm_Private_(TArray< complex<T> > const & ca)
     1439{
     1440  if (ca.NbDimensions() < 1)
     1441    throw RangeCheckError("TArray< complex<T> >::Norm2()  - Not Allocated Array ! ");
     1442  complex<T> ret= complex<T>(0., 0.);
     1443  const complex<T> * pe;
     1444  sa_size_t j,k;
     1445  if (ca.AvgStep() > 0)   {  // regularly spaced elements
     1446    sa_size_t step = ca.AvgStep();
     1447    sa_size_t maxx = ca.Size()*step;
     1448    pe = ca.Data();
     1449    for(k=0; k<maxx; k+=step )  ret += pe[k]*conj(pe[k]);
     1450  }
     1451  else {    // Non regular data spacing ...
     1452    int_4 ka = ca.MaxSizeKA();
     1453    sa_size_t step = ca.Step(ka);
     1454    sa_size_t gpas = ca.Size(ka)*step;
     1455    sa_size_t naxa = ca.Size()/ca.Size(ka);
     1456    for(j=0; j<naxa; j++)  {
     1457      pe = ca.DataBlock().Begin()+ca.Offset(ka,j);
     1458      for(k=0; k<gpas; k+=step)  ret += pe[k]*conj(pe[k]) ;
     1459    }
     1460  }
     1461  return ret;
     1462
     1463}
     1464
     1465// --- Specialisation de la methode Norm2() pour tableaux complexes ---
     1466DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
     1467complex<r_4> TArray< complex<r_4> >::Norm2() const
     1468{
     1469  return  _ComputeComplexNorm_Private_(*this);
     1470}
     1471DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
     1472complex<r_8> TArray< complex<r_8> >::Norm2() const
     1473{
     1474  return  _ComputeComplexNorm_Private_(*this);
     1475}
     1476//-------------------
    14221477
    14231478//! Return the minimum and the maximum values of the array elements

  • trunk/SophyaLib/TArray/tarray.h

    r3173 r3332  
    193193// Calcul du produit scalaire ( Somme_i (*this)(i)*a(i) )
    194194  virtual T ScalarProduct(const TArray<T>& a) const ;
    195 // Norme(^2)
    196   //! Returns the squarred of the array norm, defined as  Sum_k (*this)(k)*(*this)(k)
    197   inline T Norm2() const { return ScalarProduct(*this); }
    198195
    199196// Somme et produit des elements
    … …  
    201198  virtual T Product() const ;
    202199// Somme du carre des elements
    203   virtual T SumX2() const;
    204 // Valeur min et max des elements
     200  virtual T SumSq() const;
     201//  Norme^2 , identique a SumSq pour tableaux reels/integer , sum[el * conj(el)] pour complexe
     202  virtual T Norm2() const;
     203// Valeur min et max des elements (sauf tableaux complexes -> exception)
    205204  virtual void MinMax(T& min, T& max) const ;
    206205
Note: See TracChangeset for help on using the changeset viewer.