CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // -*- C++ -*- 00002 // CLASSDOC OFF 00003 // --------------------------------------------------------------------------- 00004 // CLASSDOC ON 00005 // 00006 // This file is a part of the CLHEP - a Class Library for High Energy Physics. 00007 // 00008 // This is the definition of the HepMatrix class. 00009 // 00010 // This software written by Nobu Katayama and Mike Smyth, Cornell University. 00011 // 00012 // 00013 // .SS Usage 00014 // The Matrix class does all the obvious things, in all the obvious ways. 00015 // You declare a Matrix by saying 00016 // 00017 // .ft B 00018 // HepMatrix m1(n, m); 00019 // 00020 // To declare a Matrix as a copy of a Matrix m2, say 00021 // 00022 // .ft B 00023 // HepMatrix m1(m2); 00024 // or 00025 // .ft B 00026 // HepMatrix m1 = m2; 00027 // 00028 // You can declare initilizations of a Matrix, by giving it a third 00029 // integer argument, either 0 or 1. 0 means initialize to 0, one means 00030 // the identity matrix. If you do not give a third argument the memory 00031 // is not initialized. 00032 // 00033 // .ft B 00034 // HepMatrix m1(n, m, 1); 00035 // 00036 // ./"This code has been written by Mike Smyth, and the algorithms used are 00037 // ./"described in the thesis "A Tracking Library for a Silicon Vertex Detector" 00038 // ./"(Mike Smyth, Cornell University, June 1993). 00039 // ./"This is file contains C++ stuff for doing things with Matrices. 00040 // ./"To turn on bound checking, define MATRIX_BOUND_CHECK before including 00041 // ./"this file. 00042 // 00043 // To find the number of rows or number of columns, say 00044 // 00045 // .ft B 00046 // nr = m1.num_row(); 00047 // 00048 // or 00049 // 00050 // .ft B 00051 // nc = m1.num_col(); 00052 // 00053 // You can print a Matrix by 00054 // 00055 // .ft B 00056 // cout << m1; 00057 // 00058 // You can add, 00059 // subtract, and multiply two Matrices. You can multiply a Matrix by a 00060 // scalar, on the left or the right. +=, *=, -= act in the obvious way. 00061 // m1 *= m2 is as m1 = m1*m2. You can also divide by a scalar on the 00062 // right, or use /= to do the same thing. 00063 // 00064 // You can read or write a Matrix element by saying 00065 // 00066 // .ft B 00067 // m(row, col) = blah. (1 <= row <= num_row(), 1 <= col <= num_col()) 00068 // 00069 // .ft B 00070 // blah = m(row, col) ... 00071 // 00072 // m(row, col) is inline, and by default does not do bound checking. 00073 // If bound checking is desired, say #define MATRIX_BOUND_CHECK before 00074 // including matrix.h. 00075 // 00076 // You can also access the element using C subscripting operator [] 00077 // 00078 // .ft B 00079 // m[row][col] = blah. (0 <= row < num_row(), 0 <= col < num_col()) 00080 // 00081 // .ft B 00082 // blah = m[row][col] ... 00083 // 00084 // m[row][col] is inline, and by default does not do bound checking. 00085 // If bound checking is desired, say #define MATRIX_BOUND_CHECK before 00086 // including matrix.h. (Notice the difference in bases in two access 00087 // methods.) 00088 // 00089 // .SS Comments on precision 00090 // 00091 // The user would normally use "Matrix" class whose precision is the same 00092 // as the other classes of CLHEP (ThreeVec, for example). However, he/she 00093 // can explicitly choose Matrix (double) or MatrixF (float) if he/she wishes. 00094 // In the former case, include "Matrix.h." In the latter case, include either 00095 // "Matrix.h," or "MatrixF.h," or both. The only operators that can talk 00096 // to each other are the copy constructor and assignment operator. 00097 // 00098 // .ft B 00099 // Matrix d(3,4,HEP_MATRIX_IDENTITY); 00100 // 00101 // .ft B 00102 // MatrixF f; 00103 // 00104 // .ft B 00105 // f = d; 00106 // 00107 // .ft B 00108 // MatrixF g(d); 00109 // 00110 // will convert from one to the other. 00111 // 00112 // .SS Other functions 00113 // 00114 // .ft B 00115 // mt = m.T(); 00116 // 00117 // returns the transpose of m. 00118 // 00119 // .ft B 00120 // ms = m2.sub(row_min, row_max, col_min, col_max); 00121 // 00122 // returns the subMatrix. 00123 // m2(row_min:row_max, col_min:col_max) in matlab terminology. 00124 // If instead you say 00125 // 00126 // .ft B 00127 // m2.sub(row, col, m1); 00128 // 00129 // then the subMatrix 00130 // m2(row:row+m1.num_row()-1, col:col+m1.num_col()-1) is replaced with m1. 00131 // 00132 // .ft B 00133 // md = dsum(m1,m2); 00134 // 00135 // returns the direct sum of the two matrices. 00136 // 00137 // Operations that do not have to be member functions or friends are listed 00138 // towards the end of this man page. 00139 // 00140 // 00141 // To invert a matrix, say 00142 // 00143 // .ft B 00144 // minv = m.inverse(ierr); 00145 // 00146 // ierr will be non-zero if the matrix is singular. 00147 // 00148 // If you can overwrite the matrix, you can use the invert function to avoid 00149 // two extra copies. Use 00150 // 00151 // .ft B 00152 // m.invert(ierr); 00153 // 00154 // Many basic linear algebra functions such as QR decomposition are defined. 00155 // In order to keep the header file a reasonable size, the functions are 00156 // defined in MatrixLinear.h. 00157 // 00158 // 00159 // .ft B 00160 // Note that Matrices run from (1,1) to (n, m), and [0][0] to 00161 // [n-1][m-1]. The users of the latter form should be careful about sub 00162 // functions. 00163 // 00164 // ./" The program that this class was orginally used with lots of small 00165 // ./" Matrices. It was very costly to malloc and free array space every 00166 // ./" time a Matrix is needed. So, a number of previously freed arrays are 00167 // ./" kept around, and when needed again one of these array is used. Right 00168 // ./" now, a set of piles of arrays with row <= row_max and col <= col_max 00169 // ./" are kept around. The piles of kept Matrices can be easily changed. 00170 // ./" Array mallocing and freeing are done only in new_m() and delete_m(), 00171 // ./" so these are the only routines that need to be rewritten. 00172 // 00173 // You can do things thinking of a Matrix as a list of numbers. 00174 // 00175 // .ft B 00176 // m = m1.apply(HEP_MATRIX_ELEMENT (*f)(HEP_MATRIX_ELEMENT, int r, int c)); 00177 // 00178 // applies f to every element of m1 and places it in m. 00179 // 00180 // .SS See Also: 00181 // SymMatrix[DF].h, GenMatrix[DF].h, DiagMatrix[DF].h Vector[DF].h 00182 // MatrixLinear[DF].h 00183 00184 #ifndef _Matrix_H_ 00185 #define _Matrix_H_ 00186 00187 #ifdef GNUPRAGMA 00188 #pragma interface 00189 #endif 00190 00191 #include <vector> 00192 00193 #include "CLHEP/Matrix/defs.h" 00194 #include "CLHEP/Matrix/GenMatrix.h" 00195 00196 namespace CLHEP { 00197 00198 class HepRandom; 00199 00200 class HepSymMatrix; 00201 class HepDiagMatrix; 00202 class HepVector; 00203 class HepRotation; 00204 00209 class HepMatrix : public HepGenMatrix { 00210 public: 00211 inline HepMatrix(); 00212 // Default constructor. Gives 0 x 0 matrix. Another Matrix can be 00213 // assigned to it. 00214 00215 HepMatrix(int p, int q); 00216 // Constructor. Gives an unitialized p x q matrix. 00217 HepMatrix(int p, int q, int i); 00218 // Constructor. Gives an initialized p x q matrix. 00219 // If i=0, it is initialized to all 0. If i=1, the diagonal elements 00220 // are set to 1.0. 00221 00222 HepMatrix(int p, int q, HepRandom &r); 00223 // Constructor with a Random object. 00224 00225 HepMatrix(const HepMatrix &m1); 00226 // Copy constructor. 00227 00228 HepMatrix(const HepSymMatrix &m1); 00229 HepMatrix(const HepDiagMatrix &m1); 00230 HepMatrix(const HepVector &m1); 00231 // Constructors from SymMatrix, DiagMatrix and Vector. 00232 00233 virtual ~HepMatrix(); 00234 // Destructor. 00235 00236 virtual int num_row() const; 00237 // Returns the number of rows. 00238 00239 virtual int num_col() const; 00240 // Returns the number of columns. 00241 00242 virtual const double & operator()(int row, int col) const; 00243 virtual double & operator()(int row, int col); 00244 // Read or write a matrix element. 00245 // ** Note that the indexing starts from (1,1). ** 00246 00247 HepMatrix & operator *= (double t); 00248 // Multiply a Matrix by a floating number. 00249 00250 HepMatrix & operator /= (double t); 00251 // Divide a Matrix by a floating number. 00252 00253 HepMatrix & operator += ( const HepMatrix &m2); 00254 HepMatrix & operator += ( const HepSymMatrix &m2); 00255 HepMatrix & operator += ( const HepDiagMatrix &m2); 00256 HepMatrix & operator += ( const HepVector &m2); 00257 HepMatrix & operator -= ( const HepMatrix &m2); 00258 HepMatrix & operator -= ( const HepSymMatrix &m2); 00259 HepMatrix & operator -= ( const HepDiagMatrix &m2); 00260 HepMatrix & operator -= ( const HepVector &m2); 00261 // Add or subtract a Matrix. 00262 // When adding/subtracting Vector, Matrix must have num_col of one. 00263 00264 HepMatrix & operator = ( const HepMatrix &m2); 00265 HepMatrix & operator = ( const HepSymMatrix &m2); 00266 HepMatrix & operator = ( const HepDiagMatrix &m2); 00267 HepMatrix & operator = ( const HepVector &m2); 00268 HepMatrix & operator = ( const HepRotation &m2); 00269 // Assignment operators. 00270 00271 HepMatrix operator- () const; 00272 // unary minus, ie. flip the sign of each element. 00273 00274 HepMatrix apply(double (*f)(double, int, int)) const; 00275 // Apply a function to all elements of the matrix. 00276 00277 HepMatrix T() const; 00278 // Returns the transpose of a Matrix. 00279 00280 HepMatrix sub(int min_row, int max_row, int min_col, int max_col) const; 00281 // Returns a sub matrix of a Matrix. 00282 // WARNING: rows and columns are numbered from 1 00283 void sub(int row, int col, const HepMatrix &m1); 00284 // Sub matrix of this Matrix is replaced with m1. 00285 // WARNING: rows and columns are numbered from 1 00286 00287 friend inline void swap(HepMatrix &m1, HepMatrix &m2); 00288 // Swap m1 with m2. 00289 00290 inline HepMatrix inverse(int& ierr) const; 00291 // Invert a Matrix. Matrix must be square and is not changed. 00292 // Returns ierr = 0 (zero) when successful, otherwise non-zero. 00293 00294 virtual void invert(int& ierr); 00295 // Invert a Matrix. Matrix must be square. 00296 // N.B. the contents of the matrix are replaced by the inverse. 00297 // Returns ierr = 0 (zero) when successful, otherwise non-zero. 00298 // This method has less overhead then inverse(). 00299 00300 inline void invert(); 00301 // Invert a matrix. Throw std::runtime_error on failure. 00302 00303 inline HepMatrix inverse() const; 00304 // Invert a matrix. Throw std::runtime_error on failure. 00305 00306 00307 double determinant() const; 00308 // calculate the determinant of the matrix. 00309 00310 double trace() const; 00311 // calculate the trace of the matrix (sum of diagonal elements). 00312 00313 class HepMatrix_row { 00314 public: 00315 inline HepMatrix_row(HepMatrix&,int); 00316 double & operator[](int); 00317 private: 00318 HepMatrix& _a; 00319 int _r; 00320 }; 00321 class HepMatrix_row_const { 00322 public: 00323 inline HepMatrix_row_const (const HepMatrix&,int); 00324 const double & operator[](int) const; 00325 private: 00326 const HepMatrix& _a; 00327 int _r; 00328 }; 00329 // helper classes for implementing m[i][j] 00330 00331 inline HepMatrix_row operator[] (int); 00332 inline const HepMatrix_row_const operator[] (int) const; 00333 // Read or write a matrix element. 00334 // While it may not look like it, you simply do m[i][j] to get an 00335 // element. 00336 // ** Note that the indexing starts from [0][0]. ** 00337 00338 protected: 00339 virtual int num_size() const; 00340 virtual void invertHaywood4(int& ierr); 00341 virtual void invertHaywood5(int& ierr); 00342 virtual void invertHaywood6(int& ierr); 00343 00344 private: 00345 friend class HepMatrix_row; 00346 friend class HepMatrix_row_const; 00347 friend class HepVector; 00348 friend class HepSymMatrix; 00349 friend class HepDiagMatrix; 00350 // Friend classes. 00351 00352 friend HepMatrix operator+(const HepMatrix &m1, const HepMatrix &m2); 00353 friend HepMatrix operator-(const HepMatrix &m1, const HepMatrix &m2); 00354 friend HepMatrix operator*(const HepMatrix &m1, const HepMatrix &m2); 00355 friend HepMatrix operator*(const HepMatrix &m1, const HepSymMatrix &m2); 00356 friend HepMatrix operator*(const HepMatrix &m1, const HepDiagMatrix &m2); 00357 friend HepMatrix operator*(const HepSymMatrix &m1, const HepMatrix &m2); 00358 friend HepMatrix operator*(const HepDiagMatrix &m1, const HepMatrix &m2); 00359 friend HepMatrix operator*(const HepVector &m1, const HepMatrix &m2); 00360 friend HepVector operator*(const HepMatrix &m1, const HepVector &m2); 00361 friend HepMatrix operator*(const HepSymMatrix &m1, const HepSymMatrix &m2); 00362 // Multiply a Matrix by a Matrix or Vector. 00363 00364 friend HepVector solve(const HepMatrix &, const HepVector &); 00365 // solve the system of linear eq 00366 friend HepVector qr_solve(HepMatrix *, const HepVector &); 00367 friend HepMatrix qr_solve(HepMatrix *, const HepMatrix &b); 00368 friend void tridiagonal(HepSymMatrix *a,HepMatrix *hsm); 00369 friend void row_house(HepMatrix *,const HepMatrix &, double, 00370 int, int, int, int); 00371 friend void row_house(HepMatrix *,const HepVector &, double, 00372 int, int); 00373 friend void back_solve(const HepMatrix &R, HepVector *b); 00374 friend void back_solve(const HepMatrix &R, HepMatrix *b); 00375 friend void col_givens(HepMatrix *A, double c, 00376 double s, int k1, int k2, 00377 int rowmin, int rowmax); 00378 // Does a column Givens update. 00379 friend void row_givens(HepMatrix *A, double c, 00380 double s, int k1, int k2, 00381 int colmin, int colmax); 00382 friend void col_house(HepMatrix *,const HepMatrix &, double, 00383 int, int, int, int); 00384 friend HepVector house(const HepMatrix &a,int row,int col); 00385 friend void house_with_update(HepMatrix *a,int row,int col); 00386 friend void house_with_update(HepMatrix *a,HepMatrix *v,int row,int col); 00387 friend void house_with_update2(HepSymMatrix *a,HepMatrix *v, 00388 int row,int col); 00389 00390 int dfact_matrix(double &det, int *ir); 00391 // factorize the matrix. If successful, the return code is 0. On 00392 // return, det is the determinant and ir[] is row-interchange 00393 // matrix. See CERNLIB's DFACT routine. 00394 00395 int dfinv_matrix(int *ir); 00396 // invert the matrix. See CERNLIB DFINV. 00397 00398 #ifdef DISABLE_ALLOC 00399 std::vector<double > m; 00400 #else 00401 std::vector<double,Alloc<double,25> > m; 00402 #endif 00403 int nrow, ncol; 00404 int size_; 00405 }; 00406 00407 // Operations other than member functions for Matrix 00408 // implemented in Matrix.cc and Matrix.icc (inline). 00409 00410 HepMatrix operator*(const HepMatrix &m1, const HepMatrix &m2); 00411 HepMatrix operator*(double t, const HepMatrix &m1); 00412 HepMatrix operator*(const HepMatrix &m1, double t); 00413 // Multiplication operators 00414 // Note that m *= m1 is always faster than m = m * m1. 00415 00416 HepMatrix operator/(const HepMatrix &m1, double t); 00417 // m = m1 / t. (m /= t is faster if you can use it.) 00418 00419 HepMatrix operator+(const HepMatrix &m1, const HepMatrix &m2); 00420 // m = m1 + m2; 00421 // Note that m += m1 is always faster than m = m + m1. 00422 00423 HepMatrix operator-(const HepMatrix &m1, const HepMatrix &m2); 00424 // m = m1 - m2; 00425 // Note that m -= m1 is always faster than m = m - m1. 00426 00427 HepMatrix dsum(const HepMatrix&, const HepMatrix&); 00428 // Direct sum of two matrices. The direct sum of A and B is the matrix 00429 // A 0 00430 // 0 B 00431 00432 HepVector solve(const HepMatrix &, const HepVector &); 00433 // solve the system of linear equations using LU decomposition. 00434 00435 std::ostream& operator<<(std::ostream &s, const HepMatrix &q); 00436 // Read in, write out Matrix into a stream. 00437 00438 // 00439 // Specialized linear algebra functions 00440 // 00441 00442 HepVector qr_solve(const HepMatrix &A, const HepVector &b); 00443 HepVector qr_solve(HepMatrix *A, const HepVector &b); 00444 HepMatrix qr_solve(const HepMatrix &A, const HepMatrix &b); 00445 HepMatrix qr_solve(HepMatrix *A, const HepMatrix &b); 00446 // Works like backsolve, except matrix does not need to be upper 00447 // triangular. For nonsquare matrix, it solves in the least square sense. 00448 00449 HepMatrix qr_inverse(const HepMatrix &A); 00450 HepMatrix qr_inverse(HepMatrix *A); 00451 // Finds the inverse of a matrix using QR decomposition. Note, often what 00452 // you really want is solve or backsolve, they can be much quicker than 00453 // inverse in many calculations. 00454 00455 00456 void qr_decomp(HepMatrix *A, HepMatrix *hsm); 00457 HepMatrix qr_decomp(HepMatrix *A); 00458 // Does a QR decomposition of a matrix. 00459 00460 void back_solve(const HepMatrix &R, HepVector *b); 00461 void back_solve(const HepMatrix &R, HepMatrix *b); 00462 // Solves R*x = b where R is upper triangular. Also has a variation that 00463 // solves a number of equations of this form in one step, where b is a matrix 00464 // with each column a different vector. See also solve. 00465 00466 void col_house(HepMatrix *a, const HepMatrix &v, double vnormsq, 00467 int row, int col, int row_start, int col_start); 00468 void col_house(HepMatrix *a, const HepMatrix &v, int row, int col, 00469 int row_start, int col_start); 00470 // Does a column Householder update. 00471 00472 void col_givens(HepMatrix *A, double c, double s, 00473 int k1, int k2, int row_min=1, int row_max=0); 00474 // do a column Givens update 00475 00476 void row_givens(HepMatrix *A, double c, double s, 00477 int k1, int k2, int col_min=1, int col_max=0); 00478 // do a row Givens update 00479 00480 void givens(double a, double b, double *c, double *s); 00481 // algorithm 5.1.5 in Golub and Van Loan 00482 00483 HepVector house(const HepMatrix &a, int row=1, int col=1); 00484 // Returns a Householder vector to zero elements. 00485 00486 void house_with_update(HepMatrix *a, int row=1, int col=1); 00487 void house_with_update(HepMatrix *a, HepMatrix *v, int row=1, int col=1); 00488 // Finds and does Householder reflection on matrix. 00489 00490 void row_house(HepMatrix *a, const HepVector &v, double vnormsq, 00491 int row=1, int col=1); 00492 void row_house(HepMatrix *a, const HepMatrix &v, double vnormsq, 00493 int row, int col, int row_start, int col_start); 00494 void row_house(HepMatrix *a, const HepMatrix &v, int row, int col, 00495 int row_start, int col_start); 00496 // Does a row Householder update. 00497 00498 } // namespace CLHEP 00499 00500 #ifdef ENABLE_BACKWARDS_COMPATIBILITY 00501 // backwards compatibility will be enabled ONLY in CLHEP 1.9 00502 using namespace CLHEP; 00503 #endif 00504 00505 #ifndef HEP_DEBUG_INLINE 00506 #include "CLHEP/Matrix/Matrix.icc" 00507 #endif 00508 00509 #endif /*_Matrix_H*/