include/Geant4/G4ErrorSymMatrix.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 // Class Description:
0028 //
0029 // Simplified version of CLHEP HepSymMatrix class.
0030
0031 // History:
0032 // - Imported from CLHEP and modified:   P. Arce    May 2007
0033 // --------------------------------------------------------------------
0034
0035 #ifndef G4ErrorSymMatrix_hh
0036 #define G4ErrorSymMatrix_hh
0037
0038 #include <vector>
0039 #include "globals.hh"
0040
0041 class G4ErrorMatrix;
0042
0043 class G4ErrorSymMatrix
0044 {
0045  public:  // with description
0046   inline G4ErrorSymMatrix();
0047   // Default constructor. Gives 0x0 symmetric matrix.
0048   // Another G4ErrorSymMatrix can be assigned to it.
0049
0050   explicit G4ErrorSymMatrix(G4int p);
0051   G4ErrorSymMatrix(G4int p, G4int);
0052   // Constructor. Gives p x p symmetric matrix.
0053   // With a second argument, the matrix is initialized. 0 means a zero
0054   // matrix, 1 means the identity matrix.
0055
0056   G4ErrorSymMatrix(const G4ErrorSymMatrix& m1);
0057   G4ErrorSymMatrix(G4ErrorSymMatrix&&) = default;
0058   // Copy and move constructor.
0059
0060   // Constructor from DiagMatrix
0061
0062   virtual ~G4ErrorSymMatrix();
0063   // Destructor.
0064
0065   inline G4int num_row() const;
0066   inline G4int num_col() const;
0067   // Returns number of rows/columns.
0068
0069   const G4double& operator()(G4int row, G4int col) const;
0070   G4double& operator()(G4int row, G4int col);
0071   // Read and write a G4ErrorSymMatrix element.
0072   // ** Note that indexing starts from (1,1). **
0073
0074   const G4double& fast(G4int row, G4int col) const;
0075   G4double& fast(G4int row, G4int col);
0076   // fast element access.
0077   // Must be row>=col;
0078   // ** Note that indexing starts from (1,1). **
0079
0080   void assign(const G4ErrorMatrix& m2);
0081   // Assigns m2 to s, assuming m2 is a symmetric matrix.
0082
0083   void assign(const G4ErrorSymMatrix& m2);
0084   // Another form of assignment. For consistency.
0085
0086   G4ErrorSymMatrix& operator*=(G4double t);
0087   // Multiply a G4ErrorSymMatrix by a floating number.
0088
0089   G4ErrorSymMatrix& operator/=(G4double t);
0090   // Divide a G4ErrorSymMatrix by a floating number.
0091
0092   G4ErrorSymMatrix& operator+=(const G4ErrorSymMatrix& m2);
0093   G4ErrorSymMatrix& operator-=(const G4ErrorSymMatrix& m2);
0094   // Add or subtract a G4ErrorSymMatrix.
0095
0096   G4ErrorSymMatrix& operator=(const G4ErrorSymMatrix& m2);
0097   G4ErrorSymMatrix& operator=(G4ErrorSymMatrix&&) = default;
0098   // Assignment and move operators.
0099   // Notice that there is no G4ErrorSymMatrix = Matrix.
0100
0101   G4ErrorSymMatrix operator-() const;
0102   // unary minus, ie. flip the sign of each element.
0103
0104   G4ErrorSymMatrix T() const;
0105   // Returns the transpose of a G4ErrorSymMatrix (which is itself).
0106
0107   G4ErrorSymMatrix apply(G4double (*f)(G4double, G4int, G4int)) const;
0108   // Apply a function to all elements of the matrix.
0109
0110   G4ErrorSymMatrix similarity(const G4ErrorMatrix& m1) const;
0111   G4ErrorSymMatrix similarity(const G4ErrorSymMatrix& m1) const;
0112   // Returns m1*s*m1.T().
0113
0114   G4ErrorSymMatrix similarityT(const G4ErrorMatrix& m1) const;
0115   // temporary. test of new similarity.
0116   // Returns m1.T()*s*m1.
0117
0118   G4ErrorSymMatrix sub(G4int min_row, G4int max_row) const;
0119   // Returns a sub matrix of a G4ErrorSymMatrix.
0120
0121   void sub(G4int row, const G4ErrorSymMatrix& m1);
0122   // Sub matrix of this G4ErrorSymMatrix is replaced with m1.
0123
0124   G4ErrorSymMatrix sub(G4int min_row, G4int max_row);
0125   // SGI CC bug. I have to have both with/without const. I should not need
0126   // one without const.
0127
0128   inline G4ErrorSymMatrix inverse(G4int& ifail) const;
0129   // Invert a Matrix. The matrix is not changed
0130   // Returns 0 when successful, otherwise non-zero.
0131
0132   void invert(G4int& ifail);
0133   // Invert a Matrix.
0134   // N.B. the contents of the matrix are replaced by the inverse.
0135   // Returns ierr = 0 when successful, otherwise non-zero.
0136   // This method has less overhead then inverse().
0137
0138   G4double determinant() const;
0139   // calculate the determinant of the matrix.
0140
0141   G4double trace() const;
0142   // calculate the trace of the matrix (sum of diagonal elements).
0143
0144   class G4ErrorSymMatrix_row
0145   {
0146    public:
0147     inline G4ErrorSymMatrix_row(G4ErrorSymMatrix&, G4int);
0148     inline G4double& operator[](G4int);
0149
0150    private:
0151     G4ErrorSymMatrix& _a;
0152     G4int _r;
0153   };
0154   class G4ErrorSymMatrix_row_const
0155   {
0156    public:
0157     inline G4ErrorSymMatrix_row_const(const G4ErrorSymMatrix&, G4int);
0158     inline const G4double& operator[](G4int) const;
0159
0160    private:
0161     const G4ErrorSymMatrix& _a;
0162     G4int _r;
0163   };
0164   // helper class to implement m[i][j]
0165
0166   inline G4ErrorSymMatrix_row operator[](G4int);
0167   inline G4ErrorSymMatrix_row_const operator[](G4int) const;
0168   // Read or write a matrix element.
0169   // While it may not look like it, you simply do m[i][j] to get an
0170   // element.
0171   // ** Note that the indexing starts from [0][0]. **
0172
0173   // Special-case inversions for 5x5 and 6x6 symmetric positive definite:
0174   // These set ifail=0 and invert if the matrix was positive definite;
0175   // otherwise ifail=1 and the matrix is left unaltered.
0176
0177   void invertCholesky5(G4int& ifail);
0178   void invertCholesky6(G4int& ifail);
0179
0180   // Inversions for 5x5 and 6x6 forcing use of specific methods:  The
0181   // behavior (though not the speed) will be identical to invert(ifail).
0182
0183   void invertHaywood4(G4int& ifail);
0184   void invertHaywood5(G4int& ifail);
0185   void invertHaywood6(G4int& ifail);
0186   void invertBunchKaufman(G4int& ifail);
0187
0188  protected:
0189   inline G4int num_size() const;
0190
0191  private:
0192   friend class G4ErrorSymMatrix_row;
0193   friend class G4ErrorSymMatrix_row_const;
0194   friend class G4ErrorMatrix;
0195
0196   friend void tridiagonal(G4ErrorSymMatrix* a, G4ErrorMatrix* hsm);
0197   friend G4double condition(const G4ErrorSymMatrix& m);
0198   friend void diag_step(G4ErrorSymMatrix* t, G4int begin, G4int end);
0199   friend void diag_step(G4ErrorSymMatrix* t, G4ErrorMatrix* u, G4int begin,
0200                         G4int end);
0201   friend G4ErrorMatrix diagonalize(G4ErrorSymMatrix* s);
0202   friend void house_with_update2(G4ErrorSymMatrix* a, G4ErrorMatrix* v,
0203                                  G4int row, G4int col);
0204
0205   friend G4ErrorSymMatrix operator+(const G4ErrorSymMatrix& m1,
0206                                     const G4ErrorSymMatrix& m2);
0207   friend G4ErrorSymMatrix operator-(const G4ErrorSymMatrix& m1,
0208                                     const G4ErrorSymMatrix& m2);
0209   friend G4ErrorMatrix operator*(const G4ErrorSymMatrix& m1,
0210                                  const G4ErrorSymMatrix& m2);
0211   friend G4ErrorMatrix operator*(const G4ErrorSymMatrix& m1,
0212                                  const G4ErrorMatrix& m2);
0213   friend G4ErrorMatrix operator*(const G4ErrorMatrix& m1,
0214                                  const G4ErrorSymMatrix& m2);
0215   // Multiply a Matrix by a Matrix or Vector.
0216
0217   // Returns v * v.T();
0218
0219   std::vector<G4double> m;
0220
0221   G4int nrow;
0222   G4int size;  // total number of elements
0223
0224   static G4ThreadLocal G4double posDefFraction5x5;
0225   static G4ThreadLocal G4double adjustment5x5;
0226   static const G4double CHOLESKY_THRESHOLD_5x5;
0227   static const G4double CHOLESKY_CREEP_5x5;
0228
0229   static G4ThreadLocal G4double posDefFraction6x6;
0230   static G4ThreadLocal G4double adjustment6x6;
0231   static const G4double CHOLESKY_THRESHOLD_6x6;
0232   static const G4double CHOLESKY_CREEP_6x6;
0233
0234   void invert4(G4int& ifail);
0235   void invert5(G4int& ifail);
0236   void invert6(G4int& ifail);
0237 };
0238
0239 //
0240 // Operations other than member functions for Matrix, G4ErrorSymMatrix,
0241 // DiagMatrix and Vectors
0242 //
0243
0244 std::ostream& operator<<(std::ostream& s, const G4ErrorSymMatrix& q);
0245 // Write out Matrix, G4ErrorSymMatrix, DiagMatrix and Vector into ostream.
0246
0247 G4ErrorMatrix operator*(const G4ErrorMatrix& m1, const G4ErrorSymMatrix& m2);
0248 G4ErrorMatrix operator*(const G4ErrorSymMatrix& m1, const G4ErrorMatrix& m2);
0249 G4ErrorMatrix operator*(const G4ErrorSymMatrix& m1, const G4ErrorSymMatrix& m2);
0250 G4ErrorSymMatrix operator*(G4double t, const G4ErrorSymMatrix& s1);
0251 G4ErrorSymMatrix operator*(const G4ErrorSymMatrix& s1, G4double t);
0252 // Multiplication operators.
0253 // Note that m *= m1 is always faster than m = m * m1
0254
0255 G4ErrorSymMatrix operator/(const G4ErrorSymMatrix& m1, G4double t);
0256 // s = s1 / t. (s /= t is faster if you can use it.)
0257
0258 G4ErrorMatrix operator+(const G4ErrorMatrix& m1, const G4ErrorSymMatrix& s2);
0259 G4ErrorMatrix operator+(const G4ErrorSymMatrix& s1, const G4ErrorMatrix& m2);
0260 G4ErrorSymMatrix operator+(const G4ErrorSymMatrix& s1,
0261                            const G4ErrorSymMatrix& s2);
0262 // Addition operators
0263
0264 G4ErrorMatrix operator-(const G4ErrorMatrix& m1, const G4ErrorSymMatrix& s2);
0265 G4ErrorMatrix operator-(const G4ErrorSymMatrix& m1, const G4ErrorMatrix& m2);
0266 G4ErrorSymMatrix operator-(const G4ErrorSymMatrix& s1,
0267                            const G4ErrorSymMatrix& s2);
0268 // subtraction operators
0269
0270 G4ErrorSymMatrix dsum(const G4ErrorSymMatrix& s1, const G4ErrorSymMatrix& s2);
0271 // Direct sum of two symmetric matrices;
0272
0273 G4double condition(const G4ErrorSymMatrix& m);
0274 // Find the conditon number of a symmetric matrix.
0275
0276 void diag_step(G4ErrorSymMatrix* t, G4int begin, G4int end);
0277 void diag_step(G4ErrorSymMatrix* t, G4ErrorMatrix* u, G4int begin, G4int end);
0278 // Implicit symmetric QR step with Wilkinson Shift
0279
0280 G4ErrorMatrix diagonalize(G4ErrorSymMatrix* s);
0281 // Diagonalize a symmetric matrix.
0282 // It returns the matrix U so that s_old = U * s_diag * U.T()
0283
0284 void house_with_update2(G4ErrorSymMatrix* a, G4ErrorMatrix* v, G4int row = 1,
0285                         G4int col = 1);
0286 // Finds and does Householder reflection on matrix.
0287
0288 void tridiagonal(G4ErrorSymMatrix* a, G4ErrorMatrix* hsm);
0289 G4ErrorMatrix tridiagonal(G4ErrorSymMatrix* a);
0290 // Does a Householder tridiagonalization of a symmetric matrix.
0291
0292 #include "G4ErrorSymMatrix.icc"
0293
0294 #endif