source: trunk/source/error_propagation/include/G4ErrorMatrix.hh @ 815

Last change on this file since 815 was 815, checked in by garnier, 16 years ago

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