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


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

File:
1 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);
Note: See TracChangeset for help on using the changeset viewer.