| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // ------------------------------------------- 27 // ------------------------------------------------------------
28 // GEANT 4 class implementation file 28 // GEANT 4 class implementation file
29 // ------------------------------------------- 29 // ------------------------------------------------------------
30 30
31 #include "globals.hh" 31 #include "globals.hh"
32 #include <iostream> 32 #include <iostream>
33 #include <cmath> 33 #include <cmath>
34 34
35 #include "G4ErrorSymMatrix.hh" 35 #include "G4ErrorSymMatrix.hh"
36 #include "G4ErrorMatrix.hh" 36 #include "G4ErrorMatrix.hh"
37 37
38 // Simple operation for all elements 38 // Simple operation for all elements
39 39
40 #define SIMPLE_UOP(OPER) << 40 #define SIMPLE_UOP(OPER) \
41 G4ErrorMatrixIter a = m.begin(); << 41 G4ErrorMatrixIter a=m.begin(); \
42 G4ErrorMatrixIter e = m.begin() + num_size() << 42 G4ErrorMatrixIter e=m.begin()+num_size(); \
43 for(; a < e; a++) << 43 for(;a<e; a++) (*a) OPER t;
44 (*a) OPER t; << 44
45 << 45 #define SIMPLE_BOP(OPER) \
46 #define SIMPLE_BOP(OPER) << 46 G4ErrorMatrixIter a=m.begin(); \
47 G4ErrorMatrixIter a = m.begin(); << 47 G4ErrorMatrixConstIter b=mat2.m.begin(); \
48 G4ErrorMatrixConstIter b = mat2.m.begin(); << 48 G4ErrorMatrixConstIter e=m.begin()+num_size(); \
49 G4ErrorMatrixConstIter e = m.begin() + num_s << 49 for(;a<e; a++, b++) (*a) OPER (*b);
50 for(; a < e; a++, b++) << 50
51 (*a) OPER(*b); << 51 #define SIMPLE_TOP(OPER) \
52 << 52 G4ErrorMatrixConstIter a=mat1.m.begin(); \
53 #define SIMPLE_TOP(OPER) << 53 G4ErrorMatrixConstIter b=mat2.m.begin(); \
54 G4ErrorMatrixConstIter a = mat1.m.begin(); << 54 G4ErrorMatrixIter t=mret.m.begin(); \
55 G4ErrorMatrixConstIter b = mat2.m.begin(); << 55 G4ErrorMatrixConstIter e=mat1.m.begin()+mat1.num_size(); \
56 G4ErrorMatrixIter t = mret.m.begin(); << 56 for( ;a<e; a++, b++, t++) (*t) = (*a) OPER (*b);
57 G4ErrorMatrixConstIter e = mat1.m.begin() + << 57
58 for(; a < e; a++, b++, t++) << 58 #define CHK_DIM_2(r1,r2,c1,c2,fun) \
59 (*t) = (*a) OPER(*b); << 59 if (r1!=r2 || c1!=c2) { \
60 << 60 G4ErrorMatrix::error("Range error in Matrix function " #fun "(1)."); \
61 #define CHK_DIM_2(r1, r2, c1, c2, fun) << 61 }
62 if(r1 != r2 || c1 != c2) << 62
63 { << 63 #define CHK_DIM_1(c1,r2,fun) \
64 G4ErrorMatrix::error("Range error in Matri << 64 if (c1!=r2) { \
65 } << 65 G4ErrorMatrix::error("Range error in Matrix function " #fun "(2)."); \
66 << 66 }
67 #define CHK_DIM_1(c1, r2, fun) <<
68 if(c1 != r2) <<
69 { <<
70 G4ErrorMatrix::error("Range error in Matri <<
71 } <<
72 67
73 // Constructors. (Default constructors are inl 68 // Constructors. (Default constructors are inlined and in .icc file)
74 69
75 G4ErrorSymMatrix::G4ErrorSymMatrix(G4int p) 70 G4ErrorSymMatrix::G4ErrorSymMatrix(G4int p)
76 : m(p * (p + 1) / 2) << 71 : m(p*(p+1)/2), nrow(p)
77 , nrow(p) <<
78 { 72 {
79 size = nrow * (nrow + 1) / 2; << 73 size = nrow * (nrow+1) / 2;
80 m.assign(size, 0); << 74 m.assign(size,0);
81 } 75 }
82 76
83 G4ErrorSymMatrix::G4ErrorSymMatrix(G4int p, G4 77 G4ErrorSymMatrix::G4ErrorSymMatrix(G4int p, G4int init)
84 : m(p * (p + 1) / 2) << 78 : m(p*(p+1)/2), nrow(p)
85 , nrow(p) <<
86 { 79 {
87 size = nrow * (nrow + 1) / 2; << 80 size = nrow * (nrow+1) / 2;
88 81
89 m.assign(size, 0); << 82 m.assign(size,0);
90 switch(init) << 83 switch(init)
91 { << 84 {
92 case 0: << 85 case 0:
93 break; 86 break;
94 << 87
95 case 1: << 88 case 1:
96 { <<
97 G4ErrorMatrixIter a = m.begin(); <<
98 for(G4int i = 1; i <= nrow; i++) <<
99 { 89 {
100 *a = 1.0; << 90 G4ErrorMatrixIter a = m.begin();
101 a += (i + 1); << 91 for(G4int i=1;i<=nrow;i++)
>> 92 {
>> 93 *a = 1.0;
>> 94 a += (i+1);
>> 95 }
>> 96 break;
102 } 97 }
103 break; << 98 default:
104 } <<
105 default: <<
106 G4ErrorMatrix::error("G4ErrorSymMatrix: 99 G4ErrorMatrix::error("G4ErrorSymMatrix: initialization must be 0 or 1.");
107 } << 100 }
108 } 101 }
109 102
110 // 103 //
111 // Destructor 104 // Destructor
112 // 105 //
113 106
114 G4ErrorSymMatrix::~G4ErrorSymMatrix() {} << 107 G4ErrorSymMatrix::~G4ErrorSymMatrix()
>> 108 {
>> 109 }
115 110
116 G4ErrorSymMatrix::G4ErrorSymMatrix(const G4Err << 111 G4ErrorSymMatrix::G4ErrorSymMatrix(const G4ErrorSymMatrix &mat1)
117 : m(mat1.size) << 112 : m(mat1.size), nrow(mat1.nrow), size(mat1.size)
118 , nrow(mat1.nrow) <<
119 , size(mat1.size) <<
120 { 113 {
121 m = mat1.m; << 114 m = mat1.m;
122 } 115 }
123 116
124 // 117 //
125 // 118 //
126 // Sub matrix 119 // Sub matrix
127 // 120 //
128 // 121 //
129 122
130 G4ErrorSymMatrix G4ErrorSymMatrix::sub(G4int m 123 G4ErrorSymMatrix G4ErrorSymMatrix::sub(G4int min_row, G4int max_row) const
131 { 124 {
132 G4ErrorSymMatrix mret(max_row - min_row + 1) << 125 G4ErrorSymMatrix mret(max_row-min_row+1);
133 if(max_row > num_row()) 126 if(max_row > num_row())
134 { << 127 { G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range"); }
135 G4ErrorMatrix::error("G4ErrorSymMatrix::su << 128 G4ErrorMatrixIter a = mret.m.begin();
136 } << 129 G4ErrorMatrixConstIter b1 = m.begin() + (min_row+2)*(min_row-1)/2;
137 G4ErrorMatrixIter a = mret.m.begin(); << 130 for(G4int irow=1; irow<=mret.num_row(); irow++)
138 G4ErrorMatrixConstIter b1 = m.begin() + (min <<
139 for(G4int irow = 1; irow <= mret.num_row(); <<
140 { 131 {
141 G4ErrorMatrixConstIter b = b1; 132 G4ErrorMatrixConstIter b = b1;
142 for(G4int icol = 1; icol <= irow; icol++) << 133 for(G4int icol=1; icol<=irow; icol++)
143 { 134 {
144 *(a++) = *(b++); 135 *(a++) = *(b++);
145 } 136 }
146 b1 += irow + min_row - 1; << 137 b1 += irow+min_row-1;
147 } 138 }
148 return mret; 139 return mret;
149 } 140 }
150 141
151 G4ErrorSymMatrix G4ErrorSymMatrix::sub(G4int m << 142 G4ErrorSymMatrix G4ErrorSymMatrix::sub(G4int min_row, G4int max_row)
152 { 143 {
153 G4ErrorSymMatrix mret(max_row - min_row + 1) << 144 G4ErrorSymMatrix mret(max_row-min_row+1);
154 if(max_row > num_row()) 145 if(max_row > num_row())
155 { << 146 { G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range"); }
156 G4ErrorMatrix::error("G4ErrorSymMatrix::su << 147 G4ErrorMatrixIter a = mret.m.begin();
157 } << 148 G4ErrorMatrixIter b1 = m.begin() + (min_row+2)*(min_row-1)/2;
158 G4ErrorMatrixIter a = mret.m.begin(); << 149 for(G4int irow=1; irow<=mret.num_row(); irow++)
159 G4ErrorMatrixIter b1 = m.begin() + (min_row <<
160 for(G4int irow = 1; irow <= mret.num_row(); <<
161 { 150 {
162 G4ErrorMatrixIter b = b1; 151 G4ErrorMatrixIter b = b1;
163 for(G4int icol = 1; icol <= irow; icol++) << 152 for(G4int icol=1; icol<=irow; icol++)
164 { 153 {
165 *(a++) = *(b++); 154 *(a++) = *(b++);
166 } 155 }
167 b1 += irow + min_row - 1; << 156 b1 += irow+min_row-1;
168 } 157 }
169 return mret; 158 return mret;
170 } 159 }
171 160
172 void G4ErrorSymMatrix::sub(G4int row, const G4 << 161 void G4ErrorSymMatrix::sub(G4int row,const G4ErrorSymMatrix &mat1)
173 { 162 {
174 if(row < 1 || row + mat1.num_row() - 1 > num << 163 if(row <1 || row+mat1.num_row()-1 > num_row() )
175 { << 164 { G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range"); }
176 G4ErrorMatrix::error("G4ErrorSymMatrix::su <<
177 } <<
178 G4ErrorMatrixConstIter a = mat1.m.begin(); 165 G4ErrorMatrixConstIter a = mat1.m.begin();
179 G4ErrorMatrixIter b1 = m.begin() + (row << 166 G4ErrorMatrixIter b1 = m.begin() + (row+2)*(row-1)/2;
180 for(G4int irow = 1; irow <= mat1.num_row(); << 167 for(G4int irow=1; irow<=mat1.num_row(); irow++)
181 { 168 {
182 G4ErrorMatrixIter b = b1; 169 G4ErrorMatrixIter b = b1;
183 for(G4int icol = 1; icol <= irow; icol++) << 170 for(G4int icol=1; icol<=irow; icol++)
184 { 171 {
185 *(b++) = *(a++); 172 *(b++) = *(a++);
186 } 173 }
187 b1 += irow + row - 1; << 174 b1 += irow+row-1;
188 } 175 }
189 } 176 }
190 177
191 // 178 //
192 // Direct sum of two matricies 179 // Direct sum of two matricies
193 // 180 //
194 181
195 G4ErrorSymMatrix dsum(const G4ErrorSymMatrix& << 182 G4ErrorSymMatrix dsum(const G4ErrorSymMatrix &mat1,
196 const G4ErrorSymMatrix& << 183 const G4ErrorSymMatrix &mat2)
197 { 184 {
198 G4ErrorSymMatrix mret(mat1.num_row() + mat2. 185 G4ErrorSymMatrix mret(mat1.num_row() + mat2.num_row(), 0);
199 mret.sub(1, mat1); << 186 mret.sub(1,mat1);
200 mret.sub(mat1.num_row() + 1, mat2); << 187 mret.sub(mat1.num_row()+1,mat2);
201 return mret; 188 return mret;
202 } 189 }
203 190
204 /* ------------------------------------------- 191 /* -----------------------------------------------------------------------
205 This section contains support routines for 192 This section contains support routines for matrix.h. This section contains
206 The two argument functions +,-. They call t 193 The two argument functions +,-. They call the copy constructor and +=,-=.
207 ------------------------------------------- 194 ----------------------------------------------------------------------- */
208 195
209 G4ErrorSymMatrix G4ErrorSymMatrix::operator-() << 196 G4ErrorSymMatrix G4ErrorSymMatrix::operator- () const
210 { 197 {
211 G4ErrorSymMatrix mat2(nrow); 198 G4ErrorSymMatrix mat2(nrow);
212 G4ErrorMatrixConstIter a = m.begin(); << 199 G4ErrorMatrixConstIter a=m.begin();
213 G4ErrorMatrixIter b = mat2.m.begin(); << 200 G4ErrorMatrixIter b=mat2.m.begin();
214 G4ErrorMatrixConstIter e = m.begin() + num_s << 201 G4ErrorMatrixConstIter e=m.begin()+num_size();
215 for(; a < e; a++, b++) << 202 for(;a<e; a++, b++) { (*b) = -(*a); }
216 { <<
217 (*b) = -(*a); <<
218 } <<
219 return mat2; 203 return mat2;
220 } 204 }
221 205
222 G4ErrorMatrix operator+(const G4ErrorMatrix& m << 206 G4ErrorMatrix operator+(const G4ErrorMatrix &mat1, const G4ErrorSymMatrix &mat2)
223 { 207 {
224 G4ErrorMatrix mret(mat1); 208 G4ErrorMatrix mret(mat1);
225 CHK_DIM_2(mat1.num_row(), mat2.num_row(), ma << 209 CHK_DIM_2(mat1.num_row(),mat2.num_row(), mat1.num_col(),mat2.num_col(),+);
226 mret += mat2; 210 mret += mat2;
227 return mret; 211 return mret;
228 } 212 }
229 213
230 G4ErrorMatrix operator+(const G4ErrorSymMatrix << 214 G4ErrorMatrix operator+(const G4ErrorSymMatrix &mat1, const G4ErrorMatrix &mat2)
231 { 215 {
232 G4ErrorMatrix mret(mat2); 216 G4ErrorMatrix mret(mat2);
233 CHK_DIM_2(mat1.num_row(), mat2.num_row(), ma << 217 CHK_DIM_2(mat1.num_row(),mat2.num_row(),mat1.num_col(),mat2.num_col(),+);
234 mret += mat1; 218 mret += mat1;
235 return mret; 219 return mret;
236 } 220 }
237 221
238 G4ErrorSymMatrix operator+(const G4ErrorSymMat << 222 G4ErrorSymMatrix operator+(const G4ErrorSymMatrix &mat1,
239 const G4ErrorSymMat << 223 const G4ErrorSymMatrix &mat2)
240 { 224 {
241 G4ErrorSymMatrix mret(mat1.nrow); 225 G4ErrorSymMatrix mret(mat1.nrow);
242 CHK_DIM_1(mat1.nrow, mat2.nrow, +); << 226 CHK_DIM_1(mat1.nrow, mat2.nrow,+);
243 SIMPLE_TOP(+) 227 SIMPLE_TOP(+)
244 return mret; 228 return mret;
245 } 229 }
246 230
247 // 231 //
248 // operator - 232 // operator -
249 // 233 //
250 234
251 G4ErrorMatrix operator-(const G4ErrorMatrix& m << 235 G4ErrorMatrix operator-(const G4ErrorMatrix &mat1, const G4ErrorSymMatrix &mat2)
252 { 236 {
253 G4ErrorMatrix mret(mat1); 237 G4ErrorMatrix mret(mat1);
254 CHK_DIM_2(mat1.num_row(), mat2.num_row(), ma << 238 CHK_DIM_2(mat1.num_row(),mat2.num_row(),mat1.num_col(),mat2.num_col(),-);
255 mret -= mat2; 239 mret -= mat2;
256 return mret; 240 return mret;
257 } 241 }
258 242
259 G4ErrorMatrix operator-(const G4ErrorSymMatrix << 243 G4ErrorMatrix operator-(const G4ErrorSymMatrix &mat1, const G4ErrorMatrix &mat2)
260 { 244 {
261 G4ErrorMatrix mret(mat1); 245 G4ErrorMatrix mret(mat1);
262 CHK_DIM_2(mat1.num_row(), mat2.num_row(), ma << 246 CHK_DIM_2(mat1.num_row(),mat2.num_row(),mat1.num_col(),mat2.num_col(),-);
263 mret -= mat2; 247 mret -= mat2;
264 return mret; 248 return mret;
265 } 249 }
266 250
267 G4ErrorSymMatrix operator-(const G4ErrorSymMat << 251 G4ErrorSymMatrix operator-(const G4ErrorSymMatrix &mat1,
268 const G4ErrorSymMat << 252 const G4ErrorSymMatrix &mat2)
269 { 253 {
270 G4ErrorSymMatrix mret(mat1.num_row()); 254 G4ErrorSymMatrix mret(mat1.num_row());
271 CHK_DIM_1(mat1.num_row(), mat2.num_row(), -) << 255 CHK_DIM_1(mat1.num_row(),mat2.num_row(),-);
272 SIMPLE_TOP(-) 256 SIMPLE_TOP(-)
273 return mret; 257 return mret;
274 } 258 }
275 259
276 /* ------------------------------------------- 260 /* -----------------------------------------------------------------------
277 This section contains support routines for 261 This section contains support routines for matrix.h. This file contains
278 The two argument functions *,/. They call c 262 The two argument functions *,/. They call copy constructor and then /=,*=.
279 ------------------------------------------- 263 ----------------------------------------------------------------------- */
280 264
281 G4ErrorSymMatrix operator/(const G4ErrorSymMat << 265 G4ErrorSymMatrix operator/(const G4ErrorSymMatrix &mat1,G4double t)
282 { 266 {
283 G4ErrorSymMatrix mret(mat1); 267 G4ErrorSymMatrix mret(mat1);
284 mret /= t; 268 mret /= t;
285 return mret; 269 return mret;
286 } 270 }
287 271
288 G4ErrorSymMatrix operator*(const G4ErrorSymMat << 272 G4ErrorSymMatrix operator*(const G4ErrorSymMatrix &mat1,G4double t)
289 { 273 {
290 G4ErrorSymMatrix mret(mat1); 274 G4ErrorSymMatrix mret(mat1);
291 mret *= t; 275 mret *= t;
292 return mret; 276 return mret;
293 } 277 }
294 278
295 G4ErrorSymMatrix operator*(G4double t, const G << 279 G4ErrorSymMatrix operator*(G4double t,const G4ErrorSymMatrix &mat1)
296 { 280 {
297 G4ErrorSymMatrix mret(mat1); 281 G4ErrorSymMatrix mret(mat1);
298 mret *= t; 282 mret *= t;
299 return mret; 283 return mret;
300 } 284 }
301 285
302 G4ErrorMatrix operator*(const G4ErrorMatrix& m << 286 G4ErrorMatrix operator*(const G4ErrorMatrix &mat1, const G4ErrorSymMatrix &mat2)
303 { 287 {
304 G4ErrorMatrix mret(mat1.num_row(), mat2.num_ << 288 G4ErrorMatrix mret(mat1.num_row(),mat2.num_col());
305 CHK_DIM_1(mat1.num_col(), mat2.num_row(), *) << 289 CHK_DIM_1(mat1.num_col(),mat2.num_row(),*);
306 G4ErrorMatrixConstIter mit1, mit2, sp, snp; << 290 G4ErrorMatrixConstIter mit1, mit2, sp,snp; //mit2=0
307 G4double temp; 291 G4double temp;
308 G4ErrorMatrixIter mir = mret.m.begin(); << 292 G4ErrorMatrixIter mir=mret.m.begin();
309 for(mit1 = mat1.m.begin(); << 293 for(mit1=mat1.m.begin();
310 mit1 < mat1.m.begin() + mat1.num_row() * << 294 mit1<mat1.m.begin()+mat1.num_row()*mat1.num_col();
311 { << 295 mit1 = mit2)
312 snp = mat2.m.begin(); << 296 {
313 for(int step = 1; step <= mat2.num_row(); << 297 snp=mat2.m.begin();
314 { << 298 for(int step=1;step<=mat2.num_row();++step)
315 mit2 = mit1; << 299 {
316 sp = snp; << 300 mit2=mit1;
317 snp += step; << 301 sp=snp;
318 temp = 0; << 302 snp+=step;
319 while(sp < snp) // Loop checking, 06.08 << 303 temp=0;
320 { << 304 while(sp<snp) // Loop checking, 06.08.2015, G.Cosmo
321 temp += *(sp++) * (*(mit2++)); << 305 { temp+=*(sp++)*(*(mit2++)); }
322 } << 306 if( step<mat2.num_row() ) { // only if we aren't on the last row
323 if(step < mat2.num_row()) << 307 sp+=step-1;
324 { // only if we aren't on the last row << 308 for(int stept=step+1;stept<=mat2.num_row();stept++)
325 sp += step - 1; << 309 {
326 for(int stept = step + 1; stept <= mat << 310 temp+=*sp*(*(mit2++));
327 { << 311 if(stept<mat2.num_row()) sp+=stept;
328 temp += *sp * (*(mit2++)); << 312 }
329 if(stept < mat2.num_row()) << 313 } // if(step
330 sp += stept; << 314 *(mir++)=temp;
331 } << 315 } // for(step
332 } // if(step << 316 } // for(mit1
333 *(mir++) = temp; <<
334 } // for(step <<
335 } // for(mit1 <<
336 return mret; 317 return mret;
337 } 318 }
338 319
339 G4ErrorMatrix operator*(const G4ErrorSymMatrix << 320 G4ErrorMatrix operator*(const G4ErrorSymMatrix &mat1, const G4ErrorMatrix &mat2)
340 { 321 {
341 G4ErrorMatrix mret(mat1.num_row(), mat2.num_ << 322 G4ErrorMatrix mret(mat1.num_row(),mat2.num_col());
342 CHK_DIM_1(mat1.num_col(), mat2.num_row(), *) << 323 CHK_DIM_1(mat1.num_col(),mat2.num_row(),*);
343 G4int step, stept; << 324 G4int step,stept;
344 G4ErrorMatrixConstIter mit1, mit2, sp, snp; << 325 G4ErrorMatrixConstIter mit1,mit2,sp,snp;
345 G4double temp; 326 G4double temp;
346 G4ErrorMatrixIter mir = mret.m.begin(); << 327 G4ErrorMatrixIter mir=mret.m.begin();
347 for(step = 1, snp = mat1.m.begin(); step <= << 328 for(step=1,snp=mat1.m.begin();step<=mat1.num_row();snp+=step++)
348 { 329 {
349 for(mit1 = mat2.m.begin(); mit1 < mat2.m.b << 330 for(mit1=mat2.m.begin();mit1<mat2.m.begin()+mat2.num_col();mit1++)
350 { 331 {
351 mit2 = mit1; << 332 mit2=mit1;
352 sp = snp; << 333 sp=snp;
353 temp = 0; << 334 temp=0;
354 while(sp < snp + step) // Loop checking << 335 while(sp<snp+step) // Loop checking, 06.08.2015, G.Cosmo
355 { 336 {
356 temp += *mit2 * (*(sp++)); << 337 temp+=*mit2*(*(sp++));
357 mit2 += mat2.num_col(); << 338 mit2+=mat2.num_col();
358 } 339 }
359 sp += step - 1; << 340 sp+=step-1;
360 for(stept = step + 1; stept <= mat1.num_ << 341 for(stept=step+1;stept<=mat1.num_row();stept++)
361 { 342 {
362 temp += *mit2 * (*sp); << 343 temp+=*mit2*(*sp);
363 mit2 += mat2.num_col(); << 344 mit2+=mat2.num_col();
364 sp += stept; << 345 sp+=stept;
365 } 346 }
366 *(mir++) = temp; << 347 *(mir++)=temp;
367 } 348 }
368 } 349 }
369 return mret; 350 return mret;
370 } 351 }
371 352
372 G4ErrorMatrix operator*(const G4ErrorSymMatrix << 353 G4ErrorMatrix operator*(const G4ErrorSymMatrix &mat1, const G4ErrorSymMatrix &mat2)
373 const G4ErrorSymMatrix <<
374 { 354 {
375 G4ErrorMatrix mret(mat1.num_row(), mat1.num_ << 355 G4ErrorMatrix mret(mat1.num_row(),mat1.num_row());
376 CHK_DIM_1(mat1.num_col(), mat2.num_row(), *) << 356 CHK_DIM_1(mat1.num_col(),mat2.num_row(),*);
377 G4int step1, stept1, step2, stept2; << 357 G4int step1,stept1,step2,stept2;
378 G4ErrorMatrixConstIter snp1, sp1, snp2, sp2; << 358 G4ErrorMatrixConstIter snp1,sp1,snp2,sp2;
379 G4double temp; 359 G4double temp;
380 G4ErrorMatrixIter mr = mret.m.begin(); 360 G4ErrorMatrixIter mr = mret.m.begin();
381 for(step1 = 1, snp1 = mat1.m.begin(); step1 << 361 for(step1=1,snp1=mat1.m.begin();step1<=mat1.num_row();snp1+=step1++)
382 snp1 += step1++) <<
383 { 362 {
384 for(step2 = 1, snp2 = mat2.m.begin(); step << 363 for(step2=1,snp2=mat2.m.begin();step2<=mat2.num_row();)
385 { 364 {
386 sp1 = snp1; << 365 sp1=snp1;
387 sp2 = snp2; << 366 sp2=snp2;
388 snp2 += step2; << 367 snp2+=step2;
389 temp = 0; << 368 temp=0;
390 if(step1 < step2) << 369 if(step1<step2)
391 { << 370 {
392 while(sp1 < snp1 + step1) // Loop che << 371 while(sp1<snp1+step1) // Loop checking, 06.08.2015, G.Cosmo
393 { << 372 { temp+=(*(sp1++))*(*(sp2++)); }
394 temp += (*(sp1++)) * (*(sp2++)); << 373 sp1+=step1-1;
395 } << 374 for(stept1=step1+1;stept1!=step2+1;sp1+=stept1++)
396 sp1 += step1 - 1; << 375 { temp+=(*sp1)*(*(sp2++)); }
397 for(stept1 = step1 + 1; stept1 != step << 376 sp2+=step2-1;
398 { << 377 for(stept2=++step2;stept2<=mat2.num_row();sp1+=stept1++,sp2+=stept2++)
399 temp += (*sp1) * (*(sp2++)); << 378 { temp+=(*sp1)*(*sp2); }
400 } <<
401 sp2 += step2 - 1; <<
402 for(stept2 = ++step2; stept2 <= mat2.n <<
403 sp1 += stept1++, sp2 += stept2++) <<
404 { <<
405 temp += (*sp1) * (*sp2); <<
406 } <<
407 } 379 }
408 else 380 else
409 { 381 {
410 while(sp2 < snp2) // Loop checking, 0 << 382 while(sp2<snp2) // Loop checking, 06.08.2015, G.Cosmo
411 { << 383 { temp+=(*(sp1++))*(*(sp2++)); }
412 temp += (*(sp1++)) * (*(sp2++)); << 384 sp2+=step2-1;
413 } << 385 for(stept2=++step2;stept2!=step1+1;sp2+=stept2++)
414 sp2 += step2 - 1; << 386 { temp+=(*(sp1++))*(*sp2); }
415 for(stept2 = ++step2; stept2 != step1 << 387 sp1+=step1-1;
416 { << 388 for(stept1=step1+1;stept1<=mat1.num_row();sp1+=stept1++,sp2+=stept2++)
417 temp += (*(sp1++)) * (*sp2); << 389 { temp+=(*sp1)*(*sp2); }
418 } <<
419 sp1 += step1 - 1; <<
420 for(stept1 = step1 + 1; stept1 <= mat1 <<
421 sp1 += stept1++, sp2 += stept2++) <<
422 { <<
423 temp += (*sp1) * (*sp2); <<
424 } <<
425 } 390 }
426 *(mr++) = temp; << 391 *(mr++)=temp;
427 } 392 }
428 } 393 }
429 return mret; 394 return mret;
430 } 395 }
431 396
432 /* ------------------------------------------- 397 /* -----------------------------------------------------------------------
433 This section contains the assignment and in 398 This section contains the assignment and inplace operators =,+=,-=,*=,/=.
434 ------------------------------------------- 399 ----------------------------------------------------------------------- */
435 400
436 G4ErrorMatrix& G4ErrorMatrix::operator+=(const << 401 G4ErrorMatrix & G4ErrorMatrix::operator+=(const G4ErrorSymMatrix &mat2)
437 { 402 {
438 CHK_DIM_2(num_row(), mat2.num_row(), num_col << 403 CHK_DIM_2(num_row(),mat2.num_row(),num_col(),mat2.num_col(),+=);
439 G4int n = num_col(); << 404 G4int n = num_col();
440 G4ErrorMatrixConstIter sjk = mat2.m.begin(); 405 G4ErrorMatrixConstIter sjk = mat2.m.begin();
441 G4ErrorMatrixIter m1j = m.begin(); << 406 G4ErrorMatrixIter m1j = m.begin();
442 G4ErrorMatrixIter mj = m.begin(); << 407 G4ErrorMatrixIter mj = m.begin();
443 // j >= k 408 // j >= k
444 for(G4int j = 1; j <= num_row(); j++) << 409 for(G4int j=1;j<=num_row();j++)
445 { 410 {
446 G4ErrorMatrixIter mjk = mj; 411 G4ErrorMatrixIter mjk = mj;
447 G4ErrorMatrixIter mkj = m1j; 412 G4ErrorMatrixIter mkj = m1j;
448 for(G4int k = 1; k <= j; k++) << 413 for(G4int k=1;k<=j;k++)
449 { 414 {
450 *(mjk++) += *sjk; 415 *(mjk++) += *sjk;
451 if(j != k) << 416 if(j!=k) *mkj += *sjk;
452 *mkj += *sjk; <<
453 sjk++; 417 sjk++;
454 mkj += n; 418 mkj += n;
455 } 419 }
456 mj += n; 420 mj += n;
457 m1j++; 421 m1j++;
458 } 422 }
459 return (*this); 423 return (*this);
460 } 424 }
461 425
462 G4ErrorSymMatrix& G4ErrorSymMatrix::operator+= << 426 G4ErrorSymMatrix & G4ErrorSymMatrix::operator+=(const G4ErrorSymMatrix &mat2)
463 { 427 {
464 CHK_DIM_2(num_row(), mat2.num_row(), num_col << 428 CHK_DIM_2(num_row(),mat2.num_row(),num_col(),mat2.num_col(),+=);
465 SIMPLE_BOP(+=) 429 SIMPLE_BOP(+=)
466 return (*this); 430 return (*this);
467 } 431 }
468 432
469 G4ErrorMatrix& G4ErrorMatrix::operator-=(const << 433 G4ErrorMatrix & G4ErrorMatrix::operator-=(const G4ErrorSymMatrix &mat2)
470 { 434 {
471 CHK_DIM_2(num_row(), mat2.num_row(), num_col << 435 CHK_DIM_2(num_row(),mat2.num_row(),num_col(),mat2.num_col(),-=);
472 G4int n = num_col(); << 436 G4int n = num_col();
473 G4ErrorMatrixConstIter sjk = mat2.m.begin(); 437 G4ErrorMatrixConstIter sjk = mat2.m.begin();
474 G4ErrorMatrixIter m1j = m.begin(); << 438 G4ErrorMatrixIter m1j = m.begin();
475 G4ErrorMatrixIter mj = m.begin(); << 439 G4ErrorMatrixIter mj = m.begin();
476 // j >= k 440 // j >= k
477 for(G4int j = 1; j <= num_row(); j++) << 441 for(G4int j=1;j<=num_row();j++)
478 { 442 {
479 G4ErrorMatrixIter mjk = mj; 443 G4ErrorMatrixIter mjk = mj;
480 G4ErrorMatrixIter mkj = m1j; 444 G4ErrorMatrixIter mkj = m1j;
481 for(G4int k = 1; k <= j; k++) << 445 for(G4int k=1;k<=j;k++)
482 { 446 {
483 *(mjk++) -= *sjk; 447 *(mjk++) -= *sjk;
484 if(j != k) << 448 if(j!=k) *mkj -= *sjk;
485 *mkj -= *sjk; <<
486 sjk++; 449 sjk++;
487 mkj += n; 450 mkj += n;
488 } 451 }
489 mj += n; 452 mj += n;
490 m1j++; 453 m1j++;
491 } 454 }
492 return (*this); 455 return (*this);
493 } 456 }
494 457
495 G4ErrorSymMatrix& G4ErrorSymMatrix::operator-= << 458 G4ErrorSymMatrix & G4ErrorSymMatrix::operator-=(const G4ErrorSymMatrix &mat2)
496 { 459 {
497 CHK_DIM_2(num_row(), mat2.num_row(), num_col << 460 CHK_DIM_2(num_row(),mat2.num_row(),num_col(),mat2.num_col(),-=);
498 SIMPLE_BOP(-=) 461 SIMPLE_BOP(-=)
499 return (*this); 462 return (*this);
500 } 463 }
501 464
502 G4ErrorSymMatrix& G4ErrorSymMatrix::operator/= << 465 G4ErrorSymMatrix & G4ErrorSymMatrix::operator/=(G4double t)
503 { 466 {
504 SIMPLE_UOP(/=) 467 SIMPLE_UOP(/=)
505 return (*this); 468 return (*this);
506 } 469 }
507 470
508 G4ErrorSymMatrix& G4ErrorSymMatrix::operator*= << 471 G4ErrorSymMatrix & G4ErrorSymMatrix::operator*=(G4double t)
509 { 472 {
510 SIMPLE_UOP(*=) 473 SIMPLE_UOP(*=)
511 return (*this); 474 return (*this);
512 } 475 }
513 476
514 G4ErrorMatrix& G4ErrorMatrix::operator=(const << 477 G4ErrorMatrix & G4ErrorMatrix::operator=(const G4ErrorSymMatrix &mat1)
515 { 478 {
516 if(mat1.nrow * mat1.nrow != size) << 479 if(mat1.nrow*mat1.nrow != size)
517 { << 480 {
518 size = mat1.nrow * mat1.nrow; << 481 size = mat1.nrow * mat1.nrow;
519 m.resize(size); << 482 m.resize(size);
520 } << 483 }
521 nrow = mat1.nrow; << 484 nrow = mat1.nrow;
522 ncol = mat1.nrow; << 485 ncol = mat1.nrow;
523 G4int n = ncol; << 486 G4int n = ncol;
524 G4ErrorMatrixConstIter sjk = mat1.m.begin(); << 487 G4ErrorMatrixConstIter sjk = mat1.m.begin();
525 G4ErrorMatrixIter m1j = m.begin(); << 488 G4ErrorMatrixIter m1j = m.begin();
526 G4ErrorMatrixIter mj = m.begin(); << 489 G4ErrorMatrixIter mj = m.begin();
527 // j >= k << 490 // j >= k
528 for(G4int j = 1; j <= num_row(); j++) << 491 for(G4int j=1;j<=num_row();j++)
529 { << 492 {
530 G4ErrorMatrixIter mjk = mj; << 493 G4ErrorMatrixIter mjk = mj;
531 G4ErrorMatrixIter mkj = m1j; << 494 G4ErrorMatrixIter mkj = m1j;
532 for(G4int k = 1; k <= j; k++) << 495 for(G4int k=1;k<=j;k++)
533 { << 496 {
534 *(mjk++) = *sjk; << 497 *(mjk++) = *sjk;
535 if(j != k) << 498 if(j!=k) *mkj = *sjk;
536 *mkj = *sjk; << 499 sjk++;
537 sjk++; << 500 mkj += n;
538 mkj += n; << 501 }
539 } << 502 mj += n;
540 mj += n; << 503 m1j++;
541 m1j++; << 504 }
542 } << 505 return (*this);
543 return (*this); << 506 }
544 } << 507
545 << 508 G4ErrorSymMatrix & G4ErrorSymMatrix::operator=(const G4ErrorSymMatrix &mat1)
546 G4ErrorSymMatrix& G4ErrorSymMatrix::operator=( << 509 {
547 { << 510 if (&mat1 == this) { return *this; }
548 if(&mat1 == this) << 511 if(mat1.nrow != nrow)
549 { << 512 {
550 return *this; << 513 nrow = mat1.nrow;
551 } << 514 size = mat1.size;
552 if(mat1.nrow != nrow) << 515 m.resize(size);
553 { << 516 }
554 nrow = mat1.nrow; << 517 m = mat1.m;
555 size = mat1.size; << 518 return (*this);
556 m.resize(size); <<
557 } <<
558 m = mat1.m; <<
559 return (*this); <<
560 } 519 }
561 520
562 // Print the Matrix. 521 // Print the Matrix.
563 522
564 std::ostream& operator<<(std::ostream& os, con << 523 std::ostream& operator<<(std::ostream &os, const G4ErrorSymMatrix &q)
565 { 524 {
566 os << G4endl; 525 os << G4endl;
567 526
568 // Fixed format needs 3 extra characters for 527 // Fixed format needs 3 extra characters for field,
569 // while scientific needs 7 528 // while scientific needs 7
570 529
571 std::size_t width; << 530 G4int width;
572 if(os.flags() & std::ios::fixed) 531 if(os.flags() & std::ios::fixed)
573 { 532 {
574 width = os.precision() + 3; << 533 width = os.precision()+3;
575 } 534 }
576 else 535 else
577 { 536 {
578 width = os.precision() + 7; << 537 width = os.precision()+7;
579 } 538 }
580 for(G4int irow = 1; irow <= q.num_row(); ++i << 539 for(G4int irow = 1; irow<= q.num_row(); irow++)
581 { 540 {
582 for(G4int icol = 1; icol <= q.num_col(); + << 541 for(G4int icol = 1; icol <= q.num_col(); icol++)
583 { 542 {
584 os.width(width); 543 os.width(width);
585 os << q(irow, icol) << " "; << 544 os << q(irow,icol) << " ";
586 } 545 }
587 os << G4endl; 546 os << G4endl;
588 } 547 }
589 return os; 548 return os;
590 } 549 }
591 550
592 G4ErrorSymMatrix G4ErrorSymMatrix::apply(G4dou << 551 G4ErrorSymMatrix G4ErrorSymMatrix::
593 << 552 apply(G4double (*f)(G4double, G4int, G4int)) const
594 { 553 {
595 G4ErrorSymMatrix mret(num_row()); 554 G4ErrorSymMatrix mret(num_row());
596 G4ErrorMatrixConstIter a = m.cbegin(); << 555 G4ErrorMatrixConstIter a = m.begin();
597 G4ErrorMatrixIter b = mret.m.begin(); << 556 G4ErrorMatrixIter b = mret.m.begin();
598 for(G4int ir = 1; ir <= num_row(); ++ir) << 557 for(G4int ir=1;ir<=num_row();ir++)
599 { 558 {
600 for(G4int ic = 1; ic <= ir; ++ic) << 559 for(G4int ic=1;ic<=ir;ic++)
601 { 560 {
602 *(b++) = (*f)(*(a++), ir, ic); 561 *(b++) = (*f)(*(a++), ir, ic);
603 } 562 }
604 } 563 }
605 return mret; 564 return mret;
606 } 565 }
607 566
608 void G4ErrorSymMatrix::assign(const G4ErrorMat << 567 void G4ErrorSymMatrix::assign (const G4ErrorMatrix &mat1)
609 { 568 {
610 if(mat1.nrow != nrow) << 569 if(mat1.nrow != nrow)
611 { << 570 {
612 nrow = mat1.nrow; << 571 nrow = mat1.nrow;
613 size = nrow * (nrow + 1) / 2; << 572 size = nrow * (nrow+1) / 2;
614 m.resize(size); << 573 m.resize(size);
615 } << 574 }
616 G4ErrorMatrixConstIter a = mat1.m.begin(); << 575 G4ErrorMatrixConstIter a = mat1.m.begin();
617 G4ErrorMatrixIter b = m.begin(); << 576 G4ErrorMatrixIter b = m.begin();
618 for(G4int r = 1; r <= nrow; r++) << 577 for(G4int r=1;r<=nrow;r++)
619 { << 578 {
620 G4ErrorMatrixConstIter d = a; << 579 G4ErrorMatrixConstIter d = a;
621 for(G4int c = 1; c <= r; c++) << 580 for(G4int c=1;c<=r;c++)
622 { << 581 {
623 *(b++) = *(d++); << 582 *(b++) = *(d++);
624 } << 583 }
625 a += nrow; << 584 a += nrow;
626 } << 585 }
627 } 586 }
628 587
629 G4ErrorSymMatrix G4ErrorSymMatrix::similarity( << 588 G4ErrorSymMatrix G4ErrorSymMatrix::similarity(const G4ErrorMatrix &mat1) const
630 { 589 {
631 G4ErrorSymMatrix mret(mat1.num_row()); 590 G4ErrorSymMatrix mret(mat1.num_row());
632 G4ErrorMatrix temp = mat1 * (*this); << 591 G4ErrorMatrix temp = mat1*(*this);
633 592
634 // If mat1*(*this) has correct dimensions, t << 593 // If mat1*(*this) has correct dimensions, then so will the mat1.T multiplication.
635 // multiplication. So there is no need to ch << 594 // So there is no need to check dimensions again.
636 595
637 G4int n = mat1.num_col(); << 596 G4int n = mat1.num_col();
638 G4ErrorMatrixIter mr = mret.m.begin(); << 597 G4ErrorMatrixIter mr = mret.m.begin();
639 G4ErrorMatrixIter tempr1 = temp.m.begin(); 598 G4ErrorMatrixIter tempr1 = temp.m.begin();
640 for(G4int r = 1; r <= mret.num_row(); r++) << 599 for(G4int r=1;r<=mret.num_row();r++)
641 { 600 {
642 G4ErrorMatrixConstIter m1c1 = mat1.m.begin 601 G4ErrorMatrixConstIter m1c1 = mat1.m.begin();
643 for(G4int c = 1; c <= r; c++) << 602 for(G4int c=1;c<=r;c++)
644 { 603 {
645 G4double tmp = 0.0; << 604 G4double tmp = 0.0;
646 G4ErrorMatrixIter tempri = tempr1; << 605 G4ErrorMatrixIter tempri = tempr1;
647 G4ErrorMatrixConstIter m1ci = m1c1; 606 G4ErrorMatrixConstIter m1ci = m1c1;
648 for(G4int i = 1; i <= mat1.num_col(); i+ << 607 for(G4int i=1;i<=mat1.num_col();i++)
649 { 608 {
650 tmp += (*(tempri++)) * (*(m1ci++)); << 609 tmp+=(*(tempri++))*(*(m1ci++));
651 } 610 }
652 *(mr++) = tmp; 611 *(mr++) = tmp;
653 m1c1 += n; 612 m1c1 += n;
654 } 613 }
655 tempr1 += n; 614 tempr1 += n;
656 } 615 }
657 return mret; 616 return mret;
658 } 617 }
659 618
660 G4ErrorSymMatrix G4ErrorSymMatrix::similarity( << 619 G4ErrorSymMatrix G4ErrorSymMatrix::similarity(const G4ErrorSymMatrix &mat1) const
661 const G4ErrorSymMatrix& mat1) const <<
662 { 620 {
663 G4ErrorSymMatrix mret(mat1.num_row()); 621 G4ErrorSymMatrix mret(mat1.num_row());
664 G4ErrorMatrix temp = mat1 * (*this); << 622 G4ErrorMatrix temp = mat1*(*this);
665 G4int n = mat1.num_col(); << 623 G4int n = mat1.num_col();
666 G4ErrorMatrixIter mr = mret.m.begin(); << 624 G4ErrorMatrixIter mr = mret.m.begin();
667 G4ErrorMatrixIter tempr1 = temp.m.begin(); 625 G4ErrorMatrixIter tempr1 = temp.m.begin();
668 for(G4int r = 1; r <= mret.num_row(); r++) << 626 for(G4int r=1;r<=mret.num_row();r++)
669 { 627 {
670 G4ErrorMatrixConstIter m1c1 = mat1.m.begin 628 G4ErrorMatrixConstIter m1c1 = mat1.m.begin();
671 G4int c; 629 G4int c;
672 for(c = 1; c <= r; c++) << 630 for(c=1;c<=r;c++)
673 { 631 {
674 G4double tmp = 0.0; << 632 G4double tmp = 0.0;
675 G4ErrorMatrixIter tempri = tempr1; << 633 G4ErrorMatrixIter tempri = tempr1;
676 G4ErrorMatrixConstIter m1ci = m1c1; 634 G4ErrorMatrixConstIter m1ci = m1c1;
677 G4int i; 635 G4int i;
678 for(i = 1; i < c; i++) << 636 for(i=1;i<c;i++)
679 { 637 {
680 tmp += (*(tempri++)) * (*(m1ci++)); << 638 tmp+=(*(tempri++))*(*(m1ci++));
681 } 639 }
682 for(i = c; i <= mat1.num_col(); i++) << 640 for(i=c;i<=mat1.num_col();i++)
683 { 641 {
684 tmp += (*(tempri++)) * (*(m1ci)); << 642 tmp+=(*(tempri++))*(*(m1ci));
685 m1ci += i; 643 m1ci += i;
686 } 644 }
687 *(mr++) = tmp; 645 *(mr++) = tmp;
688 m1c1 += c; 646 m1c1 += c;
689 } 647 }
690 tempr1 += n; 648 tempr1 += n;
691 } 649 }
692 return mret; 650 return mret;
693 } 651 }
694 652
695 G4ErrorSymMatrix G4ErrorSymMatrix::similarityT << 653 G4ErrorSymMatrix G4ErrorSymMatrix::similarityT(const G4ErrorMatrix &mat1) const
696 { 654 {
697 G4ErrorSymMatrix mret(mat1.num_col()); 655 G4ErrorSymMatrix mret(mat1.num_col());
698 G4ErrorMatrix temp = (*this) * mat1; << 656 G4ErrorMatrix temp = (*this)*mat1;
699 G4int n = mat1.num_col(); << 657 G4int n = mat1.num_col();
700 G4ErrorMatrixIter mrc = mret.m.begin(); << 658 G4ErrorMatrixIter mrc = mret.m.begin();
701 G4ErrorMatrixIter temp1r = temp.m.begin(); 659 G4ErrorMatrixIter temp1r = temp.m.begin();
702 for(G4int r = 1; r <= mret.num_row(); r++) << 660 for(G4int r=1;r<=mret.num_row();r++)
703 { 661 {
704 G4ErrorMatrixConstIter m11c = mat1.m.begin 662 G4ErrorMatrixConstIter m11c = mat1.m.begin();
705 for(G4int c = 1; c <= r; c++) << 663 for(G4int c=1;c<=r;c++)
706 { 664 {
707 G4double tmp = 0.0; << 665 G4double tmp = 0.0;
708 G4ErrorMatrixIter tempir = temp1r; << 666 G4ErrorMatrixIter tempir = temp1r;
709 G4ErrorMatrixConstIter m1ic = m11c; 667 G4ErrorMatrixConstIter m1ic = m11c;
710 for(G4int i = 1; i <= mat1.num_row(); i+ << 668 for(G4int i=1;i<=mat1.num_row();i++)
711 { 669 {
712 tmp += (*(tempir)) * (*(m1ic)); << 670 tmp+=(*(tempir))*(*(m1ic));
713 tempir += n; 671 tempir += n;
714 m1ic += n; 672 m1ic += n;
715 } 673 }
716 *(mrc++) = tmp; 674 *(mrc++) = tmp;
717 m11c++; 675 m11c++;
718 } 676 }
719 temp1r++; 677 temp1r++;
720 } 678 }
721 return mret; 679 return mret;
722 } 680 }
723 681
724 void G4ErrorSymMatrix::invert(G4int& ifail) << 682 void G4ErrorSymMatrix::invert(G4int &ifail)
725 { << 683 {
726 ifail = 0; 684 ifail = 0;
727 685
728 switch(nrow) 686 switch(nrow)
729 { 687 {
730 case 3: << 688 case 3:
731 { 689 {
732 G4double det, temp; 690 G4double det, temp;
733 G4double t1, t2, t3; 691 G4double t1, t2, t3;
734 G4double c11, c12, c13, c22, c23, c33; << 692 G4double c11,c12,c13,c22,c23,c33;
735 c11 = (*(m.begin() + 2)) * (*(m.begin() << 693 c11 = (*(m.begin()+2)) * (*(m.begin()+5))
736 (*(m.begin() + 4)) * (*(m.begin() << 694 - (*(m.begin()+4)) * (*(m.begin()+4));
737 c12 = (*(m.begin() + 4)) * (*(m.begin() << 695 c12 = (*(m.begin()+4)) * (*(m.begin()+3))
738 (*(m.begin() + 1)) * (*(m.begin() << 696 - (*(m.begin()+1)) * (*(m.begin()+5));
739 c13 = (*(m.begin() + 1)) * (*(m.begin() << 697 c13 = (*(m.begin()+1)) * (*(m.begin()+4))
740 (*(m.begin() + 2)) * (*(m.begin() << 698 - (*(m.begin()+2)) * (*(m.begin()+3));
741 c22 = (*(m.begin() + 5)) * (*m.begin()) << 699 c22 = (*(m.begin()+5)) * (*m.begin())
742 (*(m.begin() + 3)) * (*(m.begin() << 700 - (*(m.begin()+3)) * (*(m.begin()+3));
743 c23 = (*(m.begin() + 3)) * (*(m.begin() << 701 c23 = (*(m.begin()+3)) * (*(m.begin()+1))
744 (*(m.begin() + 4)) * (*m.begin()); << 702 - (*(m.begin()+4)) * (*m.begin());
745 c33 = (*m.begin()) * (*(m.begin() + 2)) << 703 c33 = (*m.begin()) * (*(m.begin()+2))
746 (*(m.begin() + 1)) * (*(m.begin() << 704 - (*(m.begin()+1)) * (*(m.begin()+1));
747 t1 = std::fabs(*m.begin()); 705 t1 = std::fabs(*m.begin());
748 t2 = std::fabs(*(m.begin() + 1)); << 706 t2 = std::fabs(*(m.begin()+1));
749 t3 = std::fabs(*(m.begin() + 3)); << 707 t3 = std::fabs(*(m.begin()+3));
750 if(t1 >= t2) << 708 if (t1 >= t2)
751 { 709 {
752 if(t3 >= t1) << 710 if (t3 >= t1)
753 { 711 {
754 temp = *(m.begin() + 3); << 712 temp = *(m.begin()+3);
755 det = c23 * c12 - c22 * c13; << 713 det = c23*c12-c22*c13;
756 } 714 }
757 else 715 else
758 { 716 {
759 temp = *m.begin(); 717 temp = *m.begin();
760 det = c22 * c33 - c23 * c23; << 718 det = c22*c33-c23*c23;
761 } 719 }
762 } 720 }
763 else if(t3 >= t2) << 721 else if (t3 >= t2)
764 { 722 {
765 temp = *(m.begin() + 3); << 723 temp = *(m.begin()+3);
766 det = c23 * c12 - c22 * c13; << 724 det = c23*c12-c22*c13;
767 } 725 }
768 else 726 else
769 { 727 {
770 temp = *(m.begin() + 1); << 728 temp = *(m.begin()+1);
771 det = c13 * c23 - c12 * c33; << 729 det = c13*c23-c12*c33;
772 } 730 }
773 if(det == 0) << 731 if (det==0)
774 { 732 {
775 ifail = 1; 733 ifail = 1;
776 return; 734 return;
777 } 735 }
778 { 736 {
779 G4double ss = temp / det; << 737 G4double ss = temp/det;
780 G4ErrorMatrixIter mq = m.begin(); 738 G4ErrorMatrixIter mq = m.begin();
781 *(mq++) = ss * c11; << 739 *(mq++) = ss*c11;
782 *(mq++) = ss * c12; << 740 *(mq++) = ss*c12;
783 *(mq++) = ss * c22; << 741 *(mq++) = ss*c22;
784 *(mq++) = ss * c13; << 742 *(mq++) = ss*c13;
785 *(mq++) = ss * c23; << 743 *(mq++) = ss*c23;
786 *(mq) = ss * c33; << 744 *(mq) = ss*c33;
787 } 745 }
788 } 746 }
789 break; 747 break;
790 case 2: << 748 case 2:
791 { 749 {
792 G4double det, temp, ss; 750 G4double det, temp, ss;
793 det = (*m.begin()) * (*(m.begin() + 2)) << 751 det = (*m.begin())*(*(m.begin()+2)) - (*(m.begin()+1))*(*(m.begin()+1));
794 (*(m.begin() + 1)) * (*(m.begin() << 752 if (det==0)
795 if(det == 0) <<
796 { 753 {
797 ifail = 1; 754 ifail = 1;
798 return; 755 return;
799 } 756 }
800 ss = 1.0 / det; << 757 ss = 1.0/det;
801 *(m.begin() + 1) *= -ss; << 758 *(m.begin()+1) *= -ss;
802 temp = ss * (*(m.begin() + 2 << 759 temp = ss*(*(m.begin()+2));
803 *(m.begin() + 2) = ss * (*m.begin()); << 760 *(m.begin()+2) = ss*(*m.begin());
804 *m.begin() = temp; << 761 *m.begin() = temp;
805 break; 762 break;
806 } 763 }
807 case 1: << 764 case 1:
808 { 765 {
809 if((*m.begin()) == 0) << 766 if ((*m.begin())==0)
810 { 767 {
811 ifail = 1; 768 ifail = 1;
812 return; 769 return;
813 } 770 }
814 *m.begin() = 1.0 / (*m.begin()); << 771 *m.begin() = 1.0/(*m.begin());
815 break; 772 break;
816 } 773 }
817 case 5: << 774 case 5:
818 { 775 {
819 invert5(ifail); 776 invert5(ifail);
820 return; 777 return;
821 } 778 }
822 case 6: << 779 case 6:
823 { 780 {
824 invert6(ifail); 781 invert6(ifail);
825 return; 782 return;
826 } 783 }
827 case 4: << 784 case 4:
828 { 785 {
829 invert4(ifail); 786 invert4(ifail);
830 return; 787 return;
831 } 788 }
832 default: << 789 default:
833 { 790 {
834 invertBunchKaufman(ifail); 791 invertBunchKaufman(ifail);
835 return; 792 return;
836 } 793 }
837 } 794 }
838 return; // inversion successful << 795 return; // inversion successful
839 } 796 }
840 797
841 G4double G4ErrorSymMatrix::determinant() const 798 G4double G4ErrorSymMatrix::determinant() const
842 { 799 {
843 static const G4int max_array = 20; 800 static const G4int max_array = 20;
844 801
845 // ir must point to an array which is ***1 l 802 // ir must point to an array which is ***1 longer than*** nrow
846 803
847 static std::vector<G4int> ir_vec(max_array + << 804 static std::vector<G4int> ir_vec (max_array+1);
848 if(ir_vec.size() <= static_cast<unsigned int << 805 if (ir_vec.size() <= static_cast<unsigned int>(nrow))
849 { 806 {
850 ir_vec.resize(nrow + 1); << 807 ir_vec.resize(nrow+1);
851 } 808 }
852 G4int* ir = &ir_vec[0]; << 809 G4int * ir = &ir_vec[0];
853 810
854 G4double det; 811 G4double det;
855 G4ErrorMatrix mt(*this); 812 G4ErrorMatrix mt(*this);
856 G4int i = mt.dfact_matrix(det, ir); 813 G4int i = mt.dfact_matrix(det, ir);
857 if(i == 0) << 814 if(i==0) { return det; }
858 { <<
859 return det; <<
860 } <<
861 return 0.0; 815 return 0.0;
862 } 816 }
863 817
864 G4double G4ErrorSymMatrix::trace() const 818 G4double G4ErrorSymMatrix::trace() const
865 { 819 {
866 G4double t = 0.0; << 820 G4double t = 0.0;
867 for(G4int i = 0; i < nrow; i++) << 821 for (G4int i=0; i<nrow; i++)
868 { << 822 { t += *(m.begin() + (i+3)*i/2); }
869 t += *(m.begin() + (i + 3) * i / 2); << 823 return t;
870 } <<
871 return t; <<
872 } 824 }
873 825
874 void G4ErrorSymMatrix::invertBunchKaufman(G4in << 826 void G4ErrorSymMatrix::invertBunchKaufman(G4int &ifail)
875 { 827 {
876 // Bunch-Kaufman diagonal pivoting method 828 // Bunch-Kaufman diagonal pivoting method
877 // It is decribed in J.R. Bunch, L. Kaufman << 829 // It is decribed in J.R. Bunch, L. Kaufman (1977).
878 // "Some Stable Methods for Calculating Iner << 830 // "Some Stable Methods for Calculating Inertia and Solving Symmetric
879 // Linear Systems", Math. Comp. 31, p. 162-1 << 831 // Linear Systems", Math. Comp. 31, p. 162-179. or in Gene H. Golub,
880 // Charles F. van Loan, "Matrix Computations << 832 // Charles F. van Loan, "Matrix Computations" (the second edition
881 // has a bug.) and implemented in "lapack" 833 // has a bug.) and implemented in "lapack"
882 // Mario Stanke, 09/97 834 // Mario Stanke, 09/97
883 835
884 G4int i, j, k, ss; 836 G4int i, j, k, ss;
885 G4int pivrow; 837 G4int pivrow;
886 838
887 // Establish the two working-space arrays ne 839 // Establish the two working-space arrays needed: x and piv are
888 // used as pointers to arrays of doubles and 840 // used as pointers to arrays of doubles and ints respectively, each
889 // of length nrow. We do not want to reallo 841 // of length nrow. We do not want to reallocate each time through
890 // unless the size needs to grow. We do not 842 // unless the size needs to grow. We do not want to leak memory, even
891 // by having a new without a delete that is 843 // by having a new without a delete that is only done once.
892 << 844
893 static const G4int max_array << 845 static const G4int max_array = 25;
894 static G4ThreadLocal std::vector<G4double>* << 846 static G4ThreadLocal std::vector<G4double> *xvec = 0;
895 if(!xvec) << 847 if (!xvec) xvec = new std::vector<G4double> (max_array) ;
896 xvec = new std::vector<G4double>(max_array << 848 static G4ThreadLocal std::vector<G4int> *pivv = 0;
897 static G4ThreadLocal std::vector<G4int>* piv << 849 if (!pivv) pivv = new std::vector<G4int> (max_array) ;
898 if(!pivv) << 850 typedef std::vector<G4int>::iterator pivIter;
899 pivv = new std::vector<G4int>(max_array); << 851 if (xvec->size() < static_cast<unsigned int>(nrow)) xvec->resize(nrow);
900 typedef std::vector<G4int>::iterator pivIter << 852 if (pivv->size() < static_cast<unsigned int>(nrow)) pivv->resize(nrow);
901 if(xvec->size() < static_cast<unsigned int>( << 853 // Note - resize should do nothing if the size is already larger than nrow,
902 xvec->resize(nrow); << 854 // but on VC++ there are indications that it does so we check.
903 if(pivv->size() < static_cast<unsigned int>( << 855 // Note - the data elements in a vector are guaranteed to be contiguous,
904 pivv->resize(nrow); << 856 // so x[i] and piv[i] are optimally fast.
905 // Note - resize should do nothing if the s << 857 G4ErrorMatrixIter x = xvec->begin();
906 // but on VC++ there are indications << 858 // x[i] is used as helper storage, needs to have at least size nrow.
907 // Note - the data elements in a vector are <<
908 // so x[i] and piv[i] are optimally f <<
909 G4ErrorMatrixIter x = xvec->begin(); <<
910 // x[i] is used as helper storage, needs to <<
911 pivIter piv = pivv->begin(); 859 pivIter piv = pivv->begin();
912 // piv[i] is used to store details of exchan << 860 // piv[i] is used to store details of exchanges
913 << 861
914 G4double temp1, temp2; 862 G4double temp1, temp2;
915 G4ErrorMatrixIter ip, mjj, iq; 863 G4ErrorMatrixIter ip, mjj, iq;
916 G4double lambda, sigma; 864 G4double lambda, sigma;
917 const G4double alpha = .6404; // = (1+sqr << 865 const G4double alpha = .6404; // = (1+sqrt(17))/8
918 const G4double epsilon = 32 * DBL_EPSILON; << 866 const G4double epsilon = 32*DBL_EPSILON;
919 // whenever a sum of two doubles is below or << 867 // whenever a sum of two doubles is below or equal to epsilon
920 // it is set to zero. << 868 // it is set to zero.
921 // this constant could be set to zero but th << 869 // this constant could be set to zero but then the algorithm
922 // doesn't neccessarily detect that a matrix << 870 // doesn't neccessarily detect that a matrix is singular
923 << 871
924 for(i = 0; i < nrow; i++) << 872 for (i = 0; i < nrow; i++)
925 { 873 {
926 piv[i] = i + 1; << 874 piv[i] = i+1;
927 } 875 }
928 << 876
929 ifail = 0; 877 ifail = 0;
930 << 878
931 // compute the factorization P*A*P^T = L * D << 879 // compute the factorization P*A*P^T = L * D * L^T
932 // L is unit lower triangular, D is direct s 880 // L is unit lower triangular, D is direct sum of 1x1 and 2x2 matrices
933 // L and D^-1 are stored in A = *this, P is 881 // L and D^-1 are stored in A = *this, P is stored in piv[]
934 << 882
935 for(j = 1; j < nrow; j += ss) // main loop << 883 for (j=1; j < nrow; j+=ss) // main loop over columns
936 { 884 {
937 mjj = m.begin() + j * (j - 1) / 2 + j - << 885 mjj = m.begin() + j*(j-1)/2 + j-1;
938 lambda = 0; // compute lambda = max of A( << 886 lambda = 0; // compute lambda = max of A(j+1:n,j)
939 pivrow = j + 1; << 887 pivrow = j+1;
940 ip = m.begin() + (j + 1) * j / 2 + j - << 888 ip = m.begin() + (j+1)*j/2 + j-1;
941 for(i = j + 1; i <= nrow; ip += i++) << 889 for (i=j+1; i <= nrow ; ip += i++)
942 { <<
943 if(std::fabs(*ip) > lambda) <<
944 { <<
945 lambda = std::fabs(*ip); <<
946 pivrow = i; <<
947 } <<
948 } <<
949 if(lambda == 0) <<
950 { <<
951 if(*mjj == 0) <<
952 { <<
953 ifail = 1; <<
954 return; <<
955 } <<
956 ss = 1; <<
957 *mjj = 1. / *mjj; <<
958 } <<
959 else <<
960 { <<
961 if(std::fabs(*mjj) >= lambda * alpha) <<
962 { <<
963 ss = 1; <<
964 pivrow = j; <<
965 } <<
966 else <<
967 { <<
968 sigma = 0; // compute sigma = max A(p <<
969 ip = m.begin() + pivrow * (pivrow - <<
970 for(k = j; k < pivrow; k++) <<
971 { <<
972 if(std::fabs(*ip) > sigma) <<
973 sigma = std::fabs(*ip); <<
974 ip++; <<
975 } <<
976 if(sigma * std::fabs(*mjj) >= alpha * <<
977 { <<
978 ss = 1; <<
979 pivrow = j; <<
980 } <<
981 else if(std::fabs(*(m.begin() + pivrow <<
982 1)) >= alpha * sig <<
983 { <<
984 ss = 1; <<
985 } <<
986 else <<
987 { <<
988 ss = 2; <<
989 } <<
990 } <<
991 if(pivrow == j) // no permutation necce <<
992 { <<
993 piv[j - 1] = pivrow; <<
994 if(*mjj == 0) <<
995 { <<
996 ifail = 1; <<
997 return; <<
998 } <<
999 temp2 = *mjj = 1. / *mjj; // invert D <<
1000 <<
1001 // update A(j+1:n, j+1,n) <<
1002 for(i = j + 1; i <= nrow; i++) <<
1003 { <<
1004 temp1 = *(m.begin() + i * (i - 1) / <<
1005 ip = m.begin() + i * (i - 1) / 2 <<
1006 for(k = j + 1; k <= i; k++) <<
1007 { 890 {
1008 *ip -= temp1 * *(m.begin() + k * << 891 if (std::fabs(*ip) > lambda)
1009 if(std::fabs(*ip) <= epsilon) <<
1010 { <<
1011 *ip = 0; <<
1012 } <<
1013 ip++; <<
1014 } <<
1015 } <<
1016 // update L <<
1017 ip = m.begin() + (j + 1) * j / 2 + j <<
1018 for(i = j + 1; i <= nrow; ip += i++) <<
1019 { <<
1020 *ip *= temp2; <<
1021 } <<
1022 } <<
1023 else if(ss == 1) // 1x1 pivot <<
1024 { <<
1025 piv[j - 1] = pivrow; <<
1026 <<
1027 // interchange rows and columns j and <<
1028 // submatrix (j:n,j:n) <<
1029 ip = m.begin() + pivrow * (pivrow - 1 <<
1030 for(i = j + 1; i < pivrow; i++, ip++) <<
1031 { <<
1032 temp1 = *(m.begin() + i * (i - 1) / <<
1033 *(m.begin() + i * (i - 1) / 2 + j - <<
1034 *ip <<
1035 } <<
1036 temp1 = *mjj; <<
1037 *mjj = *(m.begin() + pivrow * (pivro <<
1038 *(m.begin() + pivrow * (pivrow - 1) / <<
1039 ip = m.begin() + (pivrow + 1) * pivro <<
1040 iq = ip + pivrow - j; <<
1041 for(i = pivrow + 1; i <= nrow; ip += <<
1042 { <<
1043 temp1 = *iq; <<
1044 *iq = *ip; <<
1045 *ip = temp1; <<
1046 } <<
1047 <<
1048 if(*mjj == 0) <<
1049 { <<
1050 ifail = 1; <<
1051 return; <<
1052 } <<
1053 temp2 = *mjj = 1. / *mjj; // invert <<
1054 <<
1055 // update A(j+1:n, j+1:n) <<
1056 for(i = j + 1; i <= nrow; i++) <<
1057 { <<
1058 temp1 = *(m.begin() + i * (i - 1) / <<
1059 ip = m.begin() + i * (i - 1) / 2 <<
1060 for(k = j + 1; k <= i; k++) <<
1061 { <<
1062 *ip -= temp1 * *(m.begin() + k * <<
1063 if(std::fabs(*ip) <= epsilon) <<
1064 { <<
1065 *ip = 0; <<
1066 } <<
1067 ip++; <<
1068 } <<
1069 } <<
1070 // update L <<
1071 ip = m.begin() + (j + 1) * j / 2 + j <<
1072 for(i = j + 1; i <= nrow; ip += i++) <<
1073 { <<
1074 *ip *= temp2; <<
1075 } <<
1076 } <<
1077 else // ss=2, ie use a 2x2 pivot <<
1078 { <<
1079 piv[j - 1] = -pivrow; <<
1080 piv[j] = 0; // that means this i <<
1081 <<
1082 if(j + 1 != pivrow) <<
1083 { <<
1084 // interchange rows and columns j+1 <<
1085 // submatrix (j:n,j:n) <<
1086 ip = m.begin() + pivrow * (pivrow - <<
1087 for(i = j + 2; i < pivrow; i++, ip+ <<
1088 { <<
1089 temp1 = *(m.begin() + i * (i - 1) <<
1090 *(m.begin() + i * (i - 1) / 2 + j <<
1091 *ip <<
1092 } <<
1093 temp1 = *(mjj + j + 1); <<
1094 *(mjj + j + 1) = <<
1095 *(m.begin() + pivrow * (pivrow - <<
1096 *(m.begin() + pivrow * (pivrow - 1) <<
1097 temp1 <<
1098 *(mjj + j) = *(m.begin() + pivrow * <<
1099 *(m.begin() + pivrow * (pivrow - 1) <<
1100 ip = m.begin() + (pivrow + 1) * piv <<
1101 iq = ip + pivrow - (j + 1); <<
1102 for(i = pivrow + 1; i <= nrow; ip + <<
1103 { <<
1104 temp1 = *iq; <<
1105 *iq = *ip; <<
1106 *ip = temp1; <<
1107 } <<
1108 } <<
1109 // invert D(j:j+1,j:j+1) <<
1110 temp2 = *mjj * *(mjj + j + 1) - *(mjj <<
1111 if(temp2 == 0) <<
1112 { <<
1113 G4Exception("G4ErrorSymMatrix::bunc <<
1114 "GEANT4e-Notification", <<
1115 "Error in pivot choice! <<
1116 } <<
1117 temp2 = 1. / temp2; <<
1118 <<
1119 // this quotient is guaranteed to exi <<
1120 // of the pivot <<
1121 <<
1122 temp1 = *mjj; <<
1123 *mjj = *(mjj + j + 1) * tem <<
1124 *(mjj + j + 1) = temp1 * temp2; <<
1125 *(mjj + j) = -*(mjj + j) * temp2; <<
1126 <<
1127 if(j < nrow - 1) // otherwise do not <<
1128 { <<
1129 // update A(j+2:n, j+2:n) <<
1130 for(i = j + 2; i <= nrow; i++) <<
1131 { <<
1132 ip = m.begin() + i * (i - 1) / <<
1133 temp1 = *ip * *mjj + *(ip + 1) * <<
1134 if(std::fabs(temp1) <= epsilon) <<
1135 { <<
1136 temp1 = 0; <<
1137 } <<
1138 temp2 = *ip * *(mjj + j) + *(ip + <<
1139 if(std::fabs(temp2) <= epsilon) <<
1140 { <<
1141 temp2 = 0; <<
1142 } <<
1143 for(k = j + 2; k <= i; k++) <<
1144 { <<
1145 ip = m.begin() + i * (i - 1) / <<
1146 iq = m.begin() + k * (k - 1) / <<
1147 *ip -= temp1 * *iq + temp2 * *( <<
1148 if(std::fabs(*ip) <= epsilon) <<
1149 { 892 {
1150 *ip = 0; << 893 lambda = std::fabs(*ip);
>> 894 pivrow = i;
1151 } 895 }
1152 } <<
1153 } 896 }
1154 // update L << 897 if (lambda == 0 )
1155 for(i = j + 2; i <= nrow; i++) <<
1156 { <<
1157 ip = m.begin() + i * (i - 1) / <<
1158 temp1 = *ip * *mjj + *(ip + 1) * <<
1159 if(std::fabs(temp1) <= epsilon) <<
1160 { 898 {
1161 temp1 = 0; << 899 if (*mjj == 0)
>> 900 {
>> 901 ifail = 1;
>> 902 return;
>> 903 }
>> 904 ss=1;
>> 905 *mjj = 1./ *mjj;
1162 } 906 }
1163 *(ip + 1) = *ip * *(mjj + j) + *( << 907 else
1164 if(std::fabs(*(ip + 1)) <= epsilo <<
1165 { 908 {
1166 *(ip + 1) = 0; << 909 if (std::fabs(*mjj) >= lambda*alpha)
>> 910 {
>> 911 ss=1;
>> 912 pivrow=j;
>> 913 }
>> 914 else
>> 915 {
>> 916 sigma = 0; // compute sigma = max A(pivrow, j:pivrow-1)
>> 917 ip = m.begin() + pivrow*(pivrow-1)/2+j-1;
>> 918 for (k=j; k < pivrow; k++)
>> 919 {
>> 920 if (std::fabs(*ip) > sigma)
>> 921 sigma = std::fabs(*ip);
>> 922 ip++;
>> 923 }
>> 924 if (sigma * std::fabs(*mjj) >= alpha * lambda * lambda)
>> 925 {
>> 926 ss=1;
>> 927 pivrow = j;
>> 928 }
>> 929 else if (std::fabs(*(m.begin()+pivrow*(pivrow-1)/2+pivrow-1))
>> 930 >= alpha * sigma)
>> 931 { ss=1; }
>> 932 else
>> 933 { ss=2; }
>> 934 }
>> 935 if (pivrow == j) // no permutation neccessary
>> 936 {
>> 937 piv[j-1] = pivrow;
>> 938 if (*mjj == 0)
>> 939 {
>> 940 ifail=1;
>> 941 return;
>> 942 }
>> 943 temp2 = *mjj = 1./ *mjj; // invert D(j,j)
>> 944
>> 945 // update A(j+1:n, j+1,n)
>> 946 for (i=j+1; i <= nrow; i++)
>> 947 {
>> 948 temp1 = *(m.begin() + i*(i-1)/2 + j-1) * temp2;
>> 949 ip = m.begin()+i*(i-1)/2+j;
>> 950 for (k=j+1; k<=i; k++)
>> 951 {
>> 952 *ip -= temp1 * *(m.begin() + k*(k-1)/2 + j-1);
>> 953 if (std::fabs(*ip) <= epsilon)
>> 954 { *ip=0; }
>> 955 ip++;
>> 956 }
>> 957 }
>> 958 // update L
>> 959 ip = m.begin() + (j+1)*j/2 + j-1;
>> 960 for (i=j+1; i <= nrow; ip += i++)
>> 961 {
>> 962 *ip *= temp2;
>> 963 }
>> 964 }
>> 965 else if (ss==1) // 1x1 pivot
>> 966 {
>> 967 piv[j-1] = pivrow;
>> 968
>> 969 // interchange rows and columns j and pivrow in
>> 970 // submatrix (j:n,j:n)
>> 971 ip = m.begin() + pivrow*(pivrow-1)/2 + j;
>> 972 for (i=j+1; i < pivrow; i++, ip++)
>> 973 {
>> 974 temp1 = *(m.begin() + i*(i-1)/2 + j-1);
>> 975 *(m.begin() + i*(i-1)/2 + j-1)= *ip;
>> 976 *ip = temp1;
>> 977 }
>> 978 temp1 = *mjj;
>> 979 *mjj = *(m.begin()+pivrow*(pivrow-1)/2+pivrow-1);
>> 980 *(m.begin()+pivrow*(pivrow-1)/2+pivrow-1) = temp1;
>> 981 ip = m.begin() + (pivrow+1)*pivrow/2 + j-1;
>> 982 iq = ip + pivrow-j;
>> 983 for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
>> 984 {
>> 985 temp1 = *iq;
>> 986 *iq = *ip;
>> 987 *ip = temp1;
>> 988 }
>> 989
>> 990 if (*mjj == 0)
>> 991 {
>> 992 ifail = 1;
>> 993 return;
>> 994 }
>> 995 temp2 = *mjj = 1./ *mjj; // invert D(j,j)
>> 996
>> 997 // update A(j+1:n, j+1:n)
>> 998 for (i = j+1; i <= nrow; i++)
>> 999 {
>> 1000 temp1 = *(m.begin() + i*(i-1)/2 + j-1) * temp2;
>> 1001 ip = m.begin()+i*(i-1)/2+j;
>> 1002 for (k=j+1; k<=i; k++)
>> 1003 {
>> 1004 *ip -= temp1 * *(m.begin() + k*(k-1)/2 + j-1);
>> 1005 if (std::fabs(*ip) <= epsilon)
>> 1006 { *ip=0; }
>> 1007 ip++;
>> 1008 }
>> 1009 }
>> 1010 // update L
>> 1011 ip = m.begin() + (j+1)*j/2 + j-1;
>> 1012 for (i=j+1; i<=nrow; ip += i++)
>> 1013 {
>> 1014 *ip *= temp2;
>> 1015 }
>> 1016 }
>> 1017 else // ss=2, ie use a 2x2 pivot
>> 1018 {
>> 1019 piv[j-1] = -pivrow;
>> 1020 piv[j] = 0; // that means this is the second row of a 2x2 pivot
>> 1021
>> 1022 if (j+1 != pivrow)
>> 1023 {
>> 1024 // interchange rows and columns j+1 and pivrow in
>> 1025 // submatrix (j:n,j:n)
>> 1026 ip = m.begin() + pivrow*(pivrow-1)/2 + j+1;
>> 1027 for (i=j+2; i < pivrow; i++, ip++)
>> 1028 {
>> 1029 temp1 = *(m.begin() + i*(i-1)/2 + j);
>> 1030 *(m.begin() + i*(i-1)/2 + j) = *ip;
>> 1031 *ip = temp1;
>> 1032 }
>> 1033 temp1 = *(mjj + j + 1);
>> 1034 *(mjj + j + 1) =
>> 1035 *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1);
>> 1036 *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1) = temp1;
>> 1037 temp1 = *(mjj + j);
>> 1038 *(mjj + j) = *(m.begin() + pivrow*(pivrow-1)/2 + j-1);
>> 1039 *(m.begin() + pivrow*(pivrow-1)/2 + j-1) = temp1;
>> 1040 ip = m.begin() + (pivrow+1)*pivrow/2 + j;
>> 1041 iq = ip + pivrow-(j+1);
>> 1042 for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
>> 1043 {
>> 1044 temp1 = *iq;
>> 1045 *iq = *ip;
>> 1046 *ip = temp1;
>> 1047 }
>> 1048 }
>> 1049 // invert D(j:j+1,j:j+1)
>> 1050 temp2 = *mjj * *(mjj + j + 1) - *(mjj + j) * *(mjj + j);
>> 1051 if (temp2 == 0)
>> 1052 {
>> 1053 G4Exception("G4ErrorSymMatrix::bunch_invert()",
>> 1054 "GEANT4e-Notification", JustWarning,
>> 1055 "Error in pivot choice!");
>> 1056 }
>> 1057 temp2 = 1. / temp2;
>> 1058
>> 1059 // this quotient is guaranteed to exist by the choice
>> 1060 // of the pivot
>> 1061
>> 1062 temp1 = *mjj;
>> 1063 *mjj = *(mjj + j + 1) * temp2;
>> 1064 *(mjj + j + 1) = temp1 * temp2;
>> 1065 *(mjj + j) = - *(mjj + j) * temp2;
>> 1066
>> 1067 if (j < nrow-1) // otherwise do nothing
>> 1068 {
>> 1069 // update A(j+2:n, j+2:n)
>> 1070 for (i=j+2; i <= nrow ; i++)
>> 1071 {
>> 1072 ip = m.begin() + i*(i-1)/2 + j-1;
>> 1073 temp1 = *ip * *mjj + *(ip + 1) * *(mjj + j);
>> 1074 if (std::fabs(temp1 ) <= epsilon)
>> 1075 { temp1 = 0; }
>> 1076 temp2 = *ip * *(mjj + j) + *(ip + 1) * *(mjj + j + 1);
>> 1077 if (std::fabs(temp2 ) <= epsilon)
>> 1078 { temp2 = 0; }
>> 1079 for (k = j+2; k <= i ; k++)
>> 1080 {
>> 1081 ip = m.begin() + i*(i-1)/2 + k-1;
>> 1082 iq = m.begin() + k*(k-1)/2 + j-1;
>> 1083 *ip -= temp1 * *iq + temp2 * *(iq+1);
>> 1084 if (std::fabs(*ip) <= epsilon)
>> 1085 { *ip = 0; }
>> 1086 }
>> 1087 }
>> 1088 // update L
>> 1089 for (i=j+2; i <= nrow ; i++)
>> 1090 {
>> 1091 ip = m.begin() + i*(i-1)/2 + j-1;
>> 1092 temp1 = *ip * *mjj + *(ip+1) * *(mjj + j);
>> 1093 if (std::fabs(temp1) <= epsilon)
>> 1094 { temp1 = 0; }
>> 1095 *(ip+1) = *ip * *(mjj + j) + *(ip+1) * *(mjj + j + 1);
>> 1096 if (std::fabs(*(ip+1)) <= epsilon)
>> 1097 { *(ip+1) = 0; }
>> 1098 *ip = temp1;
>> 1099 }
>> 1100 }
>> 1101 }
1167 } 1102 }
1168 *ip = temp1; << 1103 } // end of main loop over columns
1169 } <<
1170 } <<
1171 } <<
1172 } <<
1173 } // end of main loop over columns <<
1174 1104
1175 if(j == nrow) // the the last pivot is 1x1 << 1105 if (j == nrow) // the the last pivot is 1x1
1176 { 1106 {
1177 mjj = m.begin() + j * (j - 1) / 2 + j - 1 << 1107 mjj = m.begin() + j*(j-1)/2 + j-1;
1178 if(*mjj == 0) << 1108 if (*mjj == 0)
1179 { 1109 {
1180 ifail = 1; 1110 ifail = 1;
1181 return; 1111 return;
1182 } 1112 }
1183 else 1113 else
1184 { 1114 {
1185 *mjj = 1. / *mjj; 1115 *mjj = 1. / *mjj;
1186 } 1116 }
1187 } // end of last pivot code << 1117 } // end of last pivot code
1188 1118
1189 // computing the inverse from the factoriza 1119 // computing the inverse from the factorization
1190 << 1120
1191 for(j = nrow; j >= 1; j -= ss) // loop ove << 1121 for (j = nrow ; j >= 1 ; j -= ss) // loop over columns
1192 { 1122 {
1193 mjj = m.begin() + j * (j - 1) / 2 + j - 1 << 1123 mjj = m.begin() + j*(j-1)/2 + j-1;
1194 if(piv[j - 1] > 0) // 1x1 pivot, compute << 1124 if (piv[j-1] > 0) // 1x1 pivot, compute column j of inverse
1195 { << 1125 {
1196 ss = 1; << 1126 ss = 1;
1197 if(j < nrow) << 1127 if (j < nrow)
1198 { << 1128 {
1199 ip = m.begin() + (j + 1) * j / 2 + j << 1129 ip = m.begin() + (j+1)*j/2 + j-1;
1200 for(i = 0; i < nrow - j; ip += 1 + j << 1130 for (i=0; i < nrow-j; ip += 1+j+i++)
1201 { << 1131 {
1202 x[i] = *ip; << 1132 x[i] = *ip;
1203 } << 1133 }
1204 for(i = j + 1; i <= nrow; i++) << 1134 for (i=j+1; i<=nrow ; i++)
1205 { << 1135 {
1206 temp2 = 0; << 1136 temp2=0;
1207 ip = m.begin() + i * (i - 1) / 2 << 1137 ip = m.begin() + i*(i-1)/2 + j;
1208 for(k = 0; k <= i - j - 1; k++) << 1138 for (k=0; k <= i-j-1; k++)
1209 { << 1139 { temp2 += *ip++ * x[k]; }
1210 temp2 += *ip++ * x[k]; << 1140 for (ip += i-1; k < nrow-j; ip += 1+j+k++)
1211 } << 1141 { temp2 += *ip * x[k]; }
1212 for(ip += i - 1; k < nrow - j; ip + << 1142 *(m.begin()+ i*(i-1)/2 + j-1) = -temp2;
1213 { << 1143 }
1214 temp2 += *ip * x[k]; << 1144 temp2 = 0;
1215 } << 1145 ip = m.begin() + (j+1)*j/2 + j-1;
1216 *(m.begin() + i * (i - 1) / 2 + j - << 1146 for (k=0; k < nrow-j; ip += 1+j+k++)
1217 } << 1147 { temp2 += x[k] * *ip; }
1218 temp2 = 0; << 1148 *mjj -= temp2;
1219 ip = m.begin() + (j + 1) * j / 2 + << 1149 }
1220 for(k = 0; k < nrow - j; ip += 1 + j << 1150 }
1221 { << 1151 else //2x2 pivot, compute columns j and j-1 of the inverse
1222 temp2 += x[k] * *ip; << 1152 {
1223 } << 1153 if (piv[j-1] != 0)
1224 *mjj -= temp2; << 1154 {
1225 } << 1155 std::ostringstream message;
1226 } << 1156 message << "Error in pivot: " << piv[j-1];
1227 else // 2x2 pivot, compute columns j and << 1157 G4Exception("G4ErrorSymMatrix::invertBunchKaufman()",
1228 { << 1158 "GEANT4e-Notification", JustWarning, message);
1229 if(piv[j - 1] != 0) << 1159 }
1230 { << 1160 ss=2;
1231 std::ostringstream message; << 1161 if (j < nrow)
1232 message << "Error in pivot: " << piv[ << 1162 {
1233 G4Exception("G4ErrorSymMatrix::invert << 1163 ip = m.begin() + (j+1)*j/2 + j-1;
1234 "GEANT4e-Notification", J << 1164 for (i=0; i < nrow-j; ip += 1+j+i++)
1235 } << 1165 {
1236 ss = 2; << 1166 x[i] = *ip;
1237 if(j < nrow) << 1167 }
1238 { << 1168 for (i=j+1; i<=nrow ; i++)
1239 ip = m.begin() + (j + 1) * j / 2 + j << 1169 {
1240 for(i = 0; i < nrow - j; ip += 1 + j << 1170 temp2 = 0;
1241 { << 1171 ip = m.begin() + i*(i-1)/2 + j;
1242 x[i] = *ip; << 1172 for (k=0; k <= i-j-1; k++)
1243 } << 1173 { temp2 += *ip++ * x[k]; }
1244 for(i = j + 1; i <= nrow; i++) << 1174 for (ip += i-1; k < nrow-j; ip += 1+j+k++)
1245 { << 1175 { temp2 += *ip * x[k]; }
1246 temp2 = 0; << 1176 *(m.begin()+ i*(i-1)/2 + j-1) = -temp2;
1247 ip = m.begin() + i * (i - 1) / 2 << 1177 }
1248 for(k = 0; k <= i - j - 1; k++) << 1178 temp2 = 0;
1249 { << 1179 ip = m.begin() + (j+1)*j/2 + j-1;
1250 temp2 += *ip++ * x[k]; << 1180 for (k=0; k < nrow-j; ip += 1+j+k++)
1251 } << 1181 { temp2 += x[k] * *ip; }
1252 for(ip += i - 1; k < nrow - j; ip + << 1182 *mjj -= temp2;
1253 { << 1183 temp2 = 0;
1254 temp2 += *ip * x[k]; << 1184 ip = m.begin() + (j+1)*j/2 + j-2;
1255 } << 1185 for (i=j+1; i <= nrow; ip += i++)
1256 *(m.begin() + i * (i - 1) / 2 + j - << 1186 { temp2 += *ip * *(ip+1); }
1257 } << 1187 *(mjj-1) -= temp2;
1258 temp2 = 0; << 1188 ip = m.begin() + (j+1)*j/2 + j-2;
1259 ip = m.begin() + (j + 1) * j / 2 + << 1189 for (i=0; i < nrow-j; ip += 1+j+i++)
1260 for(k = 0; k < nrow - j; ip += 1 + j << 1190 {
1261 { << 1191 x[i] = *ip;
1262 temp2 += x[k] * *ip; << 1192 }
1263 } << 1193 for (i=j+1; i <= nrow ; i++)
1264 *mjj -= temp2; << 1194 {
1265 temp2 = 0; << 1195 temp2 = 0;
1266 ip = m.begin() + (j + 1) * j / 2 + << 1196 ip = m.begin() + i*(i-1)/2 + j;
1267 for(i = j + 1; i <= nrow; ip += i++) << 1197 for (k=0; k <= i-j-1; k++)
1268 { << 1198 { temp2 += *ip++ * x[k]; }
1269 temp2 += *ip * *(ip + 1); << 1199 for (ip += i-1; k < nrow-j; ip += 1+j+k++)
1270 } << 1200 { temp2 += *ip * x[k]; }
1271 *(mjj - 1) -= temp2; << 1201 *(m.begin()+ i*(i-1)/2 + j-2)= -temp2;
1272 ip = m.begin() + (j + 1) * j / 2 + j << 1202 }
1273 for(i = 0; i < nrow - j; ip += 1 + j << 1203 temp2 = 0;
1274 { << 1204 ip = m.begin() + (j+1)*j/2 + j-2;
1275 x[i] = *ip; << 1205 for (k=0; k < nrow-j; ip += 1+j+k++)
1276 } << 1206 { temp2 += x[k] * *ip; }
1277 for(i = j + 1; i <= nrow; i++) << 1207 *(mjj-j) -= temp2;
1278 { << 1208 }
1279 temp2 = 0; << 1209 }
1280 ip = m.begin() + i * (i - 1) / 2 << 1210
1281 for(k = 0; k <= i - j - 1; k++) << 1211 // interchange rows and columns j and piv[j-1]
1282 { << 1212 // or rows and columns j and -piv[j-2]
1283 temp2 += *ip++ * x[k]; << 1213
1284 } << 1214 pivrow = (piv[j-1]==0)? -piv[j-2] : piv[j-1];
1285 for(ip += i - 1; k < nrow - j; ip + << 1215 ip = m.begin() + pivrow*(pivrow-1)/2 + j;
1286 { << 1216 for (i=j+1;i < pivrow; i++, ip++)
1287 temp2 += *ip * x[k]; << 1217 {
1288 } << 1218 temp1 = *(m.begin() + i*(i-1)/2 + j-1);
1289 *(m.begin() + i * (i - 1) / 2 + j - << 1219 *(m.begin() + i*(i-1)/2 + j-1) = *ip;
1290 } << 1220 *ip = temp1;
1291 temp2 = 0; << 1221 }
1292 ip = m.begin() + (j + 1) * j / 2 + << 1222 temp1 = *mjj;
1293 for(k = 0; k < nrow - j; ip += 1 + j << 1223 *mjj = *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1);
1294 { << 1224 *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1) = temp1;
1295 temp2 += x[k] * *ip; << 1225 if (ss==2)
1296 } << 1226 {
1297 *(mjj - j) -= temp2; << 1227 temp1 = *(mjj-1);
1298 } << 1228 *(mjj-1) = *( m.begin() + pivrow*(pivrow-1)/2 + j-2);
1299 } << 1229 *( m.begin() + pivrow*(pivrow-1)/2 + j-2) = temp1;
1300 << 1230 }
1301 // interchange rows and columns j and piv << 1231
1302 // or rows and columns j and -piv[j-2] << 1232 ip = m.begin() + (pivrow+1)*pivrow/2 + j-1; // &A(i,j)
1303 << 1233 iq = ip + pivrow-j;
1304 pivrow = (piv[j - 1] == 0) ? -piv[j - 2] << 1234 for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
1305 ip = m.begin() + pivrow * (pivrow - 1 << 1235 {
1306 for(i = j + 1; i < pivrow; i++, ip++) << 1236 temp1 = *iq;
1307 { << 1237 *iq = *ip;
1308 temp1 = *(m.begin() + i * (i - 1) / 2 + << 1238 *ip = temp1;
1309 *(m.begin() + i * (i - 1) / 2 + j - 1) << 1239 }
1310 *ip << 1240 } // end of loop over columns (in computing inverse from factorization)
1311 } <<
1312 temp1 = *mjj; <<
1313 *mjj = *(m.begin() + pivrow * (pivrow - <<
1314 *(m.begin() + pivrow * (pivrow - 1) / 2 + <<
1315 if(ss == 2) <<
1316 { <<
1317 temp1 = *(mjj - 1); <<
1318 *(mjj - 1) = *(m.begin() + pivrow * (pi <<
1319 *(m.begin() + pivrow * (pivrow - 1) / 2 <<
1320 } <<
1321 <<
1322 ip = m.begin() + (pivrow + 1) * pivrow / <<
1323 iq = ip + pivrow - j; <<
1324 for(i = pivrow + 1; i <= nrow; ip += i, i <<
1325 { <<
1326 temp1 = *iq; <<
1327 *iq = *ip; <<
1328 *ip = temp1; <<
1329 } <<
1330 } // end of loop over columns (in computin <<
1331 1241
1332 return; // inversion successful << 1242 return; // inversion successful
1333 } 1243 }
1334 1244
1335 G4ThreadLocal G4double G4ErrorSymMatrix::posD 1245 G4ThreadLocal G4double G4ErrorSymMatrix::posDefFraction5x5 = 1.0;
1336 G4ThreadLocal G4double G4ErrorSymMatrix::posD 1246 G4ThreadLocal G4double G4ErrorSymMatrix::posDefFraction6x6 = 1.0;
1337 G4ThreadLocal G4double G4ErrorSymMatrix::adju << 1247 G4ThreadLocal G4double G4ErrorSymMatrix::adjustment5x5 = 0.0;
1338 G4ThreadLocal G4double G4ErrorSymMatrix::adju << 1248 G4ThreadLocal G4double G4ErrorSymMatrix::adjustment6x6 = 0.0;
1339 const G4double G4ErrorSymMatrix::CHOLESKY_THR << 1249 const G4double G4ErrorSymMatrix::CHOLESKY_THRESHOLD_5x5 = .5;
1340 const G4double G4ErrorSymMatrix::CHOLESKY_THR << 1250 const G4double G4ErrorSymMatrix::CHOLESKY_THRESHOLD_6x6 = .2;
1341 const G4double G4ErrorSymMatrix::CHOLESKY_CRE << 1251 const G4double G4ErrorSymMatrix::CHOLESKY_CREEP_5x5 = .005;
1342 const G4double G4ErrorSymMatrix::CHOLESKY_CRE << 1252 const G4double G4ErrorSymMatrix::CHOLESKY_CREEP_6x6 = .002;
1343 1253
1344 // Aij are indices for a 6x6 symmetric matrix 1254 // Aij are indices for a 6x6 symmetric matrix.
1345 // The indices for 5x5 or 4x4 symmetric m << 1255 // The indices for 5x5 or 4x4 symmetric matrices are the same,
1346 // ignoring all combinations with an inde 1256 // ignoring all combinations with an index which is inapplicable.
1347 1257
1348 #define A00 0 1258 #define A00 0
1349 #define A01 1 1259 #define A01 1
1350 #define A02 3 1260 #define A02 3
1351 #define A03 6 1261 #define A03 6
1352 #define A04 10 1262 #define A04 10
1353 #define A05 15 1263 #define A05 15
1354 1264
1355 #define A10 1 1265 #define A10 1
1356 #define A11 2 1266 #define A11 2
1357 #define A12 4 1267 #define A12 4
1358 #define A13 7 1268 #define A13 7
1359 #define A14 11 1269 #define A14 11
1360 #define A15 16 1270 #define A15 16
1361 1271
1362 #define A20 3 1272 #define A20 3
1363 #define A21 4 1273 #define A21 4
1364 #define A22 5 1274 #define A22 5
1365 #define A23 8 1275 #define A23 8
1366 #define A24 12 1276 #define A24 12
1367 #define A25 17 1277 #define A25 17
1368 1278
1369 #define A30 6 1279 #define A30 6
1370 #define A31 7 1280 #define A31 7
1371 #define A32 8 1281 #define A32 8
1372 #define A33 9 1282 #define A33 9
1373 #define A34 13 1283 #define A34 13
1374 #define A35 18 1284 #define A35 18
1375 1285
1376 #define A40 10 1286 #define A40 10
1377 #define A41 11 1287 #define A41 11
1378 #define A42 12 1288 #define A42 12
1379 #define A43 13 1289 #define A43 13
1380 #define A44 14 1290 #define A44 14
1381 #define A45 19 1291 #define A45 19
1382 1292
1383 #define A50 15 1293 #define A50 15
1384 #define A51 16 1294 #define A51 16
1385 #define A52 17 1295 #define A52 17
1386 #define A53 18 1296 #define A53 18
1387 #define A54 19 1297 #define A54 19
1388 #define A55 20 1298 #define A55 20
1389 1299
1390 void G4ErrorSymMatrix::invert5(G4int& ifail) << 1300 void G4ErrorSymMatrix::invert5(G4int & ifail)
1391 { << 1301 {
1392 if(posDefFraction5x5 >= CHOLESKY_THRESHOLD_ << 1302 if (posDefFraction5x5 >= CHOLESKY_THRESHOLD_5x5)
1393 { << 1303 {
1394 invertCholesky5(ifail); << 1304 invertCholesky5(ifail);
1395 posDefFraction5x5 = .9 * posDefFraction5x << 1305 posDefFraction5x5 = .9*posDefFraction5x5 + .1*(1-ifail);
1396 if(ifail != 0) // Cholesky failed -- inv << 1306 if (ifail!=0) // Cholesky failed -- invert using Haywood
1397 { << 1307 {
1398 invertHaywood5(ifail); << 1308 invertHaywood5(ifail);
1399 } << 1309 }
1400 } << 1310 }
1401 else << 1311 else
1402 { <<
1403 if(posDefFraction5x5 + adjustment5x5 >= C <<
1404 { <<
1405 invertCholesky5(ifail); <<
1406 posDefFraction5x5 = .9 * posDefFraction <<
1407 if(ifail != 0) // Cholesky failed -- i <<
1408 { 1312 {
1409 invertHaywood5(ifail); << 1313 if (posDefFraction5x5 + adjustment5x5 >= CHOLESKY_THRESHOLD_5x5)
1410 adjustment5x5 = 0; << 1314 {
>> 1315 invertCholesky5(ifail);
>> 1316 posDefFraction5x5 = .9*posDefFraction5x5 + .1*(1-ifail);
>> 1317 if (ifail!=0) // Cholesky failed -- invert using Haywood
>> 1318 {
>> 1319 invertHaywood5(ifail);
>> 1320 adjustment5x5 = 0;
>> 1321 }
>> 1322 }
>> 1323 else
>> 1324 {
>> 1325 invertHaywood5(ifail);
>> 1326 adjustment5x5 += CHOLESKY_CREEP_5x5;
>> 1327 }
1411 } 1328 }
1412 } << 1329 return;
1413 else <<
1414 { <<
1415 invertHaywood5(ifail); <<
1416 adjustment5x5 += CHOLESKY_CREEP_5x5; <<
1417 } <<
1418 } <<
1419 return; <<
1420 } 1330 }
1421 1331
1422 void G4ErrorSymMatrix::invert6(G4int& ifail) << 1332 void G4ErrorSymMatrix::invert6(G4int & ifail)
1423 { 1333 {
1424 if(posDefFraction6x6 >= CHOLESKY_THRESHOLD_ << 1334 if (posDefFraction6x6 >= CHOLESKY_THRESHOLD_6x6)
1425 { <<
1426 invertCholesky6(ifail); <<
1427 posDefFraction6x6 = .9 * posDefFraction6x <<
1428 if(ifail != 0) // Cholesky failed -- inv <<
1429 { <<
1430 invertHaywood6(ifail); <<
1431 } <<
1432 } <<
1433 else <<
1434 { <<
1435 if(posDefFraction6x6 + adjustment6x6 >= C <<
1436 { <<
1437 invertCholesky6(ifail); <<
1438 posDefFraction6x6 = .9 * posDefFraction <<
1439 if(ifail != 0) // Cholesky failed -- i <<
1440 { 1335 {
1441 invertHaywood6(ifail); << 1336 invertCholesky6(ifail);
1442 adjustment6x6 = 0; << 1337 posDefFraction6x6 = .9*posDefFraction6x6 + .1*(1-ifail);
>> 1338 if (ifail!=0) // Cholesky failed -- invert using Haywood
>> 1339 {
>> 1340 invertHaywood6(ifail);
>> 1341 }
1443 } 1342 }
1444 } << 1343 else
1445 else << 1344 {
1446 { << 1345 if (posDefFraction6x6 + adjustment6x6 >= CHOLESKY_THRESHOLD_6x6)
1447 invertHaywood6(ifail); << 1346 {
1448 adjustment6x6 += CHOLESKY_CREEP_6x6; << 1347 invertCholesky6(ifail);
1449 } << 1348 posDefFraction6x6 = .9*posDefFraction6x6 + .1*(1-ifail);
1450 } << 1349 if (ifail!=0) // Cholesky failed -- invert using Haywood
1451 return; << 1350 {
>> 1351 invertHaywood6(ifail);
>> 1352 adjustment6x6 = 0;
>> 1353 }
>> 1354 }
>> 1355 else
>> 1356 {
>> 1357 invertHaywood6(ifail);
>> 1358 adjustment6x6 += CHOLESKY_CREEP_6x6;
>> 1359 }
>> 1360 }
>> 1361 return;
1452 } 1362 }
1453 1363
1454 void G4ErrorSymMatrix::invertHaywood5(G4int& << 1364 void G4ErrorSymMatrix::invertHaywood5 (G4int & ifail)
1455 { 1365 {
1456 ifail = 0; 1366 ifail = 0;
1457 1367
1458 // Find all NECESSARY 2x2 dets: (25 of the 1368 // Find all NECESSARY 2x2 dets: (25 of them)
1459 1369
1460 G4double Det2_23_01 = m[A20] * m[A31] - m[A << 1370 G4double Det2_23_01 = m[A20]*m[A31] - m[A21]*m[A30];
1461 G4double Det2_23_02 = m[A20] * m[A32] - m[A << 1371 G4double Det2_23_02 = m[A20]*m[A32] - m[A22]*m[A30];
1462 G4double Det2_23_03 = m[A20] * m[A33] - m[A << 1372 G4double Det2_23_03 = m[A20]*m[A33] - m[A23]*m[A30];
1463 G4double Det2_23_12 = m[A21] * m[A32] - m[A << 1373 G4double Det2_23_12 = m[A21]*m[A32] - m[A22]*m[A31];
1464 G4double Det2_23_13 = m[A21] * m[A33] - m[A << 1374 G4double Det2_23_13 = m[A21]*m[A33] - m[A23]*m[A31];
1465 G4double Det2_23_23 = m[A22] * m[A33] - m[A << 1375 G4double Det2_23_23 = m[A22]*m[A33] - m[A23]*m[A32];
1466 G4double Det2_24_01 = m[A20] * m[A41] - m[A << 1376 G4double Det2_24_01 = m[A20]*m[A41] - m[A21]*m[A40];
1467 G4double Det2_24_02 = m[A20] * m[A42] - m[A << 1377 G4double Det2_24_02 = m[A20]*m[A42] - m[A22]*m[A40];
1468 G4double Det2_24_03 = m[A20] * m[A43] - m[A << 1378 G4double Det2_24_03 = m[A20]*m[A43] - m[A23]*m[A40];
1469 G4double Det2_24_04 = m[A20] * m[A44] - m[A << 1379 G4double Det2_24_04 = m[A20]*m[A44] - m[A24]*m[A40];
1470 G4double Det2_24_12 = m[A21] * m[A42] - m[A << 1380 G4double Det2_24_12 = m[A21]*m[A42] - m[A22]*m[A41];
1471 G4double Det2_24_13 = m[A21] * m[A43] - m[A << 1381 G4double Det2_24_13 = m[A21]*m[A43] - m[A23]*m[A41];
1472 G4double Det2_24_14 = m[A21] * m[A44] - m[A << 1382 G4double Det2_24_14 = m[A21]*m[A44] - m[A24]*m[A41];
1473 G4double Det2_24_23 = m[A22] * m[A43] - m[A << 1383 G4double Det2_24_23 = m[A22]*m[A43] - m[A23]*m[A42];
1474 G4double Det2_24_24 = m[A22] * m[A44] - m[A << 1384 G4double Det2_24_24 = m[A22]*m[A44] - m[A24]*m[A42];
1475 G4double Det2_34_01 = m[A30] * m[A41] - m[A << 1385 G4double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40];
1476 G4double Det2_34_02 = m[A30] * m[A42] - m[A << 1386 G4double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40];
1477 G4double Det2_34_03 = m[A30] * m[A43] - m[A << 1387 G4double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40];
1478 G4double Det2_34_04 = m[A30] * m[A44] - m[A << 1388 G4double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40];
1479 G4double Det2_34_12 = m[A31] * m[A42] - m[A << 1389 G4double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41];
1480 G4double Det2_34_13 = m[A31] * m[A43] - m[A << 1390 G4double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41];
1481 G4double Det2_34_14 = m[A31] * m[A44] - m[A << 1391 G4double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41];
1482 G4double Det2_34_23 = m[A32] * m[A43] - m[A << 1392 G4double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42];
1483 G4double Det2_34_24 = m[A32] * m[A44] - m[A << 1393 G4double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42];
1484 G4double Det2_34_34 = m[A33] * m[A44] - m[A << 1394 G4double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43];
1485 1395
1486 // Find all NECESSARY 3x3 dets: (30 of th 1396 // Find all NECESSARY 3x3 dets: (30 of them)
1487 1397
1488 G4double Det3_123_012 = << 1398 G4double Det3_123_012 = m[A10]*Det2_23_12 - m[A11]*Det2_23_02
1489 m[A10] * Det2_23_12 - m[A11] * Det2_23_02 << 1399 + m[A12]*Det2_23_01;
1490 G4double Det3_123_013 = << 1400 G4double Det3_123_013 = m[A10]*Det2_23_13 - m[A11]*Det2_23_03
1491 m[A10] * Det2_23_13 - m[A11] * Det2_23_03 << 1401 + m[A13]*Det2_23_01;
1492 G4double Det3_123_023 = << 1402 G4double Det3_123_023 = m[A10]*Det2_23_23 - m[A12]*Det2_23_03
1493 m[A10] * Det2_23_23 - m[A12] * Det2_23_03 << 1403 + m[A13]*Det2_23_02;
1494 G4double Det3_123_123 = << 1404 G4double Det3_123_123 = m[A11]*Det2_23_23 - m[A12]*Det2_23_13
1495 m[A11] * Det2_23_23 - m[A12] * Det2_23_13 << 1405 + m[A13]*Det2_23_12;
1496 G4double Det3_124_012 = << 1406 G4double Det3_124_012 = m[A10]*Det2_24_12 - m[A11]*Det2_24_02
1497 m[A10] * Det2_24_12 - m[A11] * Det2_24_02 << 1407 + m[A12]*Det2_24_01;
1498 G4double Det3_124_013 = << 1408 G4double Det3_124_013 = m[A10]*Det2_24_13 - m[A11]*Det2_24_03
1499 m[A10] * Det2_24_13 - m[A11] * Det2_24_03 << 1409 + m[A13]*Det2_24_01;
1500 G4double Det3_124_014 = << 1410 G4double Det3_124_014 = m[A10]*Det2_24_14 - m[A11]*Det2_24_04
1501 m[A10] * Det2_24_14 - m[A11] * Det2_24_04 << 1411 + m[A14]*Det2_24_01;
1502 G4double Det3_124_023 = << 1412 G4double Det3_124_023 = m[A10]*Det2_24_23 - m[A12]*Det2_24_03
1503 m[A10] * Det2_24_23 - m[A12] * Det2_24_03 << 1413 + m[A13]*Det2_24_02;
1504 G4double Det3_124_024 = << 1414 G4double Det3_124_024 = m[A10]*Det2_24_24 - m[A12]*Det2_24_04
1505 m[A10] * Det2_24_24 - m[A12] * Det2_24_04 << 1415 + m[A14]*Det2_24_02;
1506 G4double Det3_124_123 = << 1416 G4double Det3_124_123 = m[A11]*Det2_24_23 - m[A12]*Det2_24_13
1507 m[A11] * Det2_24_23 - m[A12] * Det2_24_13 << 1417 + m[A13]*Det2_24_12;
1508 G4double Det3_124_124 = << 1418 G4double Det3_124_124 = m[A11]*Det2_24_24 - m[A12]*Det2_24_14
1509 m[A11] * Det2_24_24 - m[A12] * Det2_24_14 << 1419 + m[A14]*Det2_24_12;
1510 G4double Det3_134_012 = << 1420 G4double Det3_134_012 = m[A10]*Det2_34_12 - m[A11]*Det2_34_02
1511 m[A10] * Det2_34_12 - m[A11] * Det2_34_02 << 1421 + m[A12]*Det2_34_01;
1512 G4double Det3_134_013 = << 1422 G4double Det3_134_013 = m[A10]*Det2_34_13 - m[A11]*Det2_34_03
1513 m[A10] * Det2_34_13 - m[A11] * Det2_34_03 << 1423 + m[A13]*Det2_34_01;
1514 G4double Det3_134_014 = << 1424 G4double Det3_134_014 = m[A10]*Det2_34_14 - m[A11]*Det2_34_04
1515 m[A10] * Det2_34_14 - m[A11] * Det2_34_04 << 1425 + m[A14]*Det2_34_01;
1516 G4double Det3_134_023 = << 1426 G4double Det3_134_023 = m[A10]*Det2_34_23 - m[A12]*Det2_34_03
1517 m[A10] * Det2_34_23 - m[A12] * Det2_34_03 << 1427 + m[A13]*Det2_34_02;
1518 G4double Det3_134_024 = << 1428 G4double Det3_134_024 = m[A10]*Det2_34_24 - m[A12]*Det2_34_04
1519 m[A10] * Det2_34_24 - m[A12] * Det2_34_04 << 1429 + m[A14]*Det2_34_02;
1520 G4double Det3_134_034 = << 1430 G4double Det3_134_034 = m[A10]*Det2_34_34 - m[A13]*Det2_34_04
1521 m[A10] * Det2_34_34 - m[A13] * Det2_34_04 << 1431 + m[A14]*Det2_34_03;
1522 G4double Det3_134_123 = << 1432 G4double Det3_134_123 = m[A11]*Det2_34_23 - m[A12]*Det2_34_13
1523 m[A11] * Det2_34_23 - m[A12] * Det2_34_13 << 1433 + m[A13]*Det2_34_12;
1524 G4double Det3_134_124 = << 1434 G4double Det3_134_124 = m[A11]*Det2_34_24 - m[A12]*Det2_34_14
1525 m[A11] * Det2_34_24 - m[A12] * Det2_34_14 << 1435 + m[A14]*Det2_34_12;
1526 G4double Det3_134_134 = << 1436 G4double Det3_134_134 = m[A11]*Det2_34_34 - m[A13]*Det2_34_14
1527 m[A11] * Det2_34_34 - m[A13] * Det2_34_14 << 1437 + m[A14]*Det2_34_13;
1528 G4double Det3_234_012 = << 1438 G4double Det3_234_012 = m[A20]*Det2_34_12 - m[A21]*Det2_34_02
1529 m[A20] * Det2_34_12 - m[A21] * Det2_34_02 << 1439 + m[A22]*Det2_34_01;
1530 G4double Det3_234_013 = << 1440 G4double Det3_234_013 = m[A20]*Det2_34_13 - m[A21]*Det2_34_03
1531 m[A20] * Det2_34_13 - m[A21] * Det2_34_03 << 1441 + m[A23]*Det2_34_01;
1532 G4double Det3_234_014 = << 1442 G4double Det3_234_014 = m[A20]*Det2_34_14 - m[A21]*Det2_34_04
1533 m[A20] * Det2_34_14 - m[A21] * Det2_34_04 << 1443 + m[A24]*Det2_34_01;
1534 G4double Det3_234_023 = << 1444 G4double Det3_234_023 = m[A20]*Det2_34_23 - m[A22]*Det2_34_03
1535 m[A20] * Det2_34_23 - m[A22] * Det2_34_03 << 1445 + m[A23]*Det2_34_02;
1536 G4double Det3_234_024 = << 1446 G4double Det3_234_024 = m[A20]*Det2_34_24 - m[A22]*Det2_34_04
1537 m[A20] * Det2_34_24 - m[A22] * Det2_34_04 << 1447 + m[A24]*Det2_34_02;
1538 G4double Det3_234_034 = << 1448 G4double Det3_234_034 = m[A20]*Det2_34_34 - m[A23]*Det2_34_04
1539 m[A20] * Det2_34_34 - m[A23] * Det2_34_04 << 1449 + m[A24]*Det2_34_03;
1540 G4double Det3_234_123 = << 1450 G4double Det3_234_123 = m[A21]*Det2_34_23 - m[A22]*Det2_34_13
1541 m[A21] * Det2_34_23 - m[A22] * Det2_34_13 << 1451 + m[A23]*Det2_34_12;
1542 G4double Det3_234_124 = << 1452 G4double Det3_234_124 = m[A21]*Det2_34_24 - m[A22]*Det2_34_14
1543 m[A21] * Det2_34_24 - m[A22] * Det2_34_14 << 1453 + m[A24]*Det2_34_12;
1544 G4double Det3_234_134 = << 1454 G4double Det3_234_134 = m[A21]*Det2_34_34 - m[A23]*Det2_34_14
1545 m[A21] * Det2_34_34 - m[A23] * Det2_34_14 << 1455 + m[A24]*Det2_34_13;
1546 G4double Det3_234_234 = << 1456 G4double Det3_234_234 = m[A22]*Det2_34_34 - m[A23]*Det2_34_24
1547 m[A22] * Det2_34_34 - m[A23] * Det2_34_24 << 1457 + m[A24]*Det2_34_23;
1548 1458
1549 // Find all NECESSARY 4x4 dets: (15 of th 1459 // Find all NECESSARY 4x4 dets: (15 of them)
1550 1460
1551 G4double Det4_0123_0123 = m[A00] * Det3_123 << 1461 G4double Det4_0123_0123 = m[A00]*Det3_123_123 - m[A01]*Det3_123_023
1552 m[A02] * Det3_123 << 1462 + m[A02]*Det3_123_013 - m[A03]*Det3_123_012;
1553 G4double Det4_0124_0123 = m[A00] * Det3_124 << 1463 G4double Det4_0124_0123 = m[A00]*Det3_124_123 - m[A01]*Det3_124_023
1554 m[A02] * Det3_124 << 1464 + m[A02]*Det3_124_013 - m[A03]*Det3_124_012;
1555 G4double Det4_0124_0124 = m[A00] * Det3_124 << 1465 G4double Det4_0124_0124 = m[A00]*Det3_124_124 - m[A01]*Det3_124_024
1556 m[A02] * Det3_124 << 1466 + m[A02]*Det3_124_014 - m[A04]*Det3_124_012;
1557 G4double Det4_0134_0123 = m[A00] * Det3_134 << 1467 G4double Det4_0134_0123 = m[A00]*Det3_134_123 - m[A01]*Det3_134_023
1558 m[A02] * Det3_134 << 1468 + m[A02]*Det3_134_013 - m[A03]*Det3_134_012;
1559 G4double Det4_0134_0124 = m[A00] * Det3_134 << 1469 G4double Det4_0134_0124 = m[A00]*Det3_134_124 - m[A01]*Det3_134_024
1560 m[A02] * Det3_134 << 1470 + m[A02]*Det3_134_014 - m[A04]*Det3_134_012;
1561 G4double Det4_0134_0134 = m[A00] * Det3_134 << 1471 G4double Det4_0134_0134 = m[A00]*Det3_134_134 - m[A01]*Det3_134_034
1562 m[A03] * Det3_134 << 1472 + m[A03]*Det3_134_014 - m[A04]*Det3_134_013;
1563 G4double Det4_0234_0123 = m[A00] * Det3_234 << 1473 G4double Det4_0234_0123 = m[A00]*Det3_234_123 - m[A01]*Det3_234_023
1564 m[A02] * Det3_234 << 1474 + m[A02]*Det3_234_013 - m[A03]*Det3_234_012;
1565 G4double Det4_0234_0124 = m[A00] * Det3_234 << 1475 G4double Det4_0234_0124 = m[A00]*Det3_234_124 - m[A01]*Det3_234_024
1566 m[A02] * Det3_234 << 1476 + m[A02]*Det3_234_014 - m[A04]*Det3_234_012;
1567 G4double Det4_0234_0134 = m[A00] * Det3_234 << 1477 G4double Det4_0234_0134 = m[A00]*Det3_234_134 - m[A01]*Det3_234_034
1568 m[A03] * Det3_234 << 1478 + m[A03]*Det3_234_014 - m[A04]*Det3_234_013;
1569 G4double Det4_0234_0234 = m[A00] * Det3_234 << 1479 G4double Det4_0234_0234 = m[A00]*Det3_234_234 - m[A02]*Det3_234_034
1570 m[A03] * Det3_234 << 1480 + m[A03]*Det3_234_024 - m[A04]*Det3_234_023;
1571 G4double Det4_1234_0123 = m[A10] * Det3_234 << 1481 G4double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023
1572 m[A12] * Det3_234 << 1482 + m[A12]*Det3_234_013 - m[A13]*Det3_234_012;
1573 G4double Det4_1234_0124 = m[A10] * Det3_234 << 1483 G4double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024
1574 m[A12] * Det3_234 << 1484 + m[A12]*Det3_234_014 - m[A14]*Det3_234_012;
1575 G4double Det4_1234_0134 = m[A10] * Det3_234 << 1485 G4double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034
1576 m[A13] * Det3_234 << 1486 + m[A13]*Det3_234_014 - m[A14]*Det3_234_013;
1577 G4double Det4_1234_0234 = m[A10] * Det3_234 << 1487 G4double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034
1578 m[A13] * Det3_234 << 1488 + m[A13]*Det3_234_024 - m[A14]*Det3_234_023;
1579 G4double Det4_1234_1234 = m[A11] * Det3_234 << 1489 G4double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134
1580 m[A13] * Det3_234 << 1490 + m[A13]*Det3_234_124 - m[A14]*Det3_234_123;
1581 1491
1582 // Find the 5x5 det: 1492 // Find the 5x5 det:
1583 1493
1584 G4double det = m[A00] * Det4_1234_1234 - m[ << 1494 G4double det = m[A00]*Det4_1234_1234
1585 m[A02] * Det4_1234_0134 - m[ << 1495 - m[A01]*Det4_1234_0234
1586 m[A04] * Det4_1234_0123; << 1496 + m[A02]*Det4_1234_0134
>> 1497 - m[A03]*Det4_1234_0124
>> 1498 + m[A04]*Det4_1234_0123;
1587 1499
1588 if(det == 0) << 1500 if ( det == 0 )
1589 { << 1501 {
1590 ifail = 1; 1502 ifail = 1;
1591 return; 1503 return;
1592 } << 1504 }
1593 1505
1594 G4double oneOverDet = 1.0 / det; << 1506 G4double oneOverDet = 1.0/det;
1595 G4double mn1OverDet = -oneOverDet; << 1507 G4double mn1OverDet = - oneOverDet;
1596 1508
1597 m[A00] = Det4_1234_1234 * oneOverDet; << 1509 m[A00] = Det4_1234_1234 * oneOverDet;
1598 m[A01] = Det4_1234_0234 * mn1OverDet; << 1510 m[A01] = Det4_1234_0234 * mn1OverDet;
1599 m[A02] = Det4_1234_0134 * oneOverDet; << 1511 m[A02] = Det4_1234_0134 * oneOverDet;
1600 m[A03] = Det4_1234_0124 * mn1OverDet; << 1512 m[A03] = Det4_1234_0124 * mn1OverDet;
1601 m[A04] = Det4_1234_0123 * oneOverDet; << 1513 m[A04] = Det4_1234_0123 * oneOverDet;
1602 << 1514
1603 m[A11] = Det4_0234_0234 * oneOverDet; << 1515 m[A11] = Det4_0234_0234 * oneOverDet;
1604 m[A12] = Det4_0234_0134 * mn1OverDet; << 1516 m[A12] = Det4_0234_0134 * mn1OverDet;
1605 m[A13] = Det4_0234_0124 * oneOverDet; << 1517 m[A13] = Det4_0234_0124 * oneOverDet;
1606 m[A14] = Det4_0234_0123 * mn1OverDet; << 1518 m[A14] = Det4_0234_0123 * mn1OverDet;
1607 << 1519
1608 m[A22] = Det4_0134_0134 * oneOverDet; << 1520 m[A22] = Det4_0134_0134 * oneOverDet;
1609 m[A23] = Det4_0134_0124 * mn1OverDet; << 1521 m[A23] = Det4_0134_0124 * mn1OverDet;
1610 m[A24] = Det4_0134_0123 * oneOverDet; << 1522 m[A24] = Det4_0134_0123 * oneOverDet;
1611 1523
1612 m[A33] = Det4_0124_0124 * oneOverDet; << 1524 m[A33] = Det4_0124_0124 * oneOverDet;
1613 m[A34] = Det4_0124_0123 * mn1OverDet; << 1525 m[A34] = Det4_0124_0123 * mn1OverDet;
1614 1526
1615 m[A44] = Det4_0123_0123 * oneOverDet; << 1527 m[A44] = Det4_0123_0123 * oneOverDet;
1616 1528
1617 return; 1529 return;
1618 } 1530 }
1619 1531
1620 void G4ErrorSymMatrix::invertHaywood6(G4int& << 1532 void G4ErrorSymMatrix::invertHaywood6 (G4int & ifail)
1621 { 1533 {
1622 ifail = 0; 1534 ifail = 0;
1623 1535
1624 // Find all NECESSARY 2x2 dets: (39 of the 1536 // Find all NECESSARY 2x2 dets: (39 of them)
1625 1537
1626 G4double Det2_34_01 = m[A30] * m[A41] - m[A << 1538 G4double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40];
1627 G4double Det2_34_02 = m[A30] * m[A42] - m[A << 1539 G4double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40];
1628 G4double Det2_34_03 = m[A30] * m[A43] - m[A << 1540 G4double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40];
1629 G4double Det2_34_04 = m[A30] * m[A44] - m[A << 1541 G4double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40];
1630 G4double Det2_34_12 = m[A31] * m[A42] - m[A << 1542 G4double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41];
1631 G4double Det2_34_13 = m[A31] * m[A43] - m[A << 1543 G4double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41];
1632 G4double Det2_34_14 = m[A31] * m[A44] - m[A << 1544 G4double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41];
1633 G4double Det2_34_23 = m[A32] * m[A43] - m[A << 1545 G4double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42];
1634 G4double Det2_34_24 = m[A32] * m[A44] - m[A << 1546 G4double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42];
1635 G4double Det2_34_34 = m[A33] * m[A44] - m[A << 1547 G4double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43];
1636 G4double Det2_35_01 = m[A30] * m[A51] - m[A << 1548 G4double Det2_35_01 = m[A30]*m[A51] - m[A31]*m[A50];
1637 G4double Det2_35_02 = m[A30] * m[A52] - m[A << 1549 G4double Det2_35_02 = m[A30]*m[A52] - m[A32]*m[A50];
1638 G4double Det2_35_03 = m[A30] * m[A53] - m[A << 1550 G4double Det2_35_03 = m[A30]*m[A53] - m[A33]*m[A50];
1639 G4double Det2_35_04 = m[A30] * m[A54] - m[A << 1551 G4double Det2_35_04 = m[A30]*m[A54] - m[A34]*m[A50];
1640 G4double Det2_35_05 = m[A30] * m[A55] - m[A << 1552 G4double Det2_35_05 = m[A30]*m[A55] - m[A35]*m[A50];
1641 G4double Det2_35_12 = m[A31] * m[A52] - m[A << 1553 G4double Det2_35_12 = m[A31]*m[A52] - m[A32]*m[A51];
1642 G4double Det2_35_13 = m[A31] * m[A53] - m[A << 1554 G4double Det2_35_13 = m[A31]*m[A53] - m[A33]*m[A51];
1643 G4double Det2_35_14 = m[A31] * m[A54] - m[A << 1555 G4double Det2_35_14 = m[A31]*m[A54] - m[A34]*m[A51];
1644 G4double Det2_35_15 = m[A31] * m[A55] - m[A << 1556 G4double Det2_35_15 = m[A31]*m[A55] - m[A35]*m[A51];
1645 G4double Det2_35_23 = m[A32] * m[A53] - m[A << 1557 G4double Det2_35_23 = m[A32]*m[A53] - m[A33]*m[A52];
1646 G4double Det2_35_24 = m[A32] * m[A54] - m[A << 1558 G4double Det2_35_24 = m[A32]*m[A54] - m[A34]*m[A52];
1647 G4double Det2_35_25 = m[A32] * m[A55] - m[A << 1559 G4double Det2_35_25 = m[A32]*m[A55] - m[A35]*m[A52];
1648 G4double Det2_35_34 = m[A33] * m[A54] - m[A << 1560 G4double Det2_35_34 = m[A33]*m[A54] - m[A34]*m[A53];
1649 G4double Det2_35_35 = m[A33] * m[A55] - m[A << 1561 G4double Det2_35_35 = m[A33]*m[A55] - m[A35]*m[A53];
1650 G4double Det2_45_01 = m[A40] * m[A51] - m[A << 1562 G4double Det2_45_01 = m[A40]*m[A51] - m[A41]*m[A50];
1651 G4double Det2_45_02 = m[A40] * m[A52] - m[A << 1563 G4double Det2_45_02 = m[A40]*m[A52] - m[A42]*m[A50];
1652 G4double Det2_45_03 = m[A40] * m[A53] - m[A << 1564 G4double Det2_45_03 = m[A40]*m[A53] - m[A43]*m[A50];
1653 G4double Det2_45_04 = m[A40] * m[A54] - m[A << 1565 G4double Det2_45_04 = m[A40]*m[A54] - m[A44]*m[A50];
1654 G4double Det2_45_05 = m[A40] * m[A55] - m[A << 1566 G4double Det2_45_05 = m[A40]*m[A55] - m[A45]*m[A50];
1655 G4double Det2_45_12 = m[A41] * m[A52] - m[A << 1567 G4double Det2_45_12 = m[A41]*m[A52] - m[A42]*m[A51];
1656 G4double Det2_45_13 = m[A41] * m[A53] - m[A << 1568 G4double Det2_45_13 = m[A41]*m[A53] - m[A43]*m[A51];
1657 G4double Det2_45_14 = m[A41] * m[A54] - m[A << 1569 G4double Det2_45_14 = m[A41]*m[A54] - m[A44]*m[A51];
1658 G4double Det2_45_15 = m[A41] * m[A55] - m[A << 1570 G4double Det2_45_15 = m[A41]*m[A55] - m[A45]*m[A51];
1659 G4double Det2_45_23 = m[A42] * m[A53] - m[A << 1571 G4double Det2_45_23 = m[A42]*m[A53] - m[A43]*m[A52];
1660 G4double Det2_45_24 = m[A42] * m[A54] - m[A << 1572 G4double Det2_45_24 = m[A42]*m[A54] - m[A44]*m[A52];
1661 G4double Det2_45_25 = m[A42] * m[A55] - m[A << 1573 G4double Det2_45_25 = m[A42]*m[A55] - m[A45]*m[A52];
1662 G4double Det2_45_34 = m[A43] * m[A54] - m[A << 1574 G4double Det2_45_34 = m[A43]*m[A54] - m[A44]*m[A53];
1663 G4double Det2_45_35 = m[A43] * m[A55] - m[A << 1575 G4double Det2_45_35 = m[A43]*m[A55] - m[A45]*m[A53];
1664 G4double Det2_45_45 = m[A44] * m[A55] - m[A << 1576 G4double Det2_45_45 = m[A44]*m[A55] - m[A45]*m[A54];
1665 1577
1666 // Find all NECESSARY 3x3 dets: (65 of the 1578 // Find all NECESSARY 3x3 dets: (65 of them)
1667 1579
1668 G4double Det3_234_012 = << 1580 G4double Det3_234_012 = m[A20]*Det2_34_12 - m[A21]*Det2_34_02
1669 m[A20] * Det2_34_12 - m[A21] * Det2_34_02 << 1581 + m[A22]*Det2_34_01;
1670 G4double Det3_234_013 = << 1582 G4double Det3_234_013 = m[A20]*Det2_34_13 - m[A21]*Det2_34_03
1671 m[A20] * Det2_34_13 - m[A21] * Det2_34_03 << 1583 + m[A23]*Det2_34_01;
1672 G4double Det3_234_014 = << 1584 G4double Det3_234_014 = m[A20]*Det2_34_14 - m[A21]*Det2_34_04
1673 m[A20] * Det2_34_14 - m[A21] * Det2_34_04 << 1585 + m[A24]*Det2_34_01;
1674 G4double Det3_234_023 = << 1586 G4double Det3_234_023 = m[A20]*Det2_34_23 - m[A22]*Det2_34_03
1675 m[A20] * Det2_34_23 - m[A22] * Det2_34_03 << 1587 + m[A23]*Det2_34_02;
1676 G4double Det3_234_024 = << 1588 G4double Det3_234_024 = m[A20]*Det2_34_24 - m[A22]*Det2_34_04
1677 m[A20] * Det2_34_24 - m[A22] * Det2_34_04 << 1589 + m[A24]*Det2_34_02;
1678 G4double Det3_234_034 = << 1590 G4double Det3_234_034 = m[A20]*Det2_34_34 - m[A23]*Det2_34_04
1679 m[A20] * Det2_34_34 - m[A23] * Det2_34_04 << 1591 + m[A24]*Det2_34_03;
1680 G4double Det3_234_123 = << 1592 G4double Det3_234_123 = m[A21]*Det2_34_23 - m[A22]*Det2_34_13
1681 m[A21] * Det2_34_23 - m[A22] * Det2_34_13 << 1593 + m[A23]*Det2_34_12;
1682 G4double Det3_234_124 = << 1594 G4double Det3_234_124 = m[A21]*Det2_34_24 - m[A22]*Det2_34_14
1683 m[A21] * Det2_34_24 - m[A22] * Det2_34_14 << 1595 + m[A24]*Det2_34_12;
1684 G4double Det3_234_134 = << 1596 G4double Det3_234_134 = m[A21]*Det2_34_34 - m[A23]*Det2_34_14
1685 m[A21] * Det2_34_34 - m[A23] * Det2_34_14 << 1597 + m[A24]*Det2_34_13;
1686 G4double Det3_234_234 = << 1598 G4double Det3_234_234 = m[A22]*Det2_34_34 - m[A23]*Det2_34_24
1687 m[A22] * Det2_34_34 - m[A23] * Det2_34_24 << 1599 + m[A24]*Det2_34_23;
1688 G4double Det3_235_012 = << 1600 G4double Det3_235_012 = m[A20]*Det2_35_12 - m[A21]*Det2_35_02
1689 m[A20] * Det2_35_12 - m[A21] * Det2_35_02 << 1601 + m[A22]*Det2_35_01;
1690 G4double Det3_235_013 = << 1602 G4double Det3_235_013 = m[A20]*Det2_35_13 - m[A21]*Det2_35_03
1691 m[A20] * Det2_35_13 - m[A21] * Det2_35_03 << 1603 + m[A23]*Det2_35_01;
1692 G4double Det3_235_014 = << 1604 G4double Det3_235_014 = m[A20]*Det2_35_14 - m[A21]*Det2_35_04
1693 m[A20] * Det2_35_14 - m[A21] * Det2_35_04 << 1605 + m[A24]*Det2_35_01;
1694 G4double Det3_235_015 = << 1606 G4double Det3_235_015 = m[A20]*Det2_35_15 - m[A21]*Det2_35_05
1695 m[A20] * Det2_35_15 - m[A21] * Det2_35_05 << 1607 + m[A25]*Det2_35_01;
1696 G4double Det3_235_023 = << 1608 G4double Det3_235_023 = m[A20]*Det2_35_23 - m[A22]*Det2_35_03
1697 m[A20] * Det2_35_23 - m[A22] * Det2_35_03 << 1609 + m[A23]*Det2_35_02;
1698 G4double Det3_235_024 = << 1610 G4double Det3_235_024 = m[A20]*Det2_35_24 - m[A22]*Det2_35_04
1699 m[A20] * Det2_35_24 - m[A22] * Det2_35_04 << 1611 + m[A24]*Det2_35_02;
1700 G4double Det3_235_025 = << 1612 G4double Det3_235_025 = m[A20]*Det2_35_25 - m[A22]*Det2_35_05
1701 m[A20] * Det2_35_25 - m[A22] * Det2_35_05 << 1613 + m[A25]*Det2_35_02;
1702 G4double Det3_235_034 = << 1614 G4double Det3_235_034 = m[A20]*Det2_35_34 - m[A23]*Det2_35_04
1703 m[A20] * Det2_35_34 - m[A23] * Det2_35_04 << 1615 + m[A24]*Det2_35_03;
1704 G4double Det3_235_035 = << 1616 G4double Det3_235_035 = m[A20]*Det2_35_35 - m[A23]*Det2_35_05
1705 m[A20] * Det2_35_35 - m[A23] * Det2_35_05 << 1617 + m[A25]*Det2_35_03;
1706 G4double Det3_235_123 = << 1618 G4double Det3_235_123 = m[A21]*Det2_35_23 - m[A22]*Det2_35_13
1707 m[A21] * Det2_35_23 - m[A22] * Det2_35_13 << 1619 + m[A23]*Det2_35_12;
1708 G4double Det3_235_124 = << 1620 G4double Det3_235_124 = m[A21]*Det2_35_24 - m[A22]*Det2_35_14
1709 m[A21] * Det2_35_24 - m[A22] * Det2_35_14 << 1621 + m[A24]*Det2_35_12;
1710 G4double Det3_235_125 = << 1622 G4double Det3_235_125 = m[A21]*Det2_35_25 - m[A22]*Det2_35_15
1711 m[A21] * Det2_35_25 - m[A22] * Det2_35_15 << 1623 + m[A25]*Det2_35_12;
1712 G4double Det3_235_134 = << 1624 G4double Det3_235_134 = m[A21]*Det2_35_34 - m[A23]*Det2_35_14
1713 m[A21] * Det2_35_34 - m[A23] * Det2_35_14 << 1625 + m[A24]*Det2_35_13;
1714 G4double Det3_235_135 = << 1626 G4double Det3_235_135 = m[A21]*Det2_35_35 - m[A23]*Det2_35_15
1715 m[A21] * Det2_35_35 - m[A23] * Det2_35_15 << 1627 + m[A25]*Det2_35_13;
1716 G4double Det3_235_234 = << 1628 G4double Det3_235_234 = m[A22]*Det2_35_34 - m[A23]*Det2_35_24
1717 m[A22] * Det2_35_34 - m[A23] * Det2_35_24 << 1629 + m[A24]*Det2_35_23;
1718 G4double Det3_235_235 = << 1630 G4double Det3_235_235 = m[A22]*Det2_35_35 - m[A23]*Det2_35_25
1719 m[A22] * Det2_35_35 - m[A23] * Det2_35_25 << 1631 + m[A25]*Det2_35_23;
1720 G4double Det3_245_012 = << 1632 G4double Det3_245_012 = m[A20]*Det2_45_12 - m[A21]*Det2_45_02
1721 m[A20] * Det2_45_12 - m[A21] * Det2_45_02 << 1633 + m[A22]*Det2_45_01;
1722 G4double Det3_245_013 = << 1634 G4double Det3_245_013 = m[A20]*Det2_45_13 - m[A21]*Det2_45_03
1723 m[A20] * Det2_45_13 - m[A21] * Det2_45_03 << 1635 + m[A23]*Det2_45_01;
1724 G4double Det3_245_014 = << 1636 G4double Det3_245_014 = m[A20]*Det2_45_14 - m[A21]*Det2_45_04
1725 m[A20] * Det2_45_14 - m[A21] * Det2_45_04 << 1637 + m[A24]*Det2_45_01;
1726 G4double Det3_245_015 = << 1638 G4double Det3_245_015 = m[A20]*Det2_45_15 - m[A21]*Det2_45_05
1727 m[A20] * Det2_45_15 - m[A21] * Det2_45_05 << 1639 + m[A25]*Det2_45_01;
1728 G4double Det3_245_023 = << 1640 G4double Det3_245_023 = m[A20]*Det2_45_23 - m[A22]*Det2_45_03
1729 m[A20] * Det2_45_23 - m[A22] * Det2_45_03 << 1641 + m[A23]*Det2_45_02;
1730 G4double Det3_245_024 = << 1642 G4double Det3_245_024 = m[A20]*Det2_45_24 - m[A22]*Det2_45_04
1731 m[A20] * Det2_45_24 - m[A22] * Det2_45_04 << 1643 + m[A24]*Det2_45_02;
1732 G4double Det3_245_025 = << 1644 G4double Det3_245_025 = m[A20]*Det2_45_25 - m[A22]*Det2_45_05
1733 m[A20] * Det2_45_25 - m[A22] * Det2_45_05 << 1645 + m[A25]*Det2_45_02;
1734 G4double Det3_245_034 = << 1646 G4double Det3_245_034 = m[A20]*Det2_45_34 - m[A23]*Det2_45_04
1735 m[A20] * Det2_45_34 - m[A23] * Det2_45_04 << 1647 + m[A24]*Det2_45_03;
1736 G4double Det3_245_035 = << 1648 G4double Det3_245_035 = m[A20]*Det2_45_35 - m[A23]*Det2_45_05
1737 m[A20] * Det2_45_35 - m[A23] * Det2_45_05 << 1649 + m[A25]*Det2_45_03;
1738 G4double Det3_245_045 = << 1650 G4double Det3_245_045 = m[A20]*Det2_45_45 - m[A24]*Det2_45_05
1739 m[A20] * Det2_45_45 - m[A24] * Det2_45_05 << 1651 + m[A25]*Det2_45_04;
1740 G4double Det3_245_123 = << 1652 G4double Det3_245_123 = m[A21]*Det2_45_23 - m[A22]*Det2_45_13
1741 m[A21] * Det2_45_23 - m[A22] * Det2_45_13 << 1653 + m[A23]*Det2_45_12;
1742 G4double Det3_245_124 = << 1654 G4double Det3_245_124 = m[A21]*Det2_45_24 - m[A22]*Det2_45_14
1743 m[A21] * Det2_45_24 - m[A22] * Det2_45_14 << 1655 + m[A24]*Det2_45_12;
1744 G4double Det3_245_125 = << 1656 G4double Det3_245_125 = m[A21]*Det2_45_25 - m[A22]*Det2_45_15
1745 m[A21] * Det2_45_25 - m[A22] * Det2_45_15 << 1657 + m[A25]*Det2_45_12;
1746 G4double Det3_245_134 = << 1658 G4double Det3_245_134 = m[A21]*Det2_45_34 - m[A23]*Det2_45_14
1747 m[A21] * Det2_45_34 - m[A23] * Det2_45_14 << 1659 + m[A24]*Det2_45_13;
1748 G4double Det3_245_135 = << 1660 G4double Det3_245_135 = m[A21]*Det2_45_35 - m[A23]*Det2_45_15
1749 m[A21] * Det2_45_35 - m[A23] * Det2_45_15 << 1661 + m[A25]*Det2_45_13;
1750 G4double Det3_245_145 = << 1662 G4double Det3_245_145 = m[A21]*Det2_45_45 - m[A24]*Det2_45_15
1751 m[A21] * Det2_45_45 - m[A24] * Det2_45_15 << 1663 + m[A25]*Det2_45_14;
1752 G4double Det3_245_234 = << 1664 G4double Det3_245_234 = m[A22]*Det2_45_34 - m[A23]*Det2_45_24
1753 m[A22] * Det2_45_34 - m[A23] * Det2_45_24 << 1665 + m[A24]*Det2_45_23;
1754 G4double Det3_245_235 = << 1666 G4double Det3_245_235 = m[A22]*Det2_45_35 - m[A23]*Det2_45_25
1755 m[A22] * Det2_45_35 - m[A23] * Det2_45_25 << 1667 + m[A25]*Det2_45_23;
1756 G4double Det3_245_245 = << 1668 G4double Det3_245_245 = m[A22]*Det2_45_45 - m[A24]*Det2_45_25
1757 m[A22] * Det2_45_45 - m[A24] * Det2_45_25 << 1669 + m[A25]*Det2_45_24;
1758 G4double Det3_345_012 = << 1670 G4double Det3_345_012 = m[A30]*Det2_45_12 - m[A31]*Det2_45_02
1759 m[A30] * Det2_45_12 - m[A31] * Det2_45_02 << 1671 + m[A32]*Det2_45_01;
1760 G4double Det3_345_013 = << 1672 G4double Det3_345_013 = m[A30]*Det2_45_13 - m[A31]*Det2_45_03
1761 m[A30] * Det2_45_13 - m[A31] * Det2_45_03 << 1673 + m[A33]*Det2_45_01;
1762 G4double Det3_345_014 = << 1674 G4double Det3_345_014 = m[A30]*Det2_45_14 - m[A31]*Det2_45_04
1763 m[A30] * Det2_45_14 - m[A31] * Det2_45_04 << 1675 + m[A34]*Det2_45_01;
1764 G4double Det3_345_015 = << 1676 G4double Det3_345_015 = m[A30]*Det2_45_15 - m[A31]*Det2_45_05
1765 m[A30] * Det2_45_15 - m[A31] * Det2_45_05 << 1677 + m[A35]*Det2_45_01;
1766 G4double Det3_345_023 = << 1678 G4double Det3_345_023 = m[A30]*Det2_45_23 - m[A32]*Det2_45_03
1767 m[A30] * Det2_45_23 - m[A32] * Det2_45_03 << 1679 + m[A33]*Det2_45_02;
1768 G4double Det3_345_024 = << 1680 G4double Det3_345_024 = m[A30]*Det2_45_24 - m[A32]*Det2_45_04
1769 m[A30] * Det2_45_24 - m[A32] * Det2_45_04 << 1681 + m[A34]*Det2_45_02;
1770 G4double Det3_345_025 = << 1682 G4double Det3_345_025 = m[A30]*Det2_45_25 - m[A32]*Det2_45_05
1771 m[A30] * Det2_45_25 - m[A32] * Det2_45_05 << 1683 + m[A35]*Det2_45_02;
1772 G4double Det3_345_034 = << 1684 G4double Det3_345_034 = m[A30]*Det2_45_34 - m[A33]*Det2_45_04
1773 m[A30] * Det2_45_34 - m[A33] * Det2_45_04 << 1685 + m[A34]*Det2_45_03;
1774 G4double Det3_345_035 = << 1686 G4double Det3_345_035 = m[A30]*Det2_45_35 - m[A33]*Det2_45_05
1775 m[A30] * Det2_45_35 - m[A33] * Det2_45_05 << 1687 + m[A35]*Det2_45_03;
1776 G4double Det3_345_045 = << 1688 G4double Det3_345_045 = m[A30]*Det2_45_45 - m[A34]*Det2_45_05
1777 m[A30] * Det2_45_45 - m[A34] * Det2_45_05 << 1689 + m[A35]*Det2_45_04;
1778 G4double Det3_345_123 = << 1690 G4double Det3_345_123 = m[A31]*Det2_45_23 - m[A32]*Det2_45_13
1779 m[A31] * Det2_45_23 - m[A32] * Det2_45_13 << 1691 + m[A33]*Det2_45_12;
1780 G4double Det3_345_124 = << 1692 G4double Det3_345_124 = m[A31]*Det2_45_24 - m[A32]*Det2_45_14
1781 m[A31] * Det2_45_24 - m[A32] * Det2_45_14 << 1693 + m[A34]*Det2_45_12;
1782 G4double Det3_345_125 = << 1694 G4double Det3_345_125 = m[A31]*Det2_45_25 - m[A32]*Det2_45_15
1783 m[A31] * Det2_45_25 - m[A32] * Det2_45_15 << 1695 + m[A35]*Det2_45_12;
1784 G4double Det3_345_134 = << 1696 G4double Det3_345_134 = m[A31]*Det2_45_34 - m[A33]*Det2_45_14
1785 m[A31] * Det2_45_34 - m[A33] * Det2_45_14 << 1697 + m[A34]*Det2_45_13;
1786 G4double Det3_345_135 = << 1698 G4double Det3_345_135 = m[A31]*Det2_45_35 - m[A33]*Det2_45_15
1787 m[A31] * Det2_45_35 - m[A33] * Det2_45_15 << 1699 + m[A35]*Det2_45_13;
1788 G4double Det3_345_145 = << 1700 G4double Det3_345_145 = m[A31]*Det2_45_45 - m[A34]*Det2_45_15
1789 m[A31] * Det2_45_45 - m[A34] * Det2_45_15 << 1701 + m[A35]*Det2_45_14;
1790 G4double Det3_345_234 = << 1702 G4double Det3_345_234 = m[A32]*Det2_45_34 - m[A33]*Det2_45_24
1791 m[A32] * Det2_45_34 - m[A33] * Det2_45_24 << 1703 + m[A34]*Det2_45_23;
1792 G4double Det3_345_235 = << 1704 G4double Det3_345_235 = m[A32]*Det2_45_35 - m[A33]*Det2_45_25
1793 m[A32] * Det2_45_35 - m[A33] * Det2_45_25 << 1705 + m[A35]*Det2_45_23;
1794 G4double Det3_345_245 = << 1706 G4double Det3_345_245 = m[A32]*Det2_45_45 - m[A34]*Det2_45_25
1795 m[A32] * Det2_45_45 - m[A34] * Det2_45_25 << 1707 + m[A35]*Det2_45_24;
1796 G4double Det3_345_345 = << 1708 G4double Det3_345_345 = m[A33]*Det2_45_45 - m[A34]*Det2_45_35
1797 m[A33] * Det2_45_45 - m[A34] * Det2_45_35 << 1709 + m[A35]*Det2_45_34;
1798 1710
1799 // Find all NECESSARY 4x4 dets: (55 of the 1711 // Find all NECESSARY 4x4 dets: (55 of them)
1800 1712
1801 G4double Det4_1234_0123 = m[A10] * Det3_234 << 1713 G4double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023
1802 m[A12] * Det3_234 << 1714 + m[A12]*Det3_234_013 - m[A13]*Det3_234_012;
1803 G4double Det4_1234_0124 = m[A10] * Det3_234 << 1715 G4double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024
1804 m[A12] * Det3_234 << 1716 + m[A12]*Det3_234_014 - m[A14]*Det3_234_012;
1805 G4double Det4_1234_0134 = m[A10] * Det3_234 << 1717 G4double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034
1806 m[A13] * Det3_234 << 1718 + m[A13]*Det3_234_014 - m[A14]*Det3_234_013;
1807 G4double Det4_1234_0234 = m[A10] * Det3_234 << 1719 G4double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034
1808 m[A13] * Det3_234 << 1720 + m[A13]*Det3_234_024 - m[A14]*Det3_234_023;
1809 G4double Det4_1234_1234 = m[A11] * Det3_234 << 1721 G4double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134
1810 m[A13] * Det3_234 << 1722 + m[A13]*Det3_234_124 - m[A14]*Det3_234_123;
1811 G4double Det4_1235_0123 = m[A10] * Det3_235 << 1723 G4double Det4_1235_0123 = m[A10]*Det3_235_123 - m[A11]*Det3_235_023
1812 m[A12] * Det3_235 << 1724 + m[A12]*Det3_235_013 - m[A13]*Det3_235_012;
1813 G4double Det4_1235_0124 = m[A10] * Det3_235 << 1725 G4double Det4_1235_0124 = m[A10]*Det3_235_124 - m[A11]*Det3_235_024
1814 m[A12] * Det3_235 << 1726 + m[A12]*Det3_235_014 - m[A14]*Det3_235_012;
1815 G4double Det4_1235_0125 = m[A10] * Det3_235 << 1727 G4double Det4_1235_0125 = m[A10]*Det3_235_125 - m[A11]*Det3_235_025
1816 m[A12] * Det3_235 << 1728 + m[A12]*Det3_235_015 - m[A15]*Det3_235_012;
1817 G4double Det4_1235_0134 = m[A10] * Det3_235 << 1729 G4double Det4_1235_0134 = m[A10]*Det3_235_134 - m[A11]*Det3_235_034
1818 m[A13] * Det3_235 << 1730 + m[A13]*Det3_235_014 - m[A14]*Det3_235_013;
1819 G4double Det4_1235_0135 = m[A10] * Det3_235 << 1731 G4double Det4_1235_0135 = m[A10]*Det3_235_135 - m[A11]*Det3_235_035
1820 m[A13] * Det3_235 << 1732 + m[A13]*Det3_235_015 - m[A15]*Det3_235_013;
1821 G4double Det4_1235_0234 = m[A10] * Det3_235 << 1733 G4double Det4_1235_0234 = m[A10]*Det3_235_234 - m[A12]*Det3_235_034
1822 m[A13] * Det3_235 << 1734 + m[A13]*Det3_235_024 - m[A14]*Det3_235_023;
1823 G4double Det4_1235_0235 = m[A10] * Det3_235 << 1735 G4double Det4_1235_0235 = m[A10]*Det3_235_235 - m[A12]*Det3_235_035
1824 m[A13] * Det3_235 << 1736 + m[A13]*Det3_235_025 - m[A15]*Det3_235_023;
1825 G4double Det4_1235_1234 = m[A11] * Det3_235 << 1737 G4double Det4_1235_1234 = m[A11]*Det3_235_234 - m[A12]*Det3_235_134
1826 m[A13] * Det3_235 << 1738 + m[A13]*Det3_235_124 - m[A14]*Det3_235_123;
1827 G4double Det4_1235_1235 = m[A11] * Det3_235 << 1739 G4double Det4_1235_1235 = m[A11]*Det3_235_235 - m[A12]*Det3_235_135
1828 m[A13] * Det3_235 << 1740 + m[A13]*Det3_235_125 - m[A15]*Det3_235_123;
1829 G4double Det4_1245_0123 = m[A10] * Det3_245 << 1741 G4double Det4_1245_0123 = m[A10]*Det3_245_123 - m[A11]*Det3_245_023
1830 m[A12] * Det3_245 << 1742 + m[A12]*Det3_245_013 - m[A13]*Det3_245_012;
1831 G4double Det4_1245_0124 = m[A10] * Det3_245 << 1743 G4double Det4_1245_0124 = m[A10]*Det3_245_124 - m[A11]*Det3_245_024
1832 m[A12] * Det3_245 << 1744 + m[A12]*Det3_245_014 - m[A14]*Det3_245_012;
1833 G4double Det4_1245_0125 = m[A10] * Det3_245 << 1745 G4double Det4_1245_0125 = m[A10]*Det3_245_125 - m[A11]*Det3_245_025
1834 m[A12] * Det3_245 << 1746 + m[A12]*Det3_245_015 - m[A15]*Det3_245_012;
1835 G4double Det4_1245_0134 = m[A10] * Det3_245 << 1747 G4double Det4_1245_0134 = m[A10]*Det3_245_134 - m[A11]*Det3_245_034
1836 m[A13] * Det3_245 << 1748 + m[A13]*Det3_245_014 - m[A14]*Det3_245_013;
1837 G4double Det4_1245_0135 = m[A10] * Det3_245 << 1749 G4double Det4_1245_0135 = m[A10]*Det3_245_135 - m[A11]*Det3_245_035
1838 m[A13] * Det3_245 << 1750 + m[A13]*Det3_245_015 - m[A15]*Det3_245_013;
1839 G4double Det4_1245_0145 = m[A10] * Det3_245 << 1751 G4double Det4_1245_0145 = m[A10]*Det3_245_145 - m[A11]*Det3_245_045
1840 m[A14] * Det3_245 << 1752 + m[A14]*Det3_245_015 - m[A15]*Det3_245_014;
1841 G4double Det4_1245_0234 = m[A10] * Det3_245 << 1753 G4double Det4_1245_0234 = m[A10]*Det3_245_234 - m[A12]*Det3_245_034
1842 m[A13] * Det3_245 << 1754 + m[A13]*Det3_245_024 - m[A14]*Det3_245_023;
1843 G4double Det4_1245_0235 = m[A10] * Det3_245 << 1755 G4double Det4_1245_0235 = m[A10]*Det3_245_235 - m[A12]*Det3_245_035
1844 m[A13] * Det3_245 << 1756 + m[A13]*Det3_245_025 - m[A15]*Det3_245_023;
1845 G4double Det4_1245_0245 = m[A10] * Det3_245 << 1757 G4double Det4_1245_0245 = m[A10]*Det3_245_245 - m[A12]*Det3_245_045
1846 m[A14] * Det3_245 << 1758 + m[A14]*Det3_245_025 - m[A15]*Det3_245_024;
1847 G4double Det4_1245_1234 = m[A11] * Det3_245 << 1759 G4double Det4_1245_1234 = m[A11]*Det3_245_234 - m[A12]*Det3_245_134
1848 m[A13] * Det3_245 << 1760 + m[A13]*Det3_245_124 - m[A14]*Det3_245_123;
1849 G4double Det4_1245_1235 = m[A11] * Det3_245 << 1761 G4double Det4_1245_1235 = m[A11]*Det3_245_235 - m[A12]*Det3_245_135
1850 m[A13] * Det3_245 << 1762 + m[A13]*Det3_245_125 - m[A15]*Det3_245_123;
1851 G4double Det4_1245_1245 = m[A11] * Det3_245 << 1763 G4double Det4_1245_1245 = m[A11]*Det3_245_245 - m[A12]*Det3_245_145
1852 m[A14] * Det3_245 << 1764 + m[A14]*Det3_245_125 - m[A15]*Det3_245_124;
1853 G4double Det4_1345_0123 = m[A10] * Det3_345 << 1765 G4double Det4_1345_0123 = m[A10]*Det3_345_123 - m[A11]*Det3_345_023
1854 m[A12] * Det3_345 << 1766 + m[A12]*Det3_345_013 - m[A13]*Det3_345_012;
1855 G4double Det4_1345_0124 = m[A10] * Det3_345 << 1767 G4double Det4_1345_0124 = m[A10]*Det3_345_124 - m[A11]*Det3_345_024
1856 m[A12] * Det3_345 << 1768 + m[A12]*Det3_345_014 - m[A14]*Det3_345_012;
1857 G4double Det4_1345_0125 = m[A10] * Det3_345 << 1769 G4double Det4_1345_0125 = m[A10]*Det3_345_125 - m[A11]*Det3_345_025
1858 m[A12] * Det3_345 << 1770 + m[A12]*Det3_345_015 - m[A15]*Det3_345_012;
1859 G4double Det4_1345_0134 = m[A10] * Det3_345 << 1771 G4double Det4_1345_0134 = m[A10]*Det3_345_134 - m[A11]*Det3_345_034
1860 m[A13] * Det3_345 << 1772 + m[A13]*Det3_345_014 - m[A14]*Det3_345_013;
1861 G4double Det4_1345_0135 = m[A10] * Det3_345 << 1773 G4double Det4_1345_0135 = m[A10]*Det3_345_135 - m[A11]*Det3_345_035
1862 m[A13] * Det3_345 << 1774 + m[A13]*Det3_345_015 - m[A15]*Det3_345_013;
1863 G4double Det4_1345_0145 = m[A10] * Det3_345 << 1775 G4double Det4_1345_0145 = m[A10]*Det3_345_145 - m[A11]*Det3_345_045
1864 m[A14] * Det3_345 << 1776 + m[A14]*Det3_345_015 - m[A15]*Det3_345_014;
1865 G4double Det4_1345_0234 = m[A10] * Det3_345 << 1777 G4double Det4_1345_0234 = m[A10]*Det3_345_234 - m[A12]*Det3_345_034
1866 m[A13] * Det3_345 << 1778 + m[A13]*Det3_345_024 - m[A14]*Det3_345_023;
1867 G4double Det4_1345_0235 = m[A10] * Det3_345 << 1779 G4double Det4_1345_0235 = m[A10]*Det3_345_235 - m[A12]*Det3_345_035
1868 m[A13] * Det3_345 << 1780 + m[A13]*Det3_345_025 - m[A15]*Det3_345_023;
1869 G4double Det4_1345_0245 = m[A10] * Det3_345 << 1781 G4double Det4_1345_0245 = m[A10]*Det3_345_245 - m[A12]*Det3_345_045
1870 m[A14] * Det3_345 << 1782 + m[A14]*Det3_345_025 - m[A15]*Det3_345_024;
1871 G4double Det4_1345_0345 = m[A10] * Det3_345 << 1783 G4double Det4_1345_0345 = m[A10]*Det3_345_345 - m[A13]*Det3_345_045
1872 m[A14] * Det3_345 << 1784 + m[A14]*Det3_345_035 - m[A15]*Det3_345_034;
1873 G4double Det4_1345_1234 = m[A11] * Det3_345 << 1785 G4double Det4_1345_1234 = m[A11]*Det3_345_234 - m[A12]*Det3_345_134
1874 m[A13] * Det3_345 << 1786 + m[A13]*Det3_345_124 - m[A14]*Det3_345_123;
1875 G4double Det4_1345_1235 = m[A11] * Det3_345 << 1787 G4double Det4_1345_1235 = m[A11]*Det3_345_235 - m[A12]*Det3_345_135
1876 m[A13] * Det3_345 << 1788 + m[A13]*Det3_345_125 - m[A15]*Det3_345_123;
1877 G4double Det4_1345_1245 = m[A11] * Det3_345 << 1789 G4double Det4_1345_1245 = m[A11]*Det3_345_245 - m[A12]*Det3_345_145
1878 m[A14] * Det3_345 << 1790 + m[A14]*Det3_345_125 - m[A15]*Det3_345_124;
1879 G4double Det4_1345_1345 = m[A11] * Det3_345 << 1791 G4double Det4_1345_1345 = m[A11]*Det3_345_345 - m[A13]*Det3_345_145
1880 m[A14] * Det3_345 << 1792 + m[A14]*Det3_345_135 - m[A15]*Det3_345_134;
1881 G4double Det4_2345_0123 = m[A20] * Det3_345 << 1793 G4double Det4_2345_0123 = m[A20]*Det3_345_123 - m[A21]*Det3_345_023
1882 m[A22] * Det3_345 << 1794 + m[A22]*Det3_345_013 - m[A23]*Det3_345_012;
1883 G4double Det4_2345_0124 = m[A20] * Det3_345 << 1795 G4double Det4_2345_0124 = m[A20]*Det3_345_124 - m[A21]*Det3_345_024
1884 m[A22] * Det3_345 << 1796 + m[A22]*Det3_345_014 - m[A24]*Det3_345_012;
1885 G4double Det4_2345_0125 = m[A20] * Det3_345 << 1797 G4double Det4_2345_0125 = m[A20]*Det3_345_125 - m[A21]*Det3_345_025
1886 m[A22] * Det3_345 << 1798 + m[A22]*Det3_345_015 - m[A25]*Det3_345_012;
1887 G4double Det4_2345_0134 = m[A20] * Det3_345 << 1799 G4double Det4_2345_0134 = m[A20]*Det3_345_134 - m[A21]*Det3_345_034
1888 m[A23] * Det3_345 << 1800 + m[A23]*Det3_345_014 - m[A24]*Det3_345_013;
1889 G4double Det4_2345_0135 = m[A20] * Det3_345 << 1801 G4double Det4_2345_0135 = m[A20]*Det3_345_135 - m[A21]*Det3_345_035
1890 m[A23] * Det3_345 << 1802 + m[A23]*Det3_345_015 - m[A25]*Det3_345_013;
1891 G4double Det4_2345_0145 = m[A20] * Det3_345 << 1803 G4double Det4_2345_0145 = m[A20]*Det3_345_145 - m[A21]*Det3_345_045
1892 m[A24] * Det3_345 << 1804 + m[A24]*Det3_345_015 - m[A25]*Det3_345_014;
1893 G4double Det4_2345_0234 = m[A20] * Det3_345 << 1805 G4double Det4_2345_0234 = m[A20]*Det3_345_234 - m[A22]*Det3_345_034
1894 m[A23] * Det3_345 << 1806 + m[A23]*Det3_345_024 - m[A24]*Det3_345_023;
1895 G4double Det4_2345_0235 = m[A20] * Det3_345 << 1807 G4double Det4_2345_0235 = m[A20]*Det3_345_235 - m[A22]*Det3_345_035
1896 m[A23] * Det3_345 << 1808 + m[A23]*Det3_345_025 - m[A25]*Det3_345_023;
1897 G4double Det4_2345_0245 = m[A20] * Det3_345 << 1809 G4double Det4_2345_0245 = m[A20]*Det3_345_245 - m[A22]*Det3_345_045
1898 m[A24] * Det3_345 << 1810 + m[A24]*Det3_345_025 - m[A25]*Det3_345_024;
1899 G4double Det4_2345_0345 = m[A20] * Det3_345 << 1811 G4double Det4_2345_0345 = m[A20]*Det3_345_345 - m[A23]*Det3_345_045
1900 m[A24] * Det3_345 << 1812 + m[A24]*Det3_345_035 - m[A25]*Det3_345_034;
1901 G4double Det4_2345_1234 = m[A21] * Det3_345 << 1813 G4double Det4_2345_1234 = m[A21]*Det3_345_234 - m[A22]*Det3_345_134
1902 m[A23] * Det3_345 << 1814 + m[A23]*Det3_345_124 - m[A24]*Det3_345_123;
1903 G4double Det4_2345_1235 = m[A21] * Det3_345 << 1815 G4double Det4_2345_1235 = m[A21]*Det3_345_235 - m[A22]*Det3_345_135
1904 m[A23] * Det3_345 << 1816 + m[A23]*Det3_345_125 - m[A25]*Det3_345_123;
1905 G4double Det4_2345_1245 = m[A21] * Det3_345 << 1817 G4double Det4_2345_1245 = m[A21]*Det3_345_245 - m[A22]*Det3_345_145
1906 m[A24] * Det3_345 << 1818 + m[A24]*Det3_345_125 - m[A25]*Det3_345_124;
1907 G4double Det4_2345_1345 = m[A21] * Det3_345 << 1819 G4double Det4_2345_1345 = m[A21]*Det3_345_345 - m[A23]*Det3_345_145
1908 m[A24] * Det3_345 << 1820 + m[A24]*Det3_345_135 - m[A25]*Det3_345_134;
1909 G4double Det4_2345_2345 = m[A22] * Det3_345 << 1821 G4double Det4_2345_2345 = m[A22]*Det3_345_345 - m[A23]*Det3_345_245
1910 m[A24] * Det3_345 << 1822 + m[A24]*Det3_345_235 - m[A25]*Det3_345_234;
1911 1823
1912 // Find all NECESSARY 5x5 dets: (19 of the 1824 // Find all NECESSARY 5x5 dets: (19 of them)
1913 1825
1914 G4double Det5_01234_01234 = << 1826 G4double Det5_01234_01234 = m[A00]*Det4_1234_1234 - m[A01]*Det4_1234_0234
1915 m[A00] * Det4_1234_1234 - m[A01] * Det4_1 << 1827 + m[A02]*Det4_1234_0134 - m[A03]*Det4_1234_0124 + m[A04]*Det4_1234_0123;
1916 m[A02] * Det4_1234_0134 - m[A03] * Det4_1 << 1828 G4double Det5_01235_01234 = m[A00]*Det4_1235_1234 - m[A01]*Det4_1235_0234
1917 G4double Det5_01235_01234 = << 1829 + m[A02]*Det4_1235_0134 - m[A03]*Det4_1235_0124 + m[A04]*Det4_1235_0123;
1918 m[A00] * Det4_1235_1234 - m[A01] * Det4_1 << 1830 G4double Det5_01235_01235 = m[A00]*Det4_1235_1235 - m[A01]*Det4_1235_0235
1919 m[A02] * Det4_1235_0134 - m[A03] * Det4_1 << 1831 + m[A02]*Det4_1235_0135 - m[A03]*Det4_1235_0125 + m[A05]*Det4_1235_0123;
1920 G4double Det5_01235_01235 = << 1832 G4double Det5_01245_01234 = m[A00]*Det4_1245_1234 - m[A01]*Det4_1245_0234
1921 m[A00] * Det4_1235_1235 - m[A01] * Det4_1 << 1833 + m[A02]*Det4_1245_0134 - m[A03]*Det4_1245_0124 + m[A04]*Det4_1245_0123;
1922 m[A02] * Det4_1235_0135 - m[A03] * Det4_1 << 1834 G4double Det5_01245_01235 = m[A00]*Det4_1245_1235 - m[A01]*Det4_1245_0235
1923 G4double Det5_01245_01234 = << 1835 + m[A02]*Det4_1245_0135 - m[A03]*Det4_1245_0125 + m[A05]*Det4_1245_0123;
1924 m[A00] * Det4_1245_1234 - m[A01] * Det4_1 << 1836 G4double Det5_01245_01245 = m[A00]*Det4_1245_1245 - m[A01]*Det4_1245_0245
1925 m[A02] * Det4_1245_0134 - m[A03] * Det4_1 << 1837 + m[A02]*Det4_1245_0145 - m[A04]*Det4_1245_0125 + m[A05]*Det4_1245_0124;
1926 G4double Det5_01245_01235 = << 1838 G4double Det5_01345_01234 = m[A00]*Det4_1345_1234 - m[A01]*Det4_1345_0234
1927 m[A00] * Det4_1245_1235 - m[A01] * Det4_1 << 1839 + m[A02]*Det4_1345_0134 - m[A03]*Det4_1345_0124 + m[A04]*Det4_1345_0123;
1928 m[A02] * Det4_1245_0135 - m[A03] * Det4_1 << 1840 G4double Det5_01345_01235 = m[A00]*Det4_1345_1235 - m[A01]*Det4_1345_0235
1929 G4double Det5_01245_01245 = << 1841 + m[A02]*Det4_1345_0135 - m[A03]*Det4_1345_0125 + m[A05]*Det4_1345_0123;
1930 m[A00] * Det4_1245_1245 - m[A01] * Det4_1 << 1842 G4double Det5_01345_01245 = m[A00]*Det4_1345_1245 - m[A01]*Det4_1345_0245
1931 m[A02] * Det4_1245_0145 - m[A04] * Det4_1 << 1843 + m[A02]*Det4_1345_0145 - m[A04]*Det4_1345_0125 + m[A05]*Det4_1345_0124;
1932 G4double Det5_01345_01234 = << 1844 G4double Det5_01345_01345 = m[A00]*Det4_1345_1345 - m[A01]*Det4_1345_0345
1933 m[A00] * Det4_1345_1234 - m[A01] * Det4_1 << 1845 + m[A03]*Det4_1345_0145 - m[A04]*Det4_1345_0135 + m[A05]*Det4_1345_0134;
1934 m[A02] * Det4_1345_0134 - m[A03] * Det4_1 << 1846 G4double Det5_02345_01234 = m[A00]*Det4_2345_1234 - m[A01]*Det4_2345_0234
1935 G4double Det5_01345_01235 = << 1847 + m[A02]*Det4_2345_0134 - m[A03]*Det4_2345_0124 + m[A04]*Det4_2345_0123;
1936 m[A00] * Det4_1345_1235 - m[A01] * Det4_1 << 1848 G4double Det5_02345_01235 = m[A00]*Det4_2345_1235 - m[A01]*Det4_2345_0235
1937 m[A02] * Det4_1345_0135 - m[A03] * Det4_1 << 1849 + m[A02]*Det4_2345_0135 - m[A03]*Det4_2345_0125 + m[A05]*Det4_2345_0123;
1938 G4double Det5_01345_01245 = << 1850 G4double Det5_02345_01245 = m[A00]*Det4_2345_1245 - m[A01]*Det4_2345_0245
1939 m[A00] * Det4_1345_1245 - m[A01] * Det4_1 << 1851 + m[A02]*Det4_2345_0145 - m[A04]*Det4_2345_0125 + m[A05]*Det4_2345_0124;
1940 m[A02] * Det4_1345_0145 - m[A04] * Det4_1 << 1852 G4double Det5_02345_01345 = m[A00]*Det4_2345_1345 - m[A01]*Det4_2345_0345
1941 G4double Det5_01345_01345 = << 1853 + m[A03]*Det4_2345_0145 - m[A04]*Det4_2345_0135 + m[A05]*Det4_2345_0134;
1942 m[A00] * Det4_1345_1345 - m[A01] * Det4_1 << 1854 G4double Det5_02345_02345 = m[A00]*Det4_2345_2345 - m[A02]*Det4_2345_0345
1943 m[A03] * Det4_1345_0145 - m[A04] * Det4_1 << 1855 + m[A03]*Det4_2345_0245 - m[A04]*Det4_2345_0235 + m[A05]*Det4_2345_0234;
1944 G4double Det5_02345_01234 = << 1856 G4double Det5_12345_01234 = m[A10]*Det4_2345_1234 - m[A11]*Det4_2345_0234
1945 m[A00] * Det4_2345_1234 - m[A01] * Det4_2 << 1857 + m[A12]*Det4_2345_0134 - m[A13]*Det4_2345_0124 + m[A14]*Det4_2345_0123;
1946 m[A02] * Det4_2345_0134 - m[A03] * Det4_2 << 1858 G4double Det5_12345_01235 = m[A10]*Det4_2345_1235 - m[A11]*Det4_2345_0235
1947 G4double Det5_02345_01235 = << 1859 + m[A12]*Det4_2345_0135 - m[A13]*Det4_2345_0125 + m[A15]*Det4_2345_0123;
1948 m[A00] * Det4_2345_1235 - m[A01] * Det4_2 << 1860 G4double Det5_12345_01245 = m[A10]*Det4_2345_1245 - m[A11]*Det4_2345_0245
1949 m[A02] * Det4_2345_0135 - m[A03] * Det4_2 << 1861 + m[A12]*Det4_2345_0145 - m[A14]*Det4_2345_0125 + m[A15]*Det4_2345_0124;
1950 G4double Det5_02345_01245 = << 1862 G4double Det5_12345_01345 = m[A10]*Det4_2345_1345 - m[A11]*Det4_2345_0345
1951 m[A00] * Det4_2345_1245 - m[A01] * Det4_2 << 1863 + m[A13]*Det4_2345_0145 - m[A14]*Det4_2345_0135 + m[A15]*Det4_2345_0134;
1952 m[A02] * Det4_2345_0145 - m[A04] * Det4_2 << 1864 G4double Det5_12345_02345 = m[A10]*Det4_2345_2345 - m[A12]*Det4_2345_0345
1953 G4double Det5_02345_01345 = << 1865 + m[A13]*Det4_2345_0245 - m[A14]*Det4_2345_0235 + m[A15]*Det4_2345_0234;
1954 m[A00] * Det4_2345_1345 - m[A01] * Det4_2 << 1866 G4double Det5_12345_12345 = m[A11]*Det4_2345_2345 - m[A12]*Det4_2345_1345
1955 m[A03] * Det4_2345_0145 - m[A04] * Det4_2 << 1867 + m[A13]*Det4_2345_1245 - m[A14]*Det4_2345_1235 + m[A15]*Det4_2345_1234;
1956 G4double Det5_02345_02345 = << 1868
1957 m[A00] * Det4_2345_2345 - m[A02] * Det4_2 << 1869 // Find the determinant
1958 m[A03] * Det4_2345_0245 - m[A04] * Det4_2 << 1870
1959 G4double Det5_12345_01234 = << 1871 G4double det = m[A00]*Det5_12345_12345
1960 m[A10] * Det4_2345_1234 - m[A11] * Det4_2 << 1872 - m[A01]*Det5_12345_02345
1961 m[A12] * Det4_2345_0134 - m[A13] * Det4_2 << 1873 + m[A02]*Det5_12345_01345
1962 G4double Det5_12345_01235 = << 1874 - m[A03]*Det5_12345_01245
1963 m[A10] * Det4_2345_1235 - m[A11] * Det4_2 << 1875 + m[A04]*Det5_12345_01235
1964 m[A12] * Det4_2345_0135 - m[A13] * Det4_2 << 1876 - m[A05]*Det5_12345_01234;
1965 G4double Det5_12345_01245 = <<
1966 m[A10] * Det4_2345_1245 - m[A11] * Det4_2 <<
1967 m[A12] * Det4_2345_0145 - m[A14] * Det4_2 <<
1968 G4double Det5_12345_01345 = <<
1969 m[A10] * Det4_2345_1345 - m[A11] * Det4_2 <<
1970 m[A13] * Det4_2345_0145 - m[A14] * Det4_2 <<
1971 G4double Det5_12345_02345 = <<
1972 m[A10] * Det4_2345_2345 - m[A12] * Det4_2 <<
1973 m[A13] * Det4_2345_0245 - m[A14] * Det4_2 <<
1974 G4double Det5_12345_12345 = <<
1975 m[A11] * Det4_2345_2345 - m[A12] * Det4_2 <<
1976 m[A13] * Det4_2345_1245 - m[A14] * Det4_2 <<
1977 <<
1978 // Find the determinant <<
1979 <<
1980 G4double det = m[A00] * Det5_12345_12345 - <<
1981 m[A02] * Det5_12345_01345 - <<
1982 m[A04] * Det5_12345_01235 - <<
1983 1877
1984 if(det == 0) << 1878 if ( det == 0 )
1985 { << 1879 {
1986 ifail = 1; 1880 ifail = 1;
1987 return; 1881 return;
1988 } << 1882 }
1989 1883
1990 G4double oneOverDet = 1.0 / det; << 1884 G4double oneOverDet = 1.0/det;
1991 G4double mn1OverDet = -oneOverDet; << 1885 G4double mn1OverDet = - oneOverDet;
1992 1886
1993 m[A00] = Det5_12345_12345 * oneOverDet; << 1887 m[A00] = Det5_12345_12345*oneOverDet;
1994 m[A01] = Det5_12345_02345 * mn1OverDet; << 1888 m[A01] = Det5_12345_02345*mn1OverDet;
1995 m[A02] = Det5_12345_01345 * oneOverDet; << 1889 m[A02] = Det5_12345_01345*oneOverDet;
1996 m[A03] = Det5_12345_01245 * mn1OverDet; << 1890 m[A03] = Det5_12345_01245*mn1OverDet;
1997 m[A04] = Det5_12345_01235 * oneOverDet; << 1891 m[A04] = Det5_12345_01235*oneOverDet;
1998 m[A05] = Det5_12345_01234 * mn1OverDet; << 1892 m[A05] = Det5_12345_01234*mn1OverDet;
1999 << 1893
2000 m[A11] = Det5_02345_02345 * oneOverDet; << 1894 m[A11] = Det5_02345_02345*oneOverDet;
2001 m[A12] = Det5_02345_01345 * mn1OverDet; << 1895 m[A12] = Det5_02345_01345*mn1OverDet;
2002 m[A13] = Det5_02345_01245 * oneOverDet; << 1896 m[A13] = Det5_02345_01245*oneOverDet;
2003 m[A14] = Det5_02345_01235 * mn1OverDet; << 1897 m[A14] = Det5_02345_01235*mn1OverDet;
2004 m[A15] = Det5_02345_01234 * oneOverDet; << 1898 m[A15] = Det5_02345_01234*oneOverDet;
2005 << 1899
2006 m[A22] = Det5_01345_01345 * oneOverDet; << 1900 m[A22] = Det5_01345_01345*oneOverDet;
2007 m[A23] = Det5_01345_01245 * mn1OverDet; << 1901 m[A23] = Det5_01345_01245*mn1OverDet;
2008 m[A24] = Det5_01345_01235 * oneOverDet; << 1902 m[A24] = Det5_01345_01235*oneOverDet;
2009 m[A25] = Det5_01345_01234 * mn1OverDet; << 1903 m[A25] = Det5_01345_01234*mn1OverDet;
2010 << 1904
2011 m[A33] = Det5_01245_01245 * oneOverDet; << 1905 m[A33] = Det5_01245_01245*oneOverDet;
2012 m[A34] = Det5_01245_01235 * mn1OverDet; << 1906 m[A34] = Det5_01245_01235*mn1OverDet;
2013 m[A35] = Det5_01245_01234 * oneOverDet; << 1907 m[A35] = Det5_01245_01234*oneOverDet;
2014 1908
2015 m[A44] = Det5_01235_01235 * oneOverDet; << 1909 m[A44] = Det5_01235_01235*oneOverDet;
2016 m[A45] = Det5_01235_01234 * mn1OverDet; << 1910 m[A45] = Det5_01235_01234*mn1OverDet;
2017 1911
2018 m[A55] = Det5_01234_01234 * oneOverDet; << 1912 m[A55] = Det5_01234_01234*oneOverDet;
2019 1913
2020 return; 1914 return;
2021 } 1915 }
2022 1916
2023 void G4ErrorSymMatrix::invertCholesky5(G4int& << 1917 void G4ErrorSymMatrix::invertCholesky5 (G4int & ifail)
2024 { 1918 {
2025 // Invert by << 1919
>> 1920 // Invert by
2026 // 1921 //
2027 // a) decomposing M = G*G^T with G lower tr 1922 // a) decomposing M = G*G^T with G lower triangular
2028 // (if M is not positive definite this w 1923 // (if M is not positive definite this will fail, leaving this unchanged)
2029 // b) inverting G to form H 1924 // b) inverting G to form H
2030 // c) multiplying H^T * H to get M^-1. 1925 // c) multiplying H^T * H to get M^-1.
2031 // 1926 //
2032 // If the matrix is pos. def. it is inverte 1927 // If the matrix is pos. def. it is inverted and 1 is returned.
2033 // If the matrix is not pos. def. it remain 1928 // If the matrix is not pos. def. it remains unaltered and 0 is returned.
2034 1929
2035 G4double h10; // below-diagonal elements o << 1930 G4double h10; // below-diagonal elements of H
2036 G4double h20, h21; 1931 G4double h20, h21;
2037 G4double h30, h31, h32; 1932 G4double h30, h31, h32;
2038 G4double h40, h41, h42, h43; 1933 G4double h40, h41, h42, h43;
2039 1934
2040 G4double h00, h11, h22, h33, h44; // 1/dia << 1935 G4double h00, h11, h22, h33, h44; // 1/diagonal elements of G =
2041 // diago << 1936 // diagonal elements of H
2042 1937
2043 G4double g10; // below-diagonal elements o << 1938 G4double g10; // below-diagonal elements of G
2044 G4double g20, g21; 1939 G4double g20, g21;
2045 G4double g30, g31, g32; 1940 G4double g30, g31, g32;
2046 G4double g40, g41, g42, g43; 1941 G4double g40, g41, g42, g43;
2047 1942
2048 ifail = 1; // We start by assuing we won't 1943 ifail = 1; // We start by assuing we won't succeed...
2049 1944
2050 // Form G -- compute diagonal members of H 1945 // Form G -- compute diagonal members of H directly rather than of G
2051 //------- 1946 //-------
2052 1947
2053 // Scale first column by 1/sqrt(A00) 1948 // Scale first column by 1/sqrt(A00)
2054 1949
2055 h00 = m[A00]; << 1950 h00 = m[A00];
2056 if(h00 <= 0) << 1951 if (h00 <= 0) { return; }
2057 { <<
2058 return; <<
2059 } <<
2060 h00 = 1.0 / std::sqrt(h00); 1952 h00 = 1.0 / std::sqrt(h00);
2061 1953
2062 g10 = m[A10] * h00; 1954 g10 = m[A10] * h00;
2063 g20 = m[A20] * h00; 1955 g20 = m[A20] * h00;
2064 g30 = m[A30] * h00; 1956 g30 = m[A30] * h00;
2065 g40 = m[A40] * h00; 1957 g40 = m[A40] * h00;
2066 1958
2067 // Form G11 (actually, just h11) 1959 // Form G11 (actually, just h11)
2068 1960
2069 h11 = m[A11] - (g10 * g10); 1961 h11 = m[A11] - (g10 * g10);
2070 if(h11 <= 0) << 1962 if (h11 <= 0) { return; }
2071 { <<
2072 return; <<
2073 } <<
2074 h11 = 1.0 / std::sqrt(h11); 1963 h11 = 1.0 / std::sqrt(h11);
2075 1964
2076 // Subtract inter-column column dot product 1965 // Subtract inter-column column dot products from rest of column 1 and
2077 // scale to get column 1 of G << 1966 // scale to get column 1 of G
2078 1967
2079 g21 = (m[A21] - (g10 * g20)) * h11; 1968 g21 = (m[A21] - (g10 * g20)) * h11;
2080 g31 = (m[A31] - (g10 * g30)) * h11; 1969 g31 = (m[A31] - (g10 * g30)) * h11;
2081 g41 = (m[A41] - (g10 * g40)) * h11; 1970 g41 = (m[A41] - (g10 * g40)) * h11;
2082 1971
2083 // Form G22 (actually, just h22) 1972 // Form G22 (actually, just h22)
2084 1973
2085 h22 = m[A22] - (g20 * g20) - (g21 * g21); 1974 h22 = m[A22] - (g20 * g20) - (g21 * g21);
2086 if(h22 <= 0) << 1975 if (h22 <= 0) { return; }
2087 { <<
2088 return; <<
2089 } <<
2090 h22 = 1.0 / std::sqrt(h22); 1976 h22 = 1.0 / std::sqrt(h22);
2091 1977
2092 // Subtract inter-column column dot product 1978 // Subtract inter-column column dot products from rest of column 2 and
2093 // scale to get column 2 of G << 1979 // scale to get column 2 of G
2094 1980
2095 g32 = (m[A32] - (g20 * g30) - (g21 * g31)) 1981 g32 = (m[A32] - (g20 * g30) - (g21 * g31)) * h22;
2096 g42 = (m[A42] - (g20 * g40) - (g21 * g41)) 1982 g42 = (m[A42] - (g20 * g40) - (g21 * g41)) * h22;
2097 1983
2098 // Form G33 (actually, just h33) 1984 // Form G33 (actually, just h33)
2099 1985
2100 h33 = m[A33] - (g30 * g30) - (g31 * g31) - 1986 h33 = m[A33] - (g30 * g30) - (g31 * g31) - (g32 * g32);
2101 if(h33 <= 0) << 1987 if (h33 <= 0) { return; }
2102 { <<
2103 return; <<
2104 } <<
2105 h33 = 1.0 / std::sqrt(h33); 1988 h33 = 1.0 / std::sqrt(h33);
2106 1989
2107 // Subtract inter-column column dot product 1990 // Subtract inter-column column dot product from A43 and scale to get G43
2108 1991
2109 g43 = (m[A43] - (g30 * g40) - (g31 * g41) - 1992 g43 = (m[A43] - (g30 * g40) - (g31 * g41) - (g32 * g42)) * h33;
2110 1993
2111 // Finally form h44 - if this is possible i 1994 // Finally form h44 - if this is possible inversion succeeds
2112 1995
2113 h44 = m[A44] - (g40 * g40) - (g41 * g41) - 1996 h44 = m[A44] - (g40 * g40) - (g41 * g41) - (g42 * g42) - (g43 * g43);
2114 if(h44 <= 0) << 1997 if (h44 <= 0) { return; }
2115 { <<
2116 return; <<
2117 } <<
2118 h44 = 1.0 / std::sqrt(h44); 1998 h44 = 1.0 / std::sqrt(h44);
2119 1999
2120 // Form H = 1/G -- diagonal members of H ar 2000 // Form H = 1/G -- diagonal members of H are already correct
2121 //------------- 2001 //-------------
2122 2002
2123 // The order here is dictated by speed cons 2003 // The order here is dictated by speed considerations
2124 2004
2125 h43 = -h33 * g43 * h44; << 2005 h43 = -h33 * g43 * h44;
2126 h32 = -h22 * g32 * h33; << 2006 h32 = -h22 * g32 * h33;
2127 h42 = -h22 * (g32 * h43 + g42 * h44); 2007 h42 = -h22 * (g32 * h43 + g42 * h44);
2128 h21 = -h11 * g21 * h22; << 2008 h21 = -h11 * g21 * h22;
2129 h31 = -h11 * (g21 * h32 + g31 * h33); 2009 h31 = -h11 * (g21 * h32 + g31 * h33);
2130 h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * 2010 h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * h44);
2131 h10 = -h00 * g10 * h11; << 2011 h10 = -h00 * g10 * h11;
2132 h20 = -h00 * (g10 * h21 + g20 * h22); 2012 h20 = -h00 * (g10 * h21 + g20 * h22);
2133 h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * 2013 h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * h33);
2134 h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * 2014 h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * h43 + g40 * h44);
2135 2015
2136 // Change this to its inverse = H^T*H 2016 // Change this to its inverse = H^T*H
2137 //------------------------------------ 2017 //------------------------------------
2138 2018
2139 m[A00] = h00 * h00 + h10 * h10 + h20 * h20 2019 m[A00] = h00 * h00 + h10 * h10 + h20 * h20 + h30 * h30 + h40 * h40;
2140 m[A01] = h10 * h11 + h20 * h21 + h30 * h31 2020 m[A01] = h10 * h11 + h20 * h21 + h30 * h31 + h40 * h41;
2141 m[A11] = h11 * h11 + h21 * h21 + h31 * h31 2021 m[A11] = h11 * h11 + h21 * h21 + h31 * h31 + h41 * h41;
2142 m[A02] = h20 * h22 + h30 * h32 + h40 * h42; 2022 m[A02] = h20 * h22 + h30 * h32 + h40 * h42;
2143 m[A12] = h21 * h22 + h31 * h32 + h41 * h42; 2023 m[A12] = h21 * h22 + h31 * h32 + h41 * h42;
2144 m[A22] = h22 * h22 + h32 * h32 + h42 * h42; 2024 m[A22] = h22 * h22 + h32 * h32 + h42 * h42;
2145 m[A03] = h30 * h33 + h40 * h43; 2025 m[A03] = h30 * h33 + h40 * h43;
2146 m[A13] = h31 * h33 + h41 * h43; 2026 m[A13] = h31 * h33 + h41 * h43;
2147 m[A23] = h32 * h33 + h42 * h43; 2027 m[A23] = h32 * h33 + h42 * h43;
2148 m[A33] = h33 * h33 + h43 * h43; 2028 m[A33] = h33 * h33 + h43 * h43;
2149 m[A04] = h40 * h44; 2029 m[A04] = h40 * h44;
2150 m[A14] = h41 * h44; 2030 m[A14] = h41 * h44;
2151 m[A24] = h42 * h44; 2031 m[A24] = h42 * h44;
2152 m[A34] = h43 * h44; 2032 m[A34] = h43 * h44;
2153 m[A44] = h44 * h44; 2033 m[A44] = h44 * h44;
2154 2034
2155 ifail = 0; 2035 ifail = 0;
2156 return; 2036 return;
2157 } 2037 }
2158 2038
2159 void G4ErrorSymMatrix::invertCholesky6(G4int& << 2039 void G4ErrorSymMatrix::invertCholesky6 (G4int & ifail)
2160 { 2040 {
2161 // Invert by << 2041 // Invert by
2162 // 2042 //
2163 // a) decomposing M = G*G^T with G lower tr 2043 // a) decomposing M = G*G^T with G lower triangular
2164 // (if M is not positive definite this w 2044 // (if M is not positive definite this will fail, leaving this unchanged)
2165 // b) inverting G to form H 2045 // b) inverting G to form H
2166 // c) multiplying H^T * H to get M^-1. 2046 // c) multiplying H^T * H to get M^-1.
2167 // 2047 //
2168 // If the matrix is pos. def. it is inverte 2048 // If the matrix is pos. def. it is inverted and 1 is returned.
2169 // If the matrix is not pos. def. it remain 2049 // If the matrix is not pos. def. it remains unaltered and 0 is returned.
2170 2050
2171 G4double h10; // below-diagonal elements o << 2051 G4double h10; // below-diagonal elements of H
2172 G4double h20, h21; 2052 G4double h20, h21;
2173 G4double h30, h31, h32; 2053 G4double h30, h31, h32;
2174 G4double h40, h41, h42, h43; 2054 G4double h40, h41, h42, h43;
2175 G4double h50, h51, h52, h53, h54; 2055 G4double h50, h51, h52, h53, h54;
2176 2056
2177 G4double h00, h11, h22, h33, h44, h55; // << 2057 G4double h00, h11, h22, h33, h44, h55; // 1/diagonal elements of G =
2178 // 2058 // diagonal elements of H
2179 2059
2180 G4double g10; // below-diagonal elements o << 2060 G4double g10; // below-diagonal elements of G
2181 G4double g20, g21; 2061 G4double g20, g21;
2182 G4double g30, g31, g32; 2062 G4double g30, g31, g32;
2183 G4double g40, g41, g42, g43; 2063 G4double g40, g41, g42, g43;
2184 G4double g50, g51, g52, g53, g54; 2064 G4double g50, g51, g52, g53, g54;
2185 2065
2186 ifail = 1; // We start by assuing we won't 2066 ifail = 1; // We start by assuing we won't succeed...
2187 2067
2188 // Form G -- compute diagonal members of H 2068 // Form G -- compute diagonal members of H directly rather than of G
2189 //------- 2069 //-------
2190 2070
2191 // Scale first column by 1/sqrt(A00) 2071 // Scale first column by 1/sqrt(A00)
2192 2072
2193 h00 = m[A00]; << 2073 h00 = m[A00];
2194 if(h00 <= 0) << 2074 if (h00 <= 0) { return; }
2195 { <<
2196 return; <<
2197 } <<
2198 h00 = 1.0 / std::sqrt(h00); 2075 h00 = 1.0 / std::sqrt(h00);
2199 2076
2200 g10 = m[A10] * h00; 2077 g10 = m[A10] * h00;
2201 g20 = m[A20] * h00; 2078 g20 = m[A20] * h00;
2202 g30 = m[A30] * h00; 2079 g30 = m[A30] * h00;
2203 g40 = m[A40] * h00; 2080 g40 = m[A40] * h00;
2204 g50 = m[A50] * h00; 2081 g50 = m[A50] * h00;
2205 2082
2206 // Form G11 (actually, just h11) 2083 // Form G11 (actually, just h11)
2207 2084
2208 h11 = m[A11] - (g10 * g10); 2085 h11 = m[A11] - (g10 * g10);
2209 if(h11 <= 0) << 2086 if (h11 <= 0) { return; }
2210 { <<
2211 return; <<
2212 } <<
2213 h11 = 1.0 / std::sqrt(h11); 2087 h11 = 1.0 / std::sqrt(h11);
2214 2088
2215 // Subtract inter-column column dot product 2089 // Subtract inter-column column dot products from rest of column 1 and
2216 // scale to get column 1 of G << 2090 // scale to get column 1 of G
2217 2091
2218 g21 = (m[A21] - (g10 * g20)) * h11; 2092 g21 = (m[A21] - (g10 * g20)) * h11;
2219 g31 = (m[A31] - (g10 * g30)) * h11; 2093 g31 = (m[A31] - (g10 * g30)) * h11;
2220 g41 = (m[A41] - (g10 * g40)) * h11; 2094 g41 = (m[A41] - (g10 * g40)) * h11;
2221 g51 = (m[A51] - (g10 * g50)) * h11; 2095 g51 = (m[A51] - (g10 * g50)) * h11;
2222 2096
2223 // Form G22 (actually, just h22) 2097 // Form G22 (actually, just h22)
2224 2098
2225 h22 = m[A22] - (g20 * g20) - (g21 * g21); 2099 h22 = m[A22] - (g20 * g20) - (g21 * g21);
2226 if(h22 <= 0) << 2100 if (h22 <= 0) { return; }
2227 { <<
2228 return; <<
2229 } <<
2230 h22 = 1.0 / std::sqrt(h22); 2101 h22 = 1.0 / std::sqrt(h22);
2231 2102
2232 // Subtract inter-column column dot product 2103 // Subtract inter-column column dot products from rest of column 2 and
2233 // scale to get column 2 of G << 2104 // scale to get column 2 of G
2234 2105
2235 g32 = (m[A32] - (g20 * g30) - (g21 * g31)) 2106 g32 = (m[A32] - (g20 * g30) - (g21 * g31)) * h22;
2236 g42 = (m[A42] - (g20 * g40) - (g21 * g41)) 2107 g42 = (m[A42] - (g20 * g40) - (g21 * g41)) * h22;
2237 g52 = (m[A52] - (g20 * g50) - (g21 * g51)) 2108 g52 = (m[A52] - (g20 * g50) - (g21 * g51)) * h22;
2238 2109
2239 // Form G33 (actually, just h33) 2110 // Form G33 (actually, just h33)
2240 2111
2241 h33 = m[A33] - (g30 * g30) - (g31 * g31) - 2112 h33 = m[A33] - (g30 * g30) - (g31 * g31) - (g32 * g32);
2242 if(h33 <= 0) << 2113 if (h33 <= 0) { return; }
2243 { <<
2244 return; <<
2245 } <<
2246 h33 = 1.0 / std::sqrt(h33); 2114 h33 = 1.0 / std::sqrt(h33);
2247 2115
2248 // Subtract inter-column column dot product 2116 // Subtract inter-column column dot products from rest of column 3 and
2249 // scale to get column 3 of G << 2117 // scale to get column 3 of G
2250 2118
2251 g43 = (m[A43] - (g30 * g40) - (g31 * g41) - 2119 g43 = (m[A43] - (g30 * g40) - (g31 * g41) - (g32 * g42)) * h33;
2252 g53 = (m[A53] - (g30 * g50) - (g31 * g51) - 2120 g53 = (m[A53] - (g30 * g50) - (g31 * g51) - (g32 * g52)) * h33;
2253 2121
2254 // Form G44 (actually, just h44) 2122 // Form G44 (actually, just h44)
2255 2123
2256 h44 = m[A44] - (g40 * g40) - (g41 * g41) - 2124 h44 = m[A44] - (g40 * g40) - (g41 * g41) - (g42 * g42) - (g43 * g43);
2257 if(h44 <= 0) << 2125 if (h44 <= 0) { return; }
2258 { <<
2259 return; <<
2260 } <<
2261 h44 = 1.0 / std::sqrt(h44); 2126 h44 = 1.0 / std::sqrt(h44);
2262 2127
2263 // Subtract inter-column column dot product 2128 // Subtract inter-column column dot product from M54 and scale to get G54
2264 2129
2265 g54 = (m[A54] - (g40 * g50) - (g41 * g51) - 2130 g54 = (m[A54] - (g40 * g50) - (g41 * g51) - (g42 * g52) - (g43 * g53)) * h44;
2266 2131
2267 // Finally form h55 - if this is possible i 2132 // Finally form h55 - if this is possible inversion succeeds
2268 2133
2269 h55 = m[A55] - (g50 * g50) - (g51 * g51) - << 2134 h55 = m[A55] - (g50*g50) - (g51*g51) - (g52*g52) - (g53*g53) - (g54*g54);
2270 (g54 * g54); << 2135 if (h55 <= 0) { return; }
2271 if(h55 <= 0) <<
2272 { <<
2273 return; <<
2274 } <<
2275 h55 = 1.0 / std::sqrt(h55); 2136 h55 = 1.0 / std::sqrt(h55);
2276 2137
2277 // Form H = 1/G -- diagonal members of H ar 2138 // Form H = 1/G -- diagonal members of H are already correct
2278 //------------- 2139 //-------------
2279 2140
2280 // The order here is dictated by speed cons 2141 // The order here is dictated by speed considerations
2281 2142
2282 h54 = -h44 * g54 * h55; << 2143 h54 = -h44 * g54 * h55;
2283 h43 = -h33 * g43 * h44; << 2144 h43 = -h33 * g43 * h44;
2284 h53 = -h33 * (g43 * h54 + g53 * h55); 2145 h53 = -h33 * (g43 * h54 + g53 * h55);
2285 h32 = -h22 * g32 * h33; << 2146 h32 = -h22 * g32 * h33;
2286 h42 = -h22 * (g32 * h43 + g42 * h44); 2147 h42 = -h22 * (g32 * h43 + g42 * h44);
2287 h52 = -h22 * (g32 * h53 + g42 * h54 + g52 * 2148 h52 = -h22 * (g32 * h53 + g42 * h54 + g52 * h55);
2288 h21 = -h11 * g21 * h22; << 2149 h21 = -h11 * g21 * h22;
2289 h31 = -h11 * (g21 * h32 + g31 * h33); 2150 h31 = -h11 * (g21 * h32 + g31 * h33);
2290 h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * 2151 h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * h44);
2291 h51 = -h11 * (g21 * h52 + g31 * h53 + g41 * 2152 h51 = -h11 * (g21 * h52 + g31 * h53 + g41 * h54 + g51 * h55);
2292 h10 = -h00 * g10 * h11; << 2153 h10 = -h00 * g10 * h11;
2293 h20 = -h00 * (g10 * h21 + g20 * h22); 2154 h20 = -h00 * (g10 * h21 + g20 * h22);
2294 h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * 2155 h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * h33);
2295 h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * 2156 h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * h43 + g40 * h44);
2296 h50 = -h00 * (g10 * h51 + g20 * h52 + g30 * 2157 h50 = -h00 * (g10 * h51 + g20 * h52 + g30 * h53 + g40 * h54 + g50 * h55);
2297 2158
2298 // Change this to its inverse = H^T*H 2159 // Change this to its inverse = H^T*H
2299 //------------------------------------ 2160 //------------------------------------
2300 2161
2301 m[A00] = << 2162 m[A00] = h00 * h00 + h10 * h10 + h20 * h20 + h30 * h30 + h40 * h40 + h50*h50;
2302 h00 * h00 + h10 * h10 + h20 * h20 + h30 * <<
2303 m[A01] = h10 * h11 + h20 * h21 + h30 * h31 2163 m[A01] = h10 * h11 + h20 * h21 + h30 * h31 + h40 * h41 + h50 * h51;
2304 m[A11] = h11 * h11 + h21 * h21 + h31 * h31 2164 m[A11] = h11 * h11 + h21 * h21 + h31 * h31 + h41 * h41 + h51 * h51;
2305 m[A02] = h20 * h22 + h30 * h32 + h40 * h42 2165 m[A02] = h20 * h22 + h30 * h32 + h40 * h42 + h50 * h52;
2306 m[A12] = h21 * h22 + h31 * h32 + h41 * h42 2166 m[A12] = h21 * h22 + h31 * h32 + h41 * h42 + h51 * h52;
2307 m[A22] = h22 * h22 + h32 * h32 + h42 * h42 2167 m[A22] = h22 * h22 + h32 * h32 + h42 * h42 + h52 * h52;
2308 m[A03] = h30 * h33 + h40 * h43 + h50 * h53; 2168 m[A03] = h30 * h33 + h40 * h43 + h50 * h53;
2309 m[A13] = h31 * h33 + h41 * h43 + h51 * h53; 2169 m[A13] = h31 * h33 + h41 * h43 + h51 * h53;
2310 m[A23] = h32 * h33 + h42 * h43 + h52 * h53; 2170 m[A23] = h32 * h33 + h42 * h43 + h52 * h53;
2311 m[A33] = h33 * h33 + h43 * h43 + h53 * h53; 2171 m[A33] = h33 * h33 + h43 * h43 + h53 * h53;
2312 m[A04] = h40 * h44 + h50 * h54; 2172 m[A04] = h40 * h44 + h50 * h54;
2313 m[A14] = h41 * h44 + h51 * h54; 2173 m[A14] = h41 * h44 + h51 * h54;
2314 m[A24] = h42 * h44 + h52 * h54; 2174 m[A24] = h42 * h44 + h52 * h54;
2315 m[A34] = h43 * h44 + h53 * h54; 2175 m[A34] = h43 * h44 + h53 * h54;
2316 m[A44] = h44 * h44 + h54 * h54; 2176 m[A44] = h44 * h44 + h54 * h54;
2317 m[A05] = h50 * h55; 2177 m[A05] = h50 * h55;
2318 m[A15] = h51 * h55; 2178 m[A15] = h51 * h55;
2319 m[A25] = h52 * h55; 2179 m[A25] = h52 * h55;
2320 m[A35] = h53 * h55; 2180 m[A35] = h53 * h55;
2321 m[A45] = h54 * h55; 2181 m[A45] = h54 * h55;
2322 m[A55] = h55 * h55; 2182 m[A55] = h55 * h55;
2323 2183
2324 ifail = 0; 2184 ifail = 0;
2325 return; 2185 return;
2326 } 2186 }
2327 2187
2328 void G4ErrorSymMatrix::invert4(G4int& ifail) << 2188 void G4ErrorSymMatrix::invert4 (G4int & ifail)
2329 { 2189 {
2330 ifail = 0; 2190 ifail = 0;
2331 2191
2332 // Find all NECESSARY 2x2 dets: (14 of the 2192 // Find all NECESSARY 2x2 dets: (14 of them)
2333 2193
2334 G4double Det2_12_01 = m[A10] * m[A21] - m[A << 2194 G4double Det2_12_01 = m[A10]*m[A21] - m[A11]*m[A20];
2335 G4double Det2_12_02 = m[A10] * m[A22] - m[A << 2195 G4double Det2_12_02 = m[A10]*m[A22] - m[A12]*m[A20];
2336 G4double Det2_12_12 = m[A11] * m[A22] - m[A << 2196 G4double Det2_12_12 = m[A11]*m[A22] - m[A12]*m[A21];
2337 G4double Det2_13_01 = m[A10] * m[A31] - m[A << 2197 G4double Det2_13_01 = m[A10]*m[A31] - m[A11]*m[A30];
2338 G4double Det2_13_02 = m[A10] * m[A32] - m[A << 2198 G4double Det2_13_02 = m[A10]*m[A32] - m[A12]*m[A30];
2339 G4double Det2_13_03 = m[A10] * m[A33] - m[A << 2199 G4double Det2_13_03 = m[A10]*m[A33] - m[A13]*m[A30];
2340 G4double Det2_13_12 = m[A11] * m[A32] - m[A << 2200 G4double Det2_13_12 = m[A11]*m[A32] - m[A12]*m[A31];
2341 G4double Det2_13_13 = m[A11] * m[A33] - m[A << 2201 G4double Det2_13_13 = m[A11]*m[A33] - m[A13]*m[A31];
2342 G4double Det2_23_01 = m[A20] * m[A31] - m[A << 2202 G4double Det2_23_01 = m[A20]*m[A31] - m[A21]*m[A30];
2343 G4double Det2_23_02 = m[A20] * m[A32] - m[A << 2203 G4double Det2_23_02 = m[A20]*m[A32] - m[A22]*m[A30];
2344 G4double Det2_23_03 = m[A20] * m[A33] - m[A << 2204 G4double Det2_23_03 = m[A20]*m[A33] - m[A23]*m[A30];
2345 G4double Det2_23_12 = m[A21] * m[A32] - m[A << 2205 G4double Det2_23_12 = m[A21]*m[A32] - m[A22]*m[A31];
2346 G4double Det2_23_13 = m[A21] * m[A33] - m[A << 2206 G4double Det2_23_13 = m[A21]*m[A33] - m[A23]*m[A31];
2347 G4double Det2_23_23 = m[A22] * m[A33] - m[A << 2207 G4double Det2_23_23 = m[A22]*m[A33] - m[A23]*m[A32];
2348 2208
2349 // Find all NECESSARY 3x3 dets: (10 of th 2209 // Find all NECESSARY 3x3 dets: (10 of them)
2350 2210
2351 G4double Det3_012_012 = << 2211 G4double Det3_012_012 = m[A00]*Det2_12_12 - m[A01]*Det2_12_02
2352 m[A00] * Det2_12_12 - m[A01] * Det2_12_02 << 2212 + m[A02]*Det2_12_01;
2353 G4double Det3_013_012 = << 2213 G4double Det3_013_012 = m[A00]*Det2_13_12 - m[A01]*Det2_13_02
2354 m[A00] * Det2_13_12 - m[A01] * Det2_13_02 << 2214 + m[A02]*Det2_13_01;
2355 G4double Det3_013_013 = << 2215 G4double Det3_013_013 = m[A00]*Det2_13_13 - m[A01]*Det2_13_03
2356 m[A00] * Det2_13_13 - m[A01] * Det2_13_03 << 2216 + m[A03]*Det2_13_01;
2357 G4double Det3_023_012 = << 2217 G4double Det3_023_012 = m[A00]*Det2_23_12 - m[A01]*Det2_23_02
2358 m[A00] * Det2_23_12 - m[A01] * Det2_23_02 << 2218 + m[A02]*Det2_23_01;
2359 G4double Det3_023_013 = << 2219 G4double Det3_023_013 = m[A00]*Det2_23_13 - m[A01]*Det2_23_03
2360 m[A00] * Det2_23_13 - m[A01] * Det2_23_03 << 2220 + m[A03]*Det2_23_01;
2361 G4double Det3_023_023 = << 2221 G4double Det3_023_023 = m[A00]*Det2_23_23 - m[A02]*Det2_23_03
2362 m[A00] * Det2_23_23 - m[A02] * Det2_23_03 << 2222 + m[A03]*Det2_23_02;
2363 G4double Det3_123_012 = << 2223 G4double Det3_123_012 = m[A10]*Det2_23_12 - m[A11]*Det2_23_02
2364 m[A10] * Det2_23_12 - m[A11] * Det2_23_02 << 2224 + m[A12]*Det2_23_01;
2365 G4double Det3_123_013 = << 2225 G4double Det3_123_013 = m[A10]*Det2_23_13 - m[A11]*Det2_23_03
2366 m[A10] * Det2_23_13 - m[A11] * Det2_23_03 << 2226 + m[A13]*Det2_23_01;
2367 G4double Det3_123_023 = << 2227 G4double Det3_123_023 = m[A10]*Det2_23_23 - m[A12]*Det2_23_03
2368 m[A10] * Det2_23_23 - m[A12] * Det2_23_03 << 2228 + m[A13]*Det2_23_02;
2369 G4double Det3_123_123 = << 2229 G4double Det3_123_123 = m[A11]*Det2_23_23 - m[A12]*Det2_23_13
2370 m[A11] * Det2_23_23 - m[A12] * Det2_23_13 << 2230 + m[A13]*Det2_23_12;
2371 2231
2372 // Find the 4x4 det: 2232 // Find the 4x4 det:
2373 2233
2374 G4double det = m[A00] * Det3_123_123 - m[A0 << 2234 G4double det = m[A00]*Det3_123_123
2375 m[A02] * Det3_123_013 - m[A0 << 2235 - m[A01]*Det3_123_023
>> 2236 + m[A02]*Det3_123_013
>> 2237 - m[A03]*Det3_123_012;
2376 2238
2377 if(det == 0) << 2239 if ( det == 0 )
2378 { << 2240 {
2379 ifail = 1; 2241 ifail = 1;
2380 return; 2242 return;
2381 } << 2243 }
>> 2244
>> 2245 G4double oneOverDet = 1.0/det;
>> 2246 G4double mn1OverDet = - oneOverDet;
2382 2247
2383 G4double oneOverDet = 1.0 / det; << 2248 m[A00] = Det3_123_123 * oneOverDet;
2384 G4double mn1OverDet = -oneOverDet; << 2249 m[A01] = Det3_123_023 * mn1OverDet;
>> 2250 m[A02] = Det3_123_013 * oneOverDet;
>> 2251 m[A03] = Det3_123_012 * mn1OverDet;
2385 2252
2386 m[A00] = Det3_123_123 * oneOverDet; <<
2387 m[A01] = Det3_123_023 * mn1OverDet; <<
2388 m[A02] = Det3_123_013 * oneOverDet; <<
2389 m[A03] = Det3_123_012 * mn1OverDet; <<
2390 <<
2391 m[A11] = Det3_023_023 * oneOverDet; <<
2392 m[A12] = Det3_023_013 * mn1OverDet; <<
2393 m[A13] = Det3_023_012 * oneOverDet; <<
2394 2253
2395 m[A22] = Det3_013_013 * oneOverDet; << 2254 m[A11] = Det3_023_023 * oneOverDet;
2396 m[A23] = Det3_013_012 * mn1OverDet; << 2255 m[A12] = Det3_023_013 * mn1OverDet;
>> 2256 m[A13] = Det3_023_012 * oneOverDet;
2397 2257
2398 m[A33] = Det3_012_012 * oneOverDet; << 2258 m[A22] = Det3_013_013 * oneOverDet;
>> 2259 m[A23] = Det3_013_012 * mn1OverDet;
>> 2260
>> 2261 m[A33] = Det3_012_012 * oneOverDet;
2399 2262
2400 return; 2263 return;
2401 } 2264 }
2402 2265
2403 void G4ErrorSymMatrix::invertHaywood4(G4int& << 2266 void G4ErrorSymMatrix::invertHaywood4 (G4int & ifail)
2404 { 2267 {
2405 invert4(ifail); // For the 4x4 case, the m << 2268 invert4(ifail); // For the 4x4 case, the method we use for invert is already
2406 // the Haywood method. << 2269 // the Haywood method.
2407 } 2270 }
>> 2271
2408 2272