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

RandomObjects/CLHEP/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 hm1(n, m);
19//
20// To declare a Matrix as a copy of a Matrix hm2, say
21//
22// .ft B
23// HepMatrix hm1(hm2);
24// or
25// .ft B
26// HepMatrix hm1 = hm2;
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 hm1(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 = hm1.num_row();
47//
48// or
49//
50// .ft B
51// nc = hm1.num_col();
52//
53// You can print a Matrix by
54//
55// .ft B
56// cout << hm1;
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// hm1 *= hm2 is as hm1 = hm1*hm2. 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 = hm2.sub(row_min, row_max, col_min, col_max);
121//
122// returns the subMatrix.
123// hm2(row_min:row_max, col_min:col_max) in matlab terminology.
124// If instead you say
125//
126// .ft B
127// hm2.sub(row, col, hm1);
128//
129// then the subMatrix
130// hm2(row:row+hm1.num_row()-1, col:col+hm1.num_col()-1) is replaced with hm1.
131//
132// .ft B
133// md = dsum(hm1,hm2);
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 = hm1.apply(HEP_MATRIX_ELEMENT (*f)(HEP_MATRIX_ELEMENT, int r, int c));
177//
178// applies f to every element of hm1 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
196namespace CLHEP {
197
198class HepRandom;
199
200class HepSymMatrix;
201class HepDiagMatrix;
202class HepVector;
203class HepRotation;
204
209class HepMatrix : public HepGenMatrix {
210public:
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 &hm1);
226 // Copy constructor.
227
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
248 // Multiply a Matrix by a floating number.
249
251 // Divide a Matrix by a floating number.
252
261 // Add or subtract a Matrix.
262 // When adding/subtracting Vector, Matrix must have num_col of one.
263
269 // Assignment operators.
270
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 &hm1);
284 // Sub matrix of this Matrix is replaced with hm1.
285 // WARNING: rows and columns are numbered from 1
286
287 friend inline void swap(HepMatrix &hm1, HepMatrix &hm2);
288 // Swap hm1 with hm2.
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:
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
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
338protected:
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
344private:
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 &hm1, const HepMatrix &hm2);
353 friend HepMatrix operator-(const HepMatrix &hm1, const HepMatrix &hm2);
354 friend HepMatrix operator*(const HepMatrix &hm1, const HepMatrix &hm2);
355 friend HepMatrix operator*(const HepMatrix &hm1, const HepSymMatrix &hm2);
356 friend HepMatrix operator*(const HepMatrix &hm1, const HepDiagMatrix &hm2);
357 friend HepMatrix operator*(const HepSymMatrix &hm1, const HepMatrix &hm2);
358 friend HepMatrix operator*(const HepDiagMatrix &hm1, const HepMatrix &hm2);
359 friend HepMatrix operator*(const HepVector &hm1, const HepMatrix &hm2);
360 friend HepVector operator*(const HepMatrix &hm1, const HepVector &hm2);
361 friend HepMatrix operator*(const HepSymMatrix &hm1, const HepSymMatrix &hm2);
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);
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
410HepMatrix operator*(const HepMatrix &, const HepMatrix &);
411HepMatrix operator*(double t, const HepMatrix &);
412HepMatrix operator*(const HepMatrix &, double );
413// Multiplication operators
414// Note that m *= hm1 is always faster than m = m * hm1.
415
416HepMatrix operator/(const HepMatrix &, double );
417// m = hm1 / t. (m /= t is faster if you can use it.)
418
419HepMatrix operator+(const HepMatrix &hm1, const HepMatrix &hm2);
420// m = hm1 + hm2;
421// Note that m += hm1 is always faster than m = m + hm1.
422
423HepMatrix operator-(const HepMatrix &hm1, const HepMatrix &hm2);
424// m = hm1 - hm2;
425// Note that m -= hm1 is always faster than m = m - hm1.
426
427HepMatrix 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
432HepVector solve(const HepMatrix &, const HepVector &);
433// solve the system of linear equations using LU decomposition.
434
435std::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
442HepVector qr_solve(const HepMatrix &A, const HepVector &b);
444HepMatrix qr_solve(const 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
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
456void qr_decomp(HepMatrix *A, HepMatrix *hsm);
458// Does a QR decomposition of a matrix.
459
460void back_solve(const HepMatrix &R, HepVector *b);
461void 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
466void col_house(HepMatrix *a, const HepMatrix &v, double vnormsq,
467 int row, int col, int row_start, int col_start);
468void 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
472void 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
476void 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
480void givens(double a, double b, double *c, double *s);
481// algorithm 5.1.5 in Golub and Van Loan
482
483HepVector house(const HepMatrix &a, int row=1, int col=1);
484// Returns a Householder vector to zero elements.
485
486void house_with_update(HepMatrix *a, int row=1, int col=1);
487void house_with_update(HepMatrix *a, HepMatrix *v, int row=1, int col=1);
488// Finds and does Householder reflection on matrix.
489
490void row_house(HepMatrix *a, const HepVector &v, double vnormsq,
491 int row=1, int col=1);
492void row_house(HepMatrix *a, const HepMatrix &v, double vnormsq,
493 int row, int col, int row_start, int col_start);
494void 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
502using namespace CLHEP;
503#endif
504
505#ifndef HEP_DEBUG_INLINE
506#include "CLHEP/Matrix/Matrix.icc"
507#endif
508
509#endif /*_Matrix_H*/
HepMatrix_row_const(const HepMatrix &, int)
const double & operator[](int) const
HepMatrix_row(HepMatrix &, int)
friend HepVector house(const HepMatrix &a, int row, int col)
HepMatrix(int p, int q, int i)
friend void house_with_update2(HepSymMatrix *a, HepMatrix *v, int row, int col)
virtual void invertHaywood6(int &ierr)
HepMatrix(const HepVector &)
HepMatrix(int p, int q)
void sub(int row, int col, const HepMatrix &hm1)
friend HepMatrix operator+(const HepMatrix &hm1, const HepMatrix &hm2)
Definition: Matrix.cc:278
virtual int num_size() const
HepMatrix sub(int min_row, int max_row, int min_col, int max_col) const
virtual ~HepMatrix()
virtual int num_col() const
HepMatrix & operator=(const HepMatrix &)
Definition: Matrix.cc:417
HepMatrix(const HepMatrix &hm1)
friend void row_house(HepMatrix *, const HepMatrix &, double, int, int, int, int)
friend void tridiagonal(HepSymMatrix *a, HepMatrix *hsm)
friend HepVector solve(const HepMatrix &, const HepVector &)
Definition: Vector.cc:576
HepMatrix apply(double(*f)(double, int, int)) const
friend void col_givens(HepMatrix *A, double c, double s, int k1, int k2, int rowmin, int rowmax)
HepMatrix(const HepSymMatrix &)
HepMatrix & operator/=(double t)
Definition: Matrix.cc:405
virtual const double & operator()(int row, int col) const
friend void row_givens(HepMatrix *A, double c, double s, int k1, int k2, int colmin, int colmax)
HepMatrix & operator+=(const HepMatrix &)
Definition: Matrix.cc:391
virtual void invert(int &ierr)
HepMatrix_row operator[](int)
virtual double & operator()(int row, int col)
friend void swap(HepMatrix &hm1, HepMatrix &hm2)
HepMatrix operator-() const
Definition: Matrix.cc:261
virtual int num_row() const
HepMatrix & operator*=(double t)
Definition: Matrix.cc:411
virtual void invertHaywood4(int &ierr)
friend void house_with_update(HepMatrix *a, int row, int col)
HepMatrix(int p, int q, HepRandom &r)
virtual void invertHaywood5(int &ierr)
HepMatrix & operator-=(const HepMatrix &)
Definition: Matrix.cc:398
friend void col_house(HepMatrix *, const HepMatrix &, double, int, int, int, int)
HepMatrix T() const
double determinant() const
friend HepMatrix operator*(const HepMatrix &hm1, const HepMatrix &hm2)
Definition: Matrix.cc:351
friend HepVector qr_solve(HepMatrix *, const HepVector &)
HepMatrix(const HepDiagMatrix &)
HepMatrix inverse(int &ierr) const
HepMatrix inverse() const
friend void back_solve(const HepMatrix &R, HepVector *b)
Definition: MatrixLinear.cc:60
double trace() const
void f(void g())
Definition: excDblThrow.cc:38
void col_house(HepMatrix *a, const HepMatrix &v, double vnormsq, int row, int col, int row_start, int col_start)
HepVector qr_solve(const HepMatrix &A, const HepVector &b)
Hep3Vector operator-(const Hep3Vector &, const Hep3Vector &)
Hep3Vector operator+(const Hep3Vector &, const Hep3Vector &)
HepMatrix qr_inverse(const HepMatrix &A)
void house_with_update(HepMatrix *a, int row=1, int col=1)
HepLorentzRotation operator*(const HepRotation &r, const HepLorentzRotation &lt)
void back_solve(const HepMatrix &R, HepVector *b)
Definition: MatrixLinear.cc:60
void qr_decomp(HepMatrix *A, HepMatrix *hsm)
HepVector house(const HepMatrix &a, int row=1, int col=1)
void row_givens(HepMatrix *A, double c, double s, int k1, int k2, int col_min=1, int col_max=0)
std::ostream & operator<<(std::ostream &os, const HepAxisAngle &aa)
Definition: AxisAngle.cc:86
void givens(double a, double b, double *c, double *s)
HepDiagMatrix dsum(const HepDiagMatrix &s1, const HepDiagMatrix &s2)
Definition: DiagMatrix.cc:164
void col_givens(HepMatrix *A, double c, double s, int k1, int k2, int row_min=1, int row_max=0)
void row_house(HepMatrix *a, const HepVector &v, double vnormsq, int row=1, int col=1)
HepLorentzVector operator/(const HepLorentzVector &, double a)
HepVector solve(const HepMatrix &, const HepVector &)
Definition: Vector.cc:576
@ b
@ a