CLHEP/Matrix/Matrix.h

00001 // -*- C++ -*-
00002 // CLASSDOC OFF
00003 // $Id: Matrix.h,v 1.14 2002/04/12 15:02:44 evc 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 "CLHEP/Matrix/GenMatrix.h"
00232 #ifdef HEP_USE_RANDOM
00233 class HepRandom;
00234 #endif
00235 
00236 #ifdef HEP_NO_INLINE_IN_DECLARATION
00237 #define inline
00238 #endif
00239 
00240 class HepSymMatrix;
00241 class HepDiagMatrix;
00242 class HepVector;
00243 #ifdef HEP_USE_VECTOR_MODULE
00244 class HepRotation;
00245 #endif
00246 
00251 class HepMatrix : public HepGenMatrix {
00252 public:
00253    inline HepMatrix();
00254    // Default constructor. Gives 0 x 0 matrix. Another Matrix can be
00255    // assigned to it.
00256 
00257    HepMatrix(int p, int q);
00258    // Constructor. Gives an unitialized p x q matrix.
00259    HepMatrix(int p, int q, int i);
00260    // Constructor. Gives an initialized p x q matrix. 
00261    // If i=0, it is initialized to all 0. If i=1, the diagonal elements
00262    // are set to 1.0.
00263 
00264 #ifdef HEP_USE_RANDOM
00265    HepMatrix(int p, int q, HepRandom &r);
00266    // Constructor with a Random object.
00267 #endif
00268 
00269    HepMatrix(const HepMatrix &m1);
00270    // Copy constructor.
00271 
00272    HepMatrix(const HepSymMatrix &m1);
00273    HepMatrix(const HepDiagMatrix &m1);
00274    HepMatrix(const HepVector &m1);
00275    // Constructors from SymMatrix, DiagMatrix and Vector.
00276 
00277    virtual ~HepMatrix();
00278    // Destructor.
00279 
00280    inline virtual int num_row() const;
00281    // Returns the number of rows.
00282 
00283    inline virtual int num_col() const;
00284    // Returns the number of columns.
00285 
00286    inline virtual const double & operator()(int row, int col) const;
00287    inline virtual double & operator()(int row, int col);
00288    // Read or write a matrix element. 
00289    // ** Note that the indexing starts from (1,1). **
00290 
00291    HepMatrix & operator *= (double t);
00292    // Multiply a Matrix by a floating number.
00293 
00294    HepMatrix & operator /= (double t); 
00295    // Divide a Matrix by a floating number.
00296 
00297    HepMatrix & operator += ( const HepMatrix &m2);
00298    HepMatrix & operator += ( const HepSymMatrix &m2);
00299    HepMatrix & operator += ( const HepDiagMatrix &m2);
00300    HepMatrix & operator += ( const HepVector &m2);
00301    HepMatrix & operator -= ( const HepMatrix &m2);
00302    HepMatrix & operator -= ( const HepSymMatrix &m2);
00303    HepMatrix & operator -= ( const HepDiagMatrix &m2);
00304    HepMatrix & operator -= ( const HepVector &m2);
00305    // Add or subtract a Matrix. 
00306    // When adding/subtracting Vector, Matrix must have num_col of one.
00307 
00308    HepMatrix & operator = ( const HepMatrix &m2);
00309    HepMatrix & operator = ( const HepSymMatrix &m2);
00310    HepMatrix & operator = ( const HepDiagMatrix &m2);
00311    HepMatrix & operator = ( const HepVector &m2);
00312 #ifdef HEP_USE_VECTOR_MODULE
00313    HepMatrix & operator = ( const HepRotation &m2);
00314 #endif
00315    // Assignment operators.
00316 
00317    HepMatrix operator- () const;
00318    // unary minus, ie. flip the sign of each element.
00319 
00320    HepMatrix apply(double (*f)(double, int, int)) const;
00321    // Apply a function to all elements of the matrix.
00322 
00323    HepMatrix T() const;
00324    // Returns the transpose of a Matrix.
00325 
00326    HepMatrix sub(int min_row, int max_row, int min_col, int max_col) const;
00327    // Returns a sub matrix of a Matrix.
00328    // WARNING: rows and columns are numbered from 1
00329    void sub(int row, int col, const HepMatrix &m1);
00330    // Sub matrix of this Matrix is replaced with m1.
00331    // WARNING: rows and columns are numbered from 1
00332 
00333    friend inline void swap(HepMatrix &m1, HepMatrix &m2);
00334    // Swap m1 with m2.
00335 
00336    inline HepMatrix inverse(int& ierr) const;
00337    // Invert a Matrix. Matrix must be square and is not changed.
00338    // Returns ierr = 0 (zero) when successful, otherwise non-zero.
00339 
00340    virtual void invert(int& ierr);
00341    // Invert a Matrix. Matrix must be square.
00342    // N.B. the contents of the matrix are replaced by the inverse.
00343    // Returns ierr = 0 (zero) when successful, otherwise non-zero. 
00344    // This method has less overhead then inverse().
00345 
00346    double determinant() const;
00347    // calculate the determinant of the matrix.
00348 
00349    double trace() const;
00350    // calculate the trace of the matrix (sum of diagonal elements).
00351 
00352    class HepMatrix_row {
00353    public:
00354       inline HepMatrix_row(HepMatrix&,int);
00355       double & operator[](int);
00356    private:
00357       HepMatrix& _a;
00358       int _r;
00359    };
00360    class HepMatrix_row_const {
00361    public:
00362       inline HepMatrix_row_const (const HepMatrix&,int);
00363       const double & operator[](int) const;
00364    private:
00365       const HepMatrix& _a;
00366       int _r;
00367    };
00368    // helper classes for implementing m[i][j]
00369 
00370    inline HepMatrix_row operator[] (int);
00371    inline const HepMatrix_row_const operator[] (int) const;
00372    // Read or write a matrix element.
00373    // While it may not look like it, you simply do m[i][j] to get an
00374    // element. 
00375    // ** Note that the indexing starts from [0][0]. **
00376 
00377 protected:
00378    virtual inline int num_size() const;
00379    virtual void invertHaywood4(int& ierr);
00380    virtual void invertHaywood5(int& ierr);
00381    virtual void invertHaywood6(int& ierr);
00382 
00383 private:
00384    friend class HepMatrix_row;
00385    friend class HepMatrix_row_const;
00386    friend class HepVector;
00387    friend class HepSymMatrix;
00388    friend class HepDiagMatrix;
00389    // Friend classes.
00390 
00391    friend HepMatrix operator+(const HepMatrix &m1, const HepMatrix &m2);
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 HepSymMatrix &m2);
00395    friend HepMatrix operator*(const HepMatrix &m1, const HepDiagMatrix &m2);
00396    friend HepMatrix operator*(const HepSymMatrix &m1, const HepMatrix &m2);
00397    friend HepMatrix operator*(const HepDiagMatrix &m1, const HepMatrix &m2);
00398    friend HepMatrix operator*(const HepVector &m1, const HepMatrix &m2);
00399    friend HepVector operator*(const HepMatrix &m1, const HepVector &m2);
00400    friend HepMatrix operator*(const HepSymMatrix &m1, const HepSymMatrix &m2);
00401    // Multiply a Matrix by a Matrix or Vector.
00402 
00403    friend HepVector solve(const HepMatrix &, const HepVector &);
00404    // solve the system of linear eq
00405    friend HepVector qr_solve(HepMatrix *, const HepVector &);
00406    friend HepMatrix qr_solve(HepMatrix *, const HepMatrix &b);
00407    friend void tridiagonal(HepSymMatrix *a,HepMatrix *hsm);
00408    friend void row_house(HepMatrix *,const HepMatrix &, double,
00409                          int, int, int, int);
00410    friend void row_house(HepMatrix *,const HepVector &, double,
00411                          int, int);
00412    friend void back_solve(const HepMatrix &R, HepVector *b);
00413    friend void back_solve(const HepMatrix &R, HepMatrix *b);
00414    friend void col_givens(HepMatrix *A, double c,
00415                           double s, int k1, int k2, 
00416                           int rowmin, int rowmax);
00417    //    Does a column Givens update.
00418    friend void row_givens(HepMatrix *A, double c,
00419                           double s, int k1, int k2, 
00420                           int colmin, int colmax);
00421    friend void col_house(HepMatrix *,const HepMatrix &, double,
00422                          int, int, int, int);
00423    friend HepVector house(const HepMatrix &a,int row,int col);
00424    friend void house_with_update(HepMatrix *a,int row,int col);
00425    friend void house_with_update(HepMatrix *a,HepMatrix *v,int row,int col);
00426    friend void house_with_update2(HepSymMatrix *a,HepMatrix *v,
00427                                   int row,int col); 
00428 
00429    int dfact_matrix(double &det, int *ir);
00430    // factorize the matrix. If successful, the return code is 0. On
00431    // return, det is the determinant and ir[] is row-interchange
00432    // matrix. See CERNLIB's DFACT routine.
00433 
00434    int dfinv_matrix(int *ir);
00435    // invert the matrix. See CERNLIB DFINV.
00436 
00437    double *m;
00438    int nrow, ncol;
00439    int size;
00440 };
00441 
00442 // Operations other than member functions for Matrix
00443 // implemented in Matrix.cc and Matrix.icc (inline).
00444 
00445 HepMatrix operator*(const HepMatrix &m1, const HepMatrix &m2);
00446 HepMatrix operator*(double t, const HepMatrix &m1);
00447 HepMatrix operator*(const HepMatrix &m1, double t);
00448 // Multiplication operators
00449 // Note that m *= m1 is always faster than m = m * m1.
00450 
00451 HepMatrix operator/(const HepMatrix &m1, double t);
00452 // m = m1 / t. (m /= t is faster if you can use it.)
00453 
00454 HepMatrix operator+(const HepMatrix &m1, const HepMatrix &m2);
00455 // m = m1 + m2;
00456 // Note that m += m1 is always faster than m = m + m1.
00457 
00458 HepMatrix operator-(const HepMatrix &m1, const HepMatrix &m2);
00459 // m = m1 - m2;
00460 // Note that m -= m1 is always faster than m = m - m1.
00461 
00462 HepMatrix dsum(const HepMatrix&, const HepMatrix&);
00463 // Direct sum of two matrices. The direct sum of A and B is the matrix 
00464 //        A 0
00465 //        0 B
00466 
00467 HepVector solve(const HepMatrix &, const HepVector &);
00468 // solve the system of linear equations using LU decomposition.
00469 
00470 HepStd::ostream& operator<<(HepStd::ostream &s, const HepMatrix &q);
00471 // Read in, write out Matrix into a stream.
00472  
00473 //
00474 // Specialized linear algebra functions
00475 //
00476 
00477 HepVector qr_solve(const HepMatrix &A, const HepVector &b);
00478 HepVector qr_solve(HepMatrix *A, const HepVector &b);
00479 HepMatrix qr_solve(const HepMatrix &A, const HepMatrix &b);
00480 HepMatrix qr_solve(HepMatrix *A, const HepMatrix &b);
00481 // Works like backsolve, except matrix does not need to be upper
00482 // triangular. For nonsquare matrix, it solves in the least square sense.
00483 
00484 HepMatrix qr_inverse(const HepMatrix &A);
00485 HepMatrix qr_inverse(HepMatrix *A);
00486 // Finds the inverse of a matrix using QR decomposition.  Note, often what
00487 // you really want is solve or backsolve, they can be much quicker than
00488 // inverse in many calculations.
00489 
00490 
00491 void qr_decomp(HepMatrix *A, HepMatrix *hsm);
00492 HepMatrix qr_decomp(HepMatrix *A);
00493 // Does a QR decomposition of a matrix.
00494 
00495 void back_solve(const HepMatrix &R, HepVector *b);
00496 void back_solve(const HepMatrix &R, HepMatrix *b);
00497 // Solves R*x = b where R is upper triangular.  Also has a variation that
00498 // solves a number of equations of this form in one step, where b is a matrix
00499 // with each column a different vector. See also solve.
00500 
00501 void col_house(HepMatrix *a, const HepMatrix &v, double vnormsq,
00502                int row, int col, int row_start, int col_start);
00503 void col_house(HepMatrix *a, const HepMatrix &v, int row, int col,
00504                int row_start, int col_start);
00505 // Does a column Householder update.
00506 
00507 void col_givens(HepMatrix *A, double c, double s,
00508                 int k1, int k2, int row_min=1, int row_max=0);
00509 // do a column Givens update
00510 
00511 void row_givens(HepMatrix *A, double c, double s,
00512                 int k1, int k2, int col_min=1, int col_max=0);
00513 // do a row Givens update
00514 
00515 void givens(double a, double b, double *c, double *s);
00516 // algorithm 5.1.5 in Golub and Van Loan
00517 
00518 HepVector house(const HepMatrix &a, int row=1, int col=1);
00519 // Returns a Householder vector to zero elements.
00520 
00521 void house_with_update(HepMatrix *a, int row=1, int col=1);
00522 void house_with_update(HepMatrix *a, HepMatrix *v, int row=1, int col=1);
00523 // Finds and does Householder reflection on matrix.
00524 
00525 void row_house(HepMatrix *a, const HepVector &v, double vnormsq,
00526                int row=1, int col=1);
00527 void row_house(HepMatrix *a, const HepMatrix &v, double vnormsq,
00528                int row, int col, int row_start, int col_start);
00529 void row_house(HepMatrix *a, const HepMatrix &v, int row, int col,
00530                int row_start, int col_start);
00531 // Does a row Householder update.
00532 
00533 
00534 #ifdef HEP_NO_INLINE_IN_DECLARATION
00535 #undef inline
00536 #endif
00537 
00538 #ifdef HEP_SHORT_NAMES
00539 typedef HepMatrix Matrix;
00540 #endif
00541 
00542 #ifndef HEP_DEBUG_INLINE
00543 #include "CLHEP/Matrix/Matrix.icc"
00544 #endif
00545 
00546 #endif /*_Matrix_H*/

Class Library for High Energy Physics (version 1.8)