ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ErrorMatrix.hh
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ErrorMatrix.hh
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // Class Description:
28 //
29 // Simplified version of CLHEP HepMatrix class
30 
31 // History:
32 // - Imported from CLHEP and modified: P. Arce May 2007
33 // --------------------------------------------------------------------
34 
35 #ifndef G4ErrorMatrix_hh
36 #define G4ErrorMatrix_hh
37 
38 #include <vector>
39 #include "G4Types.hh"
40 
41 class G4ErrorSymMatrix;
42 
43 typedef std::vector<G4double >::iterator G4ErrorMatrixIter;
44 typedef std::vector<G4double >::const_iterator G4ErrorMatrixConstIter;
45 
46 class G4ErrorMatrix
47 {
48 
49  public: // with description
50 
51  G4ErrorMatrix();
52  // Default constructor. Gives 0 x 0 matrix.
53  // Another G4ErrorMatrix can be assigned to it.
54 
55  G4ErrorMatrix(G4int p, G4int q);
56  // Constructor. Gives an unitialized p x q matrix.
57 
58  G4ErrorMatrix(G4int p, G4int q, G4int i);
59  // Constructor. Gives an initialized p x q matrix.
60  // If i=0, it is initialized to all 0. If i=1, the diagonal elements
61  // are set to 1.0.
62 
63  G4ErrorMatrix(const G4ErrorMatrix &m1);
64  // Copy constructor.
65 
66  G4ErrorMatrix(const G4ErrorSymMatrix &m1);
67  // Constructors from G4ErrorSymG4ErrorMatrix, DiagG4ErrorMatrix and Vector.
68 
69  virtual ~G4ErrorMatrix();
70  // Destructor.
71 
72  inline virtual G4int num_row() const;
73  // Returns the number of rows.
74 
75  inline virtual G4int num_col() const;
76  // Returns the number of columns.
77 
78  inline virtual const G4double & operator()(G4int row, G4int col) const;
79  inline virtual G4double & operator()(G4int row, G4int col);
80  // Read or write a matrix element.
81  // ** Note that the indexing starts from (1,1). **
82 
83  G4ErrorMatrix & operator *= (G4double t);
84  // Multiply a G4ErrorMatrix by a floating number.
85 
86  G4ErrorMatrix & operator /= (G4double t);
87  // Divide a G4ErrorMatrix by a floating number.
88 
89  G4ErrorMatrix & operator += ( const G4ErrorMatrix &m2);
90  G4ErrorMatrix & operator += ( const G4ErrorSymMatrix &m2);
91  G4ErrorMatrix & operator -= ( const G4ErrorMatrix &m2);
92  G4ErrorMatrix & operator -= ( const G4ErrorSymMatrix &m2);
93  // Add or subtract a G4ErrorMatrix.
94  // When adding/subtracting Vector, G4ErrorMatrix must have num_col of one.
95 
96  G4ErrorMatrix & operator = ( const G4ErrorMatrix &m2);
97  G4ErrorMatrix & operator = ( const G4ErrorSymMatrix &m2);
98  // Assignment operators.
99 
100  G4ErrorMatrix operator- () const;
101  // unary minus, ie. flip the sign of each element.
102 
103  G4ErrorMatrix apply(G4double (*f)(G4double, G4int, G4int)) const;
104  // Apply a function to all elements of the matrix.
105 
106  G4ErrorMatrix T() const;
107  // Returns the transpose of a G4ErrorMatrix.
108 
109  G4ErrorMatrix sub(G4int min_row, G4int max_row,
110  G4int min_col, G4int max_col) const;
111  // Returns a sub matrix of a G4ErrorMatrix.
112  // WARNING: rows and columns are numbered from 1
113 
114  void sub(G4int row, G4int col, const G4ErrorMatrix &m1);
115  // Sub matrix of this G4ErrorMatrix is replaced with m1.
116  // WARNING: rows and columns are numbered from 1
117 
118  inline G4ErrorMatrix inverse(G4int& ierr) const;
119  // Invert a G4ErrorMatrix. G4ErrorMatrix must be square and is not changed.
120  // Returns ierr = 0 (zero) when successful, otherwise non-zero.
121 
122  virtual void invert(G4int& ierr);
123  // Invert a G4ErrorMatrix. G4ErrorMatrix must be square.
124  // N.B. the contents of the matrix are replaced by the inverse.
125  // Returns ierr = 0 (zero) when successful, otherwise non-zero.
126  // This method has less overhead then inverse().
127 
128  G4double determinant() const;
129  // calculate the determinant of the matrix.
130 
131  G4double trace() const;
132  // calculate the trace of the matrix (sum of diagonal elements).
133 
134  class G4ErrorMatrix_row
135  {
136  typedef std::vector<G4double >::const_iterator G4ErrorMatrixConstIter;
137  public:
138  inline G4ErrorMatrix_row(G4ErrorMatrix&,G4int);
139  G4double & operator[](G4int);
140  private:
141  G4ErrorMatrix& _a;
142  G4int _r;
143  };
144  class G4ErrorMatrix_row_const
145  {
146  public:
147  inline G4ErrorMatrix_row_const (const G4ErrorMatrix&,G4int);
148  const G4double & operator[](G4int) const;
149  private:
150  const G4ErrorMatrix& _a;
151  G4int _r;
152  };
153  // helper classes for implementing m[i][j]
154 
155  inline G4ErrorMatrix_row operator[] (G4int);
156  inline const G4ErrorMatrix_row_const operator[] (G4int) const;
157  // Read or write a matrix element.
158  // While it may not look like it, you simply do m[i][j] to get an element.
159  // ** Note that the indexing starts from [0][0]. **
160 
161  protected:
162 
163  virtual inline G4int num_size() const;
164  virtual void invertHaywood4(G4int& ierr);
165  virtual void invertHaywood5(G4int& ierr);
166  virtual void invertHaywood6(G4int& ierr);
167 
168  public:
169 
170  static void error(const char *s);
171 
172  private:
173 
174  friend class G4ErrorMatrix_row;
175  friend class G4ErrorMatrix_row_const;
176  friend class G4ErrorSymMatrix;
177  // Friend classes.
178 
179  friend G4ErrorMatrix operator+(const G4ErrorMatrix &m1,
180  const G4ErrorMatrix &m2);
181  friend G4ErrorMatrix operator-(const G4ErrorMatrix &m1,
182  const G4ErrorMatrix &m2);
183  friend G4ErrorMatrix operator*(const G4ErrorMatrix &m1,
184  const G4ErrorMatrix &m2);
185  friend G4ErrorMatrix operator*(const G4ErrorMatrix &m1,
186  const G4ErrorSymMatrix &m2);
187  friend G4ErrorMatrix operator*(const G4ErrorSymMatrix &m1,
188  const G4ErrorMatrix &m2);
189  friend G4ErrorMatrix operator*(const G4ErrorSymMatrix &m1,
190  const G4ErrorSymMatrix &m2);
191  // Multiply a G4ErrorMatrix by a G4ErrorMatrix or Vector.
192 
193  // solve the system of linear eq
194 
195  friend G4ErrorMatrix qr_solve(G4ErrorMatrix *, const G4ErrorMatrix &b);
196  friend void tridiagonal(G4ErrorSymMatrix *a,G4ErrorMatrix *hsm);
197  friend void row_house(G4ErrorMatrix *,const G4ErrorMatrix &, G4double,
198  G4int, G4int, G4int, G4int);
199  friend void back_solve(const G4ErrorMatrix &R, G4ErrorMatrix *b);
200  friend void col_givens(G4ErrorMatrix *A, G4double c,
201  G4double s, G4int k1, G4int k2,
202  G4int rowmin, G4int rowmax);
203  // Does a column Givens update.
204 
205  friend void row_givens(G4ErrorMatrix *A, G4double c,
206  G4double s, G4int k1, G4int k2,
207  G4int colmin, G4int colmax);
208  friend void col_house(G4ErrorMatrix *,const G4ErrorMatrix &, G4double,
209  G4int, G4int, G4int, G4int);
210  friend void house_with_update(G4ErrorMatrix *a,G4int row,G4int col);
211  friend void house_with_update(G4ErrorMatrix *a,G4ErrorMatrix *v,
212  G4int row, G4int col);
213  friend void house_with_update2(G4ErrorSymMatrix *a,G4ErrorMatrix *v,
214  G4int row, G4int col);
215 
216  G4int dfact_matrix(G4double &det, G4int *ir);
217  // factorize the matrix. If successful, the return code is 0. On
218  // return, det is the determinant and ir[] is row-interchange
219  // matrix. See CERNLIB's DFACT routine.
220 
221  G4int dfinv_matrix(G4int *ir);
222  // invert the matrix. See CERNLIB DFINV.
223 
224  std::vector<G4double > m;
225 
226  G4int nrow, ncol;
227  G4int size;
228 };
229 
230 // Operations other than member functions for G4ErrorMatrix
231 
232 G4ErrorMatrix operator*(const G4ErrorMatrix &m1, const G4ErrorMatrix &m2);
233 G4ErrorMatrix operator*(G4double t, const G4ErrorMatrix &m1);
234 G4ErrorMatrix operator*(const G4ErrorMatrix &m1, G4double t);
235 // Multiplication operators
236 // Note that m *= m1 is always faster than m = m * m1.
237 
238 G4ErrorMatrix operator/(const G4ErrorMatrix &m1, G4double t);
239 // m = m1 / t. (m /= t is faster if you can use it.)
240 
241 G4ErrorMatrix operator+(const G4ErrorMatrix &m1, const G4ErrorMatrix &m2);
242 // m = m1 + m2;
243 // Note that m += m1 is always faster than m = m + m1.
244 
245 G4ErrorMatrix operator-(const G4ErrorMatrix &m1, const G4ErrorMatrix &m2);
246 // m = m1 - m2;
247 // Note that m -= m1 is always faster than m = m - m1.
248 
249 G4ErrorMatrix dsum(const G4ErrorMatrix&, const G4ErrorMatrix&);
250 // Direct sum of two matrices. The direct sum of A and B is the matrix
251 // A 0
252 // 0 B
253 
254 std::ostream& operator<<(std::ostream &s, const G4ErrorMatrix &q);
255 // Read in, write out G4ErrorMatrix into a stream.
256 
257 //
258 // Specialized linear algebra functions
259 //
260 
261 G4ErrorMatrix qr_solve(const G4ErrorMatrix &A, const G4ErrorMatrix &b);
262 G4ErrorMatrix qr_solve(G4ErrorMatrix *A, const G4ErrorMatrix &b);
263 // Works like backsolve, except matrix does not need to be upper
264 // triangular. For nonsquare matrix, it solves in the least square sense.
265 
266 G4ErrorMatrix qr_inverse(const G4ErrorMatrix &A);
267 G4ErrorMatrix qr_inverse(G4ErrorMatrix *A);
268 // Finds the inverse of a matrix using QR decomposition. Note, often what
269 // you really want is solve or backsolve, they can be much quicker than
270 // inverse in many calculations.
271 
272 
273 void qr_decomp(G4ErrorMatrix *A, G4ErrorMatrix *hsm);
274 G4ErrorMatrix qr_decomp(G4ErrorMatrix *A);
275 // Does a QR decomposition of a matrix.
276 
277 void back_solve(const G4ErrorMatrix &R, G4ErrorMatrix *b);
278 // Solves R*x = b where R is upper triangular. Also has a variation that
279 // solves a number of equations of this form in one step, where b is a matrix
280 // with each column a different vector. See also solve.
281 
282 void col_house(G4ErrorMatrix *a, const G4ErrorMatrix &v, G4double vnormsq,
283  G4int row, G4int col, G4int row_start, G4int col_start);
284 void col_house(G4ErrorMatrix *a, const G4ErrorMatrix &v, G4int row, G4int col,
285  G4int row_start, G4int col_start);
286 // Does a column Householder update.
287 
288 void col_givens(G4ErrorMatrix *A, G4double c, G4double s,
289  G4int k1, G4int k2, G4int row_min=1, G4int row_max=0);
290 // do a column Givens update
291 
292 void row_givens(G4ErrorMatrix *A, G4double c, G4double s,
293  G4int k1, G4int k2, G4int col_min=1, G4int col_max=0);
294 // do a row Givens update
295 
296 // void givens(G4double a, G4double b, G4double *c, G4double *s);
297 // algorithm 5.1.5 in Golub and Van Loan
298 
299 // Returns a Householder vector to zero elements.
300 
301 void house_with_update(G4ErrorMatrix *a, G4int row=1, G4int col=1);
302 void house_with_update(G4ErrorMatrix *a, G4ErrorMatrix *v,
303  G4int row=1, G4int col=1);
304 // Finds and does Householder reflection on matrix.
305 
306 void row_house(G4ErrorMatrix *a, const G4ErrorMatrix &v, G4double vnormsq,
307  G4int row, G4int col, G4int row_start, G4int col_start);
308 void row_house(G4ErrorMatrix *a, const G4ErrorMatrix &v, G4int row, G4int col,
309  G4int row_start, G4int col_start);
310 // Does a row Householder update.
311 
312 #include "G4ErrorMatrix.icc"
313 
314 #endif