CLHEP/Matrix/SymMatrix.h

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

Class Library for High Energy Physics (version 1.8)