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
inline
void
invert
();
301
// Invert a matrix. Throw std::runtime_error on failure.
302
303
inline
HepMatrix
inverse
()
const
;
304
// Invert a matrix. Throw std::runtime_error on failure.
305
306
307
double
determinant
()
const
;
308
// calculate the determinant of the matrix.
309
310
double
trace
()
const
;
311
// calculate the trace of the matrix (sum of diagonal elements).
312
313
class
HepMatrix_row
{
314
public
:
315
inline
HepMatrix_row
(
HepMatrix
&,
int
);
316
double
&
operator[]
(
int
);
317
private
:
318
HepMatrix
& _a;
319
int
_r;
320
};
321
class
HepMatrix_row_const
{
322
public
:
323
inline
HepMatrix_row_const
(
const
HepMatrix
&,
int
);
324
const
double
&
operator[]
(
int
)
const
;
325
private
:
326
const
HepMatrix
& _a;
327
int
_r;
328
};
329
// helper classes for implementing m[i][j]
330
331
inline
HepMatrix_row
operator[]
(
int
);
332
inline
const
HepMatrix_row_const
operator[]
(
int
)
const
;
333
// Read or write a matrix element.
334
// While it may not look like it, you simply do m[i][j] to get an
335
// element.
336
// ** Note that the indexing starts from [0][0]. **
337
338
protected
:
339
virtual
int
num_size
()
const
;
340
virtual
void
invertHaywood4
(
int
& ierr);
341
virtual
void
invertHaywood5
(
int
& ierr);
342
virtual
void
invertHaywood6
(
int
& ierr);
343
344
private
:
345
friend
class
HepMatrix_row
;
346
friend
class
HepMatrix_row_const
;
347
friend
class
HepVector
;
348
friend
class
HepSymMatrix
;
349
friend
class
HepDiagMatrix
;
350
// Friend classes.
351
352
friend
HepMatrix
operator+
(
const
HepMatrix
&m1,
const
HepMatrix
&m2);
353
friend
HepMatrix
operator-
(
const
HepMatrix
&m1,
const
HepMatrix
&m2);
354
friend
HepMatrix
operator*
(
const
HepMatrix
&m1,
const
HepMatrix
&m2);
355
friend
HepMatrix
operator*
(
const
HepMatrix
&m1,
const
HepSymMatrix
&m2);
356
friend
HepMatrix
operator*
(
const
HepMatrix
&m1,
const
HepDiagMatrix
&m2);
357
friend
HepMatrix
operator*
(
const
HepSymMatrix
&m1,
const
HepMatrix
&m2);
358
friend
HepMatrix
operator*
(
const
HepDiagMatrix
&m1,
const
HepMatrix
&m2);
359
friend
HepMatrix
operator*
(
const
HepVector
&m1,
const
HepMatrix
&m2);
360
friend
HepVector
operator*
(
const
HepMatrix
&m1,
const
HepVector
&m2);
361
friend
HepMatrix
operator*
(
const
HepSymMatrix
&m1,
const
HepSymMatrix
&m2);
362
// Multiply a Matrix by a Matrix or Vector.
363
364
friend
HepVector
solve
(
const
HepMatrix
&,
const
HepVector
&);
365
// solve the system of linear eq
366
friend
HepVector
qr_solve
(
HepMatrix
*,
const
HepVector
&);
367
friend
HepMatrix
qr_solve
(
HepMatrix
*,
const
HepMatrix
&
b
);
368
friend
void
tridiagonal
(
HepSymMatrix
*
a
,
HepMatrix
*hsm);
369
friend
void
row_house
(
HepMatrix
*,
const
HepMatrix
&,
double
,
370
int
,
int
,
int
,
int
);
371
friend
void
row_house
(
HepMatrix
*,
const
HepVector
&,
double
,
372
int
,
int
);
373
friend
void
back_solve
(
const
HepMatrix
&R,
HepVector
*
b
);
374
friend
void
back_solve
(
const
HepMatrix
&R,
HepMatrix
*
b
);
375
friend
void
col_givens
(
HepMatrix
*A,
double
c,
376
double
s,
int
k1,
int
k2,
377
int
rowmin,
int
rowmax);
378
// Does a column Givens update.
379
friend
void
row_givens
(
HepMatrix
*A,
double
c,
380
double
s,
int
k1,
int
k2,
381
int
colmin,
int
colmax);
382
friend
void
col_house
(
HepMatrix
*,
const
HepMatrix
&,
double
,
383
int
,
int
,
int
,
int
);
384
friend
HepVector
house
(
const
HepMatrix
&
a
,
int
row,
int
col);
385
friend
void
house_with_update
(
HepMatrix
*
a
,
int
row,
int
col);
386
friend
void
house_with_update
(
HepMatrix
*
a
,
HepMatrix
*v,
int
row,
int
col);
387
friend
void
house_with_update2
(
HepSymMatrix
*
a
,
HepMatrix
*v,
388
int
row,
int
col);
389
390
int
dfact_matrix(
double
&det,
int
*ir);
391
// factorize the matrix. If successful, the return code is 0. On
392
// return, det is the determinant and ir[] is row-interchange
393
// matrix. See CERNLIB's DFACT routine.
394
395
int
dfinv_matrix(
int
*ir);
396
// invert the matrix. See CERNLIB DFINV.
397
398
#ifdef DISABLE_ALLOC
399
std::vector<double > m;
400
#else
401
std::vector<double,Alloc<double,25> > m;
402
#endif
403
int
nrow, ncol;
404
int
size_;
405
};
406
407
// Operations other than member functions for Matrix
408
// implemented in Matrix.cc and Matrix.icc (inline).
409
410
HepMatrix
operator*
(
const
HepMatrix &m1,
const
HepMatrix &m2);
411
HepMatrix
operator*
(
double
t,
const
HepMatrix &m1);
412
HepMatrix
operator*
(
const
HepMatrix &m1,
double
t);
413
// Multiplication operators
414
// Note that m *= m1 is always faster than m = m * m1.
415
416
HepMatrix
operator/
(
const
HepMatrix &m1,
double
t);
417
// m = m1 / t. (m /= t is faster if you can use it.)
418
419
HepMatrix
operator+
(
const
HepMatrix &m1,
const
HepMatrix &m2);
420
// m = m1 + m2;
421
// Note that m += m1 is always faster than m = m + m1.
422
423
HepMatrix
operator-
(
const
HepMatrix &m1,
const
HepMatrix &m2);
424
// m = m1 - m2;
425
// Note that m -= m1 is always faster than m = m - m1.
426
427
HepMatrix
dsum
(
const
HepMatrix&,
const
HepMatrix&);
428
// Direct sum of two matrices. The direct sum of A and B is the matrix
429
// A 0
430
// 0 B
431
432
HepVector
solve
(
const
HepMatrix &,
const
HepVector &);
433
// solve the system of linear equations using LU decomposition.
434
435
std::ostream&
operator<<
(std::ostream &s,
const
HepMatrix &q);
436
// Read in, write out Matrix into a stream.
437
438
//
439
// Specialized linear algebra functions
440
//
441
442
HepVector
qr_solve
(
const
HepMatrix &A,
const
HepVector &
b
);
443
HepVector
qr_solve
(HepMatrix *A,
const
HepVector &
b
);
444
HepMatrix
qr_solve
(
const
HepMatrix &A,
const
HepMatrix &
b
);
445
HepMatrix
qr_solve
(HepMatrix *A,
const
HepMatrix &
b
);
446
// Works like backsolve, except matrix does not need to be upper
447
// triangular. For nonsquare matrix, it solves in the least square sense.
448
449
HepMatrix
qr_inverse
(
const
HepMatrix &A);
450
HepMatrix
qr_inverse
(HepMatrix *A);
451
// Finds the inverse of a matrix using QR decomposition. Note, often what
452
// you really want is solve or backsolve, they can be much quicker than
453
// inverse in many calculations.
454
455
456
void
qr_decomp
(HepMatrix *A, HepMatrix *hsm);
457
HepMatrix
qr_decomp
(HepMatrix *A);
458
// Does a QR decomposition of a matrix.
459
460
void
back_solve
(
const
HepMatrix &R, HepVector *
b
);
461
void
back_solve
(
const
HepMatrix &R, HepMatrix *
b
);
462
// Solves R*x = b where R is upper triangular. Also has a variation that
463
// solves a number of equations of this form in one step, where b is a matrix
464
// with each column a different vector. See also solve.
465
466
void
col_house
(HepMatrix *
a
,
const
HepMatrix &v,
double
vnormsq,
467
int
row,
int
col,
int
row_start,
int
col_start);
468
void
col_house
(HepMatrix *
a
,
const
HepMatrix &v,
int
row,
int
col,
469
int
row_start,
int
col_start);
470
// Does a column Householder update.
471
472
void
col_givens
(HepMatrix *A,
double
c,
double
s,
473
int
k1,
int
k2,
int
row_min=1,
int
row_max=0);
474
// do a column Givens update
475
476
void
row_givens
(HepMatrix *A,
double
c,
double
s,
477
int
k1,
int
k2,
int
col_min=1,
int
col_max=0);
478
// do a row Givens update
479
480
void
givens
(
double
a
,
double
b
,
double
*c,
double
*s);
481
// algorithm 5.1.5 in Golub and Van Loan
482
483
HepVector
house
(
const
HepMatrix &
a
,
int
row=1,
int
col=1);
484
// Returns a Householder vector to zero elements.
485
486
void
house_with_update
(HepMatrix *
a
,
int
row=1,
int
col=1);
487
void
house_with_update
(HepMatrix *
a
, HepMatrix *v,
int
row=1,
int
col=1);
488
// Finds and does Householder reflection on matrix.
489
490
void
row_house
(HepMatrix *
a
,
const
HepVector &v,
double
vnormsq,
491
int
row=1,
int
col=1);
492
void
row_house
(HepMatrix *
a
,
const
HepMatrix &v,
double
vnormsq,
493
int
row,
int
col,
int
row_start,
int
col_start);
494
void
row_house
(HepMatrix *
a
,
const
HepMatrix &v,
int
row,
int
col,
495
int
row_start,
int
col_start);
496
// Does a row Householder update.
497
498
}
// namespace CLHEP
499
500
#ifdef ENABLE_BACKWARDS_COMPATIBILITY
501
// backwards compatibility will be enabled ONLY in CLHEP 1.9
502
using namespace
CLHEP;
503
#endif
504
505
#ifndef HEP_DEBUG_INLINE
506
#include "CLHEP/Matrix/Matrix.icc"
507
#endif
508
509
#endif
/*_Matrix_H*/
Generated on Mon May 6 2013 04:04:11 for CLHEP by
1.8.1.2