include/Geant4/G4ErrorMatrix.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 HepMatrix class
0030
0031 // History:
0032 // - Imported from CLHEP and modified:   P. Arce    May 2007
0033 // --------------------------------------------------------------------
0034
0035 #ifndef G4ErrorMatrix_hh
0036 #define G4ErrorMatrix_hh
0037
0038 #include <vector>
0039 #include "G4Types.hh"
0040
0041 class G4ErrorSymMatrix;
0042
0043 typedef std::vector<G4double>::iterator G4ErrorMatrixIter;
0044 typedef std::vector<G4double>::const_iterator G4ErrorMatrixConstIter;
0045
0046 class G4ErrorMatrix
0047 {
0048  public:  // with description
0049   G4ErrorMatrix();
0050   // Default constructor. Gives 0 x 0 matrix.
0051   // Another G4ErrorMatrix can be assigned to it.
0052
0053   G4ErrorMatrix(G4int p, G4int q);
0054   // Constructor. Gives an unitialized p x q matrix.
0055
0056   G4ErrorMatrix(G4int p, G4int q, G4int i);
0057   // Constructor. Gives an initialized p x q matrix.
0058   // If i=0, it is initialized to all 0. If i=1, the diagonal elements
0059   // are set to 1.0.
0060
0061   G4ErrorMatrix(const G4ErrorMatrix& m1);
0062   G4ErrorMatrix(G4ErrorMatrix&&) = default;
0063   // Copy and move constructor.
0064
0065   G4ErrorMatrix(const G4ErrorSymMatrix& m1);
0066   // Constructors from G4ErrorSymG4ErrorMatrix, DiagG4ErrorMatrix and Vector.
0067
0068   virtual ~G4ErrorMatrix();
0069   // Destructor.
0070
0071   inline virtual G4int num_row() const;
0072   // Returns the number of rows.
0073
0074   inline virtual G4int num_col() const;
0075   // Returns the number of columns.
0076
0077   inline virtual const G4double& operator()(G4int row, G4int col) const;
0078   inline virtual G4double& operator()(G4int row, G4int col);
0079   // Read or write a matrix element.
0080   // ** Note that the indexing starts from (1,1). **
0081
0082   G4ErrorMatrix& operator*=(G4double t);
0083   // Multiply a G4ErrorMatrix by a floating number.
0084
0085   G4ErrorMatrix& operator/=(G4double t);
0086   // Divide a G4ErrorMatrix by a floating number.
0087
0088   G4ErrorMatrix& operator+=(const G4ErrorMatrix& m2);
0089   G4ErrorMatrix& operator+=(const G4ErrorSymMatrix& m2);
0090   G4ErrorMatrix& operator-=(const G4ErrorMatrix& m2);
0091   G4ErrorMatrix& operator-=(const G4ErrorSymMatrix& m2);
0092   // Add or subtract a G4ErrorMatrix.
0093   // When adding/subtracting Vector, G4ErrorMatrix must have num_col of one.
0094
0095   G4ErrorMatrix& operator=(const G4ErrorMatrix& m2);
0096   G4ErrorMatrix& operator=(const G4ErrorSymMatrix& m2);
0097   G4ErrorMatrix& operator=(G4ErrorMatrix&&) = default;
0098   // Assignment and move operators.
0099
0100   G4ErrorMatrix operator-() const;
0101   // unary minus, ie. flip the sign of each element.
0102
0103   G4ErrorMatrix apply(G4double (*f)(G4double, G4int, G4int)) const;
0104   // Apply a function to all elements of the matrix.
0105
0106   G4ErrorMatrix T() const;
0107   // Returns the transpose of a G4ErrorMatrix.
0108
0109   G4ErrorMatrix sub(G4int min_row, G4int max_row, G4int min_col,
0110                     G4int max_col) const;
0111   // Returns a sub matrix of a G4ErrorMatrix.
0112   // WARNING: rows and columns are numbered from 1
0113
0114   void sub(G4int row, G4int col, const G4ErrorMatrix& m1);
0115   // Sub matrix of this G4ErrorMatrix is replaced with m1.
0116   // WARNING: rows and columns are numbered from 1
0117
0118   inline G4ErrorMatrix inverse(G4int& ierr) const;
0119   // Invert a G4ErrorMatrix. G4ErrorMatrix must be square and is not changed.
0120   // Returns ierr = 0 (zero) when successful, otherwise non-zero.
0121
0122   virtual void invert(G4int& ierr);
0123   // Invert a G4ErrorMatrix. G4ErrorMatrix must be square.
0124   // N.B. the contents of the matrix are replaced by the inverse.
0125   // Returns ierr = 0 (zero) when successful, otherwise non-zero.
0126   // This method has less overhead then inverse().
0127
0128   G4double determinant() const;
0129   // calculate the determinant of the matrix.
0130
0131   G4double trace() const;
0132   // calculate the trace of the matrix (sum of diagonal elements).
0133
0134   class G4ErrorMatrix_row
0135   {
0136     typedef std::vector<G4double>::const_iterator G4ErrorMatrixConstIter;
0137
0138    public:
0139     inline G4ErrorMatrix_row(G4ErrorMatrix&, G4int);
0140     G4double& operator[](G4int);
0141
0142    private:
0143     G4ErrorMatrix& _a;
0144     G4int _r;
0145   };
0146   class G4ErrorMatrix_row_const
0147   {
0148    public:
0149     inline G4ErrorMatrix_row_const(const G4ErrorMatrix&, G4int);
0150     const G4double& operator[](G4int) const;
0151
0152    private:
0153     const G4ErrorMatrix& _a;
0154     G4int _r;
0155   };
0156   // helper classes for implementing m[i][j]
0157
0158   inline G4ErrorMatrix_row operator[](G4int);
0159   inline const G4ErrorMatrix_row_const operator[](G4int) const;
0160   // Read or write a matrix element.
0161   // While it may not look like it, you simply do m[i][j] to get an element.
0162   // ** Note that the indexing starts from [0][0]. **
0163
0164  protected:
0165   virtual inline G4int num_size() const;
0166   virtual void invertHaywood4(G4int& ierr);
0167   virtual void invertHaywood5(G4int& ierr);
0168   virtual void invertHaywood6(G4int& ierr);
0169
0170  public:
0171   static void error(const char* s);
0172
0173  private:
0174   friend class G4ErrorMatrix_row;
0175   friend class G4ErrorMatrix_row_const;
0176   friend class G4ErrorSymMatrix;
0177   // Friend classes.
0178
0179   friend G4ErrorMatrix operator+(const G4ErrorMatrix& m1,
0180                                  const G4ErrorMatrix& m2);
0181   friend G4ErrorMatrix operator-(const G4ErrorMatrix& m1,
0182                                  const G4ErrorMatrix& m2);
0183   friend G4ErrorMatrix operator*(const G4ErrorMatrix& m1,
0184                                  const G4ErrorMatrix& m2);
0185   friend G4ErrorMatrix operator*(const G4ErrorMatrix& m1,
0186                                  const G4ErrorSymMatrix& m2);
0187   friend G4ErrorMatrix operator*(const G4ErrorSymMatrix& m1,
0188                                  const G4ErrorMatrix& m2);
0189   friend G4ErrorMatrix operator*(const G4ErrorSymMatrix& m1,
0190                                  const G4ErrorSymMatrix& m2);
0191   // Multiply a G4ErrorMatrix by a G4ErrorMatrix or Vector.
0192
0193   // solve the system of linear eq
0194
0195   friend G4ErrorMatrix qr_solve(G4ErrorMatrix*, const G4ErrorMatrix& b);
0196   friend void tridiagonal(G4ErrorSymMatrix* a, G4ErrorMatrix* hsm);
0197   friend void row_house(G4ErrorMatrix*, const G4ErrorMatrix&, G4double, G4int,
0198                         G4int, G4int, G4int);
0199   friend void back_solve(const G4ErrorMatrix& R, G4ErrorMatrix* b);
0200   friend void col_givens(G4ErrorMatrix* A, G4double c, G4double s, G4int k1,
0201                          G4int k2, G4int rowmin, G4int rowmax);
0202   // Does a column Givens update.
0203
0204   friend void row_givens(G4ErrorMatrix* A, G4double c, G4double s, G4int k1,
0205                          G4int k2, G4int colmin, G4int colmax);
0206   friend void col_house(G4ErrorMatrix*, const G4ErrorMatrix&, G4double, G4int,
0207                         G4int, G4int, G4int);
0208   friend void house_with_update(G4ErrorMatrix* a, G4int row, G4int col);
0209   friend void house_with_update(G4ErrorMatrix* a, G4ErrorMatrix* v, G4int row,
0210                                 G4int col);
0211   friend void house_with_update2(G4ErrorSymMatrix* a, G4ErrorMatrix* v,
0212                                  G4int row, G4int col);
0213
0214   G4int dfact_matrix(G4double& det, G4int* ir);
0215   // factorize the matrix. If successful, the return code is 0. On
0216   // return, det is the determinant and ir[] is row-interchange
0217   // matrix. See CERNLIB's DFACT routine.
0218
0219   G4int dfinv_matrix(G4int* ir);
0220   // invert the matrix. See CERNLIB DFINV.
0221
0222   std::vector<G4double> m;
0223
0224   G4int nrow, ncol;
0225   G4int size;
0226 };
0227
0228 // Operations other than member functions for G4ErrorMatrix
0229
0230 G4ErrorMatrix operator*(const G4ErrorMatrix& m1, const G4ErrorMatrix& m2);
0231 G4ErrorMatrix operator*(G4double t, const G4ErrorMatrix& m1);
0232 G4ErrorMatrix operator*(const G4ErrorMatrix& m1, G4double t);
0233 // Multiplication operators
0234 // Note that m *= m1 is always faster than m = m * m1.
0235
0236 G4ErrorMatrix operator/(const G4ErrorMatrix& m1, G4double t);
0237 // m = m1 / t. (m /= t is faster if you can use it.)
0238
0239 G4ErrorMatrix operator+(const G4ErrorMatrix& m1, const G4ErrorMatrix& m2);
0240 // m = m1 + m2;
0241 // Note that m += m1 is always faster than m = m + m1.
0242
0243 G4ErrorMatrix operator-(const G4ErrorMatrix& m1, const G4ErrorMatrix& m2);
0244 // m = m1 - m2;
0245 // Note that m -= m1 is always faster than m = m - m1.
0246
0247 G4ErrorMatrix dsum(const G4ErrorMatrix&, const G4ErrorMatrix&);
0248 // Direct sum of two matrices. The direct sum of A and B is the matrix
0249 //        A 0
0250 //        0 B
0251
0252 std::ostream& operator<<(std::ostream& s, const G4ErrorMatrix& q);
0253 // Read in, write out G4ErrorMatrix into a stream.
0254
0255 //
0256 // Specialized linear algebra functions
0257 //
0258
0259 G4ErrorMatrix qr_solve(const G4ErrorMatrix& A, const G4ErrorMatrix& b);
0260 G4ErrorMatrix qr_solve(G4ErrorMatrix* A, const G4ErrorMatrix& b);
0261 // Works like backsolve, except matrix does not need to be upper
0262 // triangular. For nonsquare matrix, it solves in the least square sense.
0263
0264 G4ErrorMatrix qr_inverse(const G4ErrorMatrix& A);
0265 G4ErrorMatrix qr_inverse(G4ErrorMatrix* A);
0266 // Finds the inverse of a matrix using QR decomposition.  Note, often what
0267 // you really want is solve or backsolve, they can be much quicker than
0268 // inverse in many calculations.
0269
0270 void qr_decomp(G4ErrorMatrix* A, G4ErrorMatrix* hsm);
0271 G4ErrorMatrix qr_decomp(G4ErrorMatrix* A);
0272 // Does a QR decomposition of a matrix.
0273
0274 void back_solve(const G4ErrorMatrix& R, G4ErrorMatrix* b);
0275 // Solves R*x = b where R is upper triangular.  Also has a variation that
0276 // solves a number of equations of this form in one step, where b is a matrix
0277 // with each column a different vector. See also solve.
0278
0279 void col_house(G4ErrorMatrix* a, const G4ErrorMatrix& v, G4double vnormsq,
0280                G4int row, G4int col, G4int row_start, G4int col_start);
0281 void col_house(G4ErrorMatrix* a, const G4ErrorMatrix& v, G4int row, G4int col,
0282                G4int row_start, G4int col_start);
0283 // Does a column Householder update.
0284
0285 void col_givens(G4ErrorMatrix* A, G4double c, G4double s, G4int k1, G4int k2,
0286                 G4int row_min = 1, G4int row_max = 0);
0287 // do a column Givens update
0288
0289 void row_givens(G4ErrorMatrix* A, G4double c, G4double s, G4int k1, G4int k2,
0290                 G4int col_min = 1, G4int col_max = 0);
0291 // do a row Givens update
0292
0293 // void givens(G4double a, G4double b, G4double *c, G4double *s);
0294 // algorithm 5.1.5 in Golub and Van Loan
0295
0296 // Returns a Householder vector to zero elements.
0297
0298 void house_with_update(G4ErrorMatrix* a, G4int row = 1, G4int col = 1);
0299 void house_with_update(G4ErrorMatrix* a, G4ErrorMatrix* v, G4int row = 1,
0300                        G4int col = 1);
0301 // Finds and does Householder reflection on matrix.
0302
0303 void row_house(G4ErrorMatrix* a, const G4ErrorMatrix& v, G4double vnormsq,
0304                G4int row, G4int col, G4int row_start, G4int col_start);
0305 void row_house(G4ErrorMatrix* a, const G4ErrorMatrix& v, G4int row, G4int col,
0306                G4int row_start, G4int col_start);
0307 // Does a row Householder update.
0308
0309 #include "G4ErrorMatrix.icc"
0310
0311 #endif