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

Matrix/CLHEP/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
77namespace CLHEP {
78
79class HepRandom;
80
81class HepMatrix;
82class HepDiagMatrix;
83class HepVector;
84
89class HepSymMatrix : public HepGenMatrix {
90public:
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);
141 HepSymMatrix & operator-=( const HepSymMatrix &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
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:
202 inline double & operator[](int);
203 private:
204 HepSymMatrix& _a;
205 int _r;
206 };
208 public:
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
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
237protected:
238 inline int num_size() const;
239
240private:
241 friend class HepSymMatrix_row;
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
296std::ostream& operator<<(std::ostream &s, const HepSymMatrix &q);
297// Write out Matrix, SymMatrix, DiagMatrix and Vector into ostream.
298
299HepMatrix operator*(const HepMatrix &hm1, const HepSymMatrix &hm2);
300HepMatrix operator*(const HepSymMatrix &hm1, const HepMatrix &hm2);
301HepMatrix operator*(const HepSymMatrix &hm1, const HepSymMatrix &hm2);
302HepSymMatrix operator*(double t, const HepSymMatrix &s1);
303HepSymMatrix operator*(const HepSymMatrix &s1, double t);
304// Multiplication operators.
305// Note that m *= hm1 is always faster than m = m * hm1
306
307HepSymMatrix operator/(const HepSymMatrix &hm1, double t);
308// s = s1 / t. (s /= t is faster if you can use it.)
309
310HepMatrix operator+(const HepMatrix &hm1, const HepSymMatrix &s2);
311HepMatrix operator+(const HepSymMatrix &s1, const HepMatrix &hm2);
312HepSymMatrix operator+(const HepSymMatrix &s1, const HepSymMatrix &s2);
313// Addition operators
314
315HepMatrix operator-(const HepMatrix &hm1, const HepSymMatrix &s2);
316HepMatrix operator-(const HepSymMatrix &hm1, const HepMatrix &hm2);
317HepSymMatrix operator-(const HepSymMatrix &s1, const HepSymMatrix &s2);
318// subtraction operators
319
320HepSymMatrix dsum(const HepSymMatrix &s1, const HepSymMatrix &s2);
321// Direct sum of two symmetric matrices;
322
323double condition(const HepSymMatrix &m);
324// Find the conditon number of a symmetric matrix.
325
326void diag_step(HepSymMatrix *t, int begin, int end);
327void diag_step(HepSymMatrix *t, HepMatrix *u, int begin, int end);
328// Implicit symmetric QR step with Wilkinson Shift
329
331// Diagonalize a symmetric matrix.
332// It returns the matrix U so that s_old = U * s_diag * U.T()
333
334HepVector house(const HepSymMatrix &a, int row=1, int col=1);
335void house_with_update2(HepSymMatrix *a, HepMatrix *v, int row=1, int col=1);
336// Finds and does Householder reflection on matrix.
337
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
346using namespace CLHEP;
347#endif
348
349#ifndef HEP_DEBUG_INLINE
350#include "CLHEP/Matrix/SymMatrix.icc"
351#endif
352
353#endif
HepSymMatrix_row_const(const HepSymMatrix &, int)
HepSymMatrix_row(HepSymMatrix &, int)
HepSymMatrix & operator*=(double t)
Definition SymMatrix.cc:614
double trace() const
Definition SymMatrix.cc:957
friend double condition(const HepSymMatrix &m)
double & fast(int row, int col)
void assign(const HepSymMatrix &hm2)
friend HepVector house(const HepSymMatrix &a, int row, int col)
void invertHaywood5(int &ifail)
double determinant() const
Definition SymMatrix.cc:943
int num_row() const
void invertCholesky6(int &ifail)
friend void house_with_update2(HepSymMatrix *a, HepMatrix *v, int row, int col)
void invertHaywood4(int &ifail)
HepSymMatrix_row operator[](int)
HepSymMatrix & operator+=(const HepSymMatrix &hm2)
Definition SymMatrix.cc:578
int num_col() const
void assign(const HepMatrix &hm2)
Definition SymMatrix.cc:718
HepSymMatrix apply(double(*f)(double, int, int)) const
Definition SymMatrix.cc:700
friend void tridiagonal(HepSymMatrix *a, HepMatrix *hsm)
int num_size() const
friend HepMatrix diagonalize(HepSymMatrix *s)
HepSymMatrix T() const
friend void diag_step(HepSymMatrix *t, int begin, int end)
HepSymMatrix similarityT(const HepMatrix &hm1) const
Definition SymMatrix.cc:816
HepSymMatrix & operator=(const HepSymMatrix &hm2)
Definition SymMatrix.cc:645
const double & operator()(int row, int col) const
void invertCholesky5(int &ifail)
HepSymMatrix sub(int min_row, int max_row) const
Definition SymMatrix.cc:134
friend HepMatrix operator*(const HepSymMatrix &hm1, const HepSymMatrix &hm2)
Definition SymMatrix.cc:437
HepSymMatrix operator-() const
Definition SymMatrix.cc:214
HepSymMatrix & operator-=(const HepSymMatrix &hm2)
Definition SymMatrix.cc:601
void invertHaywood6(int &ifail)
HepSymMatrix inverse(int &ifail) const
HepSymMatrix similarity(const HepMatrix &hm1) const
Definition SymMatrix.cc:737
friend HepSymMatrix vT_times_v(const HepVector &v)
Definition SymMatrix.cc:542
void invertBunchKaufman(int &ifail)
Definition SymMatrix.cc:964
const double & fast(int row, int col) const
double & operator()(int row, int col)
virtual ~HepSymMatrix()
Definition SymMatrix.cc:103
friend HepSymMatrix operator+(const HepSymMatrix &hm1, const HepSymMatrix &hm2)
Definition SymMatrix.cc:256
HepSymMatrix inverse() const
HepSymMatrix & operator/=(double t)
Definition SymMatrix.cc:608
void f(void g())
void tridiagonal(HepSymMatrix *a, HepMatrix *hsm)
void house_with_update2(HepSymMatrix *a, HepMatrix *v, int row=1, int col=1)
Hep3Vector operator-(const Hep3Vector &, const Hep3Vector &)
Hep3Vector operator+(const Hep3Vector &, const Hep3Vector &)
HepLorentzRotation operator*(const HepRotation &r, const HepLorentzRotation &lt)
HepVector house(const HepMatrix &a, int row=1, int col=1)
std::ostream & operator<<(std::ostream &os, const HepAxisAngle &aa)
Definition AxisAngle.cc:86
HepDiagMatrix dsum(const HepDiagMatrix &s1, const HepDiagMatrix &s2)
double condition(const HepSymMatrix &m)
HepMatrix diagonalize(HepSymMatrix *s)
HepLorentzVector operator/(const HepLorentzVector &, double a)
void diag_step(HepSymMatrix *t, int begin, int end)