Changeset 3332 in Sophya for trunk/SophyaLib/TArray/sopemtx.cc
- Timestamp:
- Oct 3, 2007, 3:04:34 PM (18 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/SophyaLib/TArray/sopemtx.cc
r3012 r3332 36 36 }; 37 37 TMatrixRC(); 38 TMatrixRC(TMatrix<T>& m, TRCKind kind, uint_4index=0);38 TMatrixRC(TMatrix<T>& m, TRCKind kind, sa_size_t index=0); 39 39 virtual ~TMatrixRC() {} 40 40 41 41 // Acces aux rangees et colonnes de matrices 42 static TMatrixRC<T> Row(TMatrix<T> & m, uint_4r);43 static TMatrixRC<T> Col(TMatrix<T> & m, uint_4c);42 static TMatrixRC<T> Row(TMatrix<T> & m, sa_size_t r); 43 static TMatrixRC<T> Col(TMatrix<T> & m, sa_size_t c); 44 44 static TMatrixRC<T> Diag(TMatrix<T> & m); 45 45 46 // ---- A virer $CHECK$ Reza 03/200047 // int_4 Next();48 // int_4 Prev();49 46 int_4 SetCol(int_4 c); 50 47 int_4 SetRow(int_4 r); 51 48 int_4 SetDiag(); 52 49 53 static uint_4Step(const TMatrix<T>& m, TRCKind rckind);54 static T* Org(const TMatrix<T>&, TRCKind rckind, uint_4ind=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); 55 52 56 53 //! Return the kind of TMatrix (line,column,diagonal) 57 54 TRCKind Kind() const { return kind; } 58 uint_4NElts() const;59 T& operator()( uint_4i);60 T operator()( uint_4i) const;55 sa_size_t NElts() const; 56 T& operator()(sa_size_t i); 57 T operator()(sa_size_t i) const; 61 58 62 59 63 60 TMatrixRC<T>& operator = (const TMatrixRC<T>& rc); 64 61 65 // ---- A virer $CHECK$ Reza 03/200066 // TVector<T> GetVect() const;67 // TMatrixRC<T>& operator += (const TMatrixRC<T>& rc);68 // TMatrixRC<T>& operator -= (const TMatrixRC<T>& rc);69 62 70 63 TMatrixRC<T>& operator *= (T x); … … 75 68 76 69 77 TMatrixRC<T>& LinComb(T a, T b, const TMatrixRC& rc, uint_4first=0);78 TMatrixRC<T>& LinComb(T b, const TMatrixRC<T>& rc, uint_4first=0);79 80 uint_4 IMaxAbs(uint_4first=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); 81 74 void Print(ostream & os) const ; 82 75 … … 93 86 //! Define Absolute value for int_8 94 87 inline static double Abs_Value(int_8 v) {return (v>0)? (double) v: (double) -v;} 95 //! Define Absolute value for uint_488 //! Define Absolute value for sa_size_t 96 89 inline static double Abs_Value(uint_4 v) {return (double) v;} 97 90 //! Define Absolute value for uint_8 … … 112 105 T* data; //!< pointer to the beginnig of interesting datas 113 106 int_4 index; //!< index of the line/column 114 uint_4step; //!< step of the line/column107 sa_size_t step; //!< step of the line/column 115 108 TRCKind kind; //!< type: line, column or diagonal 116 109 }; … … 128 121 throw(SzMismatchError("TMatrixRC::operator * type mismatch\n")); 129 122 T sum = 0; 130 for( uint_4i=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); 131 124 return sum; 132 125 } … … 138 131 */ 139 132 template <class T> 140 inline uint_4TMatrixRC<T>::Step(const TMatrix<T>& m, TRCKind rckind)133 inline sa_size_t TMatrixRC<T>::Step(const TMatrix<T>& m, TRCKind rckind) 141 134 { switch (rckind) { case TmatrixRow : return m.Step(m.ColsKA()); 142 135 case TmatrixCol : return m.Step(m.RowsKA()); … … 153 146 */ 154 147 template <class T> 155 inline T* TMatrixRC<T>::Org(const TMatrix<T>& m, TRCKind rckind, uint_4index)148 inline T* TMatrixRC<T>::Org(const TMatrix<T>& m, TRCKind rckind, sa_size_t index) 156 149 { switch (rckind) { case TmatrixRow : return const_cast<T *>(&(m(index,0))); 157 150 case TmatrixCol : return const_cast<T *>(&(m(0,index))); … … 160 153 161 154 //! return number of elements for a TMatrixRC 162 template <class T> inline uint_4TMatrixRC<T>::NElts() const155 template <class T> inline sa_size_t TMatrixRC<T>::NElts() const 163 156 { if (!matrix) return 0; 164 157 switch (kind) { case TmatrixRow : return matrix->NCols(); … … 169 162 //! access of element \b i 170 163 template <class T> 171 inline T& TMatrixRC<T>::operator()( uint_4i) {return data[i*step];}164 inline T& TMatrixRC<T>::operator()(sa_size_t i) {return data[i*step];} 172 165 //! access of element \b i 173 166 template <class T> 174 inline T TMatrixRC<T>::operator()( uint_4i) const {return data[i*step];}167 inline T TMatrixRC<T>::operator()(sa_size_t i) const {return data[i*step];} 175 168 176 169 //////////////////////////////////////////////////////////////// … … 190 183 \param ind : number of the line or column 191 184 */ 192 template <class T> TMatrixRC<T>::TMatrixRC(TMatrix<T>& m,TRCKind rckind, uint_4ind)185 template <class T> TMatrixRC<T>::TMatrixRC(TMatrix<T>& m,TRCKind rckind,sa_size_t ind) 193 186 : matrix(&m), data(Org(m,rckind,ind)), 194 187 index(ind), step(Step(m,rckind)), kind(rckind) … … 203 196 //! Return TMatrixRC for line \b r of matrix \b m 204 197 template <class T> 205 TMatrixRC<T> TMatrixRC<T>::Row(TMatrix<T> & m, uint_4r)198 TMatrixRC<T> TMatrixRC<T>::Row(TMatrix<T> & m, sa_size_t r) 206 199 { 207 200 TMatrixRC<T> rc(m, TmatrixRow, r); … … 211 204 //! Return TMatrixRC for column \b r of matrix \b m 212 205 template <class T> 213 TMatrixRC<T> TMatrixRC<T>::Col(TMatrix<T> & m, uint_4c)206 TMatrixRC<T> TMatrixRC<T>::Col(TMatrix<T> & m, sa_size_t c) 214 207 { 215 208 TMatrixRC<T> rc(m, TmatrixCol, c); … … 226 219 227 220 228 // ---- A virer $CHECK$ Reza 03/2000229 // 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 // }252 221 253 222 //! Set column \b c for this TMatrixRC … … 299 268 } 300 269 301 // ---- A virer $CHECK$ Reza 03/2000302 // template <class T> TVector<T> TMatrixRC<T>::GetVect() const303 // {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 // }328 270 329 271 //! Operator to multiply by constant \b x 330 272 template <class T> TMatrixRC<T>& TMatrixRC<T>::operator *= (T x) 331 273 { 332 for( uint_4i=0; i<NElts(); i++) (*this)(i) *= x;274 for(sa_size_t i=0; i<NElts(); i++) (*this)(i) *= x; 333 275 return *this; 334 276 } … … 337 279 template <class T> TMatrixRC<T>& TMatrixRC<T>::operator /= (T x) 338 280 { 339 for( uint_4i=0; i<NElts(); i++) (*this)(i) /= x;281 for(sa_size_t i=0; i<NElts(); i++) (*this)(i) /= x; 340 282 return *this; 341 283 } 342 284 343 285 344 // ---- A virer $CHECK$ Reza 03/2000345 // 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 // }356 286 357 287 //! Linear combination … … 361 291 */ 362 292 template <class T> 363 TMatrixRC<T>& TMatrixRC<T>::LinComb(T a, T b, const TMatrixRC<T>& rc, uint_4first)293 TMatrixRC<T>& TMatrixRC<T>::LinComb(T a, T b, const TMatrixRC<T>& rc, sa_size_t first) 364 294 { 365 295 if ( NElts() != rc.NElts() ) … … 367 297 if ( kind != rc.kind ) 368 298 throw(SzMismatchError("TMatrixRC::LinComb type mismatch\n")); 369 for( uint_4i=first; i<NElts(); i++) (*this)(i) = (*this)(i)*a + rc(i)*b;299 for(sa_size_t i=first; i<NElts(); i++) (*this)(i) = (*this)(i)*a + rc(i)*b; 370 300 return *this; 371 301 } … … 376 306 */ 377 307 template <class T> 378 TMatrixRC<T>& TMatrixRC<T>::LinComb(T b, const TMatrixRC<T>& rc, uint_4first)308 TMatrixRC<T>& TMatrixRC<T>::LinComb(T b, const TMatrixRC<T>& rc, sa_size_t first) 379 309 { 380 310 if ( NElts() != rc.NElts() ) … … 382 312 if ( kind != rc.kind ) 383 313 throw(SzMismatchError("TMatrixRC::LinComb type mismatch\n")); 384 for( uint_4i=first; i<NElts(); i++) (*this)(i) += rc(i)*b;314 for(sa_size_t i=first; i<NElts(); i++) (*this)(i) += rc(i)*b; 385 315 return *this; 386 316 } 387 317 388 318 //! Find maximum absolute value in TMatrixRC, search begin at \b first 389 template <class T> uint_4 TMatrixRC<T>::IMaxAbs(uint_4first)319 template <class T> sa_size_t TMatrixRC<T>::IMaxAbs(sa_size_t first) 390 320 { 391 321 if (first>NElts()) 392 322 throw(SzMismatchError("TMatrixRC::IMaxAbs size mismatch\n")); 393 uint_4imax = first;323 sa_size_t imax = first; 394 324 double vmax = Abs_Value((*this)(first)); 395 for( uint_4i=first+1; i<NElts(); i++) {325 for(sa_size_t i=first+1; i<NElts(); i++) { 396 326 double v = Abs_Value((*this)(i)); 397 327 if(v > vmax) {vmax = v; imax = i;} … … 406 336 os << " TMatrixRC<T>::Print(ostream & os) " << NElts() << " Kind=" 407 337 << kind << " Index=" << index << " Step= " << step << endl; 408 for( uint_4i=0; i<NElts(); i++) {338 for(sa_size_t i=0; i<NElts(); i++) { 409 339 os << (*this)(i); 410 340 if (kind == TmatrixCol) os << endl; … … 423 353 throw(SzMismatchError("TMatrixRC::Swap type mismatch\n")); 424 354 if(rc1.data == rc2.data) return; 425 for( uint_4i=0; i<rc1.NElts(); i++)355 for(sa_size_t i=0; i<rc1.NElts(); i++) 426 356 {T tmp = rc1(i); rc1(i) = rc2(i); rc2(i) = tmp;} 427 357 } … … 468 398 // {TMatrix A(a); TMatrix B(b); return (T) TMatrix::GausPiv(A,B);} 469 399 { 470 uint_4n = a.NRows();400 sa_size_t n = a.NRows(); 471 401 if(n!=b.NRows()) 472 402 throw(SzMismatchError("TMatrix::GausPiv size mismatch\n")); … … 484 414 } else if(GausPivScaling==PIV_ROW_SCALE) { 485 415 double nrm = 0.; 486 for( uint_4iii=0; iii<a.NRows(); iii++) {487 uint_4jjj; double vmax = -1.;416 for(sa_size_t iii=0; iii<a.NRows(); iii++) { 417 sa_size_t jjj; double vmax = -1.; 488 418 for(jjj=0; jjj<a.NCols(); jjj++) { 489 419 double v = TMatrixRC<T>::Abs_Value(a(iii,jjj)); … … 505 435 } else { 506 436 double vmin=MAXDOUBLE, vmax=0; 507 for( uint_4iii=0; iii<a.NRows(); iii++)508 for( uint_4jjj=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++) { 509 439 double v = TMatrixRC<T>::Abs_Value(a(iii,jjj)); 510 440 if(v>vmax) vmax = v; … … 529 459 TMatrixRC<T> pivRowb(b,TMatrixRC<T>::TmatrixRow); 530 460 531 for( uint_4k=0; k<n-1; k++) {532 uint_4iPiv = TMatrixRC<T>::Col(a, k).IMaxAbs(k);461 for(sa_size_t k=0; k<n-1; k++) { 462 sa_size_t iPiv = TMatrixRC<T>::Col(a, k).IMaxAbs(k); 533 463 if(iPiv != k) { 534 464 TMatrixRC<T> aIPiv(TMatrixRC<T>::Row(a,iPiv)); … … 544 474 pivRowa.SetRow(k); // to avoid constructors 545 475 pivRowb.SetRow(k); 546 for ( uint_4i=k+1; i<n; i++) {476 for (sa_size_t i=k+1; i<n; i++) { 547 477 T r = -a(i,k)/pivot; 548 478 TMatrixRC<T>::Row(a, i).LinComb(r, pivRowa); // + rapide que -= r * pivRowa … … 553 483 554 484 // on remonte 555 for( uint_4kk=n-1; kk>0; kk--) {485 for(sa_size_t kk=n-1; kk>0; kk--) { 556 486 T pivot = a(kk,kk); 557 487 if( TMatrixRC<T>::Abs_Value(pivot) <= 1.e-50 ) return (T) 0; 558 488 pivRowa.SetRow(kk); // to avoid constructors 559 489 pivRowb.SetRow(kk); 560 for( uint_4jj=0; jj<kk; jj++) {490 for(sa_size_t jj=0; jj<kk; jj++) { 561 491 T r = -a(jj,kk)/pivot; 562 492 TMatrixRC<T>::Row(a, jj).LinComb(r, pivRowa); … … 565 495 } 566 496 567 for( uint_4l=0; l<n; l++) {497 for(sa_size_t l=0; l<n; l++) { 568 498 if( TMatrixRC<T>::Abs_Value(a(l,l)) <= 1.e-50 ) return (T) 0; 569 499 TMatrixRC<T>::Row(b, l) /= a(l,l); … … 616 546 template <class T> 617 547 r_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_4n = x.NElts();548 sa_size_t nf, T (*f)(sa_size_t,T), TVector<T>& c) 549 { 550 sa_size_t n = x.NElts(); 621 551 if (n != y.NElts()) 622 552 throw SzMismatchError("LinFitter<T>::LinFit(x,y,nf,f,c)/Error x.NElts() <> y.Nelts() "); … … 624 554 TMatrix<T> fx(nf,n); 625 555 626 for( uint_4i=0; i<nf; i++)627 for( uint_4j=0; j<n; j++) fx(i,j) = f(i,x(j));556 for(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)); 628 558 629 559 return LinFit(fx,y,c); … … 643 573 r_8 LinFitter<T>::LinFit(const TMatrix<T>& fx, const TVector<T>& y, TVector<T>& c) 644 574 { 645 uint_4n = y.NElts();575 sa_size_t n = y.NElts(); 646 576 if (n != fx.NCol()) 647 577 throw SzMismatchError("LinFitter<T>::LinFit(fx, y, c)/Error y.NElts() <> fx.Nelts() "); 648 578 649 uint_4nf = fx.NRows();579 sa_size_t nf = fx.NRows(); 650 580 651 581 TMatrix<T> a(nf,nf); 652 582 653 for( uint_4j=0; j<nf; j++)654 for( uint_4k=j; k<nf; k++)583 for(sa_size_t j=0; j<nf; j++) 584 for(sa_size_t k=j; k<nf; k++) 655 585 a(j,k) = a(k,j) = TMatrixRC<T>::Row(const_cast<TMatrix<T> &>(fx), j) 656 586 * TMatrixRC<T>::Row(const_cast<TMatrix<T> &>(fx), k); … … 663 593 r_8 xi2 = 0., ax; 664 594 T x; 665 for( uint_4k=0; k<n; k++) {595 for(sa_size_t k=0; k<n; k++) { 666 596 x = (T) 0; 667 for( uint_4i=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); 668 598 x -= y(k); 669 599 ax = TMatrixRC<T>::Abs_Value(x); … … 690 620 template <class T> 691 621 r_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), 693 623 TVector<T>& c, TVector<T>& errC) 694 624 { 695 uint_4n = x.NElts();625 sa_size_t n = x.NElts(); 696 626 if (n != y.NElts()) 697 627 throw SzMismatchError("LinFitter<T>::LinFit(x,y,errY2,nf,f,c,errC)/Error x.NElts() <> y.Nelts() "); 698 628 699 629 TMatrix<T> fx(nf, n); 700 for( uint_4i=0; i<nf; i++)701 for( uint_4j=0; j<n; j++)630 for(sa_size_t i=0; i<nf; i++) 631 for(sa_size_t j=0; j<n; j++) 702 632 fx(i,j) = f(i,x(j)); 703 633 … … 723 653 const TVector<T>& errY2,TVector<T>& c, TVector<T>& errC) 724 654 { 725 uint_4n = y.NElts();655 sa_size_t n = y.NElts(); 726 656 if( n != errY2.NElts()) 727 657 throw SzMismatchError("LinFitter<T>::LinFit(fx,y,errY2,c,errC)/Error y.NElts() <> errY2.Nelts() "); 728 658 if( n != fx.NCol()) 729 659 throw SzMismatchError("LinFitter<T>::LinFit(fx,y,errY2,c,errC)/Error y.NElts() <> fx.NCols() "); 730 uint_4nf = fx.NRows();660 sa_size_t nf = fx.NRows(); 731 661 732 662 TMatrix<T> a(nf,nf); … … 735 665 errC.Realloc(nf); 736 666 737 for( uint_4j=0; j<nf; j++)738 for( uint_4k=j; k<nf; k++) {667 for(sa_size_t j=0; j<nf; j++) 668 for(sa_size_t k=j; k<nf; k++) { 739 669 T x=0; 740 670 // Matrice a inverser 741 for( uint_4l=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); 742 672 a(j,k) = a(k,j) = x; 743 673 } 744 674 745 675 TMatrix<T> d(nf,nf+1); 746 for( uint_4k=0; k<nf; k++) {676 for(sa_size_t k=0; k<nf; k++) { 747 677 T x = (T) 0; 748 678 // Second membre 1ere colonne 749 for( uint_4l=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); 750 680 d(k,0) = x; 751 681 // Reste second membre = Id. 752 for( uint_4m=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; 753 683 } 754 684 … … 757 687 758 688 759 for( uint_4l=0; l<nf; l++) {689 for(sa_size_t l=0; l<nf; l++) { 760 690 c(l) = d(l,0); // Parametres = 1ere colonne 761 691 errC(l) = d(l,l+1); // Erreurs = diag inverse. … … 764 694 r_8 xi2 = 0., ax; 765 695 T x; 766 for( uint_4jj=0; jj<n; jj++) {696 for(sa_size_t jj=0; jj<n; jj++) { 767 697 x = (T) 0; 768 for( uint_4ii=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); 769 699 x -= y(jj); 770 700 ax = TMatrixRC<T>::Abs_Value(x);
Note:
See TracChangeset
for help on using the changeset viewer.