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 HepSymMatrix 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 // .SS Usage 00041 // 00042 // This is very much like the Matrix, except of course it is restricted to 00043 // Symmetric Matrix. All the operations for Matrix can also be done here 00044 // (except for the +=,-=,*= that don't yield a symmetric matrix. e.g. 00045 // +=(const Matrix &) is not defined) 00046 00047 // The matrix is stored as a lower triangular matrix. 00048 // In addition to the (row, col) method of finding element, fast(row, col) 00049 // returns an element with checking to see if row and col need to be 00050 // interchanged so that row >= col. 00051 00052 // New operations are: 00053 // 00054 // .ft B 00055 // sym = s.similarity(m); 00056 // 00057 // This returns m*s*m.T(). This is a similarity 00058 // transform when m is orthogonal, but nothing 00059 // restricts m to be orthogonal. It is just 00060 // a more efficient way to calculate m*s*m.T, 00061 // and it realizes that this should be a 00062 // HepSymMatrix (the explicit operation m*s*m.T 00063 // will return a Matrix, not realizing that 00064 // it is symmetric). 00065 // 00066 // .ft B 00067 // sym = similarityT(m); 00068 // 00069 // This returns m.T()*s*m. 00070 // 00071 // .ft B 00072 // s << m; 00073 // 00074 // This takes the matrix m, and treats it 00075 // as symmetric matrix that is copied to s. 00076 // This is useful for operations that yield 00077 // symmetric matrix, but which the computer 00078 // is too dumb to realize. 00079 // 00080 // .ft B 00081 // s = vT_times_v(const HepVector v); 00082 // 00083 // calculates v.T()*v. 00084 // 00085 // ./"This code has been written by Mike Smyth, and the algorithms used are 00086 // ./"described in the thesis "A Tracking Library for a Silicon Vertex Detector" 00087 // ./"(Mike Smyth, Cornell University, June 1993). 00088 // ./"Copyright (C) Cornell University 1993. Permission is granted to copy and 00089 // ./"distribute this code, provided this copyright is not changed or deleted. 00090 // ./"You may modify your copy, providing that you cause the modified file to 00091 // ./"carry prominent notices stating that you changed the files, and the date 00092 // ./"of any change. This code may not be sold, nor may it be contained in 00093 // ./"programs that are to be sold without the written permission of the author. 00094 // ./"You may, however, charge a fee for the physical act of transferring a 00095 // ./"copy of this code. The code is offered "as is" without warranty of any 00096 // ./"kind, either expressed or implied. By copying, distributing, or 00097 // ./"modifying this code you indicate your acceptance of this license to 00098 // ./"do so, and all its terms and conditions. 00099 // ./"This is file contains C++ stuff for doing things with Matrixes. 00100 // ./"To turn on bound checking, define MATRIX_BOUND_CHECK before including 00101 // ./"this file. 00102 // 00103 00104 #ifndef _SYMMatrix_H_ 00105 #define _SYMMatrix_H_ 00106 00107 #ifdef GNUPRAGMA 00108 #pragma interface 00109 #endif 00110 00111 #include <vector> 00112 00113 #include "CLHEP/Matrix/defs.h" 00114 #include "CLHEP/Matrix/GenMatrix.h" 00115 00116 namespace CLHEP { 00117 00118 class HepRandom; 00119 00120 class HepMatrix; 00121 class HepDiagMatrix; 00122 class HepVector; 00123 00128 class HepSymMatrix : public HepGenMatrix { 00129 public: 00130 inline HepSymMatrix(); 00131 // Default constructor. Gives 0x0 symmetric matrix. 00132 // Another SymMatrix can be assigned to it. 00133 00134 explicit HepSymMatrix(int p); 00135 HepSymMatrix(int p, int); 00136 // Constructor. Gives p x p symmetric matrix. 00137 // With a second argument, the matrix is initialized. 0 means a zero 00138 // matrix, 1 means the identity matrix. 00139 00140 HepSymMatrix(int p, HepRandom &r); 00141 00142 HepSymMatrix(const HepSymMatrix &m1); 00143 // Copy constructor. 00144 00145 HepSymMatrix(const HepDiagMatrix &m1); 00146 // Constructor from DiagMatrix 00147 00148 virtual ~HepSymMatrix(); 00149 // Destructor. 00150 00151 inline int num_row() const; 00152 inline int num_col() const; 00153 // Returns number of rows/columns. 00154 00155 const double & operator()(int row, int col) const; 00156 double & operator()(int row, int col); 00157 // Read and write a SymMatrix element. 00158 // ** Note that indexing starts from (1,1). ** 00159 00160 const double & fast(int row, int col) const; 00161 double & fast(int row, int col); 00162 // fast element access. 00163 // Must be row>=col; 00164 // ** Note that indexing starts from (1,1). ** 00165 00166 void assign(const HepMatrix &m2); 00167 // Assigns m2 to s, assuming m2 is a symmetric matrix. 00168 00169 void assign(const HepSymMatrix &m2); 00170 // Another form of assignment. For consistency. 00171 00172 HepSymMatrix & operator*=(double t); 00173 // Multiply a SymMatrix by a floating number. 00174 00175 HepSymMatrix & operator/=(double t); 00176 // Divide a SymMatrix by a floating number. 00177 00178 HepSymMatrix & operator+=( const HepSymMatrix &m2); 00179 HepSymMatrix & operator+=( const HepDiagMatrix &m2); 00180 HepSymMatrix & operator-=( const HepSymMatrix &m2); 00181 HepSymMatrix & operator-=( const HepDiagMatrix &m2); 00182 // Add or subtract a SymMatrix. 00183 00184 HepSymMatrix & operator=( const HepSymMatrix &m2); 00185 HepSymMatrix & operator=( const HepDiagMatrix &m2); 00186 // Assignment operators. Notice that there is no SymMatrix = Matrix. 00187 00188 HepSymMatrix operator- () const; 00189 // unary minus, ie. flip the sign of each element. 00190 00191 HepSymMatrix T() const; 00192 // Returns the transpose of a SymMatrix (which is itself). 00193 00194 HepSymMatrix apply(double (*f)(double, int, int)) const; 00195 // Apply a function to all elements of the matrix. 00196 00197 HepSymMatrix similarity(const HepMatrix &m1) const; 00198 HepSymMatrix similarity(const HepSymMatrix &m1) const; 00199 // Returns m1*s*m1.T(). 00200 00201 HepSymMatrix similarityT(const HepMatrix &m1) const; 00202 // temporary. test of new similarity. 00203 // Returns m1.T()*s*m1. 00204 00205 double similarity(const HepVector &v) const; 00206 // Returns v.T()*s*v (This is a scaler). 00207 00208 HepSymMatrix sub(int min_row, int max_row) const; 00209 // Returns a sub matrix of a SymMatrix. 00210 void sub(int row, const HepSymMatrix &m1); 00211 // Sub matrix of this SymMatrix is replaced with m1. 00212 HepSymMatrix sub(int min_row, int max_row); 00213 // SGI CC bug. I have to have both with/without const. I should not need 00214 // one without const. 00215 00216 inline HepSymMatrix inverse(int &ifail) const; 00217 // Invert a Matrix. The matrix is not changed 00218 // Returns 0 when successful, otherwise non-zero. 00219 00220 void invert(int &ifail); 00221 // Invert a Matrix. 00222 // N.B. the contents of the matrix are replaced by the inverse. 00223 // Returns ierr = 0 when successful, otherwise non-zero. 00224 // This method has less overhead then inverse(). 00225 00226 double determinant() const; 00227 // calculate the determinant of the matrix. 00228 00229 double trace() const; 00230 // calculate the trace of the matrix (sum of diagonal elements). 00231 00232 class HepSymMatrix_row { 00233 public: 00234 inline HepSymMatrix_row(HepSymMatrix&,int); 00235 inline double & operator[](int); 00236 private: 00237 HepSymMatrix& _a; 00238 int _r; 00239 }; 00240 class HepSymMatrix_row_const { 00241 public: 00242 inline HepSymMatrix_row_const(const HepSymMatrix&,int); 00243 inline const double & operator[](int) const; 00244 private: 00245 const HepSymMatrix& _a; 00246 int _r; 00247 }; 00248 // helper class to implement m[i][j] 00249 00250 inline HepSymMatrix_row operator[] (int); 00251 inline HepSymMatrix_row_const operator[] (int) const; 00252 // Read or write a matrix element. 00253 // While it may not look like it, you simply do m[i][j] to get an 00254 // element. 00255 // ** Note that the indexing starts from [0][0]. ** 00256 00257 // Special-case inversions for 5x5 and 6x6 symmetric positive definite: 00258 // These set ifail=0 and invert if the matrix was positive definite; 00259 // otherwise ifail=1 and the matrix is left unaltered. 00260 void invertCholesky5 (int &ifail); 00261 void invertCholesky6 (int &ifail); 00262 00263 // Inversions for 5x5 and 6x6 forcing use of specific methods: The 00264 // behavior (though not the speed) will be identical to invert(ifail). 00265 void invertHaywood4 (int & ifail); 00266 void invertHaywood5 (int &ifail); 00267 void invertHaywood6 (int &ifail); 00268 void invertBunchKaufman (int &ifail); 00269 00270 protected: 00271 inline int num_size() const; 00272 00273 private: 00274 friend class HepSymMatrix_row; 00275 friend class HepSymMatrix_row_const; 00276 friend class HepMatrix; 00277 friend class HepDiagMatrix; 00278 00279 friend void tridiagonal(HepSymMatrix *a,HepMatrix *hsm); 00280 friend double condition(const HepSymMatrix &m); 00281 friend void diag_step(HepSymMatrix *t,int begin,int end); 00282 friend void diag_step(HepSymMatrix *t,HepMatrix *u,int begin,int end); 00283 friend HepMatrix diagonalize(HepSymMatrix *s); 00284 friend HepVector house(const HepSymMatrix &a,int row,int col); 00285 friend void house_with_update2(HepSymMatrix *a,HepMatrix *v,int row,int col); 00286 00287 friend HepSymMatrix operator+(const HepSymMatrix &m1, 00288 const HepSymMatrix &m2); 00289 friend HepSymMatrix operator-(const HepSymMatrix &m1, 00290 const HepSymMatrix &m2); 00291 friend HepMatrix operator*(const HepSymMatrix &m1, const HepSymMatrix &m2); 00292 friend HepMatrix operator*(const HepSymMatrix &m1, const HepMatrix &m2); 00293 friend HepMatrix operator*(const HepMatrix &m1, const HepSymMatrix &m2); 00294 friend HepVector operator*(const HepSymMatrix &m1, const HepVector &m2); 00295 // Multiply a Matrix by a Matrix or Vector. 00296 00297 friend HepSymMatrix vT_times_v(const HepVector &v); 00298 // Returns v * v.T(); 00299 00300 #ifdef DISABLE_ALLOC 00301 std::vector<double > m; 00302 #else 00303 std::vector<double,Alloc<double,25> > m; 00304 #endif 00305 int nrow; 00306 int size; // total number of elements 00307 00308 static double posDefFraction5x5; 00309 static double adjustment5x5; 00310 static const double CHOLESKY_THRESHOLD_5x5; 00311 static const double CHOLESKY_CREEP_5x5; 00312 00313 static double posDefFraction6x6; 00314 static double adjustment6x6; 00315 static const double CHOLESKY_THRESHOLD_6x6; 00316 static const double CHOLESKY_CREEP_6x6; 00317 00318 void invert4 (int & ifail); 00319 void invert5 (int & ifail); 00320 void invert6 (int & ifail); 00321 00322 }; 00323 00324 // 00325 // Operations other than member functions for Matrix, SymMatrix, DiagMatrix 00326 // and Vectors implemented in Matrix.cc and Matrix.icc (inline). 00327 // 00328 00329 std::ostream& operator<<(std::ostream &s, const HepSymMatrix &q); 00330 // Write out Matrix, SymMatrix, DiagMatrix and Vector into ostream. 00331 00332 HepMatrix operator*(const HepMatrix &m1, const HepSymMatrix &m2); 00333 HepMatrix operator*(const HepSymMatrix &m1, const HepMatrix &m2); 00334 HepMatrix operator*(const HepSymMatrix &m1, const HepSymMatrix &m2); 00335 HepSymMatrix operator*(double t, const HepSymMatrix &s1); 00336 HepSymMatrix operator*(const HepSymMatrix &s1, double t); 00337 // Multiplication operators. 00338 // Note that m *= m1 is always faster than m = m * m1 00339 00340 HepSymMatrix operator/(const HepSymMatrix &m1, double t); 00341 // s = s1 / t. (s /= t is faster if you can use it.) 00342 00343 HepMatrix operator+(const HepMatrix &m1, const HepSymMatrix &s2); 00344 HepMatrix operator+(const HepSymMatrix &s1, const HepMatrix &m2); 00345 HepSymMatrix operator+(const HepSymMatrix &s1, const HepSymMatrix &s2); 00346 // Addition operators 00347 00348 HepMatrix operator-(const HepMatrix &m1, const HepSymMatrix &s2); 00349 HepMatrix operator-(const HepSymMatrix &m1, const HepMatrix &m2); 00350 HepSymMatrix operator-(const HepSymMatrix &s1, const HepSymMatrix &s2); 00351 // subtraction operators 00352 00353 HepSymMatrix dsum(const HepSymMatrix &s1, const HepSymMatrix &s2); 00354 // Direct sum of two symmetric matrices; 00355 00356 double condition(const HepSymMatrix &m); 00357 // Find the conditon number of a symmetric matrix. 00358 00359 void diag_step(HepSymMatrix *t, int begin, int end); 00360 void diag_step(HepSymMatrix *t, HepMatrix *u, int begin, int end); 00361 // Implicit symmetric QR step with Wilkinson Shift 00362 00363 HepMatrix diagonalize(HepSymMatrix *s); 00364 // Diagonalize a symmetric matrix. 00365 // It returns the matrix U so that s_old = U * s_diag * U.T() 00366 00367 HepVector house(const HepSymMatrix &a, int row=1, int col=1); 00368 void house_with_update2(HepSymMatrix *a, HepMatrix *v, int row=1, int col=1); 00369 // Finds and does Householder reflection on matrix. 00370 00371 void tridiagonal(HepSymMatrix *a, HepMatrix *hsm); 00372 HepMatrix tridiagonal(HepSymMatrix *a); 00373 // Does a Householder tridiagonalization of a symmetric matrix. 00374 00375 } // namespace CLHEP 00376 00377 #ifdef ENABLE_BACKWARDS_COMPATIBILITY 00378 // backwards compatibility will be enabled ONLY in CLHEP 1.9 00379 using namespace CLHEP; 00380 #endif 00381 00382 #ifndef HEP_DEBUG_INLINE 00383 #include "CLHEP/Matrix/SymMatrix.icc" 00384 #endif 00385 00386 #endif