Geant4  10.02.p01
G4ErrorMatrix.hh
Go to the documentation of this file.
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 // $Id: G4ErrorMatrix.hh 66892 2013-01-17 10:57:59Z gunter $
27 //
28 // Class Description:
29 //
30 // Simplified version of CLHEP HepMatrix class
31 
32 // History:
33 // - Imported from CLHEP and modified: P. Arce May 2007
34 // --------------------------------------------------------------------
35 
36 #ifndef G4ErrorMatrix_hh
37 #define G4ErrorMatrix_hh
38 
39 #include <vector>
40 
42 
43 typedef std::vector<G4double >::iterator G4ErrorMatrixIter;
44 typedef std::vector<G4double >::const_iterator G4ErrorMatrixConstIter;
45 
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 
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 
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 
84  // Multiply a G4ErrorMatrix by a floating number.
85 
87  // Divide a G4ErrorMatrix by a floating number.
88 
93  // Add or subtract a G4ErrorMatrix.
94  // When adding/subtracting Vector, G4ErrorMatrix must have num_col of one.
95 
98  // Assignment operators.
99 
100  G4ErrorMatrix operator- () const;
101  // unary minus, ie. flip the sign of each element.
102 
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 
135  {
136  typedef std::vector<G4double >::const_iterator G4ErrorMatrixConstIter;
137  public:
140  private:
143  };
145  {
146  public:
148  const G4double & operator[](G4int) const;
149  private:
152  };
153  // helper classes for implementing m[i][j]
154 
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;
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 
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);
212  G4int row, G4int col);
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 
228 };
229 
230 // Operations other than member functions for G4ErrorMatrix
231 
235 // Multiplication operators
236 // Note that m *= m1 is always faster than m = m * m1.
237 
239 // m = m1 / t. (m /= t is faster if you can use it.)
240 
242 // m = m1 + m2;
243 // Note that m += m1 is always faster than m = m + m1.
244 
246 // m = m1 - m2;
247 // Note that m -= m1 is always faster than m = m - m1.
248 
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 
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 
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 
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 
289  G4int k1, G4int k2, G4int row_min=1, G4int row_max=0);
290 // do a column Givens update
291 
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);
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
void house_with_update(G4ErrorMatrix *a, G4int row=1, G4int col=1)
virtual G4int num_row() const
G4ErrorMatrix dsum(const G4ErrorMatrix &, const G4ErrorMatrix &)
G4ErrorMatrix operator*(const G4ErrorMatrix &m1, const G4ErrorMatrix &m2)
void back_solve(const G4ErrorMatrix &R, G4ErrorMatrix *b)
const G4double & operator[](G4int) const
G4ErrorMatrix inverse(G4int &ierr) const
G4ErrorMatrix operator/(const G4ErrorMatrix &m1, G4double t)
virtual G4int num_col() const
void givens(G4double a, G4double b, G4double *c, G4double *s)
void col_givens(G4ErrorMatrix *A, G4double c, G4double s, G4int k1, G4int k2, G4int row_min=1, G4int row_max=0)
friend void back_solve(const G4ErrorMatrix &R, G4ErrorMatrix *b)
virtual ~G4ErrorMatrix()
G4double a
Definition: TRTMaterials.hh:39
friend void col_house(G4ErrorMatrix *, const G4ErrorMatrix &, G4double, G4int, G4int, G4int, G4int)
friend void tridiagonal(G4ErrorSymMatrix *a, G4ErrorMatrix *hsm)
int G4int
Definition: G4Types.hh:78
G4ErrorMatrix & operator*=(G4double t)
G4ErrorMatrix & operator=(const G4ErrorMatrix &m2)
std::vector< G4double >::const_iterator G4ErrorMatrixConstIter
G4ErrorMatrix_row_const(const G4ErrorMatrix &, G4int)
G4ErrorMatrix sub(G4int min_row, G4int max_row, G4int min_col, G4int max_col) const
static const double s
Definition: G4SIunits.hh:168
G4ErrorMatrix operator-(const G4ErrorMatrix &m1, const G4ErrorMatrix &m2)
friend void house_with_update2(G4ErrorSymMatrix *a, G4ErrorMatrix *v, G4int row, G4int col)
friend G4ErrorMatrix operator*(const G4ErrorMatrix &m1, const G4ErrorMatrix &m2)
virtual void invert(G4int &ierr)
double A(double temperature)
friend void row_givens(G4ErrorMatrix *A, G4double c, G4double s, G4int k1, G4int k2, G4int colmin, G4int colmax)
friend G4ErrorMatrix operator+(const G4ErrorMatrix &m1, const G4ErrorMatrix &m2)
static const double m2
Definition: G4SIunits.hh:129
std::vector< G4double >::const_iterator G4ErrorMatrixConstIter
void row_givens(G4ErrorMatrix *A, G4double c, G4double s, G4int k1, G4int k2, G4int col_min=1, G4int col_max=0)
void col_house(G4ErrorMatrix *a, const G4ErrorMatrix &v, G4double vnormsq, G4int row, G4int col, G4int row_start, G4int col_start)
G4ErrorMatrix & operator/=(G4double t)
void row_house(G4ErrorMatrix *a, const G4ErrorMatrix &v, G4double vnormsq, G4int row, G4int col, G4int row_start, G4int col_start)
std::vector< G4double > m
G4int dfact_matrix(G4double &det, G4int *ir)
static void error(const char *s)
G4ErrorMatrix T() const
G4double trace() const
std::vector< G4double >::iterator G4ErrorMatrixIter
G4ErrorMatrix apply(G4double(*f)(G4double, G4int, G4int)) const
G4int dfinv_matrix(G4int *ir)
friend void row_house(G4ErrorMatrix *, const G4ErrorMatrix &, G4double, G4int, G4int, G4int, G4int)
virtual G4int num_size() const
virtual void invertHaywood4(G4int &ierr)
void qr_decomp(G4ErrorMatrix *A, G4ErrorMatrix *hsm)
G4double determinant() const
friend void col_givens(G4ErrorMatrix *A, G4double c, G4double s, G4int k1, G4int k2, G4int rowmin, G4int rowmax)
virtual void invertHaywood6(G4int &ierr)
G4ErrorMatrix_row operator[](G4int)
virtual void invertHaywood5(G4int &ierr)
friend G4ErrorMatrix qr_solve(G4ErrorMatrix *, const G4ErrorMatrix &b)
G4ErrorMatrix operator-() const
virtual const G4double & operator()(G4int row, G4int col) const
double G4double
Definition: G4Types.hh:76
G4ErrorMatrix qr_inverse(const G4ErrorMatrix &A)
G4ErrorMatrix & operator+=(const G4ErrorMatrix &m2)
G4ErrorMatrix & operator-=(const G4ErrorMatrix &m2)
friend void house_with_update(G4ErrorMatrix *a, G4int row, G4int col)
G4ErrorMatrix operator+(const G4ErrorMatrix &m1, const G4ErrorMatrix &m2)
G4ErrorMatrix_row(G4ErrorMatrix &, G4int)
G4ErrorMatrix qr_solve(const G4ErrorMatrix &A, const G4ErrorMatrix &b)
std::ostream & operator<<(std::ostream &s, const G4ErrorMatrix &q)