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