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