G4ErrorSymMatrix.hh

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id: G4ErrorSymMatrix.hh 69766 2013-05-14 14:33:55Z gcosmo $
00027 //
00028 // Class Description:
00029 //
00030 // Simplified version of CLHEP HepSymMatrix class.
00031 
00032 // History:
00033 // - Imported from CLHEP and modified:   P. Arce    May 2007
00034 // --------------------------------------------------------------------
00035 
00036 #ifndef G4ErrorSymMatrix_hh
00037 #define G4ErrorSymMatrix_hh
00038 
00039 #include <vector>
00040 
00041 class G4ErrorMatrix;
00042 
00043 class G4ErrorSymMatrix
00044 {
00045  public:  //with description
00046 
00047    inline G4ErrorSymMatrix();
00048    // Default constructor. Gives 0x0 symmetric matrix.
00049    // Another G4ErrorSymMatrix can be assigned to it.
00050 
00051    explicit G4ErrorSymMatrix(G4int p);
00052    G4ErrorSymMatrix(G4int p, G4int);
00053    // Constructor. Gives p x p symmetric matrix.
00054    // With a second argument, the matrix is initialized. 0 means a zero
00055    // matrix, 1 means the identity matrix.
00056 
00057    G4ErrorSymMatrix(const G4ErrorSymMatrix &m1);
00058    // Copy constructor.
00059 
00060    // Constructor from DiagMatrix
00061 
00062    virtual ~G4ErrorSymMatrix();
00063    // Destructor.
00064 
00065    inline G4int num_row() const;
00066    inline G4int num_col() const;
00067    // Returns number of rows/columns.
00068 
00069    const G4double & operator()(G4int row, G4int col) const; 
00070    G4double & operator()(G4int row, G4int col);
00071    // Read and write a G4ErrorSymMatrix element.
00072    // ** Note that indexing starts from (1,1). **
00073 
00074    const G4double & fast(G4int row, G4int col) const;
00075    G4double & fast(G4int row, G4int col);
00076    // fast element access.
00077    // Must be row>=col;
00078    // ** Note that indexing starts from (1,1). **
00079 
00080    void assign(const G4ErrorMatrix &m2);
00081    // Assigns m2 to s, assuming m2 is a symmetric matrix.
00082 
00083    void assign(const G4ErrorSymMatrix &m2);
00084    // Another form of assignment. For consistency.
00085 
00086    G4ErrorSymMatrix & operator*=(G4double t);
00087    // Multiply a G4ErrorSymMatrix by a floating number.
00088 
00089    G4ErrorSymMatrix & operator/=(G4double t); 
00090    // Divide a G4ErrorSymMatrix by a floating number.
00091 
00092    G4ErrorSymMatrix & operator+=( const G4ErrorSymMatrix &m2);
00093    G4ErrorSymMatrix & operator-=( const G4ErrorSymMatrix &m2);
00094    // Add or subtract a G4ErrorSymMatrix.
00095 
00096    G4ErrorSymMatrix & operator=( const G4ErrorSymMatrix &m2);
00097    // Assignment operators. Notice that there is no G4ErrorSymMatrix = Matrix.
00098 
00099    G4ErrorSymMatrix operator- () const;
00100    // unary minus, ie. flip the sign of each element.
00101 
00102    G4ErrorSymMatrix T() const;
00103    // Returns the transpose of a G4ErrorSymMatrix (which is itself).
00104 
00105    G4ErrorSymMatrix apply(G4double (*f)(G4double, G4int, G4int)) const;
00106    // Apply a function to all elements of the matrix.
00107 
00108    G4ErrorSymMatrix similarity(const G4ErrorMatrix &m1) const;
00109    G4ErrorSymMatrix similarity(const G4ErrorSymMatrix &m1) const;
00110    // Returns m1*s*m1.T().
00111 
00112    G4ErrorSymMatrix similarityT(const G4ErrorMatrix &m1) const;
00113    // temporary. test of new similarity.
00114    // Returns m1.T()*s*m1.
00115 
00116    G4ErrorSymMatrix sub(G4int min_row, G4int max_row) const;
00117    // Returns a sub matrix of a G4ErrorSymMatrix.
00118 
00119    void sub(G4int row, const G4ErrorSymMatrix &m1);
00120    // Sub matrix of this G4ErrorSymMatrix is replaced with m1.
00121 
00122    G4ErrorSymMatrix sub(G4int min_row, G4int max_row);
00123    // SGI CC bug. I have to have both with/without const. I should not need
00124    // one without const.
00125 
00126    inline G4ErrorSymMatrix inverse(G4int &ifail) const;
00127    // Invert a Matrix. The matrix is not changed
00128    // Returns 0 when successful, otherwise non-zero.
00129 
00130    void invert(G4int &ifail);
00131    // Invert a Matrix.
00132    // N.B. the contents of the matrix are replaced by the inverse.
00133    // Returns ierr = 0 when successful, otherwise non-zero. 
00134    // This method has less overhead then inverse().
00135 
00136    G4double determinant() const;
00137    // calculate the determinant of the matrix.
00138 
00139    G4double trace() const;
00140    // calculate the trace of the matrix (sum of diagonal elements).
00141 
00142    class G4ErrorSymMatrix_row
00143    {
00144      public:
00145       inline G4ErrorSymMatrix_row(G4ErrorSymMatrix&,G4int);
00146       inline G4double & operator[](G4int);
00147      private:
00148       G4ErrorSymMatrix& _a;
00149       G4int _r;
00150    };
00151    class G4ErrorSymMatrix_row_const
00152    {
00153      public:
00154       inline G4ErrorSymMatrix_row_const(const G4ErrorSymMatrix&,G4int);
00155       inline const G4double & operator[](G4int) const;
00156      private:
00157       const G4ErrorSymMatrix& _a;
00158       G4int _r;
00159    };
00160    // helper class to implement m[i][j]
00161 
00162    inline G4ErrorSymMatrix_row operator[] (G4int);
00163    inline G4ErrorSymMatrix_row_const operator[] (G4int) const;
00164    // Read or write a matrix element.
00165    // While it may not look like it, you simply do m[i][j] to get an
00166    // element. 
00167    // ** Note that the indexing starts from [0][0]. **
00168 
00169    // Special-case inversions for 5x5 and 6x6 symmetric positive definite:
00170    // These set ifail=0 and invert if the matrix was positive definite;
00171    // otherwise ifail=1 and the matrix is left unaltered.
00172 
00173    void invertCholesky5 (G4int &ifail);  
00174    void invertCholesky6 (G4int &ifail);
00175 
00176    // Inversions for 5x5 and 6x6 forcing use of specific methods:  The
00177    // behavior (though not the speed) will be identical to invert(ifail).
00178 
00179    void invertHaywood4 (G4int & ifail);  
00180    void invertHaywood5 (G4int &ifail);  
00181    void invertHaywood6 (G4int &ifail);
00182    void invertBunchKaufman (G4int &ifail);  
00183 
00184  protected:
00185 
00186    inline G4int num_size() const;
00187   
00188  private:
00189 
00190    friend class G4ErrorSymMatrix_row;
00191    friend class G4ErrorSymMatrix_row_const;
00192    friend class G4ErrorMatrix;
00193 
00194    friend void tridiagonal(G4ErrorSymMatrix *a, G4ErrorMatrix *hsm);
00195    friend G4double condition(const G4ErrorSymMatrix &m);
00196    friend void diag_step(G4ErrorSymMatrix *t, G4int begin, G4int end);
00197    friend void diag_step(G4ErrorSymMatrix *t, G4ErrorMatrix *u,
00198                          G4int begin, G4int end);
00199    friend G4ErrorMatrix diagonalize(G4ErrorSymMatrix *s);
00200    friend void house_with_update2(G4ErrorSymMatrix *a, G4ErrorMatrix *v,
00201                                   G4int row, G4int col);
00202 
00203    friend G4ErrorSymMatrix operator+(const G4ErrorSymMatrix &m1, 
00204                                      const G4ErrorSymMatrix &m2);
00205    friend G4ErrorSymMatrix operator-(const G4ErrorSymMatrix &m1, 
00206                                      const G4ErrorSymMatrix &m2);
00207    friend G4ErrorMatrix operator*(const G4ErrorSymMatrix &m1,
00208                                   const G4ErrorSymMatrix &m2);
00209    friend G4ErrorMatrix operator*(const G4ErrorSymMatrix &m1,
00210                                   const G4ErrorMatrix &m2);
00211    friend G4ErrorMatrix operator*(const G4ErrorMatrix &m1,
00212                                   const G4ErrorSymMatrix &m2);
00213    // Multiply a Matrix by a Matrix or Vector.
00214    
00215    // Returns v * v.T();
00216 
00217    std::vector<G4double > m;
00218 
00219    G4int nrow;
00220    G4int size;   // total number of elements
00221 
00222    static G4double posDefFraction5x5;
00223    static G4double adjustment5x5;
00224    static const  G4double CHOLESKY_THRESHOLD_5x5;
00225    static const  G4double CHOLESKY_CREEP_5x5;
00226 
00227    static G4double posDefFraction6x6;
00228    static G4double adjustment6x6;
00229    static const G4double CHOLESKY_THRESHOLD_6x6;
00230    static const G4double CHOLESKY_CREEP_6x6;
00231 
00232    void invert4  (G4int & ifail);
00233    void invert5  (G4int & ifail);
00234    void invert6  (G4int & ifail);
00235 };
00236 
00237 //
00238 // Operations other than member functions for Matrix, G4ErrorSymMatrix,
00239 // DiagMatrix and Vectors
00240 //
00241 
00242 std::ostream& operator<<(std::ostream &s, const G4ErrorSymMatrix &q);
00243 // Write out Matrix, G4ErrorSymMatrix, DiagMatrix and Vector into ostream.
00244 
00245 G4ErrorMatrix operator*(const G4ErrorMatrix &m1,
00246                         const G4ErrorSymMatrix &m2);
00247 G4ErrorMatrix operator*(const G4ErrorSymMatrix &m1,
00248                         const G4ErrorMatrix &m2);
00249 G4ErrorMatrix operator*(const G4ErrorSymMatrix &m1,
00250                         const G4ErrorSymMatrix &m2);
00251 G4ErrorSymMatrix operator*(G4double t, const G4ErrorSymMatrix &s1);
00252 G4ErrorSymMatrix operator*(const G4ErrorSymMatrix &s1, G4double t);
00253 // Multiplication operators.
00254 // Note that m *= m1 is always faster than m = m * m1
00255 
00256 G4ErrorSymMatrix operator/(const G4ErrorSymMatrix &m1, G4double t);
00257 // s = s1 / t. (s /= t is faster if you can use it.)
00258 
00259 G4ErrorMatrix operator+(const G4ErrorMatrix &m1,
00260                         const G4ErrorSymMatrix &s2);
00261 G4ErrorMatrix operator+(const G4ErrorSymMatrix &s1,
00262                         const G4ErrorMatrix &m2);
00263 G4ErrorSymMatrix operator+(const G4ErrorSymMatrix &s1,
00264                            const G4ErrorSymMatrix &s2);
00265 // Addition operators
00266 
00267 G4ErrorMatrix operator-(const G4ErrorMatrix &m1,
00268                         const G4ErrorSymMatrix &s2);
00269 G4ErrorMatrix operator-(const G4ErrorSymMatrix &m1,
00270                         const G4ErrorMatrix &m2);
00271 G4ErrorSymMatrix operator-(const G4ErrorSymMatrix &s1,
00272                            const G4ErrorSymMatrix &s2);
00273 // subtraction operators
00274 
00275 G4ErrorSymMatrix dsum(const G4ErrorSymMatrix &s1,
00276                       const G4ErrorSymMatrix &s2);
00277 // Direct sum of two symmetric matrices;
00278 
00279 G4double condition(const G4ErrorSymMatrix &m);
00280 // Find the conditon number of a symmetric matrix.
00281 
00282 void diag_step(G4ErrorSymMatrix *t, G4int begin, G4int end);
00283 void diag_step(G4ErrorSymMatrix *t, G4ErrorMatrix *u, G4int begin, G4int end);
00284 // Implicit symmetric QR step with Wilkinson Shift
00285 
00286 G4ErrorMatrix diagonalize(G4ErrorSymMatrix *s);
00287 // Diagonalize a symmetric matrix.
00288 // It returns the matrix U so that s_old = U * s_diag * U.T()
00289 
00290 void house_with_update2(G4ErrorSymMatrix *a, G4ErrorMatrix *v,
00291                         G4int row=1, G4int col=1);
00292 // Finds and does Householder reflection on matrix.
00293 
00294 void tridiagonal(G4ErrorSymMatrix *a, G4ErrorMatrix *hsm);
00295 G4ErrorMatrix tridiagonal(G4ErrorSymMatrix *a);
00296 // Does a Householder tridiagonalization of a symmetric matrix.
00297 
00298 #include "G4ErrorSymMatrix.icc"
00299 
00300 #endif

Generated on Mon May 27 17:48:12 2013 for Geant4 by  doxygen 1.4.7