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