| 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 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 public: // with description
49 G4ErrorMatrix();
50 // Default constructor. Gives 0 x 0 matrix.
51 // Another G4ErrorMatrix can be assigned to it.
52
53 G4ErrorMatrix(G4int p, G4int q);
54 // Constructor. Gives an unitialized p x q matrix.
55
56 G4ErrorMatrix(G4int p, G4int q, G4int i);
57 // Constructor. Gives an initialized p x q matrix.
58 // If i=0, it is initialized to all 0. If i=1, the diagonal elements
59 // are set to 1.0.
60
61 G4ErrorMatrix(const G4ErrorMatrix& m1);
62 G4ErrorMatrix(G4ErrorMatrix&&) = default;
63 // Copy and move constructor.
64
65 G4ErrorMatrix(const G4ErrorSymMatrix& m1);
66 // Constructors from G4ErrorSymG4ErrorMatrix, DiagG4ErrorMatrix and Vector.
67
68 virtual ~G4ErrorMatrix();
69 // Destructor.
70
71 inline virtual G4int num_row() const;
72 // Returns the number of rows.
73
74 inline virtual G4int num_col() const;
75 // Returns the number of columns.
76
77 inline virtual const G4double& operator()(G4int row, G4int col) const;
78 inline virtual G4double& operator()(G4int row, G4int col);
79 // Read or write a matrix element.
80 // ** Note that the indexing starts from (1,1). **
81
82 G4ErrorMatrix& operator*=(G4double t);
83 // Multiply a G4ErrorMatrix by a floating number.
84
85 G4ErrorMatrix& operator/=(G4double t);
86 // Divide a G4ErrorMatrix by a floating number.
87
88 G4ErrorMatrix& operator+=(const G4ErrorMatrix& m2);
89 G4ErrorMatrix& operator+=(const G4ErrorSymMatrix& m2);
90 G4ErrorMatrix& operator-=(const G4ErrorMatrix& m2);
91 G4ErrorMatrix& operator-=(const G4ErrorSymMatrix& m2);
92 // Add or subtract a G4ErrorMatrix.
93 // When adding/subtracting Vector, G4ErrorMatrix must have num_col of one.
94
95 G4ErrorMatrix& operator=(const G4ErrorMatrix& m2);
96 G4ErrorMatrix& operator=(const G4ErrorSymMatrix& m2);
97 G4ErrorMatrix& operator=(G4ErrorMatrix&&) = default;
98 // Assignment and move 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, G4int min_col,
110 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
138 public:
139 inline G4ErrorMatrix_row(G4ErrorMatrix&, G4int);
140 G4double& operator[](G4int);
141
142 private:
143 G4ErrorMatrix& _a;
144 G4int _r;
145 };
146 class G4ErrorMatrix_row_const
147 {
148 public:
149 inline G4ErrorMatrix_row_const(const G4ErrorMatrix&, G4int);
150 const G4double& operator[](G4int) const;
151
152 private:
153 const G4ErrorMatrix& _a;
154 G4int _r;
155 };
156 // helper classes for implementing m[i][j]
157
158 inline G4ErrorMatrix_row operator[](G4int);
159 inline const G4ErrorMatrix_row_const operator[](G4int) const;
160 // Read or write a matrix element.
161 // While it may not look like it, you simply do m[i][j] to get an element.
162 // ** Note that the indexing starts from [0][0]. **
163
164 protected:
165 virtual inline G4int num_size() const;
166 virtual void invertHaywood4(G4int& ierr);
167 virtual void invertHaywood5(G4int& ierr);
168 virtual void invertHaywood6(G4int& ierr);
169
170 public:
171 static void error(const char* s);
172
173 private:
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, G4int,
198 G4int, G4int, G4int);
199 friend void back_solve(const G4ErrorMatrix& R, G4ErrorMatrix* b);
200 friend void col_givens(G4ErrorMatrix* A, G4double c, G4double s, G4int k1,
201 G4int k2, G4int rowmin, G4int rowmax);
202 // Does a column Givens update.
203
204 friend void row_givens(G4ErrorMatrix* A, G4double c, G4double s, G4int k1,
205 G4int k2, G4int colmin, G4int colmax);
206 friend void col_house(G4ErrorMatrix*, const G4ErrorMatrix&, G4double, G4int,
207 G4int, G4int, G4int);
208 friend void house_with_update(G4ErrorMatrix* a, G4int row, G4int col);
209 friend void house_with_update(G4ErrorMatrix* a, G4ErrorMatrix* v, G4int row,
210 G4int col);
211 friend void house_with_update2(G4ErrorSymMatrix* a, G4ErrorMatrix* v,
212 G4int row, G4int col);
213
214 G4int dfact_matrix(G4double& det, G4int* ir);
215 // factorize the matrix. If successful, the return code is 0. On
216 // return, det is the determinant and ir[] is row-interchange
217 // matrix. See CERNLIB's DFACT routine.
218
219 G4int dfinv_matrix(G4int* ir);
220 // invert the matrix. See CERNLIB DFINV.
221
222 std::vector<G4double> m;
223
224 G4int nrow, ncol;
225 G4int size;
226 };
227
228 // Operations other than member functions for G4ErrorMatrix
229
230 G4ErrorMatrix operator*(const G4ErrorMatrix& m1, const G4ErrorMatrix& m2);
231 G4ErrorMatrix operator*(G4double t, const G4ErrorMatrix& m1);
232 G4ErrorMatrix operator*(const G4ErrorMatrix& m1, G4double t);
233 // Multiplication operators
234 // Note that m *= m1 is always faster than m = m * m1.
235
236 G4ErrorMatrix operator/(const G4ErrorMatrix& m1, G4double t);
237 // m = m1 / t. (m /= t is faster if you can use it.)
238
239 G4ErrorMatrix operator+(const G4ErrorMatrix& m1, const G4ErrorMatrix& m2);
240 // m = m1 + m2;
241 // Note that m += m1 is always faster than m = m + m1.
242
243 G4ErrorMatrix operator-(const G4ErrorMatrix& m1, const G4ErrorMatrix& m2);
244 // m = m1 - m2;
245 // Note that m -= m1 is always faster than m = m - m1.
246
247 G4ErrorMatrix dsum(const G4ErrorMatrix&, const G4ErrorMatrix&);
248 // Direct sum of two matrices. The direct sum of A and B is the matrix
249 // A 0
250 // 0 B
251
252 std::ostream& operator<<(std::ostream& s, const G4ErrorMatrix& q);
253 // Read in, write out G4ErrorMatrix into a stream.
254
255 //
256 // Specialized linear algebra functions
257 //
258
259 G4ErrorMatrix qr_solve(const G4ErrorMatrix& A, const G4ErrorMatrix& b);
260 G4ErrorMatrix qr_solve(G4ErrorMatrix* A, const G4ErrorMatrix& b);
261 // Works like backsolve, except matrix does not need to be upper
262 // triangular. For nonsquare matrix, it solves in the least square sense.
263
264 G4ErrorMatrix qr_inverse(const G4ErrorMatrix& A);
265 G4ErrorMatrix qr_inverse(G4ErrorMatrix* A);
266 // Finds the inverse of a matrix using QR decomposition. Note, often what
267 // you really want is solve or backsolve, they can be much quicker than
268 // inverse in many calculations.
269
270 void qr_decomp(G4ErrorMatrix* A, G4ErrorMatrix* hsm);
271 G4ErrorMatrix qr_decomp(G4ErrorMatrix* A);
272 // Does a QR decomposition of a matrix.
273
274 void back_solve(const G4ErrorMatrix& R, G4ErrorMatrix* b);
275 // Solves R*x = b where R is upper triangular. Also has a variation that
276 // solves a number of equations of this form in one step, where b is a matrix
277 // with each column a different vector. See also solve.
278
279 void col_house(G4ErrorMatrix* a, const G4ErrorMatrix& v, G4double vnormsq,
280 G4int row, G4int col, G4int row_start, G4int col_start);
281 void col_house(G4ErrorMatrix* a, const G4ErrorMatrix& v, G4int row, G4int col,
282 G4int row_start, G4int col_start);
283 // Does a column Householder update.
284
285 void col_givens(G4ErrorMatrix* A, G4double c, G4double s, G4int k1, G4int k2,
286 G4int row_min = 1, G4int row_max = 0);
287 // do a column Givens update
288
289 void row_givens(G4ErrorMatrix* A, G4double c, G4double s, G4int k1, G4int k2,
290 G4int col_min = 1, G4int col_max = 0);
291 // do a row Givens update
292
293 // void givens(G4double a, G4double b, G4double *c, G4double *s);
294 // algorithm 5.1.5 in Golub and Van Loan
295
296 // Returns a Householder vector to zero elements.
297
298 void house_with_update(G4ErrorMatrix* a, G4int row = 1, G4int col = 1);
299 void house_with_update(G4ErrorMatrix* a, G4ErrorMatrix* v, G4int row = 1,
300 G4int col = 1);
301 // Finds and does Householder reflection on matrix.
302
303 void row_house(G4ErrorMatrix* a, const G4ErrorMatrix& v, G4double vnormsq,
304 G4int row, G4int col, G4int row_start, G4int col_start);
305 void row_house(G4ErrorMatrix* a, const G4ErrorMatrix& v, G4int row, G4int col,
306 G4int row_start, G4int col_start);
307 // Does a row Householder update.
308
309 #include "G4ErrorMatrix.icc"
310
311 #endif
312