00001 // -*- C++ -*- 00002 // CLASSDOC OFF 00003 // $Id: Matrix.h,v 1.14 2002/04/12 15:02:44 evc 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 "CLHEP/Matrix/GenMatrix.h" 00232 #ifdef HEP_USE_RANDOM 00233 class HepRandom; 00234 #endif 00235 00236 #ifdef HEP_NO_INLINE_IN_DECLARATION 00237 #define inline 00238 #endif 00239 00240 class HepSymMatrix; 00241 class HepDiagMatrix; 00242 class HepVector; 00243 #ifdef HEP_USE_VECTOR_MODULE 00244 class HepRotation; 00245 #endif 00246 00251 class HepMatrix : public HepGenMatrix { 00252 public: 00253 inline HepMatrix(); 00254 // Default constructor. Gives 0 x 0 matrix. Another Matrix can be 00255 // assigned to it. 00256 00257 HepMatrix(int p, int q); 00258 // Constructor. Gives an unitialized p x q matrix. 00259 HepMatrix(int p, int q, int i); 00260 // Constructor. Gives an initialized p x q matrix. 00261 // If i=0, it is initialized to all 0. If i=1, the diagonal elements 00262 // are set to 1.0. 00263 00264 #ifdef HEP_USE_RANDOM 00265 HepMatrix(int p, int q, HepRandom &r); 00266 // Constructor with a Random object. 00267 #endif 00268 00269 HepMatrix(const HepMatrix &m1); 00270 // Copy constructor. 00271 00272 HepMatrix(const HepSymMatrix &m1); 00273 HepMatrix(const HepDiagMatrix &m1); 00274 HepMatrix(const HepVector &m1); 00275 // Constructors from SymMatrix, DiagMatrix and Vector. 00276 00277 virtual ~HepMatrix(); 00278 // Destructor. 00279 00280 inline virtual int num_row() const; 00281 // Returns the number of rows. 00282 00283 inline virtual int num_col() const; 00284 // Returns the number of columns. 00285 00286 inline virtual const double & operator()(int row, int col) const; 00287 inline virtual double & operator()(int row, int col); 00288 // Read or write a matrix element. 00289 // ** Note that the indexing starts from (1,1). ** 00290 00291 HepMatrix & operator *= (double t); 00292 // Multiply a Matrix by a floating number. 00293 00294 HepMatrix & operator /= (double t); 00295 // Divide a Matrix by a floating number. 00296 00297 HepMatrix & operator += ( const HepMatrix &m2); 00298 HepMatrix & operator += ( const HepSymMatrix &m2); 00299 HepMatrix & operator += ( const HepDiagMatrix &m2); 00300 HepMatrix & operator += ( const HepVector &m2); 00301 HepMatrix & operator -= ( const HepMatrix &m2); 00302 HepMatrix & operator -= ( const HepSymMatrix &m2); 00303 HepMatrix & operator -= ( const HepDiagMatrix &m2); 00304 HepMatrix & operator -= ( const HepVector &m2); 00305 // Add or subtract a Matrix. 00306 // When adding/subtracting Vector, Matrix must have num_col of one. 00307 00308 HepMatrix & operator = ( const HepMatrix &m2); 00309 HepMatrix & operator = ( const HepSymMatrix &m2); 00310 HepMatrix & operator = ( const HepDiagMatrix &m2); 00311 HepMatrix & operator = ( const HepVector &m2); 00312 #ifdef HEP_USE_VECTOR_MODULE 00313 HepMatrix & operator = ( const HepRotation &m2); 00314 #endif 00315 // Assignment operators. 00316 00317 HepMatrix operator- () const; 00318 // unary minus, ie. flip the sign of each element. 00319 00320 HepMatrix apply(double (*f)(double, int, int)) const; 00321 // Apply a function to all elements of the matrix. 00322 00323 HepMatrix T() const; 00324 // Returns the transpose of a Matrix. 00325 00326 HepMatrix sub(int min_row, int max_row, int min_col, int max_col) const; 00327 // Returns a sub matrix of a Matrix. 00328 // WARNING: rows and columns are numbered from 1 00329 void sub(int row, int col, const HepMatrix &m1); 00330 // Sub matrix of this Matrix is replaced with m1. 00331 // WARNING: rows and columns are numbered from 1 00332 00333 friend inline void swap(HepMatrix &m1, HepMatrix &m2); 00334 // Swap m1 with m2. 00335 00336 inline HepMatrix inverse(int& ierr) const; 00337 // Invert a Matrix. Matrix must be square and is not changed. 00338 // Returns ierr = 0 (zero) when successful, otherwise non-zero. 00339 00340 virtual void invert(int& ierr); 00341 // Invert a Matrix. Matrix must be square. 00342 // N.B. the contents of the matrix are replaced by the inverse. 00343 // Returns ierr = 0 (zero) when successful, otherwise non-zero. 00344 // This method has less overhead then inverse(). 00345 00346 double determinant() const; 00347 // calculate the determinant of the matrix. 00348 00349 double trace() const; 00350 // calculate the trace of the matrix (sum of diagonal elements). 00351 00352 class HepMatrix_row { 00353 public: 00354 inline HepMatrix_row(HepMatrix&,int); 00355 double & operator[](int); 00356 private: 00357 HepMatrix& _a; 00358 int _r; 00359 }; 00360 class HepMatrix_row_const { 00361 public: 00362 inline HepMatrix_row_const (const HepMatrix&,int); 00363 const double & operator[](int) const; 00364 private: 00365 const HepMatrix& _a; 00366 int _r; 00367 }; 00368 // helper classes for implementing m[i][j] 00369 00370 inline HepMatrix_row operator[] (int); 00371 inline const HepMatrix_row_const operator[] (int) const; 00372 // Read or write a matrix element. 00373 // While it may not look like it, you simply do m[i][j] to get an 00374 // element. 00375 // ** Note that the indexing starts from [0][0]. ** 00376 00377 protected: 00378 virtual inline int num_size() const; 00379 virtual void invertHaywood4(int& ierr); 00380 virtual void invertHaywood5(int& ierr); 00381 virtual void invertHaywood6(int& ierr); 00382 00383 private: 00384 friend class HepMatrix_row; 00385 friend class HepMatrix_row_const; 00386 friend class HepVector; 00387 friend class HepSymMatrix; 00388 friend class HepDiagMatrix; 00389 // Friend classes. 00390 00391 friend HepMatrix operator+(const HepMatrix &m1, const HepMatrix &m2); 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 HepSymMatrix &m2); 00395 friend HepMatrix operator*(const HepMatrix &m1, const HepDiagMatrix &m2); 00396 friend HepMatrix operator*(const HepSymMatrix &m1, const HepMatrix &m2); 00397 friend HepMatrix operator*(const HepDiagMatrix &m1, const HepMatrix &m2); 00398 friend HepMatrix operator*(const HepVector &m1, const HepMatrix &m2); 00399 friend HepVector operator*(const HepMatrix &m1, const HepVector &m2); 00400 friend HepMatrix operator*(const HepSymMatrix &m1, const HepSymMatrix &m2); 00401 // Multiply a Matrix by a Matrix or Vector. 00402 00403 friend HepVector solve(const HepMatrix &, const HepVector &); 00404 // solve the system of linear eq 00405 friend HepVector qr_solve(HepMatrix *, const HepVector &); 00406 friend HepMatrix qr_solve(HepMatrix *, const HepMatrix &b); 00407 friend void tridiagonal(HepSymMatrix *a,HepMatrix *hsm); 00408 friend void row_house(HepMatrix *,const HepMatrix &, double, 00409 int, int, int, int); 00410 friend void row_house(HepMatrix *,const HepVector &, double, 00411 int, int); 00412 friend void back_solve(const HepMatrix &R, HepVector *b); 00413 friend void back_solve(const HepMatrix &R, HepMatrix *b); 00414 friend void col_givens(HepMatrix *A, double c, 00415 double s, int k1, int k2, 00416 int rowmin, int rowmax); 00417 // Does a column Givens update. 00418 friend void row_givens(HepMatrix *A, double c, 00419 double s, int k1, int k2, 00420 int colmin, int colmax); 00421 friend void col_house(HepMatrix *,const HepMatrix &, double, 00422 int, int, int, int); 00423 friend HepVector house(const HepMatrix &a,int row,int col); 00424 friend void house_with_update(HepMatrix *a,int row,int col); 00425 friend void house_with_update(HepMatrix *a,HepMatrix *v,int row,int col); 00426 friend void house_with_update2(HepSymMatrix *a,HepMatrix *v, 00427 int row,int col); 00428 00429 int dfact_matrix(double &det, int *ir); 00430 // factorize the matrix. If successful, the return code is 0. On 00431 // return, det is the determinant and ir[] is row-interchange 00432 // matrix. See CERNLIB's DFACT routine. 00433 00434 int dfinv_matrix(int *ir); 00435 // invert the matrix. See CERNLIB DFINV. 00436 00437 double *m; 00438 int nrow, ncol; 00439 int size; 00440 }; 00441 00442 // Operations other than member functions for Matrix 00443 // implemented in Matrix.cc and Matrix.icc (inline). 00444 00445 HepMatrix operator*(const HepMatrix &m1, const HepMatrix &m2); 00446 HepMatrix operator*(double t, const HepMatrix &m1); 00447 HepMatrix operator*(const HepMatrix &m1, double t); 00448 // Multiplication operators 00449 // Note that m *= m1 is always faster than m = m * m1. 00450 00451 HepMatrix operator/(const HepMatrix &m1, double t); 00452 // m = m1 / t. (m /= t is faster if you can use it.) 00453 00454 HepMatrix operator+(const HepMatrix &m1, const HepMatrix &m2); 00455 // m = m1 + m2; 00456 // Note that m += m1 is always faster than m = m + m1. 00457 00458 HepMatrix operator-(const HepMatrix &m1, const HepMatrix &m2); 00459 // m = m1 - m2; 00460 // Note that m -= m1 is always faster than m = m - m1. 00461 00462 HepMatrix dsum(const HepMatrix&, const HepMatrix&); 00463 // Direct sum of two matrices. The direct sum of A and B is the matrix 00464 // A 0 00465 // 0 B 00466 00467 HepVector solve(const HepMatrix &, const HepVector &); 00468 // solve the system of linear equations using LU decomposition. 00469 00470 HepStd::ostream& operator<<(HepStd::ostream &s, const HepMatrix &q); 00471 // Read in, write out Matrix into a stream. 00472 00473 // 00474 // Specialized linear algebra functions 00475 // 00476 00477 HepVector qr_solve(const HepMatrix &A, const HepVector &b); 00478 HepVector qr_solve(HepMatrix *A, const HepVector &b); 00479 HepMatrix qr_solve(const HepMatrix &A, const HepMatrix &b); 00480 HepMatrix qr_solve(HepMatrix *A, const HepMatrix &b); 00481 // Works like backsolve, except matrix does not need to be upper 00482 // triangular. For nonsquare matrix, it solves in the least square sense. 00483 00484 HepMatrix qr_inverse(const HepMatrix &A); 00485 HepMatrix qr_inverse(HepMatrix *A); 00486 // Finds the inverse of a matrix using QR decomposition. Note, often what 00487 // you really want is solve or backsolve, they can be much quicker than 00488 // inverse in many calculations. 00489 00490 00491 void qr_decomp(HepMatrix *A, HepMatrix *hsm); 00492 HepMatrix qr_decomp(HepMatrix *A); 00493 // Does a QR decomposition of a matrix. 00494 00495 void back_solve(const HepMatrix &R, HepVector *b); 00496 void back_solve(const HepMatrix &R, HepMatrix *b); 00497 // Solves R*x = b where R is upper triangular. Also has a variation that 00498 // solves a number of equations of this form in one step, where b is a matrix 00499 // with each column a different vector. See also solve. 00500 00501 void col_house(HepMatrix *a, const HepMatrix &v, double vnormsq, 00502 int row, int col, int row_start, int col_start); 00503 void col_house(HepMatrix *a, const HepMatrix &v, int row, int col, 00504 int row_start, int col_start); 00505 // Does a column Householder update. 00506 00507 void col_givens(HepMatrix *A, double c, double s, 00508 int k1, int k2, int row_min=1, int row_max=0); 00509 // do a column Givens update 00510 00511 void row_givens(HepMatrix *A, double c, double s, 00512 int k1, int k2, int col_min=1, int col_max=0); 00513 // do a row Givens update 00514 00515 void givens(double a, double b, double *c, double *s); 00516 // algorithm 5.1.5 in Golub and Van Loan 00517 00518 HepVector house(const HepMatrix &a, int row=1, int col=1); 00519 // Returns a Householder vector to zero elements. 00520 00521 void house_with_update(HepMatrix *a, int row=1, int col=1); 00522 void house_with_update(HepMatrix *a, HepMatrix *v, int row=1, int col=1); 00523 // Finds and does Householder reflection on matrix. 00524 00525 void row_house(HepMatrix *a, const HepVector &v, double vnormsq, 00526 int row=1, int col=1); 00527 void row_house(HepMatrix *a, const HepMatrix &v, double vnormsq, 00528 int row, int col, int row_start, int col_start); 00529 void row_house(HepMatrix *a, const HepMatrix &v, int row, int col, 00530 int row_start, int col_start); 00531 // Does a row Householder update. 00532 00533 00534 #ifdef HEP_NO_INLINE_IN_DECLARATION 00535 #undef inline 00536 #endif 00537 00538 #ifdef HEP_SHORT_NAMES 00539 typedef HepMatrix Matrix; 00540 #endif 00541 00542 #ifndef HEP_DEBUG_INLINE 00543 #include "CLHEP/Matrix/Matrix.icc" 00544 #endif 00545 00546 #endif /*_Matrix_H*/