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

RandomObjects/CLHEP/Matrix/Matrix.h
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 // CLASSDOC OFF
00003 // ---------------------------------------------------------------------------
00004 // CLASSDOC ON
00005 //
00006 // This file is a part of the CLHEP - a Class Library for High Energy Physics.
00007 // 
00008 // This is the definition of the HepMatrix class.
00009 // 
00010 // This software written by Nobu Katayama and Mike Smyth, Cornell University.
00011 //
00012 //
00013 // .SS Usage
00014 // The Matrix class does all the obvious things, in all the obvious ways.
00015 // You declare a Matrix by saying
00016 //
00017 // .ft B
00018 //       HepMatrix m1(n, m);
00019 //
00020 //  To declare a Matrix as a copy of a Matrix m2, say
00021 //
00022 // .ft B
00023 //       HepMatrix m1(m2);
00024 // or
00025 // .ft B
00026 //       HepMatrix m1 = m2;
00027 // 
00028 // You can declare initilizations of a Matrix, by giving it a third
00029 // integer argument, either 0 or 1. 0 means initialize to 0, one means
00030 // the identity matrix. If you do not give a third argument the memory
00031 // is not initialized.
00032 //
00033 // .ft B
00034 //       HepMatrix m1(n, m, 1);
00035 //
00036 // ./"This code has been written by Mike Smyth, and the algorithms used are
00037 // ./"described in the thesis "A Tracking Library for a Silicon Vertex Detector"
00038 // ./"(Mike Smyth, Cornell University, June 1993).
00039 // ./"This is file contains C++ stuff for doing things with Matrices.
00040 // ./"To turn on bound checking, define MATRIX_BOUND_CHECK before including
00041 // ./"this file.
00042 // 
00043 //  To find the number of rows or number of columns, say 
00044 //
00045 // .ft B
00046 // nr = m1.num_row();
00047 //
00048 // or
00049 //
00050 // .ft B
00051 // nc = m1.num_col();
00052 //
00053 // You can print a Matrix by
00054 //
00055 // .ft B
00056 // cout << m1;
00057 //
00058 //  You can add,
00059 //  subtract, and multiply two Matrices.  You can multiply a Matrix by a
00060 //  scalar, on the left or the right.  +=, *=, -= act in the obvious way.
00061 //  m1 *= m2 is as m1 = m1*m2. You can also divide by a scalar on the
00062 //  right, or use /= to do the same thing.  
00063 // 
00064 //  You can read or write a Matrix element by saying
00065 //
00066 // .ft B
00067 //  m(row, col) = blah. (1 <= row <= num_row(), 1 <= col <= num_col())
00068 //
00069 // .ft B
00070 //  blah = m(row, col) ...
00071 //
00072 //  m(row, col) is inline, and by default does not do bound checking. 
00073 //  If bound checking is desired, say #define MATRIX_BOUND_CHECK before
00074 //  including matrix.h.
00075 // 
00076 //  You can also access the element using C subscripting operator []
00077 //
00078 // .ft B
00079 //  m[row][col] = blah. (0 <= row < num_row(), 0 <= col < num_col())
00080 //
00081 // .ft B
00082 //  blah = m[row][col] ...
00083 //
00084 //  m[row][col] is inline, and by default does not do bound checking. 
00085 //  If bound checking is desired, say #define MATRIX_BOUND_CHECK before
00086 //  including matrix.h. (Notice the difference in bases in two access
00087 //  methods.) 
00088 //
00089 // .SS Comments on precision
00090 //
00091 //  The user would normally use "Matrix" class whose precision is the same
00092 //  as the other classes of CLHEP (ThreeVec, for example). However, he/she
00093 //  can explicitly choose Matrix (double) or MatrixF (float) if he/she wishes.
00094 //  In the former case, include "Matrix.h." In the latter case, include either
00095 //  "Matrix.h," or "MatrixF.h," or both. The only operators that can talk
00096 //  to each other are the copy constructor and assignment operator.
00097 //
00098 // .ft B
00099 //  Matrix d(3,4,HEP_MATRIX_IDENTITY);
00100 //
00101 // .ft B
00102 //  MatrixF f;
00103 //
00104 // .ft B
00105 //  f = d;
00106 //
00107 // .ft B
00108 //  MatrixF g(d);
00109 //
00110 //  will convert from one to the other.
00111 //
00112 // .SS Other functions
00113 //
00114 // .ft B
00115 //  mt = m.T();
00116 //
00117 //  returns the transpose of m. 
00118 //
00119 // .ft B
00120 //  ms = m2.sub(row_min, row_max, col_min, col_max);
00121 //
00122 //  returns the subMatrix.
00123 //  m2(row_min:row_max, col_min:col_max) in matlab terminology.
00124 //  If instead you say
00125 //
00126 // .ft B
00127 //  m2.sub(row, col, m1);
00128 //
00129 //  then the subMatrix
00130 //  m2(row:row+m1.num_row()-1, col:col+m1.num_col()-1) is replaced with m1.
00131 //
00132 // .ft B
00133 // md = dsum(m1,m2);
00134 //
00135 // returns the direct sum of the two matrices.
00136 //
00137 // Operations that do not have to be member functions or friends are listed
00138 // towards the end of this man page.
00139 //
00140 //
00141 // To invert a matrix, say
00142 //
00143 // .ft B
00144 // minv = m.inverse(ierr);
00145 //
00146 // ierr will be non-zero if the matrix is singular.
00147 //
00148 // If you can overwrite the matrix, you can use the invert function to avoid
00149 // two extra copies. Use
00150 //
00151 // .ft B
00152 // m.invert(ierr);
00153 //
00154 // Many basic linear algebra functions such as QR decomposition are defined.
00155 // In order to keep the header file a reasonable size, the functions are
00156 // defined in MatrixLinear.h.
00157 //
00158 // 
00159 // .ft B 
00160 //  Note that Matrices run from (1,1) to (n, m), and [0][0] to
00161 //  [n-1][m-1]. The users of the latter form should be careful about sub
00162 //  functions.
00163 //
00164 // ./" The program that this class was orginally used with lots of small
00165 // ./" Matrices.  It was very costly to malloc and free array space every
00166 // ./" time a Matrix is needed.  So, a number of previously freed arrays are
00167 // ./" kept around, and when needed again one of these array is used.  Right
00168 // ./" now, a set of piles of arrays with row <= row_max and col <= col_max
00169 // ./" are kept around.  The piles of kept Matrices can be easily changed.
00170 // ./" Array mallocing and freeing are done only in new_m() and delete_m(),
00171 // ./" so these are the only routines that need to be rewritten.
00172 // 
00173 //  You can do things thinking of a Matrix as a list of numbers.  
00174 //
00175 // .ft B
00176 //  m = m1.apply(HEP_MATRIX_ELEMENT (*f)(HEP_MATRIX_ELEMENT, int r, int c));
00177 // 
00178 //  applies f to every element of m1 and places it in m.
00179 //
00180 // .SS See Also:
00181 // SymMatrix[DF].h, GenMatrix[DF].h, DiagMatrix[DF].h Vector[DF].h
00182 // MatrixLinear[DF].h 
00183 
00184 #ifndef _Matrix_H_
00185 #define _Matrix_H_
00186 
00187 #ifdef GNUPRAGMA
00188 #pragma interface
00189 #endif
00190 
00191 #include <vector>
00192 
00193 #include "CLHEP/Matrix/defs.h"
00194 #include "CLHEP/Matrix/GenMatrix.h"
00195 
00196 namespace CLHEP {
00197 
00198 class HepRandom;
00199 
00200 class HepSymMatrix;
00201 class HepDiagMatrix;
00202 class HepVector;
00203 class HepRotation;
00204 
00209 class HepMatrix : public HepGenMatrix {
00210 public:
00211    inline HepMatrix();
00212    // Default constructor. Gives 0 x 0 matrix. Another Matrix can be
00213    // assigned to it.
00214 
00215    HepMatrix(int p, int q);
00216    // Constructor. Gives an unitialized p x q matrix.
00217    HepMatrix(int p, int q, int i);
00218    // Constructor. Gives an initialized p x q matrix. 
00219    // If i=0, it is initialized to all 0. If i=1, the diagonal elements
00220    // are set to 1.0.
00221 
00222    HepMatrix(int p, int q, HepRandom &r);
00223    // Constructor with a Random object.
00224 
00225    HepMatrix(const HepMatrix &m1);
00226    // Copy constructor.
00227 
00228    HepMatrix(const HepSymMatrix &m1);
00229    HepMatrix(const HepDiagMatrix &m1);
00230    HepMatrix(const HepVector &m1);
00231    // Constructors from SymMatrix, DiagMatrix and Vector.
00232 
00233    virtual ~HepMatrix();
00234    // Destructor.
00235 
00236    virtual int num_row() const;
00237    // Returns the number of rows.
00238 
00239    virtual int num_col() const;
00240    // Returns the number of columns.
00241 
00242    virtual const double & operator()(int row, int col) const;
00243    virtual double & operator()(int row, int col);
00244    // Read or write a matrix element. 
00245    // ** Note that the indexing starts from (1,1). **
00246 
00247    HepMatrix & operator *= (double t);
00248    // Multiply a Matrix by a floating number.
00249 
00250    HepMatrix & operator /= (double t); 
00251    // Divide a Matrix by a floating number.
00252 
00253    HepMatrix & operator += ( const HepMatrix &m2);
00254    HepMatrix & operator += ( const HepSymMatrix &m2);
00255    HepMatrix & operator += ( const HepDiagMatrix &m2);
00256    HepMatrix & operator += ( const HepVector &m2);
00257    HepMatrix & operator -= ( const HepMatrix &m2);
00258    HepMatrix & operator -= ( const HepSymMatrix &m2);
00259    HepMatrix & operator -= ( const HepDiagMatrix &m2);
00260    HepMatrix & operator -= ( const HepVector &m2);
00261    // Add or subtract a Matrix. 
00262    // When adding/subtracting Vector, Matrix must have num_col of one.
00263 
00264    HepMatrix & operator = ( const HepMatrix &m2);
00265    HepMatrix & operator = ( const HepSymMatrix &m2);
00266    HepMatrix & operator = ( const HepDiagMatrix &m2);
00267    HepMatrix & operator = ( const HepVector &m2);
00268    HepMatrix & operator = ( const HepRotation &m2);
00269    // Assignment operators.
00270 
00271    HepMatrix operator- () const;
00272    // unary minus, ie. flip the sign of each element.
00273 
00274    HepMatrix apply(double (*f)(double, int, int)) const;
00275    // Apply a function to all elements of the matrix.
00276 
00277    HepMatrix T() const;
00278    // Returns the transpose of a Matrix.
00279 
00280    HepMatrix sub(int min_row, int max_row, int min_col, int max_col) const;
00281    // Returns a sub matrix of a Matrix.
00282    // WARNING: rows and columns are numbered from 1
00283    void sub(int row, int col, const HepMatrix &m1);
00284    // Sub matrix of this Matrix is replaced with m1.
00285    // WARNING: rows and columns are numbered from 1
00286 
00287    friend inline void swap(HepMatrix &m1, HepMatrix &m2);
00288    // Swap m1 with m2.
00289 
00290    inline HepMatrix inverse(int& ierr) const;
00291    // Invert a Matrix. Matrix must be square and is not changed.
00292    // Returns ierr = 0 (zero) when successful, otherwise non-zero.
00293 
00294    virtual void invert(int& ierr);
00295    // Invert a Matrix. Matrix must be square.
00296    // N.B. the contents of the matrix are replaced by the inverse.
00297    // Returns ierr = 0 (zero) when successful, otherwise non-zero. 
00298    // This method has less overhead then inverse().
00299 
00300    double determinant() const;
00301    // calculate the determinant of the matrix.
00302 
00303    double trace() const;
00304    // calculate the trace of the matrix (sum of diagonal elements).
00305 
00306    class HepMatrix_row {
00307    public:
00308       inline HepMatrix_row(HepMatrix&,int);
00309       double & operator[](int);
00310    private:
00311       HepMatrix& _a;
00312       int _r;
00313    };
00314    class HepMatrix_row_const {
00315    public:
00316       inline HepMatrix_row_const (const HepMatrix&,int);
00317       const double & operator[](int) const;
00318    private:
00319       const HepMatrix& _a;
00320       int _r;
00321    };
00322    // helper classes for implementing m[i][j]
00323 
00324    inline HepMatrix_row operator[] (int);
00325    inline const HepMatrix_row_const operator[] (int) const;
00326    // Read or write a matrix element.
00327    // While it may not look like it, you simply do m[i][j] to get an
00328    // element. 
00329    // ** Note that the indexing starts from [0][0]. **
00330 
00331 protected:
00332    virtual  int num_size() const;
00333    virtual void invertHaywood4(int& ierr);
00334    virtual void invertHaywood5(int& ierr);
00335    virtual void invertHaywood6(int& ierr);
00336 
00337 private:
00338    friend class HepMatrix_row;
00339    friend class HepMatrix_row_const;
00340    friend class HepVector;
00341    friend class HepSymMatrix;
00342    friend class HepDiagMatrix;
00343    // Friend classes.
00344 
00345    friend HepMatrix operator+(const HepMatrix &m1, const HepMatrix &m2);
00346    friend HepMatrix operator-(const HepMatrix &m1, const HepMatrix &m2);
00347    friend HepMatrix operator*(const HepMatrix &m1, const HepMatrix &m2);
00348    friend HepMatrix operator*(const HepMatrix &m1, const HepSymMatrix &m2);
00349    friend HepMatrix operator*(const HepMatrix &m1, const HepDiagMatrix &m2);
00350    friend HepMatrix operator*(const HepSymMatrix &m1, const HepMatrix &m2);
00351    friend HepMatrix operator*(const HepDiagMatrix &m1, const HepMatrix &m2);
00352    friend HepMatrix operator*(const HepVector &m1, const HepMatrix &m2);
00353    friend HepVector operator*(const HepMatrix &m1, const HepVector &m2);
00354    friend HepMatrix operator*(const HepSymMatrix &m1, const HepSymMatrix &m2);
00355    // Multiply a Matrix by a Matrix or Vector.
00356 
00357    friend HepVector solve(const HepMatrix &, const HepVector &);
00358    // solve the system of linear eq
00359    friend HepVector qr_solve(HepMatrix *, const HepVector &);
00360    friend HepMatrix qr_solve(HepMatrix *, const HepMatrix &b);
00361    friend void tridiagonal(HepSymMatrix *a,HepMatrix *hsm);
00362    friend void row_house(HepMatrix *,const HepMatrix &, double,
00363                          int, int, int, int);
00364    friend void row_house(HepMatrix *,const HepVector &, double,
00365                          int, int);
00366    friend void back_solve(const HepMatrix &R, HepVector *b);
00367    friend void back_solve(const HepMatrix &R, HepMatrix *b);
00368    friend void col_givens(HepMatrix *A, double c,
00369                           double s, int k1, int k2, 
00370                           int rowmin, int rowmax);
00371    //    Does a column Givens update.
00372    friend void row_givens(HepMatrix *A, double c,
00373                           double s, int k1, int k2, 
00374                           int colmin, int colmax);
00375    friend void col_house(HepMatrix *,const HepMatrix &, double,
00376                          int, int, int, int);
00377    friend HepVector house(const HepMatrix &a,int row,int col);
00378    friend void house_with_update(HepMatrix *a,int row,int col);
00379    friend void house_with_update(HepMatrix *a,HepMatrix *v,int row,int col);
00380    friend void house_with_update2(HepSymMatrix *a,HepMatrix *v,
00381                                   int row,int col); 
00382 
00383    int dfact_matrix(double &det, int *ir);
00384    // factorize the matrix. If successful, the return code is 0. On
00385    // return, det is the determinant and ir[] is row-interchange
00386    // matrix. See CERNLIB's DFACT routine.
00387 
00388    int dfinv_matrix(int *ir);
00389    // invert the matrix. See CERNLIB DFINV.
00390 
00391 #ifdef DISABLE_ALLOC
00392    std::vector<double > m;
00393 #else
00394    std::vector<double,Alloc<double,25> > m;
00395 #endif
00396    int nrow, ncol;
00397    int size_;
00398 };
00399 
00400 // Operations other than member functions for Matrix
00401 // implemented in Matrix.cc and Matrix.icc (inline).
00402 
00403 HepMatrix operator*(const HepMatrix &m1, const HepMatrix &m2);
00404 HepMatrix operator*(double t, const HepMatrix &m1);
00405 HepMatrix operator*(const HepMatrix &m1, double t);
00406 // Multiplication operators
00407 // Note that m *= m1 is always faster than m = m * m1.
00408 
00409 HepMatrix operator/(const HepMatrix &m1, double t);
00410 // m = m1 / t. (m /= t is faster if you can use it.)
00411 
00412 HepMatrix operator+(const HepMatrix &m1, const HepMatrix &m2);
00413 // m = m1 + m2;
00414 // Note that m += m1 is always faster than m = m + m1.
00415 
00416 HepMatrix operator-(const HepMatrix &m1, const HepMatrix &m2);
00417 // m = m1 - m2;
00418 // Note that m -= m1 is always faster than m = m - m1.
00419 
00420 HepMatrix dsum(const HepMatrix&, const HepMatrix&);
00421 // Direct sum of two matrices. The direct sum of A and B is the matrix 
00422 //        A 0
00423 //        0 B
00424 
00425 HepVector solve(const HepMatrix &, const HepVector &);
00426 // solve the system of linear equations using LU decomposition.
00427 
00428 std::ostream& operator<<(std::ostream &s, const HepMatrix &q);
00429 // Read in, write out Matrix into a stream.
00430  
00431 //
00432 // Specialized linear algebra functions
00433 //
00434 
00435 HepVector qr_solve(const HepMatrix &A, const HepVector &b);
00436 HepVector qr_solve(HepMatrix *A, const HepVector &b);
00437 HepMatrix qr_solve(const HepMatrix &A, const HepMatrix &b);
00438 HepMatrix qr_solve(HepMatrix *A, const HepMatrix &b);
00439 // Works like backsolve, except matrix does not need to be upper
00440 // triangular. For nonsquare matrix, it solves in the least square sense.
00441 
00442 HepMatrix qr_inverse(const HepMatrix &A);
00443 HepMatrix qr_inverse(HepMatrix *A);
00444 // Finds the inverse of a matrix using QR decomposition.  Note, often what
00445 // you really want is solve or backsolve, they can be much quicker than
00446 // inverse in many calculations.
00447 
00448 
00449 void qr_decomp(HepMatrix *A, HepMatrix *hsm);
00450 HepMatrix qr_decomp(HepMatrix *A);
00451 // Does a QR decomposition of a matrix.
00452 
00453 void back_solve(const HepMatrix &R, HepVector *b);
00454 void back_solve(const HepMatrix &R, HepMatrix *b);
00455 // Solves R*x = b where R is upper triangular.  Also has a variation that
00456 // solves a number of equations of this form in one step, where b is a matrix
00457 // with each column a different vector. See also solve.
00458 
00459 void col_house(HepMatrix *a, const HepMatrix &v, double vnormsq,
00460                int row, int col, int row_start, int col_start);
00461 void col_house(HepMatrix *a, const HepMatrix &v, int row, int col,
00462                int row_start, int col_start);
00463 // Does a column Householder update.
00464 
00465 void col_givens(HepMatrix *A, double c, double s,
00466                 int k1, int k2, int row_min=1, int row_max=0);
00467 // do a column Givens update
00468 
00469 void row_givens(HepMatrix *A, double c, double s,
00470                 int k1, int k2, int col_min=1, int col_max=0);
00471 // do a row Givens update
00472 
00473 void givens(double a, double b, double *c, double *s);
00474 // algorithm 5.1.5 in Golub and Van Loan
00475 
00476 HepVector house(const HepMatrix &a, int row=1, int col=1);
00477 // Returns a Householder vector to zero elements.
00478 
00479 void house_with_update(HepMatrix *a, int row=1, int col=1);
00480 void house_with_update(HepMatrix *a, HepMatrix *v, int row=1, int col=1);
00481 // Finds and does Householder reflection on matrix.
00482 
00483 void row_house(HepMatrix *a, const HepVector &v, double vnormsq,
00484                int row=1, int col=1);
00485 void row_house(HepMatrix *a, const HepMatrix &v, double vnormsq,
00486                int row, int col, int row_start, int col_start);
00487 void row_house(HepMatrix *a, const HepMatrix &v, int row, int col,
00488                int row_start, int col_start);
00489 // Does a row Householder update.
00490 
00491 }  // namespace CLHEP
00492 
00493 #ifdef ENABLE_BACKWARDS_COMPATIBILITY
00494 //  backwards compatibility will be enabled ONLY in CLHEP 1.9
00495 using namespace CLHEP;
00496 #endif
00497 
00498 #ifndef HEP_DEBUG_INLINE
00499 #include "CLHEP/Matrix/Matrix.icc"
00500 #endif
00501 
00502 #endif /*_Matrix_H*/