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 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*/