CLHEP VERSION Reference Documentation
CLHEP Home Page
CLHEP Documentation
CLHEP Bug Reports
Main Page
Namespaces
Classes
Files
File List
File Members
Matrix
Matrix
Matrix/Matrix/Matrix.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 HepMatrix class.
9
//
10
// This software written by Nobu Katayama and Mike Smyth, Cornell University.
11
//
12
//
13
// .SS Usage
14
// The Matrix class does all the obvious things, in all the obvious ways.
15
// You declare a Matrix by saying
16
//
17
// .ft B
18
// HepMatrix m1(n, m);
19
//
20
// To declare a Matrix as a copy of a Matrix m2, say
21
//
22
// .ft B
23
// HepMatrix m1(m2);
24
// or
25
// .ft B
26
// HepMatrix m1 = m2;
27
//
28
// You can declare initilizations of a Matrix, by giving it a third
29
// integer argument, either 0 or 1. 0 means initialize to 0, one means
30
// the identity matrix. If you do not give a third argument the memory
31
// is not initialized.
32
//
33
// .ft B
34
// HepMatrix m1(n, m, 1);
35
//
36
// ./"This code has been written by Mike Smyth, and the algorithms used are
37
// ./"described in the thesis "A Tracking Library for a Silicon Vertex Detector"
38
// ./"(Mike Smyth, Cornell University, June 1993).
39
// ./"This is file contains C++ stuff for doing things with Matrices.
40
// ./"To turn on bound checking, define MATRIX_BOUND_CHECK before including
41
// ./"this file.
42
//
43
// To find the number of rows or number of columns, say
44
//
45
// .ft B
46
// nr = m1.num_row();
47
//
48
// or
49
//
50
// .ft B
51
// nc = m1.num_col();
52
//
53
// You can print a Matrix by
54
//
55
// .ft B
56
// cout << m1;
57
//
58
// You can add,
59
// subtract, and multiply two Matrices. You can multiply a Matrix by a
60
// scalar, on the left or the right. +=, *=, -= act in the obvious way.
61
// m1 *= m2 is as m1 = m1*m2. You can also divide by a scalar on the
62
// right, or use /= to do the same thing.
63
//
64
// You can read or write a Matrix element by saying
65
//
66
// .ft B
67
// m(row, col) = blah. (1 <= row <= num_row(), 1 <= col <= num_col())
68
//
69
// .ft B
70
// blah = m(row, col) ...
71
//
72
// m(row, col) is inline, and by default does not do bound checking.
73
// If bound checking is desired, say #define MATRIX_BOUND_CHECK before
74
// including matrix.h.
75
//
76
// You can also access the element using C subscripting operator []
77
//
78
// .ft B
79
// m[row][col] = blah. (0 <= row < num_row(), 0 <= col < num_col())
80
//
81
// .ft B
82
// blah = m[row][col] ...
83
//
84
// m[row][col] is inline, and by default does not do bound checking.
85
// If bound checking is desired, say #define MATRIX_BOUND_CHECK before
86
// including matrix.h. (Notice the difference in bases in two access
87
// methods.)
88
//
89
// .SS Comments on precision
90
//
91
// The user would normally use "Matrix" class whose precision is the same
92
// as the other classes of CLHEP (ThreeVec, for example). However, he/she
93
// can explicitly choose Matrix (double) or MatrixF (float) if he/she wishes.
94
// In the former case, include "Matrix.h." In the latter case, include either
95
// "Matrix.h," or "MatrixF.h," or both. The only operators that can talk
96
// to each other are the copy constructor and assignment operator.
97
//
98
// .ft B
99
// Matrix d(3,4,HEP_MATRIX_IDENTITY);
100
//
101
// .ft B
102
// MatrixF f;
103
//
104
// .ft B
105
// f = d;
106
//
107
// .ft B
108
// MatrixF g(d);
109
//
110
// will convert from one to the other.
111
//
112
// .SS Other functions
113
//
114
// .ft B
115
// mt = m.T();
116
//
117
// returns the transpose of m.
118
//
119
// .ft B
120
// ms = m2.sub(row_min, row_max, col_min, col_max);
121
//
122
// returns the subMatrix.
123
// m2(row_min:row_max, col_min:col_max) in matlab terminology.
124
// If instead you say
125
//
126
// .ft B
127
// m2.sub(row, col, m1);
128
//
129
// then the subMatrix
130
// m2(row:row+m1.num_row()-1, col:col+m1.num_col()-1) is replaced with m1.
131
//
132
// .ft B
133
// md = dsum(m1,m2);
134
//
135
// returns the direct sum of the two matrices.
136
//
137
// Operations that do not have to be member functions or friends are listed
138
// towards the end of this man page.
139
//
140
//
141
// To invert a matrix, say
142
//
143
// .ft B
144
// minv = m.inverse(ierr);
145
//
146
// ierr will be non-zero if the matrix is singular.
147
//
148
// If you can overwrite the matrix, you can use the invert function to avoid
149
// two extra copies. Use
150
//
151
// .ft B
152
// m.invert(ierr);
153
//
154
// Many basic linear algebra functions such as QR decomposition are defined.
155
// In order to keep the header file a reasonable size, the functions are
156
// defined in MatrixLinear.h.
157
//
158
//
159
// .ft B
160
// Note that Matrices run from (1,1) to (n, m), and [0][0] to
161
// [n-1][m-1]. The users of the latter form should be careful about sub
162
// functions.
163
//
164
// ./" The program that this class was orginally used with lots of small
165
// ./" Matrices. It was very costly to malloc and free array space every
166
// ./" time a Matrix is needed. So, a number of previously freed arrays are
167
// ./" kept around, and when needed again one of these array is used. Right
168
// ./" now, a set of piles of arrays with row <= row_max and col <= col_max
169
// ./" are kept around. The piles of kept Matrices can be easily changed.
170
// ./" Array mallocing and freeing are done only in new_m() and delete_m(),
171
// ./" so these are the only routines that need to be rewritten.
172
//
173
// You can do things thinking of a Matrix as a list of numbers.
174
//
175
// .ft B
176
// m = m1.apply(HEP_MATRIX_ELEMENT (*f)(HEP_MATRIX_ELEMENT, int r, int c));
177
//
178
// applies f to every element of m1 and places it in m.
179
//
180
// .SS See Also:
181
// SymMatrix[DF].h, GenMatrix[DF].h, DiagMatrix[DF].h Vector[DF].h
182
// MatrixLinear[DF].h
183
184
#ifndef _Matrix_H_
185
#define _Matrix_H_
186
187
#ifdef GNUPRAGMA
188
#pragma interface
189
#endif
190
191
#include <vector>
192
193
#include "CLHEP/Matrix/defs.h"
194
#include "CLHEP/Matrix/GenMatrix.h"
195
196
namespace
CLHEP {
197
198
class
HepRandom;
199
200
class
HepSymMatrix;
201
class
HepDiagMatrix;
202
class
HepVector;
203
class
HepRotation;
204
209
class
HepMatrix :
public
HepGenMatrix {
210
public
:
211
inline
HepMatrix
();
212
// Default constructor. Gives 0 x 0 matrix. Another Matrix can be
213
// assigned to it.
214
215
HepMatrix
(
int
p,
int
q);
216
// Constructor. Gives an unitialized p x q matrix.
217
HepMatrix
(
int
p,
int
q,
int
i);
218
// Constructor. Gives an initialized p x q matrix.
219
// If i=0, it is initialized to all 0. If i=1, the diagonal elements
220
// are set to 1.0.
221
222
HepMatrix
(
int
p,
int
q, HepRandom &r);
223
// Constructor with a Random object.
224
225
HepMatrix
(
const
HepMatrix
&m1);
226
// Copy constructor.
227
228
HepMatrix
(
const
HepSymMatrix
&m1);
229
HepMatrix
(
const
HepDiagMatrix
&m1);
230
HepMatrix
(
const
HepVector
&m1);
231
// Constructors from SymMatrix, DiagMatrix and Vector.
232
233
virtual
~HepMatrix
();
234
// Destructor.
235
236
virtual
int
num_row
()
const
;
237
// Returns the number of rows.
238
239
virtual
int
num_col
()
const
;
240
// Returns the number of columns.
241
242
virtual
const
double
&
operator()
(
int
row,
int
col)
const
;
243
virtual
double
&
operator()
(
int
row,
int
col);
244
// Read or write a matrix element.
245
// ** Note that the indexing starts from (1,1). **
246
247
HepMatrix
&
operator *=
(
double
t);
248
// Multiply a Matrix by a floating number.
249
250
HepMatrix
&
operator /=
(
double
t);
251
// Divide a Matrix by a floating number.
252
253
HepMatrix
&
operator +=
(
const
HepMatrix
&m2);
254
HepMatrix
&
operator +=
(
const
HepSymMatrix
&m2);
255
HepMatrix
&
operator +=
(
const
HepDiagMatrix
&m2);
256
HepMatrix
&
operator +=
(
const
HepVector
&m2);
257
HepMatrix
&
operator -=
(
const
HepMatrix
&m2);
258
HepMatrix
&
operator -=
(
const
HepSymMatrix
&m2);
259
HepMatrix
&
operator -=
(
const
HepDiagMatrix
&m2);
260
HepMatrix
&
operator -=
(
const
HepVector
&m2);
261
// Add or subtract a Matrix.
262
// When adding/subtracting Vector, Matrix must have num_col of one.
263
264
HepMatrix
&
operator =
(
const
HepMatrix
&m2);
265
HepMatrix
&
operator =
(
const
HepSymMatrix
&m2);
266
HepMatrix
&
operator =
(
const
HepDiagMatrix
&m2);
267
HepMatrix
&
operator =
(
const
HepVector
&m2);
268
HepMatrix
&
operator =
(
const
HepRotation &m2);
269
// Assignment operators.
270
271
HepMatrix
operator-
()
const
;
272
// unary minus, ie. flip the sign of each element.
273
274
HepMatrix
apply
(
double
(*
f
)(
double
,
int
,
int
))
const
;
275
// Apply a function to all elements of the matrix.
276
277
HepMatrix
T
()
const
;
278
// Returns the transpose of a Matrix.
279
280
HepMatrix
sub
(
int
min_row,
int
max_row,
int
min_col,
int
max_col)
const
;
281
// Returns a sub matrix of a Matrix.
282
// WARNING: rows and columns are numbered from 1
283
void
sub
(
int
row,
int
col,
const
HepMatrix
&m1);
284
// Sub matrix of this Matrix is replaced with m1.
285
// WARNING: rows and columns are numbered from 1
286
287
friend
inline
void
swap
(
HepMatrix
&m1,
HepMatrix
&m2);
288
// Swap m1 with m2.
289
290
inline
HepMatrix
inverse
(
int
& ierr)
const
;
291
// Invert a Matrix. Matrix must be square and is not changed.
292
// Returns ierr = 0 (zero) when successful, otherwise non-zero.
293
294
virtual
void
invert
(
int
& ierr);
295
// Invert a Matrix. Matrix must be square.
296
// N.B. the contents of the matrix are replaced by the inverse.
297
// Returns ierr = 0 (zero) when successful, otherwise non-zero.
298
// This method has less overhead then inverse().
299
300
double
determinant
()
const
;
301
// calculate the determinant of the matrix.
302
303
double
trace
()
const
;
304
// calculate the trace of the matrix (sum of diagonal elements).
305
306
class
HepMatrix_row
{
307
public
:
308
inline
HepMatrix_row
(
HepMatrix
&,
int
);
309
double
&
operator[]
(
int
);
310
private
:
311
HepMatrix
& _a;
312
int
_r;
313
};
314
class
HepMatrix_row_const
{
315
public
:
316
inline
HepMatrix_row_const
(
const
HepMatrix
&,
int
);
317
const
double
&
operator[]
(
int
)
const
;
318
private
:
319
const
HepMatrix
& _a;
320
int
_r;
321
};
322
// helper classes for implementing m[i][j]
323
324
inline
HepMatrix_row
operator[]
(
int
);
325
inline
const
HepMatrix_row_const
operator[]
(
int
)
const
;
326
// Read or write a matrix element.
327
// While it may not look like it, you simply do m[i][j] to get an
328
// element.
329
// ** Note that the indexing starts from [0][0]. **
330
331
protected
:
332
virtual
int
num_size
()
const
;
333
virtual
void
invertHaywood4
(
int
& ierr);
334
virtual
void
invertHaywood5
(
int
& ierr);
335
virtual
void
invertHaywood6
(
int
& ierr);
336
337
private
:
338
friend
class
HepMatrix_row
;
339
friend
class
HepMatrix_row_const
;
340
friend
class
HepVector
;
341
friend
class
HepSymMatrix
;
342
friend
class
HepDiagMatrix
;
343
// Friend classes.
344
345
friend
HepMatrix
operator+
(
const
HepMatrix
&m1,
const
HepMatrix
&m2);
346
friend
HepMatrix
operator-
(
const
HepMatrix
&m1,
const
HepMatrix
&m2);
347
friend
HepMatrix
operator*
(
const
HepMatrix
&m1,
const
HepMatrix
&m2);
348
friend
HepMatrix
operator*
(
const
HepMatrix
&m1,
const
HepSymMatrix
&m2);
349
friend
HepMatrix
operator*
(
const
HepMatrix
&m1,
const
HepDiagMatrix
&m2);
350
friend
HepMatrix
operator*
(
const
HepSymMatrix
&m1,
const
HepMatrix
&m2);
351
friend
HepMatrix
operator*
(
const
HepDiagMatrix
&m1,
const
HepMatrix
&m2);
352
friend
HepMatrix
operator*
(
const
HepVector
&m1,
const
HepMatrix
&m2);
353
friend
HepVector
operator*
(
const
HepMatrix
&m1,
const
HepVector
&m2);
354
friend
HepMatrix
operator*
(
const
HepSymMatrix
&m1,
const
HepSymMatrix
&m2);
355
// Multiply a Matrix by a Matrix or Vector.
356
357
friend
HepVector
solve
(
const
HepMatrix
&,
const
HepVector
&);
358
// solve the system of linear eq
359
friend
HepVector
qr_solve
(
HepMatrix
*,
const
HepVector
&);
360
friend
HepMatrix
qr_solve
(
HepMatrix
*,
const
HepMatrix
&
b
);
361
friend
void
tridiagonal
(
HepSymMatrix
*
a
,
HepMatrix
*hsm);
362
friend
void
row_house
(
HepMatrix
*,
const
HepMatrix
&,
double
,
363
int
,
int
,
int
,
int
);
364
friend
void
row_house
(
HepMatrix
*,
const
HepVector
&,
double
,
365
int
,
int
);
366
friend
void
back_solve
(
const
HepMatrix
&R,
HepVector
*
b
);
367
friend
void
back_solve
(
const
HepMatrix
&R,
HepMatrix
*
b
);
368
friend
void
col_givens
(
HepMatrix
*A,
double
c,
369
double
s,
int
k1,
int
k2,
370
int
rowmin,
int
rowmax);
371
// Does a column Givens update.
372
friend
void
row_givens
(
HepMatrix
*A,
double
c,
373
double
s,
int
k1,
int
k2,
374
int
colmin,
int
colmax);
375
friend
void
col_house
(
HepMatrix
*,
const
HepMatrix
&,
double
,
376
int
,
int
,
int
,
int
);
377
friend
HepVector
house
(
const
HepMatrix
&
a
,
int
row,
int
col);
378
friend
void
house_with_update
(
HepMatrix
*
a
,
int
row,
int
col);
379
friend
void
house_with_update
(
HepMatrix
*
a
,
HepMatrix
*v,
int
row,
int
col);
380
friend
void
house_with_update2
(
HepSymMatrix
*
a
,
HepMatrix
*v,
381
int
row,
int
col);
382
383
int
dfact_matrix(
double
&det,
int
*ir);
384
// factorize the matrix. If successful, the return code is 0. On
385
// return, det is the determinant and ir[] is row-interchange
386
// matrix. See CERNLIB's DFACT routine.
387
388
int
dfinv_matrix(
int
*ir);
389
// invert the matrix. See CERNLIB DFINV.
390
391
#ifdef DISABLE_ALLOC
392
std::vector<double > m;
393
#else
394
std::vector<double,Alloc<double,25> > m;
395
#endif
396
int
nrow, ncol;
397
int
size_;
398
};
399
400
// Operations other than member functions for Matrix
401
// implemented in Matrix.cc and Matrix.icc (inline).
402
403
HepMatrix
operator*
(
const
HepMatrix &m1,
const
HepMatrix &m2);
404
HepMatrix
operator*
(
double
t,
const
HepMatrix &m1);
405
HepMatrix
operator*
(
const
HepMatrix &m1,
double
t);
406
// Multiplication operators
407
// Note that m *= m1 is always faster than m = m * m1.
408
409
HepMatrix
operator/
(
const
HepMatrix &m1,
double
t);
410
// m = m1 / t. (m /= t is faster if you can use it.)
411
412
HepMatrix
operator+
(
const
HepMatrix &m1,
const
HepMatrix &m2);
413
// m = m1 + m2;
414
// Note that m += m1 is always faster than m = m + m1.
415
416
HepMatrix
operator-
(
const
HepMatrix &m1,
const
HepMatrix &m2);
417
// m = m1 - m2;
418
// Note that m -= m1 is always faster than m = m - m1.
419
420
HepMatrix
dsum
(
const
HepMatrix&,
const
HepMatrix&);
421
// Direct sum of two matrices. The direct sum of A and B is the matrix
422
// A 0
423
// 0 B
424
425
HepVector
solve
(
const
HepMatrix &,
const
HepVector &);
426
// solve the system of linear equations using LU decomposition.
427
428
std::ostream&
operator<<
(std::ostream &s,
const
HepMatrix &q);
429
// Read in, write out Matrix into a stream.
430
431
//
432
// Specialized linear algebra functions
433
//
434
435
HepVector
qr_solve
(
const
HepMatrix &A,
const
HepVector &
b
);
436
HepVector
qr_solve
(HepMatrix *A,
const
HepVector &
b
);
437
HepMatrix
qr_solve
(
const
HepMatrix &A,
const
HepMatrix &
b
);
438
HepMatrix
qr_solve
(HepMatrix *A,
const
HepMatrix &
b
);
439
// Works like backsolve, except matrix does not need to be upper
440
// triangular. For nonsquare matrix, it solves in the least square sense.
441
442
HepMatrix
qr_inverse
(
const
HepMatrix &A);
443
HepMatrix
qr_inverse
(HepMatrix *A);
444
// Finds the inverse of a matrix using QR decomposition. Note, often what
445
// you really want is solve or backsolve, they can be much quicker than
446
// inverse in many calculations.
447
448
449
void
qr_decomp
(HepMatrix *A, HepMatrix *hsm);
450
HepMatrix
qr_decomp
(HepMatrix *A);
451
// Does a QR decomposition of a matrix.
452
453
void
back_solve
(
const
HepMatrix &R, HepVector *
b
);
454
void
back_solve
(
const
HepMatrix &R, HepMatrix *
b
);
455
// Solves R*x = b where R is upper triangular. Also has a variation that
456
// solves a number of equations of this form in one step, where b is a matrix
457
// with each column a different vector. See also solve.
458
459
void
col_house
(HepMatrix *
a
,
const
HepMatrix &v,
double
vnormsq,
460
int
row,
int
col,
int
row_start,
int
col_start);
461
void
col_house
(HepMatrix *
a
,
const
HepMatrix &v,
int
row,
int
col,
462
int
row_start,
int
col_start);
463
// Does a column Householder update.
464
465
void
col_givens
(HepMatrix *A,
double
c,
double
s,
466
int
k1,
int
k2,
int
row_min=1,
int
row_max=0);
467
// do a column Givens update
468
469
void
row_givens
(HepMatrix *A,
double
c,
double
s,
470
int
k1,
int
k2,
int
col_min=1,
int
col_max=0);
471
// do a row Givens update
472
473
void
givens
(
double
a
,
double
b
,
double
*c,
double
*s);
474
// algorithm 5.1.5 in Golub and Van Loan
475
476
HepVector
house
(
const
HepMatrix &
a
,
int
row=1,
int
col=1);
477
// Returns a Householder vector to zero elements.
478
479
void
house_with_update
(HepMatrix *
a
,
int
row=1,
int
col=1);
480
void
house_with_update
(HepMatrix *
a
, HepMatrix *v,
int
row=1,
int
col=1);
481
// Finds and does Householder reflection on matrix.
482
483
void
row_house
(HepMatrix *
a
,
const
HepVector &v,
double
vnormsq,
484
int
row=1,
int
col=1);
485
void
row_house
(HepMatrix *
a
,
const
HepMatrix &v,
double
vnormsq,
486
int
row,
int
col,
int
row_start,
int
col_start);
487
void
row_house
(HepMatrix *
a
,
const
HepMatrix &v,
int
row,
int
col,
488
int
row_start,
int
col_start);
489
// Does a row Householder update.
490
491
}
// namespace CLHEP
492
493
#ifdef ENABLE_BACKWARDS_COMPATIBILITY
494
// backwards compatibility will be enabled ONLY in CLHEP 1.9
495
using namespace
CLHEP;
496
#endif
497
498
#ifndef HEP_DEBUG_INLINE
499
#include "CLHEP/Matrix/Matrix.icc"
500
#endif
501
502
#endif
/*_Matrix_H*/
Generated on Sun Jun 17 2012 08:08:26 for CLHEP by
1.8.1.1