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

Matrix/Matrix/SymMatrix.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 is the definition of the HepSymMatrix class.
9 //
10 // This software written by Nobu Katayama and Mike Smyth, Cornell University.
11 //
12 // .SS Usage
13 //
14 // This is very much like the Matrix, except of course it is restricted to
15 // Symmetric Matrix. All the operations for Matrix can also be done here
16 // (except for the +=,-=,*= that don't yield a symmetric matrix. e.g.
17 // +=(const Matrix &) is not defined)
18 
19 // The matrix is stored as a lower triangular matrix.
20 // In addition to the (row, col) method of finding element, fast(row, col)
21 // returns an element with checking to see if row and col need to be
22 // interchanged so that row >= col.
23 
24 // New operations are:
25 //
26 // .ft B
27 // sym = s.similarity(m);
28 //
29 // This returns m*s*m.T(). This is a similarity
30 // transform when m is orthogonal, but nothing
31 // restricts m to be orthogonal. It is just
32 // a more efficient way to calculate m*s*m.T,
33 // and it realizes that this should be a
34 // HepSymMatrix (the explicit operation m*s*m.T
35 // will return a Matrix, not realizing that
36 // it is symmetric).
37 //
38 // .ft B
39 // sym = similarityT(m);
40 //
41 // This returns m.T()*s*m.
42 //
43 // .ft B
44 // s << m;
45 //
46 // This takes the matrix m, and treats it
47 // as symmetric matrix that is copied to s.
48 // This is useful for operations that yield
49 // symmetric matrix, but which the computer
50 // is too dumb to realize.
51 //
52 // .ft B
53 // s = vT_times_v(const HepVector v);
54 //
55 // calculates v.T()*v.
56 //
57 // ./"This code has been written by Mike Smyth, and the algorithms used are
58 // ./"described in the thesis "A Tracking Library for a Silicon Vertex Detector"
59 // ./"(Mike Smyth, Cornell University, June 1993).
60 // ./"This is file contains C++ stuff for doing things with Matrixes.
61 // ./"To turn on bound checking, define MATRIX_BOUND_CHECK before including
62 // ./"this file.
63 //
64 
65 #ifndef _SYMMatrix_H_
66 #define _SYMMatrix_H_
67 
68 #ifdef GNUPRAGMA
69 #pragma interface
70 #endif
71 
72 #include <vector>
73 
74 #include "CLHEP/Matrix/defs.h"
75 #include "CLHEP/Matrix/GenMatrix.h"
76 
77 namespace CLHEP {
78 
79 class HepRandom;
80 
81 class HepMatrix;
82 class HepDiagMatrix;
83 class HepVector;
84 
89 class HepSymMatrix : public HepGenMatrix {
90 public:
91  inline HepSymMatrix();
92  // Default constructor. Gives 0x0 symmetric matrix.
93  // Another SymMatrix can be assigned to it.
94 
95  explicit HepSymMatrix(int p);
96  HepSymMatrix(int p, int);
97  // Constructor. Gives p x p symmetric matrix.
98  // With a second argument, the matrix is initialized. 0 means a zero
99  // matrix, 1 means the identity matrix.
100 
101  HepSymMatrix(int p, HepRandom &r);
102 
103  HepSymMatrix(const HepSymMatrix &hm1);
104  // Copy constructor.
105 
106  HepSymMatrix(const HepDiagMatrix &hm1);
107  // Constructor from DiagMatrix
108 
109  virtual ~HepSymMatrix();
110  // Destructor.
111 
112  inline int num_row() const;
113  inline int num_col() const;
114  // Returns number of rows/columns.
115 
116  const double & operator()(int row, int col) const;
117  double & operator()(int row, int col);
118  // Read and write a SymMatrix element.
119  // ** Note that indexing starts from (1,1). **
120 
121  const double & fast(int row, int col) const;
122  double & fast(int row, int col);
123  // fast element access.
124  // Must be row>=col;
125  // ** Note that indexing starts from (1,1). **
126 
127  void assign(const HepMatrix &hm2);
128  // Assigns hm2 to s, assuming hm2 is a symmetric matrix.
129 
130  void assign(const HepSymMatrix &hm2);
131  // Another form of assignment. For consistency.
132 
133  HepSymMatrix & operator*=(double t);
134  // Multiply a SymMatrix by a floating number.
135 
136  HepSymMatrix & operator/=(double t);
137  // Divide a SymMatrix by a floating number.
138 
139  HepSymMatrix & operator+=( const HepSymMatrix &hm2);
140  HepSymMatrix & operator+=( const HepDiagMatrix &hm2);
141  HepSymMatrix & operator-=( const HepSymMatrix &hm2);
142  HepSymMatrix & operator-=( const HepDiagMatrix &hm2);
143  // Add or subtract a SymMatrix.
144 
145  HepSymMatrix & operator=( const HepSymMatrix &hm2);
146  HepSymMatrix & operator=( const HepDiagMatrix &hm2);
147  // Assignment operators. Notice that there is no SymMatrix = Matrix.
148 
149  HepSymMatrix operator- () const;
150  // unary minus, ie. flip the sign of each element.
151 
152  HepSymMatrix T() const;
153  // Returns the transpose of a SymMatrix (which is itself).
154 
155  HepSymMatrix apply(double (*f)(double, int, int)) const;
156  // Apply a function to all elements of the matrix.
157 
158  HepSymMatrix similarity(const HepMatrix &hm1) const;
159  HepSymMatrix similarity(const HepSymMatrix &hm1) const;
160  // Returns hm1*s*hm1.T().
161 
162  HepSymMatrix similarityT(const HepMatrix &hm1) const;
163  // temporary. test of new similarity.
164  // Returns hm1.T()*s*hm1.
165 
166  double similarity(const HepVector &v) const;
167  // Returns v.T()*s*v (This is a scaler).
168 
169  HepSymMatrix sub(int min_row, int max_row) const;
170  // Returns a sub matrix of a SymMatrix.
171  void sub(int row, const HepSymMatrix &hm1);
172  // Sub matrix of this SymMatrix is replaced with hm1.
173  HepSymMatrix sub(int min_row, int max_row);
174  // SGI CC bug. I have to have both with/without const. I should not need
175  // one without const.
176 
177  inline HepSymMatrix inverse(int &ifail) const;
178  // Invert a Matrix. The matrix is not changed
179  // Returns 0 when successful, otherwise non-zero.
180 
181  void invert(int &ifail);
182  // Invert a Matrix.
183  // N.B. the contents of the matrix are replaced by the inverse.
184  // Returns ierr = 0 when successful, otherwise non-zero.
185  // This method has less overhead then inverse().
186 
187  inline void invert();
188  // Invert a matrix. Throw std::runtime_error on failure.
189 
190  inline HepSymMatrix inverse() const;
191  // Invert a matrix. Throw std::runtime_error on failure.
192 
193  double determinant() const;
194  // calculate the determinant of the matrix.
195 
196  double trace() const;
197  // calculate the trace of the matrix (sum of diagonal elements).
198 
199  class HepSymMatrix_row {
200  public:
201  inline HepSymMatrix_row(HepSymMatrix&,int);
202  inline double & operator[](int);
203  private:
204  HepSymMatrix& _a;
205  int _r;
206  };
207  class HepSymMatrix_row_const {
208  public:
209  inline HepSymMatrix_row_const(const HepSymMatrix&,int);
210  inline const double & operator[](int) const;
211  private:
212  const HepSymMatrix& _a;
213  int _r;
214  };
215  // helper class to implement m[i][j]
216 
217  inline HepSymMatrix_row operator[] (int);
218  inline HepSymMatrix_row_const operator[] (int) const;
219  // Read or write a matrix element.
220  // While it may not look like it, you simply do m[i][j] to get an
221  // element.
222  // ** Note that the indexing starts from [0][0]. **
223 
224  // Special-case inversions for 5x5 and 6x6 symmetric positive definite:
225  // These set ifail=0 and invert if the matrix was positive definite;
226  // otherwise ifail=1 and the matrix is left unaltered.
227  void invertCholesky5 (int &ifail);
228  void invertCholesky6 (int &ifail);
229 
230  // Inversions for 5x5 and 6x6 forcing use of specific methods: The
231  // behavior (though not the speed) will be identical to invert(ifail).
232  void invertHaywood4 (int & ifail);
233  void invertHaywood5 (int &ifail);
234  void invertHaywood6 (int &ifail);
235  void invertBunchKaufman (int &ifail);
236 
237 protected:
238  inline int num_size() const;
239 
240 private:
241  friend class HepSymMatrix_row;
242  friend class HepSymMatrix_row_const;
243  friend class HepMatrix;
244  friend class HepDiagMatrix;
245 
246  friend void tridiagonal(HepSymMatrix *a,HepMatrix *hsm);
247  friend double condition(const HepSymMatrix &m);
248  friend void diag_step(HepSymMatrix *t,int begin,int end);
249  friend void diag_step(HepSymMatrix *t,HepMatrix *u,int begin,int end);
251  friend HepVector house(const HepSymMatrix &a,int row,int col);
252  friend void house_with_update2(HepSymMatrix *a,HepMatrix *v,int row,int col);
253 
254  friend HepSymMatrix operator+(const HepSymMatrix &hm1,
255  const HepSymMatrix &hm2);
256  friend HepSymMatrix operator-(const HepSymMatrix &hm1,
257  const HepSymMatrix &hm2);
258  friend HepMatrix operator*(const HepSymMatrix &hm1, const HepSymMatrix &hm2);
259  friend HepMatrix operator*(const HepSymMatrix &hm1, const HepMatrix &hm2);
260  friend HepMatrix operator*(const HepMatrix &hm1, const HepSymMatrix &hm2);
261  friend HepVector operator*(const HepSymMatrix &hm1, const HepVector &hm2);
262  // Multiply a Matrix by a Matrix or Vector.
263 
264  friend HepSymMatrix vT_times_v(const HepVector &v);
265  // Returns v * v.T();
266 
267 #ifdef DISABLE_ALLOC
268  std::vector<double > m;
269 #else
270  std::vector<double,Alloc<double,25> > m;
271 #endif
272  int nrow;
273  int size_; // total number of elements
274 
275  static double posDefFraction5x5;
276  static double adjustment5x5;
277  static const double CHOLESKY_THRESHOLD_5x5;
278  static const double CHOLESKY_CREEP_5x5;
279 
280  static double posDefFraction6x6;
281  static double adjustment6x6;
282  static const double CHOLESKY_THRESHOLD_6x6;
283  static const double CHOLESKY_CREEP_6x6;
284 
285  void invert4 (int & ifail);
286  void invert5 (int & ifail);
287  void invert6 (int & ifail);
288 
289 };
290 
291 //
292 // Operations other than member functions for Matrix, SymMatrix, DiagMatrix
293 // and Vectors implemented in Matrix.cc and Matrix.icc (inline).
294 //
295 
296 std::ostream& operator<<(std::ostream &s, const HepSymMatrix &q);
297 // Write out Matrix, SymMatrix, DiagMatrix and Vector into ostream.
298 
299 HepMatrix operator*(const HepMatrix &hm1, const HepSymMatrix &hm2);
300 HepMatrix operator*(const HepSymMatrix &hm1, const HepMatrix &hm2);
301 HepMatrix operator*(const HepSymMatrix &hm1, const HepSymMatrix &hm2);
302 HepSymMatrix operator*(double t, const HepSymMatrix &s1);
303 HepSymMatrix operator*(const HepSymMatrix &s1, double t);
304 // Multiplication operators.
305 // Note that m *= hm1 is always faster than m = m * hm1
306 
307 HepSymMatrix operator/(const HepSymMatrix &hm1, double t);
308 // s = s1 / t. (s /= t is faster if you can use it.)
309 
310 HepMatrix operator+(const HepMatrix &hm1, const HepSymMatrix &s2);
311 HepMatrix operator+(const HepSymMatrix &s1, const HepMatrix &hm2);
312 HepSymMatrix operator+(const HepSymMatrix &s1, const HepSymMatrix &s2);
313 // Addition operators
314 
315 HepMatrix operator-(const HepMatrix &hm1, const HepSymMatrix &s2);
316 HepMatrix operator-(const HepSymMatrix &hm1, const HepMatrix &hm2);
317 HepSymMatrix operator-(const HepSymMatrix &s1, const HepSymMatrix &s2);
318 // subtraction operators
319 
320 HepSymMatrix dsum(const HepSymMatrix &s1, const HepSymMatrix &s2);
321 // Direct sum of two symmetric matrices;
322 
323 double condition(const HepSymMatrix &m);
324 // Find the conditon number of a symmetric matrix.
325 
326 void diag_step(HepSymMatrix *t, int begin, int end);
327 void diag_step(HepSymMatrix *t, HepMatrix *u, int begin, int end);
328 // Implicit symmetric QR step with Wilkinson Shift
329 
330 HepMatrix diagonalize(HepSymMatrix *s);
331 // Diagonalize a symmetric matrix.
332 // It returns the matrix U so that s_old = U * s_diag * U.T()
333 
334 HepVector house(const HepSymMatrix &a, int row=1, int col=1);
335 void house_with_update2(HepSymMatrix *a, HepMatrix *v, int row=1, int col=1);
336 // Finds and does Householder reflection on matrix.
337 
338 void tridiagonal(HepSymMatrix *a, HepMatrix *hsm);
339 HepMatrix tridiagonal(HepSymMatrix *a);
340 // Does a Householder tridiagonalization of a symmetric matrix.
341 
342 } // namespace CLHEP
343 
344 #ifdef ENABLE_BACKWARDS_COMPATIBILITY
345 // backwards compatibility will be enabled ONLY in CLHEP 1.9
346 using namespace CLHEP;
347 #endif
348 
349 #ifndef HEP_DEBUG_INLINE
350 #include "CLHEP/Matrix/SymMatrix.icc"
351 #endif
352 
353 #endif
friend void diag_step(HepSymMatrix *t, int begin, int end)
Hep3Vector operator+(const Hep3Vector &, const Hep3Vector &)
int num_col() const
void invertBunchKaufman(int &ifail)
Definition: SymMatrix.cc:964
void house_with_update2(HepSymMatrix *a, HepMatrix *v, int row=1, int col=1)
double condition(const HepSymMatrix &m)
HepSymMatrix_row(HepSymMatrix &, int)
HepLorentzVector operator/(const HepLorentzVector &, double a)
HepSymMatrix_row operator[](int)
HepSymMatrix apply(double(*f)(double, int, int)) const
Definition: SymMatrix.cc:700
void invertHaywood5(int &ifail)
HepSymMatrix inverse() const
const double & fast(int row, int col) const
HepSymMatrix T() const
friend HepMatrix operator*(const HepSymMatrix &hm1, const HepSymMatrix &hm2)
Definition: SymMatrix.cc:437
friend void tridiagonal(HepSymMatrix *a, HepMatrix *hsm)
double determinant() const
Definition: SymMatrix.cc:943
virtual ~HepSymMatrix()
Definition: SymMatrix.cc:103
HepVector house(const HepMatrix &a, int row=1, int col=1)
HepLorentzRotation operator*(const HepRotation &r, const HepLorentzRotation &lt)
double trace() const
Definition: SymMatrix.cc:957
HepDiagMatrix dsum(const HepDiagMatrix &s1, const HepDiagMatrix &s2)
Definition: DiagMatrix.cc:164
int num_size() const
friend HepSymMatrix operator+(const HepSymMatrix &hm1, const HepSymMatrix &hm2)
Definition: SymMatrix.cc:256
void invertCholesky6(int &ifail)
HepSymMatrix & operator*=(double t)
Definition: SymMatrix.cc:614
int num_row() const
HepSymMatrix similarityT(const HepMatrix &hm1) const
Definition: SymMatrix.cc:816
HepSymMatrix & operator-=(const HepSymMatrix &hm2)
Definition: SymMatrix.cc:601
const double & operator()(int row, int col) const
HepSymMatrix_row_const(const HepSymMatrix &, int)
HepSymMatrix similarity(const HepMatrix &hm1) const
Definition: SymMatrix.cc:737
friend HepVector house(const HepSymMatrix &a, int row, int col)
HepSymMatrix & operator+=(const HepSymMatrix &hm2)
Definition: SymMatrix.cc:578
void invertHaywood6(int &ifail)
HepSymMatrix & operator/=(double t)
Definition: SymMatrix.cc:608
void diag_step(HepSymMatrix *t, int begin, int end)
void f(void g())
Definition: excDblThrow.cc:38
friend HepSymMatrix vT_times_v(const HepVector &v)
Definition: SymMatrix.cc:542
Hep3Vector operator-(const Hep3Vector &, const Hep3Vector &)
friend void house_with_update2(HepSymMatrix *a, HepMatrix *v, int row, int col)
friend HepMatrix diagonalize(HepSymMatrix *s)
const double & operator[](int) const
void assign(const HepMatrix &hm2)
Definition: SymMatrix.cc:718
HepSymMatrix sub(int min_row, int max_row) const
Definition: SymMatrix.cc:134
friend double condition(const HepSymMatrix &m)
HepSymMatrix operator-() const
Definition: SymMatrix.cc:214
void invertCholesky5(int &ifail)
HepMatrix diagonalize(HepSymMatrix *s)
std::ostream & operator<<(std::ostream &os, const HepAxisAngle &aa)
Definition: AxisAngle.cc:86
void tridiagonal(HepSymMatrix *a, HepMatrix *hsm)
void invertHaywood4(int &ifail)
HepSymMatrix & operator=(const HepSymMatrix &hm2)
Definition: SymMatrix.cc:645