CLHEP/Matrix/Matrix.h

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

Class Library for High Energy Physics (version 1.8)