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
CLHEP::HepSymMatrix::operator[]
HepSymMatrix_row operator[](int)
CLHEP::house
HepVector house(const HepMatrix &a, int row=1, int col=1)
Definition: MatrixLinear.cc:368
a
@ a
Definition: testCategories.cc:125
CLHEP::HepSymMatrix::operator+
friend HepSymMatrix operator+(const HepSymMatrix &hm1, const HepSymMatrix &hm2)
Definition: SymMatrix.cc:256
CLHEP::HepSymMatrix::invertCholesky6
void invertCholesky6(int &ifail)
Definition: SymMatrixInvert.cc:806
CLHEP::HepSymMatrix::operator/=
HepSymMatrix & operator/=(double t)
Definition: SymMatrix.cc:608
CLHEP::HepSymMatrix::diag_step
friend void diag_step(HepSymMatrix *t, int begin, int end)
Definition: MatrixLinear.cc:224
CLHEP::condition
double condition(const HepSymMatrix &m)
Definition: MatrixLinear.cc:198
CLHEP::HepSymMatrix::invert
void invert()
CLHEP::HepSymMatrix::operator*=
HepSymMatrix & operator*=(double t)
Definition: SymMatrix.cc:614
CLHEP::HepSymMatrix::HepSymMatrix
HepSymMatrix()
CLHEP::HepSymMatrix::sub
HepSymMatrix sub(int min_row, int max_row) const
Definition: SymMatrix.cc:134
CLHEP::HepSymMatrix::invertHaywood5
void invertHaywood5(int &ifail)
Definition: SymMatrixInvert.cc:124
CLHEP::HepSymMatrix::house_with_update2
friend void house_with_update2(HepSymMatrix *a, HepMatrix *v, int row, int col)
Definition: MatrixLinear.cc:459
CLHEP::house_with_update2
void house_with_update2(HepSymMatrix *a, HepMatrix *v, int row=1, int col=1)
Definition: MatrixLinear.cc:459
CLHEP::HepSymMatrix::HepMatrix
friend class HepMatrix
Definition: Matrix/CLHEP/Matrix/SymMatrix.h:243
CLHEP::HepSymMatrix::assign
void assign(const HepMatrix &hm2)
Definition: SymMatrix.cc:718
CLHEP::HepSymMatrix::HepSymMatrix_row_const
friend class HepSymMatrix_row_const
Definition: Matrix/CLHEP/Matrix/SymMatrix.h:242
CLHEP::HepSymMatrix::HepSymMatrix_row_const::HepSymMatrix_row_const
HepSymMatrix_row_const(const HepSymMatrix &, int)
CLHEP::tridiagonal
void tridiagonal(HepSymMatrix *a, HepMatrix *hsm)
Definition: MatrixLinear.cc:774
CLHEP::HepSymMatrix::trace
double trace() const
Definition: SymMatrix.cc:957
CLHEP::HepSymMatrix::operator=
HepSymMatrix & operator=(const HepSymMatrix &hm2)
Definition: SymMatrix.cc:645
CLHEP::HepSymMatrix::inverse
HepSymMatrix inverse() const
CLHEP::HepSymMatrix::vT_times_v
friend HepSymMatrix vT_times_v(const HepVector &v)
Definition: SymMatrix.cc:542
CLHEP::HepSymMatrix::house
friend HepVector house(const HepSymMatrix &a, int row, int col)
Definition: MatrixLinear.cc:350
CLHEP::HepSymMatrix::operator-
HepSymMatrix operator-() const
Definition: SymMatrix.cc:214
CLHEP::HepSymMatrix::operator+=
HepSymMatrix & operator+=(const HepSymMatrix &hm2)
Definition: SymMatrix.cc:578
CLHEP::diagonalize
HepMatrix diagonalize(HepSymMatrix *s)
Definition: MatrixLinear.cc:315
CLHEP::HepSymMatrix::~HepSymMatrix
virtual ~HepSymMatrix()
Definition: SymMatrix.cc:103
CLHEP::operator/
HepLorentzVector operator/(const HepLorentzVector &, double a)
Definition: LorentzVector.cc:162
f
void f(void g())
Definition: excDblThrow.cc:38
CLHEP
Definition: ClhepVersion.h:13
CLHEP::operator-
Hep3Vector operator-(const Hep3Vector &, const Hep3Vector &)
CLHEP::HepSymMatrix::invertHaywood4
void invertHaywood4(int &ifail)
Definition: SymMatrixInvert.cc:1038
CLHEP::HepSymMatrix::HepSymMatrix_row_const::operator[]
const double & operator[](int) const
v
they are gone ZOOM Features Discontinued The following features of the ZOOM package were felt to be extreme overkill These have been after checking that no existing user code was utilizing as in SpaceVector v
Definition: keyMergeIssues.doc:324
CLHEP::HepSymMatrix::HepSymMatrix_row::HepSymMatrix_row
HepSymMatrix_row(HepSymMatrix &, int)
CLHEP::HepSymMatrix::HepSymMatrix_row::operator[]
double & operator[](int)
CLHEP::HepSymMatrix::operator-=
HepSymMatrix & operator-=(const HepSymMatrix &hm2)
Definition: SymMatrix.cc:601
CLHEP::HepSymMatrix::invertCholesky5
void invertCholesky5(int &ifail)
Definition: SymMatrixInvert.cc:683
CLHEP::HepSymMatrix::operator()
const double & operator()(int row, int col) const
CLHEP::HepSymMatrix::num_size
int num_size() const
CLHEP::HepSymMatrix::HepSymMatrix_row
friend class HepSymMatrix_row
Definition: Matrix/CLHEP/Matrix/SymMatrix.h:241
CLHEP::HepSymMatrix::num_row
int num_row() const
CLHEP::HepSymMatrix::determinant
double determinant() const
Definition: SymMatrix.cc:943
CLHEP::HepSymMatrix::fast
const double & fast(int row, int col) const
CLHEP::HepSymMatrix::condition
friend double condition(const HepSymMatrix &m)
Definition: MatrixLinear.cc:198
CLHEP::operator+
Hep3Vector operator+(const Hep3Vector &, const Hep3Vector &)
CLHEP::HepSymMatrix::similarityT
HepSymMatrix similarityT(const HepMatrix &hm1) const
Definition: SymMatrix.cc:816
CLHEP::HepSymMatrix::HepDiagMatrix
friend class HepDiagMatrix
Definition: Matrix/CLHEP/Matrix/SymMatrix.h:244
CLHEP::HepSymMatrix::apply
HepSymMatrix apply(double(*f)(double, int, int)) const
Definition: SymMatrix.cc:700
CLHEP::HepSymMatrix::T
HepSymMatrix T() const
CLHEP::HepSymMatrix::similarity
HepSymMatrix similarity(const HepMatrix &hm1) const
Definition: SymMatrix.cc:737
CLHEP::HepSymMatrix::operator*
friend HepMatrix operator*(const HepSymMatrix &hm1, const HepSymMatrix &hm2)
Definition: SymMatrix.cc:437
CLHEP::operator<<
std::ostream & operator<<(std::ostream &os, const HepAxisAngle &aa)
Definition: AxisAngle.cc:86
CLHEP::HepSymMatrix::invertBunchKaufman
void invertBunchKaufman(int &ifail)
Definition: SymMatrix.cc:964
CLHEP::HepSymMatrix::invertHaywood6
void invertHaywood6(int &ifail)
Definition: SymMatrixInvert.cc:295
CLHEP::operator*
HepLorentzRotation operator*(const HepRotation &r, const HepLorentzRotation &lt)
Definition: LorentzRotation.cc:262
n_constructors::m
int m
Definition: testSharedPtr.cc:370
CLHEP::HepSymMatrix::tridiagonal
friend void tridiagonal(HepSymMatrix *a, HepMatrix *hsm)
Definition: MatrixLinear.cc:774
CLHEP::dsum
HepDiagMatrix dsum(const HepDiagMatrix &s1, const HepDiagMatrix &s2)
Definition: DiagMatrix.cc:164
CLHEP::diag_step
void diag_step(HepSymMatrix *t, int begin, int end)
Definition: MatrixLinear.cc:224
CLHEP::HepSymMatrix::diagonalize
friend HepMatrix diagonalize(HepSymMatrix *s)
Definition: MatrixLinear.cc:315
CLHEP::HepSymMatrix::num_col
int num_col() const