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

Matrix/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 
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 &hm1);
226  // Copy constructor.
227 
228  HepMatrix(const HepSymMatrix &);
229  HepMatrix(const HepDiagMatrix &);
230  HepMatrix(const HepVector &);
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 &);
254  HepMatrix & operator += ( const HepSymMatrix &);
256  HepMatrix & operator += ( const HepVector &);
257  HepMatrix & operator -= ( const HepMatrix &);
258  HepMatrix & operator -= ( const HepSymMatrix &);
260  HepMatrix & operator -= ( const HepVector &);
261  // Add or subtract a Matrix.
262  // When adding/subtracting Vector, Matrix must have num_col of one.
263 
264  HepMatrix & operator = ( const HepMatrix &);
265  HepMatrix & operator = ( const HepSymMatrix &);
266  HepMatrix & operator = ( const HepDiagMatrix &);
267  HepMatrix & operator = ( const HepVector &);
268  HepMatrix & operator = ( const HepRotation &);
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 &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 
314  public:
315  inline HepMatrix_row(HepMatrix&,int);
316  double & operator[](int);
317  private:
318  HepMatrix& _a;
319  int _r;
320  };
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 &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);
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 &, const HepMatrix &);
411 HepMatrix operator*(double t, const HepMatrix &);
412 HepMatrix operator*(const HepMatrix &, double );
413 // Multiplication operators
414 // Note that m *= hm1 is always faster than m = m * hm1.
415 
416 HepMatrix operator/(const HepMatrix &, double );
417 // m = hm1 / t. (m /= t is faster if you can use it.)
418 
419 HepMatrix operator+(const HepMatrix &hm1, const HepMatrix &hm2);
420 // m = hm1 + hm2;
421 // Note that m += hm1 is always faster than m = m + hm1.
422 
423 HepMatrix operator-(const HepMatrix &hm1, const HepMatrix &hm2);
424 // m = hm1 - hm2;
425 // Note that m -= hm1 is always faster than m = m - hm1.
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);
444 HepMatrix 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 
449 HepMatrix qr_inverse(const 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);
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*/
HepMatrix apply(double(*f)(double, int, int)) const
Definition: Matrix.cc:476
HepVector solve(const HepMatrix &, const HepVector &)
Definition: Vector.cc:576
HepMatrix operator-() const
Definition: Matrix.cc:261
Hep3Vector operator+(const Hep3Vector &, const Hep3Vector &)
HepMatrix_row(HepMatrix &, int)
void row_givens(HepMatrix *A, double c, double s, int k1, int k2, int col_min=1, int col_max=0)
friend void col_givens(HepMatrix *A, double c, double s, int k1, int k2, int rowmin, int rowmax)
friend HepVector house(const HepMatrix &a, int row, int col)
virtual void invertHaywood5(int &ierr)
void row_house(HepMatrix *a, const HepVector &v, double vnormsq, int row=1, int col=1)
void house_with_update(HepMatrix *a, int row=1, int col=1)
HepMatrix qr_inverse(const HepMatrix &A)
HepMatrix & operator=(const HepMatrix &)
Definition: Matrix.cc:417
virtual int num_col() const
Definition: Matrix.cc:122
virtual const double & operator()(int row, int col) const
Definition: Matrix.cc:137
HepMatrix_row_const(const HepMatrix &, int)
void col_house(HepMatrix *a, const HepMatrix &v, double vnormsq, int row, int col, int row_start, int col_start)
friend void tridiagonal(HepSymMatrix *a, HepMatrix *hsm)
HepLorentzVector operator/(const HepLorentzVector &, double a)
friend void row_givens(HepMatrix *A, double c, double s, int k1, int k2, int colmin, int colmax)
friend HepMatrix operator+(const HepMatrix &hm1, const HepMatrix &hm2)
Definition: Matrix.cc:278
void qr_decomp(HepMatrix *A, HepMatrix *hsm)
friend HepMatrix operator*(const HepMatrix &hm1, const HepMatrix &hm2)
Definition: Matrix.cc:351
friend void house_with_update2(HepSymMatrix *a, HepMatrix *v, int row, int col)
double determinant() const
Definition: Matrix.cc:815
HepVector house(const HepMatrix &a, int row=1, int col=1)
friend void row_house(HepMatrix *, const HepMatrix &, double, int, int, int, int)
HepLorentzRotation operator*(const HepRotation &r, const HepLorentzRotation &lt)
friend void col_house(HepMatrix *, const HepMatrix &, double, int, int, int, int)
virtual void invertHaywood6(int &ierr)
virtual ~HepMatrix()
Definition: Matrix.cc:108
HepMatrix inverse() const
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)
const double & operator[](int) const
HepMatrix & operator/=(double t)
Definition: Matrix.cc:405
void back_solve(const HepMatrix &R, HepVector *b)
Definition: MatrixLinear.cc:60
HepMatrix & operator*=(double t)
Definition: Matrix.cc:411
HepMatrix_row operator[](int)
void givens(double a, double b, double *c, double *s)
friend HepVector solve(const HepMatrix &, const HepVector &)
Definition: Vector.cc:576
friend void swap(HepMatrix &hm1, HepMatrix &hm2)
void f(void g())
Definition: excDblThrow.cc:38
HepVector qr_solve(const HepMatrix &A, const HepVector &b)
virtual void invertHaywood4(int &ierr)
Hep3Vector operator-(const Hep3Vector &, const Hep3Vector &)
HepMatrix sub(int min_row, int max_row, int min_col, int max_col) const
Definition: Matrix.cc:195
HepMatrix & operator-=(const HepMatrix &)
Definition: Matrix.cc:398
HepMatrix T() const
Definition: Matrix.cc:456
friend void house_with_update(HepMatrix *a, int row, int col)
HepMatrix & operator+=(const HepMatrix &)
Definition: Matrix.cc:391
virtual int num_row() const
Definition: Matrix.cc:120
virtual int num_size() const
Definition: Matrix.cc:124
std::ostream & operator<<(std::ostream &os, const HepAxisAngle &aa)
Definition: AxisAngle.cc:86
double trace() const
Definition: Matrix.cc:832
friend void back_solve(const HepMatrix &R, HepVector *b)
Definition: MatrixLinear.cc:60
friend HepVector qr_solve(HepMatrix *, const HepVector &)