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