![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
Changeset 508 in Sophya for trunk/SophyaLib/NTools/matrix.cc
- Timestamp:
- Oct 25, 1999, 12:36:22 PM (26 years ago)
- File:
-
- 1 edited
-
trunk/SophyaLib/NTools/matrix.cc (modified) (67 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/SophyaLib/NTools/matrix.cc
r490 r508 1 // $Id: matrix.cc,v 1. 5 1999-10-21 15:25:49 ansari Exp $1 // $Id: matrix.cc,v 1.6 1999-10-25 10:36:09 ansari Exp $ 2 2 3 3 #include "machdefs.h" … … 26 26 27 27 //++ 28 // Class Matrix28 // Class OMatrix 29 29 // Lib Outils++ 30 30 // include matrix.h … … 38 38 39 39 //++ 40 Matrix::Matrix()40 OMatrix::OMatrix() 41 41 // 42 42 // Construit une matrice 1x1 (pour ppersist). … … 50 50 51 51 //++ 52 Matrix::Matrix(int r, int c)52 OMatrix::OMatrix(int r, int c) 53 53 // 54 54 // Construit une matrice de r lignes et c colonnes. … … 62 62 63 63 //++ 64 Matrix::Matrix(int r, int c, double* values)64 OMatrix::OMatrix(int r, int c, double* values) 65 65 // 66 66 // Construit une matrice de r lignes et c colonnes. On fournit … … 73 73 74 74 //++ 75 Matrix::Matrix(constMatrix& a)75 OMatrix::OMatrix(const OMatrix& a) 76 76 // 77 77 // Constructeur par copie. … … 85 85 86 86 //++ 87 Matrix::Matrix(const TMatrix<r_8>& a)87 OMatrix::OMatrix(const TMatrix<r_8>& a) 88 88 // 89 89 // Constructeur par copie a partir d'une TMatrix<r_8>. … … 99 99 100 100 101 Matrix::~Matrix()101 OMatrix::~OMatrix() 102 102 { 103 103 DBASSERT(ndata == nr*nc); … … 111 111 112 112 //++ 113 void Matrix::Zero()113 void OMatrix::Zero() 114 114 // 115 115 // Remise à zero de tous les éléments … … 121 121 122 122 //++ 123 void Matrix::Realloc(int r, int c, bool force)123 void OMatrix::Realloc(int r, int c, bool force) 124 124 // 125 125 // Change la taille de la matrice. Réallocation physique seulement si … … 145 145 146 146 //++ 147 Matrix& Matrix::operator = (constMatrix& a)147 OMatrix& OMatrix::operator = (const OMatrix& a) 148 148 // 149 149 // Opérateur d'affectation. … … 158 158 159 159 //++ 160 Matrix&Matrix::operator = (double x)160 OMatrix& OMatrix::operator = (double x) 161 161 // 162 162 // Opérateur d'affectation depuis scalaire : identité * scalaire. … … 176 176 177 177 //++ 178 ostream& operator << (ostream& s, const Matrix& a)178 ostream& operator << (ostream& s, const OMatrix& a) 179 179 // 180 180 // Impression … … 194 194 195 195 //++ 196 Matrix&Matrix::operator *= (double b)196 OMatrix& OMatrix::operator *= (double b) 197 197 // 198 198 //-- … … 208 208 209 209 //++ 210 Matrix&Matrix::operator += (double b)210 OMatrix& OMatrix::operator += (double b) 211 211 // 212 212 //-- … … 222 222 223 223 //++ 224 Matrix&Matrix::operator -= (double b)224 OMatrix& OMatrix::operator -= (double b) 225 225 // 226 226 //-- … … 236 236 237 237 //++ 238 Matrix&Matrix::operator /= (double b)238 OMatrix& OMatrix::operator /= (double b) 239 239 // 240 240 //-- … … 251 251 252 252 //++ 253 Matrix operator * (constMatrix& a, double b)254 // 255 //-- 256 { 257 Matrix result(a);253 OMatrix operator * (const OMatrix& a, double b) 254 // 255 //-- 256 { 257 OMatrix result(a); 258 258 return (result *= b); 259 259 } 260 260 261 261 //++ 262 Matrix operator * (double b, constMatrix& a)263 // 264 //-- 265 { 266 Matrix result(a);262 OMatrix operator * (double b, const OMatrix& a) 263 // 264 //-- 265 { 266 OMatrix result(a); 267 267 return (result *= b); 268 268 } 269 269 270 270 //++ 271 Matrix operator + (constMatrix& a, double b)272 // 273 //-- 274 { 275 Matrix result(a);271 OMatrix operator + (const OMatrix& a, double b) 272 // 273 //-- 274 { 275 OMatrix result(a); 276 276 return (result += b); 277 277 } 278 278 279 279 //++ 280 Matrix operator + (double b, constMatrix& a)281 // 282 //-- 283 { 284 Matrix result(a);280 OMatrix operator + (double b, const OMatrix& a) 281 // 282 //-- 283 { 284 OMatrix result(a); 285 285 return (result += b); 286 286 } 287 287 288 288 //++ 289 Matrix operator - (constMatrix& a, double b)290 // 291 //-- 292 { 293 Matrix result(a);289 OMatrix operator - (const OMatrix& a, double b) 290 // 291 //-- 292 { 293 OMatrix result(a); 294 294 return (result -= b); 295 295 } 296 296 297 297 //++ 298 Matrix operator - (double b, constMatrix& a)299 // 300 //-- 301 { 302 Matrix result(a);298 OMatrix operator - (double b, const OMatrix& a) 299 // 300 //-- 301 { 302 OMatrix result(a); 303 303 result *= -1; 304 304 return (result += b); … … 306 306 307 307 //++ 308 Matrix operator / (constMatrix& a, double b)309 // 310 //-- 311 { 312 Matrix result(a);308 OMatrix operator / (const OMatrix& a, double b) 309 // 310 //-- 311 { 312 OMatrix result(a); 313 313 return (result /= b); 314 314 } 315 315 316 316 //++ 317 Matrix operator * (constMatrix& a, int b)318 // 319 //-- 320 { 321 Matrix result(a);317 OMatrix operator * (const OMatrix& a, int b) 318 // 319 //-- 320 { 321 OMatrix result(a); 322 322 return (result *= b); 323 323 } 324 324 325 325 //++ 326 Matrix operator * (int b, constMatrix& a)327 // 328 //-- 329 { 330 Matrix result(a);326 OMatrix operator * (int b, const OMatrix& a) 327 // 328 //-- 329 { 330 OMatrix result(a); 331 331 return (result *= b); 332 332 } 333 333 334 334 //++ 335 Matrix operator + (constMatrix& a, int b)336 // 337 //-- 338 { 339 Matrix result(a);335 OMatrix operator + (const OMatrix& a, int b) 336 // 337 //-- 338 { 339 OMatrix result(a); 340 340 return (result += b); 341 341 } 342 342 343 343 //++ 344 Matrix operator + (int b, constMatrix& a)345 // 346 //-- 347 { 348 Matrix result(a);344 OMatrix operator + (int b, const OMatrix& a) 345 // 346 //-- 347 { 348 OMatrix result(a); 349 349 return (result += b); 350 350 } 351 351 352 352 //++ 353 Matrix operator - (constMatrix& a, int b)354 // 355 //-- 356 { 357 Matrix result(a);353 OMatrix operator - (const OMatrix& a, int b) 354 // 355 //-- 356 { 357 OMatrix result(a); 358 358 return (result -= b); 359 359 } 360 360 361 361 //++ 362 Matrix operator - (int b, constMatrix& a)363 // 364 //-- 365 { 366 Matrix result(a);362 OMatrix operator - (int b, const OMatrix& a) 363 // 364 //-- 365 { 366 OMatrix result(a); 367 367 result *= -1; 368 368 return (result += b); … … 370 370 371 371 //++ 372 Matrix operator / (constMatrix& a, int b)373 // 374 //-- 375 { 376 Matrix result(a);372 OMatrix operator / (const OMatrix& a, int b) 373 // 374 //-- 375 { 376 OMatrix result(a); 377 377 return (result /= b); 378 378 } … … 382 382 383 383 //++ 384 Matrix& Matrix::operator += (constMatrix& a)384 OMatrix& OMatrix::operator += (const OMatrix& a) 385 385 // 386 386 //-- … … 399 399 400 400 //++ 401 Matrix& Matrix::operator -= (constMatrix& a)401 OMatrix& OMatrix::operator -= (const OMatrix& a) 402 402 // 403 403 //-- … … 417 417 418 418 //++ 419 Matrix& Matrix::operator *= (constMatrix& a)419 OMatrix& OMatrix::operator *= (const OMatrix& a) 420 420 // 421 421 //-- … … 448 448 449 449 //++ 450 Matrix operator + (const Matrix& a, constMatrix& b)450 OMatrix operator + (const OMatrix& a, const OMatrix& b) 451 451 // 452 452 //-- … … 454 454 if (b.nc != a.nc || b.nr != a.nr) THROW(sizeMismatchErr); 455 455 456 Matrix c(a);456 OMatrix c(a); 457 457 return (c += b); 458 458 } 459 459 460 460 //++ 461 Matrix operator - (const Matrix& a, constMatrix& b)461 OMatrix operator - (const OMatrix& a, const OMatrix& b) 462 462 // 463 463 //-- … … 465 465 if (b.nc != a.nc || b.nr != a.nr) THROW(sizeMismatchErr); 466 466 467 Matrix c(a);467 OMatrix c(a); 468 468 return (c -= b); 469 469 } … … 471 471 #if 0 472 472 473 Matrix operator * (const Matrix& a, constMatrix& b)473 OMatrix operator * (const OMatrix& a, const OMatrix& b) 474 474 { 475 475 if (a.nc != b.nr) THROW(sizeMismatchErr); 476 476 477 Matrix c(a.nr, b.nc);477 OMatrix c(a.nr, b.nc); 478 478 479 479 double* pc = c.data; … … 505 505 506 506 //++ 507 Matrix operator * (const Matrix& a, constMatrix& b)507 OMatrix operator * (const OMatrix& a, const OMatrix& b) 508 508 // 509 509 //-- … … 511 511 if (a.nc != b.nr) THROW(sizeMismatchErr); 512 512 513 Matrix c(a.nr, b.nc);513 OMatrix c(a.nr, b.nc); 514 514 515 515 for (int rc = 0; rc < c.nr; rc++) … … 527 527 528 528 //++ 529 MatrixMatrix::Transpose() const529 OMatrix OMatrix::Transpose() const 530 530 // 531 531 // Retourne la transposée. … … 536 536 { 537 537 #if !HAS_NAMED_RETURN 538 Matrix a(nc,nr);538 OMatrix a(nc,nr); 539 539 #endif 540 540 for (int i=0; i<nr; i++) … … 548 548 549 549 //++ 550 MatrixMatrix::Inverse() const550 OMatrix OMatrix::Inverse() const 551 551 // 552 552 // Retourne la matrice inverse. … … 559 559 { 560 560 #if !HAS_NAMED_RETURN 561 Matrix b(nc,nr);562 #endif 563 Matrix a(*this);561 OMatrix b(nc,nr); 562 #endif 563 OMatrix a(*this); 564 564 b = 1.0; 565 if (fabs( Matrix::GausPiv(a,b)) < 1.e-50) THROW(singMatxErr);565 if (fabs(OMatrix::GausPiv(a,b)) < 1.e-50) THROW(singMatxErr); 566 566 #if !HAS_NAMED_RETURN 567 567 return b; … … 570 570 571 571 572 void Matrix::WriteSelf(POutPersist& s) const572 void OMatrix::WriteSelf(POutPersist& s) const 573 573 { 574 574 ASSERT(ndata == nr*nc); … … 577 577 } 578 578 579 void Matrix::ReadSelf(PInPersist& s)579 void OMatrix::ReadSelf(PInPersist& s) 580 580 { 581 581 int_4 r,c; … … 587 587 588 588 589 MatrixRC::MatrixRC()589 OMatrixRC::OMatrixRC() 590 590 : matrix(0), data(0), index(0), step(0) 591 591 {} 592 592 593 MatrixRC::MatrixRC(Matrix& m, RCKind rckind, int ind)593 OMatrixRC::OMatrixRC(OMatrix& m, RCKind rckind, int ind) 594 594 : matrix(&m), data(Org(m,rckind,ind)), 595 595 index(ind), step(Step(m,rckind)), kind(rckind) … … 599 599 600 600 601 int MatrixRC::Next()601 int OMatrixRC::Next() 602 602 { 603 603 if (!matrix) return -1; // Failure ? … … 621 621 622 622 623 int MatrixRC::Prev()623 int OMatrixRC::Prev() 624 624 { 625 625 if (!matrix) return -1; // Failure ? … … 638 638 639 639 640 int MatrixRC::SetCol(int c)640 int OMatrixRC::SetCol(int c) 641 641 { 642 642 if (!matrix) return -1; // Failure ? … … 650 650 651 651 652 int MatrixRC::SetRow(int r)652 int OMatrixRC::SetRow(int r) 653 653 { 654 654 if (!matrix) return -1; // Failure ? … … 662 662 663 663 664 int MatrixRC::SetDiag()664 int OMatrixRC::SetDiag() 665 665 { 666 666 if (!matrix) return -1; // Failure ? … … 674 674 675 675 676 MatrixRC& MatrixRC::operator = (constMatrixRC& rc)676 OMatrixRC& OMatrixRC::operator = (const OMatrixRC& rc) 677 677 { 678 678 matrix = rc.matrix; … … 684 684 } 685 685 686 // MatrixRC::operatorVector() const687 VectorMatrixRC::GetVect() const686 //OMatrixRC::operator OVector() const 687 OVector OMatrixRC::GetVect() const 688 688 #if HAS_NAMED_RETURN 689 689 return v(NElts()) … … 691 691 { 692 692 #if !HAS_NAMED_RETURN 693 Vector v(NElts());693 OVector v(NElts()); 694 694 #endif 695 695 int n = NElts(); … … 702 702 703 703 704 MatrixRC& MatrixRC::operator += (constMatrixRC& rc)704 OMatrixRC& OMatrixRC::operator += (const OMatrixRC& rc) 705 705 { 706 706 int n = NElts(); … … 716 716 717 717 718 MatrixRC& MatrixRC::operator -= (constMatrixRC& rc)718 OMatrixRC& OMatrixRC::operator -= (const OMatrixRC& rc) 719 719 { 720 720 int n = NElts(); … … 729 729 730 730 731 MatrixRC&MatrixRC::operator *= (double x)731 OMatrixRC& OMatrixRC::operator *= (double x) 732 732 { 733 733 int n = NElts(); … … 740 740 741 741 742 MatrixRC&MatrixRC::operator /= (double x)742 OMatrixRC& OMatrixRC::operator /= (double x) 743 743 { 744 744 int n = NElts(); … … 751 751 752 752 753 MatrixRC&MatrixRC::operator -= (double x)753 OMatrixRC& OMatrixRC::operator -= (double x) 754 754 { 755 755 int n = NElts(); … … 762 762 763 763 764 MatrixRC&MatrixRC::operator += (double x)764 OMatrixRC& OMatrixRC::operator += (double x) 765 765 { 766 766 int n = NElts(); … … 774 774 775 775 776 double operator * (const MatrixRC& a, constMatrixRC& b)776 double operator * (const OMatrixRC& a, const OMatrixRC& b) 777 777 { 778 778 int n = a.NElts(); … … 786 786 787 787 788 MatrixRC& MatrixRC::LinComb(double a, double b, constMatrixRC& rc, int first)788 OMatrixRC& OMatrixRC::LinComb(double a, double b, const OMatrixRC& rc, int first) 789 789 { 790 790 int n = NElts(); … … 799 799 800 800 801 MatrixRC& MatrixRC::LinComb(double b, constMatrixRC& rc, int first)801 OMatrixRC& OMatrixRC::LinComb(double b, const OMatrixRC& rc, int first) 802 802 { 803 803 int n = NElts(); … … 813 813 814 814 815 MatrixRCMatrix::Row(int r) const815 OMatrixRC OMatrix::Row(int r) const 816 816 #if HAS_NAMED_RETURN 817 return rc(( Matrix&)*this, matrixRow, r)817 return rc((OMatrix&)*this, matrixRow, r) 818 818 #endif 819 819 { 820 820 #if !HAS_NAMED_RETURN 821 MatrixRC rc((Matrix&)*this, matrixRow, r);821 OMatrixRC rc((OMatrix&)*this, matrixRow, r); 822 822 return rc; 823 823 #endif … … 825 825 826 826 827 MatrixRCMatrix::Col(int c) const827 OMatrixRC OMatrix::Col(int c) const 828 828 #if HAS_NAMED_RETURN 829 return rc(( Matrix&)*this, matrixCol, c)829 return rc((OMatrix&)*this, matrixCol, c) 830 830 #endif 831 831 { 832 832 #if !HAS_NAMED_RETURN 833 MatrixRC rc((Matrix&)*this, matrixCol, c);833 OMatrixRC rc((OMatrix&)*this, matrixCol, c); 834 834 return rc; 835 835 #endif … … 837 837 838 838 839 MatrixRCMatrix::Diag() const839 OMatrixRC OMatrix::Diag() const 840 840 #if HAS_NAMED_RETURN 841 return rc(( Matrix&)*this, matrixDiag)841 return rc((OMatrix&)*this, matrixDiag) 842 842 #endif 843 843 { 844 844 #if !HAS_NAMED_RETURN 845 MatrixRC rc((Matrix&)*this, matrixDiag);845 OMatrixRC rc((OMatrix&)*this, matrixDiag); 846 846 return rc; 847 847 #endif … … 849 849 850 850 851 int MatrixRC::IMaxAbs(int first)851 int OMatrixRC::IMaxAbs(int first) 852 852 { 853 853 int n=NElts(); … … 865 865 866 866 867 void MatrixRC::Swap(MatrixRC& rc1,MatrixRC& rc2)867 void OMatrixRC::Swap(OMatrixRC& rc1, OMatrixRC& rc2) 868 868 { 869 869 int n=rc1.NElts(); … … 877 877 878 878 879 double Matrix::GausPiv(Matrix& a,Matrix& b)879 double OMatrix::GausPiv(OMatrix& a, OMatrix& b) 880 880 { 881 881 int n = a.NRows(); … … 902 902 double ld = a.NRows() * log(nrm); 903 903 if (ld <= LN_MINDOUBLE || ld >= LN_MAXDOUBLE) { 904 // cerr << " Matrix warning, overflow for det" << endl;904 // cerr << "OMatrix warning, overflow for det" << endl; 905 905 } else { 906 906 det = exp(ld); … … 908 908 } 909 909 910 MatrixRC pivRowa(a,matrixRow);911 MatrixRC pivRowb(b,matrixRow);910 OMatrixRC pivRowa(a,matrixRow); 911 OMatrixRC pivRowb(b,matrixRow); 912 912 913 913 for (int k=0; k<n-1; k++) { 914 914 int iPiv = a.Col(k).IMaxAbs(k); 915 915 if (iPiv != k) { 916 MatrixRC aIPiv(a.Row(iPiv));917 MatrixRC aK(a.Row(k));918 MatrixRC::Swap(aIPiv,aK);919 MatrixRC bIPiv(b.Row(iPiv));920 MatrixRC bK(b.Row(k));921 MatrixRC::Swap(bIPiv,bK);916 OMatrixRC aIPiv(a.Row(iPiv)); 917 OMatrixRC aK(a.Row(k)); 918 OMatrixRC::Swap(aIPiv,aK); 919 OMatrixRC bIPiv(b.Row(iPiv)); 920 OMatrixRC bK(b.Row(k)); 921 OMatrixRC::Swap(bIPiv,bK); 922 922 } 923 923 double pivot = a(k,k); … … 958 958 959 959 //++ 960 double Matrix::Norm1()960 double OMatrix::Norm1() 961 961 // 962 962 // Norme 1 : somme des valeurs absolues. … … 970 970 971 971 //++ 972 double Matrix::Norm2()972 double OMatrix::Norm2() 973 973 // 974 974 // Norme 2, euclidienne. … … 983 983 ////////////////////////////////////////////////////////// 984 984 //++ 985 MatrixMatrix::FitResidus(GeneralFit& gfit985 OMatrix OMatrix::FitResidus(GeneralFit& gfit 986 986 ,double xorg,double yorg,double dx,double dy) 987 987 // … … 992 992 { 993 993 if(NCol()<=0||NRows()<=0) 994 throw(SzMismatchError(" Matrix::FitResidus: size mismatch\n"));994 throw(SzMismatchError("OMatrix::FitResidus: size mismatch\n")); 995 995 GeneralFunction* f = gfit.GetFunction(); 996 996 if(f==NULL) 997 throw(NullPtrError(" Matrix::FitResidus: NULL pointer\n"));998 Vector par = gfit.GetParm();999 Matrix m(*this);997 throw(NullPtrError("OMatrix::FitResidus: NULL pointer\n")); 998 OVector par = gfit.GetParm(); 999 OMatrix m(*this); 1000 1000 for(int i=0;i<NRows();i++) for(int j=0;j<NCol();j++) { 1001 1001 double x[2] = {xorg+j*dx,yorg+i*dy}; … … 1006 1006 1007 1007 //++ 1008 MatrixMatrix::FitFunction(GeneralFit& gfit1008 OMatrix OMatrix::FitFunction(GeneralFit& gfit 1009 1009 ,double xorg,double yorg,double dx,double dy) 1010 1010 // … … 1015 1015 { 1016 1016 if(NCol()<=0||NRows()<=0) 1017 throw(SzMismatchError(" Matrix::FitFunction: size mismatch\n"));1017 throw(SzMismatchError("OMatrix::FitFunction: size mismatch\n")); 1018 1018 GeneralFunction* f = gfit.GetFunction(); 1019 1019 if(f==NULL) 1020 throw(NullPtrError(" Matrix::FitFunction: NULL pointer\n"));1021 Vector par = gfit.GetParm();1022 Matrix m(*this);1020 throw(NullPtrError("OMatrix::FitFunction: NULL pointer\n")); 1021 OVector par = gfit.GetParm(); 1022 OMatrix m(*this); 1023 1023 for(int i=0;i<NRows();i++) for(int j=0;j<NCol();j++) { 1024 1024 double x[2] = {xorg+j*dx,yorg+i*dy};
Note:
See TracChangeset
for help on using the changeset viewer.
