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 software written by Nobu Katayama and Mike Smyth, Cornell University. 00009 // 00010 // DiagMatrix is a class for diagonal matrix. This is useful for a covariance 00011 // matrix of measured quantities since they are uncorrelated to each other 00012 // and therefore diagonal. It is obviously smaller and faster to manipulate. 00013 // 00014 00015 #ifndef _DIAGMatrix_H_ 00016 #define _DIAGMatrix_H_ 00017 00018 #ifdef GNUPRAGMA 00019 #pragma interface 00020 #endif 00021 00022 #include <vector> 00023 00024 #include "CLHEP/Matrix/defs.h" 00025 #include "CLHEP/Matrix/GenMatrix.h" 00026 00027 namespace CLHEP { 00028 00029 class HepRandom; 00030 00031 class HepMatrix; 00032 class HepSymMatrix; 00033 class HepVector; 00034 00039 class HepDiagMatrix: public HepGenMatrix { 00040 public: 00041 inline HepDiagMatrix(); 00042 // Default constructor. Gives 0x0 matrix. Another Matrix can be assigned 00043 // to it. 00044 00045 explicit HepDiagMatrix(int p); 00046 HepDiagMatrix(int p, int); 00047 // Constructor. Gives p x p diagonal matrix. 00048 // With a second argument, either 0 or 1, the matrix is initialized. 00049 // 0 means a zero matrix, 1 means the identity matrix. 00050 00051 HepDiagMatrix(int p, HepRandom &r); 00052 00053 HepDiagMatrix(const HepDiagMatrix &m1); 00054 // Copy constructor. 00055 00056 virtual ~HepDiagMatrix(); 00057 // Destructor. 00058 00059 inline int num_row() const; 00060 inline int num_col() const; 00061 // Returns the number of rows/columns. (Should be equal.) 00062 00063 double &operator()(int row, int col); 00064 const double &operator()(int row, int col) const; 00065 // Read or write a matrix element. row must be equal to col. 00066 // ** Note that indexing starts from (1,1). ** 00067 00068 double &fast(int row, int col); 00069 const double &fast(int row, int col) const; 00070 // fast element access. 00071 // Must be row>=col; 00072 // ** Note that indexing starts from (1,1). ** 00073 00074 void assign(const HepMatrix &m2); 00075 // Assigns m2 to d, assuming m2 is a diagnal matrix. 00076 00077 void assign(const HepSymMatrix &m2); 00078 // Assigns m2 to d, assuming m2 is a diagnal matrix. 00079 00080 void assign(const HepDiagMatrix &m2); 00081 // Another form of assignment. For consistency. 00082 00083 HepDiagMatrix & operator*=(double t); 00084 // Multiply a DiagMatrix by a floating number 00085 00086 HepDiagMatrix & operator/=(double t); 00087 // Divide a DiagMatrix by a floating number 00088 00089 HepDiagMatrix & operator+=( const HepDiagMatrix &m2); 00090 HepDiagMatrix & operator-=( const HepDiagMatrix &m2); 00091 // Add or subtract a DiagMatrix. 00092 00093 HepDiagMatrix & operator=( const HepDiagMatrix &m2); 00094 // Assignment operator. To assign SymMatrix to DiagMatrix, use d<<s. 00095 00096 HepDiagMatrix operator- () const; 00097 // unary minus, ie. flip the sign of each element. 00098 00099 HepDiagMatrix T() const; 00100 // Returns the transpose of a DiagMatrix (which is itself). 00101 00102 HepDiagMatrix apply(double (*f)(double, 00103 int, int)) const; 00104 // Apply a function to all elements of the matrix. 00105 00106 HepSymMatrix similarity(const HepMatrix &m1) const; 00107 // Returns m1*s*m1.T(). 00108 HepSymMatrix similarityT(const HepMatrix &m1) const; 00109 // Returns m1.T()*s*m1. 00110 00111 double similarity(const HepVector &) const; 00112 // Returns v.T()*s*v (This is a scaler). 00113 00114 HepDiagMatrix sub(int min_row, int max_row) const; 00115 // Returns a sub matrix of a SymMatrix. 00116 HepDiagMatrix sub(int min_row, int max_row); 00117 // SGI CC bug. I have to have both with/without const. I should not need 00118 // one without const. 00119 00120 void sub(int row, const HepDiagMatrix &m1); 00121 // Sub matrix of this SymMatrix is replaced with m1. 00122 00123 HepDiagMatrix inverse(int&ierr) const; 00124 // Invert a Matrix. The matrix is not changed 00125 // Returns 0 when successful, otherwise non-zero. 00126 00127 void invert(int&ierr); 00128 // Invert a Matrix. 00129 // N.B. the contents of the matrix are replaced by the inverse. 00130 // Returns ierr = 0 when successful, otherwise non-zero. 00131 // This method has less overhead then inverse(). 00132 00133 inline void invert(); 00134 // Invert a matrix. Throw std::runtime_error on failure. 00135 00136 inline HepDiagMatrix inverse() const; 00137 // Invert a matrix. Throw std::runtime_error on failure. 00138 00139 double determinant() const; 00140 // calculate the determinant of the matrix. 00141 00142 double trace() const; 00143 // calculate the trace of the matrix (sum of diagonal elements). 00144 00145 class HepDiagMatrix_row { 00146 public: 00147 inline HepDiagMatrix_row(HepDiagMatrix&,int); 00148 inline double & operator[](int); 00149 private: 00150 HepDiagMatrix& _a; 00151 int _r; 00152 }; 00153 class HepDiagMatrix_row_const { 00154 public: 00155 inline HepDiagMatrix_row_const(const HepDiagMatrix&,int); 00156 inline const double & operator[](int) const; 00157 private: 00158 const HepDiagMatrix& _a; 00159 int _r; 00160 }; 00161 // helper classes to implement m[i][j] 00162 00163 inline HepDiagMatrix_row operator[] (int); 00164 inline HepDiagMatrix_row_const operator[] (int) const; 00165 // Read or write a matrix element. 00166 // While it may not look like it, you simply do m[i][j] to get an 00167 // element. 00168 // ** Note that the indexing starts from [0][0]. ** 00169 00170 protected: 00171 inline int num_size() const; 00172 00173 private: 00174 friend class HepDiagMatrix_row; 00175 friend class HepDiagMatrix_row_const; 00176 friend class HepMatrix; 00177 friend class HepSymMatrix; 00178 00179 friend HepDiagMatrix operator*(const HepDiagMatrix &m1, 00180 const HepDiagMatrix &m2); 00181 friend HepDiagMatrix operator+(const HepDiagMatrix &m1, 00182 const HepDiagMatrix &m2); 00183 friend HepDiagMatrix operator-(const HepDiagMatrix &m1, 00184 const HepDiagMatrix &m2); 00185 friend HepMatrix operator*(const HepDiagMatrix &m1, const HepMatrix &m2); 00186 friend HepMatrix operator*(const HepMatrix &m1, const HepDiagMatrix &m2); 00187 friend HepVector operator*(const HepDiagMatrix &m1, const HepVector &m2); 00188 00189 #ifdef DISABLE_ALLOC 00190 std::vector<double > m; 00191 #else 00192 std::vector<double,Alloc<double,25> > m; 00193 #endif 00194 int nrow; 00195 #if defined(__sun) || !defined(__GNUG__) 00196 // 00197 // Sun CC 4.0.1 has this bug. 00198 // 00199 static double zero; 00200 #else 00201 static const double zero; 00202 #endif 00203 }; 00204 00205 std::ostream& operator<<(std::ostream &s, const HepDiagMatrix &q); 00206 // Write out Matrix, SymMatrix, DiagMatrix and Vector into ostream. 00207 00208 HepMatrix operator*(const HepMatrix &m1, const HepDiagMatrix &m2); 00209 HepMatrix operator*(const HepDiagMatrix &m1, const HepMatrix &m2); 00210 HepDiagMatrix operator*(double t, const HepDiagMatrix &d1); 00211 HepDiagMatrix operator*(const HepDiagMatrix &d1, double t); 00212 // Multiplication operators 00213 // Note that m *= m1 is always faster than m = m * m1 00214 00215 HepDiagMatrix operator/(const HepDiagMatrix &m1, double t); 00216 // d = d1 / t. (d /= t is faster if you can use it.) 00217 00218 HepMatrix operator+(const HepMatrix &m1, const HepDiagMatrix &d2); 00219 HepMatrix operator+(const HepDiagMatrix &d1, const HepMatrix &m2); 00220 HepDiagMatrix operator+(const HepDiagMatrix &m1, const HepDiagMatrix &d2); 00221 HepSymMatrix operator+(const HepSymMatrix &s1, const HepDiagMatrix &d2); 00222 HepSymMatrix operator+(const HepDiagMatrix &d1, const HepSymMatrix &s2); 00223 // Addition operators 00224 00225 HepMatrix operator-(const HepMatrix &m1, const HepDiagMatrix &d2); 00226 HepMatrix operator-(const HepDiagMatrix &d1, const HepMatrix &m2); 00227 HepDiagMatrix operator-(const HepDiagMatrix &d1, const HepDiagMatrix &d2); 00228 HepSymMatrix operator-(const HepSymMatrix &s1, const HepDiagMatrix &d2); 00229 HepSymMatrix operator-(const HepDiagMatrix &d1, const HepSymMatrix &s2); 00230 // Subtraction operators 00231 00232 HepDiagMatrix dsum(const HepDiagMatrix &s1, const HepDiagMatrix &s2); 00233 // Direct sum of two diagonal matricies; 00234 00235 } // namespace CLHEP 00236 00237 #ifdef ENABLE_BACKWARDS_COMPATIBILITY 00238 // backwards compatibility will be enabled ONLY in CLHEP 1.9 00239 using namespace CLHEP; 00240 #endif 00241 00242 #ifndef HEP_DEBUG_INLINE 00243 #include "CLHEP/Matrix/DiagMatrix.icc" 00244 #endif 00245 00246 #endif