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

Matrix/CLHEP/Matrix/SymMatrix.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 HepSymMatrix class.
00009 // 
00010 // This software written by Nobu Katayama and Mike Smyth, Cornell University.
00011 //
00012 // .SS Usage
00013 //
00014 //   This is very much like the Matrix, except of course it is restricted to
00015 //   Symmetric Matrix.  All the operations for Matrix can also be done here
00016 //   (except for the +=,-=,*= that don't yield a symmetric matrix.  e.g.
00017 //    +=(const Matrix &) is not defined)
00018    
00019 //   The matrix is stored as a lower triangular matrix.
00020 //   In addition to the (row, col) method of finding element, fast(row, col)
00021 //   returns an element with checking to see if row and col need to be 
00022 //   interchanged so that row >= col.
00023 
00024 //   New operations are:
00025 //
00026 // .ft B
00027 //  sym = s.similarity(m);
00028 //
00029 //  This returns m*s*m.T(). This is a similarity
00030 //  transform when m is orthogonal, but nothing
00031 //  restricts m to be orthogonal.  It is just
00032 //  a more efficient way to calculate m*s*m.T,
00033 //  and it realizes that this should be a 
00034 //  HepSymMatrix (the explicit operation m*s*m.T
00035 //  will return a Matrix, not realizing that 
00036 //  it is symmetric).
00037 //
00038 // .ft B
00039 //  sym =  similarityT(m);
00040 //
00041 // This returns m.T()*s*m.
00042 //
00043 // .ft B
00044 // s << m;
00045 //
00046 // This takes the matrix m, and treats it
00047 // as symmetric matrix that is copied to s.
00048 // This is useful for operations that yield
00049 // symmetric matrix, but which the computer
00050 // is too dumb to realize.
00051 //
00052 // .ft B
00053 // s = vT_times_v(const HepVector v);
00054 //
00055 //  calculates v.T()*v.
00056 //
00057 // ./"This code has been written by Mike Smyth, and the algorithms used are
00058 // ./"described in the thesis "A Tracking Library for a Silicon Vertex Detector"
00059 // ./"(Mike Smyth, Cornell University, June 1993).
00060 // ./"This is file contains C++ stuff for doing things with Matrixes.
00061 // ./"To turn on bound checking, define MATRIX_BOUND_CHECK before including
00062 // ./"this file.
00063 //
00064 
00065 #ifndef _SYMMatrix_H_
00066 #define _SYMMatrix_H_
00067 
00068 #ifdef GNUPRAGMA
00069 #pragma interface
00070 #endif
00071 
00072 #include <vector>
00073 
00074 #include "CLHEP/Matrix/defs.h"
00075 #include "CLHEP/Matrix/GenMatrix.h"
00076 
00077 namespace CLHEP {
00078 
00079 class HepRandom;
00080 
00081 class HepMatrix;
00082 class HepDiagMatrix;
00083 class HepVector;
00084 
00089 class HepSymMatrix : public HepGenMatrix {
00090 public:
00091    inline HepSymMatrix();
00092    // Default constructor. Gives 0x0 symmetric matrix.
00093    // Another SymMatrix can be assigned to it.
00094 
00095    explicit HepSymMatrix(int p);
00096    HepSymMatrix(int p, int);
00097    // Constructor. Gives p x p symmetric matrix.
00098    // With a second argument, the matrix is initialized. 0 means a zero
00099    // matrix, 1 means the identity matrix.
00100 
00101    HepSymMatrix(int p, HepRandom &r);
00102 
00103    HepSymMatrix(const HepSymMatrix &m1);
00104    // Copy constructor.
00105 
00106    HepSymMatrix(const HepDiagMatrix &m1);
00107    // Constructor from DiagMatrix
00108 
00109    virtual ~HepSymMatrix();
00110    // Destructor.
00111 
00112    inline int num_row() const;
00113    inline int num_col() const;
00114    // Returns number of rows/columns.
00115 
00116    const double & operator()(int row, int col) const; 
00117    double & operator()(int row, int col);
00118    // Read and write a SymMatrix element.
00119    // ** Note that indexing starts from (1,1). **
00120 
00121    const double & fast(int row, int col) const;
00122    double & fast(int row, int col);
00123    // fast element access.
00124    // Must be row>=col;
00125    // ** Note that indexing starts from (1,1). **
00126 
00127    void assign(const HepMatrix &m2);
00128    // Assigns m2 to s, assuming m2 is a symmetric matrix.
00129 
00130    void assign(const HepSymMatrix &m2);
00131    // Another form of assignment. For consistency.
00132 
00133    HepSymMatrix & operator*=(double t);
00134    // Multiply a SymMatrix by a floating number.
00135 
00136    HepSymMatrix & operator/=(double t); 
00137    // Divide a SymMatrix by a floating number.
00138 
00139    HepSymMatrix & operator+=( const HepSymMatrix &m2);
00140    HepSymMatrix & operator+=( const HepDiagMatrix &m2);
00141    HepSymMatrix & operator-=( const HepSymMatrix &m2);
00142    HepSymMatrix & operator-=( const HepDiagMatrix &m2);
00143    // Add or subtract a SymMatrix.
00144 
00145    HepSymMatrix & operator=( const HepSymMatrix &m2);
00146    HepSymMatrix & operator=( const HepDiagMatrix &m2);
00147    // Assignment operators. Notice that there is no SymMatrix = Matrix.
00148 
00149    HepSymMatrix operator- () const;
00150    // unary minus, ie. flip the sign of each element.
00151 
00152    HepSymMatrix T() const;
00153    // Returns the transpose of a SymMatrix (which is itself).
00154 
00155    HepSymMatrix apply(double (*f)(double, int, int)) const;
00156    // Apply a function to all elements of the matrix.
00157 
00158    HepSymMatrix similarity(const HepMatrix &m1) const;
00159    HepSymMatrix similarity(const HepSymMatrix &m1) const;
00160    // Returns m1*s*m1.T().
00161 
00162    HepSymMatrix similarityT(const HepMatrix &m1) const;
00163    // temporary. test of new similarity.
00164    // Returns m1.T()*s*m1.
00165 
00166    double similarity(const HepVector &v) const;
00167    // Returns v.T()*s*v (This is a scaler).
00168 
00169    HepSymMatrix sub(int min_row, int max_row) const;
00170    // Returns a sub matrix of a SymMatrix.
00171    void sub(int row, const HepSymMatrix &m1);
00172    // Sub matrix of this SymMatrix is replaced with m1.
00173    HepSymMatrix sub(int min_row, int max_row);
00174    // SGI CC bug. I have to have both with/without const. I should not need
00175    // one without const.
00176 
00177    inline HepSymMatrix inverse(int &ifail) const;
00178    // Invert a Matrix. The matrix is not changed
00179    // Returns 0 when successful, otherwise non-zero.
00180 
00181    void invert(int &ifail);
00182    // Invert a Matrix.
00183    // N.B. the contents of the matrix are replaced by the inverse.
00184    // Returns ierr = 0 when successful, otherwise non-zero. 
00185    // This method has less overhead then inverse().
00186 
00187    double determinant() const;
00188    // calculate the determinant of the matrix.
00189 
00190    double trace() const;
00191    // calculate the trace of the matrix (sum of diagonal elements).
00192 
00193    class HepSymMatrix_row {
00194    public:
00195       inline HepSymMatrix_row(HepSymMatrix&,int);
00196       inline double & operator[](int);
00197    private:
00198       HepSymMatrix& _a;
00199       int _r;
00200    };
00201    class HepSymMatrix_row_const {
00202    public:
00203       inline HepSymMatrix_row_const(const HepSymMatrix&,int);
00204       inline const double & operator[](int) const;
00205    private:
00206       const HepSymMatrix& _a;
00207       int _r;
00208    };
00209    // helper class to implement m[i][j]
00210 
00211    inline HepSymMatrix_row operator[] (int);
00212    inline HepSymMatrix_row_const operator[] (int) const;
00213    // Read or write a matrix element.
00214    // While it may not look like it, you simply do m[i][j] to get an
00215    // element. 
00216    // ** Note that the indexing starts from [0][0]. **
00217 
00218    // Special-case inversions for 5x5 and 6x6 symmetric positive definite:
00219    // These set ifail=0 and invert if the matrix was positive definite;
00220    // otherwise ifail=1 and the matrix is left unaltered.
00221    void invertCholesky5 (int &ifail);  
00222    void invertCholesky6 (int &ifail);
00223 
00224    // Inversions for 5x5 and 6x6 forcing use of specific methods:  The
00225    // behavior (though not the speed) will be identical to invert(ifail).
00226    void invertHaywood4 (int & ifail);  
00227    void invertHaywood5 (int &ifail);  
00228    void invertHaywood6 (int &ifail);
00229    void invertBunchKaufman (int &ifail);  
00230 
00231 protected:
00232    inline int num_size() const;
00233   
00234 private:
00235    friend class HepSymMatrix_row;
00236    friend class HepSymMatrix_row_const;
00237    friend class HepMatrix;
00238    friend class HepDiagMatrix;
00239 
00240    friend void tridiagonal(HepSymMatrix *a,HepMatrix *hsm);
00241    friend double condition(const HepSymMatrix &m);
00242    friend void diag_step(HepSymMatrix *t,int begin,int end);
00243    friend void diag_step(HepSymMatrix *t,HepMatrix *u,int begin,int end);
00244    friend HepMatrix diagonalize(HepSymMatrix *s);
00245    friend HepVector house(const HepSymMatrix &a,int row,int col);
00246    friend void house_with_update2(HepSymMatrix *a,HepMatrix *v,int row,int col);
00247 
00248    friend HepSymMatrix operator+(const HepSymMatrix &m1, 
00249                                   const HepSymMatrix &m2);
00250    friend HepSymMatrix operator-(const HepSymMatrix &m1, 
00251                                   const HepSymMatrix &m2);
00252    friend HepMatrix operator*(const HepSymMatrix &m1, const HepSymMatrix &m2);
00253    friend HepMatrix operator*(const HepSymMatrix &m1, const HepMatrix &m2);
00254    friend HepMatrix operator*(const HepMatrix &m1, const HepSymMatrix &m2);
00255    friend HepVector operator*(const HepSymMatrix &m1, const HepVector &m2);
00256    // Multiply a Matrix by a Matrix or Vector.
00257    
00258    friend HepSymMatrix vT_times_v(const HepVector &v);
00259    // Returns v * v.T();
00260 
00261 #ifdef DISABLE_ALLOC
00262    std::vector<double > m;
00263 #else
00264    std::vector<double,Alloc<double,25> > m;
00265 #endif
00266    int nrow;
00267    int size_;                                // total number of elements
00268 
00269    static double posDefFraction5x5;
00270    static double adjustment5x5;
00271    static const  double CHOLESKY_THRESHOLD_5x5;
00272    static const  double CHOLESKY_CREEP_5x5;
00273 
00274    static double posDefFraction6x6;
00275    static double adjustment6x6;
00276    static const double CHOLESKY_THRESHOLD_6x6;
00277    static const double CHOLESKY_CREEP_6x6;
00278 
00279    void invert4  (int & ifail);
00280    void invert5  (int & ifail);
00281    void invert6  (int & ifail);
00282 
00283 };
00284 
00285 //
00286 // Operations other than member functions for Matrix, SymMatrix, DiagMatrix
00287 // and Vectors implemented in Matrix.cc and Matrix.icc (inline).
00288 //
00289 
00290 std::ostream& operator<<(std::ostream &s, const HepSymMatrix &q);
00291 // Write out Matrix, SymMatrix, DiagMatrix and Vector into ostream.
00292 
00293 HepMatrix operator*(const HepMatrix &m1, const HepSymMatrix &m2);
00294 HepMatrix operator*(const HepSymMatrix &m1, const HepMatrix &m2);
00295 HepMatrix operator*(const HepSymMatrix &m1, const HepSymMatrix &m2);
00296 HepSymMatrix operator*(double t, const HepSymMatrix &s1);
00297 HepSymMatrix operator*(const HepSymMatrix &s1, double t);
00298 // Multiplication operators.
00299 // Note that m *= m1 is always faster than m = m * m1
00300 
00301 HepSymMatrix operator/(const HepSymMatrix &m1, double t);
00302 // s = s1 / t. (s /= t is faster if you can use it.)
00303 
00304 HepMatrix operator+(const HepMatrix &m1, const HepSymMatrix &s2);
00305 HepMatrix operator+(const HepSymMatrix &s1, const HepMatrix &m2);
00306 HepSymMatrix operator+(const HepSymMatrix &s1, const HepSymMatrix &s2);
00307 // Addition operators
00308 
00309 HepMatrix operator-(const HepMatrix &m1, const HepSymMatrix &s2);
00310 HepMatrix operator-(const HepSymMatrix &m1, const HepMatrix &m2);
00311 HepSymMatrix operator-(const HepSymMatrix &s1, const HepSymMatrix &s2);
00312 // subtraction operators
00313 
00314 HepSymMatrix dsum(const HepSymMatrix &s1, const HepSymMatrix &s2);
00315 // Direct sum of two symmetric matrices;
00316 
00317 double condition(const HepSymMatrix &m);
00318 // Find the conditon number of a symmetric matrix.
00319 
00320 void diag_step(HepSymMatrix *t, int begin, int end);
00321 void diag_step(HepSymMatrix *t, HepMatrix *u, int begin, int end);
00322 // Implicit symmetric QR step with Wilkinson Shift
00323 
00324 HepMatrix diagonalize(HepSymMatrix *s);
00325 // Diagonalize a symmetric matrix.
00326 // It returns the matrix U so that s_old = U * s_diag * U.T()
00327 
00328 HepVector house(const HepSymMatrix &a, int row=1, int col=1);
00329 void house_with_update2(HepSymMatrix *a, HepMatrix *v, int row=1, int col=1);
00330 // Finds and does Householder reflection on matrix.
00331 
00332 void tridiagonal(HepSymMatrix *a, HepMatrix *hsm);
00333 HepMatrix tridiagonal(HepSymMatrix *a);
00334 // Does a Householder tridiagonalization of a symmetric matrix.
00335 
00336 }  // namespace CLHEP
00337 
00338 #ifdef ENABLE_BACKWARDS_COMPATIBILITY
00339 //  backwards compatibility will be enabled ONLY in CLHEP 1.9
00340 using namespace CLHEP;
00341 #endif
00342 
00343 #ifndef HEP_DEBUG_INLINE
00344 #include "CLHEP/Matrix/SymMatrix.icc"
00345 #endif
00346 
00347 #endif