CLHEP VERSION Reference Documentation
CLHEP Home Page
CLHEP Documentation
CLHEP Bug Reports
Main Page
Namespaces
Classes
Files
File List
File Members
RandomObjects
CLHEP
Matrix
RandomObjects/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
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
&m1);
104
// Copy constructor.
105
106
HepSymMatrix
(
const
HepDiagMatrix
&m1);
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
&m2);
128
// Assigns m2 to s, assuming m2 is a symmetric matrix.
129
130
void
assign
(
const
HepSymMatrix
&m2);
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
&m2);
140
HepSymMatrix
&
operator+=
(
const
HepDiagMatrix
&m2);
141
HepSymMatrix
&
operator-=
(
const
HepSymMatrix
&m2);
142
HepSymMatrix
&
operator-=
(
const
HepDiagMatrix
&m2);
143
// Add or subtract a SymMatrix.
144
145
HepSymMatrix
&
operator=
(
const
HepSymMatrix
&m2);
146
HepSymMatrix
&
operator=
(
const
HepDiagMatrix
&m2);
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
&m1)
const
;
159
HepSymMatrix
similarity
(
const
HepSymMatrix
&m1)
const
;
160
// Returns m1*s*m1.T().
161
162
HepSymMatrix
similarityT
(
const
HepMatrix
&m1)
const
;
163
// temporary. test of new similarity.
164
// Returns m1.T()*s*m1.
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
&m1);
172
// Sub matrix of this SymMatrix is replaced with m1.
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
double
determinant
()
const
;
188
// calculate the determinant of the matrix.
189
190
double
trace
()
const
;
191
// calculate the trace of the matrix (sum of diagonal elements).
192
193
class
HepSymMatrix_row
{
194
public
:
195
inline
HepSymMatrix_row
(
HepSymMatrix
&,
int
);
196
inline
double
&
operator[]
(
int
);
197
private
:
198
HepSymMatrix
& _a;
199
int
_r;
200
};
201
class
HepSymMatrix_row_const
{
202
public
:
203
inline
HepSymMatrix_row_const
(
const
HepSymMatrix
&,
int
);
204
inline
const
double
&
operator[]
(
int
)
const
;
205
private
:
206
const
HepSymMatrix
& _a;
207
int
_r;
208
};
209
// helper class to implement m[i][j]
210
211
inline
HepSymMatrix_row
operator[]
(
int
);
212
inline
HepSymMatrix_row_const
operator[]
(
int
)
const
;
213
// Read or write a matrix element.
214
// While it may not look like it, you simply do m[i][j] to get an
215
// element.
216
// ** Note that the indexing starts from [0][0]. **
217
218
// Special-case inversions for 5x5 and 6x6 symmetric positive definite:
219
// These set ifail=0 and invert if the matrix was positive definite;
220
// otherwise ifail=1 and the matrix is left unaltered.
221
void
invertCholesky5
(
int
&ifail);
222
void
invertCholesky6
(
int
&ifail);
223
224
// Inversions for 5x5 and 6x6 forcing use of specific methods: The
225
// behavior (though not the speed) will be identical to invert(ifail).
226
void
invertHaywood4
(
int
& ifail);
227
void
invertHaywood5
(
int
&ifail);
228
void
invertHaywood6
(
int
&ifail);
229
void
invertBunchKaufman
(
int
&ifail);
230
231
protected
:
232
inline
int
num_size
()
const
;
233
234
private
:
235
friend
class
HepSymMatrix_row
;
236
friend
class
HepSymMatrix_row_const
;
237
friend
class
HepMatrix
;
238
friend
class
HepDiagMatrix
;
239
240
friend
void
tridiagonal
(
HepSymMatrix
*
a
,
HepMatrix
*hsm);
241
friend
double
condition
(
const
HepSymMatrix
&m);
242
friend
void
diag_step
(
HepSymMatrix
*t,
int
begin,
int
end);
243
friend
void
diag_step
(
HepSymMatrix
*t,
HepMatrix
*u,
int
begin,
int
end);
244
friend
HepMatrix
diagonalize
(
HepSymMatrix
*s);
245
friend
HepVector
house
(
const
HepSymMatrix
&
a
,
int
row,
int
col);
246
friend
void
house_with_update2
(
HepSymMatrix
*
a
,
HepMatrix
*v,
int
row,
int
col);
247
248
friend
HepSymMatrix
operator+
(
const
HepSymMatrix
&m1,
249
const
HepSymMatrix
&m2);
250
friend
HepSymMatrix
operator-
(
const
HepSymMatrix
&m1,
251
const
HepSymMatrix
&m2);
252
friend
HepMatrix
operator*
(
const
HepSymMatrix
&m1,
const
HepSymMatrix
&m2);
253
friend
HepMatrix
operator*
(
const
HepSymMatrix
&m1,
const
HepMatrix
&m2);
254
friend
HepMatrix
operator*
(
const
HepMatrix
&m1,
const
HepSymMatrix
&m2);
255
friend
HepVector
operator*
(
const
HepSymMatrix
&m1,
const
HepVector &m2);
256
// Multiply a Matrix by a Matrix or Vector.
257
258
friend
HepSymMatrix
vT_times_v
(
const
HepVector &v);
259
// Returns v * v.T();
260
261
#ifdef DISABLE_ALLOC
262
std::vector<double > m;
263
#else
264
std::vector<double,Alloc<double,25> > m;
265
#endif
266
int
nrow;
267
int
size_;
// total number of elements
268
269
static
double
posDefFraction5x5;
270
static
double
adjustment5x5;
271
static
const
double
CHOLESKY_THRESHOLD_5x5;
272
static
const
double
CHOLESKY_CREEP_5x5;
273
274
static
double
posDefFraction6x6;
275
static
double
adjustment6x6;
276
static
const
double
CHOLESKY_THRESHOLD_6x6;
277
static
const
double
CHOLESKY_CREEP_6x6;
278
279
void
invert4 (
int
& ifail);
280
void
invert5 (
int
& ifail);
281
void
invert6 (
int
& ifail);
282
283
};
284
285
//
286
// Operations other than member functions for Matrix, SymMatrix, DiagMatrix
287
// and Vectors implemented in Matrix.cc and Matrix.icc (inline).
288
//
289
290
std::ostream&
operator<<
(std::ostream &s,
const
HepSymMatrix &q);
291
// Write out Matrix, SymMatrix, DiagMatrix and Vector into ostream.
292
293
HepMatrix
operator*
(
const
HepMatrix &m1,
const
HepSymMatrix &m2);
294
HepMatrix
operator*
(
const
HepSymMatrix &m1,
const
HepMatrix &m2);
295
HepMatrix
operator*
(
const
HepSymMatrix &m1,
const
HepSymMatrix &m2);
296
HepSymMatrix
operator*
(
double
t,
const
HepSymMatrix &s1);
297
HepSymMatrix
operator*
(
const
HepSymMatrix &s1,
double
t);
298
// Multiplication operators.
299
// Note that m *= m1 is always faster than m = m * m1
300
301
HepSymMatrix
operator/
(
const
HepSymMatrix &m1,
double
t);
302
// s = s1 / t. (s /= t is faster if you can use it.)
303
304
HepMatrix
operator+
(
const
HepMatrix &m1,
const
HepSymMatrix &s2);
305
HepMatrix
operator+
(
const
HepSymMatrix &s1,
const
HepMatrix &m2);
306
HepSymMatrix
operator+
(
const
HepSymMatrix &s1,
const
HepSymMatrix &s2);
307
// Addition operators
308
309
HepMatrix
operator-
(
const
HepMatrix &m1,
const
HepSymMatrix &s2);
310
HepMatrix
operator-
(
const
HepSymMatrix &m1,
const
HepMatrix &m2);
311
HepSymMatrix
operator-
(
const
HepSymMatrix &s1,
const
HepSymMatrix &s2);
312
// subtraction operators
313
314
HepSymMatrix
dsum
(
const
HepSymMatrix &s1,
const
HepSymMatrix &s2);
315
// Direct sum of two symmetric matrices;
316
317
double
condition
(
const
HepSymMatrix &m);
318
// Find the conditon number of a symmetric matrix.
319
320
void
diag_step
(HepSymMatrix *t,
int
begin,
int
end);
321
void
diag_step
(HepSymMatrix *t, HepMatrix *u,
int
begin,
int
end);
322
// Implicit symmetric QR step with Wilkinson Shift
323
324
HepMatrix
diagonalize
(HepSymMatrix *s);
325
// Diagonalize a symmetric matrix.
326
// It returns the matrix U so that s_old = U * s_diag * U.T()
327
328
HepVector
house
(
const
HepSymMatrix &
a
,
int
row=1,
int
col=1);
329
void
house_with_update2
(HepSymMatrix *
a
, HepMatrix *v,
int
row=1,
int
col=1);
330
// Finds and does Householder reflection on matrix.
331
332
void
tridiagonal
(HepSymMatrix *
a
, HepMatrix *hsm);
333
HepMatrix
tridiagonal
(HepSymMatrix *
a
);
334
// Does a Householder tridiagonalization of a symmetric matrix.
335
336
}
// namespace CLHEP
337
338
#ifdef ENABLE_BACKWARDS_COMPATIBILITY
339
// backwards compatibility will be enabled ONLY in CLHEP 1.9
340
using namespace
CLHEP;
341
#endif
342
343
#ifndef HEP_DEBUG_INLINE
344
#include "CLHEP/Matrix/SymMatrix.icc"
345
#endif
346
347
#endif
Generated on Sun Jun 17 2012 08:08:27 for CLHEP by
1.8.1.1