source: trunk/source/error_propagation/include/G4ErrorSymMatrix.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: 11.1 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: G4ErrorSymMatrix.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 HepSymMatrix class.
32
33// History:
34// - Imported from CLHEP and modified:   P. Arce    May 2007
35// --------------------------------------------------------------------
36
37#ifndef G4ErrorSymMatrix_hh
38#define G4ErrorSymMatrix_hh
39
40#include <vector>
41
42class G4ErrorMatrix;
43
44class G4ErrorSymMatrix
45{
46 public:  //with description
47
48   inline G4ErrorSymMatrix();
49   // Default constructor. Gives 0x0 symmetric matrix.
50   // Another G4ErrorSymMatrix can be assigned to it.
51
52   explicit G4ErrorSymMatrix(G4int p);
53   G4ErrorSymMatrix(G4int p, G4int);
54   // Constructor. Gives p x p symmetric matrix.
55   // With a second argument, the matrix is initialized. 0 means a zero
56   // matrix, 1 means the identity matrix.
57
58   G4ErrorSymMatrix(const G4ErrorSymMatrix &m1);
59   // Copy constructor.
60
61   // Constructor from DiagMatrix
62
63   virtual ~G4ErrorSymMatrix();
64   // Destructor.
65
66   inline G4int num_row() const;
67   inline G4int num_col() const;
68   // Returns number of rows/columns.
69
70   const G4double & operator()(G4int row, G4int col) const; 
71   G4double & operator()(G4int row, G4int col);
72   // Read and write a G4ErrorSymMatrix element.
73   // ** Note that indexing starts from (1,1). **
74
75   const G4double & fast(G4int row, G4int col) const;
76   G4double & fast(G4int row, G4int col);
77   // fast element access.
78   // Must be row>=col;
79   // ** Note that indexing starts from (1,1). **
80
81   void assign(const G4ErrorMatrix &m2);
82   // Assigns m2 to s, assuming m2 is a symmetric matrix.
83
84   void assign(const G4ErrorSymMatrix &m2);
85   // Another form of assignment. For consistency.
86
87   G4ErrorSymMatrix & operator*=(G4double t);
88   // Multiply a G4ErrorSymMatrix by a floating number.
89
90   G4ErrorSymMatrix & operator/=(G4double t); 
91   // Divide a G4ErrorSymMatrix by a floating number.
92
93   G4ErrorSymMatrix & operator+=( const G4ErrorSymMatrix &m2);
94   G4ErrorSymMatrix & operator-=( const G4ErrorSymMatrix &m2);
95   // Add or subtract a G4ErrorSymMatrix.
96
97   G4ErrorSymMatrix & operator=( const G4ErrorSymMatrix &m2);
98   // Assignment operators. Notice that there is no G4ErrorSymMatrix = Matrix.
99
100   G4ErrorSymMatrix operator- () const;
101   // unary minus, ie. flip the sign of each element.
102
103   G4ErrorSymMatrix T() const;
104   // Returns the transpose of a G4ErrorSymMatrix (which is itself).
105
106   G4ErrorSymMatrix apply(G4double (*f)(G4double, G4int, G4int)) const;
107   // Apply a function to all elements of the matrix.
108
109   G4ErrorSymMatrix similarity(const G4ErrorMatrix &m1) const;
110   G4ErrorSymMatrix similarity(const G4ErrorSymMatrix &m1) const;
111   // Returns m1*s*m1.T().
112
113   G4ErrorSymMatrix similarityT(const G4ErrorMatrix &m1) const;
114   // temporary. test of new similarity.
115   // Returns m1.T()*s*m1.
116
117   G4ErrorSymMatrix sub(G4int min_row, G4int max_row) const;
118   // Returns a sub matrix of a G4ErrorSymMatrix.
119
120   void sub(G4int row, const G4ErrorSymMatrix &m1);
121   // Sub matrix of this G4ErrorSymMatrix is replaced with m1.
122
123   G4ErrorSymMatrix sub(G4int min_row, G4int max_row);
124   // SGI CC bug. I have to have both with/without const. I should not need
125   // one without const.
126
127   inline G4ErrorSymMatrix inverse(G4int &ifail) const;
128   // Invert a Matrix. The matrix is not changed
129   // Returns 0 when successful, otherwise non-zero.
130
131   void invert(G4int &ifail);
132   // Invert a Matrix.
133   // N.B. the contents of the matrix are replaced by the inverse.
134   // Returns ierr = 0 when successful, otherwise non-zero.
135   // This method has less overhead then inverse().
136
137   G4double determinant() const;
138   // calculate the determinant of the matrix.
139
140   G4double trace() const;
141   // calculate the trace of the matrix (sum of diagonal elements).
142
143   class G4ErrorSymMatrix_row
144   {
145     public:
146      inline G4ErrorSymMatrix_row(G4ErrorSymMatrix&,G4int);
147      inline G4double & operator[](G4int);
148     private:
149      G4ErrorSymMatrix& _a;
150      G4int _r;
151   };
152   class G4ErrorSymMatrix_row_const
153   {
154     public:
155      inline G4ErrorSymMatrix_row_const(const G4ErrorSymMatrix&,G4int);
156      inline const G4double & operator[](G4int) const;
157     private:
158      const G4ErrorSymMatrix& _a;
159      G4int _r;
160   };
161   // helper class to implement m[i][j]
162
163   inline G4ErrorSymMatrix_row operator[] (G4int);
164   inline G4ErrorSymMatrix_row_const operator[] (G4int) const;
165   // Read or write a matrix element.
166   // While it may not look like it, you simply do m[i][j] to get an
167   // element.
168   // ** Note that the indexing starts from [0][0]. **
169
170   // Special-case inversions for 5x5 and 6x6 symmetric positive definite:
171   // These set ifail=0 and invert if the matrix was positive definite;
172   // otherwise ifail=1 and the matrix is left unaltered.
173
174   void invertCholesky5 (G4int &ifail); 
175   void invertCholesky6 (G4int &ifail);
176
177   // Inversions for 5x5 and 6x6 forcing use of specific methods:  The
178   // behavior (though not the speed) will be identical to invert(ifail).
179
180   void invertHaywood4 (G4int & ifail); 
181   void invertHaywood5 (G4int &ifail); 
182   void invertHaywood6 (G4int &ifail);
183   void invertBunchKaufman (G4int &ifail); 
184
185 protected:
186
187   inline G4int num_size() const;
188 
189 private:
190
191   friend class G4ErrorSymMatrix_row;
192   friend class G4ErrorSymMatrix_row_const;
193   friend class G4ErrorMatrix;
194
195   friend void tridiagonal(G4ErrorSymMatrix *a, G4ErrorMatrix *hsm);
196   friend G4double condition(const G4ErrorSymMatrix &m);
197   friend void diag_step(G4ErrorSymMatrix *t, G4int begin, G4int end);
198   friend void diag_step(G4ErrorSymMatrix *t, G4ErrorMatrix *u,
199                         G4int begin, G4int end);
200   friend G4ErrorMatrix diagonalize(G4ErrorSymMatrix *s);
201   friend void house_with_update2(G4ErrorSymMatrix *a, G4ErrorMatrix *v,
202                                  G4int row, G4int col);
203
204   friend G4ErrorSymMatrix operator+(const G4ErrorSymMatrix &m1, 
205                                     const G4ErrorSymMatrix &m2);
206   friend G4ErrorSymMatrix operator-(const G4ErrorSymMatrix &m1, 
207                                     const G4ErrorSymMatrix &m2);
208   friend G4ErrorMatrix operator*(const G4ErrorSymMatrix &m1,
209                                  const G4ErrorSymMatrix &m2);
210   friend G4ErrorMatrix operator*(const G4ErrorSymMatrix &m1,
211                                  const G4ErrorMatrix &m2);
212   friend G4ErrorMatrix operator*(const G4ErrorMatrix &m1,
213                                  const G4ErrorSymMatrix &m2);
214   // Multiply a Matrix by a Matrix or Vector.
215   
216   // Returns v * v.T();
217
218   std::vector<G4double > m;
219
220   G4int nrow;
221   G4int size;   // total number of elements
222
223   static G4double posDefFraction5x5;
224   static G4double adjustment5x5;
225   static const  G4double CHOLESKY_THRESHOLD_5x5;
226   static const  G4double CHOLESKY_CREEP_5x5;
227
228   static G4double posDefFraction6x6;
229   static G4double adjustment6x6;
230   static const G4double CHOLESKY_THRESHOLD_6x6;
231   static const G4double CHOLESKY_CREEP_6x6;
232
233   void invert4  (G4int & ifail);
234   void invert5  (G4int & ifail);
235   void invert6  (G4int & ifail);
236};
237
238//
239// Operations other than member functions for Matrix, G4ErrorSymMatrix,
240// DiagMatrix and Vectors
241//
242
243std::ostream& operator<<(std::ostream &s, const G4ErrorSymMatrix &q);
244// Write out Matrix, G4ErrorSymMatrix, DiagMatrix and Vector into ostream.
245
246G4ErrorMatrix operator*(const G4ErrorMatrix &m1,
247                        const G4ErrorSymMatrix &m2);
248G4ErrorMatrix operator*(const G4ErrorSymMatrix &m1,
249                        const G4ErrorMatrix &m2);
250G4ErrorMatrix operator*(const G4ErrorSymMatrix &m1,
251                        const G4ErrorSymMatrix &m2);
252G4ErrorSymMatrix operator*(G4double t, const G4ErrorSymMatrix &s1);
253G4ErrorSymMatrix operator*(const G4ErrorSymMatrix &s1, G4double t);
254// Multiplication operators.
255// Note that m *= m1 is always faster than m = m * m1
256
257G4ErrorSymMatrix operator/(const G4ErrorSymMatrix &m1, G4double t);
258// s = s1 / t. (s /= t is faster if you can use it.)
259
260G4ErrorMatrix operator+(const G4ErrorMatrix &m1,
261                        const G4ErrorSymMatrix &s2);
262G4ErrorMatrix operator+(const G4ErrorSymMatrix &s1,
263                        const G4ErrorMatrix &m2);
264G4ErrorSymMatrix operator+(const G4ErrorSymMatrix &s1,
265                           const G4ErrorSymMatrix &s2);
266// Addition operators
267
268G4ErrorMatrix operator-(const G4ErrorMatrix &m1,
269                        const G4ErrorSymMatrix &s2);
270G4ErrorMatrix operator-(const G4ErrorSymMatrix &m1,
271                        const G4ErrorMatrix &m2);
272G4ErrorSymMatrix operator-(const G4ErrorSymMatrix &s1,
273                           const G4ErrorSymMatrix &s2);
274// subtraction operators
275
276G4ErrorSymMatrix dsum(const G4ErrorSymMatrix &s1,
277                      const G4ErrorSymMatrix &s2);
278// Direct sum of two symmetric matrices;
279
280G4double condition(const G4ErrorSymMatrix &m);
281// Find the conditon number of a symmetric matrix.
282
283void diag_step(G4ErrorSymMatrix *t, G4int begin, G4int end);
284void diag_step(G4ErrorSymMatrix *t, G4ErrorMatrix *u, G4int begin, G4int end);
285// Implicit symmetric QR step with Wilkinson Shift
286
287G4ErrorMatrix diagonalize(G4ErrorSymMatrix *s);
288// Diagonalize a symmetric matrix.
289// It returns the matrix U so that s_old = U * s_diag * U.T()
290
291void house_with_update2(G4ErrorSymMatrix *a, G4ErrorMatrix *v,
292                        G4int row=1, G4int col=1);
293// Finds and does Householder reflection on matrix.
294
295void tridiagonal(G4ErrorSymMatrix *a, G4ErrorMatrix *hsm);
296G4ErrorMatrix tridiagonal(G4ErrorSymMatrix *a);
297// Does a Householder tridiagonalization of a symmetric matrix.
298
299#include "G4ErrorSymMatrix.icc"
300
301#endif
Note: See TracBrowser for help on using the repository browser.