CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

Matrix/CLHEP/Matrix/DiagMatrix.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 // CLASSDOC OFF
3 // ---------------------------------------------------------------------------
4 // CLASSDOC ON
5 //
6 // This file is a part of the CLHEP - a Class Library for High Energy Physics.
7 //
8 // This software written by Nobu Katayama and Mike Smyth, Cornell University.
9 //
10 // DiagMatrix is a class for diagonal matrix. This is useful for a covariance
11 // matrix of measured quantities since they are uncorrelated to each other
12 // and therefore diagonal. It is obviously smaller and faster to manipulate.
13 //
14 
15 #ifndef _DIAGMatrix_H_
16 #define _DIAGMatrix_H_
17 
18 #ifdef GNUPRAGMA
19 #pragma interface
20 #endif
21 
22 #include <vector>
23 
24 #include "CLHEP/Matrix/defs.h"
25 #include "CLHEP/Matrix/GenMatrix.h"
26 
27 namespace CLHEP {
28 
29 class HepRandom;
30 
31 class HepMatrix;
32 class HepSymMatrix;
33 class HepVector;
34 
39 class HepDiagMatrix: public HepGenMatrix {
40 public:
41  inline HepDiagMatrix();
42  // Default constructor. Gives 0x0 matrix. Another Matrix can be assigned
43  // to it.
44 
45  explicit HepDiagMatrix(int p);
46  HepDiagMatrix(int p, int);
47  // Constructor. Gives p x p diagonal matrix.
48  // With a second argument, either 0 or 1, the matrix is initialized.
49  // 0 means a zero matrix, 1 means the identity matrix.
50 
51  HepDiagMatrix(int p, HepRandom &r);
52 
53  HepDiagMatrix(const HepDiagMatrix &m1);
54  // Copy constructor.
55 
56  virtual ~HepDiagMatrix();
57  // Destructor.
58 
59  inline int num_row() const;
60  inline int num_col() const;
61  // Returns the number of rows/columns. (Should be equal.)
62 
63  double &operator()(int row, int col);
64  const double &operator()(int row, int col) const;
65  // Read or write a matrix element. row must be equal to col.
66  // ** Note that indexing starts from (1,1). **
67 
68  double &fast(int row, int col);
69  const double &fast(int row, int col) const;
70  // fast element access.
71  // Must be row>=col;
72  // ** Note that indexing starts from (1,1). **
73 
74  void assign(const HepMatrix &m2);
75  // Assigns m2 to d, assuming m2 is a diagnal matrix.
76 
77  void assign(const HepSymMatrix &m2);
78  // Assigns m2 to d, assuming m2 is a diagnal matrix.
79 
80  void assign(const HepDiagMatrix &m2);
81  // Another form of assignment. For consistency.
82 
83  HepDiagMatrix & operator*=(double t);
84  // Multiply a DiagMatrix by a floating number
85 
86  HepDiagMatrix & operator/=(double t);
87  // Divide a DiagMatrix by a floating number
88 
91  // Add or subtract a DiagMatrix.
92 
94  // Assignment operator. To assign SymMatrix to DiagMatrix, use d<<s.
95 
96  HepDiagMatrix operator- () const;
97  // unary minus, ie. flip the sign of each element.
98 
99  HepDiagMatrix T() const;
100  // Returns the transpose of a DiagMatrix (which is itself).
101 
102  HepDiagMatrix apply(double (*f)(double,
103  int, int)) const;
104  // Apply a function to all elements of the matrix.
105 
106  HepSymMatrix similarity(const HepMatrix &m1) const;
107  // Returns m1*s*m1.T().
108  HepSymMatrix similarityT(const HepMatrix &m1) const;
109  // Returns m1.T()*s*m1.
110 
111  double similarity(const HepVector &) const;
112  // Returns v.T()*s*v (This is a scaler).
113 
114  HepDiagMatrix sub(int min_row, int max_row) const;
115  // Returns a sub matrix of a SymMatrix.
116  HepDiagMatrix sub(int min_row, int max_row);
117  // SGI CC bug. I have to have both with/without const. I should not need
118  // one without const.
119 
120  void sub(int row, const HepDiagMatrix &m1);
121  // Sub matrix of this SymMatrix is replaced with m1.
122 
123  HepDiagMatrix inverse(int&ierr) const;
124  // Invert a Matrix. The matrix is not changed
125  // Returns 0 when successful, otherwise non-zero.
126 
127  void invert(int&ierr);
128  // Invert a Matrix.
129  // N.B. the contents of the matrix are replaced by the inverse.
130  // Returns ierr = 0 when successful, otherwise non-zero.
131  // This method has less overhead then inverse().
132 
133  double determinant() const;
134  // calculate the determinant of the matrix.
135 
136  double trace() const;
137  // calculate the trace of the matrix (sum of diagonal elements).
138 
140  public:
141  inline HepDiagMatrix_row(HepDiagMatrix&,int);
142  inline double & operator[](int);
143  private:
144  HepDiagMatrix& _a;
145  int _r;
146  };
148  public:
149  inline HepDiagMatrix_row_const(const HepDiagMatrix&,int);
150  inline const double & operator[](int) const;
151  private:
152  const HepDiagMatrix& _a;
153  int _r;
154  };
155  // helper classes to implement m[i][j]
156 
157  inline HepDiagMatrix_row operator[] (int);
158  inline HepDiagMatrix_row_const operator[] (int) const;
159  // Read or write a matrix element.
160  // While it may not look like it, you simply do m[i][j] to get an
161  // element.
162  // ** Note that the indexing starts from [0][0]. **
163 
164 protected:
165  inline int num_size() const;
166 
167 private:
168  friend class HepDiagMatrix_row;
170  friend class HepMatrix;
171  friend class HepSymMatrix;
172 
173  friend HepDiagMatrix operator*(const HepDiagMatrix &m1,
174  const HepDiagMatrix &m2);
175  friend HepDiagMatrix operator+(const HepDiagMatrix &m1,
176  const HepDiagMatrix &m2);
177  friend HepDiagMatrix operator-(const HepDiagMatrix &m1,
178  const HepDiagMatrix &m2);
179  friend HepMatrix operator*(const HepDiagMatrix &m1, const HepMatrix &m2);
180  friend HepMatrix operator*(const HepMatrix &m1, const HepDiagMatrix &m2);
181  friend HepVector operator*(const HepDiagMatrix &m1, const HepVector &m2);
182 
183 #ifdef DISABLE_ALLOC
184  std::vector<double > m;
185 #else
186  std::vector<double,Alloc<double,25> > m;
187 #endif
188  int nrow;
189 #if defined(__sun) || !defined(__GNUG__)
190 //
191 // Sun CC 4.0.1 has this bug.
192 //
193  static double zero;
194 #else
195  static const double zero;
196 #endif
197 };
198 
199 std::ostream& operator<<(std::ostream &s, const HepDiagMatrix &q);
200 // Write out Matrix, SymMatrix, DiagMatrix and Vector into ostream.
201 
202 HepMatrix operator*(const HepMatrix &m1, const HepDiagMatrix &m2);
203 HepMatrix operator*(const HepDiagMatrix &m1, const HepMatrix &m2);
204 HepDiagMatrix operator*(double t, const HepDiagMatrix &d1);
205 HepDiagMatrix operator*(const HepDiagMatrix &d1, double t);
206 // Multiplication operators
207 // Note that m *= m1 is always faster than m = m * m1
208 
209 HepDiagMatrix operator/(const HepDiagMatrix &m1, double t);
210 // d = d1 / t. (d /= t is faster if you can use it.)
211 
212 HepMatrix operator+(const HepMatrix &m1, const HepDiagMatrix &d2);
213 HepMatrix operator+(const HepDiagMatrix &d1, const HepMatrix &m2);
214 HepDiagMatrix operator+(const HepDiagMatrix &m1, const HepDiagMatrix &d2);
215 HepSymMatrix operator+(const HepSymMatrix &s1, const HepDiagMatrix &d2);
216 HepSymMatrix operator+(const HepDiagMatrix &d1, const HepSymMatrix &s2);
217 // Addition operators
218 
219 HepMatrix operator-(const HepMatrix &m1, const HepDiagMatrix &d2);
220 HepMatrix operator-(const HepDiagMatrix &d1, const HepMatrix &m2);
221 HepDiagMatrix operator-(const HepDiagMatrix &d1, const HepDiagMatrix &d2);
222 HepSymMatrix operator-(const HepSymMatrix &s1, const HepDiagMatrix &d2);
223 HepSymMatrix operator-(const HepDiagMatrix &d1, const HepSymMatrix &s2);
224 // Subtraction operators
225 
226 HepDiagMatrix dsum(const HepDiagMatrix &s1, const HepDiagMatrix &s2);
227 // Direct sum of two diagonal matricies;
228 
229 } // namespace CLHEP
230 
231 #ifdef ENABLE_BACKWARDS_COMPATIBILITY
232 // backwards compatibility will be enabled ONLY in CLHEP 1.9
233 using namespace CLHEP;
234 #endif
235 
236 #ifndef HEP_DEBUG_INLINE
237 #include "CLHEP/Matrix/DiagMatrix.icc"
238 #endif
239 
240 #endif